Different results of area of a deformed mesh in parallel

[chenyongxin]

I deformed a mesh by appending a deformation array to the mesh with mesh.geometry.x[:, d] += deformation_array.x.array[:], where d is the dimension of mesh to be moved. I got the correct result with a serial run. The result was wrong with parallel runs but the visualization was correct.

To brief the problem: I have two non-matching meshes in a FSI problem. Every time I need to interpolate the fields from the background (bigger) domain to the foreground (smaller) domain and then form the deformation array. In the following code, serial run returns the expected result 1, but it gives 1.30419921875 with
mpirun -n 4 python parallel_area.py
The doflinx version is 0.6.0.

MWE:

# parallel_area.py
# Test to measure area of mesh after deformation in parallel 

from mpi4py import MPI
from dolfinx import fem, io, mesh
from ufl import dx

#  Dimension of the problem
ndims = 2

#  Foreground mesh and function space
mesh_square = mesh.create_unit_square(MPI.COMM_WORLD, 32, 32)
mesh_square.geometry.x[:, :2] -= 0.5                  # move center to (0, 0) with size (-0.5, -0.5) -> (0.5, 0.5)
V_square = fem.VectorFunctionSpace(mesh_square, ("Lagrange", 1))
Q_square =       fem.FunctionSpace(mesh_square, ("Lagrange", 1))

#  Background mesh and function space
mesh_background = mesh.create_unit_square(MPI.COMM_WORLD, 64, 64)
mesh_background.geometry.x[:, :2] -= 0.5              # move center to (0, 0)
mesh_background.geometry.x[:, :2] *= 2                # double the size: (-1, -1) -> (1, 1)
V_background = fem.VectorFunctionSpace(mesh_background, ("Lagrange", 1))
Q_background =       fem.FunctionSpace(mesh_background, ("Lagrange", 1))

#  Original area
one = fem.Function(Q_square)
one.x.array[:] = 1
f = one*dx
area0 = MPI.COMM_WORLD.allreduce(fem.assemble_scalar(fem.form(one*dx)), op = MPI.SUM)
print("Original area:", area0)                        # always returns 1 as expected
with io.XDMFFile(MPI.COMM_WORLD, "original_area.xdmf", "w") as file:
    file.write_mesh(mesh_square)


## Deformation occurs
#  Define background fields
u_background = fem.Function(V_background)
u_background.sub(0).interpolate(lambda x: x[1])       # like a shear flow

#  Interpolate to foreground fields
u_square = fem.Function(V_square)
for d in range(ndims):
    u_square.sub(d).interpolate(u_background.sub(d))

#  Deform the foreground mesh
dt = 1                                                # an coefficient which mimics the one in the real problem
for d in range(ndims):
    mesh_square.geometry.x[:, d] += u_square.sub(d).collapse().x.array[:] * dt

# Compute new area
one = fem.Function(Q_square)
one.x.array[:] = 1
f = one*dx
area1 = MPI.COMM_WORLD.allreduce(fem.assemble_scalar(fem.form(one*dx)), op = MPI.SUM)
with io.XDMFFile(MPI.COMM_WORLD, "changed_area.xdmf", "w") as file:
    file.write_mesh(mesh_square)

print("Changed area:", area1)                         # mpirun -n 1: 1; -n 2: 1; -n 3: 1.45068359375; -n 4: 1.30419921875

[dokken]

Running:

with io.VTXWriter(MPI.COMM_SELF, f"changed_area_{MPI.COMM_WORLD.rank}.bp", mesh_square, engine="BP4") as file:
    file.write(0.0)

you can see that the mapping actually isnt correct in parallel.

To fix this consider:

# parallel_area.py
# Test to measure area of mesh after deformation in parallel

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

#  Dimension of the problem
ndims = 2

#  Foreground mesh and function space
mesh_square = mesh.create_unit_square(MPI.COMM_WORLD, 32, 32)
# move center to (0, 0) with size (-0.5, -0.5) -> (0.5, 0.5)
mesh_square.geometry.x[:, :2] -= 0.5
V_square = fem.VectorFunctionSpace(mesh_square, ("Lagrange", 1))
Q_square = fem.FunctionSpace(mesh_square, ("Lagrange", 1))

