Couple Escript with ESyS-Particle and Suggestions to install software properly

I have forwarded you a copy of the paper.

You should be able to install escript and esys-particle on top of each other into the same diretory "esys/". The order might be critical but if there is a problem this would be in esys/__init__.py. We haven't done this for a while and I assume you need to switch off MPI for escript. Please don;t ask me for a script how to link the codes.

Thanks Lutz Gross, that solved my question.

Hi Lutz Gross,

Currently I'm running coupled Escript-ESyS-Particle simulation similar with https://yade-dem.org/doc/FEMxDEM.html.
Multiprocessing is used in the Escript to arrange the independent ESyS-Particle simulation by call('mpirun -np 2 `which esysparticle` load.py', shell=True).
Now I want to run the coupled simulation on more than one nodes.
I have asked this question https://groups.google.com/forum/#!topic/mpi4py/YBr70v17bsg.
Howevery, if I used the method method mentioned above using mpi4py, then I have to turn MPI for Escript, which will contradict
with the MPI of ESyS-Particle.
Could you please give me some suggestion?

Best regards,
Jiadun

In principle you can launch escript using mpirun and mpi4py. The run-escript shell is just for convenience to deal with different mph flavors, threading, etc. However it is vital that mpi4py, esys-particle and encrypt are linked against the same mpi shared libraries.

The otehr thing the launcher does is to add the underlying C++ libraries and escript python libs to LD_LIBRARY_PATH and PYTHONPATH.
(Usually esys/ and lib/ under the install location).

So make sure you've made the appropriate adjustments.

Thank you very much for your help.
I'll try it.

This question was expired because it remained in the 'Open' state without activity for the last 15 days.

Hi Lutz Gross and Joel Fenwick,

I think that mpi is properly turned on for escript. However I have problem where I need your help.

Details to show that MPI option is correctly installed. The problem is at the end of this thread.

(1) [bwjliu@int2 Escript_MPI]$ run-escript
Python 2.7.9 (default, Mar 9 2015, 11:04:54)
[GCC 4.4.7 20120313 (Red Hat 4.4.7-11)] on linux2
Type "help", "copyright", "credits" or "license" for more information.
>>> from esys.lsm import *
>>> from esys.escript import *
>>> from mpi4py import MPI
>>>
(2)[bwjliu@int2 Escript_MPI]$ run-escript -p 24 example01a.py
Expected Number of time outputs is: 200.0
time step 1 at t=9.131063e+01 days completed. total energy = 1.554331e+14.
time step 2 at t=1.826213e+02 days completed. total energy = 1.554331e+14.
time step 3 at t=2.739319e+02 days completed. total energy = 1.554331e+14.
time step 4 at t=3.652425e+02 days completed. total energy = 1.554331e+14.
time step 5 at t=4.565531e+02 days completed. total energy = 1.554331e+14.
time step 6 at t=5.478637e+02 days completed. total energy = 1.554331e+14.
time step 7 at t=6.39174
(3) refer to https://groups.google.com/forum/#!topic/mpi4py/YBr70v17bsg
Here totalFormulationDynamicQuasiStaticFDEM_Biaxial_2.py is the Escript python script wherein multiprocessing and call(,shell = True) is used to invoke the separate ESyS-Particle simulation.
[bwjliu@int2 FDEM_2_T_MPI]$ run-escript totalFormulationDynamicQuasiStaticFDEM_Biaxial_2.py
8
-1.99206367725
CSubLatticeControler::initMPI()
CSubLatticeControler::initMPI()
CSubLatticeControler::initMPI()
CSubLatticeControler::initMPI()
CSubLatticeControler::initMPI()
CSubLatticeControler::initMPI()
CSubLatticeControler::initMPI()
CSubLatticeControler::initMPI()
slave started at local/global rank 0 / 1
slave started at local/global rank 0 / 1
slave started at local/global rank 0 / 1
slave started at local/global rank 0 / 1
set time step to 0
set time step to 0
set time step to 0
set time step to 0
reading reading 502 RotFriction interactions
reading 503 RotFriction interactions
reading 504 RotFriction interactions
slave started at local/global rank 0 / 1
set time step to 0
reading 504 RotFriction interactions
504 RotFriction interactions
slave started at local/global rank 0 / 1
constructing FieldMaster for field e_kin
constructing FieldMaster for field e_kin
set time step to 0
constructing FieldMaster for field e_kin
reading 499 RotFriction interactions
4e+02 days completed. total energy = 1.554331e+14.

