Hi,

I defining ParticleSwelling ()(increase particle radiues over time) by using pyrunner in the engine and define the iterperiod=10 000. The simulation does not run the command every 10 000 iteration, as defined. Instead is seems to run ParticleSwelling() for realPeriod=1, which is A period I defined for commands checkForce()/Particlswelling() in my old (saved) simulation.
So my question, how can I get the command run for each 10 000 simulation.

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?

Best regards,
Mithu

Here is my code:
from __future__ import print_function
from yade import utils, plot, timing
import pandas as pd
import numpy as np
from PIL import Image
from scipy.interpolate import interp1d
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
#elif save==0:
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

Cyl_height=0.015

##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]

#*************************************
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)

r_save=[]
initial_save=[]
size_save=[]
swell_t=np.zeros((1,no_p))
i=0
for b in O.bodies:
if isinstance(b.shape, Sphere):
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='ParticleSwelling()',iterPeriod=10000)
#DeformControl(label="DefControl")
]

def ParticleSwelling():
time_current=(o.iter-initial_save)*o.dt
Liq_pos=0.24*(time_current**0.574) #m
Liq_pos=Liq_pos*1e-3 #convert to m
for b in O.bodies:
if isinstance(b.shape, Sphere):
par_pos=(initial_save-b.state.pos)
k=b.id
if Liq_pos>=par_pos:
r_0=r_save[k+1]
if swell_t[k]==0:
swell_t[k]=time_current
continue
time=time_current-swell_t[k]
t = np.linspace(0,time)
r = odeint(model,r_0,t,args=(P,P,P,r_0,P,))
r_new=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:
size_current=[time_current,utils.aabbDim()]
size_save.append(size_current)
if time_current>20:
o.pause()
path_save='/home/mithushan/Swelling'
base_filename='PH101_swelling_data.csv'
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==0:
initial_save.append(O.iter)
initial_save.append(utils.aabbExtrema()+r1)

Question information

Language:
English Edit question
Status:
For:
Assignee:
No assignee Edit question
Last query:
 Revision history for this message Mithushan Soundaranathan (mithushan93) said on 2021-11-12 #1

Here is code for saved simulation, which loaded into the new simulation:
from __future__ import print_function
from yade import utils, plot, timing
import pandas as pd
import numpy as np
from PIL import Image
from scipy.interpolate import interp1d
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

Cyl_height=0.015

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]

#*************************************
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)

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

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 and isinstance(b.shape, Sphere):
highSphere = b.state.pos
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
for i in upper_punch:
body= O.bodies[i]
body.state.vel = (0,0,-2)

force_up=0
for i in upper_punch:
body= O.bodies[i]
force_up=force_up+abs(O.forces.f(body.id))
# if abs(O.forces.f(plate.id)) > 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

for i in upper_punch:
body=O.bodies[i]
pos_up=body.state.pos
if pos_up > 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()+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()

 Revision history for this message Jan Stránský (honzik) said on 2021-11-12: #2

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

 Revision history for this message Mithushan Soundaranathan (mithushan93) said on 2021-11-12: #3

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

 Revision history for this message Jan Stránský (honzik) said on 2021-11-13: #4

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

it is possible to save them, it is easy.
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

 Revision history for this message Mithushan Soundaranathan (mithushan93) said on 2021-11-13 #5

Hi Jan,

>#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:
for i in upper_punch:
body=O.bodies[i]
pos_up=body.state.pos
if pos_up > 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()+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)*o.dt
Liq_pos=0.24*(time_current**0.574)
Liq_pos=Liq_pos*1e-3 #convert to m
for b in O.bodies:
if isinstance(b.shape, Sphere):
par_pos=(initial_save-b.state.pos)
k=b.id
if Liq_pos>=par_pos:
r_0=r_save[k+1]
if swell_t[k]==0:
swell_t[k]=time_current
continue
time=time_current-swell_t[k]
t = np.linspace(0,time)
r = odeint(model,r_0,t,args=(P,P,P,r_0,P,))
r_new=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:
size_current=[time_current,utils.aabbDim()]
size_save.append(size_current)
if time_current>20:
o.pause()
path_save='/home/mithushan/Swelling'
base_filename='PH101_swelling_data.csv'
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.

 Revision history for this message Mithushan Soundaranathan (mithushan93) said on 2021-11-13: #6

Here is overall new code:
from __future__ import print_function
from yade import utils, plot, timing
import pandas as pd
import numpy as np
from PIL import Image
from scipy.interpolate import interp1d
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

Cyl_height=0.015

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]

#*************************************
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)

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

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 and isinstance(b.shape, Sphere):
highSphere = b.state.pos
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
for i in upper_punch:
body= O.bodies[i]
body.state.vel = (0,0,-2)

force_up=0
for i in upper_punch:
body= O.bodies[i]
force_up=force_up+abs(O.forces.f(body.id))
# if abs(O.forces.f(plate.id)) > 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

for i in upper_punch:
body=O.bodies[i]
pos_up=body.state.pos
if pos_up > 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()+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)*o.dt
Liq_pos=0.24*(time_current**0.574)
Liq_pos=Liq_pos*1e-3 #convert to m
for b in O.bodies:
if isinstance(b.shape, Sphere):
par_pos=(initial_save-b.state.pos)
k=b.id
if Liq_pos>=par_pos:
r_0=r_save[k+1]
if swell_t[k]==0:
swell_t[k]=time_current
continue
time=time_current-swell_t[k]
t = np.linspace(0,time)
r = odeint(model,r_0,t,args=(P,P,P,r_0,P,))
r_new=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:
size_current=[time_current,utils.aabbDim()]
size_save.append(size_current)
if time_current>20:
o.pause()
path_save='/home/mithushan/Swelling'
base_filename='PH101_swelling_data.csv'
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.dt=1e-7
initial_save.append(O.iter)
initial_save.append(utils.aabbExtrema()+r1)

 Revision history for this message Jan Stránský (honzik) said on 2021-11-13: #7

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 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)]