Boundary conditions on decoupled system

Meryem · December 31, 2025, 12:28pm

Hello,

I am trying to solve a decoupled system. I first solve for pressure and saturation (mixed space) (non linear system) then I solve for the displacement (linear system). I define separate dirichlet boundary conditions on each problem. Here is a minimal working example. For the first problem, the unknown is the increments of the pressure and the saturation. and for the second, the unknown is the displacement. The pressure and displacement results take huge values since the first iteration and when I visualise the results I notice that the boundary conditions are not respected in the external boundary during the simulation. I already tried to reproduce the mixed space bc in the fenicsX tutorial. Could you please tell me what I am doing wrong with the boundary condition definitions ? is it the update of the variables? Thank you.

from mpi4py import MPI
import basix
from dolfinx import mesh, fem, default_scalar_type, log
from dolfinx.fem.petsc import NonlinearProblem, LinearProblem
from dolfinx.nls.petsc import NewtonSolver
import petsc4py
import adios4dolfinx
from pathlib import Path
import time
from dolfinx.io import XDMFFile
import numpy as np
import ufl
import time


start = time.time()

##### FUNCTIONS #####

def deform(u):
    return ufl.sym(ufl.grad(u))
def lambda_s(E,nu):
        return E*nu/((1+nu)*(1-2*nu))
def mu(E,nu):
	return E/(2*(1+nu))

def centre_region(x):
    # x has shape (3, N)
    dx = x.T - center
    return np.linalg.norm(dx, axis=1) <= r

dx = ufl.Measure("dx", domain=domain)

##### FUNCTIONS #####

comm = MPI.COMM_WORLD
domain = mesh.create_unit_cube(comm, 16, 16, 16)

P1 = basix.ufl.element("Lagrange", domain.topology.cell_name(),  degree=1)
P2 = basix.ufl.element("Lagrange", domain.topology.cell_name(), degree=2, shape=(3,))
MS = basix.ufl.mixed_element([P1, P1])

MS_func         = fem.functionspace(domain, MS)
P1_func         = fem.functionspace(domain, P1)
P2_func         = fem.functionspace(domain, P2)


# subdomain centre and radius
center = np.array([0.5, 0.5, 0.5], dtype=np.float64)
r = 0.2



#----------------------------------------------------------------------
# Read meshtags 
metadata = {"quadrature_degree": 4}
dx       = ufl.Measure("dx", domain=domain, metadata=metadata)
print("metadata and dx ok")
#----------------------------------------------------------------------
#Locate exterior facets
fdim = domain.topology.dim - 1

domain.topology.create_connectivity(fdim, domain.topology.dim)
boundary_facets = mesh.exterior_facet_indices(domain.topology)

# #boundary conditions facets
facet_indices = boundary_facets
facet_markers = 3*np.ones_like(facet_indices, dtype=np.int32)  # mark all boundary facets as 3
facet_tags = mesh.meshtags(domain, fdim, facet_indices, facet_markers)
print("exterior facets located")

with XDMFFile(comm, "facet_tags_minimal_example.xdmf", "w") as xdmf:
    xdmf.write_mesh(domain)                 # Write mesh geometry and topology
    xdmf.write_meshtags(facet_tags,domain.geometry)  # visualise outside boundary, where the bcs are defined.

#----------------------------------------------------------------------

#Functions
#Fluid problem 
#Test and trial functions
(ql,qt) = ufl.TestFunctions(MS_func) 
ddw_fluid = ufl.TrialFunction(MS_func)
#Primary variables (solution increment)
dw_fluid = fem.Function(MS_func)
dpl, dSt= ufl.split(dw_fluid)

#previous solution (initial solution)
w0 = fem.Function(MS_func)
pl_n,St_n     = ufl.split(w0)

#initial conditions
pl_n_, pl_n_to_MS = MS_func.sub(0).collapse()
St_n_, St_n_to_MS = MS_func.sub(1).collapse()

# For pressure
Fpl_n_ = fem.Function(pl_n_)
Fpl_n_.x.array[:] = 460  # set all DOFs to zero
Fpl_n_.x.scatter_forward()


