Zero Displacement in Coupled system solve

dolfinx
1

Hello I am new to FEniCsx (running it via the WSL and Conda, dolfinx version: 0.9.0).
I am trying to solve the poroelastic system from “Computational Geomechanics with Special Reference to Earthquake Engineering, O. C. Zienkiewicz” given by :

\begin{align*} \overline{\sigma} &= \sigma - \alpha p \, \mathbf{Id} \\ \rho \frac{\partial^2 u}{\partial t^2} &- \nabla \cdot (\overline{\sigma}) - \rho b = 0 \\ \nabla \cdot \left(k (-\nabla p + \rho_f b )\right) &+ \nabla \cdot \left(\frac{\partial u}{\partial t} \right) + \frac{1}{Q} \frac{\partial p}{\partial t} = 0 \end{align*}

I tried to adapt the FEniCsx Tutorial for solving Navier Stokes using Taylor Hood and while the pressure is changing, u seems to be consistently zero

msh = create_rectangle( MPI.COMM_WORLD, [np.array([0, 0]), np.array([1.0, 1.0])], [30, 30], CellType.triangle)
def noslip_boundary(x):
    return np.isclose(x[1], 0.0)

P2 = element("Lagrange", msh.basix_cell(), 2, shape=(msh.geometry.dim,), dtype=default_real_type)
P1 = element("Lagrange", msh.basix_cell(), 1, dtype=default_real_type)
U, P = functionspace(msh, P2), functionspace(msh, P1)

# No-slip condition on boundaries where y = 0
u_noslip = np.array([0.0, 0.0], dtype=PETSc.ScalarType)
facets = locate_entities_boundary(msh, 1, noslip_boundary)
bc_u = dirichletbc(u_noslip, locate_dofs_topological(U, 1, facets), U)
outflow_dofs = fem.locate_dofs_geometrical(P, noslip_boundary)
bc_p = dirichletbc(PETSc.ScalarType(0), outflow_dofs, P)
bcs = [bc_u, bc_p]

I built the variational formulation like this

(u, p) = ufl.TrialFunction(U), ufl.TrialFunction(P)
u_dt = fem.Function(U)
u_ddt = fem.Function(U)
p_dt = fem.Function(P)
(u_test, p_test) = ufl.TestFunction(U),  ufl.TestFunction(P)

#set initial values
u_dt.x.array[:] = 0.0
u_ddt.x.array[:] = 0.0
p_dt.x.array[:] = 0.0 


dt_val = 1e-3
dt = fem.Constant(msh, dt_val)
b = fem.Constant(msh, PETSc.ScalarType((0.0, -9.81)))
k = 1e-10/0.00093
rho_f = 0.997
rho_s = 0.230
alpha = 1.0
phi = 0.77
KF = 1e2
KS= 1e4
Q_inv = phi/KF + (alpha-phi)/KS
rho = phi*rho_f + (1-phi)*rho_s
flux = -k * ufl.grad(p)
grav = rho_f * b

def epsilon( u):
        return 0.5*ufl.grad(u)+ufl.transpose((ufl.grad(u)))
def sigma(u):
        epsilon = epsilon(u)
        sigma = 2* 0.05*epsilon + 0.1 * ufl.tr(epsilon)* ufl.Identity(len(u))
        return sigma
sigma = sigma(u)

a00 = (rho * ufl.inner(u_test, (u/ dt**2)) * ufl.dx  + ufl.inner(sigma, ufl.grad(u_test)) * ufl.dx)
a01 =  - alpha *ufl.div(u_test)* p * ufl.dx 
a10 = p_test * ufl.div(u/dt) *ufl.dx
a11 = (ufl.inner(ufl.grad(p_test) , flux )* ufl.dx + Q_inv * p_test * (p/dt)* ufl.dx)
a = fem.form([ [a00, a01], [a10, a11]]) 

L0 = (-rho *  ufl.inner(u_test, ((-2*u_dt + u_ddt)/ dt**2) )* ufl.dx + rho *  ufl.inner(u_test, b) * ufl.dx)
L1 = (- p_test*(-p_dt/dt) * Q_inv * ufl.dx +  ufl.inner(ufl.grad(p_test) , (k* grav)) *ufl.dx - p_test * ufl.div(-u_dt/dt )* ufl.dx)
L = fem.form( [L0, L1])

a_p = fem.form([[a00, None], [None, a11]])

