Periodic Boundary Conditons for Mixed Functionspace with dolfinx_mpc

[tuderic]

Hello everybody,

I was wondering if there is an opportunity to implement periodic boundary conditions for every sub space in a mixed function space.

When I am trying to do this, I got the following error message: RuntimeError: Periodic conditions for vector valued spaces are not implemented. However, I have seen that there was already an issue , so it seems to me that it should work.

Here is my minimal example:

# Load required libraries
#########################
from dolfinx      import fem, mesh
from dolfinx.fem  import FunctionSpace, Function
from dolfinx.mesh import locate_entities_boundary, meshtags
from dolfinx_mpc import MultiPointConstraint
import numpy as np
from mpi4py import MPI  
from ufl import TrialFunction, TestFunction, FiniteElement, MixedElement, split
from typing import List


##############################################################################
#                                     Define mesh                            #
##############################################################################

domain = mesh.create_unit_square(MPI.COMM_WORLD, 8, 8, mesh.CellType.quadrilateral)
dim = domain.topology.dim                                                       # Dimension


##############################################################################
#               Definition of elements, functions and spaces                 #
##############################################################################

# Elements
P1 = FiniteElement("CG", domain.ufl_cell(), 1)                               
mel = MixedElement([P1, P1, P1])
# Define FunctionSpace for coupled mixed space
V = FunctionSpace(domain, mel)

# Define TestFunctions
(dv_1, dv_2, dv_3) = TestFunction(V)                                            # test functions
# Define Function for solution and split into subspaces
q_sol = Function(V)                                                             # solution vector
(v_1, v_2, v_3) = split(q_sol)                                          # split solution vector in primary variables

# Define TrialFunctions
dq_sol = TrialFunction(V)                                                       # for the Jacobian (nonlinear problems)


##############################################################################
#                            periodic boundary conditions                    #
##############################################################################

def dirichlet_and_periodic_bcs(domain, functionspace, boundary_condition: List[str] = ["dirichlet", "periodic"], dbc_value = 0):
    """
    Function to set either dirichlet or periodic boundary conditions
    ----------
    boundary_condition
        First item describes b.c. on {x=0} and {x=1}
        Second item describes b.c. on {y=0} and {y=1}
    """
    
    fdim = domain.topology.dim - 1

    bcs = []
    pbc_directions = []
    pbc_slave_tags = []
    pbc_is_slave = []
    pbc_is_master = []
    pbc_meshtags = []
    pbc_slave_to_master_maps = []

    def generate_pbc_slave_to_master_map(i):
        def pbc_slave_to_master_map(x):
            out_x = x.copy() # was ist x.copy()?
            out_x[i] = x[i] - domain.geometry.x.max()
            return out_x
        return pbc_slave_to_master_map

    def generate_pbc_is_slave(i):
        return lambda x: np.isclose(x[i], domain.geometry.x.max())

    def generate_pbc_is_master(i):
        return lambda x: np.isclose(x[i], domain.geometry.x.min())

    # Parse boundary conditions
    for i, bc_type in enumerate(boundary_condition):
        
        if bc_type == "dirichlet":
            u_bc = fem.Function(functionspace)
            u_bc.x.array[:] = dbc_value # value of dirichlet bc needs to be passed into this function!

            def dirichletboundary(x):
                return np.logical_or(np.isclose(x[i], domain.geometry.x.min()), np.isclose(x[i], domain.geometry.x.max()))
            facets = locate_entities_boundary(domain, fdim, dirichletboundary)
            topological_dofs = fem.locate_dofs_topological(functionspace, fdim, facets)
            bcs.append(fem.dirichletbc(u_bc, topological_dofs))
        
        elif bc_type == "periodic":
            pbc_directions.append(i)
            pbc_slave_tags.append(i + 2)
            pbc_is_slave.append(generate_pbc_is_slave(i))
            pbc_is_master.append(generate_pbc_is_master(i))
            pbc_slave_to_master_maps.append(generate_pbc_slave_to_master_map(i))

            facets = locate_entities_boundary(domain, fdim, pbc_is_slave[-1])
            arg_sort = np.argsort(facets)
            pbc_meshtags.append(meshtags(domain,
                                          fdim,
                                          facets[arg_sort],
                                          np.full(len(facets), pbc_slave_tags[-1], dtype=np.int32)))

    # Create MultiPointConstraint object
    mpc = MultiPointConstraint(functionspace)
    
    N_pbc = len(pbc_directions)
    for i in range(N_pbc):
        if N_pbc > 1:
            def pbc_slave_to_master_map(x):
                out_x = pbc_slave_to_master_maps[i](x)
                idx = pbc_is_slave[(i + 1) % N_pbc](x)
                out_x[pbc_directions[i]][idx] = np.nan
                return out_x
        else:
            pbc_slave_to_master_map = pbc_slave_to_master_maps[i]

        mpc.create_periodic_constraint_topological(functionspace, pbc_meshtags[i],
                                                    pbc_slave_tags[i],
                                                    pbc_slave_to_master_map,
                                                    bcs)
        
    if len(pbc_directions) > 1:
        # Map intersection(slaves_x, slaves_y) to intersection(masters_x, masters_y),
        # i.e. map the slave dof at (1, 1) to the master dof at (0, 0)
        def pbc_slave_to_master_map(x):
            out_x = x.copy()
            out_x[0] = x[0] - domain.geometry.x.max()
            out_x[1] = x[1] - domain.geometry.x.max()
            idx = np.logical_and(pbc_is_slave[0](x), pbc_is_slave[1](x))
            out_x[0][~idx] = np.nan
            out_x[1][~idx] = np.nan
            return out_x
        mpc.create_periodic_constraint_topological(functionspace, pbc_meshtags[1],
                                                    pbc_slave_tags[1],
                                                    pbc_slave_to_master_map,
                                                    bcs)
        
    mpc.finalize()
    
    return mpc, bcs


mpc, bcs = dirichlet_and_periodic_bcs(domain, V, ["periodic", "periodic"]) 

Please note that I have only made small changes the function for the boundary conditions from here .

Thank you for your help.

[dokken]

You should constrain each sub element of your mixed element, as for instance done for the 0th component in: https://github.com/jorgensd/dolfinx_mpc/blob/main/python/demos/demo_periodic3d_topological.py#L76-L79
i.e.

@tuderic

     mpc.create_periodic_constraint_topological(functionspace, pbc_meshtags[1],
                                                    pbc_slave_tags[1],
                                                    pbc_slave_to_master_map,
                                                    bcs)

should be:

for i in range(functionspace.num_sub_spaces):
     mpc.create_periodic_constraint_topological(functionspace.sub(i), pbc_meshtags[1],
                                                pbc_slave_tags[1],
                                                pbc_slave_to_master_map,
                                                bcs)

[tuderic]

Thank you very much for your reply!

I figured it out by myself, but I wanted to verify the solution before posting it to this entry. However, an error occurs applying the class NonlinearMPCProblem with TypeError: NonlinearProblem.__init__() got an unexpected keyword argument 'form_compiler_options'.

