Yade

Run several unloading tests together

Asked by Yu Tian

I want to do a series of biaxial unloading tests.
Before the unloading, I do a biaxial loading test, and save the model in many ".bz2" files when the vertical strain reaches different values. The obtained files are: save0.00010yade.bz2, save0.00020yade.bz2, save0.00030yade.bz2...
Now, I want to unload the sample from these values of vertical strain, until the vertical stress reduces to the confining pressure. Because there are too many .bz2 files, I hope the code could call these files one by one, and after one unloading test is completed, another test would be done automatically.
I first write a code that can do an unloading test once at a time, and it works. Then, I put the code in a “for” loop. However, it seems that the “O.engines” does not start when the unloading test is to be done, because no information about the running time is shown in the terminal window and no file is outputted (I introduce two PyRunners into O.engines).
Could you help me solve this problem?

##########################################
# unicode: UTF-8
# For 2D biaxial simulation
# 19/09/2022

#########################################
### Defining parameters and variables ###
#########################################

#Material constants
FrictionAngle = 35.000000000000000
Damp = 0.5
Porosity = 0.0

#Confining variables
ConfPress = -1.0e5

#Loading control
LoadRate = 0.01

#time calculation
startT = O.time
endT = O.time
timeSpent = endT - startT

########################################
#import necessary packages
from yade import pack,plot,os,timing

#record of e22
recorde = ['0.00010','0.00020','0.00031','0.00041','0.00051','0.00061','0.00071','0.00082','0.00092','0.00102',\
           '0.00112','0.00122','0.00132','0.00143','0.00153','0.00163','0.00173','0.00183','0.00194','0.00204',\
           '0.00214','0.00224','0.00234','0.00245','0.00255','0.00265','0.00275','0.00286','0.00296','0.00306',\
           '0.00316','0.00326','0.00337','0.00347','0.00357','0.00367','0.00377','0.00388','0.00398','0.00408',\
           '0.00418','0.00428','0.00439','0.00449','0.00459','0.00469','0.00479','0.00490','0.00500']

