### Problems with mpi

251

views

-1

Dear Community,

I have two questions: The first one is if there is anyway to parallel only the solver instead of use * mpirun -np 2 python Twist_CLSM.py*. Second is: What is the reason of the following reason when I apply the aforementioned command?

```
Traceback (most recent call last):
File "Twist_CLSM.py", line 134, in <module>
Traceback (most recent call last):
File "Twist_CLSM.py", line 134, in <module>
xc = centroid[0.0]; yc = centroid[1.0]
File "/usr/lib/python2.7/dist-packages/dolfin/cpp/la.py", line 1465, in __getitem__
xc = centroid[0.0]; yc = centroid[1.0]
File "/usr/lib/python2.7/dist-packages/dolfin/cpp/la.py", line 1465, in __getitem__
indices = self._check_indices(indices)
File "/usr/lib/python2.7/dist-packages/dolfin/cpp/la.py", line 1422, in _check_indices
indices = self._check_indices(indices)
File "/usr/lib/python2.7/dist-packages/dolfin/cpp/la.py", line 1422, in _check_indices
raise TypeError("expected an int or a list or numpy array of "\
TypeError: expected an int or a list or numpy array of integers or a boolean numpy array as indices.
raise TypeError("expected an int or a list or numpy array of "\
TypeError: expected an int or a list or numpy array of integers or a boolean numpy array as indices.
--------------------------------------------------------------------------
mpirun noticed that process rank 0 with PID 0 on node ubuntu exited on signal 6 (Aborted).
```

Thanks so much in advance.

```
from fenics import*
import time
import ufl
import numpy as np
import os
from math import floor, ceil
from sympy import symbols
dt = 0.001 # Time Step
k = Constant(dt) # Time Step
t_end = 1.0 # Total simulation time
theta = Constant(0.5) # Interpolation schema
Dx = 150
Dy = Dx
mesh = RectangleMesh(Point(0.0,0.0), Point(1.0,1.0), Dx, Dy,'crossed')
V = VectorFunctionSpace(mesh,"Lagrange",2,dim=2)
L = FunctionSpace(mesh,"Lagrange",2)
n = FacetNormal(mesh)
h = CellSize(mesh)
beta=0.5
center = Point(0.5,0.75)
epsilon = (beta*((1.0/Dx)**(0.9)))
radius = 0.15
phi = Expression('1.0/(1.0+exp((sqrt((x[0]-A)*(x[0]-A) + (x[1]-B)*(x[1]-B))-r)/(eps)))',degree=2, eps=epsilon, A=center[0], B=center[1],r=radius)
l = Function(L); l0 = Function(L)
l.assign(interpolate(phi,L)); l0.assign(interpolate(phi,L))
v = Function(V); v0 = Function(V); ff=Function(V)
ff=(Expression(("2.0*sin(2.0*pi*x[1])*pow(sin(pi*x[0]),2)*cos(pi*t)", "-2.0*sin(2.0*pi*x[0])*pow(sin(pi*x[1]),2)*cos(pi*t)"),degree=2,pi=np.pi,t=0.0,element=V.ufl_element()))
v.assign(ff)
plot(l)
def level_set(v,v0,l0,l):
l_ = TestFunction(L)
nb = sqrt(dot(v,v))
tau = h*pow(2.0*nb,-1.0)
r = ( ((l-l0)/k) + theta*(inner(v,grad(l)))+ (1.0-theta)*(inner(v0,grad(l0))) )*tau*inner(v,grad(l_))*dx
def LS(l,v,l_):
return(inner(v,grad(l))*l_*dx)
l_ = TestFunction(L)
F = inner((l-l0)/k,l_)*dx + theta*LS(l,v,l_) + (1.0-theta)*LS(l0,v0,l_) + r
bc = []
begin("level_set")
J = derivative(F,l)
problem=NonlinearVariationalProblem(F,l,bc,J)
solver=NonlinearVariationalSolver(problem)
solver.solve()
end()
return(l)
def reinit(l,epsilon,beta,Dx,mesh):
# time-step
dtau = Constant(1.0/(beta*((1.0/Dx)**(1.1))))
# space definition
V = VectorFunctionSpace(mesh,"Lagrange",1,dim=2)
FE = FunctionSpace(mesh,"Lagrange",2)
# functions setup
phi = Function(FE); phi0 = Function(FE)
w = TestFunction(FE)
# intial value
phi0.assign(interpolate(l,FE))
# Unit normal vector (does not change during this process)
grad_n = project(grad(l),V)
n = grad_n/(sqrt(dot(grad_n,grad_n)))
# FEM linearization
F = dtau*(phi-phi0)*w*dx - 0.5*(phi+phi0)*dot(grad(w),n)*dx + phi*phi0*dot(grad(w),n)*dx + epsilon*0.5*dot(grad(phi+phi0),n)*dot(grad(w),n)*dx
# Newton-Raphson parameters
bc = []
E = 1e10; E_old = 1e10
cont = 0; num_steps = 10
for n in range(num_steps):
begin("Reinitialization")
solve(F == 0, phi, bc)
end()
error = (((phi - phi0)*dtau)**2)*dx
E = sqrt(abs(assemble(error)))
fail = 0
if (E_old < E ):
fail = 1
print('fail',"at:", cont)
break
phi0.assign(phi)
cont +=1
E_old = E
return phi
vfile = File("Twist10/velocity.pvd")
lfile = File("Twist10/level.pvd")
R = VectorFunctionSpace(mesh,"R",0,dim=2)
position = Function(V)
ve= Function(V)
position.assign(Expression(["x[0]","x[1]"], element=V.ufl_element()))
c = TestFunction(R)
out_dt = 0.1; count = 0
t = dt
if __name__ == "__main__":
while t < t_end:
l1 = level_set(v,v0,l0,l)
l2 = reinit(l1,epsilon,beta,Dx,mesh)
l.assign(interpolate(l2,L))
# Mass
Vo=assemble(conditional(gt(l,0.5),1.0,0.0)*dx)
# Center of mass
volume = assemble(conditional(ge(l,0.5),Constant(1.0),0.0)*dx)
centroid = assemble(conditional(ge(l,0.5),dot(c,position),0.0)*dx)
centroid /= volume
xc = centroid[0.0]; yc = centroid[1.0]
# Velocity
u_c = assemble(conditional(ge(l,0.5),dot(c,v),0.0)*dx)
u_c /= volume
vxc = u_c[0][0]; vyc = u_c[1][0]
#myfile3.write('%e %e' '\n' % (vyc, t))
print("%e %e %e %e" %(Vo, yc, vyc, t))
# Save solution
if (t >= float(count)*out_dt):
count+=1
v.rename("velocity", "velocity")
l.rename("level", "level")
vfile << v
lfile << l
t += dt
l0.assign(l)
v0.assign(v)
ff.t=t
v=interpolate(ff,V)
v.rename("velocity", "velocity")
l.rename("level", "level_set")
vfile << v
lfile << l
```

