Convert script to Openmpi

Dear all,

I am new on running Yade on multiple nodes using openmpi.
Still checking how to convert the code below to run in mpi.

Runnning the code below I get the messsage:
Friction: 1 porosity: 1.0 None
Master: single-core, fall back to O.run() 

Probably I am not using O.run and mpirun correctly.
If anyone with more experimence see what I am duing wrong, please let me know.

Cheers,

#!/usr/bin/python
# -*- encoding=utf-8 -*-
#*************************************************************************
# Copyright (C) 2010 by Bruno Chareyre *
# bruno.chareyre_at_grenoble-inp.fr *
# *
# This program is free software; it is licensed under the terms of the *
# GNU General Public License v2 or later. See file LICENSE for details. *
#*************************************************************************/
import os
from yade import mpy as mp
from yade import pack
from yade import bodiesHandling
from yade import export
from yade import utils
from yade import ymport
import math

############################################
### DEFINING VARIABLES AND MATERIALS ###
############################################

# The following 5 lines will be used later for batch execution
nRead=readParamsFromTable(
    num_spheres=3000,# number of spheres
    compFricDegree = 1, # contact friction during the confining phase
    key='_triax_base_', # put you simulation's name here
    unknownOk=True
)
from yade.params import table

num_spheres=table.num_spheres# number of spheres
key=table.key

targetPorosity = 0.35 #the porosity we want for the packing
compFricDegree = table.compFricDegree # initial contact friction during the confining phase (will be decreased during the REFD compaction process)
finalFricDegree = 30 # contact friction during the deviatoric loading
rate=0 # loading rate (strain rate)
damp=0.8 # damping coefficient
stabilityThreshold=0.01 # we test unbalancedForce against this value in different loops (see below)
#2e4+70e4medio 1e4+70e4bom 1e4+60e4bom 3e4+90e4+w3,1,-1-the best
young=20e5 # contact stiffness200e4
young2=20e5
youngcoat=20e5
bondstr=0.3e3#2e7
bondstr2=0.3e3
bondstrcoat=1e3

## create materials for spheres and plates
mat=O.materials.append(JCFpmMat(type=1,young=young,poisson=0.3,frictionAngle=radians(compFricDegree),density=2000,tensileStrength=bondstr,cohesion=bondstr,jointNormalStiffness=0,jointShearStiffness=0,jointCohesion=bondstr,jointFrictionAngle=radians(0),jointDilationAngle=0.0,label='spheres'))
O.materials.append(JCFpmMat(type=1,young=20e7,poisson=0.3,frictionAngle=radians(0),density=2600,tensileStrength=0,cohesion=0,jointNormalStiffness=0,jointShearStiffness=0,jointCohesion=0,jointFrictionAngle=radians(0),jointDilationAngle=0.0,label='walls'))
O.materials.append(JCFpmMat(type=1,young=youngcoat,poisson=0.3,frictionAngle=radians(1),density=1500,tensileStrength=bondstrcoat,cohesion=bondstrcoat,jointNormalStiffness=0,jointShearStiffness=0,jointCohesion=bondstrcoat,jointFrictionAngle=radians(0),jointDilationAngle=0.0,label='spherescoat'))

## create walls around the packing

mn,mx=Vector3(0,0,0),Vector3(0.0101,0.0101,0.0101)
mnbox,mxbox=Vector3(0,0,0),Vector3(0.0101,0.0115,0.0101)
walls=aabbWalls([mnbox,mxbox],thickness=0,material='walls')
wallIds=O.bodies.append(walls)

O.bodies.append(ymport.textExt("matrix_vtest.txt", format='x_y_z_r', shift=Vector3(0,0,0), scale=1.0,material='spherescoat',color=(0,1,1)))

################Particle substitution by large aggregate######################################################################

bodid=[]
a=[]
for b in O.bodies:
 if b and isinstance(b.shape,Sphere):
# print (b.shape.radius)
  if b.shape.radius==0.0005:
   bodid.append(b.id)
   a.append(b.state.pos)

i=0
for p in bodid:
 t=a[i]
 f1=O.bodies.append(ymport.textExt("agg5e4_5e5.txt", format='x_y_z_r', shift=t-Vector3(0,0,0), scale=1.0,material='spheres',color=(0,1,1)))
 O.bodies.erase(bodid[i])
 i=i+1

bodidd=[]
aa=[]
for bb in O.bodies:# in sp:
 if bb and isinstance(bb.shape,Sphere):
# print (bb.shape.radius)
  if bb.shape.radius==0.0003:
   bodidd.append(bb.id)
   aa.append(bb.state.pos)

ii=0
for pp in bodidd:
 tt=aa[ii]
 f2=O.bodies.append(ymport.textExt("agg3e4_3e5.txt", format='x_y_z_r', shift=tt-Vector3(0,0,0), scale=1.0,material='spheres',color=(0,1,1)))
 O.bodies.erase(bodidd[ii])
 ii=ii+1