The detail of the problem is

[bwjliu@int2 FDEM_2_T_MPI]$ run-escript -p 16 totalFormulationDynamicQuasiStaticFDEM_Biaxial_2.py
4
Traceback (most recent call last):
  File "totalFormulationDynamicQuasiStaticFDEM_Biaxial_2.py", line 153, in <module>
Traceback (most recent call last):
  File "totalFormulationDynamicQuasiStaticFDEM_Biaxial_2.py", line 153, in <module>
Traceback (most recent call last):
  File "totalFormulationDynamicQuasiStaticFDEM_Biaxial_2.py", line 153, in <module>
Traceback (most recent call last):
  File "totalFormulationDynamicQuasiStaticFDEM_Biaxial_2.py", line 153, in <module>
Traceback (most recent call last):
Traceback (most recent call last):
  File "totalFormulationDynamicQuasiStaticFDEM_Biaxial_2.py", line 153, in <module>
Traceback (most recent call last):
  File "totalFormulationDynamicQuasiStaticFDEM_Biaxial_2.py", line 153, in <module>
Traceback (most recent call last):
  File "totalFormulationDynamicQuasiStaticFDEM_Biaxial_2.py", line 153, in <module>
Traceback (most recent call last):
  File "totalFormulationDynamicQuasiStaticFDEM_Biaxial_2.py", line 153, in <module>
Traceback (most recent call last):
  File "totalFormulationDynamicQuasiStaticFDEM_Biaxial_2.py", line 153, in <module>
Traceback (most recent call last):
  File "totalFormulationDynamicQuasiStaticFDEM_Biaxial_2.py", line 153, in <module>
Traceback (most recent call last):
  File "totalFormulationDynamicQuasiStaticFDEM_Biaxial_2.py", line 153, in <module>
Traceback (most recent call last):
  File "totalFormulationDynamicQuasiStaticFDEM_Biaxial_2.py", line 153, in <module>
    lengthHeightInitial += (1.0/(topEndNode - topBeginNode + 1))*d.getTupleForDataPoint(i)[1]
    lengthHeightInitial += (1.0/(topEndNode - topBeginNode + 1))*d.getTupleForDataPoint(i)[1]
RuntimeError: Error - Data::getValueOfDataPointAsTuple: invalid sampleNo.
    lengthHeightInitial += (1.0/(topEndNode - topBeginNode + 1))*d.getTupleForDataPoint(i)[1]
RuntimeError: Error - Data::getValueOfDataPointAsTuple: invalid sampleNo.
    lengthHeightInitial += (1.0/(topEndNode - topBeginNode + 1))*d.getTupleForDataPoint(i)[1]
RuntimeError: Error - Data::getValueOfDataPointAsTuple: invalid sampleNo.
    lengthHeightInitial += (1.0/(topEndNode - topBeginNode + 1))*d.getTupleForDataPoint(i)[1]
RuntimeError: Error - Data::getValueOfDataPointAsTuple: invalid sampleNo.
    lengthHeightInitial += (1.0/(topEndNode - topBeginNode + 1))*d.getTupleForDataPoint(i)[1]
RuntimeError: Error - Data::getValueOfDataPointAsTuple: invalid sampleNo.
RuntimeError: Error - Data::getValueOfDataPointAsTuple: invalid sampleNo.
    lengthHeightInitial += (1.0/(topEndNode - topBeginNode + 1))*d.getTupleForDataPoint(i)[1]
RuntimeError: Error - Data::getValueOfDataPointAsTuple: invalid sampleNo.
Traceback (most recent call last):
  File "totalFormulationDynamicQuasiStaticFDEM_Biaxial_2.py", line 153, in <module>
  File "totalFormulationDynamicQuasiStaticFDEM_Biaxial_2.py", line 153, in <module>
Traceback (most recent call last):
  File "totalFormulationDynamicQuasiStaticFDEM_Biaxial_2.py", line 153, in <module>
    lengthHeightInitial += (1.0/(topEndNode - topBeginNode + 1))*d.getTupleForDataPoint(i)[1]
    lengthHeightInitial += (1.0/(topEndNode - topBeginNode + 1))*d.getTupleForDataPoint(i)[1]
