Help with the dof maps between two meshes in mpi

1

Hi,

I have code that has two identical meshes. I move one of them and solve a problem on it. Now I need to map back the solution to the mesh that I didn’t move. I have a simple code attached below that works fine in series. But when I use a parallel run even with 2 processes it gives me the attached error. I think it is because the process break the two meshes differently. Could you please help me fix this code to work for the parallel case?
from fenics import *
import numpy as np

#Defining two identical meshes
mesh1 = UnitSquareMesh(20,20)  #we move this one
mesh2 = UnitSquareMesh(20,20)  #we keep this one as the reference mesh

#See what the 
File('test/ref_mesh.pvd')<<mesh1

#function spaces on the two meshes
V_1 = VectorFunctionSpace(mesh1,'P',1)
V_2 = VectorFunctionSpace(mesh2,'P',1)
W_1 = FunctionSpace(mesh1,'P',1)
W_2 = FunctionSpace(mesh2,'P',1)


#defining a vector to move mesh1 with
u = Expression(('sin(x[0])','cos(x[0])'),degree=0)

#Move mesh 1
ALE.move(mesh1,u)

#see what the moved mesh looks like
File('test/moved_mesh.pvd')<<mesh1


#Defining dof to vertex and vertex to dof maps for moved and reference regions
dof2v_W_1 = np.array(dof_to_vertex_map(W_1), dtype=int)
dof2v_W_2 = np.array(dof_to_vertex_map(W_2), dtype=int)    
v2dof_W_1 = np.array(vertex_to_dof_map(W_1), dtype=int)
v2dof_W_2 = np.array(vertex_to_dof_map(W_2), dtype=int)


#defining the problem variables
p1 = Function(W_1)
v1 = TestFunction(W_1)
p0 = Expression('x[1]*x[1]+x[0]*x[0]',degree=0)   #initial condition
p1 = interpolate(p0,W_1)

#the weak form and solution in the moved mesh
F = dot(grad(p1),grad(v1))*dx-p1*v1*dx
solve(F==0,p1)

#the solution in the moved mesh
File('test/P_in_moved_mesh.pvd')<< p1


#Define the function p2 on the reference domain and put its values in an auxillary vector Array 
p2 = Function(W_2)
Array = p2.vector().get_local()

#set the values of the auxillary vector equal to the values of p1 from the moved mesh using the dof maps
Array[v2dof_W_2[dof2v_W_1]] = p1.vector().get_local()

#Let p2 be the solution p1 but mapped backed into the reference domain          
p2.vector()[:]=Array

#see the result
File('test/P_in_ref_mesh.pvd')<< p2

And the error is:
Traceback (most recent call last):
File “test.py”, line 54, in
Traceback (most recent call last):
File “test.py”, line 54, in
Array[v2dof_W_2[dof2v_W_1]] = p1.vector().get_local()
Array[v2dof_W_2[dof2v_W_1]] = p1.vector().get_local()
ValueError: shape mismatch: value array of shape (223,) could not be broadcast to indexing result of shape (233,)
ValueError: shape mismatch: value array of shape (218,) could not be broadcast to indexing result of shape (229,)
--------------------------------------------------------------------------
mpirun noticed that process rank 0 with PID 8367 on node user1-HP-EliteBook-8470p exited on signal 6 (Aborted).

2

Why not simply save the initial state of the mesh, and revert your mesh changes after solving on the perturbed mesh? I.e.

mesh =UnitSquareMesh(10,10)
V=FunctionSpace(mesh, “CG”,1)
movement = project(Expression(("x[0]","x[1]"), degree=1),V)
ALE.move(mesh, movement)
# solve on perturbed mesh
#...
movement.vector()[:]*=-1
ALE.move(mesh, movement)
1 Like
3

That’s actually a very good idea and that might solve my problem. But I am genuinely curious about the compatibility of the dofmaps with the mpi. There is not a lot about it I can read. Maybe guide me to a book or document I can read would be extremely helpful. Thank you again.

4

As MPI can only transfer integers or vectors of i integers, you need to transform the dofmap into a transferrable format. See for instance https://mpi4py.readthedocs.io/en/stable/

5

Thank you dokken. This was actually the reason I asked the other question here
I think I am going to use an np array to save them and then use a dolfin vector to gather them. I have been asking many questions lately because I am learning MPI and simultaneously doing a project about it. But you have been such a great help so far to clarify things.

6

Do we have an update or a fix for this mapping ? @dokken @nami.

Unfortunately, in my case - I cannot revert my mesh changes back to the original mesh using ALE.

I would appreciate any alternative approach/help.

7

As legacy DOLFIN is no longer actively developed, I’m not adding patches for things such as this.

Why can’t you do this? Do you have a minimal reproducible example that highlights why you cannot follow this approach?

8

Hi @dokken,

