# Solving coupling equation with MixedElement form

Hello I am now stuck in the implementation of the coupling equation with FEniCSx.
It may be quite a long question, but it would be a great pleasure if you comment on this problem. I will explain detailedly as much as I can.

1) Coupling equations

The problem consists of two coupling equations, formulated with u(displacement) and p(pressure).

The goal is to find the u and p, by solving these coupling equations with FEM.

2) Boundary conditions

For a simpler statement, here the domain is defined as rectangular mesh.
Here, it is imposed as :

• u : Dirichlet condition(=0) for left / right / bottom boundary
• p : Dicirchlet condition(=0) for top boundary
• n dot sigma (called traction) : Dicirchlet condition(=t0) only for top boundary

3) Code implementation

Here, I will show you the code. If you find any mistakes, please notify me.

• Defining Domain
``````domain = mesh.create_unit_square(MPI.COMM_WORLD, 10, 10)
x = SpatialCoordinate(domain)

dim = domain.topology.dim
bdim = dim - 1
``````
• Defining function space
``````# Define function spaces
u_elem     = VectorElement('CG', domain.ufl_cell(), 2, dim=2) # displacement(2D)
p_elem     = FiniteElement('CG', domain.ufl_cell(), 1) # pore water pressure(1D)
V    = FunctionSpace(domain, MixedElement([u_elem, p_elem])) #Define vector function space for [u,p]
``````
• Defining boundary conditions
``````# Define boundary conditions
def left(x):
return np.isclose(x, 0) #x=0

def right(x):
return np.isclose(x, 1) #x=1

def bottom(x):
return np.isclose(x, 0) #y=0

def top(x):
return np.isclose(x, 0) #y=1

# Constrained boundaries

#Dirichlet BC for u : left, right, bottom

U, _ = V.sub(0).collapse()
u_bc = Function(U)
u_bc.x.set(0)

#left
boundary_facets_left = locate_entities_boundary(domain, bdim, left)
boundary_dofs_left = locate_dofs_topological((V.sub(0), U), bdim, boundary_facets_left)
BC_u_left = dirichletbc(u_bc, boundary_dofs_left, V.sub(0))

#right
boundary_facets_right = locate_entities_boundary(domain, bdim, left)
boundary_dofs_right = locate_dofs_topological((V.sub(0), U), bdim, boundary_facets_right)
BC_u_right = dirichletbc(u_bc, boundary_dofs_right, V.sub(0))

#bottom
boundary_facets_bottom = locate_entities_boundary(domain, bdim, bottom)
boundary_dofs_bottom = locate_dofs_topological((V.sub(0), U), bdim, boundary_facets_bottom)
BC_u_bottom = dirichletbc(u_bc, boundary_dofs_bottom, V.sub(0))

#Dirichlet BC for p : top
P, _ = V.sub(1).collapse()
p_bc = Function(P)
p_bc.x.set(0)

boundary_facets_top = locate_entities_boundary(domain, bdim, top)
boundary_dofs_top = locate_dofs_topological((V.sub(1), P), bdim, boundary_facets_top)
BC_p_top = dirichletbc(p_bc, boundary_dofs_top, V.sub(1))

#Total BC
BC = [BC_u_left, BC_u_right, BC_u_bottom, BC_p_top]

trac = Constant(domain, ScalarType((0,trac0))) #Traction applied on top boundary, in downward direction

# Mark boundaries : 1, 2, 3, 4
marked_facets = np.hstack([boundary_facets_left, boundary_facets_right, boundary_facets_bottom, ])
marked_values = np.hstack([np.full_like(boundary_facets_left, 1), np.full_like(boundary_facets_right, 2), np.full_like(boundary_facets_bottom, 3), np.full_like(boundary_facets_top, 4)])
sorted_facets = np.argsort(marked_facets)
facet_tag = mesh.meshtags(domain, bdim, marked_facets[sorted_facets], marked_values[sorted_facets])

#Boundary (ds)
ds = Measure("ds", domain=domain, subdomain_data=facet_tag)
``````

For this part, I mainly followed codes in here: Boundary conditions about Mixed space - #2 by dokken

• Variational formulation and solver
``````del_x     = TestFunction(V)
eta, psi = split(del_x)

x_n1 = Function(V)
u_n1, p_n1 = split(x_n1)

#initial condition
x_n = Function(V)
u_n, p_n = split(x_n)

Res = inner(grad(eta), sigma_e(u_n1)) * dx - inner(p_n1, div(eta)) * dx - dot(eta, trac) * ds(4) \
+ inner(psi, div((u_n1 - u_n)/dt)) * dx - inner(grad(psi), w_Darcy(p_n1)) * dx

Jac = derivative(Res, x_n1) # jacobian

Prob   = NonlinearProblem(Res, x_n1, BC, Jac)
solver = nls.petsc.NewtonSolver(MPI.COMM_WORLD, Prob)
``````