Traceback (most recent call last):
  File "totalFormulationDynamicQuasiStaticFDEM_Biaxial_2.py", line 153, in <module>
RuntimeError: Error - Data::getValueOfDataPointAsTuple: invalid sampleNo.
RuntimeError: Error - Data::getValueOfDataPointAsTuple: invalid sampleNo.
    lengthHeightInitial += (1.0/(topEndNode - topBeginNode + 1))*d.getTupleForDataPoint(i)[1]
RuntimeError: Error - Data::getValueOfDataPointAsTuple: invalid sampleNo.
    lengthHeightInitial += (1.0/(topEndNode - topBeginNode + 1))*d.getTupleForDataPoint(i)[1]
    lengthHeightInitial += (1.0/(topEndNode - topBeginNode + 1))*d.getTupleForDataPoint(i)[1]
RuntimeError: Error - Data::getValueOfDataPointAsTuple: invalid sampleNo.
RuntimeError: Error - Data::getValueOfDataPointAsTuple: invalid sampleNo.
    lengthHeightInitial += (1.0/(topEndNode - topBeginNode + 1))*d.getTupleForDataPoint(i)[1]
RuntimeError: Error - Data::getValueOfDataPointAsTuple: invalid sampleNo.
    lengthHeightInitial += (1.0/(topEndNode - topBeginNode + 1))*d.getTupleForDataPoint(i)[1]
RuntimeError: Error - Data::getValueOfDataPointAsTuple: invalid sampleNo.
    lengthHeightInitial += (1.0/(topEndNode - topBeginNode + 1))*d.getTupleForDataPoint(i)[1]
RuntimeError: Error - Data::getValueOfDataPointAsTuple: invalid sampleNo.
    lengthHeightInitial += (1.0/(topEndNode - topBeginNode + 1))*d.getTupleForDataPoint(i)[1]
RuntimeError: Error - Data::getValueOfDataPointAsTuple: invalid sampleNo.
application called MPI_Abort(MPI_COMM_WORLD, 1) - process 6
application called MPI_Abort(MPI_COMM_WORLD, 1) - process 15
application called MPI_Abort(MPI_COMM_WORLD, 1) - process 14
application called MPI_Abort(MPI_COMM_WORLD, 1) - process 11
application called MPI_Abort(MPI_COMM_WORLD, 1) - process 8
application called MPI_Abort(MPI_COMM_WORLD, 1) - process 9
application called MPI_Abort(MPI_COMM_WORLD, 1) - process 10
application called MPI_Abort(MPI_COMM_WORLD, 1) - process 13

Best regards,
Jiadun

Without knowing your code in details, I would guess that problem is in assumed numbering of data points. Data objects are distributed over the available process, ie. when you object has N data points then under np ranks each rank holds roughly N/np data points where for some data objects a data overlap is introduced. So my guess is that the i loop around the d.getTupleForDataPoint(i) call does not consider this?

Thanks Lutz Gross, that solved my question.

Hi Lutz Gross,

It seems that there are some functions which are not mpi safe when running Escript by using run-escript -p NP some.py.

From example
(Function(dom).getX()).getNumberOfDataPoints()
d.getTupleForDataPoint(i)
d.setValueOfDataPoint()

How to avoid it?

Correct, these functions are not MPI-safe and hence users are not encouraged to use them (eg. they are not listed in the users guide). An appropriate replacement depends on the context?

Hi Lutz Gross,

I'm running a coupled FEM-DEM simulation similar to https://yade-dem.org/doc/FEMxDEM.html.
Coefficients need to be set for Gauss Points like https://answers.launchpad.net/escript-finley/+question/168042.

Best regards,
Jiadun

The question is what <everything you want to do with in_loc > means. Can you point me to the point where this happens in your code?

Hi Lutz Gross,
Here is beginning part of the code, the problem comes up at the for loop where displacement is set.
from esys.escript import *
from esys.escript.linearPDEs import LinearPDE,LinearPDESystem, SolverOptions
from esys.finley import ReadMesh,Rectangle
from numpy import zeros,ones
from esys.escript import util
from esys.pycad import * #domain constructor
from esys.pycad.gmsh import Design #Finite Element meshing package
import numpy,os,shutil,time
from esys.weipa import saveVTK
from esys.escript.pdetools import Projector
from subprocess import call #call shell command
import multiprocessing #run DEM in turn
from mpi4py import MPI
comm = MPI.COMM_WORLD
size = comm.Get_size()
rank = comm.Get_rank()
## record start time
start = time.time()

