Load and run new simulation

Here is code for saved simulation, which loaded into the new simulation:
from __future__ import print_function
from yade import utils, plot, timing
from yade import pack
import pandas as pd
import numpy as np
from PIL import Image
from yade import pack, export
from scipy.interpolate import interp1d
readParamsFromTable(k1=2e5,kp=2e6)
from yade.params.table import *
from scipy.integrate import odeint
import matplotlib.pyplot as plt
import os
import csv
from matplotlib.pyplot import figure
from scipy.interpolate import interp1d
from pylab import *
from scipy.optimize import curve_fit

o = Omega()
save=1
# Physical parameters
fr = 0.41
rho = 1561
Diameter = 0.0012/2
D=Diameter
r1 = Diameter/2
#r2 = Diameter/2
k1 = 1.26e5
kp = 12*k1
kc = k1 * 0.1
ks = k1 * 0.1
DeltaPMax = Diameter/3.0
Chi1 = 0.34
PhiF1=0.999
#PhiF1 = DeltaPMax*(kp-k1)*(r1+r2)/(kp*2*r1*r2)

if save==1:
    o.dt = 1.0e-7
else:
    o.dt=1.0e-5

particleMass = 4.0/3.0*math.pi*r1*r1*r1*rho
Vi1 = math.sqrt(k1/particleMass)*DeltaPMax*Chi1
m_tot=0.0005
#no_p=m_tot/particleMass
no_p=1000

Tab_rad=0.005
Cyl_height=0.015
cross_area=math.pi*(Tab_rad**2)

Comp_press= 1.2e8
Comp_force=Comp_press*cross_area

##Single particle swelling model
def model(r,t,Q_max,rho_t,rho_w,r_0,D):
    Q=((rho_w*(r**3))/(rho_t*(r_0**3)))-(rho_w/rho_t)+1;
    drdt =((D*rho_t)/(r*rho_w))*((Q_max-Q)/Q);
    return drdt
P=[1.45,rho,1000,396.39e-12]

#*************************************
# Add material
mat1 = o.materials.append(LudingMat(frictionAngle=fr, density=rho, k1=k1, kp=kp, ks=ks, kc=kc, PhiF=PhiF1, G0 = 0.0))

# Spheres for compression and walls
sp=pack.SpherePack()
sp.makeCloud((-5.0*Diameter,-5.0*Diameter,-12*Diameter),(5.0*Diameter,5.0*Diameter,7.0*Diameter), rMean=Diameter/2.0,rRelFuzz=0.18,num=1000)
sp.toSimulation(material=mat1)
walls=o.bodies.append(yade.geom.facetCylinder((0,0,0),radius=Tab_rad,height=Cyl_height,segmentsNumber=20,wallMask=6,material=mat1))

r_save=[]
radius=[]
initial_save=[]
size_save=[]
radius.append(0)
swell_t=np.zeros((1,no_p))
i=0
for b in O.bodies:
    if isinstance(b.shape, Sphere):
        radius.append(b.shape.radius)
r_save.append(radius)

# Add engines
o.engines = [
  ForceResetter(),
  InsertionSortCollider([Bo1_Sphere_Aabb(aabbEnlargeFactor=1.05),
                         Bo1_Wall_Aabb(),
                         Bo1_Facet_Aabb()
                         ]),
  InteractionLoop(
    [Ig2_Sphere_Sphere_ScGeom(interactionDetectionFactor=1.05),
     Ig2_Facet_Sphere_ScGeom(),
     Ig2_Wall_Sphere_ScGeom()],
    [Ip2_LudingMat_LudingMat_LudingPhys()],
    [Law2_ScGeom_LudingPhys_Basic()]
  ),
  NewtonIntegrator(damping=0.1, gravity=[0, 0, -9.81]),
  PyRunner(command='checkForce()', realPeriod=1, label="fCheck"),
  #DeformControl(label="DefControl")
]

