Periodic BCs in parallel not working

mesh
1

I want to apply Periodic BC in all directions (2D). I followed the method given here.

This works well in serial, but not in parallel. Below is how it looks in parallel (2 procs).

image
image852×855 94.1 KB
.

Here is the MWE:

# Import Libraries
from mpi4py import MPI
from dolfinx import mesh, fem, io, log, default_real_type
import numpy as np
import ufl
from basix.ufl import element, mixed_element
from dolfinx.fem.petsc import NonlinearProblem
from dolfinx.nls.petsc import NewtonSolver
from petsc4py import PETSc
from pathlib import Path
from script import create_periodic_mesh

# Create mesh
domain = mesh.create_unit_square(MPI.COMM_WORLD, 100, 100, cell_type=mesh.CellType.triangle)

# Apply Periodic BCs
# Extract the dimension of the mesh
L_min = [domain.comm.allreduce(np.min(domain.geometry.x[:, i]), op=MPI.MIN) for i in range(3)]
L_max = [domain.comm.allreduce(np.max(domain.geometry.x[:, i]), op=MPI.MAX) for i in range(3)]
# Define the periodic boundary condition
def i_x(x):
    return np.isclose(x[0], L_min[0])

def i_y(x):
    return np.isclose(x[1], L_min[1])

def indicator(x):
    return i_x(x) | i_y(x)

def mapping(x):
    values = x.copy()
    values[0] += i_x(x) * (L_max[0] - L_min[0])
    values[1] += i_y(x) * (L_max[1] - L_min[1])
    return values

# domain, replaced_vertices, replacement_map = create_periodic_mesh(domain, indicator, mapping)
fdim = domain.topology.dim - 1
domain.topology.create_entities(fdim)
domain, replaced_vertices, replacement_map = create_periodic_mesh(domain, indicator, mapping)
domain.topology.create_entities(fdim)
# Locate facets for boundary conditions and create meshtags
domain.topology.create_connectivity(fdim, fdim + 1)

# Create FunctionSpace
V = fem.functionspace(domain, ("Lagrange", 1))

# Define variational problem
# Define trial and test functions
u = fem.Function(V)

# Initial condition
u.x.array[:] = 0.0
rng = np.random.default_rng(2)
u.interpolate(lambda x: 0.63 + 0.1 * (0.5 - rng.random(x.shape[1])))
u.x.scatter_forward()

results_folder = Path("results")
results_folder.mkdir(exist_ok=True, parents=True)
filename = results_folder / "out_mwe.bp"
with io.VTXWriter(domain.comm, filename, [u], engine="BP4") as vtx:
    vtx.write(0)
2

This is due to the fact that the mesh is actually periodic and you have to use a different method for visualization (and setting initial conditions).
Here is an mwe that is adjusted to illustrate that it works:

# Import Libraries
from mpi4py import MPI
from dolfinx import mesh, fem, io, la
import dolfinx.fem.petsc
import numpy as np

from pathlib import Path
from script import create_periodic_mesh

# Create mesh
domain = mesh.create_unit_square(
    MPI.COMM_WORLD,
    100,
    100,
    cell_type=mesh.CellType.triangle,
    ghost_mode=mesh.GhostMode.shared_facet,
)

# Apply Periodic BCs
# Extract the dimension of the mesh
L_min = [
    domain.comm.allreduce(np.min(domain.geometry.x[:, i]), op=MPI.MIN) for i in range(3)
]
L_max = [
    domain.comm.allreduce(np.max(domain.geometry.x[:, i]), op=MPI.MAX) for i in range(3)
]


# Define the periodic boundary condition
def i_x(x):
    return np.isclose(x[0], L_min[0], atol=1e-8)


def i_y(x):
    return np.isclose(x[1], L_min[1], atol=1e-8)


def indicator(x):
    return i_x(x) | i_y(x)


def mapping(x):
    values = x.copy()
    values[0] += i_x(x) * (L_max[0] - L_min[0])
    values[1] += i_y(x) * (L_max[1] - L_min[1])
    return values


# domain, replaced_vertices, replacement_map = create_periodic_mesh(domain, indicator, mapping)
fdim = domain.topology.dim - 1
domain.topology.create_entities(fdim)
domain, replaced_vertices, replacement_map = create_periodic_mesh(
    domain, indicator, mapping
)
domain.topology.create_entities(fdim)
# Locate facets for boundary conditions and create meshtags
domain.topology.create_connectivity(fdim, fdim + 1)

print(domain.topology.index_map(2).size_global)

# Create FunctionSpace
V = fem.functionspace(domain, ("Lagrange", 1))

# Define variational problem
# Define trial and test functions
uh = fem.Function(V)

import ufl

u = ufl.TrialFunction(V)
v = ufl.TestFunction(V)
a = ufl.inner(u, v) * ufl.dx
x = ufl.SpatialCoordinate(domain)
L = ufl.inner(x[0] + 2 * x[1] ** 2, v) * ufl.dx
fem.petsc.LinearProblem(
    a,
    L,
    u=uh,
    bcs=[],
    petsc_options={
        "ksp_type": "preonly",
        "pc_type": "lu",
        "pc_factor_mat_solver_type": "mumps",
        "ksp_error_if_not_converged": True,
    },
).solve()


print(la.norm(uh.x))


V_out = fem.functionspace(domain, ("DG", 1))
v_out = fem.Function(V_out)
v_out.interpolate(uh)
v_out.x.scatter_forward()


results_folder = Path("results")
results_folder.mkdir(exist_ok=True, parents=True)
filename = results_folder / "out_mwe.bp"
with io.VTXWriter(domain.comm, filename, [v_out], engine="BP4") as vtx:
    vtx.write(0)

First of all, using random noise in serial and parallel will give different results (as the length of the arrays of x coordinates are different).
Secondly, using interpolation on a periodic mesh can be troublesome:

  • If you have a spatial dependency in the interpolation, it will depend on what process owns the dof (and what is the last cell on the process that one interpolated with respect to).
    Finally, one cannot visualize the result with a standard Lagrange space. If you compare the output of your code (in serial) when you use a periodic mesh and when you comment out that command, you will observe that the domain in the bp-file is missing a set of cells.
    This is because you need to interpolate it into an appropriate DG space to capture the geometry on both sides of the domain.
1 Like