# make sure path exists
save_path= os.path.join("Results")
mkDir(save_path)

def runOdem(num):
 os.chdir("Odem_{0:0d}".format(num))
 call('mpirun -np 2 `which esysparticle` PLoad2.py',shell=True)
 os.chdir(base_dir)

def runIdem(num):
 os.chdir("Idem_{0:0d}".format(num))
 call('mpirun -np 2 `which esysparticle` PLoadRestart2.py',shell=True)
 os.chdir(base_dir)

rho = 0.00169254697901
lx = 400.0
ly = 800.0
dom = Rectangle(l0=lx,l1=ly,n0=1,n1=2,order=1,integrationOrder=2)
#dom = ReadMesh("lin.fly",1)
base_dir = os.getcwd()
#print (Function(dom).getX()).getNumberOfDataPoints()
GaussNum = 8
ProcessesNum = GaussNum

###copy converged stress fromlast step
for i in range(GaussNum):
 call('cp Cdem_{0:0d}/{1:0d}_Stress_Cauchy.dat ./'.format(i,i),shell=True)
 call('cp Cdem_{0:0d}/{1:0d}_Stress_Nominal.dat ./'.format(i,i),shell=True)

#dom = ReadMesh("1.fly",integrationOrder=1)
pde = LinearPDESystem(dom)
#pde = LinearPDE(dom,numEquations=dom.getDim(),numSolutions=dom.getDim())
pde.getSolverOptions().setSolverMethod(SolverOptions.HRZ_LUMPING)
pde.setSymmetryOn()
pde.setValue(D=rho*kronecker(dom))

local_damping_a = Vector(0,Solution(dom))
CauchyStress=Tensor(0.,Function(dom))
NominalStress=Tensor(0.,Function(dom))
stressOld = Tensor(0.,Function(dom))
R = Tensor(0.0,Function(dom))
U = Tensor(0.0,Function(dom))
F = Tensor(0.0,Function(dom))
F = Tensor(0.0,Function(dom))
E = Tensor(0.0,Function(dom))
macroRotation = Scalar(0.0,Function(dom))
averageRotation = Scalar(0.0,Function(dom))

#stress=Tensor(0.,ReducedFunction(dom))
xInitial = dom.getX()

##From biaxialSmooth YADE FEMXDEM
proj = Projector(dom)
sig = proj(CauchyStress) # project Gauss point value to nodal value
sig_bounda = interpolate(sig,FunctionOnBoundary(dom)) # interpolate
traction = matrix_mult(sig_bounda,dom.getNormal()) # boundary traction
Bx = FunctionOnBoundary(dom).getX()

bottomSurf = whereZero(Bx[1])
topSurf = whereZero(Bx[1]-sup(Bx[1]))
rightSurf = whereZero(Bx[0]-sup(Bx[0]))
leftSurf = whereZero(Bx[0])

tractTop = traction*topSurf # traction at top surface
forceTop = integrate(tractTop,where=FunctionOnBoundary(dom)) # resultant force at top
lengthTop = integrate(topSurf,where=FunctionOnBoundary(dom)) # length of top surface

v_last = Vector(0.,Solution(dom))
d_last = Vector(0.,Solution(dom))
d = Vector(0.,Solution(dom))

#####################
## DR parameter
DR= True
alpha = 2.0
XI = 0.005
safeFactor = 0.9

tol = 1.0e-2
rtol = 1.0e-2

totalSteps = 0
loadingSteps = 0

openfilename = "runTime.data"
openfile = open(openfilename,"r")
run_time = openfile.readlines()
openfile.close()

for i in range(len(run_time)):
 totalSteps += float(run_time[i].split()[1])

################################################################################################################

###########################################################################################################################
###read displacement from Isotropic compression to calculate delt_d

bottomBeginNode = 0
bottomEndNode = 1
topBeginNode = 4
topEndNode = 5

openfilename = "iso_disp.data"
openfile = open(openfilename,"r")
disp_iso = openfile.readlines()
openfile.close()