def checkForce():
    # at the very start, unbalanced force can be low as there is only few
    # contacts, but it does not mean the packing is stable
    if O.iter < 400000:
        return
    # the rest will be run only if unbalanced is < .1 (stabilized packing)
    timing.reset()
    if unbalancedForce() > 1:
        return
    # add plate at upper box side

    highSphere = 0.0
    for b in O.bodies:
        if highSphere < b.state.pos[2] and isinstance(b.shape, Sphere):
            highSphere = b.state.pos[2]
        else:
            pass

    #o.bodies.append(wall(highSphere, axis=2, sense=-1, material=mat1))
    #global plate
    #plate = o.bodies[-1] # the last particles is the plat
    #plate.state.vel = (0, 0, -2)
    global upper_punch
    upper_punch=O.bodies.append(geom.facetCylinder((0,0,((-Cyl_height/2)+0.0001)+utils.aabbDim()[2]),Tab_rad-.00001,0,segmentsNumber=50,wallMask=1))
    for i in upper_punch:
        body= O.bodies[i]
        body.state.vel = (0,0,-2)
    fCheck.command = 'unloadPlate()'

def unloadPlate():
    force_up=0
    for i in upper_punch:
        body= O.bodies[i]
        force_up=force_up+abs(O.forces.f(body.id)[2])
    # if the force on plate exceeds maximum load, start unloading
    # if abs(O.forces.f(plate.id)[2]) > 5e-2:
    if force_up > Comp_force:
        for i in upper_punch:
            body= O.bodies[i]
            body.state.vel = (0,0,2)
        #plate.state.vel *= -1
        # next time, do not call this function anymore, but the next one
        # (stopUnloading) instead
        fCheck.command = 'stopUnloading()'

def stopUnloading():
    for i in upper_punch:
        body=O.bodies[i]
        pos_up=body.state.pos
    if pos_up[2] > Cyl_height/2:
        for i in upper_punch:
            body= O.bodies[i]
            body.state.vel = (0,0,2)
        #plate.state.vel = (0, 0, 0)
        initial_save.append(utils.aabbExtrema()[1][2]+r1)
        initial_save.append(O.iter)
        #o.engines = o.engines+[PyRunner(command='ParticleSwelling()', iterPeriod=100000)]
        #for b in O.bodies:
            #if isinstance(b.shape,Wall): O.bodies.erase(b.id) # erase punch
        for j in upper_punch: O.bodies.erase(j)
        #o.pause()
        #if save==1:
            #utils.saveVars('somethingtest',plate=plate)
            #o.save('Tablet_swelling.xml')
        fCheck.command = 'ParticleSwelling()'

def ParticleSwelling():
    if save==1:
        utils.saveVars('somethingtest',plate=upper_punch)
        o.save('Tablet_swelling.xml')
        o.pause()
O.run();
O.wait()

Hi,

> I running a new simulation on loaded data.

no. With O.load, you just continue where the O.saved simulation ended.

If your new simulation is not too much related to the old one (which seems to be the case), instead of O.save/O.load, I suggest e.g. to just export particle positions and import them in the new simulation.

O.save / O.load for "new" simulation is not recommended, for the reasons you are encountering

> The simulation does not run the command every 10 000 iteration, as defined.
> Instead is seems to run ParticleSwelling() for realPeriod=1,

The two sentences are not necessarily contradicting. It is possible it is running every 10000 iteractions AND every 1 real second.

If you define iterPeriod=10000, then it is run each 10000 iterations.

In the case of O.load, the saved PyRunners are still active, so maybe this is the confusion?

> realPeriod=1 which is A period I defined for commands checkForce()/Particlswelling() in my old (saved) simulation.

In the code provided, you defined iterPeriod=10000 for both cases...

> One more question, in my saved simulation the timestep is 1e-7. But defined a new timestep in loaded simulation to 1e-5. The new simulation still use the old time step of 1e-5. How can I define a new timestep?

