# Assembly of matrices for use in external solver

**URL:** <https://fenicsproject.discourse.group/t/assembly-of-matrices-for-use-in-external-solver/16900>\
**Category:** General\
**Tags:** boundary-conditions, assembly, dolfinx\
**Created:** [February 2, 2025, 2:25pm UTC](https://fenicsproject.discourse.group/t/assembly-of-matrices-for-use-in-external-solver/16900 "2025-02-02T14:25:19Z")\
**Posts on this page:** 5\
**Page:** 1

<div class="post-metadata">

**Author:** ![Fruitpunchsamurai](https://avatars.discourse-cdn.com/v4/letter/f/9dc877/32.png) [@Fruitpunchsamurai](https://fenicsproject.discourse.group/u/Fruitpunchsamurai)\
**Post date:** [February 2, 2025, 2:25pm UTC](https://fenicsproject.discourse.group/t/assembly-of-matrices-for-use-in-external-solver/16900/1 "2025-02-02T14:25:19Z")

</div>

I am trying to implement the example in this article [Lombardi, D., Pagliantini, C. (2024). Conformal variational discretisation of infinite dimensional Hamiltonian systems with gradient flow dissipation], specifically the one about solving Navier-Stokes in vorticity formulation on the torus.  
For simplicity, for now I just want to solve M \dot{y}=-\nu K y where M is the mass matrix and K is the stiffness matrix. My aim is to build the M and K matrices using `dolfinx_mpc` (imposing pure periodic boundary conditions), then use an external time-stepping method for the ODE above (for now I just use `solve_ivp`, but ideally I would like to test my own method on this problem).  
The issue I believe I am encountering is that the solution I obtain does not satisfy the periodic boundary conditions. Below a MWE.  
I am quite new to Fenics, so any kind of insight is greatly appreciated.

```auto
from mpi4py import MPI
from dolfinx import mesh
from dolfinx import fem
from dolfinx import plot
from dolfinx.fem import functionspace
from dolfinx_mpc import assemble_matrix
import dolfinx_mpc
import dolfinx_mpc.utils
import ufl
import numpy as np
import scipy as sp
import scipy.integrate as integ
import pyvista

# =============================================================================
# Function for the right hand side
# =============================================================================

def FOM(t, y):
    out = np.zeros(N)
    out = -nu * K @ y
    out = sp.sparse.linalg.spsolve(M, out)
    
    uh = fem.Function(mpc.function_space) # Create function to store solution
    uh.x.array[:] = out
    uh.x.scatter_forward()
    mpc.backsubstitution(uh)
    return uh.x.array

# Parameters
nu = 1e-2
lmb = 25
nx = ny = 20

domain = mesh.create_rectangle(MPI.COMM_WORLD, [np.array([0, 0]), np.array([2*np.pi, 2*np.pi])],
                               [nx, ny], mesh.CellType.triangle)
 

V = functionspace(domain, ("Lagrange", 1))
fdim = domain.topology.dim - 1

tdim = domain.topology.dim
fdim = tdim - 1
domain.topology.create_connectivity(fdim, tdim)
boundary_facets = mesh.exterior_facet_indices(domain.topology)
boundary_dofs = fem.locate_dofs_topological(V, fdim, boundary_facets)

# No Dirichlet boundary conditions
bcs = []

# Periodic boundary conditions
def top_boundary(x):
    return np.isclose(x[1], 2*np.pi)

def right_boundary(x):
    return np.isclose(x[0], 2*np.pi)

# Periodic boundary conditions
def periodic_relation_left_right(x):
    out_x = np.zeros(x.shape)
    out_x[0] = x[0] - 2*np.pi
    out_x[1] = x[1]
    return out_x

def periodic_relation_bottom_top(x):
    out_x = np.zeros(x.shape)
    out_x[0] = x[0]
    out_x[1] = x[1] - 2*np.pi
    return out_x

mpc = dolfinx_mpc.MultiPointConstraint(V)
mpc.create_periodic_constraint_geometrical(
    V, right_boundary, periodic_relation_left_right, bcs
)
mpc.create_periodic_constraint_geometrical(
    V, top_boundary, periodic_relation_bottom_top, bcs
)
mpc.finalize()

# Mass and stiffness matrices assembly
phi_i = ufl.TrialFunction(V)
phi_j = ufl.TestFunction(V)

k = ufl.inner(ufl.grad(phi_i), ufl.grad(phi_j)) * ufl.dx
m = ufl.inner(phi_i, phi_j) * ufl.dx
k_compiled = fem.form(k)
m_compiled = fem.form(m)
# Get Scipy analogues for K and M
K = dolfinx_mpc.utils.gather_PETScMatrix(assemble_matrix(k_compiled, mpc, bcs=bcs))
M = dolfinx_mpc.utils.gather_PETScMatrix(assemble_matrix(m_compiled, mpc, bcs=bcs))

# Define initial stream function
psi0 = fem.Function(V)
psi0.interpolate(lambda x: 0.25*np.cos(3*x[0])*np.sin(4*x[1]) - 0.2*np.cos(5*x[1]) - 0.2*np.sin(5*x[0]))
y0 = - lmb*psi0.x.array
    
# =============================================================================
# Solving the FOM
# =============================================================================

# Time stepping grid
T = 10
h = 0.01
t0 = np.arange(0,T,h)
L = t0.size
N = V.dofmap.index_map.size_global * V.dofmap.index_map_bs

Y = np.zeros((N, L)) # snapshot matrix
    
# solve the system with initial condition y0
sol = integ.solve_ivp(FOM, (0,T), y0, atol=1e-10, rtol=1e-10, t_eval=t0, method='LSODA')
Y = sol.y

# =============================================================================
# Visualizing the solution on the mesh
# =============================================================================

# Plot mesh
u_topology, u_cell_types, u_geometry = plot.vtk_mesh(V)
u_grid = pyvista.UnstructuredGrid(u_topology, u_cell_types, u_geometry)
u_grid.point_data["u_ref"] = Y[:,-1]
u_grid.set_active_scalars("u_ref")
u_plotter = pyvista.Plotter()
u_plotter.add_mesh(u_grid, show_edges=True)
u_plotter.view_xy()
if not pyvista.OFF_SCREEN:
    u_plotter.show()
# 3D plot
warped = u_grid.warp_by_scalar()
plotter2 = pyvista.Plotter()
plotter2.add_mesh(warped, show_edges=True, show_scalar_bar=True)
if not pyvista.OFF_SCREEN:
    plotter2.show()

```

---

<div class="post-metadata">

**Author:** ![dokken](https://yyz2.discourse-cdn.com/free1/user_avatar/fenicsproject.discourse.group/dokken/32/1560_2.png) [@dokken](https://fenicsproject.discourse.group/u/dokken)\
**Post date:** [February 2, 2025, 3:11pm UTC](https://fenicsproject.discourse.group/t/assembly-of-matrices-for-use-in-external-solver/16900/2 "2025-02-02T15:11:32Z")

</div>

> [@Fruitpunchsamurai](#):
>
> `out = -nu * K @ y`

This to me doesn’t look quite right. The assembled matrix K is made with

> [@Fruitpunchsamurai](#):
>
> `K = dolfinx_mpc.utils.gather_PETScMatrix(assemble_matrix(k_compiled, mpc, bcs=bcs))`

which ends up being P^T K ( P y) where P is a restriction operator.

As I’m not at a computer, I can’t give further guidance atm. I would suggest having a look at:

> <https://github.com/jorgensd/dolfinx_mpc/blob/main/python/tests/test_mpc_pipeline.py#L170>

That compares the global restriction matrix approach with the one in DOLFINx MPC.  
This also might be a relevant test: [dolfinx\_mpc/python/tests/test\_vector\_assembly.py at main · jorgensd/dolfinx\_mpc · GitHub](https://github.com/jorgensd/dolfinx_mpc/blob/main/python/tests/test_vector_assembly.py#L60)  
[dolfinx\_mpc/python/dolfinx\_mpc/utils/test.py at main · jorgensd/dolfinx\_mpc · GitHub](https://github.com/jorgensd/dolfinx_mpc/blob/main/python/dolfinx_mpc/utils/test.py#L281)

---

<div class="post-metadata">

**Author:** ![Fruitpunchsamurai](https://avatars.discourse-cdn.com/v4/letter/f/9dc877/32.png) [@Fruitpunchsamurai](https://fenicsproject.discourse.group/u/Fruitpunchsamurai)\
**Post date:** [February 12, 2025, 12:58pm UTC](https://fenicsproject.discourse.group/t/assembly-of-matrices-for-use-in-external-solver/16900/3 "2025-02-12T12:58:33Z")

</div>

Thank you for your suggestions! If I understand correctly, the suggested strategy is to compute the restricted operators (in my case K\_r=P^TKP and M\_r=P^TMP) which contain values only for the non-slave dofs, and solve the restricted ODE M\_r\dot{y}\_r=-\nu K\_r y\_r. The final solution is then recovered as y=Py\_r.  
I tried implementing this changes, but I still think there is something wrong.  
For comparison, I tried to solve the same ODE using the implicit midpoint method both by constructing the restricted operators and by following this [tutorial](https://jsdokken.com/dolfinx-tutorial/chapter2/diffusion_code.html). I obtain reasonable results only in the second way. Below you can find a MWE. Any suggestion is appreciated!

```auto
from mpi4py import MPI
from dolfinx import mesh
from dolfinx import fem
from dolfinx import plot
from dolfinx.fem import functionspace
import dolfinx_mpc
from petsc4py import PETSc
from dolfinx_mpc import assemble_matrix, assemble_vector
import dolfinx_mpc.utils
import ufl
import numpy as np
import scipy as sp
import scipy.integrate as integ
import pyvista

# =============================================================================
# Function for the right hand side
# =============================================================================

def implicit_midpoint(interval, y0, h):
    t = np.arange(interval[0], interval[1], h)
    y = np.zeros((len(y0),len(t)))
    y[:,0] = y0;
    for i in np.arange(1,len(t)):
        rhs = (M_mpc - 0.5*h*nu*K_mpc) @ y[:,i-1]
        y[:,i] = sp.sparse.linalg.spsolve(M_mpc + 0.5*h*nu*K_mpc, rhs)
    return y
    

# Parameters
nu = 1e-2
lmb = 25
nx = ny = 20

domain = mesh.create_rectangle(MPI.COMM_WORLD, [np.array([0, 0]), np.array([2*np.pi, 2*np.pi])],
                               [nx, ny], mesh.CellType.triangle)
 

V = functionspace(domain, ("Lagrange", 1))
fdim = domain.topology.dim - 1

tdim = domain.topology.dim
fdim = tdim - 1
domain.topology.create_connectivity(fdim, tdim)
boundary_facets = mesh.exterior_facet_indices(domain.topology)
boundary_dofs = fem.locate_dofs_topological(V, fdim, boundary_facets)

bcs = []

# Periodic boundary conditions
def top_boundary(x):
    return np.isclose(x[1], 2*np.pi)

def right_boundary(x):
    return np.isclose(x[0], 2*np.pi)

# Periodic boundary conditions
def periodic_relation_left_right(x):
    out_x = np.zeros(x.shape)
    out_x[0] = x[0] - 2*np.pi
    out_x[1] = x[1]
    return out_x

def periodic_relation_bottom_top(x):
    out_x = np.zeros(x.shape)
    out_x[0] = x[0]
    out_x[1] = x[1] - 2*np.pi
    return out_x

mpc = dolfinx_mpc.MultiPointConstraint(V)
mpc.create_periodic_constraint_geometrical(
    V, right_boundary, periodic_relation_left_right, bcs
)
mpc.create_periodic_constraint_geometrical(
    V, top_boundary, periodic_relation_bottom_top, bcs
)
mpc.finalize()

# Mass and stiffness matrices assembly
phi_i = ufl.TrialFunction(V)
phi_j = ufl.TestFunction(V)

k = ufl.inner(ufl.grad(phi_i), ufl.grad(phi_j)) * ufl.dx
m = ufl.inner(phi_i, phi_j) * ufl.dx
k_compiled = fem.form(k)
m_compiled = fem.form(m)
# Get Scipy analogues for K and M
K = fem.petsc.assemble_matrix(k_compiled, bcs)
K.assemble()
M = fem.petsc.assemble_matrix(m_compiled, bcs)
M.assemble()
K = dolfinx_mpc.utils.gather_PETScMatrix(K)
M = dolfinx_mpc.utils.gather_PETScMatrix(M)

# Compute trasformation matrix
P = dolfinx_mpc.utils.gather_transformation_matrix(mpc)

# Compute the restricted matrices
M_mpc = P.T @ M @ P
K_mpc = P.T @ K @ P

# Define initial stream function
psi0 = fem.Function(V)
psi0.interpolate(lambda x: 0.25*np.cos(3*x[0])*np.sin(4*x[1]) - 0.2*np.cos(5*x[1]) - 0.2*np.sin(5*x[0]))
y0 = - lmb*psi0.x.array
    
# =============================================================================
# Solving the FOM
# =============================================================================

# Time stepping grid
T = 1
h = 0.01
t0 = np.arange(0,T,h)
L = t0.size
N = V.dofmap.index_map.size_global * V.dofmap.index_map_bs

Y = np.zeros((N, L)) # snapshot matrix

# solve with implicit midpoint
Y = implicit_midpoint([0,T], P.T @ y0, h)
y_mpc_imp = P @ Y[:,-1]

# =============================================================================
# Solution with the classical Fenics approach
# =============================================================================
u_n = fem.Function(V)
u_n.name = "u_n"
u_n.interpolate(lambda x: 0.25*np.cos(3*x[0])*np.sin(4*x[1]) - 0.2*np.cos(5*x[1]) - 0.2*np.sin(5*x[0]))

u, v = ufl.TrialFunction(V), ufl.TestFunction(V)
a = ufl.inner(u, v)* ufl.dx + 0.5*h*nu * ufl.dot(ufl.grad(u), ufl.grad(v)) * ufl.dx
L = u_n * v * ufl.dx - 0.5*h*nu * ufl.dot(ufl.grad(u_n), ufl.grad(v)) * ufl.dx

bilinear_form = fem.form(a)
linear_form = fem.form(L)
A = assemble_matrix(bilinear_form, mpc)
b = assemble_vector(linear_form, mpc)
b.ghostUpdate(addv=PETSc.InsertMode.ADD_VALUES, mode=PETSc.ScatterMode.REVERSE)

solver = PETSc.KSP().create(domain.comm)
solver.setOperators(A)
solver.setType(PETSc.KSP.Type.PREONLY)
solver.getPC().setType(PETSc.PC.Type.LU)

uh = fem.Function(V)
uh.name = "uh"
uh.interpolate(lambda x: 0.25*np.cos(3*x[0])*np.sin(4*x[1]) - 0.2*np.cos(5*x[1]) - 0.2*np.sin(5*x[0]))

t_range = np.arange(0.0, T, h)
num_steps = len(t_range)
for i in range(num_steps):

    # Update the right hand side reusing the initial vector
    with b.localForm() as loc_b:
        loc_b.set(0)
    b = assemble_vector(linear_form, mpc)

    # Apply Dirichlet boundary condition to the vector
    b.ghostUpdate(addv=PETSc.InsertMode.ADD_VALUES, mode=PETSc.ScatterMode.REVERSE)

    # Solve linear problem
    solver.solve(b, uh.x.petsc_vec)
    uh.x.scatter_forward()
    mpc.backsubstitution(uh)

    # Update solution at previous time step (u_n)
    u_n.x.array[:] = uh.x.array
# =============================================================================
# Visualizing the solution on the mesh
# =============================================================================
print(sp.linalg.norm(uh.x.array - y_mpc_imp))

# Plot mesh
u_topology, u_cell_types, u_geometry = plot.vtk_mesh(V)
u_grid = pyvista.UnstructuredGrid(u_topology, u_cell_types, u_geometry)
u_grid.point_data["u_ref"] = uh.x.array
u_grid.set_active_scalars("u_ref")
u_plotter = pyvista.Plotter()
u_plotter.add_mesh(u_grid, show_edges=True)
u_plotter.view_xy()
if not pyvista.OFF_SCREEN:
    u_plotter.show()
# 3D plot
warped = u_grid.warp_by_scalar()
plotter2 = pyvista.Plotter()
plotter2.add_mesh(warped, show_edges=True, show_scalar_bar=True)
if not pyvista.OFF_SCREEN:
    plotter2.show()

```

---

<div class="post-metadata">

**Author:** ![dokken](https://yyz2.discourse-cdn.com/free1/user_avatar/fenicsproject.discourse.group/dokken/32/1560_2.png) [@dokken](https://fenicsproject.discourse.group/u/dokken)\
**Post date:** [February 12, 2025, 1:05pm UTC](https://fenicsproject.discourse.group/t/assembly-of-matrices-for-use-in-external-solver/16900/4 "2025-02-12T13:05:00Z")

</div>

Could you first ensure that it works for a time-independent problem, i.e. consider: [dolfinx\_mpc/python/tests/test\_linear\_problem.py at main · jorgensd/dolfinx\_mpc · GitHub](https://github.com/jorgensd/dolfinx_mpc/blob/main/python/tests/test_linear_problem.py#L88-L101) which checks all the reduction matrices.

---

<div class="post-metadata">

**Author:** ![Fruitpunchsamurai](https://avatars.discourse-cdn.com/v4/letter/f/9dc877/32.png) [@Fruitpunchsamurai](https://fenicsproject.discourse.group/u/Fruitpunchsamurai)\
**Post date:** [February 18, 2025, 10:17am UTC](https://fenicsproject.discourse.group/t/assembly-of-matrices-for-use-in-external-solver/16900/5 "2025-02-18T10:17:37Z")

</div>

Thank you for the answer. I have tried to solve the Poisson equation both with Dirichlet + periodic and with pure periodic boundary conditions (in this second case I also imposed Dirichlet conditions on the four corners of the domain in order to make the stiffness matrix invertible). In both cases I get the expected solution. For reference here is the code in the case of pure periodic conditions:

```auto
from mpi4py import MPI
from dolfinx import mesh
from dolfinx import fem
from dolfinx import plot
from dolfinx.fem import functionspace
import dolfinx_mpc
import dolfinx_mpc.utils
from petsc4py import PETSc
import ufl
import numpy as np
import scipy as sp
import pyvista
from typing import List

# Parameters
nx = ny = 50

domain = mesh.create_rectangle(MPI.COMM_WORLD, [np.array([0, 0]), np.array([2*np.pi, 2*np.pi])],
                               [nx, ny], mesh.CellType.triangle)
 
V = functionspace(domain, ("Lagrange", 1))
gdim = domain.topology.dim
fdim = domain.topology.dim - 1

# Dirichlet boundary conditions
point_dof1 = fem.locate_dofs_geometrical(V, lambda x: np.logical_and(np.isclose(x[0], 0), np.isclose(x[1], 0)))
point_dof2 = fem.locate_dofs_geometrical(V, lambda x: np.logical_and(np.isclose(x[0], 0), np.isclose(x[1], 2*np.pi)))
point_dof3 = fem.locate_dofs_geometrical(V, lambda x: np.logical_and(np.isclose(x[0], 2*np.pi), np.isclose(x[1], 0)))
point_dof4 = fem.locate_dofs_geometrical(V, lambda x: np.logical_and(np.isclose(x[0], 2*np.pi), np.isclose(x[1], 2*np.pi)))
uD = fem.Constant(domain, 1.0)
bcs = [fem.dirichletbc(uD, point_dof1, V), fem.dirichletbc(uD, point_dof2, V),
       fem.dirichletbc(uD, point_dof3, V), fem.dirichletbc(uD, point_dof4, V)]

# Periodic boundary conditions
boundary_condition: List[str] = ["periodic", "periodic"]

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

# Create MultiPointConstraint object
mpc = dolfinx_mpc.MultiPointConstraint(V)

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())

for i, bc_type in enumerate(boundary_condition):

    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)))

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(V, pbc_meshtags[i],
                                               pbc_slave_tags[i],
                                               pbc_slave_to_master_map,
                                               bcs)
mpc.finalize()

# Poisson equation
u = ufl.TrialFunction(V)
v = ufl.TestFunction(V)
a = ufl.inner(ufl.grad(u), ufl.grad(v)) * ufl.dx
x = ufl.SpatialCoordinate(domain)
f = ufl.cos(x[0]) + ufl.sin(x[1])
rhs = ufl.inner(f, v) * ufl.dx

# Solve with fenics solver
problem = dolfinx_mpc.LinearProblem(a, rhs, mpc, bcs=bcs, petsc_options={"ksp_rtol": 1e-12})
uh0 = problem.solve()

# Solve with custom solver
a_compiled = fem.form(a)
b_compiled = fem.form(rhs)

# Get Scipy analogues for K and M
A = fem.petsc.assemble_matrix(a_compiled, bcs=bcs)
A.assemble()
b = fem.petsc.assemble_vector(b_compiled)
fem.petsc.apply_lifting(b, [a_compiled], [bcs])
b.ghostUpdate(addv=PETSc.InsertMode.ADD_VALUES, mode=PETSc.ScatterMode.REVERSE)
fem.petsc.set_bc(b, bcs)
A = dolfinx_mpc.utils.gather_PETScMatrix(A)
b = dolfinx_mpc.utils.gather_PETScVector(b)

# Compute trasformation matrix
P = dolfinx_mpc.utils.gather_transformation_matrix(mpc)

A_mpc = P.T @ A @ P
b_mpc = P.T @ b
uh_arr = sp.sparse.linalg.spsolve(A_mpc, b_mpc)
uh_arr = P @ uh_arr
    
# Define exact solution and print errors
u_e = fem.Function(V)
u_e.interpolate(lambda x: np.cos(x[0]) + np.sin(x[1]))
    
print(sp.linalg.norm(u_e.x.array - uh_arr))
print(sp.linalg.norm(u_e.x.array - uh0.x.array))

```

I also modified the time-dependent code by imposing the periodic conditions as done [here](https://github.com/jorgensd/dolfinx_mpc/blob/1dc49f5f2b2c209a056a7eed6821838afe278333/python/demos/demo_periodic_gep.py#L239), but the results are still not reasonable.