for iii in recorde:
    #folders for post-processing
    try:
        os.mkdir('contact-unload'+str(iii))
    except:
        pass

    #restart
    O.load('save'+str(iii)+'yade.bz2') # obtained from "work calculation_dt=0.00001"
    for k in O.bodies:
        k.state.vel=(0,0,0) # set the initial speed to be 0
        k.state.angVel=(0,0,0)

    triax1=O.engines[4]
    #O.run(100,1)
    def Output():
        e22=-triax1.strain[1]
        #particle info
        g = open('./contact-unload'+str(iii)+'/pInfo'+'{:.5f}'.format(e22)+'.txt', 'w')
        print('particle information at Iter %d' % O.iter, file=g)

        g.write('left right bottom top\n')
        print(O.bodies[0].state.pos[0],O.bodies[1].state.pos[0],O.bodies[2].state.pos[1],O.bodies[3].state.pos[1],file=g)
        g.write('ID x y radius mass fx fy disx disy rotz velx vely omega\n')
        for b in O.bodies:
            if isinstance(b.shape,Sphere):
                print(b.id, b.state.pos[0], b.state.pos[1], b.shape.radius, b.state.mass, O.forces.f(b.id)[0], O.forces.f(b.id)[1], b.state.displ()[0], b.state.displ()[1], b.state.rot()[2], b.state.vel[0], b.state.vel[1], b.state.angVel[2], file=g)

        g.close()

    Output()

    #################################
    #####Defining triaxil engines####
    #################################

    triax1=TriaxialStressController(
        width0=0.9994224497343608,
        height0=1.4991128144248458,
        depth0=0.0400000000000000,
        wall_back_activated = False, #for 2d simulation, z-axis is along the direction of front and back
        wall_front_activated = False,
        thickness = 0.001,
        #maxMultiplier = 1.002,
        internalCompaction = False, # If true the confining pressure is generated by growing particles
        max_vel = 0.1,
        stressMask = 5,
        computeStressStrainInterval = 10,
        goal1 = ConfPress,
        goal2 = LoadRate,
        #goal3 = ConfPress,
    )

    ini_e22 = -triax1.strain[1]
    ini_e22i = -triax1.strain[1]
    kkk = 1
    # print('ini_e22',ini_e22) # The initial strain is always 0, before the engine runs

    O.usesTimeStepper=True
    # O.dt=2e-6
    newton=NewtonIntegrator(damping=Damp)
    setContactFriction(radians(FrictionAngle))
    O.trackEnergy=True
    O.timingEnabled=True

    ###engine
    O.engines=[
        ForceResetter(),
        InsertionSortCollider([Bo1_Sphere_Aabb(),Bo1_Box_Aabb()]),
        InteractionLoop(
         [Ig2_Sphere_Sphere_ScGeom(),Ig2_Box_Sphere_ScGeom()],
         [Ip2_FrictMat_FrictMat_FrictPhys()],
         [Law2_ScGeom_FrictPhys_CundallStrack()]
        ),
        GlobalStiffnessTimeStepper(active=1,timeStepUpdateInterval=100,timestepSafetyCoefficient=0.8),
        triax1,
        newton,
        PyRunner(iterPeriod=100,command='endCheck()',label='check'),
        PyRunner(command='addPlotData()',iterPeriod=10,label='record'),
        # VTKRecorder(iterPeriod=10000,recorders=['all'],fileName='./post/p1-')
        # TriaxialStateRecorder(iterPeriod=2000,file='WallStresses')
    ]

    print('ini_e22',ini_e22)

    # Simulation stop conditions defination
    def endCheck():
        #unb=unbalancedForce()
        e22=-triax1.strain[1]
        if e22<0.0 or triax1.stress(triax1.wall_top_id)[1]>ConfPress:
           O.pause()
           endT = O.time
           timeSpent = endT - startT
           print('iii',iii)
           print('Time using:',timeSpent)
           plot.saveDataTxt('unloadinglog'+str(iii)+'.txt')

    # collect history of data
    def addPlotData():
        global ini_e22i,ini_e22
        global kkk

        e22=-triax1.strain[1]
        if kkk:
            ini_e22i = -triax1.strain[1]
            ini_e22 = -triax1.strain[1]
            kkk = 0
        #print('e22',e22,'ini_e22',ini_e22,'ini_e22i',ini_e22i)

        if abs(ini_e22i) - abs(e22) < 0.00025:
            test = abs(ini_e22) - abs(e22) > 0.000001
        elif abs(ini_e22i) - abs(e22) < 0.08:
            test = abs(ini_e22) - abs(e22) > 0.0001

        if test:
            unb = unbalancedForce()
            mStress = -(triax1.stress(triax1.wall_right_id)[0]+triax1.stress(triax1.wall_top_id)[1])/2.0
            area = (O.bodies[1].state.pos[0]-O.bodies[0].state.pos[0])*(O.bodies[3].state.pos[1]-O.bodies[2].state.pos[1])
            Porosity = 1.0-sum(pi*b.shape.radius**2 for b in O.bodies if isinstance(b.shape,Sphere))/area

            plot.addData(e11=-triax1.strain[0], e22=-triax1.strain[1], e33=-triax1.strain[2],
                 ev=-triax1.strain[0]-triax1.strain[1]-triax1.strain[2],
                 s11=-triax1.stress(triax1.wall_right_id)[0],
                 s22=-triax1.stress(triax1.wall_top_id)[1],
                 s33=-triax1.stress(triax1.wall_front_id)[2],
                 ub=unbalancedForce(),
                 dstress=(-triax1.stress(triax1.wall_top_id)[1])-(-triax1.stress(triax1.wall_right_id)[0]),
                 disx=triax1.width,
                 disy=triax1.height,
                 disz=triax1.depth,
                 i=O.iter,
                 porosity=Porosity,
                 cnumber=avgNumInteractions(),
                 total=O.energy.total(),
                 **O.energy
                 )

            print('axial strain',e22,'unbalanced force:',unb,'mean stress: ',mStress,'s11',-triax1.stress(triax1.wall_right_id)[0],'s22',-triax1.stress(triax1.wall_top_id)[1],'coordination number: ',avgNumInteractions(),'porosity: ',Porosity)
            Output()
            ini_e22 = e22

    if iii==0:
        O.run(); #O.wait()