Here, by running the cell, error occurs :

``````---------------------------------------------------------------------------
ValueError                                Traceback (most recent call last)
Cell In, line 22
17 Res = inner(grad(eta), sigma_e(u_n1)) * dx - inner(p_n1, div(eta)) * dx - dot(eta, trac) * ds(4) \
18         + inner(psi, div((u_n1 - u_n)/dt)) * dx - inner(grad(psi), w_Darcy(p_n1)) * dx
20 Jac = derivative(Res, x_n1) # jacobian
---> 22 Prob   = NonlinearProblem(Res, x_n1, BC, Jac)
23 solver = nls.petsc.NewtonSolver(MPI.COMM_WORLD, Prob)
25 # Set nonlinear solver parameters

File /usr/local/dolfinx-real/lib/python3.10/dist-packages/dolfinx/fem/petsc.py:683, in NonlinearProblem.__init__(self, F, u, bcs, J, form_compiler_options, jit_options)
661 def __init__(self, F: ufl.form.Form, u: _Function, bcs: typing.List[DirichletBCMetaClass] = [],
662              J: ufl.form.Form = None, form_compiler_options={}, jit_options={}):
663     """Initialize solver for solving a non-linear problem using Newton's method, :math:`dF/du(u) du = -F(u)`.
664
665     Args:
(...)
681
682     """
--> 683     self._L = _create_form(F, form_compiler_options=form_compiler_options,
684                            jit_options=jit_options)
686     # Create the Jacobian matrix, dF/du
687     if J is None:

File /usr/local/dolfinx-real/lib/python3.10/dist-packages/dolfinx/fem/forms.py:176, in form(form, dtype, form_compiler_options, jit_options)
173         return list(map(lambda sub_form: _create_form(sub_form), form))
174     return form
--> 176 return _create_form(form)

File /usr/local/dolfinx-real/lib/python3.10/dist-packages/dolfinx/fem/forms.py:171, in form.<locals>._create_form(form)
168 """Recursively convert ufl.Forms to dolfinx.fem.Form, otherwise
169 return form argument"""
170 if isinstance(form, ufl.Form):
--> 171     return _form(form)
172 elif isinstance(form, collections.abc.Iterable):
173     return list(map(lambda sub_form: _create_form(sub_form), form))

File /usr/local/dolfinx-real/lib/python3.10/dist-packages/dolfinx/fem/forms.py:132, in form.<locals>._form(form)
130 # Extract subdomain data from UFL form
131 sd = form.subdomain_data()
--> 132 domain, = list(sd.keys())  # Assuming single domain
133 # Get subdomain data for each integral type
134 subdomains = {}

ValueError: too many values to unpack (expected 1)
``````

I am really not sure about this part. Actually, I am now fixing the code written in FEniCS, so I am highly suspicious on some mismatch in grammar.

1 Like

Hi there,

I had a quick glance but figure out the error. However, it could be useful to check out GitHub - Matt-L-McLean/Computational-Geomechanics: poroelastic numerical simulations with open source software Fenicsx/Dolfinx, which implements similar problems.

Cheers

Hello Mr. Driley, first, I’d like to say thank you for your reply on old post. I really waited for some replies Actually, after posting the question, I revised a lot of code and could eliminate errors in that part.
Below is the revised code :

``````domain = mesh.create_unit_square(MPI.COMM_WORLD, 10, 10)
x = SpatialCoordinate(domain)

dim = domain.topology.dim
bdim = dim - 1
``````
``````u_elem     = VectorElement('CG', domain.ufl_cell(), 2) # displacement(2D)
p_elem     = FiniteElement('CG', domain.ufl_cell(), 1) # pore water pressure(scalar)
V    = FunctionSpace(domain, MixedElement([u_elem, p_elem])) #Define vector function space for [u,p]
``````
``````# Define boundary conditions
def left(x):
return np.isclose(x, 0) #x=0

def right(x):
return np.isclose(x, 1) #x=1

def bottom(x):
return np.isclose(x, 0) #y=0

def top(x):
return np.isclose(x, 1) #y=1

U, _ = V.sub(0).collapse()
u_bc = Function(U)
u_bc.x.set(0)

#left
boundary_facets_left = locate_entities_boundary(domain, bdim, left)
boundary_dofs_left = locate_dofs_topological((V.sub(0), U), bdim, boundary_facets_left)
BC_u_left = dirichletbc(u_bc, boundary_dofs_left, V.sub(0))

#right
boundary_facets_right = locate_entities_boundary(domain, bdim, right)
boundary_dofs_right = locate_dofs_topological((V.sub(0), U), bdim, boundary_facets_right)
BC_u_right = dirichletbc(u_bc, boundary_dofs_right, V.sub(0))