Did you mean to write it like this? New time step defined 1e-5 and new simulation uses 1e-5 (which perfectly makes sense)?

Cheers
Jan

Hi Jan,

>If your new simulation is not too much related to the old one (which seems to be the case), instead of O.save/O.load, I suggest e.g. to just export particle positions and import them in the new simulation.

Is it possible to also save interparticle forces as well and particle wall interactions?

>The two sentences are not necessarily contradicting. It is possible it is running every 10000 iteractions AND every 1 real second.

You are right that they necessarily contradicting. But in my simulation ParticleSwelling() command have been run 8 times within iterationperiod of 4000. The command run each 500 iteration.

>Did you mean to write it like this? New time step defined 1e-5 and new simulation uses 1e-5 (which perfectly makes sense)?

Sorry is was a typo. I meant the time step is still 1e-7, even through I defined as 1e-5s.

Best regards,
Mithu

> Is it possible to also save interparticle forces as well and particle wall interactions?

it is possible to save them, it is easy.
The loading and applying would/could be more difficult :-D
E.g. for simple contact law, the normal forces are just there defined by position.
If your concerns are also shear forces, friction sliding, normal forces defined by some non-linear law, then this approach is probably not suitable and O.save/O.load is much more suitable.

> But in my simulation ParticleSwelling() command have been run 8 times within iterationperiod of 4000. The command run each 500 iteration.

Please try to reduce your codes. At least remove the unnecessary commented lines.
E.g. the line
#o.engines = o.engines+[PyRunner(command='ParticleSwelling()', iterPeriod=100000)]
is really distracting and confusing if one of the main problems is PyRunner with ParticleSwelling and iterPeriod 100000 :-)

Another tip (I don' have time and mood to test myself):
PyRunner(command='checkForce()', realPeriod=1, label="fCheck"),
fCheck.command = 'ParticleSwelling()'
So if fCheck.command is actually set like this, it runs every 1 real second.
As you load your simulation, the original O.engines is applied (see below), so maybe it is this case?

> Sorry is was a typo. I meant the time step is still 1e-7, even through I defined as 1e-5s.

You have O.load after setting O.dt, right? (I am not sure, but I got the feeling from the codes)
If so, O.load overwrites everything it can (O.dt, O.engines, ...) with the saved state (just as what load should do).

If you want to change anything (O.dt, PyRunners), do it AFTER O.load()

Cheers
Jan

Hi Jan,
Thank you for your reply.

>#o.engines = o.engines+[PyRunner(command='ParticleSwelling()', iterPeriod=100000)]
is really distracting and confusing if one of the main problems is PyRunner with ParticleSwelling and iterPeriod 100000 :-)

I combined the particle swelling and the other code together, due to issue with loading and applying interparticle force.
I want to run the ParticleSwelling() command every 10000 iteration. Still the command does not run every 10000 iteration and I don't understand why. I also tried with setting iterperiod=1. Still it chose it own iterperiod and seems to be same, irrelevant from my input as iterperiod.
Here is the code:
def stopUnloading():
    for i in upper_punch:
        body=O.bodies[i]
        pos_up=body.state.pos
    if pos_up[2] > Cyl_height/2:
        for i in upper_punch:
            body= O.bodies[i]
            body.state.vel = (0,0,2)
        #plate.state.vel = (0, 0, 0)
        initial_save.append(utils.aabbExtrema()[1][2]+r1)
        initial_save.append(O.iter)
        for j in upper_punch: O.bodies.erase(j)
        fCheck.command = 'Savecheck()'

def Savecheck():
    if save==1:
        o.save('Tablet_swelling.xml')
        o.pause()
    if save==0:
        o.engines = o.engines+[PyRunner(command='ParticleSwelling()', iterPeriod=10000)]

