Interpolation from submesh to mesh error by parallel computing?

dolfinx
,
1

Hello everyone,

I am trying to interpolate the result from submesh 1 to mesh 0 by the method from test_interpolation. When I run the code by a single core, like python3 test.py, everything is good. However, if I try it on multiple cores, like OMP_NUM_THREADS=1 mpirun -np 64 python3 test.py, I found the result on the boundary is not perfectly right as shown in the picture below. We have an original mesh 0 and it’s splitted into submesh 1 and submesh 2. Between submesh 1 and submesh 2, it should be a straight line, not like this. Can anyone help me with this problem? I want the result on the boundary to be perfectly right as a straight boundary line. This is important because I want to use the result (a function) of submesh 1 to be the boundary condition of submesh 2 (submesh 1 to mesh 0 to submesh 2). Thanks a lot for your time and help!

interpolation
interpolation1359×1343 75.9 KB

The code is attached here. Parallel interpolation from mesh to submesh is good. Parallel interpolation from submesh to mesh is not good, which is the exported function “r1” in the code.

import random

from mpi4py import MPI

import numpy as np

import basix
import ufl
from basix.ufl import blocked_element, custom_element, element, enriched_element, mixed_element
from dolfinx import default_real_type, default_scalar_type
from dolfinx.fem import (
    Expression,
    Function,
    assemble_scalar,
    create_interpolation_data,
    form,
    functionspace,
)
from dolfinx.geometry import bb_tree, compute_collisions_points
from dolfinx.mesh import (
    CellType,
    create_mesh,
    create_rectangle,
    create_submesh,
    create_unit_cube,
    create_unit_square,
    locate_entities,
    locate_entities_boundary,
    meshtags,
)

# OMP_NUM_THREADS=1 mpirun -np 64 python3

"""Test interpolation of a function between a sub-mesh and its parent mesh."""
mesh = create_unit_square(MPI.COMM_WORLD, 200, 300)

def left_locator(x):
    return x[0] <= 0.5 + 1e-14

def ref_func(x):
    return x[0] + x[1] ** 2

tdim = mesh.topology.dim
cells = locate_entities(mesh, tdim, left_locator)
submesh, entity_map, _, _ = create_submesh(mesh, tdim, cells)

u0 = Function(functionspace(mesh, ("Lagrange", 2)))
u0.interpolate(ref_func)

V1 = functionspace(submesh, ("DG", 3))
u1 = Function(V1)

smsh_cell_imap = submesh.topology.index_map(tdim)
smsh_cells = np.arange(smsh_cell_imap.size_local + smsh_cell_imap.num_ghosts)
parent_cells = entity_map.sub_topology_to_topology(smsh_cells, inverse=False)

# Interpolate u0 (defined on 'full' mesh) into u0 (defined on
# 'sub'0mesh)
u1.interpolate(u0, cells0=parent_cells, cells1=smsh_cells)

u1_exact = Function(V1)
u1_exact.interpolate(ref_func)
atol = 5 * np.finfo(default_scalar_type).resolution
np.testing.assert_allclose(u1_exact.x.array, u1.x.array, atol=atol)

# Map from sub to parent
W = functionspace(mesh, ("DG", 4))
w = Function(W)

# Interpolate Function defined on sub-mesh (u1_exact) to the part of
# a Function on the full mesh (w)
w.interpolate(u1_exact, cells0=np.arange(len(parent_cells)), cells1=parent_cells)
w_exact = Function(W)
w_exact.interpolate(ref_func, cells0=cells)
np.testing.assert_allclose(w.x.array, w_exact.x.array, atol=atol)


from dolfinx import io

xdmf_t1 = io.XDMFFile(mesh.comm, "./test_result/result_t1.xdmf", "w")
xdmf_t1.write_mesh(mesh)

xdmf_r1 = io.XDMFFile(mesh.comm, "./test_result/result_r1.xdmf", "w")
xdmf_r1.write_mesh(mesh)

xdmf_t2 = io.XDMFFile(submesh.comm, "./test_result/result_t2.xdmf", "w")
xdmf_t2.write_mesh(submesh)

xdmf_r2 = io.XDMFFile(submesh.comm, "./test_result/result_r2.xdmf", "w")
xdmf_r2.write_mesh(submesh)

