yade-users team mailing list archive

Thread
Date
Re: [Question #703247]: Run several unloading tests together

To: yade-users@xxxxxxxxxxxxxxxxxxx
From: Yu Tian <question703247@xxxxxxxxxxxxxxxxxxxxx>
Date: Sun, 09 Oct 2022 13:30:39 -0000
Reply-to: question703247@xxxxxxxxxxxxxxxxxxxxx
Sender: noreply@xxxxxxxxxxxxx
Question #703247 on Yade changed:
https://answers.launchpad.net/yade/+question/703247

    Status: Solved => Open

Yu Tian is still having a problem:
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()

-- 
You received this question notification because your team yade-users is
an answer contact for Yade.