def ParticleSwelling():
    time_current=(O.iter-initial_save[1])*o.dt
    Liq_pos=0.24*(time_current**0.574)
    Liq_pos=Liq_pos*1e-3 #convert to m
    radius=[]
    radius.append(time_current)
    for b in O.bodies:
        if isinstance(b.shape, Sphere):
            par_pos=(initial_save[0]-b.state.pos[2])
            k=b.id
            r_now=b.shape.radius
            if Liq_pos>=par_pos:
                r_0=r_save[0][k+1]
                if swell_t[0][k]==0:
                    swell_t[0][k]=time_current
                    radius.append(r_save[0][k+1])
                    continue
                time=time_current-swell_t[0][k]
                t = np.linspace(0,time)
                r = odeint(model,r_0,t,args=(P[0],P[1],P[2],r_0,P[3],))
                r_new=float(r[-1])
                radius.append(r_new)
                b.shape.radius = float(r[-1])
                f=float(r[-1])/r_new
                mcurrent=b.state.mass
                mnew=mcurrent*(f*f*f)
                b.state.mass=mnew
                Icurrent=b.state.inertia
                Inew=Icurrent*(f*f*f*f*f)
                b.state.inertia=Inew
            elif Liq_pos<par_pos:
                radius.append(r_save[0][k+1])
    r_save.append(radius)
    size_current=[time_current,utils.aabbDim()[2]]
    size_save.append(size_current)
    if time_current>20:
        o.pause()
        radius_data=pd.DataFrame(r_save)
        path_save='/home/mithushan/Swelling'
        base_filename='PH101_swelling_data.csv'
        radius_data.to_csv(os.path.join(path_save,base_filename))
        size_data=pd.DataFrame(size_save, columns=['time','height'])
        base_filename='PH101_height.csv'
        size_data.to_csv(os.path.join(path_save,base_filename))
>You have O.load after setting O.dt, right? (I am not sure, but I got the feeling from the codes)
If so, O.load overwrites everything it can (O.dt, O.engines, ...) with the saved state (just as what load should do).

Thank you this helped, it changed the o.dt to 1e-5. But simulation becomes unstablet and spheres just flyes away at the new timestep.

Here is overall new code:
from __future__ import print_function
from yade import utils, plot, timing
from yade import pack
import pandas as pd
import numpy as np
from PIL import Image
from yade import pack, export
from scipy.interpolate import interp1d
readParamsFromTable(k1=2e5,kp=2e6)
from yade.params.table import *
from scipy.integrate import odeint
import matplotlib.pyplot as plt
import os
import csv
from matplotlib.pyplot import figure
from scipy.interpolate import interp1d
from pylab import *
from scipy.optimize import curve_fit

o = Omega()
save=0
# Physical parameters
fr = 0.41
rho = 1561
Diameter = 0.0012/2
D=Diameter
r1 = Diameter/2
#r2 = Diameter/2
k1 = 1.26e5
kp = 12*k1
kc = k1 * 0.1
ks = k1 * 0.1
DeltaPMax = Diameter/3.0
Chi1 = 0.34
PhiF1=0.999
#PhiF1 = DeltaPMax*(kp-k1)*(r1+r2)/(kp*2*r1*r2)

if save==1:
    o.dt = 1.0e-7
else:
    o.dt=1.0e-5

particleMass = 4.0/3.0*math.pi*r1*r1*r1*rho
Vi1 = math.sqrt(k1/particleMass)*DeltaPMax*Chi1
m_tot=0.0005
#no_p=m_tot/particleMass
no_p=1000

Tab_rad=0.005
Cyl_height=0.015
cross_area=math.pi*(Tab_rad**2)

Comp_press= 1.2e8
Comp_force=Comp_press*cross_area

##Single particle swelling model
def model(r,t,Q_max,rho_t,rho_w,r_0,D):
    Q=((rho_w*(r**3))/(rho_t*(r_0**3)))-(rho_w/rho_t)+1;
    drdt =((D*rho_t)/(r*rho_w))*((Q_max-Q)/Q);
    return drdt