# For saturation
FSt_n_ = fem.Function(St_n_)
FSt_n_.x.array[:] = 0  # set all DOFs to zero
FSt_n_.x.scatter_forward()
dofs_centre_mixed = fem.locate_dofs_geometrical((St_n_, P1_func), centre_region)
FSt_n_.x.array[dofs_centre_mixed] = 0.5
FSt_n_.x.scatter_forward()


w0.x.array[pl_n_to_MS] = Fpl_n_.x.array
w0.x.array[St_n_to_MS] = FSt_n_.x.array
w0.x.scatter_forward()



print("Fluid Initial variables written")
#Create boundary conditions

#dpl=0
bcs_fluid       = []
PL, _ = MS_func.sub(0).collapse()
pl_func = fem.Function(PL)
pl_dofs   = fem.locate_dofs_topological((MS_func.sub(0), PL), fdim, boundary_facets)
bcs_fluid.append(fem.dirichletbc(pl_func, pl_dofs, MS_func.sub(0)))
print("bcs ux uy uz pl",len(bcs_fluid))
print(bcs_fluid)


print("pl bc ok")

#Mech problem
#Test and trial functions
v = ufl.TestFunction(P2_func)
u = ufl.TrialFunction(P2_func) 
#Primary variables
u_s  = fem.Function(P2_func) # solution for linear problem
#previous solution (initial solution)
u_n = fem.Function(P2_func) 

#U
u_s.x.array[:] = 0
u_n.x.array[:] = 0

u_n.x.scatter_forward()
u_s.x.scatter_forward()

#adios4dolfinx.read_function(filename, u_s, time = 1530, name="Solid displacement")
print("Mech Initial variables written")



u_bc = np.array((0,) * domain.geometry.dim, dtype=default_scalar_type)
bc_dofs = fem.locate_dofs_topological(P2_func, fdim, boundary_facets)
bcs_mech = [fem.dirichletbc(u_bc, bc_dofs, P2_func)]

#Also tested: 
# uy=0
#dofs   = fem.locate_dofs_topological(P2_func.sub(1), fdim, boundary_facets)
#bcs_mech.append(fem.dirichletbc(default_scalar_type(0), dofs, P2_func.sub(1)))
# ux=0
#dofs   = fem.locate_dofs_topological(P2_func.sub(0), fdim, boundary_facets)
#bcs_mech.append(fem.dirichletbc(default_scalar_type(0), dofs, P2_func.sub(0)))
# uz=0
#dofs   = fem.locate_dofs_topological(P2_func.sub(2), fdim, boundary_facets)
#bcs_mech.append(fem.dirichletbc(default_scalar_type(0), dofs, P2_func.sub(2)))

print(bcs_mech)

print("u bc ok")

#----------------------------------------------------------------------
#Constants
num_steps      = int(default_scalar_type(5))
dt             = fem.Constant(domain, default_scalar_type(0.1))
t              = default_scalar_type(0)
mu_l           = default_scalar_type(3.5e-3)
nu             = default_scalar_type(0.49)
a_             = default_scalar_type(500)
mu_t           = default_scalar_type(100)
alpha          = default_scalar_type(1e-4)
poro           = default_scalar_type(0.2)
k              = default_scalar_type(1e-15)
w_bs           = default_scalar_type(0.5)
E              = default_scalar_type(400)
ps             = default_scalar_type(200)

#----------------------------------------------------------------------
#----------------------------------------------------------------------
# Variational form


# Define form F 
F_fluid = (1/dt)*(1-(St_n+dSt))*ufl.nabla_div(u_s-u_n)*ql*dx-(1/dt)*poro*(dSt)*ql*dx+(0.551*k/(mu_l))*ufl.dot(ufl.grad(pl_n+dpl),ufl.grad(ql))*dx+ alpha *poro*(St_n+dSt)*(1-St_n-dSt)*ql*dx
#*(1-St_n-dSt)*#x

# add to mass conservation 
F_fluid += (1/dt)*(St_n+dSt)*ufl.nabla_div(u_s-u_n)*qt*dx+(1/dt)*poro*(dSt)*qt*dx+(0.551*k/(mu_t))*ufl.dot(ufl.grad(pl_n+dpl+a_*ufl.tan((np.pi/2)*(St_n+dSt))),ufl.grad(qt))*dx -alpha*poro*(St_n+dSt)*(1-St_n-dSt)*qt*dx


# Jacobian of the fluid problem
J_fluid = ufl.derivative(F_fluid, dw_fluid, ddw_fluid)


 #Momentum conservation
 #linear problem 