bodiddd=[]
aaa=[]
for bbb in O.bodies:# in sp:
 if bbb and isinstance(bbb.shape,Sphere):
# print (bbb.shape.radius)
  if bbb.shape.radius==0.0002:
   bodiddd.append(bbb.id)
   aaa.append(bbb.state.pos)

iii=0
for ppp in bodiddd:
 ttt=aaa[iii]
 f3=O.bodies.append(ymport.textExt("agg2e4_2e5.txt", format='x_y_z_r', shift=ttt-Vector3(0,0,0), scale=1.0,material='spheres',color=(0,1,1)))
 O.bodies.erase(bodiddd[iii])
 iii=iii+1

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

#====================================================================================================================================================================

############################
### DEFINING ENGINES ###
############################

triax=TriaxialStressController(
    ## TriaxialStressController will be used to control stress and strain. It controls particles size and plates positions.
    ## this control of boundary conditions was used for instance in http://dx.doi.org/10.1016/j.ijengsci.2008.07.002
    maxMultiplier=1.+2e4/young, # spheres growing factor (fast growth)
    finalMaxMultiplier=1.+2e3/young, # spheres growing factor (slow growth)
    thickness = 0,
    ## switch stress/strain control using a bitmask. What is a bitmask, huh?!
    ## Say x=1 if stess is controlled on x, else x=0. Same for for y and z, which are 1 or 0.
    ## Then an integer uniquely defining the combination of all these tests is: mask = x*1 + y*2 + z*4
    ## to put it differently, the mask is the integer whose binary representation is xyz, i.e.
    ## "100" (1) means "x", "110" (3) means "x and y", "111" (7) means "x and y and z", etc.
    stressMask = 0,
    internalCompaction=False, # If true the confining pressure is generated by growing particles
    wall_front_activated=True,
    wall_back_activated=True,
    wall_top_activated=True,
    wall_bottom_activated=True,
    wall_left_activated=True,
    wall_right_activated=True,
    goal1=-5,
    goal2=-10,
    goal3=-5,
)

newton=NewtonIntegrator(damping=damp)

O.engines=[

 ForceResetter(),
 InsertionSortCollider([Bo1_Sphere_Aabb(),Bo1_Box_Aabb()]),
 InteractionLoop(
  [Ig2_Sphere_Sphere_ScGeom(),Ig2_Box_Sphere_ScGeom()],
  [Ip2_JCFpmMat_JCFpmMat_JCFpmPhys(cohesiveTresholdIteration=-1,label='Physicspheres')],#,xSectionWeibullShapeParameter=1.5, xSectionWeibullScaleParameter=1
  [Law2_ScGeom_JCFpmPhys_JointedCohesiveFrictionalPM(smoothJoint=False)]
 ),
 GlobalStiffnessTimeStepper(active=1,timeStepUpdateInterval=100,timestepSafetyCoefficient=0.3),
 triax,
# VTKRecorder(iterPeriod=2000,recorders=['all','cracks'],fileName="/home/user/Área de Trabalho/Paper_Biopore/SandBox_Laser/Video/Videotest"),
# newton
# NewtonIntegrator(damping=0.4,gravity=[0,0,0]),
 newton

]

#Display spheres with 2 colors for seeing rotations better
#Gl1_Sphere.stripes=0
#if nRead==0: yade.qt.Controller(), yade.qt.View()

###################################################
### REACHING A SPECIFIED POROSITY PRECISELY ###
###################################################

## We will reach a prescribed value of porosity with the REFD algorithm
## (see http://dx.doi.org/10.2516/ogst/2012032 and
## http://www.geosyntheticssociety.org/Resources/Archive/GI/src/V9I2/GI-V9-N2-Paper1.pdf)

import sys #this is only for the flush() below

#poro=utils.voxelPorosityTriaxial(triax, 200, offset=0.0018)
while triax.porosity>targetPorosity:
    # we decrease friction value and apply it to all the bodies and contacts
# compFricDegree = 0.95*compFricDegree
    setContactFriction(radians(compFricDegree))
    print ("\r Friction: ",compFricDegree," porosity:",triax.porosity,
    sys.stdout.flush())
    # while we run steps, triax will tend to grow particles as the packing
    # keeps shrinking as a consequence of decreasing friction. Consequently
    # porosity will decrease
    mp.mpirun(500,1)

#O.save('compactedState'+key+'.yade.gz')
mp.mprint("### Compacted state saved ###")
#print(triax.stress(3)[1])

##############################
### DEVIATORIC LOADING ###
##############################

#We move to deviatoric loading, let us turn internal compaction off to keep particles sizes constant
#triax.internalCompaction=False