Does anyone have an idea what I am doing wrong?

Here is my MWE (the mesh files can be found here):

# Load required libraries
#########################
import dolfinx
from dolfinx      import fem
from dolfinx.fem  import FunctionSpace, Function, locate_dofs_geometrical, dirichletbc, Constant
from dolfinx.mesh import locate_entities_boundary, meshtags
from dolfinx.io   import XDMFFile
import dolfinx_mpc

import numpy as np
from mpi4py import MPI  
import ufl
from typing import List
import pytest

from petsc4py import PETSc

##############################################################################
#              Classes for nonlinear MPC problems                            #
##############################################################################

class NonlinearMPCProblem(dolfinx.fem.petsc.NonlinearProblem):

    def __init__(self, F, u, mpc, bcs=[], J=None, form_compiler_options={},
                  jit_options={}):
        self.mpc = mpc
        super().__init__(F, u, bcs=bcs, J=J,
                          form_compiler_options=form_compiler_options,
                          jit_options=jit_options)
        
    def F(self, x: PETSc.Vec, F: PETSc.Vec):
        with F.localForm() as F_local:
            F_local.set(0.0)
        dolfinx_mpc.assemble_vector(self._L, self.mpc, b=F)

        # Apply boundary condition
        dolfinx_mpc.apply_lifting(F, [self._a], bcs=[self.bcs],
                                  constraint=self.mpc, x0=[x], scale=-1.0)
        F.ghostUpdate(addv=PETSc.InsertMode.ADD, mode=PETSc.ScatterMode.REVERSE)
        dolfinx.fem.petsc.set_bc(F, self.bcs, x, -1.0)

    def J(self, x: PETSc.Vec, A: PETSc.Mat):
        A.zeroEntries()
        dolfinx_mpc.assemble_matrix(self._a, self.mpc, bcs=self.bcs, A=A)
        A.assemble()


class NewtonSolverMPC(dolfinx.cpp.nls.petsc.NewtonSolver):
    def __init__(self, comm: MPI.Intracomm, problem: NonlinearMPCProblem,
                 mpc: dolfinx_mpc.MultiPointConstraint):
        """A Newton solver for non-linear MPC problems."""
        super().__init__(comm)
        self.mpc = mpc
        self.u_mpc = dolfinx.fem.Function(mpc.function_space)

        # Create matrix and vector to be used for assembly of the non-linear
        # MPC problem
        self._A = dolfinx_mpc.cpp.mpc.create_matrix(
            problem.a, mpc._cpp_object)
        self._b = dolfinx.cpp.la.petsc.create_vector(
            mpc.function_space.dofmap.index_map,
            mpc.function_space.dofmap.index_map_bs)

        self.setF(problem.F, self._b)
        self.setJ(problem.J, self._A)
        self.set_form(problem.form)
        self.set_update(self.update)

    def update(self, solver: dolfinx.nls.petsc.NewtonSolver,
               dx: PETSc.Vec, x: PETSc.Vec):
        # We need to use a vector created on the MPC's space to update ghosts
        self.u_mpc.vector.array = x.array_r
        self.u_mpc.vector.axpy(-1.0, dx)
        self.u_mpc.vector.ghostUpdate(addv=PETSc.InsertMode.INSERT,
                                      mode=PETSc.ScatterMode.FORWARD)
        self.mpc.homogenize(self.u_mpc.vector)
        self.mpc.backsubstitution(self.u_mpc.vector)
        x.array = self.u_mpc.vector.array_r
        x.ghostUpdate(addv=PETSc.InsertMode.INSERT,
                      mode=PETSc.ScatterMode.FORWARD)

    def solve(self, u: dolfinx.fem.Function):
        """Solve non-linear problem into function u. Returns the number
        of iterations and if the solver converged."""
        n, converged = super().solve(u.vector)
        u.x.scatter_forward()
        return n, converged

    @property
    def A(self) -> PETSc.Mat:
        """Jacobian matrix"""
        return self._A

    @property
    def b(self) -> PETSc.Vec:
        """Residual vector"""
        return self._b


@pytest.mark.skipif(np.issubdtype(PETSc.ScalarType, np.complexfloating),
                    reason="This test does not work in complex mode.")
@pytest.mark.parametrize("poly_order", [1, 2, 3])
def test_nonlinear_possion(poly_order):
    # Solve a standard Poisson problem with known solution which has
    # rotational symmetry of pi/2 at (x, y) = (0.5, 0.5). Therefore we may
    # impose MPCs on those DoFs which lie on the symmetry plane(s) and test
    # our numerical approximation. We do not impose any constraints at the
    # rotationally degenerate point (x, y) = (0.5, 0.5).

    N_vals = np.array([4, 8, 16], dtype=np.int32)
    l2_error = np.zeros_like(N_vals, dtype=np.double)
    for run_no, N in enumerate(N_vals):
        mesh = dolfinx.mesh.create_unit_square(MPI.COMM_WORLD, N, N)

        V = dolfinx.fem.FunctionSpace(mesh, ("Lagrange", poly_order))

        def boundary(x):
            return np.ones_like(x[0], dtype=np.int8)

        u_bc = dolfinx.fem.Function(V)
        with u_bc.vector.localForm() as u_local:
            u_local.set(0.0)

        facets = dolfinx.mesh.locate_entities_boundary(mesh, 1, boundary)
        topological_dofs = dolfinx.fem.locate_dofs_topological(V, 1, facets)
        zero = np.array(0, dtype=PETSc.ScalarType)
        bc = dolfinx.fem.dirichletbc(zero, topological_dofs, V)
        bcs = [bc]

        # Define variational problem
        u = dolfinx.fem.Function(V)
        v = ufl.TestFunction(V)
        x = ufl.SpatialCoordinate(mesh)

        u_soln = ufl.sin(ufl.pi * x[0]) * ufl.sin(ufl.pi * x[1])
        f = -ufl.div((1 + u_soln**2) * ufl.grad(u_soln))

        F = ufl.inner((1 + u**2) * ufl.grad(u), ufl.grad(v)) * ufl.dx \
            - ufl.inner(f, v) * ufl.dx
        J = ufl.derivative(F, u)

        # -- Impose the pi/2 rotational symmetry of the solution as a constraint,
        # -- except at the centre DoF
        def periodic_boundary(x):
            eps = 1e-10
            return np.isclose(x[0], 0.5) & (
                (x[1] < 0.5 - eps) | (x[1] > 0.5 + eps))

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

        mpc = dolfinx_mpc.MultiPointConstraint(V)
        mpc.create_periodic_constraint_geometrical(
            V, periodic_boundary, periodic_relation, bcs)
        mpc.finalize()

        # Sanity check that the MPC class has some constraints to impose
        num_slaves_global = mesh.comm.allreduce(len(mpc.slaves), op=MPI.SUM)
        num_masters_global = mesh.comm.allreduce(
            len(mpc.masters.array), op=MPI.SUM)

        assert num_slaves_global > 0
        assert num_masters_global == num_slaves_global

        problem = NonlinearMPCProblem(F, u, mpc, bcs=bcs, J=J)
        solver = NewtonSolverMPC(mesh.comm, problem, mpc)

        # Ensure the solver works with nonzero initial guess
        u.interpolate(lambda x: x[0]**2 * x[1]**2)
        solver.solve(u)

        l2_error_local = dolfinx.fem.assemble_scalar(
            dolfinx.fem.form((u - u_soln) ** 2 * ufl.dx))
        l2_error_global = mesh.comm.allreduce(l2_error_local, op=MPI.SUM)

        l2_error[run_no] = l2_error_global**0.5

    rates = np.log(l2_error[:-1] / l2_error[1:]) / np.log(2.0)

    assert np.all(rates > poly_order + 0.9)