F_mech = (1/dt)*(ufl.inner(2*mu(E, nu)*deform(u), deform(v))*dx
        + lambda_s(E, nu)*ufl.nabla_div(u)*ufl.nabla_div(v)*dx
        - ps * ufl.nabla_div(v) * dx)


F_mech += -(1/dt)*(ufl.inner(2*mu(E, nu)*deform(u_n), deform(v))*dx
          + lambda_s(E, nu)*ufl.nabla_div(u_n)*ufl.nabla_div(v)*dx
          - ps * ufl.nabla_div(v) * dx)


a = ufl.lhs(F_mech)
L = ufl.rhs(F_mech)

print("All defined. Now solving")


#Solver parameters 
problem_fluid = NonlinearProblem(F_fluid, dw_fluid, bcs = bcs_fluid, J = J_fluid)
problem_solid = LinearProblem(a, L, bcs=bcs_mech, petsc_options={"ksp_type": "preonly", "pc_type": "lu", "pc_factor_mat_solver_type": "mumps"})
print("problem created")

solver_fluid  = NewtonSolver(domain.comm, problem_fluid)
# Maximum iterations
solver_fluid.max_it                = 15
# Absolute tolerance
solver_fluid.atol                  = 1e-10
# Relative tolerance
solver_fluid.rtol                  = 1e-10
# Convergence criterion
solver_fluid.convergence_criterion = "incremental"
solver_fluid.report = True

print('problem_solid.solve() written and param ok')

# Solver Pre-requisites
ksp_fluid                                              = solver_fluid.krylov_solver
opts                                                   = petsc4py.PETSc.Options()
option_prefix1                                         = ksp_fluid.getOptionsPrefix()
print(option_prefix1)

opts[f"{option_prefix1}ksp_type"]                      = "preonly" #"gmres" 
opts[f"{option_prefix1}pc_type"]                       = "lu"      #"hypre" 
opts[f"{option_prefix1}pc_factor_mat_solver_type"]     = "mumps"   #"boomeramg" 
opts[f"mat_mumps_icntl_2"]                            = 2
opts[f"mat_mumps_icntl_4"]                            = 2
ksp_fluid.setFromOptions()



log_solve=True
if log_solve:
    log.set_log_level(log.LogLevel.INFO)



xdmf = XDMFFile(comm, "Decoupled_results_minimal_example_test_bc.xdmf", "w") 
xdmf.write_mesh(domain)
 

for N in range(num_steps):
    # Dynamic adjustment of `dt`
    if   N == 10:
        dt.value = 1

    ## Advance time
    t += dt.value
    if comm.rank == 0:
        print(f"Time: {t} sec   Iteration: {N}")
    try:
        dw_fluid.x.array[:] = 0.0
        dw_fluid.x.scatter_forward()
        num_its_fluid, converged_fluid = solver_fluid.solve(dw_fluid)
        assert (converged_fluid)
    except:
        if comm.rank == 0:
            print("*************") 
            print("Fluid subproblem solver failed")
            print("*************") 
            break    
    dw_fluid.x.scatter_forward() 

    # ---- Step 2: Solve mechanics subproblem: displacement ----
    
    u_s = problem_solid.solve()
    u_s.x.scatter_forward()
    #update previous displacement 
    u_n.x.array[:] = u_s.x.array
    u_n.x.scatter_forward()

    # update previous pressure and saturation
    w0.x.array[:] += dw_fluid.x.array 
    w0.x.scatter_forward()

    # ---- Post-processing: visualization and output ----

    u_vis = fem.Function(P1_vec_space)
    u_vis.interpolate(u_s)
    u_vis.name = "Solid displacement"
    __dpl, __dSt = dw_fluid.split()
    __dpl.name = "pressure inc"
    __dSt.name = "Saturation inc"
    xdmf.write_function(__dSt,t)
    xdmf.write_function(__dpl,t)
    __pl_n, __St_n = w0.split()
    __pl_n.name = "pressure "
    __St_n.name = "Saturation "
    xdmf.write_function(__pl_n,t)
    xdmf.write_function(__St_n,t)
    xdmf.write_function(u_vis,t)

    if comm.rank == 0:
        print(f"Time step {t}, Fluid its {num_its_fluid}")