P=[1.45,rho,1000,396.39e-12]

#*************************************
# Add material
mat1 = o.materials.append(LudingMat(frictionAngle=fr, density=rho, k1=k1, kp=kp, ks=ks, kc=kc, PhiF=PhiF1, G0 = 0.0))

# Spheres for compression and walls
sp=pack.SpherePack()
sp.makeCloud((-5.0*Diameter,-5.0*Diameter,-12*Diameter),(5.0*Diameter,5.0*Diameter,7.0*Diameter), rMean=Diameter/2.0,rRelFuzz=0.18,num=1000)
sp.toSimulation(material=mat1)
walls=o.bodies.append(yade.geom.facetCylinder((0,0,0),radius=Tab_rad,height=Cyl_height,segmentsNumber=20,wallMask=6,material=mat1))

r_save=[]
radius=[]
initial_save=[]
size_save=[]
radius.append(0)
swell_t=np.zeros((1,no_p))
i=0
for b in O.bodies:
    if isinstance(b.shape, Sphere):
        radius.append(b.shape.radius)
r_save.append(radius)

# Add engines
o.engines = [
  ForceResetter(),
  InsertionSortCollider([Bo1_Sphere_Aabb(aabbEnlargeFactor=1.05),
                         Bo1_Wall_Aabb(),
                         Bo1_Facet_Aabb()
                         ]),
  InteractionLoop(
    [Ig2_Sphere_Sphere_ScGeom(interactionDetectionFactor=1.05),
     Ig2_Facet_Sphere_ScGeom(),
     Ig2_Wall_Sphere_ScGeom()],
    [Ip2_LudingMat_LudingMat_LudingPhys()],
    [Law2_ScGeom_LudingPhys_Basic()]
  ),
  NewtonIntegrator(damping=0.1, gravity=[0, 0, -9.81]),
  PyRunner(command='checkForce()', realPeriod=1, label="fCheck"),
  #DeformControl(label="DefControl")
]

def checkForce():
    # at the very start, unbalanced force can be low as there is only few
    # contacts, but it does not mean the packing is stable
    if O.iter < 400000:
        return
    # the rest will be run only if unbalanced is < .1 (stabilized packing)
    timing.reset()
    if unbalancedForce() > 1:
        return
    # add plate at upper box side

    highSphere = 0.0
    for b in O.bodies:
        if highSphere < b.state.pos[2] and isinstance(b.shape, Sphere):
            highSphere = b.state.pos[2]
        else:
            pass

    #o.bodies.append(wall(highSphere, axis=2, sense=-1, material=mat1))
    #global plate
    #plate = o.bodies[-1] # the last particles is the plat
    #plate.state.vel = (0, 0, -2)
    global upper_punch
    upper_punch=O.bodies.append(geom.facetCylinder((0,0,((-Cyl_height/2)+0.0001)+utils.aabbDim()[2]),Tab_rad-.00001,0,segmentsNumber=50,wallMask=1))
    for i in upper_punch:
        body= O.bodies[i]
        body.state.vel = (0,0,-2)
    fCheck.command = 'unloadPlate()'

def unloadPlate():
    force_up=0
    for i in upper_punch:
        body= O.bodies[i]
        force_up=force_up+abs(O.forces.f(body.id)[2])
    # if the force on plate exceeds maximum load, start unloading
    # if abs(O.forces.f(plate.id)[2]) > 5e-2:
    if force_up > Comp_force:
        for i in upper_punch:
            body= O.bodies[i]
            body.state.vel = (0,0,2)
        #plate.state.vel *= -1
        # next time, do not call this function anymore, but the next one
        # (stopUnloading) instead
        fCheck.command = 'stopUnloading()'

