Mixing dolfinx_mpc and scifem packages for a nonlinear problem

Ergonsa · July 5, 2025, 12:51am

Hello everyone,
the goal of this post is to ask whether it is possible to use all the different packages together in the same problem. For instance, I am working in a very ambitious project regarding homogenisation of magneto-viscohyperelastic materials with extreme deformation patterns. In order to implement such a problem in dolfinx, I would need MPC for imposing periodicity, lagrange multipliers per element (for incompressibility conditions) and per mesh (global multiplers for imposing integral constraints of the total stress) and nonlinear problem/solver classes. I have been reading through all the discourse group and I have seen that everything has been more or less implemented in different auxiiliary packages (dolfinx_mpc, scifem, multiphenics or github code). However, I have doubts regarding the compatibility of the global lagrange multiplier with dolfinx_mpc (in fenics legacy, it was done with R elements, but now,as far as I know, one would have to consider a nullspace and assemble the matrix directly for the linear problem and develop the entire nonlinear iteration scheme). With all this:

(1) Is it possible to implement such a problem in dolfinx by combining all the auxiliary packages or are there incompatibilities between them?

(2) If possible, how could “dolfinx_mpc” be combined with the real spaces of “scifem” to consider the global lagrange multipliers? Is there any example already considering this?

Thank you so much in advance. Any suggestion would be appreciated.

Stein · July 7, 2025, 9:49am

It sounds to me like all of what you aim to do is possible. A Real space doesn’t require knowledge on periodicity, so I wouldn’t expect problems there…
In general, I’d advice to simply build up complexity, and test to see if you can get the packages to work together on a very simple PDE first (e.g., implement a global Lagrane multiplier for a 2D Poisson equation on periodic grid). When you run into problems creating such minimal working examples people can give you more targeted advice.

Ergonsa · July 8, 2025, 6:20pm

Thank you for your suggestion @Stein. I have been lately trying to provide an example but I was not capable of making it work. At first I thought of a 2D Poisson equation, but that problem is not complex enough for imposing such BCs. Therefore, I have moved to the Navier-Stokes equations, following the implementation in: Test problem 1: Channel flow (Poiseuille flow). I wanted to try to substitute the top and bottom wall condition (u=0 at walls) by a vertical periodic condition with just u=0 at y=0 (which should yield the same results). The problem is that I really don’t know when the function “assemble_matrix” should be the one of dolfinx or the one of dolfinx_mpc in this case. This is my code (note that I have introduced a flag at the beginning of the code, MPC_FLAG, which switches between MPC active or not):

from mpi4py import MPI
from petsc4py import PETSc
import numpy as np
import pyvista

import dolfinx
from dolfinx.fem import Constant, Function, functionspace, assemble_scalar, dirichletbc, form, locate_dofs_geometrical
from dolfinx.fem.petsc import assemble_matrix, assemble_vector, apply_lifting, create_vector, set_bc
from dolfinx.io import VTXWriter
from dolfinx.mesh import create_unit_square
from dolfinx.plot import vtk_mesh
from dolfinx.common import Timer
from basix.ufl import element
from ufl import (FacetNormal, Identity, TestFunction, TrialFunction,
                 div, dot, ds, dx, inner, lhs, nabla_grad, rhs, sym)

from dolfinx_mpc import MultiPointConstraint, LinearProblem
from dolfinx_mpc import apply_lifting as apply_lifting_mpc
from dolfinx_mpc import assemble_matrix as assemble_matrix_mpc
from dolfinx_mpc import assemble_vector as assemble_vector_mpc


MPC_FLAG = True

### CREATE MESH ##
mesh = create_unit_square(MPI.COMM_WORLD, 10, 10)
t = 0
T = 10
num_steps = 500
dt = T / num_steps
tol = 250 * np.finfo(dolfinx.default_scalar_type).resolution
###
### CREATE ELEMENTS & FUNCTIONSPACES
v_cg2 = element("Lagrange", mesh.topology.cell_name(), 2, shape=(mesh.geometry.dim, ))
s_cg1 = element("Lagrange", mesh.topology.cell_name(), 1)
V = functionspace(mesh, v_cg2)
Q = functionspace(mesh, s_cg1)
###
### CREATE TEST & TRIAL FUNCTIONS
u = TrialFunction(V)
v = TestFunction(V)
p = TrialFunction(Q)
q = TestFunction(Q)
###
### CREATE BOUNDARY CONDITIONS
def walls(x):
    return np.logical_or(np.isclose(x[1], 0), np.isclose(x[1], 1))

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

if MPC_FLAG:
    wall_dofs = locate_dofs_geometrical(V, walls_mpc)