xdmf.close()

end = time.time()

dokken · January 5, 2026, 2:42pm

Dear @Meryem,

This code is far away from minimal, which means that it will take quite alot of time for anyone to debug this.
For the sake of clarity, could you:
Try to minimize the problem to:

Only defining the dirichlet conditions for the two problems, not including the variational forms.
Use dolfinx.fem.petsc.set_bc, to apply the boundary condition to an existing function, and store the result, see for example: DirichletBC with mesh tags in dolfinx - weird behavior - #2 by dokken for a use-case of this.

Meryem · January 6, 2026, 9:26am

Dear @dokken ,

Thank you for you answer. As I am using fenicsX v0.9.0, I followed this tutorial: Mixed finite element problems — FEniCS Workshop here is an example of both boundary conditions defined:

from mpi4py import MPI
import basix
from dolfinx import mesh, fem, default_scalar_type, log
from dolfinx.fem.petsc import NonlinearProblem, LinearProblem
from dolfinx.nls.petsc import NewtonSolver
import petsc4py
import adios4dolfinx
from pathlib import Path
import time
from dolfinx.io import XDMFFile
import numpy as np
import ufl
import time


start = time.time()


##### FUNCTIONS #####

comm = MPI.COMM_WORLD
domain = mesh.create_unit_cube(comm, 16, 16, 16)

P1 = basix.ufl.element("Lagrange", domain.topology.cell_name(),  degree=1)
P2 = basix.ufl.element("Lagrange", domain.topology.cell_name(), degree=2, shape=(3,))
MS = basix.ufl.mixed_element([P1, P1])

MS_func         = fem.functionspace(domain, MS)
P1_func         = fem.functionspace(domain, P1)
P2_func         = fem.functionspace(domain, P2)



dx = ufl.Measure("dx", domain=domain)

#----------------------------------------------------------------------
# Read meshtags 
metadata = {"quadrature_degree": 4}
dx       = ufl.Measure("dx", domain=domain, metadata=metadata)
print("metadata and dx ok")
#----------------------------------------------------------------------
#Locate exterior facets
fdim = domain.topology.dim - 1

domain.topology.create_connectivity(fdim, domain.topology.dim)
boundary_facets = mesh.exterior_facet_indices(domain.topology)

# #boundary conditions facets
facet_indices = boundary_facets
facet_markers = 3*np.ones_like(facet_indices, dtype=np.int32)  # mark all boundary facets as 3
facet_tags = mesh.meshtags(domain, fdim, facet_indices, facet_markers)
print("exterior facets located")
#----------------------------------------------------------------------

#Functions
#Fluid problem 
#Test and trial functions
(ql,qt) = ufl.TestFunctions(MS_func) 
ddw_fluid = ufl.TrialFunction(MS_func)
#Primary variables (solution increment)
dw_fluid = fem.Function(MS_func)
dpl, dSt= ufl.split(dw_fluid)


#Create boundary conditions

#dpl=0

PL = MS_func.sub(0)
P, P_to_PL = PL.collapse()
pl_func = fem.Function(P)
pl_func.x.array[:]=0

combined_dofs = fem.locate_dofs_topological((PL, P), fdim, boundary_facets)

bcs_fluid = fem.dirichletbc(pl_func, combined_dofs, PL)
new_dw = fem.Function(MS_func)
bcs_fluid.set(new_dw.x.array)
new_dw.x.scatter_forward()

xdmf = XDMFFile(comm, "function_with_bc.xdmf", "w") 
xdmf.write_mesh(domain)
xdmf.write_function(new_dw.sub(0))


#Mech problem
#Test and trial functions
v = ufl.TestFunction(P2_func)
u = ufl.TrialFunction(P2_func) 
#Primary variables
u_s  = fem.Function(P2_func) # solution for linear problem
#U
u_s.x.array[:] = 0  
u_s.x.scatter_forward()



u_D = fem.Function(P2_func)
u_D.x.array[:] = 0

uh = fem.Function(P2_func)
dofs = fem.locate_dofs_topological(P2_func, fdim, boundary_facets)
bcs_mech = fem.dirichletbc(u_D, dofs)
bcs_mech.set(uh.x.array)
uh.x.scatter_forward()

print("Mech boundary conditions written")

Then I call them in the solver this way:

problem_fluid = NonlinearProblem(F_fluid, dw_fluid, bcs = [bcs_fluid], J = J_fluid)
problem_solid = LinearProblem(a, L, bcs=[bcs_mech], petsc_options={"ksp_type": "preonly", "pc_type": "lu", "pc_factor_mat_solver_type": "mumps"})

When I solve my system using these boundary conditions and print values of displacement at the boundaries, I still observe that the 0 bc is not respected.

Thank you for your help !

Meryem

dokken · January 6, 2026, 2:55pm

So your question is really just about the mechanics problem , where you are applying:

Meryem:

u_bc = np.array((0,) * domain.geometry.dim, dtype=default_scalar_type)
bc_dofs = fem.locate_dofs_topological(P2_func, fdim, boundary_facets)
bcs_mech = [fem.dirichletbc(u_bc, bc_dofs, P2_func)]

To me this looks all fine.

Could you try removing the fluid solver from the problem, as:

Meryem:

F_mech = (1/dt)*(ufl.inner(2*mu(E, nu)*deform(u), deform(v))*dx
        + lambda_s(E, nu)*ufl.nabla_div(u)*ufl.nabla_div(v)*dx
        - ps * ufl.nabla_div(v) * dx)


F_mech += -(1/dt)*(ufl.inner(2*mu(E, nu)*deform(u_n), deform(v))*dx
          + lambda_s(E, nu)*ufl.nabla_div(u_n)*ufl.nabla_div(v)*dx
          - ps * ufl.nabla_div(v) * dx)

to me seems well defined as long as you set a known u_n ( i see no dependence on the fluid problem in here, please correct me if Im wrong).

Meryem · January 9, 2026, 8:58am

That was a mistake from my side, ps is a function of pl_n and St_n and u_n is set as 0.

When I solve both problems, I print values of u_s at the boundaries this way, and I also try to force them to take 0 as a value:

    dofs_u_all = fem.locate_dofs_topological(P2_func, fdim, boundary_facets)
    vals_u_bnd = u_s.x.array[dofs_u_all]
    if comm.rank == 0:
        print(f"[Step {N}] u_s boundary min/max: {vals_u_bnd.min()}, {vals_u_bnd.max()}")

    u_s.x.array[dofs_u_all]=0
    vals_u_bnd = u_s.x.array[dofs_u_all]
    if comm.rank == 0:
        print(f"[Step {N}] u_s boundary min/max: {vals_u_bnd.min()}, {vals_u_bnd.max()}")


    u_norm = np.linalg.norm(u_s.x.array)
    if comm.rank == 0:
        print(f"[Step {N}] max|u| approx={np.abs(u_s.x.array).max()}, ||u||={u_norm}")

These are the values I get already at 3 iterations

Time: 0.1 sec   Iteration: 0
[Step 0] u_s boundary min/max: -2.502683575141684e-07, 2.4070067549498826e-07
[Step 0] u_s boundary min/max: 0.0, 0.0
[Step 0] max|u| approx=2.502683575159242e-07, ||u||=1.4729756787369482e-05
Time step 0.1, Fluid its 2
Time: 0.2 sec   Iteration: 1
[Step 1] u_s boundary min/max: -22.218458537230884, 22.09838081525254
[Step 1] u_s boundary min/max: 0.0, 0.0
[Step 1] max|u| approx=22.430651772795912, ||u||=2177.283650518729
Time step 0.2, Fluid its 2
Time: 0.30000000000000004 sec   Iteration: 2
[Step 2] u_s boundary min/max: -3019219683.750289, 2966692677.8965783
[Step 2] u_s boundary min/max: 0.0, 0.0
[Step 2] max|u| approx=3025991809.6255302, ||u||=354391067006.7157
Time step 0.30000000000000004, Fluid its 2

Topic		Replies	Views
Please help me to check whether boundary condition are right set dolfinx	13	1818	August 28, 2022
Mixed dimensional Diffusion Problem General dolfinx	9	150	September 19, 2025
Boundary conditions (all direction or only x) General	8	332	February 20, 2024
Struggling with vector values bcdirichlet dolfinx boundary-conditions	21	306	August 28, 2024
Boundary condition to single direction in the subspace of Mixed function space dolfinx boundary-conditions	9	446	January 3, 2024

Boundary conditions on decoupled system

Related topics