Periodic boundary conditions not being respected with mixed functionspace

[dthillaithevan]

Hi,

I am attempting to solve a doubly periodic stokes flow problem (adapted from this post) but am running into issues with enforcing PBCs for all function subspaces. For a simple stokes flow problem it looks like PBCs are being respected for the $y$-velocity and pressure but not for $x$-velocity. I have tried to construct this problem by setting up the $x$ and $y$ velocity functionspaces individually and run into the same problem.

I do not know if I have applied the PBC constraints incorrectly or if this a bug, please see MWE below and gmsh file attached (zip file as I cannot upload .msh directly - mesh.msh.zip).

Thanks!

import dolfinx
import dolfinx_mpc
from dolfinx.fem import FunctionSpace, Function, Constant
from dolfinx.io import XDMFFile
from dolfinx.mesh import locate_entities_boundary, meshtags
from dolfinx.io.gmshio import read_from_msh

from ufl import grad, div, inner, Measure, TrialFunctions, TestFunctions, FiniteElement, VectorElement
import ufl
from petsc4py import PETSc
import numpy as np
from mpi4py import MPI

def create_pbc(mesh, V, bcs):

    # Taken from demo 
    fdim = mesh.topology.dim - 1
    
    ## periodic boundary conditions
    pbc_directions = []
    pbc_slave_tags = []
    pbc_is_slave = []
    pbc_is_master = []
    pbc_meshtags = []
    pbc_slave_to_master_maps = []

    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] - L
            return out_x
        return pbc_slave_to_master_map

    def generate_pbc_is_slave(i):
        return lambda x: np.isclose(x[i], L)

    def generate_pbc_is_master(i):
        return lambda x: np.isclose(x[i], 0.0)

    for i in range(gdim):

        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(mesh, fdim, pbc_is_slave[-1])
        arg_sort = np.argsort(facets)
        pbc_meshtags.append(meshtags(mesh,
                                     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
                print(out_x[pbc_directions[i]][idx])
                print(idx)
                return out_x
        else:
            pbc_slave_to_master_map = pbc_slave_to_master_maps[i]

        for ij in range(gdim):
            mpc.create_periodic_constraint_topological(V.sub(0).sub(ij), pbc_meshtags[i],
                                                       pbc_slave_tags[i],
                                                       pbc_slave_to_master_map,
                                                       bcs)

        mpc.create_periodic_constraint_topological(V.sub(1), pbc_meshtags[i],
                                                       pbc_slave_tags[i],
                                                       pbc_slave_to_master_map,
                                                       bcs)

    # 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] - mesh.geometry.x[:, 0].max()
        out_x[1] = x[1] - mesh.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

    for ij in range(gdim):
        mpc.create_periodic_constraint_topological(V.sub(0).sub(ij), pbc_meshtags[1],
                                                   pbc_slave_tags[1],
                                                   pbc_slave_to_master_map,
                                                   bcs)

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

    mpc.finalize()

    
    return mpc


# Fluid domain is outside circle with radius r centered on [0.5, 0.5]
def inside_fluid_domain(x):
    return (x[0] + 0.5)**2 + (x[1]+0.5)**2 >= r

def corner_nodes(x):
    return np.logical_or(
        np.logical_and(
            np.isclose(x[0], 0.),
            np.logical_or(
                np.isclose(x[1], .0),
                np.isclose(x[1], L),
            )
        ),
        np.logical_and(
            np.isclose(x[0], L),
            np.logical_or(
                np.isclose(x[1], 0.),
                np.isclose(x[1], L),
            )
        ),
    )


# MPI initialization
comm = MPI.COMM_WORLD
rank = MPI.COMM_WORLD.rank

# =============================================================================
# Inputs
# =============================================================================
# Geometry dimension
# Geometry = square - circle of radius=r centered at [0.5, 0.5]
gdim = 2
L, H = 1, 1

# Radius of cylinder
r = 0.3

# Load mesh
mesh, ct, ft = read_from_msh("mesh.msh", MPI.COMM_WORLD, 0, gdim=gdim)

solid_domain_marker = 2
fluid_domain_marker = 1
interface_marker = 3

# Integration variable
dx = Measure('dx', domain=mesh, subdomain_data=ct)#,submesh = create_submesh(mesh_1, edim, entities)[0]
num_cells = mesh.topology.index_map(mesh.topology.dim).size_global

# Taylor-Hood elements and Mixed Functionspace
P1 = FiniteElement("CG", mesh.ufl_cell(), 1)
P2 = VectorElement("CG", mesh.ufl_cell(), 2)
TH = ufl.MixedElement([P2, P1])
V = FunctionSpace(mesh, TH)

# define function, trial function and test function
(ve, pr) = TrialFunctions(V)
(dve, dpr) = TestFunctions(V)

# =============================================================================
# Boundary conditions
# =============================================================================
## Set no slip condition for boundary of where fluid domain meets circle (solid)
# create collapsed functionspace
V_coll0, _ = V.sub(0).collapse()
v_bc = Function(V_coll0)
v_bc.x.set(0)
v_bc.x.scatter_forward()

# Create no slip BC at interface
dofs_v_interf = dolfinx.fem.locate_dofs_topological((V.sub(0), V_coll0)  , ft.dim, ft.find(interface_marker))  # interface tag: 9
bcs_v_interf = dolfinx.fem.dirichletbc(v_bc, dofs_v_interf, V.sub(0))

##P = 0 at corners since we have PBC
boundary_facets = locate_entities_boundary(mesh, 0, corner_nodes) 
dofs_corners = dolfinx.fem.locate_dofs_topological(V.sub(1), 0, boundary_facets)
bc_corners = dolfinx.fem.dirichletbc(PETSc.ScalarType(0.0), dofs_corners, V.sub(1))

bcs = [bcs_v_interf, bc_corners]

# Create multipoint constraint
mpc = create_pbc(mesh, V, bcs)

# =============================================================================
# # Variational problem
# =============================================================================
f_constant = Constant(mesh, PETSc.ScalarType((1,0)))

# Stokes weak variational form
a = inner(grad(ve), grad(dve)) * dx + pr * div(dve) * dx + div(ve) * dpr * dx
L = inner(f_constant, dve) * dx

# Solve linear problem
solution = Function(V)
problem = dolfinx_mpc.LinearProblem(a, L, mpc, bcs=bcs, petsc_options={
    "ksp_type": "preonly", "pc_type": "lu", "pc_factor_mat_solver_type": "mumps"})  # define linear problem
solution = problem.solve()

# Apply PBC back substitution
mpc.backsubstitution(solution)

ve, pr = solution.sub(0).collapse(), solution.sub(1).collapse()

ve.name = 'Velocity'
pr.name = 'Pressure'

# Save output to XDMF
with XDMFFile(mesh.comm, "output.xdmf", "w") as xdmf:
    xdmf.write_mesh(mesh)
    xdmf.write_function(ve)
    xdmf.write_function(pr)

Pressure:

$x$-velocity

$y$-velocity