@pytest.mark.parametrize("element", [ufl.FiniteElement, ufl.VectorElement])
@pytest.mark.parametrize("poly_order", [1, 2, 3])
def test_homogenize(element, poly_order):
    mesh = dolfinx.mesh.create_unit_square(MPI.COMM_WORLD, 8, 8)

    V = dolfinx.fem.FunctionSpace(
        mesh, element("CG", mesh.ufl_cell(), poly_order))

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

    def periodic_relation(x):
        out_x = np.zeros(x.shape)
        out_x[0] = 1.0 - x[0]
        out_x[1] = x[1]
        out_x[2] = x[2]
        return out_x

    mpc = dolfinx_mpc.MultiPointConstraint(V)
    mpc.create_periodic_constraint_geometrical(
        V, periodic_boundary, periodic_relation, [])
    mpc.finalize()

    # Sanity check that the MPC class has some constraints to impose
    num_slaves_global = mesh.comm.allreduce(len(mpc.slaves), op=MPI.SUM)
    assert num_slaves_global > 0

    u = dolfinx.fem.Function(V)
    u.vector.set(1.0)

    assert np.isclose(u.vector.min()[1], u.vector.max()[1])
    assert np.isclose(u.vector.array_r[0], 1.0)

    mpc.homogenize(u.vector)

    with u.vector.localForm() as u_:
        for i in range(V.dofmap.index_map.size_local * V.dofmap.index_map_bs):
            if i in mpc.slaves:
                assert np.isclose(u_.array_r[i], 0.0)
            else:
                assert np.isclose(u_.array_r[i], 1.0)


##############################################################################
#                                     Define mesh                            #
##############################################################################
comm = MPI.COMM_WORLD

meshFile = "2D_OneInclusion_v3"                                                 # mesh file name

# read meshfile and define the mesh(domain), meshtags of the domain (cell_tag) and of the facets (facet_tag)
with XDMFFile(comm, "mesh/"+meshFile+"_mesh.xdmf", "r") as xdmf:
    domain = xdmf.read_mesh(name="Grid")                                        # reads the mesh
    cell_tag = xdmf.read_meshtags(domain, name="Grid")                          # cell tags
domain.topology.create_connectivity(domain.topology.dim, domain.topology.dim-1) 

with XDMFFile(comm, "mesh/"+meshFile+"_facets.xdmf", "r") as xdmf:
    facet_tag = xdmf.read_meshtags(domain, name="Grid") 
dim = domain.topology.dim                                                       # Dimension


##  indices of the mesh
# number of cells
num_cells = domain.topology.index_map(domain.topology.dim).size_local
print("Number of cells",num_cells)
# number of vertices
print(f"Number of global vertices: {domain.topology.index_map(0).size_global}")
                                

# Normal vector
n = ufl.FacetNormal(domain)
# Different measures
dx   = ufl.Measure('dx', domain=domain, subdomain_data=cell_tag)                    # surfaces
ds   = ufl.Measure('ds', domain=domain, subdomain_data=facet_tag)                   # boundaries
dInt = ufl.Measure("dS", domain, subdomain_data=facet_tag)                          # interfaces
# spatial coordinates (can also be used in variational formulation)
x_coor = ufl.SpatialCoordinate(domain)


##############################################################################
#               Definition of elements, functions and spaces                 #
##############################################################################

# Elements
P1 = ufl.FiniteElement("CG", domain.ufl_cell(), 1)                               
mel = ufl.MixedElement([P1, P1])
# Define FunctionSpace for coupled mixed space
V = FunctionSpace(domain, mel)
# Define TestFunctions
(dv_1, dv_2) = ufl.TestFunction(V)                                                
# Define Function for solution and split into subspaces
v_sol = Function(V)	
(v_1, v_2) = ufl.split(v_sol)  
# Trial function for Jacobian
dv_sol = ufl.TrialFunction(V) 


##############################################################################
#                       Material Properties                                  #
##############################################################################

V_DG = FunctionSpace(domain, ("DG", 0))  

k_1 = Function(V_DG)
k_2 = Function(V_DG)

for i in range(num_cells):
    if cell_tag.values[i] == 1: # matrix
        k_1.x.array[i] = 1 # arbitrary material parameters
        k_2.x.array[i] = 3
    elif cell_tag.values[i] == 2: # inclusions
        k_1.x.array[i] = 5
        k_1.x.array[i] = 7
        
##############################################################################
#                            periodic boundary conditions                    #
##############################################################################