Question information

Language:
English Edit question
Status:
Solved
For:
Yade Edit question
Assignee:
No assignee Edit question
Solved by:
Jan Stránský
Solved:
Last query:
Last reply:

This question was reopened

Revision history for this message
Jan Stránský (honzik) said : 		#1

Hello,
another option, different from for loop in python, is to run your single-unloading script in a loop or use yade-batch.
Cheers
Jan

Revision history for this message
Jérôme Duriez (jduriez) said : 		#2

Hi,

Regarding the yade-batch idea proposed by Jan, see https://yade-dem.org/doc/user.html#batch-queuing-and-execution-yade-batch where your .table would e.g. look like

saveIdx
0.00010
0.00020
...

and your .py would include something like
O.load('save'+str(saveIdx)+'.yade.bz2') # I added a first "." in the extension of your saved file names, I guess it is good practice

Revision history for this message
Yu Tian (yutian92) said : 		#3

Thanks, Jan and Jerome!
I am trying to run my single-unloading script in a loop first. It partly works, but the obtained stress-strain curves during unloading are not correct. I think I know where the bug is.
I will try the yade-batch later if the first method doesn't work.
Yu

Revision history for this message
Yu Tian (yutian92) said : 		#4

Hi,
I have tried the two methods mentioned by Jan.
First, I package my single-unloading script into a self-defined function, and call it in a loop. I encounter the following problem:
(1) If I put just one .bz2 file into the folder “input” (which is used to contain the .bz2 file to be called during the unloading test), my script works well;
(2) If I put two or more .bz2 files into the folder “input”, only the last .bz2 file can be called, while the following fatal error appears when the first file is called:
<FATAL ERROR> ThreadRunner:35 void yade::ThreadRunner::run(): Exception occured:
Invalid argument

Then, I try the yade-batch. Even through I can get some results, I find the following results on the terminal window:
http://localhost:9080 shows batch summary
#0 (recorde=0.00031) started on Fri Oct 7 15:04:20 2022
#0 (recorde=0.00031) FAILED (exit status 35584), duration 00:05:10, log unloading2.py.recorde=0.00031.log
#1 (recorde=0.00041) started on Fri Oct 7 15:09:31 2022
#1 (recorde=0.00041) FAILED (exit status 35584), duration 00:05:30, log unloading2.py.recorde=0.00041.log
#2 (recorde=0.00051) started on Fri Oct 7 15:15:02 2022
#2 (recorde=0.00051) FAILED (exit status 35584), duration 00:06:53, log unloading2.py.recorde=0.00051.log
All jobs finished, total time 00:17:36