Community: FEniCS Project

### 2 Answers

3

The essence of parallel FEM computation is mesh partitioning (that happens in you line

You can do "parallelise only the solver" if you construct the mesh on

`... = RectangleMesh()`

). Parallelise "only the solver" doesn't make much sense. FEniCS assembles all the objects into your linear algebra's backend (e.g. PETSc) and parallel computation takes place there.You can do "parallelise only the solver" if you construct the mesh on

`MPI_COMM_SELF`

(so every process holds whole mesh), assemble vectors and matrices using FEniCS methods and then do your custom parallel solve magic and communicate results.
1

Concerning your second question, as the error suggests an integer is expected

`xc = centroid[0]; yc = centroid[1]`

I have tried the syntax you posted before and still have the same problem. Thanks for your time.

written
4 months ago by
Cassia

2

Well, the error is different now (

`IndexError`

vs. `TypeError`

). When you run your code in parallel, `centroid`

is seemingly distributed only to the last process (only there, it has the right size). You can check it with `print(MPI.rank(mpi_comm_world()), centroid.local_size())`

. Not sure what to do about it though...
written
4 months ago by
Adam Janecka

1

In parallel, it is also usually required to set (for example):

You do not have this error message, but it may be related. An issue like this would work in serial since the function dofs are not distributed. You could also make each mpi process execute the offending statement serially (using a barrier) and explicitly check the size of the centroid result. If I understand the error message correctly, centroid is returned of size < 2.

See below for set_allow_extrapolation issues:

https://www.allanswered.com/post/lobkw/#qvowl

and

https://www.allanswered.com/local/search/page/?q=set_allow_extrapolation

```
u1 = interpolate(expr,V1)
# allow extrapolation on that function
u1.set_allow_extrapolation(True)
```

You do not have this error message, but it may be related. An issue like this would work in serial since the function dofs are not distributed. You could also make each mpi process execute the offending statement serially (using a barrier) and explicitly check the size of the centroid result. If I understand the error message correctly, centroid is returned of size < 2.

See below for set_allow_extrapolation issues:

https://www.allanswered.com/post/lobkw/#qvowl

and

https://www.allanswered.com/local/search/page/?q=set_allow_extrapolation

written
4 months ago by
jwinkle

Is it possible your centroid conditional is not passing?

Try printing the value of the centroid (or checking its size is non-zero) in each time step for debugging.

`centroid = assemble(conditional(ge(l,0.5),dot(c,position),0.0)*dx)`

Try printing the value of the centroid (or checking its size is non-zero) in each time step for debugging.

written
4 months ago by
jwinkle

It is passing and it is evaluated just fine for one processor.

written
4 months ago by
Cassia

Please login to add an answer/comment or follow this question.