Solving coupling equation with MixedElement form

[psj1866]

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], 0) #x=0

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

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

def top(x):
    return np.isclose(x[1], 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]

# Define the loading as an expression depending on t (?)
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[304], 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.

Thanks for reading.

[driley]

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

[psj1866]

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], 0) #x=0

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

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

def top(x):
    return np.isclose(x[1], 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]

# Define the loading as an expression depending on t (?)
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[1],), 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.

[Mamadou]

Hi @psj1866.

Have you fixed your code?

Because I am trying to approximate the coupling equation with the MixedElement form.

Does the following link https://github.com/Matt-L-McLean/Computational-Geomechanics work, dear @driley?

Regards.

[driley]

Hi @Mamadou,

It seems they changed the name of repository. It is now at GitHub - Matt-L-McLean/poromechanics

Cheers

[Mamadou]

Thank you @driley for the information.