#bottom
boundary_facets_bottom = locate_entities_boundary(domain, bdim, bottom)
boundary_dofs_bottom = locate_dofs_topological((V.sub(0), U), bdim, boundary_facets_bottom)
BC_u_bottom = dirichletbc(u_bc, boundary_dofs_bottom, V.sub(0))

#Dirichlet BC for p : top
P, _ = V.sub(1).collapse()
p_bc = Function(P)
p_bc.x.set(0)

boundary_facets_top = locate_entities_boundary(domain, bdim, top)
boundary_dofs_top = locate_dofs_topological((V.sub(1), P), bdim, boundary_facets_top)
BC_p_top = dirichletbc(p_bc, boundary_dofs_top, V.sub(1))

#Total BC
BC = [BC_u_left, BC_u_right, BC_u_bottom, BC_p_top]

trac = Constant(domain, ScalarType((0,trac0))) #Traction applied on top boundary, in downward direction

# Mark boundaries : 1, 2, 3, 4
marked_facets = np.hstack([boundary_facets_left, boundary_facets_right, boundary_facets_bottom, boundary_facets_top])
marked_values = np.hstack([np.full_like(boundary_facets_left, 1), np.full_like(boundary_facets_right, 2), np.full_like(boundary_facets_bottom, 3), np.full_like(boundary_facets_top, 4)])
sorted_facets = np.argsort(marked_facets)
facet_tag = mesh.meshtags(domain, bdim, marked_facets[sorted_facets], marked_values[sorted_facets])

#Boundary (ds)
ds = Measure("ds", domain=domain, subdomain_data=facet_tag)
n  = FacetNormal(domain)
``````

I revised the boundary condition for pressure on top (originally it was bottom).

``````#Test function
eta, psi = ufl.TestFunctions(V)

#Trial function
xt = Function(V)
ut, pt = split(xt)

#Define function for step
x_n = Function(V)
#x_n.sub(1).interpolate(lambda x : np.full((x.shape,), 0))
print(x_n.vector.array)
u_n, p_n = split(x_n)

Res = inner(grad(eta), sigma_e(ut)) * dx - inner(pt, div(eta)) * dx  - dot(eta, trac) * ds(4) \
+ inner(psi, div((ut - u_n)/dt)) * dx - inner(grad(psi), w_Darcy(pt)) * dx

Jac = derivative(Res, xt) # jacobian
Prob   = NonlinearProblem(Res, xt, BC, Jac)
solver = nls.petsc.NewtonSolver(MPI.COMM_WORLD, Prob)

# Set nonlinear solver parameters
solver.rtol = 1e-7 #relative_tolerance
solver.atol = 1e-7 #absolute_tolerance
solver.max_it = 100 #maximum_iterations
solver.report = True
``````

Now, there is no error. However, when I tried to solve it with defined solver…

``````# Solve system & output results
# ----------------------------------------------------------------
# time stepping
time = np.linspace(t_i, t_f, Nsteps+1)

# output file
#xdmf = XDMFFile(domain.comm, "terzaghi.xdmf", "w")
#xdmf.write_mesh(domain)

with io.XDMFFile(MPI.COMM_WORLD, "ux.xdmf", "w") as xdmf1:
xdmf1.write_mesh(domain)
with io.XDMFFile(MPI.COMM_WORLD, "uy.xdmf", "w") as xdmf2:
xdmf2.write_mesh(domain)
with io.XDMFFile(MPI.COMM_WORLD, "p.xdmf", "w") as xdmf3:
xdmf3.write_mesh(domain)

t = 0
u, p = xt.split()
u.name = "u"
p.name = "p"

dx = u.sub(0).collapse()
dy = u.sub(1).collapse()

xdmf1.write_function(dx, t)
xdmf2.write_function(dy, t)
xdmf3.write_function(p, t)

for (i, dt) in enumerate(np.diff(time)):

t = time[i+1]
print('-----------------------------------------')
print('>> t =', t, '[sec]')

# solve system
solver.solve(xt)
#xt.x.scatter_forward()

# update internal variable
x_n.vector.array[:] = xt.vector.array[:]
#u_n, p_n = split(x_n)

# data for visualization
u, p = xt.split()
#u.name = "u"
#p.name = "p"

dx = u.sub(0).collapse()
dy = u.sub(1).collapse()

#poro = project((1.+ epsilon_v(u_n1))*poro0, W) # update porosity

xdmf1.write_function(dx, t)
xdmf2.write_function(dy, t)
xdmf3.write_function(p, t)

#xdmf_file.write(poro, t)

xdmf1.close()
xdmf2.close()
xdmf3.close()
``````

I could get the result. HOWEVER. with time loop, I could get the solution for 1st step (just after initial condition), but from 2nd step, the u, p are no longer updated so my solution (u, p) is always the same through time loop.

I think you are quite familiar with such problems, so could you help me to find out what is wrong with my code? I strongly think that the last 2 paragraphs would have some errors

Thanks a lot.

1 Like