#  Background mesh and function space
mesh_background = mesh.create_unit_square(MPI.COMM_WORLD, 64, 64)
mesh_background.geometry.x[:, :2] -= 0.5              # move center to (0, 0)
# double the size: (-1, -1) -> (1, 1)
mesh_background.geometry.x[:, :2] *= 2
V_background = fem.VectorFunctionSpace(mesh_background, ("Lagrange", 1))
Q_background = fem.FunctionSpace(mesh_background, ("Lagrange", 1))

#  Original area
one = fem.Function(Q_square)
one.x.array[:] = 1
f = one*dx
area0 = MPI.COMM_WORLD.allreduce(
    fem.assemble_scalar(fem.form(one*dx)), op=MPI.SUM)
# always returns 1 as expected
print("Original area:", area0)
with io.XDMFFile(MPI.COMM_WORLD, "original_area.xdmf", "w") as file:
    file.write_mesh(mesh_square)


# Deformation occurs
#  Define background fields
u_background = fem.Function(V_background)
u_background.sub(0).interpolate(lambda x: x[1])       # like a shear flow

#  Interpolate to foreground fields
u_square = fem.Function(V_square)
for d in range(ndims):
    nms = fem.create_nonmatching_meshes_interpolation_data(
        u_square.sub(d).function_space.mesh._cpp_object,
        u_square.sub(d).function_space.element,
        u_background.sub(d).function_space.mesh._cpp_object, padding=1e-6)
    u_square.sub(d).interpolate(
        u_background.sub(d), nmm_interpolation_data=nms)
u_square.x.scatter_forward()
#  Deform the foreground mesh
# an coefficient which mimics the one in the real problem
dt = 1
for d in range(ndims):
    V_sub, sub_to_parent = u_square.sub(d).function_space.collapse()
    sub_to_parent = np.asarray(sub_to_parent)
    sub_to_parent = sub_to_parent // ndims
    mesh_square.geometry.x[sub_to_parent, d] += u_square.sub(
        d).collapse().x.array * dt
# Compute new area
one = fem.Function(Q_square)
one.x.array[:] = 1
f = one*dx
area1 = MPI.COMM_WORLD.allreduce(
    fem.assemble_scalar(fem.form(one*dx)), op=MPI.SUM)
with io.VTXWriter(MPI.COMM_SELF, f"changed_area_{MPI.COMM_WORLD.rank}.bp", mesh_square, engine="BP4") as file:
    file.write(0.0)
print("Changed area:", area1)

which I’ve executed with DOLFINx v0.7.2
yielding:

for the 0th process.

[chenyongxin]

I see. I thought the >1 result is because I added the contribution of the ghost cells for serveral times in parallel. The real problem in my code is the wrong mapping between mesh and function. So I replace with the following part and I can get 1 now with DOLFINx 0.6.0:

@dokken

for d in range(ndims):
    V_sub, sub_to_parent = u_square.sub(d).function_space.collapse()
    sub_to_parent = np.asarray(sub_to_parent)
    sub_to_parent = sub_to_parent // ndims
    mesh_square.geometry.x[sub_to_parent, d] += u_square.sub(
        d).collapse().x.array * dt

Some follow-up questions:
(1) VTXWriter vs XDMFFile. Why XDMFFile still gives correct visualization even though the mapping is actually not correct in parallel but VTXWriter can dectect that?
(2) The result of area 1 with mpirun -n 1 or mpirun -n 2 in my original code is evaluated coincidentally？The wrong mapping should be independent of number of processes.
(3) I see that you use sub_to_parent = sub_to_parent // ndims is because I used a (Vector) Function whose DOFs are ordered. What if I use a MixedElement and FunctionSpace instead of VectorFunctionSpace here? The dofs are most usually ordered for a mixed element.

[dokken]

@chenyongxin

VTXWriter vs XDMFFile. Why XDMFFile still gives correct visualization even though the mapping is actually not correct in parallel but VTXWriter can dectect that?

There are subtle differences in the implementation of the two. Also note that your XDMFFile is using

@chenyongxin

with io.XDMFFile(MPI.COMM_WORLD, "changed_area.xdmf", "w") as file:

while VTXWriter is in this case using