else:
    wall_dofs = locate_dofs_geometrical(V, walls)
u_noslip = np.array((0,) * mesh.geometry.dim, dtype=PETSc.ScalarType)
bc_noslip = dirichletbc(u_noslip, wall_dofs, V)

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

inflow_dofs = locate_dofs_geometrical(Q, inflow)
bc_inflow = dirichletbc(PETSc.ScalarType(8), inflow_dofs, Q)

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

outflow_dofs = locate_dofs_geometrical(Q, outflow)
bc_outflow = dirichletbc(PETSc.ScalarType(0), outflow_dofs, Q)
bcu = [bc_noslip]
bcp = [bc_inflow, bc_outflow]
###
### CREATE MPC FOR TOP WALL
def periodic_boundary(x):
    return np.isclose(x[1], 1, atol=tol)

def periodic_relation(x):
    out_x = np.zeros_like(x)
    out_x[0] = out_x[0]
    out_x[1] = 0.0
    out_x[2] = out_x[2]
    return out_x

with Timer("~PERIODIC: Initialize MPC"):
    mpc = MultiPointConstraint(V)
    mpc.create_periodic_constraint_geometrical(V, periodic_boundary, periodic_relation, bcu)
    mpc.finalize()
###
if MPC_FLAG:
    u_n = Function(mpc.function_space)
else:
    u_n = Function(V)
u_n.name = "u_n"
U = 0.5 * (u_n + u)
n = FacetNormal(mesh)
f = Constant(mesh, PETSc.ScalarType((0, 0)))
k = Constant(mesh, PETSc.ScalarType(dt))
mu = Constant(mesh, PETSc.ScalarType(1))
rho = Constant(mesh, PETSc.ScalarType(1))

# Define strain-rate tensor
def epsilon(u):
    return sym(nabla_grad(u))

# Define stress tensor
def sigma(u, p):
    return 2 * mu * epsilon(u) - p * Identity(len(u))

# Define the variational problem for the first step
p_n = Function(Q)
p_n.name = "p_n"
F1 = rho * dot((u - u_n) / k, v) * dx
F1 += rho * dot(dot(u_n, nabla_grad(u_n)), v) * dx
F1 += inner(sigma(U, p_n), epsilon(v)) * dx
F1 += dot(p_n * n, v) * ds - dot(mu * nabla_grad(U) * n, v) * ds
F1 -= dot(f, v) * dx
a1 = form(lhs(F1))
L1 = form(rhs(F1))
if MPC_FLAG:
    A1 = assemble_matrix_mpc(a1, constraint=mpc, bcs=bcu)
else:
    A1 = assemble_matrix(a1, bcs=bcu)
A1.assemble()
b1 = create_vector(L1)

# Define variational problem for step 2
if MPC_FLAG:
    u_ = Function(mpc.function_space)
else:
    u_ = Function(V)
a2 = form(dot(nabla_grad(p), nabla_grad(q)) * dx)
L2 = form(dot(nabla_grad(p_n), nabla_grad(q)) * dx - (rho / k) * div(u_) * q * dx)
A2 = assemble_matrix(a2, bcs=bcp)
A2.assemble()
b2 = create_vector(L2)

# Define variational problem for step 3
p_ = Function(Q)
u__ = Function(V)
a3 = form(rho * dot(u, v) * dx)
L3 = form(rho * dot(u__, v) * dx - k * dot(nabla_grad(p_ - p_n), v) * dx)
A3 = assemble_matrix(a3)
A3.assemble()
b3 = create_vector(L3)

# Solver for step 1
solver1 = PETSc.KSP().create(mesh.comm)
solver1.setOperators(A1)
solver1.setType(PETSc.KSP.Type.BCGS)
pc1 = solver1.getPC()
pc1.setType(PETSc.PC.Type.HYPRE)
pc1.setHYPREType("boomeramg")

# Solver for step 2
solver2 = PETSc.KSP().create(mesh.comm)
solver2.setOperators(A2)
solver2.setType(PETSc.KSP.Type.BCGS)
pc2 = solver2.getPC()
pc2.setType(PETSc.PC.Type.HYPRE)
pc2.setHYPREType("boomeramg")

# Solver for step 3
solver3 = PETSc.KSP().create(mesh.comm)
solver3.setOperators(A3)
solver3.setType(PETSc.KSP.Type.CG)
pc3 = solver3.getPC()
pc3.setType(PETSc.PC.Type.SOR)

from pathlib import Path
folder = Path("results")
folder.mkdir(exist_ok=True, parents=True)
vtx_u = VTXWriter(mesh.comm, folder / "poiseuille_u.bp", u_n, engine="BP4")
vtx_p = VTXWriter(mesh.comm, folder / "poiseuille_p.bp", p_n, engine="BP4")
vtx_u.write(t)
vtx_p.write(t)