def dirichlet_and_periodic_bcs(domain, funcspace, boundary_condition: List[str] = ["dirichlet", "periodic"], dbc_value = 0):
    """
    Function to set either dirichlet or periodic boundary conditions
    ----------
    boundary_condition
        First item describes b.c. on {x=0} and {x=1}
        Second item describes b.c. on {y=0} and {y=1}
    """
    
    fdim = domain.topology.dim - 1

    bcs = []
    pbc_directions = []
    pbc_slave_tags = []
    pbc_is_slave = []
    pbc_is_master = []
    pbc_meshtags = []
    pbc_slave_to_master_maps = []
    
    # number of subspaces in functionspace
    N_subspaces = funcspace.num_sub_spaces
    
    # Create MultiPointConstraint object
    mpc = dolfinx_mpc.MultiPointConstraint(funcspace)

    def generate_pbc_slave_to_master_map(i):
        def pbc_slave_to_master_map(x):
            out_x = x.copy() 
            out_x[i] = x[i] - domain.geometry.x.max()
            return out_x
        return pbc_slave_to_master_map

    def generate_pbc_is_slave(i):
        return lambda x: np.isclose(x[i], domain.geometry.x.max())

    def generate_pbc_is_master(i):
        return lambda x: np.isclose(x[i], domain.geometry.x.min())

    # Parse boundary conditions
    def parse_bcs():
        for i, bc_type in enumerate(boundary_condition):
            
            if bc_type == "dirichlet":
                u_bc = fem.Function(functionspace)
                u_bc.x.array[:] = dbc_value # value of dirichlet bc needs to be passed into this function!
    
                def dirichletboundary(x):
                    return np.logical_or(np.isclose(x[i], domain.geometry.x.min()), np.isclose(x[i], domain.geometry.x.max()))
                facets = locate_entities_boundary(domain, fdim, dirichletboundary)
                topological_dofs = fem.locate_dofs_topological(functionspace, fdim, facets)
                bcs.append(fem.dirichletbc(u_bc, topological_dofs))
            
            elif bc_type == "periodic":
                pbc_directions.append(i)
                pbc_slave_tags.append(i + 2)
                pbc_is_slave.append(generate_pbc_is_slave(i))
                pbc_is_master.append(generate_pbc_is_master(i))
                pbc_slave_to_master_maps.append(generate_pbc_slave_to_master_map(i))
    
                facets = locate_entities_boundary(domain, fdim, pbc_is_slave[-1])
                arg_sort = np.argsort(facets)
                pbc_meshtags.append(meshtags(domain,
                                              fdim,
                                              facets[arg_sort],
                                              np.full(len(facets), pbc_slave_tags[-1], dtype=np.int32)))
        
        N_pbc = len(pbc_directions)
        for i in range(N_pbc):
            if N_pbc > 1:
                def pbc_slave_to_master_map(x):
                    out_x = pbc_slave_to_master_maps[i](x)
                    idx = pbc_is_slave[(i + 1) % N_pbc](x)
                    out_x[pbc_directions[i]][idx] = np.nan
                    return out_x
            else:
                pbc_slave_to_master_map = pbc_slave_to_master_maps[i]
    
            mpc.create_periodic_constraint_topological(functionspace, pbc_meshtags[i],
                                                        pbc_slave_tags[i],
                                                        pbc_slave_to_master_map,
                                                        bcs)
            
        if len(pbc_directions) > 1:
            # Map intersection(slaves_x, slaves_y) to intersection(masters_x, masters_y),
            # i.e. map the slave dof at (1, 1) to the master dof at (0, 0)
            def pbc_slave_to_master_map(x):
                out_x = x.copy()
                out_x[0] = x[0] - domain.geometry.x.max()
                out_x[1] = x[1] - domain.geometry.x.max()
                idx = np.logical_and(pbc_is_slave[0](x), pbc_is_slave[1](x))
                out_x[0][~idx] = np.nan
                out_x[1][~idx] = np.nan
                return out_x
            mpc.create_periodic_constraint_topological(functionspace, pbc_meshtags[1],
                                                        pbc_slave_tags[1],
                                                        pbc_slave_to_master_map,
                                                        bcs)

    if N_subspaces == 0:
        functionspace = funcspace
        parse_bcs()
    else:
        for j in range(N_subspaces):
            functionspace = funcspace.sub(j)
            parse_bcs()
            
    mpc.finalize()
    
    return mpc, bcs


mpc, bcs = dirichlet_and_periodic_bcs(domain, V, ["periodic", "periodic"]) 

### fix fluctuation at one point 

def boundary_node(x):
    return np.logical_and( np.isclose(x[0], domain.geometry.x.min()), np.isclose(x[1], domain.geometry.x.min()))

dofs_node_sub0 = locate_dofs_geometrical(V.sub(0).collapse()[0], boundary_node)
dofs_node_sub1 = locate_dofs_geometrical(V.sub(1).collapse()[0], boundary_node)
bc_node_sub0 = dirichletbc(PETSc.ScalarType(0.0), dofs_node_sub0, V.sub(0))
bc_node_sub1 = dirichletbc(PETSc.ScalarType(0.0), dofs_node_sub1, V.sub(1))

bcs.append(bc_node_sub0), bcs.append(bc_node_sub1)

### constant gradients for computational homogenization approach
gradient  = Constant(domain, PETSc.ScalarType((1,0)))
gradient2 = Constant(domain, PETSc.ScalarType((0,1)))

##############################################################################
#                               Variational problem                          #
##############################################################################
# weak formulation
eqn1 = k_1 * ufl.dot(ufl.grad(ufl.dot(gradient , x_coor) + v_1), ufl.grad(dv_1)) * dx 
eqn2 = k_2 * ufl.dot(ufl.grad(ufl.dot(gradient2, x_coor) + v_2), ufl.grad(dv_2)) * dx 

eqn = eqn1 + eqn2

# Jacobi determinant
Jac = ufl.derivative(eqn, v_sol, dv_sol)

##############################################################################
#                                      Solve                                 #
##############################################################################

problem = NonlinearMPCProblem(eqn, v_sol, mpc, bcs=bcs, J = Jac)                # define nonlinear problem
solver  = NewtonSolverMPC(comm, problem, mpc)                                   # define solver

# Set Solver Configurations
solver.convergence_criterion = "incremental"
solver.rtol = 1e-8
ksp = solver.krylov_solver
opts = PETSc.Options()
option_prefix = ksp.getOptionsPrefix()
opts[f"{option_prefix}ksp_type"] = "preonly"
opts[f"{option_prefix}pc_type"] = "lu"
ksp.setFromOptions()


# solve 
number, converged = solver.solve(v_sol)                                         # solve the set of PDEs
assert(converged)

[dokken]

What version of DOLFINx and DOLFINx_mpc are you running?
the *_options prefix was changed after the release of v0.5.2 of DOLFINx. Thus change *_options, as seen in: replace 'jit_paramters' with 'jit_option', and 'form_compiler_parameters' with 'form_compiler_options' by SarahRo · Pull Request #2310 · FEniCS/dolfinx · GitHub

[tuderic]

I am working with docker and the version dolfinx_mpc:v0.5.0. However, the error occurs as described.

I will try to figure out how to change the prefix in docker as you suggested.

[dokken]

@tuderic

I am working with docker and the version dolfinx_mpc:v0.5.0. However, the error occurs as described.

As I said, these changes was introduced into DOLFINx after the 0.5.0 version. You are probably using the source code from the main branch, which is compatible with the DOLFINx main branch.

[tuderic]

Thank you again for your help!

I just changed form_compiler_options to form_compiler_params and jit_options to jit_params in order to have the same syntax as in the older version.

However, the problem just does not converge even though it is a very simple problem. Is there any mistake in my code or does it not converge due to the classes for nonlinear MPC problems that seem to be in the testing phase?

Here is the full code (only slight changes from the code above):

# Load required libraries
#########################
import dolfinx
from dolfinx      import fem, log
from dolfinx.fem  import FunctionSpace, Function, locate_dofs_geometrical, dirichletbc, Constant
from dolfinx.mesh import locate_entities_boundary, meshtags
from dolfinx.io   import XDMFFile
import dolfinx_mpc

import numpy as np
from mpi4py import MPI  
import ufl
from typing import List
import pytest

from petsc4py import PETSc


##############################################################################
#              Classes for nonlinear MPC problems                            #
##############################################################################
# copied from https://github.com/jorgensd/dolfinx_mpc/blob/main/python/tests/test_nonlinear_assembly.py

