Getting [inf inf ... inf] in this coupled displacement-pressure problem implemented with MixedElements

[gitable]

I am studying a soil model with coupling between its solid skeleton displacement \vec{u} and pore pressure p. The simplest form I can reduce the problem is with the form:

\nabla \cdot (\sigma(\vec{u}) - I p) + \rho \vec{g} = \vec{0},
\nabla^2 p= 0,

with linear elastic constitutive relation. Following the tutorials, I am trying to simultaneously solve the equations with MixedElements. Here is the MWE:

from mpi4py import MPI
from dolfinx import mesh
from dolfinx.fem import FunctionSpace, Function, dirichletbc, locate_dofs_topological, Constant
from dolfinx.fem.petsc import LinearProblem
from petsc4py.PETSc import ScalarType
import ufl
import numpy as np

def bottom_boundary(x):
    return np.isclose(x[1],0)

def left_boundary(x):
    return np.isclose(x[0],0)

def right_boundary(x):
    return np.isclose(x[0],0)

def epsilon(u):
    return ufl.sym(ufl.grad(u))

def sigma(u,Lame_params):
    return Lame_params[0]*ufl.nabla_div(u)*ufl.Identity(u.geometric_dimension()) + 2*Lame_params[1]*epsilon(u)


#Main simuluation
rho = 2000   #density
Ym = 2e7     #Young modulus
Pr = 0.3     #Poisson ratio
g = 9.8      #accel due to gravity
rho_w = 1000 #density of water
Lame_params = (Ym*Pr)/((1+Pr)*(1-2*Pr)), 0.5*Ym/(1+Pr)

domain = mesh.create_unit_square(MPI.COMM_WORLD, 10, 10, mesh.CellType.quadrilateral)

V_el = ufl.VectorElement("CG", domain.ufl_cell(), 1)
X_el = ufl.FiniteElement("CG", domain.ufl_cell(), 1)
M_el = ufl.MixedElement([V_el, X_el])
VX = FunctionSpace(domain,M_el)

(u,p) = ufl.TrialFunctions(VX)
(v,q) = ufl.TestFunctions(VX)

bcs = []
#boundary conditions on displacement
V, _ = VX.sub(0).collapse()
X, _ = VX.sub(1).collapse()

fdim = domain.topology.dim - 1
p0 = Function(X)
p0.interpolate(lambda x: -rho_w*g*x[0])
boundary_facets = mesh.locate_entities_boundary(domain, fdim, bottom_boundary)
bcs.append(dirichletbc(ScalarType((0,0)),locate_dofs_topological(V, fdim, boundary_facets), V))
bcs.append(dirichletbc(p0,locate_dofs_topological(X, fdim, boundary_facets)))

boundary_facets = mesh.locate_entities_boundary(domain, fdim, left_boundary)
bcs.append(dirichletbc(ScalarType((0,0)),locate_dofs_topological(V, fdim, boundary_facets), V))
bcs.append(dirichletbc(p0,locate_dofs_topological(X, fdim, boundary_facets)))

boundary_facets = mesh.locate_entities_boundary(domain, fdim, right_boundary)
bcs.append(dirichletbc(p0,locate_dofs_topological(X, fdim, boundary_facets)))

G = Constant(domain,ScalarType((0,-g)))
a = ufl.inner(sigma(u,Lame_params),epsilon(v))*ufl.dx
a += -ufl.inner(p*ufl.Identity(u.geometric_dimension()),epsilon(v))*ufl.dx
a += ufl.dot(ufl.grad(p),ufl.grad(q))*ufl.dx
L = ufl.dot(rho*G,v)*ufl.dx + Constant(domain,ScalarType(0))*q*ufl.dx

problem = LinearProblem(a,L,bcs,petsc_options={"ksp_type":"preonly","pc_type":"lu"})
up = problem.solve()

u,p = up.split()
print(u.x.array)
print(p.x.array)

to which I obtain u = [inf inf…] and p = [inf inf …].

Now, solving the equations by first solving p from the 2nd equation and then \vec{u} from the 1st equation with the obtained p yields a solution. The displacement (left) and pressure (right) plots are shown below:

Here is the MWE with this approach:

from mpi4py import MPI
from dolfinx import mesh
from dolfinx.fem import (FunctionSpace, Function, dirichletbc, locate_dofs_topological, Constant,
                        VectorFunctionSpace)