def u_exact(x):
    values = np.zeros((2, x.shape[1]), dtype=PETSc.ScalarType)
    values[0] = 4 * x[1] * (1.0 - x[1])
    return values


u_ex = Function(V)
u_ex.interpolate(u_exact)

L2_error = form(dot(u_ - u_ex, u_ - u_ex) * dx)

for i in range(num_steps):
    # Update current time step
    t += dt

    # Step 1: Tentative veolcity step
    with b1.localForm() as loc_1:
        loc_1.set(0)

    if MPC_FLAG:
        assemble_vector_mpc(L1, mpc, b1)
        apply_lifting_mpc(b1, [a1], [bcu], mpc)
    else:
        assemble_vector(b1, L1)
        apply_lifting(b1, [a1], [bcu])
    b1.ghostUpdate(addv=PETSc.InsertMode.ADD_VALUES, mode=PETSc.ScatterMode.REVERSE)
    set_bc(b1, bcu)
    solver1.solve(b1, u_.x.petsc_vec)
    u_.x.scatter_forward()
    if MPC_FLAG:
        mpc.backsubstitution(u_)
    u__.x.array[:] = u_.x.array[:]

    # Step 2: Pressure corrrection step
    with b2.localForm() as loc_2:
        loc_2.set(0)
    assemble_vector(b2, L2)
    apply_lifting(b2, [a2], [bcp])
    b2.ghostUpdate(addv=PETSc.InsertMode.ADD_VALUES, mode=PETSc.ScatterMode.REVERSE)
    set_bc(b2, bcp)
    solver2.solve(b2, p_.x.petsc_vec)
    p_.x.scatter_forward()

    # Step 3: Velocity correction step
    with b3.localForm() as loc_3:
        loc_3.set(0)
    assemble_vector(b3, L3)
    b3.ghostUpdate(addv=PETSc.InsertMode.ADD_VALUES, mode=PETSc.ScatterMode.REVERSE)
    solver3.solve(b3, u__.x.petsc_vec)
    u__.x.scatter_forward()
    # Update variable with solution form this time step
    u_n.x.array[:] = u__.x.array[:]
    p_n.x.array[:] = p_.x.array[:]

    # Write solutions to file
    vtx_u.write(t)
    vtx_p.write(t)

    # Compute error at current time-step
    error_L2 = np.sqrt(mesh.comm.allreduce(assemble_scalar(L2_error), op=MPI.SUM))
    error_max = mesh.comm.allreduce(np.max(u_.x.petsc_vec.array - u_ex.x.petsc_vec.array), op=MPI.MAX)
    # Print error only every 20th step and at the last step
    if (i % 20 == 0) or (i == num_steps - 1):
        print(f"Time {t:.2f}, L2-error {error_L2:.2e}, Max error {error_max:.2e}")
# Close xmdf file
vtx_u.close()
vtx_p.close()
b1.destroy()
b2.destroy()
b3.destroy()
solver1.destroy()
solver2.destroy()
solver3.destroy()

#pyvista.start_xvfb()
topology, cell_types, geometry = vtk_mesh(V)
values = np.zeros((geometry.shape[0], 3), dtype=np.float64)
values[:, :len(u_n)] = u_n.x.array.real.reshape((geometry.shape[0], len(u_n)))

# Create a point cloud of glyphs
function_grid = pyvista.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
mesh.topology.create_connectivity(mesh.topology.dim, mesh.topology.dim)
grid = pyvista.UnstructuredGrid(*vtk_mesh(mesh, mesh.topology.dim))

# Create plotter
plotter = pyvista.Plotter()
plotter.add_mesh(grid, style="wireframe", color="k")
plotter.add_mesh(glyphs)
plotter.view_xy()
if not pyvista.OFF_SCREEN:
    plotter.show()
else:
    fig_as_array = plotter.screenshot("glyphs.png")

Topic		Replies	Views
Periodic Boundary Conditons for Mixed Functionspace with dolfinx_mpc General	10	1302	January 4, 2023
Unequal results for linear and nonlinear solver with dolfinx_mpc dolfinx	5	1039	February 1, 2023
Using dolfinx_mpc in generalised plain strain formulation dolfinx dolfinx	6	51	July 12, 2025
NonlinearMPCProblem: error with periodic structure and nonlinear material dolfinx	5	70	June 23, 2025
Explicitly assembling matrices with dolfinx_mpc General	31	1052	August 11, 2023

Mixing dolfinx_mpc and scifem packages for a nonlinear problem

Related topics