def stopUnloading():
    for i in upper_punch:
        body=O.bodies[i]
        pos_up=body.state.pos
    if pos_up[2] > Cyl_height/2:
        for i in upper_punch:
            body= O.bodies[i]
            body.state.vel = (0,0,2)
        #plate.state.vel = (0, 0, 0)
        initial_save.append(utils.aabbExtrema()[1][2]+r1)
        initial_save.append(O.iter)
        for j in upper_punch: O.bodies.erase(j)
        fCheck.command = 'Savecheck()'

def Savecheck():
    if save==1:
        o.save('Tablet_swelling.xml')
        o.pause()
    if save==0:
        o.engines = o.engines+[PyRunner(command='ParticleSwelling()', iterPeriod=10000)]

def ParticleSwelling():
    time_current=(O.iter-initial_save[1])*o.dt
    Liq_pos=0.24*(time_current**0.574)
    Liq_pos=Liq_pos*1e-3 #convert to m
    radius=[]
    radius.append(time_current)
    for b in O.bodies:
        if isinstance(b.shape, Sphere):
            par_pos=(initial_save[0]-b.state.pos[2])
            k=b.id
            r_now=b.shape.radius
            if Liq_pos>=par_pos:
                r_0=r_save[0][k+1]
                if swell_t[0][k]==0:
                    swell_t[0][k]=time_current
                    radius.append(r_save[0][k+1])
                    continue
                time=time_current-swell_t[0][k]
                t = np.linspace(0,time)
                r = odeint(model,r_0,t,args=(P[0],P[1],P[2],r_0,P[3],))
                r_new=float(r[-1])
                radius.append(r_new)
                b.shape.radius = float(r[-1])
                f=float(r[-1])/r_new
                mcurrent=b.state.mass
                mnew=mcurrent*(f*f*f)
                b.state.mass=mnew
                Icurrent=b.state.inertia
                Inew=Icurrent*(f*f*f*f*f)
                b.state.inertia=Inew
            elif Liq_pos<par_pos:
                radius.append(r_save[0][k+1])
    r_save.append(radius)
    size_current=[time_current,utils.aabbDim()[2]]
    size_save.append(size_current)
    if time_current>20:
        o.pause()
        radius_data=pd.DataFrame(r_save)
        path_save='/home/mithushan/Swelling'
        base_filename='PH101_swelling_data.csv'
        radius_data.to_csv(os.path.join(path_save,base_filename))
        size_data=pd.DataFrame(size_save, columns=['time','height'])
        base_filename='PH101_height.csv'
        size_data.to_csv(os.path.join(path_save,base_filename))
if save==1:
    O.run();
#if save==1:
    #o.save('Tablet_swelling.xml')
if save==0:
    o.load('Tablet_swelling.xml')
    o.dt=1e-7
    #utils.loadVars('somethingtest')
    #from yade.params.somethingtest import *
    initial_save.append(O.iter)
    initial_save.append(utils.aabbExtrema()[1][2]+r1)

Sherlock Holmes approach:
If you have suspicion about realPeriod=1, focus on this and what is going on there.
It is used here: PyRunner(command='checkForce()', realPeriod=1, label="fCheck"),
It calls checkForce. In checkforce: fCheck.command = 'unloadPlate()'
in unloadPlate: fCheck.command = 'stopUnloading()'
in stopUnloading: fCheck.command = 'Savecheck()'
in SaveCheck: o.engines = o.engines+[PyRunner(command='ParticleSwelling()', iterPeriod=10000)]
Still done every 1 real second.

So every 1 real second you add a NEW PyRunner with command='ParticleSwelling()' and iterPeriod=10000.
Now, if you add a new PyRunner in the middle of simulation and O.iter > iterPeriod, it is called immediately the next iteration. Just then it would wait 10000 iteration for next call.

So your PyRunner works as expected, calling ParticleSwelling() every 10000 iteration. Just there is many of them..

To fix it, you can do e.g.:
###
def Savecheck():
    ...
    if save==0:
        o.engines = o.engines+[PyRunner(command='ParticleSwelling()', iterPeriod=10000)]
        fCheck.dead = True # (!!!)
###

Cheers
Jan