@dokken

with io.VTXWriter(MPI.COMM_SELF, f"changed_area_{MPI.COMM_WORLD.rank}.bp", mesh_square, engine="BP4") as file:

writing data per process, which is possible due to how VTXWriter is created.

@chenyongxin

(2) The result of area 1 with mpirun -n 1 or mpirun -n 2 in my original code is evaluated coincidentally？The wrong mapping should be independent of number of processes.

The issue is with Ghost degrees of freedom. These depend on the number of processes.
Thus the mapping depend on the number of processes.

@chenyongxin

3. I see that you use sub_to_parent = sub_to_parent // ndims is because I used a (Vector) Function whose DOFs are ordered. What if I use a MixedElement and FunctionSpace instead of VectorFunctionSpace here?

Then things become more complicated. Please make a minimal reproducible example

[chenyongxin]

@dokken

@chenyongxin

3. I see that you use sub_to_parent = sub_to_parent // ndims is because I used a (Vector) Function whose DOFs are ordered. What if I use a MixedElement and FunctionSpace instead of VectorFunctionSpace here?

Then things become more complicated. Please make a minimal reproducible example

In the real applications, I use MixedElement instead of (Vector)Function. As the DOFs of a MixedElement are not always ordered, using the same way with a MixedElement gives wrong result as follows:

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

#  Dimension of the problem and mesh
ndims = 2
mesh_square = mesh.create_unit_square(MPI.COMM_WORLD, 32, 32)

#  Original area
area0 = MPI.COMM_WORLD.allreduce(fem.assemble_scalar(fem.form(1*ufl.dx(mesh_square))), 
                                 op = MPI.SUM)
print("Original area:", area0)                # always returns 1 as expected

## Deformation occurs: replace (Vector)Function(Space) to a MixedElement with P1-P1 element
FE       = ufl.FiniteElement("Lagrange", mesh_square.ufl_cell(), 1)
ME       = ufl.MixedElement([FE, FE])
V_ME     = fem.FunctionSpace(mesh_square, ME)
u_square = fem.Function(V_ME)
u_square.sub(0).interpolate(lambda x: x[1])   # like a shear flow

#  Deform the mesh with a deformation array
dt = 1              # mimics the coefficient in the real applications
for d in range(ndims):

    V_sub, sub_to_parent = u_square.sub(d).function_space.collapse()  
    sub_to_parent = np.asarray(sub_to_parent)
    sub_to_parent = sub_to_parent // ndims
    mesh_square.geometry.x[sub_to_parent, d] += u_square.sub(
        d).collapse().x.array * dt

# Compute new area
area1 = MPI.COMM_WORLD.allreduce(fem.assemble_scalar(fem.form(1*ufl.dx(mesh_square))), 
                                 op = MPI.SUM)
print("Changed area:", area1)       
# mpirun -n 1: 0.99951171875
# mpirun -n 2: 1.02880859375

My understanding: for mesh, the vertices are ordered with vertex_id: mesh.geometry.x[:, 0] → [vertex_0_x, vertex_1_x, vertex_2_x, ...]. For Function, function.x.array[:] → [dof_0, dof_1, dof_2, ...]. For a general usage, I guess here I need an equivalent vertex_to_dof functionality instead of sub_to_parent = sub_to_parent // ndims.

[dokken]

@chenyongxin

In the real applications, I use MixedElement instead of (Vector)Function. As the DOFs of a MixedElement are not always ordered, using the same way with a MixedElement gives wrong result as follows:

Is there a specific motivation for using a Mixed element? I get that if you use it, you end up with this issue, but why not use a vector-element?

[chenyongxin]

It is a fully-coupled FSI problem. All variables are solved simultaneously but not solved in a multi-stage decouple model, so I need to use MixedElement to include all variables (ux-uy-uz-p).

[dokken]

You could use a Mixed element of a vector element + pressure element. I think that would simplify the code alot.
This is similar to what is done in the stokes demo in dolfinx

[dokken]

Irregardless, ill see if I have time to code up a dof to vertex map for you tonight.

[chenyongxin]

Application of point forces, mapping vertex indices to corresponding DOFs

