Problems with mpi


251
views
-1
4 months ago by
Cassia  

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
4 months ago by
The essence of parallel FEM computation is mesh partitioning (that happens in you line ... = 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
4 months ago by
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. 
xc = centroid[0]; yc = centroid[1]
  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 1447, in _check_indices
    self.local_size()))
IndexError: expected indices to be in [0..0]
​
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):
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?
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.

Similar posts:
Search »