class NonlinearMPCProblem(dolfinx.fem.petsc.NonlinearProblem):

    def __init__(self, F, u, mpc, bcs=[], J=None, form_compiler_params={},
                  jit_params={}):
        self.mpc = mpc
        super().__init__(F, u, bcs=bcs, J=J,
                          form_compiler_params=form_compiler_params,
                          jit_params=jit_params)
        
    def F(self, x: PETSc.Vec, F: PETSc.Vec):
        with F.localForm() as F_local:
            F_local.set(0.0)
        dolfinx_mpc.assemble_vector(self._L, self.mpc, b=F)

        # Apply boundary condition
        dolfinx_mpc.apply_lifting(F, [self._a], bcs=[self.bcs],
                                  constraint=self.mpc, x0=[x], scale=-1.0)
        F.ghostUpdate(addv=PETSc.InsertMode.ADD, mode=PETSc.ScatterMode.REVERSE)
        dolfinx.fem.petsc.set_bc(F, self.bcs, x, -1.0)

    def J(self, x: PETSc.Vec, A: PETSc.Mat):
        A.zeroEntries()
        dolfinx_mpc.assemble_matrix(self._a, self.mpc, bcs=self.bcs, A=A)
        A.assemble()


class NewtonSolverMPC(dolfinx.cpp.nls.petsc.NewtonSolver):
    def __init__(self, comm: MPI.Intracomm, problem: NonlinearMPCProblem,
                 mpc: dolfinx_mpc.MultiPointConstraint):
        """A Newton solver for non-linear MPC problems."""
        super().__init__(comm)
        self.mpc = mpc
        self.u_mpc = dolfinx.fem.Function(mpc.function_space)

        # Create matrix and vector to be used for assembly of the non-linear
        # MPC problem
        self._A = dolfinx_mpc.cpp.mpc.create_matrix(
            problem.a, mpc._cpp_object)
        self._b = dolfinx.cpp.la.petsc.create_vector(
            mpc.function_space.dofmap.index_map,
            mpc.function_space.dofmap.index_map_bs)

        self.setF(problem.F, self._b)
        self.setJ(problem.J, self._A)
        self.set_form(problem.form)
        self.set_update(self.update)

    def update(self, solver: dolfinx.nls.petsc.NewtonSolver,
               dx: PETSc.Vec, x: PETSc.Vec):
        # We need to use a vector created on the MPC's space to update ghosts
        self.u_mpc.vector.array = x.array_r
        self.u_mpc.vector.axpy(-1.0, dx)
        self.u_mpc.vector.ghostUpdate(addv=PETSc.InsertMode.INSERT,
                                      mode=PETSc.ScatterMode.FORWARD)
        self.mpc.homogenize(self.u_mpc.vector)
        self.mpc.backsubstitution(self.u_mpc.vector)
        x.array = self.u_mpc.vector.array_r
        x.ghostUpdate(addv=PETSc.InsertMode.INSERT,
                      mode=PETSc.ScatterMode.FORWARD)

    def solve(self, u: dolfinx.fem.Function):
        """Solve non-linear problem into function u. Returns the number
        of iterations and if the solver converged."""
        n, converged = super().solve(u.vector)
        u.x.scatter_forward()
        return n, converged

    @property
    def A(self) -> PETSc.Mat:
        """Jacobian matrix"""
        return self._A

    @property
    def b(self) -> PETSc.Vec:
        """Residual vector"""
        return self._b

##############################################################################
#                                     Define mesh                                                                                               #
##############################################################################
comm = MPI.COMM_WORLD

meshFile = "2D_OneInclusion_v3"                                                 # mesh file name

# read meshfile and define the mesh(domain), meshtags of the domain (cell_tag) and of the facets (facet_tag)
with XDMFFile(comm, "mesh/"+meshFile+"_mesh.xdmf", "r") as xdmf:
    domain = xdmf.read_mesh(name="Grid")                                        # reads the mesh
    cell_tag = xdmf.read_meshtags(domain, name="Grid")                          # cell tags
domain.topology.create_connectivity(domain.topology.dim, domain.topology.dim-1) 

with XDMFFile(comm, "mesh/"+meshFile+"_facets.xdmf", "r") as xdmf:
    facet_tag = xdmf.read_meshtags(domain, name="Grid") 
dim = domain.topology.dim                                                       # Dimension


##  indices of the mesh
# number of cells
num_cells = domain.topology.index_map(domain.topology.dim).size_local
print("Number of cells",num_cells)
# number of vertices
print(f"Number of global vertices: {domain.topology.index_map(0).size_global}")
                                

# Normal vector
n = ufl.FacetNormal(domain)
# Different measures
dx   = ufl.Measure('dx', domain=domain, subdomain_data=cell_tag)                    # surfaces
ds   = ufl.Measure('ds', domain=domain, subdomain_data=facet_tag)                   # boundaries
dInt = ufl.Measure("dS", domain, subdomain_data=facet_tag)                          # interfaces
# spatial coordinates (can also be used in variational formulation)
x_coor = ufl.SpatialCoordinate(domain)


##############################################################################
#               Definition of elements, functions and spaces                 #
##############################################################################

# Elements
P1 = ufl.FiniteElement("CG", domain.ufl_cell(), 1)                               
mel = ufl.MixedElement([P1, P1])
# Define FunctionSpace for coupled mixed space
V = FunctionSpace(domain, mel)

# Define TestFunctions
(dv_1, dv_2) = ufl.TestFunction(V)                                                
# # Define Function for solution and split into subspaces
v_sol = Function(V)	
(v_1, v_2) = ufl.split(v_sol)  
# # Trial function for Jacobian
dv_sol = ufl.TrialFunction(V) 


##############################################################################
#                       Material Properties                                  #
##############################################################################

V_DG = FunctionSpace(domain, ("DG", 0))  

k_1 = Function(V_DG)
k_2 = Function(V_DG)

for i in range(num_cells):
    if cell_tag.values[i] == 1: # matrix
        k_1.x.array[i] = 1 # arbitrary material parameters
        k_2.x.array[i] = 3
    elif cell_tag.values[i] == 2: # inclusions
        k_1.x.array[i] = 5
        k_1.x.array[i] = 7
        
##############################################################################
#                            periodic boundary conditions                    #
##############################################################################
# only slight changes compared to function "assemble_and_solve" in https://github.com/jorgensd/dolfinx_mpc/blob/main/python/demos/demo_periodic_gep.py