# Change contact friction (remember that decreasing it would generate instantaneous instabilities)
setContactFriction(radians(finalFricDegree))
#O.materials.append(JCFpmMat(type=1,young=young,poisson=0.3,frictionAngle=radians(finalFricDegree),density=2000,tensileStrength=bondstr2,cohesion=bondstr2,jointNormalStiffness=0,jointShearStiffness=0,jointCohesion=bondstr2,jointFrictionAngle=radians(0),jointDilationAngle=0.0,label='spheres'))

#=======================================

#set stress control on x and z, we will impose strain rate on y
triax.stressMask = 2
triax.wall_bottom_activated=0
#now goal2 is the target strain rate
triax.goal1=rate#triax.stress(1)[0]
triax.goal3=rate#triax.stress(5)[2]
triax.goal2=triax.stress(3)[1]

##Save temporary state in live memory. This state will be reloaded from the interface with the "reload" button.
#O.saveTmp()

#####################################################
### Example of how to record and plot data ###
#####################################################

#from yade import plot
from yade import plot
mp.mpirun(10,True)

#strain is logarithmic strain or true strain which is ls=(ln1+e) where e=dl/L (strain)
ei0=-triax.strain[0];ei1=-triax.strain[1];ei2=-triax.strain[2]
si0=-triax.stress(0)[0];si1=-triax.stress(2)[1];si2=-triax.stress(4)[2]

## a function saving variables
def history():
    plot.addData(e11=-triax.strain[0]-ei0, e22=-triax.strain[1]-ei1, e33=-triax.strain[2]-ei2,
            ev=-triax.strain[0]-triax.strain[1]-triax.strain[2],
            s11=-triax.stress(triax.wall_right_id)[0]-si0,
            s22=-triax.stress(triax.wall_top_id)[1]-si1,
            s33=-triax.stress(triax.wall_front_id)[2]-si2,
            e=math.exp(-triax.strain[1]-ei1)-1,
            pc=-unsat.bndCondValue[2],
            sw=unsat.getSaturation(False),
            z1=O.bodies[3].state.pos[1],
            i=O.iter)

if 1:
  ## include a periodic engine calling that function in the simulation loop
  O.engines=O.engines[0:5]+[PyRunner(iterPeriod=20,command='history()',label='recorder')]+O.engines[5:7]

#plot.plots={'pc':('sw',None,'e22')}
#plot.plot()
#print("test",math.exp(2))

#######################################################
## Drainage Test under oedometer conditions ###
#######################################################
##Instantiate a two-phase engine
unsat=TwoPhaseFlowEngine()
#meanDiameter=(O.bodies[-1].shape.radius + O.bodies[6].shape.radius) / 2.

##set boundary conditions, the drainage is controlled by decreasing W-phase pressure and keeping NW-phase pressure constant
unsat.bndCondIsPressure=[0,0,1,1,0,0]
unsat.bndCondIsWaterReservoir=[0,0,1,0,0,0]
unsat.bndCondValue=[0,0,-1e8,0,0,0]
unsat.isPhaseTrapped=True #the W-phase can be disconnected from its reservoir
unsat.initialization()
unsat.surfaceTension = 0.0728

#start invasion, the data of normalized pc-sw-strain will be written into pcSwStrain.txt

f5=open('SwPc3.txt',"w")

ts=O.dt
pgstep= 45000000*ts #30Pa/s
print (pgstep)
pgmax= 9000#9316 #Pa

for pg in arange(1.0e-8,pgmax,pgstep):
  unsat.bndCondValue=[0,0,(-1.0)*pg,0,0,0]
  unsat.invasion()
  unsat.computeCapillaryForce()
  unsat.meshUpdateInterval=500
  unsat.defTolerance=-1
  unsat.updateTriangulation=True
  print(unsat.getSaturation(False),pg,-triax.strain[1])

  for b in O.bodies:
    O.forces.setPermF(b.id, unsat.fluidForce(b.id))

  while 1:
    mp.mpirun(100,True)
    unb=unbalancedForce()
    if unb<0.1:
      break
  f5.write(str(pg)+" "+str(unsat.getSaturation(False))+" "+str(triax.strain[1])+"\n")
f5.close()
plot.plots={'pc':('sw',None,'e22')}
plot.plot(noShow=True).savefig('Fig.png')

mp.ERASE_REMOTE_MASTER = True #keep remote bodies in master?
mp.DOMAIN_DECOMPOSITION= True #automatic splitting/domain decomposition
mp.mpirun(20000) #passive mode run
mp.MERGE_W_INTERACTIONS = False
#mp.mpirun(NSTEPS,numThreads,withMerge=True) # interactive run, numThreads is the number of workers to be initialized, see below for withMerge explanation.
mp.mergeScene() #merge scene after run.