from dolfinx.fem.petsc import LinearProblem
from petsc4py.PETSc import ScalarType
import ufl
import numpy as np
from dolfinx.io import XDMFFile

def bottom_boundary(x):
    return np.isclose(x[1],0)

def left_boundary(x):
    return np.isclose(x[0],0)

def right_boundary(x):
    return np.isclose(x[0],1)

def epsilon(u):
    return ufl.sym(ufl.grad(u))

def sigma(u,Lame_params):
    return Lame_params[0]*ufl.nabla_div(u)*ufl.Identity(u.geometric_dimension()) + 2*Lame_params[1]*epsilon(u)

#Main simuluation
rho = 2000   #density
Ym = 2e6     #Young modulus
Pr = 0.3     #Poisson ratio
g = 9.8      #accel due to gravity
rho_w = 1000 #density of water
Lame_params = (Ym*Pr)/((1+Pr)*(1-2*Pr)), 0.5*Ym/(1+Pr)


domain = mesh.create_unit_square(MPI.COMM_WORLD, 10, 10, mesh.CellType.quadrilateral)

V = VectorFunctionSpace(domain,("CG",1))
u, v = ufl.TrialFunction(V), ufl.TestFunction(V)
X = FunctionSpace(domain,("CG",1))
p, q = ufl.TrialFunction(X), ufl.TestFunction(X)

bc_u = []
bc_p = []

fdim = domain.topology.dim - 1
p0 = Function(X)
p0.interpolate(lambda x: -rho_w*g*x[1])
boundary_facets = mesh.locate_entities_boundary(domain, fdim, bottom_boundary)
bc_u.append(dirichletbc(ScalarType((0,0)),locate_dofs_topological(V, fdim, boundary_facets), V))
bc_p.append(dirichletbc(p0,locate_dofs_topological(X, fdim, boundary_facets)))

boundary_facets = mesh.locate_entities_boundary(domain, fdim, left_boundary)
bc_u.append(dirichletbc(ScalarType((0,0)),locate_dofs_topological(V, fdim, boundary_facets), V))
bc_p.append(dirichletbc(p0,locate_dofs_topological(X, fdim, boundary_facets)))

boundary_facets = mesh.locate_entities_boundary(domain, fdim, right_boundary)
bc_p.append(dirichletbc(p0,locate_dofs_topological(X, fdim, boundary_facets)))

ap = ufl.dot(ufl.grad(p),ufl.grad(q))*ufl.dx
Lp = Constant(domain,ScalarType(0))*q*ufl.dx
problem = LinearProblem(ap,Lp,bcs=bc_p,petsc_options={"ksp_type":"preonly","pc_type":"lu"})
ph = problem.solve()

rho_g = Constant(domain,ScalarType((0,-rho*g)))
au = ufl.inner(sigma(u,Lame_params),epsilon(v))*ufl.dx
Lu = ufl.inner(ph*ufl.Identity(u.geometric_dimension()),epsilon(v))*ufl.dx + ufl.dot(rho_g,v)*ufl.dx
problem = LinearProblem(au,Lu,bcs=bc_u,petsc_options={"ksp_type":"preonly","pc_type":"lu"})
uh = problem.solve()

xdmf_p = XDMFFile(domain.comm, "pore_pressures.xdmf", "w")
xdmf_p.write_mesh(domain)
xdmf_p.write_function(ph,0)

xdmf_u = XDMFFile(domain.comm, "displacements.xdmf", "w")
xdmf_u.write_mesh(domain)
xdmf_u.write_function(uh,0)

Please check my scripts and explain why no solution is obtained in the MixedElements approach. Any comments will be appreciated.

[conpierce8]

You should not apply your bcs to the collapsed subspaces V and X, see:

Please help me to check whether boundary condition are right set

@skatulla Ah! You caught me not paying close attention. It looks like you should not apply your Dirichlet bcs to the collapsed subspaces, i.e. replacing bc1 = fem.dirichletbc(fixed_displacement, left_u_dofs, V0) bc2 = fem.dirichletbc(incremented_displacement, right_u_dofs, V0) bc3 = fem.dirichletbc(fixed_director, top_wf_dofs, V1) bc4 = fem.dirichletbc(fixed_director, bottom_wf_dofs,V1) with bc1 = fem.dirichletbc(fixed_displacement, left_u_dofs, VY.sub(0)) bc2 = fem.dirichletbc…