def dirichlet_and_periodic_bcs(domain, funcspace, boundary_condition: List[str] = ["dirichlet", "periodic"], dbc_value = 0):
    """
    Function to set either dirichlet or periodic boundary conditions
    ----------
    boundary_condition
        First item describes b.c. on {x=0} and {x=1}
        Second item describes b.c. on {y=0} and {y=1}
    """
    
    fdim = domain.topology.dim - 1

    bcs = []
    pbc_directions = []
    pbc_slave_tags = []
    pbc_is_slave = []
    pbc_is_master = []
    pbc_meshtags = []
    pbc_slave_to_master_maps = []
    
    # number of subspaces in functionspace
    N_subspaces = funcspace.num_sub_spaces
    
    # Create MultiPointConstraint object
    mpc = dolfinx_mpc.MultiPointConstraint(funcspace)

    def generate_pbc_slave_to_master_map(i):
        def pbc_slave_to_master_map(x):
            out_x = x.copy() 
            out_x[i] = x[i] - domain.geometry.x.max()
            return out_x
        return pbc_slave_to_master_map

    def generate_pbc_is_slave(i):
        return lambda x: np.isclose(x[i], domain.geometry.x.max())

    def generate_pbc_is_master(i):
        return lambda x: np.isclose(x[i], domain.geometry.x.min())

    # Parse boundary conditions
    def parse_bcs():
        for i, bc_type in enumerate(boundary_condition):
            
            if bc_type == "dirichlet":
                u_bc = fem.Function(functionspace)
                u_bc.x.array[:] = dbc_value # value of dirichlet bc needs to be passed into this function!
    
                def dirichletboundary(x):
                    return np.logical_or(np.isclose(x[i], domain.geometry.x.min()), np.isclose(x[i], domain.geometry.x.max()))
                facets = locate_entities_boundary(domain, fdim, dirichletboundary)
                topological_dofs = fem.locate_dofs_topological(functionspace, fdim, facets)
                bcs.append(fem.dirichletbc(u_bc, topological_dofs))
            
            elif bc_type == "periodic":
                pbc_directions.append(i)
                pbc_slave_tags.append(i + 2)
                pbc_is_slave.append(generate_pbc_is_slave(i))
                pbc_is_master.append(generate_pbc_is_master(i))
                pbc_slave_to_master_maps.append(generate_pbc_slave_to_master_map(i))
    
                facets = locate_entities_boundary(domain, fdim, pbc_is_slave[-1])
                arg_sort = np.argsort(facets)
                pbc_meshtags.append(meshtags(domain,
                                              fdim,
                                              facets[arg_sort],
                                              np.full(len(facets), pbc_slave_tags[-1], dtype=np.int32)))
        
        N_pbc = len(pbc_directions)
        for i in range(N_pbc):
            if N_pbc > 1:
                def pbc_slave_to_master_map(x):
                    out_x = pbc_slave_to_master_maps[i](x)
                    idx = pbc_is_slave[(i + 1) % N_pbc](x)
                    out_x[pbc_directions[i]][idx] = np.nan
                    return out_x
            else:
                pbc_slave_to_master_map = pbc_slave_to_master_maps[i]
    
            mpc.create_periodic_constraint_topological(functionspace, pbc_meshtags[i],
                                                        pbc_slave_tags[i],
                                                        pbc_slave_to_master_map,
                                                        bcs)
            
        if len(pbc_directions) > 1:
            # Map intersection(slaves_x, slaves_y) to intersection(masters_x, masters_y),
            # i.e. map the slave dof at (1, 1) to the master dof at (0, 0)
            def pbc_slave_to_master_map(x):
                out_x = x.copy()
                out_x[0] = x[0] - domain.geometry.x.max()
                out_x[1] = x[1] - domain.geometry.x.max()
                idx = np.logical_and(pbc_is_slave[0](x), pbc_is_slave[1](x))
                out_x[0][~idx] = np.nan
                out_x[1][~idx] = np.nan
                return out_x
            mpc.create_periodic_constraint_topological(functionspace, pbc_meshtags[1],
                                                        pbc_slave_tags[1],
                                                        pbc_slave_to_master_map,
                                                        bcs)

    if N_subspaces == 0:
        functionspace = funcspace
        parse_bcs()
    else:
        for j in range(N_subspaces):
            functionspace = funcspace.sub(j)
            parse_bcs()
            
    mpc.finalize()
    
    return mpc, bcs


mpc, bcs = dirichlet_and_periodic_bcs(domain, V, ["periodic", "periodic"]) 

### fix fluctuation at one point 

def boundary_node(x):
    return np.logical_and( np.isclose(x[0], domain.geometry.x.min()), np.isclose(x[1], domain.geometry.x.min()))

dofs_node_sub0 = locate_dofs_geometrical(V.sub(0).collapse()[0], boundary_node)
dofs_node_sub1 = locate_dofs_geometrical(V.sub(1).collapse()[0], boundary_node)
bc_node_sub0 = dirichletbc(PETSc.ScalarType(0.0), dofs_node_sub0, V.sub(0))
bc_node_sub1 = dirichletbc(PETSc.ScalarType(0.0), dofs_node_sub1, V.sub(1))

bcs.append(bc_node_sub0), bcs.append(bc_node_sub1)


### constant gradients for computational homogenization approach
gradient  = Constant(domain, PETSc.ScalarType((1,0)))
gradient2 = Constant(domain, PETSc.ScalarType((0,1)))

##############################################################################
#                               Variational problem                          #
##############################################################################
# weak formulation
eqn1 = - k_1 * ufl.dot(ufl.grad(ufl.dot(gradient , x_coor) + v_1), ufl.grad(dv_1)) * dx 
eqn2 = - k_2 * ufl.dot(ufl.grad(ufl.dot(gradient2, x_coor) + v_2), ufl.grad(dv_2)) * dx 

eqn = eqn1 + eqn2

# Jacobi determinant
Jac = ufl.derivative(eqn, v_sol, dv_sol)

##############################################################################
#                                      Solve                                 #
##############################################################################

problem = NonlinearMPCProblem(eqn, v_sol, mpc, bcs=bcs, J = Jac)                # define nonlinear problem
solver  = NewtonSolverMPC(comm, problem, mpc)                                   # define solver

# solve 
log.set_log_level(log.LogLevel.INFO) 
number, converged = solver.solve(v_sol)                                         # solve the set of PDEs
assert(converged)
log.set_log_level(log.LogLevel.OFF) 

[dokken]

@tuderic

does it not converge due to the classes for nonlinear MPC problems that seem to be in the testing phase?

It is not in a testing phase, the implementation should solve it (@nate please correct me if Im wrong).

You say that your problem is very simple, but you have supplied 460 lines of code, without any mathematical description of your problem.
This is way too lengthy for me to parse and interpret to look for potential errors.
I would suggest you start with a simpler example, look at the test in DOLFINX_MPC, and extend it bit by bit, until it stops working.

[nate]

Since your example is so long, doesn’t provide a mesh, and for some reason contains copied tests from dolfinx_mpc, I doubt anyone here has time to help.

Nonlinear problems are not straightforward. I agree with @dokken regarding taking things slow and building up your problem ensuring at each step that not only does the code complete execution, but that the numerical result is verifiably correct. Consider discussing with collaborators if you cannot create a minimal working example to post here.

[tuderic]

Thanks to both of you for giving me the advice to create a simple example first. I try to provide a minimal working example:

In this example, we have a quadratic domain with one inclusion in the center. These two subdomains have different material parameters.
I want to solve the steady state heat equation (which is actually in this case the poisson equation). However, in the context of computational homogenization, we split the primary variable (the temperature) into a constant part (a prescribed temperature gradient * x) and a fluctuation part, which is the new primary variable. We can write:

\kappa \nabla^2(\overline{\nabla T} \cdot x + \widetilde{T}) = 0

I want to use periodic boundary conditions (function modified from here). Therefore, I adapted this function. When I solve this linear problem with multipoint constraints, the result seems to be right. However, I tried to solve the same problem with the nonlinear MPC classes. I do not understand two things:

1. When I try to use a boundary condition (in addition to periodic boundary conditions) in order to fix one degree of freedom (as you prevent “rigid body movements” in the mechanical case), the solver diverges. This behavior cannot be seen for the linear solver.
2. When I do not use an additional boundary condition for one DOF, the nonlinear solver converges with one iteration. This should be the case for linear problems. However, the results between the nonlinear and linear solver are totally different.

The code is provided below. The mesh files can be found here.

# Load required libraries
#########################
import dolfinx
from dolfinx import fem, mesh, io, log
import dolfinx_mpc

import numpy as np
from mpi4py import MPI
import ufl
from petsc4py import PETSc
from typing import List
##############################################################################
#                                     Define mesh                            #
##############################################################################
comm = MPI.COMM_WORLD
rank = comm.rank

meshFile = "2D_OneInclusion_v4"                                                 # mesh file name

with io.XDMFFile(comm, "mesh/"+meshFile+"_mesh.xdmf", "r") as xdmf:
    domain = xdmf.read_mesh(name="Grid")                                        # reads the mesh
    ct = xdmf.read_meshtags(domain, name="Grid")                                # cell tags
domain.topology.create_connectivity(domain.topology.dim, domain.topology.dim-1) 

with io.XDMFFile(comm, "mesh/"+meshFile+"_facets.xdmf", "r") as xdmf:
    ft = xdmf.read_meshtags(domain, name="Grid")                                # facet tags

num_cells = domain.topology.index_map(domain.topology.dim).size_local           # number of cells

dx = ufl.Measure('dx', domain=domain, subdomain_data=ct) 
x_coor = ufl.SpatialCoordinate(domain)

##############################################################################
#                                   Function space                           #
##############################################################################

# Elementformulierung
P1 = ufl.FiniteElement("CG", domain.ufl_cell(), 1)  
V = fem.FunctionSpace(domain, P1)                                              

# for linear problem:
T_lin  = ufl.TrialFunction(V)                                                       # fluctuation: primary variable linear problem
# for nonlinear problem:
T_nonl  = fem.Function(V)                                                           # fluctuation: primary variable nonlinear problem
dT_nonl  = ufl.TrialFunction(V)  

# testfunction      
v = ufl.TestFunction(V)                                               

##############################################################################
#                                   Randbedinungen                           #
##############################################################################

### periodic boundary conditions 
if(1):
    mpc, bcs = dirichlet_and_periodic_bcs(domain, V, ["periodic", "periodic"]) 

### Gradient
T_gradient = fem.Constant(domain, PETSc.ScalarType((1,0)))

### fix fluctuation on one point
def boundary_node(x):
    return np.logical_and( np.isclose(x[0], domain.geometry.x.min()), np.isclose(x[1], domain.geometry.x.min()))
dofs_node = fem.locate_dofs_geometrical(V, boundary_node)
bc_node = fem.dirichletbc(PETSc.ScalarType(0.0), dofs_node, V)

# this cause divergence for nonlinear solving ... Why?
if(1):    
    bcs.append(bc_node)
    
#############################################################################
#                               Definition von Parametern                    #
##############################################################################

S = fem.FunctionSpace(domain, ("DG", 0))
kappa = fem.Function(S)

for i in range(num_cells):
    if ct.values[i] == 1: # matrix
        kappa.x.array[i] = 1
    elif ct.values[i] == 2: # inclusions
        kappa.x.array[i] = 5
    
##############################################################################
#                                 Variationsproblem                          #
##############################################################################

# linear form
eqn_lin  = kappa * ufl.dot(ufl.grad(ufl.dot(T_gradient , x_coor) + T_lin), ufl.grad(v)) * dx
# nonlinear form
eqn_nonl = kappa * ufl.dot(ufl.grad(ufl.dot(T_gradient , x_coor) + T_nonl), ufl.grad(v)) * dx
Jac = ufl.derivative(eqn_nonl, T_nonl, dT_nonl)

##############################################################################
#                                  Lösungsschleife                           #
##############################################################################

### linear solver
a_form = ufl.lhs(eqn_lin)
L_form = ufl.rhs(eqn_lin)
problem_lin = dolfinx_mpc.LinearProblem(a_form, L_form, mpc, bcs=bcs)

T_lin = problem_lin.solve()

### nonlinear solver
problem_nonl = NonlinearMPCProblem(eqn_nonl, T_nonl, mpc, bcs=bcs, J = Jac)
solver_nonl  = NewtonSolverMPC(domain.comm, problem_nonl, mpc)

log.set_log_level(log.LogLevel.INFO) 
number, converged = solver_nonl.solve(T_nonl)
assert(converged)
log.set_log_level(log.LogLevel.OFF) 

if(0):
    with io.XDMFFile(domain.comm, "results.xdmf", "w") as xdmf:
        xdmf.write_mesh(domain)
        xdmf.write_function(T_lin) 
    
    with io.XDMFFile(domain.comm, "results_nonl.xdmf", "w") as xdmf:
        xdmf.write_mesh(domain)
        xdmf.write_function(T_nonl) 

This function is used to create periodic boundary conditions. It should work since reasonable results appear for linear solving.

