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