I understand that legacy fenics is not being actively developed anymore. I am using Yu’s algorithm for my FSI solver (vanDANA) built in legacy FEniCS (2019.2.0.dev). Here if you look at Scheme 1 (Eq’s 46 - 49) or Scheme 2, this approach requires you to solve the solid momentum equation (Eq. 47) on the reference configuration and lagrange mulitplier (Eq. 48) on the current configuration of the mesh. Hence the need arises to map variables back and forth and also retain the mesh on either configuration at each time step. I hope that makes sense.

So here is a similar and more generalized MWE, as compared to @nami 's

To make the mesh :

from dolfin import *

m=40

mesh1 = UnitCubeMesh(m,m,m)

File('test/ref_mesh.pvd')<<mesh1

hdf = HDF5File(MPI.comm_world, "test/file_s.h5", "w")
hdf.write(mesh1, "/mesh")

To test the mapping : Note that this code also runs in parallel (for eg: -n 32)

from dolfin import *

mpi_comm = MPI.comm_world

# Read meshes
mesh1 = Mesh(mpi_comm)                                  #we move this one
hdf1 = HDF5File(mpi_comm, "test/file_s.h5", "r")
hdf1.read(mesh1, "/mesh", False)

mesh2 = Mesh(mpi_comm)                                  #we keep this one as the reference mesh
hdf2 = HDF5File(mpi_comm, "test/file_s.h5", "r")
hdf2.read(mesh2, "/mesh", False)

# --------------------------------------------
# Function to assign-vector in MPI
def vector_assign_in_parallel(x, y):

    x.vector().set_local(y.vector().get_local()[:])
    x.vector().apply("insert")

# --------------------------------------------
V_1 = VectorFunctionSpace(mesh1,'P',2)
p1 = Function(V_1); v1 = TestFunction(V_1)

print("DOF for the mesh variable: ", MPI.sum(mpi_comm, p1.vector().get_local().size), "\n", flush=True)

#defining a vector to move mesh1 with
u = Expression(('sin(x[0])','cos(x[0])','sin(x[0])+cos(x[0])'),degree=0)

#Move mesh 1
ALE.move(mesh1,u)

#see what the moved mesh looks like
File('test/moved_mesh.pvd')<<mesh1

#defining the problem variables
p0 = Expression((('x[2]*x[2]', 'x[0]*x[0]', 'x[1]*x[1]')), degree=0)   #initial condition
p1 = interpolate(p0, V_1)

#the weak form and solution in the moved mesh
F = inner(grad(p1),grad(v1))*dx - dot(p1,v1)*dx
solve(F==0, p1, solver_parameters={"newton_solver":{"linear_solver": "bicgstab", "preconditioner":'hypre_amg'}})

#the solution in the moved mesh
File('test/P_in_moved_mesh.pvd')<< p1

# ---------------------------------------------
#Define the function p2 on the reference domain and put its values in an auxillary vector Array
V_2 = VectorFunctionSpace(mesh2,'P',2)
p2 = Function(V_2)

vector_assign_in_parallel(p2, p1)

#see the result
File('test/P_in_ref_mesh.pvd')<< p2

# ---------------------------------------------
#Erase p1 vector and map p2 back to moved mesh
p1.vector().zero()                                      # Feel free to comment these 3 lines while testing the rest of the code
vector_assign_in_parallel(p1, p2)
File('test/P_remapped_in_moved_mesh.pvd')<< p1

Here is the error message:

Traceback (most recent call last):
  File "new.py", line 52, in <module>
    vector_assign_in_parallel(p2, p1)
  File "new.py", line 18, in vector_assign_in_parallel
    x.vector().set_local(y.vector().get_local()[:])
RuntimeError: 

*** -------------------------------------------------------------------------
*** DOLFIN encountered an error. If you are not able to resolve this issue
*** using the information listed below, you can ask for help at
***
***     fenics-support@googlegroups.com
***
*** Remember to include the error message listed below and, if possible,
*** include a *minimal* running example to reproduce the error.
***
*** -------------------------------------------------------------------------
*** Error:   Unable to successfully call PETSc function 'VecSetValuesLocal'.
*** Reason:  PETSc error code is: 63 (Argument out of range).
*** Where:   This error was encountered inside /tmp/src/dolfin/dolfin/la/PETScVector.cpp.
*** Process: 17
*** 
*** DOLFIN version: 2019.2.0.dev0
*** Git changeset:  f0a90b8fffc958ed612484cd3465e59994b695f9
*** -------------------------------------------------------------------------

However this MWE/code works fine (in parallel) for a smaller mesh (for eg: m=20). The interesting thing to notice here is that since I am using the same mpi_comm to read both the meshes, I expect the partitioning to remain the same, however this approach doesnt always seem to work (especially for larger mesh sizes).

Hence, here I am looking for an alternative approach to define vector_assign_in_parallel(x, y), so that it works at higher DOF’s.

9

Why cant you just map the mesh back and forth between the solves, instead of having two mesh objects, functionspaces etc. it is alot cheaper to keep to vectors in memory than the whole mesh twice.