I then recreated the solve function from the tutorial (using the same “is_u” and “is_p” marker:

def solve():
    A, b, Pr = block_operators()
    ksp = PETSc.KSP().create(msh.comm)
    ksp.setOperators(A, Pr)
    ksp.setType("gmres")
    ksp.setTolerances(rtol=1.0e-6)
    #fieldsplit
    pc = ksp.getPC()
    pc.setType("fieldsplit")
    pc.setFieldSplitType(PETSc.PC.CompositeType.ADDITIVE)
    pc.setFieldSplitIS(("u", is_u), ("p", is_p))
    ksp_u, ksp_p = pc.getFieldSplitSubKSP()
    ksp_u.setType("preonly")
    ksp_u.getPC().setType("gamg")
    ksp_p.setType("preonly")
    ksp_p.getPC().setType("jacobi")
    pc.setUp()
    x = A.createVecRight()
    ksp.solve(b, x)
    return x

Which I use in the time loop:

uh = fem.Function(U)
ph = fem.Function(P)
t = 0.0
T = 1.
step = 0
while t < T:
    x = solve()
    # Extract displacement and pressure from solution vector
    ndofs_u = U.dofmap.index_map.size_local * U.dofmap.index_map_bs
    ndofs_p = P.dofmap.index_map.size_local
    uh.x.array[:ndofs_u] = x.array_r[:ndofs_u]
    ph.x.array[:ndofs_p] = x.array_r[ndofs_u:ndofs_u + ndofs_p]
    print("Max displacement:", la.norm(uh.x))

    # Update previous steps
    u_ddt.x.array[:] = u_dt.x.array  # u(t-2) = u(t-1)
    u_dt.x.array[:] = uh.x.array     # u(t-1) = u(t)
    p_dt.x.array[:] = ph.x.array     # p(t-1) = p(t)
    t += dt_val
    step += 1

While I see changes in p, u seems to be consistently zero. I tried other boundary conditions, and parameters but I cannot figure out what is wrong. Can someone help me?

Full code:

from petsc4py import PETSc
import numpy as np
import ufl
from basix.ufl import element
from dolfinx import default_real_type, fem, la
from dolfinx.fem import (
    dirichletbc,
    functionspace,
    locate_dofs_topological,
)
from dolfinx.fem.petsc import assemble_matrix_block, assemble_vector_block
from dolfinx.mesh import CellType, locate_entities_boundary, create_rectangle
import pyvista as pv
from dolfinx.plot import vtk_mesh
import dolfinx
print(dolfinx.__version__)

msh = create_rectangle( MPI.COMM_WORLD, [np.array([0, 0]), np.array([1.0, 1.0])], [30, 30], CellType.triangle)
def noslip_boundary(x):
    return np.isclose(x[1], 0.0)

P2 = element("Lagrange", msh.basix_cell(), 2, shape=(msh.geometry.dim,), dtype=default_real_type)
P1 = element("Lagrange", msh.basix_cell(), 1, dtype=default_real_type)
U, P = functionspace(msh, P2), functionspace(msh, P1)

# No-slip condition on boundaries where y = 0
u_noslip = np.array([0.0, 0.0], dtype=PETSc.ScalarType)
facets = locate_entities_boundary(msh, 1, noslip_boundary)
bc_u = dirichletbc(u_noslip, locate_dofs_topological(U, 1, facets), U)
outflow_dofs = fem.locate_dofs_geometrical(P, noslip_boundary)
bc_p = dirichletbc(PETSc.ScalarType(0), outflow_dofs, P)
bcs = [bc_u, bc_p]

(u, p) = ufl.TrialFunction(U), ufl.TrialFunction(P)
u_dt = fem.Function(U)
u_ddt = fem.Function(U)
p_dt = fem.Function(P)
(u_test, p_test) = ufl.TestFunction(U),  ufl.TestFunction(P)

#set initial values
u_dt.x.array[:] = 0.0
u_ddt.x.array[:] = 0.0
p_dt.x.array[:] = 0.0 

dt_val = 1e-3
dt = fem.Constant(msh, dt_val)
b = fem.Constant(msh, PETSc.ScalarType((0.0, -9.81)))
k = 1e-10/0.00093
rho_f = 0.997
rho_s = 0.230
alpha = 1.0
phi = 0.77
KF = 1e2
KS = 1e4
Q_inv = phi/KF + (alpha-phi)/KS
rho = phi* rho_f + (1-phi)* rho_s
flux = -k * ufl.grad(p)
grav = rho_f *b

def epsilon( u):
        return 0.5*ufl.grad(u)+ufl.transpose((ufl.grad(u)))
def sigma(u, epsilon):
        sigma = 2* 0.05*epsilon + 0.1 * ufl.tr(epsilon)* ufl.Identity(len(u))
        return sigma
epsilon = epsilon(u)
sigma = sigma(u, epsilon)
a00 = (rho * ufl.inner(u_test, (u/ dt**2)) * ufl.dx  + ufl.inner(sigma, ufl.grad(u_test)) * ufl.dx)
a01 =  - alpha *ufl.div(u_test)* p * ufl.dx 
a10 = p_test * ufl.div(u/dt) *ufl.dx
a11 = (ufl.inner(ufl.grad(p_test) , flux )* ufl.dx + Q_inv * p_test * (p/dt)* ufl.dx)
a = fem.form([ [a00, a01], [a10, a11]]) 

L0 = (-rho *  ufl.inner(u_test, ((-2*u_dt + u_ddt)/ dt**2) )* ufl.dx + rho *  ufl.inner(u_test, b) * ufl.dx)
L1 = (- p_test*(-p_dt/dt) * Q_inv * ufl.dx +  ufl.inner(ufl.grad(p_test) , (k* grav)) *ufl.dx - p_test * ufl.div(-u_dt/dt )* ufl.dx)
L = fem.form( [L0, L1])

a_p = fem.form([[a00, None], [None, a11]])


def block_operators():
    """Return block operators and block RHS vector for the Stokes
    problem"""
    A = assemble_matrix_block(a, bcs=bcs)
    A.assemble()
    Pr = assemble_matrix_block(a_p, bcs=bcs)
    Pr.assemble()
    b = assemble_vector_block(L, a, bcs=bcs)
    return A, b, Pr

IS = PETSc.IS()
U_map = U.dofmap.index_map
P_map = P.dofmap.index_map
offset_u = U_map.local_range[0] * U.dofmap.index_map_bs + P_map.local_range[0]
offset_p = offset_u + U_map.size_local * U.dofmap.index_map_bs
is_u = IS.createStride(
    U_map.size_local *  U.dofmap.index_map_bs, offset_u, 1, comm= PETSc.COMM_SELF
)
is_p = IS.createStride(
    P_map.size_local, offset_p, 1, comm= PETSc.COMM_SELF
)

ksp = PETSc.KSP().create(msh.comm)
ksp.setType("gmres")
ksp.setTolerances(rtol=1.0e-6)
#fieldsplit
pc = ksp.getPC()
pc.setType("fieldsplit")
pc.setFieldSplitType(PETSc.PC.CompositeType.ADDITIVE)
pc.setFieldSplitIS(("u", is_u), ("p", is_p))
ksp_u, ksp_p = pc.getFieldSplitSubKSP()
ksp_u.setType("preonly")
ksp_u.getPC().setType("gamg")
ksp_p.setType("preonly")
ksp_p.getPC().setType("jacobi")


def solve():
    A, b, Pr = block_operators()
    ksp.setOperators(A, Pr)
    x = A.createVecRight()
    ksp.solve(b, x)
    pc.setUp()
    return x 

uh = fem.Function(U)
ph = fem.Function(P)
t = 0.0
T = 1.
step = 0
while t < T:
    print(f"\nStep {step}, time {t:.4f}")
    x = solve()
    # Extract displacement and pressure from solution vector
    ndofs_u = U.dofmap.index_map.size_local * U.dofmap.index_map_bs
    ndofs_p = P.dofmap.index_map.size_local
    uh.x.array[:ndofs_u] = x.array_r[:ndofs_u]
    ph.x.array[:ndofs_p] = x.array_r[ndofs_u:ndofs_u + ndofs_p]

    # Update previous steps
    u_ddt.x.array[:] = u_dt.x.array  # u(t-2) = u(t-1)
    u_dt.x.array[:] = uh.x.array     # u(t-1) = u(t)
    p_dt.x.array[:] = ph.x.array     # p(t-1) = p(t)
    t += dt_val
    step += 1

print("Max displacement:", la.norm(uh.x))
print("Max pressure:", la.norm(ph.x))

topology, cell_types, geometry = vtk_mesh(U)
values = np.zeros((geometry.shape[0], 3), dtype=np.float64)
values[:, :len(uh)] = uh.x.array.real.reshape((geometry.shape[0], len(uh)))

# Create a point cloud of glyphs
function_grid = pv.UnstructuredGrid(topology, cell_types, geometry)
function_grid["u"] = values
glyphs = function_grid.glyph(orient="u", factor=0.2)

# Create a pyvista-grid for the mesh
msh.topology.create_connectivity(msh.topology.dim, msh.topology.dim)
grid = pv.UnstructuredGrid(*vtk_mesh(msh, msh.topology.dim))

# Create plotter
plotter = pv.Plotter()
plotter.add_mesh(grid, style="wireframe", color="k")
plotter.add_mesh(glyphs)
plotter.view_xy()

if pv.OFF_SCREEN:
    pv.start_xvfb(wait=0.1)
    plotter.screenshot("Deformation.png")
else:
    plotter.show()```
2

Could you rework your post to make sure all match renders correctly, and to create a runnable code snippet that people can copy and paste to test and help?

Some points I am noticing:

  • Most of the code in your solve loop can (and should) be outside of that loop as it does not change in time. In particular, creating the ksp and setting all options is not something you want to do over and over again.
  • You’re using minres as a solver, but your system is not symmetric. That shouldn’t work (and might be the issue). Try changing to gmres
3

Hi, thanks for checking. I updated the post with the full pasteable code. I tried changing the type to gmres, sadly that did not help.