# adapted from https://github.com/jorgensd/dolfinx_mpc/blob/main/python/demos/demo_periodic_gep.py
def dirichlet_and_periodic_bcs(domain, funcspace, boundary_condition: List[str] = ["dirichlet", "periodic"], dbc_value = 0):
    """
    - Function to set either dirichlet or periodic boundary conditions.
    - This function is was taken from https://github.com/jorgensd/dolfinx_mpc/blob/main/python/demos/demo_periodic_gep.py and adapted.
    ----------
    boundary_condition
        First item describes b.c. on {x=0} and {x=1}
        Second item describes b.c. on {y=0} and {y=1} 
    """
    
    fdim = domain.topology.dim - 1

    bcs = []
    pbc_directions = []
    pbc_slave_tags = []
    pbc_is_slave = []
    pbc_is_master = []
    pbc_meshtags = []
    pbc_slave_to_master_maps = []
    
    # number of subspaces in functionspace
    N_subspaces = funcspace.num_sub_spaces
    
    # Create MultiPointConstraint object
    mpc = dolfinx_mpc.MultiPointConstraint(funcspace)

    def generate_pbc_slave_to_master_map(i):
        def pbc_slave_to_master_map(x):
            out_x = x.copy() 
            out_x[i] = x[i] - domain.geometry.x[:,i].max()
            return out_x
        return pbc_slave_to_master_map

    def generate_pbc_is_slave(i):
        return lambda x: np.isclose(x[i], domain.geometry.x[:,i].max())

    def generate_pbc_is_master(i):
        return lambda x: np.isclose(x[i], domain.geometry.x[:,i].min())

    # Parse boundary conditions
    def parse_bcs():
        for i, bc_type in enumerate(boundary_condition):
            
            if bc_type == "dirichlet":
               
    
                def dirichletboundary(x):
                    return np.logical_or(np.isclose(x[i], domain.geometry.x[:,i].min()), np.isclose(x[i], domain.geometry.x[:,i].max()))
                
                if N_subspaces == 0:
                    u_bc = fem.Function(functionspace)
                    u_bc.x.array[:] = dbc_value # value of dirichlet bc needs to be passed into this function!
                    facets = mesh.locate_entities_boundary(domain, fdim, dirichletboundary)
                    topological_dofs = fem.locate_dofs_topological(functionspace, fdim, facets)
                    bcs.append(fem.dirichletbc(u_bc, topological_dofs))
                else: 
                    geometrical_dofs = fem.locate_dofs_geometrical(functionspace.collapse()[0], dirichletboundary)
                    bcs.append(fem.dirichletbc(PETSc.ScalarType(dbc_value), geometrical_dofs, functionspace))
            
            elif bc_type == "periodic":
                pbc_directions.append(i)
                pbc_slave_tags.append(i + 2)
                pbc_is_slave.append(generate_pbc_is_slave(i))
                pbc_is_master.append(generate_pbc_is_master(i))
                pbc_slave_to_master_maps.append(generate_pbc_slave_to_master_map(i))
    
                facets = mesh.locate_entities_boundary(domain, fdim, pbc_is_slave[-1])
                arg_sort = np.argsort(facets)
                pbc_meshtags.append(mesh.meshtags(domain,
                                              fdim,
                                              facets[arg_sort],
                                              np.full(len(facets), pbc_slave_tags[-1], dtype=np.int32)))
                
            elif bc_type == "neumann":
                pass
        
        N_pbc = len(pbc_directions)
        for i in range(N_pbc):
            if N_pbc > 1:
                def pbc_slave_to_master_map(x):
                    out_x = pbc_slave_to_master_maps[i](x)
                    idx = pbc_is_slave[(i + 1) % N_pbc](x)
                    out_x[pbc_directions[i]][idx] = np.nan
                    return out_x
            else:
                pbc_slave_to_master_map = pbc_slave_to_master_maps[i]
    
            mpc.create_periodic_constraint_topological(functionspace, pbc_meshtags[i],
                                                        pbc_slave_tags[i],
                                                        pbc_slave_to_master_map,
                                                        bcs)
            
        if len(pbc_directions) > 1:
            # Map intersection(slaves_x, slaves_y) to intersection(masters_x, masters_y),
            # i.e. map the slave dof at (1, 1) to the master dof at (0, 0)
            def pbc_slave_to_master_map(x):
                out_x = x.copy()
                out_x[0] = x[0] - domain.geometry.x[:,0].max()
                out_x[1] = x[1] - domain.geometry.x[:,1].max()
                idx = np.logical_and(pbc_is_slave[0](x), pbc_is_slave[1](x))
                out_x[0][~idx] = np.nan
                out_x[1][~idx] = np.nan
                return out_x
            mpc.create_periodic_constraint_topological(functionspace, pbc_meshtags[1],
                                                        pbc_slave_tags[1],
                                                        pbc_slave_to_master_map,
                                                        bcs)

    if N_subspaces == 0:
        functionspace = funcspace
        parse_bcs()
    else:
        for j in range(N_subspaces):
            functionspace = funcspace.sub(j)
            parse_bcs()
            
    mpc.finalize()
    
    return mpc, bcs

These classes are used for NonlinearMPC problems.

# copied from https://github.com/jorgensd/dolfinx_mpc/blob/main/python/tests/test_nonlinear_assembly.py
class NonlinearMPCProblem(dolfinx.fem.petsc.NonlinearProblem):

    def __init__(self, F, u, mpc, bcs=[], J=None, form_compiler_params={},
                  jit_params={}):
        self.mpc = mpc
        super().__init__(F, u, bcs=bcs, J=J,
                          form_compiler_params=form_compiler_params,
                          jit_params=jit_params)
        
    def F(self, x: PETSc.Vec, F: PETSc.Vec):
        with F.localForm() as F_local:
            F_local.set(0.0)
        dolfinx_mpc.assemble_vector(self._L, self.mpc, b=F)

        # Apply boundary condition
        dolfinx_mpc.apply_lifting(F, [self._a], bcs=[self.bcs],
                                  constraint=self.mpc, x0=[x], scale=-1.0)
        F.ghostUpdate(addv=PETSc.InsertMode.ADD, mode=PETSc.ScatterMode.REVERSE)
        dolfinx.fem.petsc.set_bc(F, self.bcs, x, -1.0)

    def J(self, x: PETSc.Vec, A: PETSc.Mat):
        A.zeroEntries()
        dolfinx_mpc.assemble_matrix(self._a, self.mpc, bcs=self.bcs, A=A)
        A.assemble()


class NewtonSolverMPC(dolfinx.cpp.nls.petsc.NewtonSolver):
    def __init__(self, comm: MPI.Intracomm, problem: NonlinearMPCProblem,
                 mpc: dolfinx_mpc.MultiPointConstraint):
        """A Newton solver for non-linear MPC problems."""
        super().__init__(comm)
        self.mpc = mpc
        self.u_mpc = dolfinx.fem.Function(mpc.function_space)

        # Create matrix and vector to be used for assembly of the non-linear
        # MPC problem
        self._A = dolfinx_mpc.cpp.mpc.create_matrix(
            problem.a, mpc._cpp_object)
        self._b = dolfinx.cpp.la.petsc.create_vector(
            mpc.function_space.dofmap.index_map,
            mpc.function_space.dofmap.index_map_bs)

        self.setF(problem.F, self._b)
        self.setJ(problem.J, self._A)
        self.set_form(problem.form)
        self.set_update(self.update)

    def update(self, solver: dolfinx.nls.petsc.NewtonSolver,
               dx: PETSc.Vec, x: PETSc.Vec):
        # We need to use a vector created on the MPC's space to update ghosts
        self.u_mpc.vector.array = x.array_r
        self.u_mpc.vector.axpy(-1.0, dx)
        self.u_mpc.vector.ghostUpdate(addv=PETSc.InsertMode.INSERT,
                                      mode=PETSc.ScatterMode.FORWARD)
        self.mpc.homogenize(self.u_mpc.vector)
        self.mpc.backsubstitution(self.u_mpc.vector)
        x.array = self.u_mpc.vector.array_r
        x.ghostUpdate(addv=PETSc.InsertMode.INSERT,
                      mode=PETSc.ScatterMode.FORWARD)

    def solve(self, u: dolfinx.fem.Function):
        """Solve non-linear problem into function u. Returns the number
        of iterations and if the solver converged."""
        n, converged = super().solve(u.vector)
        u.x.scatter_forward()
        return n, converged

    @property
    def A(self) -> PETSc.Mat:
        """Jacobian matrix"""
        return self._A

    @property
    def b(self) -> PETSc.Vec:
        """Residual vector"""
        return self._b

[tuderic]

Hi everybody,
I still do not understand why the heat equation with periodic boundary conditions does not converge (to the right solution) with the solver for nonlinear MPC problems. Does anyone have any ideas? Alternatively, how can I simplify the problem even more that I can find the mistake?