MPI and tripcolor

[qlyan913]

Hi,

I am trying to solve a PDE with a 2D piecewise constant coefficient. One step is define the coefficient function and plot it.
Here, I want to define the function on [0,4]\times[0,4] :

When I use mpiexec -n 2 to run the same code, it gives:

which is clearly incorrect.

Here is the sample code:

import matplotlib.pyplot as plt
from firedrake import *
from firedrake.petsc import PETSc
from firedrake.__future__ import interpolate
import numpy as np
L=4
nx=L
ny=L
a0=1   # pc constant range from [a0, a1]
a1=10 
coefplotfile='coef.png'
mesh = SquareMesh(nx,ny,L,quadrilateral=True)
aelt = 'DG'
adeg = 0
V=FunctionSpace(mesh,aelt,adeg)
aexpr = Function(V)
av=aexpr.vector()
nc2=nx*ny
np.random.seed(5)
aval=a0+np.random.rand(nc2)*(a1-a0)
print("aval: ", aval)
n_min,n_max=av.local_range()
print("common world rank", COMM_WORLD.rank, "set", n_min, "to ", n_max)
av.set_local(aval[n_min:n_max])
print("common world rank", COMM_WORLD.rank,aval[n_min:n_max])
A = assemble(interpolate(aexpr, V))
# evaluate coefficient, save to file and plot
print("integral of A", assemble(A*dx))
plt.clf()
fig, axes = plt.subplots()
collection = tripcolor(A, axes=axes)
fig.colorbar(collection);
plt.title('coefficient')
plt.savefig(coefplotfile, dpi=500)
PETSc.Sys.Print("> coefficient plotted to {}".format(coefplotfile))

Any helps are greatly appreciated.

If I run with N=1 process saved to h5 file and load it with M=2 then plot the loaded function, the results is:

If I save function with M=2 process and load it with N=1 then plot the functions, then:

But if I run with 1 process and save the figure it should be:

So, there is a kind of permutation.

[ksagiyam]

Associating DoFs of a Function distributed over N processes (say N=1) to those of the same Function distributed over M processes (say M=2) is not trivial.
You could create a randomised Function with N=1, save it with CheckpointFile (https://www.firedrakeproject.org/checkpointing.html), later load it with M=2, and then plot the loaded Function to do this comparison.

[qlyan913]

Thanks. When I save function with N=1 and load it with M=2 then plot the loaded function, the results is still incorrect. On the other hand, if I save function with M=2 and load it with N=1 then plot the functions, then I can get what expected.

[ksagiyam]

Could you post what you get for those cases?

[qlyan913]

I just updated the figure in the post.

[ksagiyam]

The followings worked fine for me:

Run with mpiexec -n 1 python ...:

from firedrake import *
import matplotlib.pyplot as plt


mesh = UnitSquareMesh(4, 4, quadrilateral=True, name="m")
x, y = SpatialCoordinate(mesh)
V = FunctionSpace(mesh, "DG", 0)
f = Function(V, name="f").interpolate(x + 2 * y)
with CheckpointFile("random.h5", "w") as afile:
    afile.save_function(f)
#
fig = plt.gcf()
ax = plt.gca()
ax.set_xlim([0, 1])
ax.set_ylim([0, 1])
t = tripcolor(f, axes=ax, vmin=0.0, vmax=2.5)
fig.colorbar(t)
plt.savefig("saved_with_1.pdf")

saved_with_1.pdf

Run with mpiexec -n 2 python ...:

from firedrake import *
import matplotlib.pyplot as plt


with CheckpointFile("random.h5", "r") as afile:
    mesh = afile.load_mesh("m")
    f = afile.load_function(mesh, "f")
#
fig = plt.gcf()
ax = plt.gca()
ax.set_xlim([0, 1])
ax.set_ylim([0, 1])
t = tripcolor(f, axes=ax, vmin=0.0, vmax=2.5)
fig.colorbar(t)
plt.savefig(f"loaded_with_2_{mesh.comm.rank}.pdf")

loaded_with_2_0.pdf
loaded_with_2_1.pdf

[qlyan913]

Thanks. But this is not what I want. I want to generate function with MPI -n 2, and get a correct plot. I also updated the post. It seems that there is a permutation of function after using the checkpoint.

[dham]

The basic problem is that you are relying on the DoF numbering in order to set your Function values. However, the DoF numbering isn't a part of the mathematical properties of the finite element space, it's an implementation detail. It generally won't be the same if you run on a different number of processors. We also don't promise that it won't change in future Firedrake or PETSc versions. This means that what you are observing is expected behaviour, it's not wrong.

If you want to recover the same forcing function regardless of parallel decomposition, then you need to ensure that your function values depend on the coordinate location of the cell, not the numbering of the DoFs. The procedure for interpolating external data provides a mechanism for doing this that does not depend on the DoF ordering: https://www.firedrakeproject.org/interpolation.html#interpolation-from-external-data

[qlyan913]

Thanks for your suggestion. But when I need to solve the eigenvalue problem with slepcy4py, I still need

eigenfun.vector().set_local(v)

to assign the value of eigenvector to eigenfunction in fem space. So, it still depends on the DoF in the end.

Do you know how to avoid this issue in such a situation if I want use mpi to solve eigenvalue problem and do the correct visualization? Thank you.