In your case, changing your first MWE to the following code should allow you to solve using MixedElements. Note that I corrected the definitions of right_boundary and p0 to match the definitions in your second example.

from mpi4py import MPI
from dolfinx import mesh
from dolfinx.fem import FunctionSpace, Function, dirichletbc, locate_dofs_topological, Constant
from dolfinx.fem.petsc import LinearProblem
from petsc4py.PETSc import ScalarType
import ufl
import numpy as np

def bottom_boundary(x):
    return np.isclose(x[1],0)

def left_boundary(x):
    return np.isclose(x[0],0)

def right_boundary(x):
    return np.isclose(x[0],1)

def epsilon(u):
    return ufl.sym(ufl.grad(u))

def sigma(u,Lame_params):
    return Lame_params[0]*ufl.nabla_div(u)*ufl.Identity(u.geometric_dimension()) + 2*Lame_params[1]*epsilon(u)


#Main simuluation
rho = 2000   #density
Ym = 2e7     #Young modulus
Pr = 0.3     #Poisson ratio
g = 9.8      #accel due to gravity
rho_w = 1000 #density of water
Lame_params = (Ym*Pr)/((1+Pr)*(1-2*Pr)), 0.5*Ym/(1+Pr)

domain = mesh.create_unit_square(MPI.COMM_WORLD, 10, 10, mesh.CellType.quadrilateral)

V_el = ufl.VectorElement("CG", domain.ufl_cell(), 1)
X_el = ufl.FiniteElement("CG", domain.ufl_cell(), 1)
M_el = ufl.MixedElement([V_el, X_el])
VX = FunctionSpace(domain,M_el)

(u,p) = ufl.TrialFunctions(VX)
(v,q) = ufl.TestFunctions(VX)

bcs = []
#boundary conditions on displacement
V, V_dofs = VX.sub(0).collapse()
X, X_dofs = VX.sub(1).collapse()

fdim = domain.topology.dim - 1
p0 = Function(X)
p0.interpolate(lambda x: -rho_w*g*x[1])
u0 = Function(V)
boundary_facets = mesh.locate_entities_boundary(domain, fdim, bottom_boundary)
bcs.append(dirichletbc(u0,locate_dofs_topological((VX.sub(0), V), fdim, boundary_facets), VX.sub(0)))
bcs.append(dirichletbc(p0,locate_dofs_topological((VX.sub(1), X), fdim, boundary_facets), VX.sub(1)))

boundary_facets = mesh.locate_entities_boundary(domain, fdim, left_boundary)
bcs.append(dirichletbc(u0,locate_dofs_topological((VX.sub(0), V), fdim, boundary_facets), VX.sub(0)))
bcs.append(dirichletbc(p0,locate_dofs_topological((VX.sub(1), X), fdim, boundary_facets), VX.sub(1)))

boundary_facets = mesh.locate_entities_boundary(domain, fdim, right_boundary)
bcs.append(dirichletbc(p0,locate_dofs_topological((VX.sub(1), X), fdim, boundary_facets), VX.sub(1)))

G = Constant(domain,ScalarType((0,-g)))
a = ufl.inner(sigma(u,Lame_params),epsilon(v))*ufl.dx
a += -ufl.inner(p*ufl.Identity(u.geometric_dimension()),epsilon(v))*ufl.dx
a += ufl.dot(ufl.grad(p),ufl.grad(q))*ufl.dx
L = ufl.dot(rho*G,v)*ufl.dx + Constant(domain,ScalarType(0))*q*ufl.dx

problem = LinearProblem(a,L,bcs,petsc_options={"ksp_type":"preonly","pc_type":"lu"})
up = problem.solve()

u,p = up.split()
print(u.x.array)
print(p.x.array)

from dolfinx import io
with io.XDMFFile(MPI.COMM_WORLD, "u.xdmf", "w") as xdmf:
    xdmf.write_mesh(domain)
    xdmf.write_function(u)
with io.XDMFFile(MPI.COMM_WORLD, "p.xdmf", "w") as xdmf:
    xdmf.write_mesh(domain)
    xdmf.write_function(p)

[gitable]

This solves my problem. Thanks!