Incorrect parallel result of finding DOFs of a phase excepting interface in a two-phase problem

General
,
1

For a two-phase system, the domain can be categoried into three parts: phase1, phase2 and interface. I try to find DOFs of one phase excluding those on the interface for applying boundary conditions. The following my attempt gives different results in serial and parallel runs, where the serial run is correct. The result of a simplified problem is illustrated by the attached figure.

marked_result
marked_result2088×1146 263 KB

I first seperate cells by two phases and then get their vertices of each phase. The vertices on the interface can be obtained from the intersection. The vertices of one phase excluding those on the interface can be obtained by the difference between phase vertices and those on the interface. Finally, get corresponding DOFs by using fem.locate_dofs_topological as

fluid_vertices = mesh.compute_incident_entities(_mesh.topology, Omega_f, tdim, 0)
solid_vertices = mesh.compute_incident_entities(_mesh.topology, Omega_s, tdim, 0)
interface_vertices = np.intersect1d(fluid_vertices, solid_vertices)
Omega_s_inside_vertices = np.setdiff1d(solid_vertices, interface_vertices)

The full MWE with dolfinx v0.8.0:

from mpi4py import MPI
import numpy as np
from dolfinx import fem, mesh, io

comm = MPI.COMM_WORLD
_mesh = mesh.create_unit_square(comm, 50, 50)

tdim, fdim = 2, 1

def solid_domain(x):
    return x[1] < 0.5 + 1e-6
    
num_cells = (_mesh.topology.index_map(tdim).size_local + 
             _mesh.topology.index_map(tdim).num_ghosts)
Omega_s = mesh.locate_entities(_mesh, tdim, solid_domain)
Omega_f = np.setdiff1d(np.arange(num_cells), Omega_s)

fluid_vertices = mesh.compute_incident_entities(_mesh.topology, Omega_f, tdim, 0)
solid_vertices = mesh.compute_incident_entities(_mesh.topology, Omega_s, tdim, 0)
interface_vertices = np.intersect1d(fluid_vertices, solid_vertices)
Omega_s_inside_vertices = np.setdiff1d(solid_vertices, interface_vertices)

V = fem.functionspace(_mesh, ("Lagrange", 1))
f = fem.Function(V)
_mesh.topology.create_connectivity(0, tdim)
f.x.array[fem.locate_dofs_topological(V, 0, Omega_s_inside_vertices)] = 1
f.x.scatter_forward()

with io.XDMFFile(MPI.COMM_WORLD, f"test.xdmf", "w") as xdmf_s:
        xdmf_s.write_mesh(_mesh)
        xdmf_s.write_function(f, 0)
2

This is not so simple, as vertices might be ghosted on processes where they do not share cell information.

Some processes might have a facet on the surface, but not have a cell on either side of the surface.
You therefore need to communicate some information between these processes:

from mpi4py import MPI
import numpy as np
from dolfinx import fem, mesh, io, la

comm = MPI.COMM_WORLD
_mesh = mesh.create_unit_square(comm, 10, 10, ghost_mode=mesh.GhostMode.shared_facet)

tdim, fdim = 2, 1

def solid_domain(x):
    return x[1] <= 0.5 + 1e-6
 
num_cells = (_mesh.topology.index_map(tdim).size_local + 
             _mesh.topology.index_map(tdim).num_ghosts)
Omega_s = mesh.locate_entities(_mesh, tdim, solid_domain)

# Scatter ghosted_entities
def scatter_ghosts(index_map, local_entities):
    """
    Given a set of owned entities, distribute them to processes that has them as ghosts
    """      
    vector = la.vector(index_map, 1, dtype=np.int32)
    vector.array[local_entities] = 1
    vector.scatter_forward()
    return np.flatnonzero(vector.array).astype(np.int32)

def accumulate_ghosts(index_map, entities):
    """
    Given a set of entities (owned or ghosted), accumulate them on the owning process
    """      
    vector = la.vector(index_map, 1, dtype=np.int32)
    vector.array[:] = entities
    vector.scatter_reverse(la.InsertMode.add)
    vector.scatter_forward()
    return vector.array



Omega_s = scatter_ghosts(_mesh.topology.index_map(tdim), Omega_s)

Omega_f = np.setdiff1d(np.arange(num_cells), Omega_s)

fluid_vertices = mesh.compute_incident_entities(_mesh.topology, Omega_f, tdim, 0)
solid_vertices = mesh.compute_incident_entities(_mesh.topology, Omega_s, tdim, 0)

num_vertices_on_process = _mesh.topology.index_map(0).size_local + _mesh.topology.index_map(0).num_ghosts

fluid_vertex_marker = np.zeros(num_vertices_on_process, dtype=np.int32)
fluid_vertex_marker[fluid_vertices] = 1
fluid_vertex_marker = accumulate_ghosts(_mesh.topology.index_map(0), fluid_vertex_marker)

solid_vertex_marker = np.zeros(num_vertices_on_process, dtype=np.int32)
solid_vertex_marker[solid_vertices] = 1
solid_vertex_marker = accumulate_ghosts(_mesh.topology.index_map(0), solid_vertex_marker)




fluid_vertices = scatter_ghosts(_mesh.topology.index_map(0), np.flatnonzero(fluid_vertex_marker).astype(np.int32))
solid_vertices = scatter_ghosts(_mesh.topology.index_map(0), np.flatnonzero(solid_vertex_marker).astype(np.int32))


interface_vertices = np.intersect1d(fluid_vertices, solid_vertices)
Omega_s_inside_vertices = np.setdiff1d(solid_vertices, interface_vertices)

V = fem.functionspace(_mesh, ("Lagrange", 1))
f = fem.Function(V)
_mesh.topology.create_connectivity(0, tdim)
f.x.array[fem.locate_dofs_topological(V, 0, Omega_s_inside_vertices)] = 1
f.x.scatter_forward()

with io.XDMFFile(MPI.COMM_WORLD, f"test.xdmf", "w") as xdmf_s:
        xdmf_s.write_mesh(_mesh)
        xdmf_s.write_function(f, 0)
3

I see the communication need for vertices between processes. I am wondering if this is also necessary for cells? i.e. if the following line is necessary? I did not see difference of cell id before and after this.

4

Locate entities only returns owned entities

I don’t think it is needed, as all ghosts are found with locate entities (dolfinx/cpp/dolfinx/mesh/utils.h at 201be830056d20ff0bdd678dc95a82cd9c9ecc77 · FEniCS/dolfinx · GitHub). I was being overly cautious