Consider the following code, that can map from any mesh node to any Lagrange space: import ufl import dolfinx from mpi4py import MPI import numpy as np import scipy def f(x): return x[0] + 3 * x[1] N = 4 P = 1 Q = 2 mesh = dolfinx.mesh.create_unit_cube(MPI.COMM_WORLD, N, N, N) el0 = ufl.VectorElement("Lagrange", mesh.ufl_cell(), P) el1 = ufl.FiniteElement("Lagrange", mesh.ufl_cell(), Q) el = ufl.MixedElement([el0, el1]) V = dolfinx.fem.FunctionSpace(mesh, el) u = dolfinx.fem.Function(V…

Based on the above your post, I still didn’t get it how to deal with the DOFs of a mixed element. Will be waiting for your code of map method here

I think I can work around the mixed element DOF problem with the following code, which is based on the facts:
(1) The DOFs of VectorFunctionSpace are ordered;
(2) The DOFs of a mixed element are not always ordered;
(3) The vertex-to-DOF mapping is only valid with P1 element.
So firstly I interpolate a higher order mixed element function (P2) to a vector function (P1) and then move mesh. Any comments and/or suggestions?

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

#  Dimension of the problem and mesh
ndims = 2
mesh_square = mesh.create_unit_square(MPI.COMM_WORLD, 32, 32)

#  Original area
area0 = MPI.COMM_WORLD.allreduce(fem.assemble_scalar(fem.form(1*ufl.dx(mesh_square))), 
                                op = MPI.SUM)
print("Original area:", area0)                # always returns 1 as expected

## Deformation occurs: replace (Vector)Function(Space) to a MixedElement with P1-P1 element
deg      = 2
FE       = ufl.FiniteElement("Lagrange", mesh_square.ufl_cell(), deg)
ME       = ufl.MixedElement([FE, FE])
V_ME     = fem.FunctionSpace(mesh_square, ME)
u_square = fem.Function(V_ME)
u_square.sub(0).interpolate(lambda x: x[1])   # like a shear flow

#  To interpolate to a Vector Function with P1 element as a deformation array
V = fem.VectorFunctionSpace(mesh_square, ("Lagrange", 1))
u = fem.Function(V)
for d in range(ndims):
   u.sub(d).interpolate(u_square.sub(d))

#  Deform the mesh with a deformation array
for d in range(ndims):

   V_sub, sub_to_parent = u.sub(d).function_space.collapse()  
   sub_to_parent = np.asarray(sub_to_parent)
   sub_to_parent = sub_to_parent // ndims
   mesh_square.geometry.x[sub_to_parent, d] += u.sub(d).collapse().x.array

# Compute new area
area1 = MPI.COMM_WORLD.allreduce(fem.assemble_scalar(fem.form(1*ufl.dx(mesh_square))), 
                                op = MPI.SUM)
print("Changed area:", area1)       

[dokken]

@chenyongxin

So firstly I interpolate a higher order mixed element function (P2) to a vector function (P1) and then move mesh. Any comments and/or suggestions?

As long as you work with function spaces that has a point evaluation dof at the vertices, you can skip the interpolation and do the following:

# parallel_area.py
# Test to measure area of mesh after deformation in parallel

from IPython import embed
import numpy as np
from mpi4py import MPI
from dolfinx import fem, io, mesh, cpp
import ufl

#  Dimension of the problem
ndims = 2

#  Foreground mesh and function space
mesh_square = mesh.create_unit_square(MPI.COMM_WORLD, 32, 32)
# move center to (0, 0) with size (-0.5, -0.5) -> (0.5, 0.5)
mesh_square.geometry.x[:, :2] -= 0.5
el = ufl.FiniteElement("Lagrange", mesh_square.ufl_cell(), 2)
mixed_el = ufl.MixedElement([el, el])
V_square = fem.FunctionSpace(mesh_square, mixed_el)

#  Background mesh and function space
mesh_background = mesh.create_unit_square(MPI.COMM_WORLD, 64, 64)
mesh_background.geometry.x[:, :2] -= 0.5              # move center to (0, 0)
# double the size: (-1, -1) -> (1, 1)
mesh_background.geometry.x[:, :2] *= 2
V_background = fem.FunctionSpace(mesh_background, mixed_el)

#  Original area
f = 1*ufl.dx(domain=mesh_square)
area0 = MPI.COMM_WORLD.allreduce(
    fem.assemble_scalar(fem.form(f)), op=MPI.SUM)

# always returns 1 as expected
print("Original area:", area0)


# Deformation occurs
#  Define background fields
u_background = fem.Function(V_background)
u_background.sub(0).interpolate(lambda x: x[1])       # like a shear flow

#  Interpolate to foreground fields
u_square = fem.Function(V_square)
for d in range(ndims):
    nms = fem.create_nonmatching_meshes_interpolation_data(
        u_square.sub(d).function_space.mesh._cpp_object,
        u_square.sub(d).function_space.element,
        u_background.sub(d).function_space.mesh._cpp_object, padding=1e-6)
    u_square.sub(d).interpolate(
        u_background.sub(d), nmm_interpolation_data=nms)
u_square.x.scatter_forward()


with io.VTXWriter(u_square.function_space.mesh.comm, f"u_square.bp", [u_square.sub(0).collapse()], engine="BP4") as file:
    file.write(0.0)
with io.VTXWriter(u_background.function_space.mesh.comm, f"u_background.bp", [u_background.sub(0).collapse()], engine="BP4") as file:
    file.write(0.0)

with io.VTXWriter(MPI.COMM_SELF, f"orginal_area_{MPI.COMM_WORLD.rank}.bp", mesh_square, engine="BP4") as file:
    file.write(0.0)

# Create a vertex to dof map for the mixed element
num_vertices = mesh_square.topology.index_map(
    0).size_local + mesh_square.topology.index_map(0).num_ghosts
tdim = mesh_square.topology.dim
num_cells = mesh_square.topology.index_map(
    tdim).size_local + mesh_square.topology.index_map(tdim).num_ghosts
vertex_to_dofmaps = np.zeros(
    (V_square.num_sub_spaces, num_vertices), dtype=np.int32)
num_vertices_per_cell = cpp.mesh.cell_num_vertices(
    mesh_square.topology.cell_types[0])
c_to_v = mesh_square.topology.connectivity(tdim, 0)
vertices_to_mesh_nodes = cpp.mesh.entities_to_geometry(
    mesh_square._cpp_object, 0, np.arange(num_vertices, dtype=np.int32), False).reshape(-1)
for d in range(vertex_to_dofmaps.shape[0]):
    dofmap = V_square.sub(d).dofmap
    local_vertex_map = np.zeros(num_vertices_per_cell, dtype=np.int32)
    # For given sub space, find the local vertex to dof map for a cell
    for v in range(num_vertices_per_cell):
        v_to_d = dofmap.dof_layout.entity_dofs(0, v)
        assert len(v_to_d) == 1
        local_vertex_map[v] = v_to_d[0]

    for cell in range(num_cells):
        vertices = c_to_v.links(cell)
        for v in range(num_vertices_per_cell):
            vertex_to_dofmaps[d, vertices[v]] = dofmap.cell_dofs(
                cell)[local_vertex_map[v]]

# Deform the foreground mesh
# an coefficient which mimics the one in the real problem


dt = 1
for d in range(ndims):
    mesh_square.geometry.x[vertices_to_mesh_nodes,
                           d] += u_square.x.array[vertex_to_dofmaps[d]]

# Compute new area
f = 1*ufl.dx(domain=mesh_square)
area1 = MPI.COMM_WORLD.allreduce(
    fem.assemble_scalar(fem.form(f)), op=MPI.SUM)
with io.VTXWriter(MPI.COMM_SELF, f"changed_area_{MPI.COMM_WORLD.rank}.bp", mesh_square, engine="BP4") as file:
    file.write(0.0)
print("Changed area:", area1)

Some neat things in this code snippet:

1. For a time dependent problem, you can create vertices_to_mesh_nodes and vertex_to_dofmaps once, and use them in a loop.
2. No interpolation or collapsing is required inside the mesh update loop,i.e.

for d in range(ndims):
    mesh_square.geometry.x[vertices_to_mesh_nodes,
                           d] += u_square.x.array[vertex_to_dofmaps[d]]

as we used the un-collapsed spaces to make the vertex to dof-map.