Hi Floriane, On 05/09/17 12:43, Floriane Gidel [RPG] wrote:
Hi David,
This is a big piece of code where I solve nonlinear equations in a 2D domain with mixed systems involving vector functions.
I'd like to run the whole code in parallel except from some lines such as saving txt files, to avoid overwriting it or saving in several files. Is there a command to select a specific process? Andrew once mentioned using op2.MPI.comm.rank. Is that a solution? How can I use it?
The terminology that MPI uses is that you have a *communicator* comprised of some number of processes, each of these is indexed by its *rank*. All Firedrake objects should have the communicator they are defined on accessible as the .comm attribute (if you find one that is not, please report a bug). When you build a mesh, and don't provide a communicator, it is built on the default COMM_WORLD communicator. Objects built on top of the mesh inherit the mesh's communicator. So if you have a mesh, and only want to print something on rank 0. You can write: if mesh.comm.rank == 0: print("I am on rank 0") An important thing to note: Most operations that you do in firedrake are *collective* over the communicator of the participating objects. That means that all ranks need to participate, otherwise you will get errors or your code might hang. So for example you can't do this: mesh = Mesh(...) V = FunctionSpace(mesh, ...) v = TestFunction(V) f = assemble(v*dx) if f.comm.rank == 0: # This line collectively accesses the data in f # it therefore will hang. print(f.dat.data) You can do: # Everyone accesses data, rank 0 prints something fdata = f.dat.data if f.comm.rank == 0: print(fdata)
The error when appending the array is not such an issue as I need it only for visualisation purpose.
For instance, I want to save a 1D function h_1D
mesh_1D = IntervalMesh(Nx,Lx)
V_1D = FunctionSpace(mesh_1D, "CG", 1)
h_1D = Function(V_1D)
on a 2D surface. So I define a 2D mesh as follows
mesh_2D = RectangleMesh(Nx,1,Lx,1.0,quadrilateral=True)
V_2D = FunctionSpace(mesh_2D, "CG", 1)
h_2D = Function(V_2D)
Then, I save in Indx[i] the indices for which x_2D[i] = x_1D[i]:
Indx = []
fori inrange(Nx+1):
Indx.append([item foritem inrange(len(mesh_2D.coordinates.dat.data[:,0])) ifmesh_2D.coordinates.dat.data[item,0] == mesh_1D.coordinates.dat.data[i]])
The idea is that I can then project h_1D to the 2D mesh through h_2D:
fori inrange(Nx+1):
h_2D.dat.data[Indx[i]] = h_1D.dat.data[i]
This works fine when I run the code on one core, but when running in parallel I get an index error:
File "main.py", line 932, in <listcomp>
Indx.append([item for item in range(len(mesh_2D.coordinates.dat.data[:,0])) if mesh_2D.coordinates.dat.data[item,0] == mesh_1D.coordinates.dat.data[j]])
IndexError: index 125 is out of bounds for axis 0 with size 125
Both the 2D mesh and the 1D mesh are distributed amongst the MPI processes, so you have no guarantee that there are Nx+1 entries in the 1D coordinate array. You probably want: for i in range(len(mesh_1D.coordinates.data.data)): Indx.append(...) The same is true of writing: Once you have made the Indx array, I would just write: h_2D.dat.data[Index, 0] = h_1D.dat.data_ro[:]
I can comment it, but then I get other (Firedrake-related) errors that I don't get otherwise, for example:
File "main.py", line 970, in <module>
mesh_WM= CubeMesh(1,1,1,1.0)
...
raise RuntimeError("Mesh must have at least one cell on every process")
RuntimeError: Mesh must have at least one cell on every process
I don't understand while running in parallel causes errors for the mesh definition?
You made a 1x1x1 cube mesh which has 6 cells. If you tried to run on more than 6 processes (mpiexec -n 7, for example), then you would get this error, because every process must own at least one cell. Possibly even with fewer than 6 processes you might get this error depending on the exact mesh decomposition. Cheers, Lawrence