Interior facets indexing on a GMSH mesh: unexpected results when solving in parallel

dolfinx
, , ,
1

Hi everyone,

I’m having an issue when trying to solve in parallel the Poisson equation with DG elements using a GMSH mesh.

Interior facet tags are required to write the DG formulation.
Unfortunatly, the gmshio.model_to_mesh function only returns exterior facets.
Here is my attempt to index interior facets (script below).

It works fine in serial computing but not in parallel (screenshots below).

Serial: expected results

poisson_serial
poisson_serial650×639 22.7 KB

Parallel: unexpected results

poisson_parallel
poisson_parallel660×653 30.2 KB

Complete script:


import gmsh
####################################################################
# Create unit square GMSH mesh
gmsh.initialize()
gdim = 2
Lx, Ly = 1.0, 1.0
rectangle = gmsh.model.occ.addRectangle(
    0.0, 0.0, 0.0, Lx, Ly
)
gmsh.model.occ.synchronize()

# Cell tag = 10
gmsh.model.addPhysicalGroup(gdim, [rectangle], tag=10)
gmsh.model.setPhysicalName(2, 10, "all")

# Boundary facets tag = 1
facet_entities = [entity[1] for entity in gmsh.model.getEntities(1)]
gmsh.model.addPhysicalGroup(1, facet_entities, tag=1)
gmsh.model.setPhysicalName(1, 1, 'boundary')

gmsh.option.setNumber("Mesh.CharacteristicLengthMin", 0.05)
gmsh.option.setNumber("Mesh.CharacteristicLengthMax", 0.05)
gmsh.model.mesh.generate(gdim)

####################################################################
import numpy as np
from dolfinx import mesh, fem
import dolfinx.fem.petsc as dfp
from dolfinx.io import gmshio, VTXWriter
import ufl
import basix.ufl
from petsc4py import PETSc
from mpi4py import MPI

####################################################################
# Convert GMSH mesh to dolfinx mesh
comm = MPI.COMM_WORLD
model_rank = 0

msh, cell_tags, boundary_facet_tags = gmshio.model_to_mesh(
    gmsh.model,
    comm=comm, rank=model_rank, gdim=gdim,
)

####################################################################
# Create meshtags including interior facets
tdim = msh.topology.dim
fdim = tdim - 1
msh.topology.create_connectivity(fdim, tdim)
facet_imap = msh.topology.index_map(tdim - 1)

num_facets_local = facet_imap.size_local
num_facets_ghost = facet_imap.num_ghosts

num_facets_local = facet_imap.size_local
num_facets_ghost = facet_imap.num_ghosts
num_facets = num_facets_local + num_facets_ghost

indices = np.arange(0, num_facets)
values = np.zeros(indices.shape, dtype=np.intc)  # all facets are tagged with zero
values[boundary_facet_tags.indices] = 1          # boundary facets tag = 1

# New facet tags including interior facets
facet_tags = mesh.meshtags(msh, tdim - 1, indices, values)

####################################################################
# Variational form

# Measure of the domain
dx = ufl.Measure('dx', domain=msh, subdomain_data=cell_tags, subdomain_id=(10))
## exterior facets
ds = ufl.Measure("ds", domain=msh, subdomain_data=facet_tags, subdomain_id=(1))
## interior facets
dS = ufl.Measure("dS", domain=msh, subdomain_data=facet_tags, subdomain_id=(0))

# Discontinuous Lagrange functionspace
V = fem.functionspace(
    msh,
    basix.ufl.element("DG", msh.topology.cell_name(), 1)
)

u = ufl.TrialFunction(V)
v = ufl.TestFunction(V)

# Function to store the solution
uh = fem.Function(V)

# Define mesh cell diameter and facet normal
h = ufl.CellDiameter(msh)
n = ufl.FacetNormal(msh)

# Bilinear form
LHS = ufl.inner(ufl.grad(v), ufl.grad(u)) * dx
LHS += - ufl.inner(ufl.avg(ufl.grad(v)), ufl.jump(u, n)) * dS
LHS += - ufl.inner(ufl.jump(v, n), ufl.avg(ufl.grad(u))) * dS

alpha = 1e3
LHS += alpha/ufl.avg(h) * ufl.inner(ufl.jump(v, n), ufl.jump(u, n)) * dS

# Linear form
f = fem.Constant(msh, PETSc.ScalarType(-6.0))
RHS = v * f * dx 

# Dirichlet BC (weak imposition)
uD = fem.Function(V)
uD.interpolate((lambda x: 1 + x[0]**2 + 2 * x[1]**2))
gamma_BC = 1e3
LHS += gamma_BC/h * u * v * ds
RHS += gamma_BC/h * uD * v * ds

# Forms
BF = fem.form(LHS)
LF = fem.form(RHS)

####################################################################
# Compute the solution
A = dfp.create_matrix(BF)
b = dfp.create_vector(LF)

solver = PETSc.KSP().create(comm)
solver.setOperators(A)
solver.setType('preonly')
solver.getPC().setType('lu')

# Assemble LHS
A.zeroEntries()
dfp.assemble_matrix(A, BF)
A.assemble()

# Assemble RHS
with b.localForm() as loc_b:
    loc_b.set(0)
dfp.assemble_vector(b, LF)

# Solve
solver.solve(b, uh.vector)
uh.x.scatter_forward()

solver.destroy()
A.destroy()
b.destroy()

####################################################################
# Export the solution
## Create continuous Lagrange functionspace for visualization
__V = fem.functionspace(
    msh,
    basix.ufl.element("CG", msh.topology.cell_name(), 1)
)
__u = fem.Function(__V)
__u.interpolate(uh)

## Write the solution
with VTXWriter(comm, "res_u.bp", __u, engine="BP4") as file:
    file.write(0.0)
    file.close()

conda environment
fenicsx: 0.8.0
gmsh: 4.12.2
petsc: 3.21.4
mpi: 3.1.6

2

Can you try adding a custom partitioner argument to that call, where the partitioner is created with partitioner = dolfinx.mesh.create_cell_partitioner(dolfinx.mesh.GhostMode.shared_facet)?
dS requires a non default ghost mode when running in parallel.

1 Like
3

Now it works! Thank you for the swift reply