Periodic boundary condition with non-periodic mesh using dolfinx_mpc

General
1

I am trying to solve a Poisson problem with periodic boundary conditions. The meshes I am importing are 3D boxes, but generally these meshes will not be periodic (ie. nodes and/or element edges do not match identically on opposite faces). As a simplified example, I have tried solving the following “manufactured” problem on a 2D rectangular mesh:

\Delta u = -\frac{4\pi}{L^2}\sin\left(\frac{2\pi y}{L}\right) ~~\textrm{in }\Omega=[0,1] \times [0,1] \\ u = \sin\left(\frac{2\pi y}{L}\right)~~\textrm{on }\Gamma_D = \left\{(0,y) \cup (1,y) \right\} \\ u(x,0) = u(x,1) ~~\textrm{on }\Gamma_P = \left\{(x,0) \cup (x,1) \right\} \\

for which the exact solution is u_{e} = \sin\left(\frac{2\pi y}{L}\right). Below, I attempt to solve this problem twice: first with matching nodes on the upper and lower portions of boundary \Gamma_P, and then with non-matching nodes.

from math import pi
import os
import numpy as np
from mpi4py import MPI
import ufl
from dolfinx import mesh, fem, io, common
import dolfinx_mpc
import gmsh
import meshio

comm = MPI.COMM_WORLD
comm_rank = comm.rank

gmsh.initialize()
gmsh.option.setNumber("General.Terminal", 0)

# Domain size and problem definition
L = 1
def u_exact(mod): # mod is either np or ufl
    return lambda x : mod.sin(2*pi*x[1]/L)

# Solve Poisson problem with periodic BC's
def solve_periodic(L, h1, h2, i):
    # L - size of square domain
    # h1 - mesh size 1
    # h2 - mesh size 2
    # i - counter for file

    # Create mesh
    if comm_rank == 0:
        gmsh.model.add("mesh_{}".format(i))
        p1 = gmsh.model.geo.addPoint(0, 0, 0, h1)
        p2 = gmsh.model.geo.addPoint(L, 0, 0, h2)
        p3 = gmsh.model.geo.addPoint(0, L, 0, h2)
        p4 = gmsh.model.geo.addPoint(L, L, 0, h1)
        c1 = gmsh.model.geo.addLine(p1,p2)
        c2 = gmsh.model.geo.addLine(p1,p3)
        c3 = gmsh.model.geo.addLine(p2,p4)
        c4 = gmsh.model.geo.addLine(p3,p4)
        cl1 = gmsh.model.geo.addCurveLoop([c1,c3,-c4,-c2])
        s1 = gmsh.model.geo.addPlaneSurface([cl1])
        pg1 = gmsh.model.addPhysicalGroup(2, [s1])
        gmsh.model.geo.synchronize()
        gmsh.model.mesh.generate(2)
        gmsh.write("mesh.msh")
        msh = meshio.read("mesh.msh")
        meshio.write("mesh.xdmf", meshio.Mesh(points=msh.points[:,0:2], cells={"triangle":msh.get_cells_type("triangle")})) # add cell data if physical groups specified
    with io.XDMFFile(comm, "mesh.xdmf", "r") as infile:
        domain = infile.read_mesh(name="Grid")

    # Boundary definitions
    b_Left = lambda x : np.isclose(x[0], 0)
    b_Right = lambda x : np.isclose(x[0], L)
    b_Bot = lambda x : np.isclose(x[1], 0)
    b_Top = lambda x : np.isclose(x[1], L)

    # Function space
    V = fem.FunctionSpace(domain, ("CG", 1))

    # Dirichlet boundary conditions - sinusoidal on left and right ends
    dofs_L = fem.locate_dofs_geometrical(V, b_Left)
    u_L = fem.Function(V)
    u_L.interpolate(u_exact(np))
    bc_L = fem.dirichletbc(u_L, dofs_L)
    dofs_R = fem.locate_dofs_geometrical(V, b_Right)
    u_R = fem.Function(V)
    u_R.interpolate(u_exact(np))
    bc_R = fem.dirichletbc(u_R, dofs_R)
    bcs = [bc_L, bc_R]

    periodic_boundary = b_Top # specify the "master" surface as bottom (note - this is function pointer)

    # Map master (bottom) surface to slave (top) surface (note: this is opposite of old fenics)
    def periodic_relation(x):
        out_x = np.zeros(x.shape)
        out_x[0] = x[0]
        out_x[1] = L - x[1]
        out_x[2] = x[2]
        return out_x

    # Periodic boundary - topological
    facets = mesh.locate_entities_boundary(domain, domain.topology.dim - 1, periodic_boundary)
    arg_sort = np.argsort(facets)
    mt = mesh.meshtags(domain, domain.topology.dim - 1, facets[arg_sort], np.full(len(facets), 2, dtype=np.int32))
    with common.Timer("~~Periodic: Compute mpc condition"):
        mpc = dolfinx_mpc.MultiPointConstraint(V)
        mpc.create_periodic_constraint_topological(V, mt, 2, periodic_relation, bcs)
        mpc.finalize()

    # Variational form
    u = ufl.TrialFunction(V)
    w = ufl.TestFunction(V)
    x = ufl.SpatialCoordinate(domain)
    a = ufl.dot(ufl.grad(w), ufl.grad(u))*ufl.dx
    f = -4*pi**2/L**2 * ufl.sin(2*pi*x[1]/L)
    L_ = w*-f*ufl.dx

    # Solve
    problem = dolfinx_mpc.LinearProblem(a, L_, mpc, bcs, petsc_options={"ksp_type": "preonly", "pc_type": "lu"})
    uh = problem.solve()
    uh.name = "u"

    # Output
    with io.XDMFFile(comm, "results_periodic_{}.xdmf".format(i), 'w') as file_results:
        file_results.write_mesh(domain)
        file_results.write_function(uh,1)

    # Evaluate norm
    error = fem.form((uh - u_exact(ufl)(x))**2 * ufl.dx)
    l2_error = np.sqrt(comm.allreduce(fem.assemble_scalar(error), MPI.SUM))
    print(l2_error)

# Solve for different mesh sizes
meshsize_1 = [L/10, L/10]
meshsize_2 = [L/10, L/30]
for i, (h1, h2) in enumerate(zip(meshsize_1, meshsize_2)):
    solve_periodic(L, h1, h2, i)

On the periodic mesh, the solution is periodic as expected. However, for the non-periodic mesh shown below, there is some error:

mesh_2
mesh_2801×644 132 KB

plot_2
plot_2801×644 9.12 KB

The plot shows the solution along the lower (red) and upper (blue) periodic boundaries, with large deviation between the solutions that grows worse near the ends where the nodes don’t overlap. Any suggestions on how to properly enforce this multi-point constraint using dolfinx_mpc? Thank you in advance.

2

You will have a problem when the mesh sizes varies as much as they do in your case, as the constraints are enforced per dof coordinate (vertex in your case).

There is no way that one can exactly enforce the relation of the bottom left corner and bottom right corner with the same constraint, as the a coarse element restricted to a set of multiple finer elements is not possible to represent in the coarse finite element space.