t1 = Function(functionspace(mesh, ("Lagrange", 1)))
r1 = Function(functionspace(mesh, ("Lagrange", 1)))
t2 = Function(functionspace(submesh, ("Lagrange", 1)))
r2 = Function(functionspace(submesh, ("Lagrange", 1)))

t1.interpolate(u0)
r1.interpolate(w_exact)
t2.interpolate(u1)
r2.interpolate(u1_exact)

xdmf_t1.write_function(t1)
xdmf_r1.write_function(r1)
xdmf_t2.write_function(t2)
xdmf_r2.write_function(r2)

xdmf_t1.close()
xdmf_r1.close()
xdmf_t2.close()
xdmf_r2.close()
2

It seems like using interpolate_nonmatching can solve this problem. I will do more tests. If it really works, I will update the code.

3 
padding = 1e-14
cell_map0 = mesh.topology.index_map(mesh.topology.dim)
num_cells_on_proc = cell_map0.size_local + cell_map0.num_ghosts
cells0 = np.arange(num_cells_on_proc, dtype=np.int32)
interpolation_data1 = create_interpolation_data(W, V1, cells0, padding=padding)
w.interpolate_nonmatching(u1_exact, cells0, interpolation_data1)
w.x.scatter_forward()
w_exact = Function(W)
w_exact.interpolate(ref_func, cells0=cells)
w_exact.x.scatter_forward()
np.testing.assert_allclose(w.x.array, w_exact.x.array, atol=atol)
4

There are two things you should notice:

  1. When you interpolate on subsets of cells, such as one should always call u1.x.scatter_forward() afterwards, as you might need to update ghost degrees of freedom.
  2. You are going from a discontinuous space (W) to outputting on a continuous space t1,r1.
    This is not a well defined operation.

Thus, for all these exports, you should use a format that supportr discontinuous finite elements (VTKFile or VTXWriter, and call


from mpi4py import MPI

import numpy as np

from dolfinx import default_scalar_type
from dolfinx.fem import (
    Function,
    functionspace,
)
from dolfinx.mesh import (
    create_submesh,
    create_unit_square,
    locate_entities,
)

"""Test interpolation of a function between a sub-mesh and its parent mesh."""
mesh = create_unit_square(MPI.COMM_WORLD, 200, 300)

def left_locator(x):
    return x[0] <= 0.5 + 1e-14

def ref_func(x):
    return x[0] + x[1] ** 2

tdim = mesh.topology.dim
cells = locate_entities(mesh, tdim, left_locator)
submesh, entity_map, _, _ = create_submesh(mesh, tdim, cells)

u0 = Function(functionspace(mesh, ("Lagrange", 2)), name="u0")
u0.interpolate(ref_func)

V1 = functionspace(submesh, ("DG", 3))
u1 = Function(V1, name="u1")
    
smsh_cell_imap = submesh.topology.index_map(tdim)
smsh_cells = np.arange(smsh_cell_imap.size_local + smsh_cell_imap.num_ghosts)
parent_cells = entity_map.sub_topology_to_topology(smsh_cells, inverse=False)

# Interpolate u0 (defined on 'full' mesh) into u1 (defined on
# 'sub' mesh)
u1.interpolate(u0, cells0=parent_cells, cells1=smsh_cells)
u1.x.scatter_forward()

u1_exact = Function(V1, name="u1_exact")
u1_exact.interpolate(ref_func)

atol = 5 * np.finfo(default_scalar_type).resolution
np.testing.assert_allclose(u1_exact.x.array, u1.x.array, atol=atol)

# Map from sub to parent
W = functionspace(mesh, ("DG", 4))
w = Function(W, name="w")

# Interpolate Function defined on sub-mesh (u1_exact) to the part of
# a Function on the full mesh (w)
w.interpolate(u1_exact, cells0=np.arange(len(parent_cells)), cells1=parent_cells)
w_exact = Function(W, name="w_exact")
w_exact.interpolate(ref_func, cells0=cells)
np.testing.assert_allclose(w.x.array, w_exact.x.array, atol=atol)


from dolfinx import io

with io.VTXWriter(mesh.comm, "results.bp", [w, w_exact]) as bp:
    bp.write(0.0)

yielding:

image
image1354×1445 58 KB

1 Like
6

That’s great to have a better standing on the interpolation. Thanks a lot for your time and help!