And according to the .log file, the error is:
Error in atexit._run_exitfuncs:
Traceback (most recent call last):
  File "/usr/lib/python3/dist-packages/IPython/core/history.py", line 780, in writeout_cache
    self._writeout_input_cache(conn)
  File "/usr/lib/python3/dist-packages/IPython/core/history.py", line 763, in _writeout_input_cache
    conn.execute("INSERT INTO history VALUES (?, ?, ?, ?)",
sqlite3.ProgrammingError: SQLite objects created in a thread can only be used in that same thread. The object was created in thread id 139639369226048 and this is thread id 139638768092928.

I don’t know why. Could you help me?

Here is my script of running yade in a loop:
#import necessary packages
from yade import pack,plot,os,timing
import re,sys

#Material constants
FrictionAngle = 35.000000000000000
Damp = 0.5
Porosity = 0.0

#Confining variables
ConfPress = -1.0e5

#Loading control
LoadRate = 0.01

#time calculation
startT = O.time
endT = O.time
timeSpent = endT - startT

#sample size after isotropic consolidation, find them from 'loadinglog.txt.bz2'
Width00=0.9994224138086599
Height00=1.4991129353301802
Depth00=0.0400000000000000

#folders for post-processing
try:
    os.mkdir('output')
except:
    pass

########################################
# get the name of all the .bz2 files
def get_filename():
    files = os.listdir("/home/yu/Downloads/code/unload/input")
    files.sort()
    return files

# output the information about contacts and particles
def Output():
    e22=-triax1.strain[1]
    #particle info
    g = open(path2+'/{:.5f}'.format(e22)+'pInfo-unload.txt', 'w')
    print('particle information at Iter %d' % O.iter, file=g)
    g.write('left right bottom top\n')
    print(O.bodies[0].state.pos[0],O.bodies[1].state.pos[0],O.bodies[2].state.pos[1],O.bodies[3].state.pos[1],file=g)

    g.write('ID x y radius mass fx fy disx disy rotz velx vely omega\n')
    for b in O.bodies:
        if isinstance(b.shape,Sphere):
            print(b.id, b.state.pos[0], b.state.pos[1], b.shape.radius, b.state.mass, O.forces.f(b.id)[0], O.forces.f(b.id)[1], b.state.displ()[0], b.state.displ()[1], b.state.rot()[2], b.state.vel[0], b.state.vel[1], b.state.angVel[2], file=g)
    g.close()

# Simulation stop conditions definition
def endCheck():
    e22=-triax1.strain[1]
    if e22<0.0 or triax1.stress(triax1.wall_top_id)[1]>ConfPress:
        O.pause()
        endT = O.time
        timeSpent = endT - startT
        print('unload from e22=',fileid)
        print('Time using:',timeSpent)
        #sys.exit()

# collect history of data
def addPlotData():
    global ini_e22
    e22=-triax1.strain[1]
    #print('e22=',e22,'ini_e22=',ini_e22)

    if float(fileid) - abs(e22) < 0.0002:
        test = abs(ini_e22) - abs(e22) > 0.00001 and abs(e22) > 1e-8
        # e22 is always 0 at the first iteration after O.load
    elif float(fileid) - abs(e22) < 0.5:
        test = abs(ini_e22) - abs(e22) > 0.0001 and abs(e22) > 1e-8

    if test:
        unb = unbalancedForce()
        mStress = -(triax1.stress(triax1.wall_right_id)[0]+triax1.stress(triax1.wall_top_id)[1])/2.0
        area = (O.bodies[1].state.pos[0]-O.bodies[0].state.pos[0])*(O.bodies[3].state.pos[1]-O.bodies[2].state.pos[1])
        Porosity = 1.0-sum(pi*b.shape.radius**2 for b in O.bodies if isinstance(b.shape,Sphere))/area

        plot.addData(e11=-triax1.strain[0], e22=-triax1.strain[1], e33=-triax1.strain[2],
             ev=-triax1.strain[0]-triax1.strain[1]-triax1.strain[2],
             s11=-triax1.stress(triax1.wall_right_id)[0],
             s22=-triax1.stress(triax1.wall_top_id)[1],
             s33=-triax1.stress(triax1.wall_front_id)[2],
             ub=unbalancedForce(),
             dstress=(-triax1.stress(triax1.wall_top_id)[1])-(-triax1.stress(triax1.wall_right_id)[0]),
             disx=triax1.width,
             disy=triax1.height,
             disz=triax1.depth,
             i=O.iter,
             porosity=Porosity,
             cnumber=avgNumInteractions(),
             total=O.energy.total(),
             **O.energy
             )
        print('axial strain',e22,'unbalanced force:',unb,'mean stress: ',mStress,'s11',-triax1.stress(triax1.wall_right_id)[0],'s22',-triax1.stress(triax1.wall_top_id)[1],'coordination number: ',avgNumInteractions(),'porosity: ',Porosity)
        Output()
        plot.saveDataTxt('unloadinglog'+fileid+'.txt')
        ini_e22 = e22

########################################
# unloading
def Yade_unload(filename):
    global triax1,fileid,ini_e22,path2

    # restart
    path1='/home/yu/Downloads/code/unload/input/'+filename
    O.load(path1)
    for k in O.bodies:
        k.state.vel=(0,0,0) # set the initial speed to be 0
        k.state.angVel=(0,0,0)

    # get epsilon2 from the filename
    fileid = re.findall(r"\d+\.?\d*",filename)[0]
    ini_e22 = float(fileid)
    # establish a subfolder for output
    path2='/home/yu/Downloads/code/unload/output/'+fileid+'-unload'
    if not os.path.exists(path2):
        os.makedirs(path2)

    triax1=O.engines[4]

    #################################
    #####Defining triaxil engines####
    #################################
    triax1=TriaxialStressController(
        width0=Width00,
        height0=Height00,
        depth0=Depth00,
        wall_back_activated = False, #for 2d simulation, z-axis is along the direction of front and back
        wall_front_activated = False,
        thickness = 0.001,
        #maxMultiplier = 1.002,
        internalCompaction = False, # If true the confining pressure is generated by growing particles
        max_vel = 0.1,
        stressMask = 5,
        computeStressStrainInterval = 10,

        goal1 = ConfPress,
        goal2 = LoadRate,
        #goal3 = ConfPress,

    )

    O.usesTimeStepper=True
    #O.dt=2e-6
    newton=NewtonIntegrator(damping=Damp)
    setContactFriction(radians(FrictionAngle))
    O.trackEnergy=True
    O.timingEnabled=True

    ###engine
    O.engines=[
        ForceResetter(),
        InsertionSortCollider([Bo1_Sphere_Aabb(),Bo1_Box_Aabb()]),
        InteractionLoop(
         [Ig2_Sphere_Sphere_ScGeom(),Ig2_Box_Sphere_ScGeom()],
         [Ip2_FrictMat_FrictMat_FrictPhys()],
         [Law2_ScGeom_FrictPhys_CundallStrack()]
        ),
        GlobalStiffnessTimeStepper(active=1,timeStepUpdateInterval=100,timestepSafetyCoefficient=0.8),
        triax1,
        newton,
        PyRunner(realPeriod=100,command='endCheck()',label='check'),
        PyRunner(command='addPlotData()',iterPeriod=10,label='record')
    ]

    O.run()

# main program
files = get_filename()
for iii in files:
    Yade_unload(iii)
    print('{} is called'.format(iii))

#######################################
Here is my script of yade-batch:
#import necessary packages
readParamsFromTable(recorde = 0.00031)
from yade.params import table
from yade import pack,plot,os,timing

#Material constants
FrictionAngle = 35.000000000000000
Damp = 0.5
Porosity = 0.0

#Confining variables
ConfPress = -1.0e5

#Loading control
LoadRate = 0.01

#time calculation
startT = O.time
endT = O.time
timeSpent = endT - startT

#sample size after isotropic consolidation, find them from 'loadinglog.txt.bz2'
Width00=0.9994224138086599
Height00=1.4991129353301802
Depth00=0.0400000000000000

#folders for post-processing
try:
    os.mkdir(str(table.recorde)+'-unload')
except:
    pass

#restart
O.load(str(table.recorde)+'-load.bz2')
for k in O.bodies:
    k.state.vel=(0,0,0) # set the initial speed to be 0
    k.state.angVel=(0,0,0)
ini_e22 = table.recorde

triax1=O.engines[4]
def Output():
    e22=-triax1.strain[1]
    #contact info
    f = open(str(table.recorde)+'-unload/{:.5f}'.format(e22)+'cInfo-unload.txt', 'w')
    print('contact information at Iter %d' % O.iter, file=f)

    f.write('left right bottom top\n')
    print(O.bodies[0].state.pos[0],O.bodies[1].state.pos[0],O.bodies[2].state.pos[1],O.bodies[3].state.pos[1],file=f)

    f.write('ctype id1 id2 nfx nfy tfx tfy xc yc rad1 x1 y1 rot1 rad2 x2 y2 rot2\n')
    for k in O.interactions:
        if isinstance(O.bodies[k.id1].shape,Sphere) and isinstance(O.bodies[k.id2].shape,Box):
            print('01',k.id1,k.id2,k.phys.normalForce[0], k.phys.normalForce[1], k.phys.shearForce[0], k.phys.shearForce[1], k.geom.contactPoint[0], k.geom.contactPoint[1], O.bodies[k.id1].shape.radius, O.bodies[k.id1].state.pos[0], O.bodies[k.id1].state.pos[1], O.bodies[k.id1].state.rot()[2], '111', O.bodies[k.id2].state.pos[0], O.bodies[k.id2].state.pos[1], O.bodies[k.id2].state.rot()[2], file=f)
        elif isinstance(O.bodies[k.id1].shape,Box) and isinstance(O.bodies[k.id2].shape,Sphere):
            print('10',k.id1,k.id2,k.phys.normalForce[0], k.phys.normalForce[1], k.phys.shearForce[0], k.phys.shearForce[1], k.geom.contactPoint[0], k.geom.contactPoint[1], '111', O.bodies[k.id1].state.pos[0], O.bodies[k.id1].state.pos[1], O.bodies[k.id1].state.rot()[2], O.bodies[k.id2].shape.radius, O.bodies[k.id2].state.pos[0], O.bodies[k.id2].state.pos[1], O.bodies[k.id2].state.rot()[2], file=f)
        elif isinstance(O.bodies[k.id1].shape,Sphere) and isinstance(O.bodies[k.id2].shape,Sphere):
            print('00',k.id1,k.id2,k.phys.normalForce[0], k.phys.normalForce[1], k.phys.shearForce[0], k.phys.shearForce[1], k.geom.contactPoint[0], k.geom.contactPoint[1], O.bodies[k.id1].shape.radius, O.bodies[k.id1].state.pos[0], O.bodies[k.id1].state.pos[1], O.bodies[k.id1].state.rot()[2], O.bodies[k.id2].shape.radius, O.bodies[k.id2].state.pos[0], O.bodies[k.id2].state.pos[1], O.bodies[k.id2].state.rot()[2], file=f)

    f.close()

    #particle info
    g = open(str(table.recorde)+'-unload/{:.5f}'.format(e22)+'pInfo-unload.txt', 'w')
    print('particle information at Iter %d' % O.iter, file=g)

    g.write('left right bottom top\n')
    print(O.bodies[0].state.pos[0],O.bodies[1].state.pos[0],O.bodies[2].state.pos[1],O.bodies[3].state.pos[1],file=g)

    g.write('ID x y radius mass fx fy disx disy rotz velx vely omega\n')
    for b in O.bodies:
        if isinstance(b.shape,Sphere):
            print(b.id,b.state.pos[0],b.state.pos[1],b.shape.radius,b.state.mass,O.forces.f(b.id)[0],O.forces.f(b.id)[1],b.state.displ()[0],b.state.displ()[1],b.state.rot()[2],b.state.vel[0],b.state.vel[1],b.state.angVel[2],file=g)

    g.close()

#################################
#####Defining triaxil engines####
#################################

triax1=TriaxialStressController(
    width0=Width00,
    height0=Height00,
    depth0=Depth00,
    wall_back_activated = False, #for 2d simulation, z-axis is along the direction of front and back
    wall_front_activated = False,
    thickness = 0.001,
    #maxMultiplier = 1.002,
    internalCompaction = False, # If true the confining pressure is generated by growing particles
    max_vel = 0.1,
    stressMask = 5,
    computeStressStrainInterval = 10,

    goal1 = ConfPress,
    goal2 = LoadRate,
    #goal3 = ConfPress,

)

O.usesTimeStepper=True
#O.dt=2e-6
newton=NewtonIntegrator(damping=Damp)
setContactFriction(radians(FrictionAngle))
O.trackEnergy=True
O.timingEnabled=True

###engine
O.engines=[
    ForceResetter(),
    InsertionSortCollider([Bo1_Sphere_Aabb(),Bo1_Box_Aabb()]),
    InteractionLoop(
     [Ig2_Sphere_Sphere_ScGeom(),Ig2_Box_Sphere_ScGeom()],
     [Ip2_FrictMat_FrictMat_FrictPhys()],
     [Law2_ScGeom_FrictPhys_CundallStrack()]
    ),
    GlobalStiffnessTimeStepper(active=1,timeStepUpdateInterval=100,timestepSafetyCoefficient=0.8),
    triax1,
    newton,
    PyRunner(realPeriod=100,command='endCheck()',label='check'),
    PyRunner(command='addPlotData()',iterPeriod=10,label='record')
]

# Simulation stop conditions defination
def endCheck():
    e22=-triax1.strain[1]
    s22=-triax1.stress(triax1.wall_top_id)[1]
    #print('e22=',e22,'s22=',s22)
    if e22<0.0 or triax1.stress(triax1.wall_top_id)[1]>ConfPress:
       O.pause()
       endT = O.time
       timeSpent = endT - startT
       print('unload from e22=',table.recorde)
       print('Time using:',timeSpent)
       import sys
       sys.exit()
       #O.exitNoBacktrace()

# collect history of data
def addPlotData():
    global ini_e22
    e22=-triax1.strain[1]
    #print('e22=',e22,'ini_e22=',ini_e22)

    if table.recorde - abs(e22) < 0.0002:
        test = abs(ini_e22) - abs(e22) > 0.00001 and abs(e22) > 1e-8
        # e22 is always 0 at the first iteration after O.load
    elif table.recorde - abs(e22) < 0.5:
        test = abs(ini_e22) - abs(e22) > 0.0001 and abs(e22) > 1e-8

    if test:
        unb = unbalancedForce()
        mStress = -(triax1.stress(triax1.wall_right_id)[0]+triax1.stress(triax1.wall_top_id)[1])/2.0
        area = (O.bodies[1].state.pos[0]-O.bodies[0].state.pos[0])*(O.bodies[3].state.pos[1]-O.bodies[2].state.pos[1])
        Porosity = 1.0-sum(pi*b.shape.radius**2 for b in O.bodies if isinstance(b.shape,Sphere))/area

        plot.addData(e11=-triax1.strain[0], e22=-triax1.strain[1], e33=-triax1.strain[2],
                 ev=-triax1.strain[0]-triax1.strain[1]-triax1.strain[2],
                 s11=-triax1.stress(triax1.wall_right_id)[0],
                 s22=-triax1.stress(triax1.wall_top_id)[1],
                 s33=-triax1.stress(triax1.wall_front_id)[2],
                 ub=unbalancedForce(),
                 dstress=(-triax1.stress(triax1.wall_top_id)[1])-(-triax1.stress(triax1.wall_right_id)[0]),
                 disx=triax1.width,
                 disy=triax1.height,
                 disz=triax1.depth,
                 i=O.iter,
                 porosity=Porosity,
                 cnumber=avgNumInteractions(),
                 total=O.energy.total(),
                 **O.energy
                 )
        print('axial strain',e22,'unbalanced force:',unb,'mean stress: ',mStress,'s11',-triax1.stress(triax1.wall_right_id)[0],'s22',-triax1.stress(triax1.wall_top_id)[1],'coordination number: ',avgNumInteractions(),'porosity: ',Porosity)
        plot.saveDataTxt('unloadinglog'+str(table.recorde)+'.txt')
        Output()
        ini_e22 = e22

O.run()
waitIfBatch()

Revision history for this message
Best Jan Stránský (honzik) said : 		#5

> (2) If I put two or more .bz2 files into the folder “input”, only the last .bz2 file can be called
>
> def Yade_unload(filename):
> ...
> O.run()

After O.run(), next commands (next loop) are executed immediately, therefore you get only the last run.
Add O.wait() or some other stop/wait condition.

> I try the yade-batch

with what command you try it?

> And according to the .log file, the error is:

seems like some IPython error, I have no idea what it is or how to fix it

Cheers
Jan

Revision history for this message
Jérôme Duriez (jduriez) said : 		#6

Hi, your error look to relate with an attempted use of parallel computing, which would be impossible. Did you use -j someNumber option when launching your yade simulation combining all unloads in one ?..

Generally speaking, I would advise you first try your intended "running several tests together from previous saves" workflow on much simpler YADE simulations and scripts (even if it has nothing to do with your final goal, e.g. just 1 or 2 spheres). It would much increase your / our chances to see exactly where the problem is.

Revision history for this message
Yu Tian (yutian92) said : 		#7

Thanks Jan Stránský, that solved my question.