for i in xrange(6):
 d.setValueOfDataPoint(i,(float(disp_iso[i].split()[1]),float(disp_iso[i].split()[2])))

 run-escript -p 6 totalFormulationDynamicQuasiStaticFDEM_Biaxial_2.py
Traceback (most recent call last):
  File "totalFormulationDynamicQuasiStaticFDEM_Biaxial_2.py", line 151, in <module>
Traceback (most recent call last):
  File "totalFormulationDynamicQuasiStaticFDEM_Biaxial_2.py", line 151, in <module>
Traceback (most recent call last):
  File "totalFormulationDynamicQuasiStaticFDEM_Biaxial_2.py", line 151, in <module>
Traceback (most recent call last):
  File "totalFormulationDynamicQuasiStaticFDEM_Biaxial_2.py", line 151, in <module>
Traceback (most recent call last):
  File "totalFormulationDynamicQuasiStaticFDEM_Biaxial_2.py", line 151, in <module>
Traceback (most recent call last):
  File "totalFormulationDynamicQuasiStaticFDEM_Biaxial_2.py", line 151, in <module>
    d.setValueOfDataPoint(i,(float(disp_iso[i].split()[1]),float(disp_iso[i].split()[2])))
    d.setValueOfDataPoint(i,(float(disp_iso[i].split()[1]),float(disp_iso[i].split()[2])))
    d.setValueOfDataPoint(i,(float(disp_iso[i].split()[1]),float(disp_iso[i].split()[2])))
RuntimeError: RuntimeError: DataExpanded::copyDataPoint: invalid sampleNo.
RuntimeError: DataExpanded::copyDataPoint: invalid sampleNo.
DataExpanded::copyDataPoint: invalid sampleNo.
    d.setValueOfDataPoint(i,(float(disp_iso[i].split()[1]),float(disp_iso[i].split()[2])))
RuntimeError: d.setValueOfDataPoint(i,(float(disp_iso[i].split()[1]),float(disp_iso[i].split()[2])))
DataExpanded::copyDataPoint: invalid sampleNo.
RuntimeError: DataExpanded::copyDataPoint: invalid sampleNo.
    d.setValueOfDataPoint(i,(float(disp_iso[i].split()[1]),float(disp_iso[i].split()[2])))
RuntimeError: DataExpanded::copyDataPoint: invalid sampleNo.
application called MPI_Abort(MPI_COMM_WORLD, 1) - process 0
application called MPI_Abort(MPI_COMM_WORLD, 1) - process 1
application called MPI_Abort(MPI_COMM_WORLD, 1) - process 3
application called MPI_Abort(MPI_COMM_WORLD, 1) - process 4
application called MPI_Abort(MPI_COMM_WORLD, 1) - process 5
application called MPI_Abort(MPI_COMM_WORLD, 1) - process 2

Best regards,
Jiadun

The relevant code fragment is
 ###########################################################################################################################
###read displacement from Isotropic compression to calculate delt_d
openfilename = "iso_disp.data"
openfile = open(openfilename,"r")
disp_iso = openfile.readlines()
openfile.close()

for i in xrange(6):
 d.setValueOfDataPoint(i,(float(disp_iso[i].split()[1]),float(disp_iso[i].split()[2])))

Question: where is the 6 in xrange(6) coming from?

dom = Rectangle(l0=lx,l1=ly,n0=1,n1=2,order=1,integrationOrder=2)
So there are 6 nodes.

This is tricky business as it is not really possible to predict on which processor a node or element will end up (plus there can be overlapping). However there is a so called reference ID which is assigned to each node or element and which is unique across all ranks. This is number is set by the first column in the node as well as the element sections of the fly file (see users guide.)

for i in xrange(6):

    d.setValueOfDataPoint(i,(float(disp_iso[i].split()[1]),float(disp_iso[i].split()[2])))

getReferenceIDFromDataPointNo(i)

This is tricky business as it is not really possible to predict on which processor a node or element will end up (plus there can be overlapping). However there is a so called reference ID which is assigned to each node or element and which is unique across all ranks. This is number is set by the first column in the node as well as the element sections of the fly file (see users guide.). If you use Rectangle IDs are generated by counting along the x axis.

assuming that the first column in "iso_disp.data" defines the reference number then your code would like like:

v={}
for l in disp_iso:
     v[int(l.split()[0]) = (float(disp_iso[i].split()[1]),float(disp_iso[i].split()[2]))

for i in xrange(6):
    id=d.getFunctionSpace().getReferenceIDFromDataPointNo(i)
    d.setValueOfDataPoint(i, v[id])

A similar technique you can apply to elements.

I tried the method above. Here is the result.

[bwjliu@int1 setEscriptDisp]$ run-escript -p 2 setDisp.py
Traceback (most recent call last):
  File "setDisp.py", line 22, in <module>
Traceback (most recent call last):
  File "setDisp.py", line 22, in <module>
    d.setValueOfDataPoint(i, v[Id])
    d.setValueOfDataPoint(i, v[Id])
RuntimeError: RuntimeError: DataExpanded::copyDataPoint: invalid sampleNo.DataExpanded::copyDataPoint: invalid sampleNo.

application called MPI_Abort(MPI_COMM_WORLD, 1) - process 1

File setDisp.py
from esys.escript import *
from esys.finley import Rectangle

lx = 400.0
ly = 800.0
dom = Rectangle(l0=lx,l1=ly,n0=1,n1=2,order=1,integrationOrder=2)

## displacement
d = Vector(0.,Solution(dom))

openfilename = "iso_disp.data"
openfile = open(openfilename,"r")
disp_iso = openfile.readlines()
openfile.close()

##set disp
v={}
for i in range(len(disp_iso)):
 v[int(disp_iso[i].split()[0])] = (float(disp_iso[i].split()[1]),float(disp_iso[i].split()[2]))
for i in xrange(6):
 Id=d.getFunctionSpace().getReferenceIDFromDataPointNo(i)
 d.setValueOfDataPoint(i, v[Id])
 #print Id

##write disp to check
DISP = FileWriter('disp_check.data')
for i in range(d.getNumberOfDataPoints()):
 DISP.write(str(int(i))+' '
   +str(float(d.getTupleForDataPoint(i)[0]))+' '
   +str(float(d.getTupleForDataPoint(i)[1]))+'\n')
DISP.close()

exit()

File iso_disp.data
0 0.0 0.0
1 -1.80853237932 0.0
2 0.286888846683 -1.40356468889
3 -1.51239344602 -1.77001262807
4 0.80758326687 -2.92594370955
5 -0.982979988092 -3.42311449366

The line "for i in xrange(6):" should be read as "for i in xrange(d.getNumberOfDataPoints()):"

After using for i in xrange(d.getNumberOfDataPoints()):, only part of the nodes displacement are set.

[bwjliu@int2 Escript_MPI]$ more setDisp.py
from esys.escript import *
from esys.finley import Rectangle

lx = 400.0
ly = 800.0
dom = Rectangle(l0=lx,l1=ly,n0=1,n1=2,order=1,integrationOrder=2)

## displacement
d = Vector(0.,Solution(dom))

openfilename = "iso_disp.data"
openfile = open(openfilename,"r")
disp_iso = openfile.readlines()
openfile.close()

v={}
for i in range(len(disp_iso)):
 v[int(disp_iso[i].split()[0])] = (float(disp_iso[i].split()[1]),float(di
sp_iso[i].split()[2]))
for i in xrange(d.getNumberOfDataPoints()):
 Id=d.getFunctionSpace().getReferenceIDFromDataPointNo(i)
 d.setValueOfDataPoint(i, v[Id])
 print i, Id

DISP = FileWriter('disp_check.data')
for i in xrange(d.getNumberOfDataPoints()):
 DISP.write(str(int(i))+' '
   +str(float(d.getTupleForDataPoint(i)[0]))+' '
   +str(float(d.getTupleForDataPoint(i)[1]))+'\n')
DISP.close()

exit()
[bwjliu@int2 Escript_MPI]$ run-escript -p 2 setDisp.py
0 0
1 1
2 2
[bwjliu@int2 Escript_MPI]$ more disp_check.data
0 0.0 0.0
1 -1.80853237932 0.0
2 0.286888846683 -1.40356468889

FileWriter writes data from rank 0 only!

Is it possible to write all node displacements to a file when running simulation with more than 1 processes?

Thank you so much for helping me solve the setValueOfDataPoint-related problem.

Thanks Lutz Gross, that solved my question.