Yade

O.save/O.load problem

Asked by Roxana Saghafian Larijani

Hi all,
I am trying to continue my simulation in a rotating drum after I stop it. To do this I save my first simulation using O.save at the end of my first script. I also export vtu files for visualisation in Paraview. This works correctly. Then I O.load the saved file (test.bz2) at the beginning of the second script. The simulation goes on and gives me vtu files. However, in the vtu files some points are Nan and I do not get a proper visualisation in Paraview (it does not show the particles in paraview). I was wondering if you could help me with this issue.

I have attached my two scripts here as well :
Script 1:

from yade import ymport
from yade import pack
from yade import utils, plot,wrapper
from yade import export

nCG=1
fr = 0.38
rho = 3700
En = 0.22
Et = 0.22
poi=0.3
yoM=200e6
r = 0.0005

Gamma = 0.072
Theta = 0
vB = 6.03 * 1e-9 * (nCG**3)
Eta= 0.001
CapType="Rabinovich"
Tc=0.00045

##defining material

mat=O.materials.append(
        ViscElCapMat(frictionAngle=fr, density=rho, poisson=poi, young= yoM, Vb=vB, gamma=Gamma, eta=Eta, theta=Theta, Capillar=True, CapillarType=CapType,en=En, et=Et,tc=Tc)
)
mat2=O.materials.append(
        ViscElCapMat(frictionAngle=fr, density=rho, poisson=poi, young= yoM, Vb=vB, gamma=Gamma, eta=0, theta=90, Capillar=True, CapillarType=CapType, en=En, et=Et,tc=Tc)
)

#defining the spheres

sp=pack.SpherePack()
sp.makeCloud((-0.035,-0.035,-0.01),(0.035,0.035,0.01),psdSizes=[0.001,0.002,0.003,0.004], psdCumm=[0.0,0.8,0.95,1.0])
sp.toSimulation(material=mat)
Nprtcl=len(O.bodies)
print(Nprtcl)
Tt= utils.PWaveTimeStep()
O.dt = 0.5*Tt
print(O.dt)

#liquidMigration
VV=0.03
Vmin=0.0
for s in O.bodies:
        if not type(s.shape)==wrapper.Sphere:
                continue
        s.state.Vf=VV * (4/3) * 3.14*(s.shape.radius)**3
        s.state.Vmin=Vmin
##

Drum=geom.facetCylinder(material=mat2,center=(0.0,0.0,0.0),radius=0.05,height=0.03,orientation=Quaternion(Vector3(0,0,1),(pi/2.0)))
walls = O.bodies.append(Drum)

##engine
O.engines = [
        ForceResetter(),
        InsertionSortCollider([Bo1_Sphere_Aabb(), Bo1_Facet_Aabb()]),
        InteractionLoop(
                [Ig2_Sphere_Sphere_ScGeom(),Ig2_Facet_Sphere_ScGeom()],
                [Ip2_ViscElCapMat_ViscElCapMat_ViscElCapPhys()],
                [Law2_ScGeom_ViscElCapPhys_Basic()],

        ),

        NewtonIntegrator(gravity=[0, -9.8, 0]),
        RotationEngine(ids=walls,rotationAxis=[0,0,1],rotateAroundZero=True, zeroPoint=[0.0,0.0,0.0], angularVelocity=2.0),
        LiqControl(label='lqc'),
        VTKRecorder(iterPeriod=1000,recorders=['spheres','facets','colors'],fileName='test-')]

O.run( 500000,True)

O.save('test.bz2')

Script 2:

from yade import ymport
from yade import pack
from yade import utils, plot,wrapper
from yade import export

O.load('test.bz2')
Nprtcl=len(O.bodies)
print(Nprtcl)

O.run(15000,True)

###
BTW, I am using this version yade-2022-04-13.git-5bd3ade.

Thanks in advance for your respopnse!

Question information

Language:
English Edit question
Status:
Solved
For:
Yade Edit question
Assignee:
No assignee Edit question
Solved by:
Karol Brzezinski
Solved:
Last query:
Last reply:
Revision history for this message
Karol Brzezinski (kbrzezinski) said : 		#1

Hi,

Are you sure it is a working example? ViscElCapMat has no attribute 'eta'.

Cheers,
Karol

Revision history for this message
Roxana Saghafian Larijani (r-saghafian) said : 		#2

Sorry, eta is an attribute that I have added to ViscElCapMat in the source code, so please consider the following script without eta:

from yade import ymport
from yade import pack
from yade import utils, plot,wrapper
from yade import export

fr = 0.38
rho = 3700
En = 0.22
Et = 0.22
poi=0.3
yoM=200e6
r = 0.0005

Gamma = 0.072
Theta = 0
vB = 6.03 * 1e-9
CapType="Rabinovich"
Tc=0.00045

##defining material

mat=O.materials.append(
        ViscElCapMat(frictionAngle=fr, density=rho, poisson=poi, young= yoM, Vb=vB, gamma=Gamma, theta=Theta, Capillar=True, CapillarType=CapType,en=En, et=Et,tc=Tc)
)
mat2=O.materials.append(
        ViscElCapMat(frictionAngle=fr, density=rho, poisson=poi, young= yoM, Vb=vB, gamma=Gamma, theta=90, Capillar=True, CapillarType=CapType, en=En, et=Et,tc=Tc)
)

#defining the spheres

sp=pack.SpherePack()
sp.makeCloud((-0.035,-0.035,-0.01),(0.035,0.035,0.01),psdSizes=[0.001,0.002,0.003,0.004], psdCumm=[0.0,0.8,0.95,1.0])
sp.toSimulation(material=mat)
Nprtcl=len(O.bodies)
print(Nprtcl)
Tt= utils.PWaveTimeStep()
O.dt = 0.5*Tt
print(O.dt)

#liquidMigration
VV=0.03
Vmin=0.0
for s in O.bodies:
        if not type(s.shape)==wrapper.Sphere:
                continue
        s.state.Vf=VV * (4/3) * 3.14*(s.shape.radius)**3
        s.state.Vmin=Vmin
##

Drum=geom.facetCylinder(material=mat2,center=(0.0,0.0,0.0),radius=0.05,height=0.03,orientation=Quaternion(Vector3(0,0,1),(pi/2.0)))
walls = O.bodies.append(Drum)

##engine
O.engines = [
        ForceResetter(),
        InsertionSortCollider([Bo1_Sphere_Aabb(), Bo1_Facet_Aabb()]),
        InteractionLoop(
                [Ig2_Sphere_Sphere_ScGeom(),Ig2_Facet_Sphere_ScGeom()],
                [Ip2_ViscElCapMat_ViscElCapMat_ViscElCapPhys()],
                [Law2_ScGeom_ViscElCapPhys_Basic()],

        ),

        NewtonIntegrator(gravity=[0, -9.8, 0]),
        RotationEngine(ids=walls,rotationAxis=[0,0,1],rotateAroundZero=True, zeroPoint=[0.0,0.0,0.0], angularVelocity=2.0),
        LiqControl(label='lqc'),
        VTKRecorder(iterPeriod=1000,recorders=['spheres','facets','colors'],fileName='test-')]

O.run( 500000,True)

O.save('test.bz2')

Regards,
Roxana Saghafian

Revision history for this message
Karol Brzezinski (kbrzezinski) said : 		#3

Hi,

In my opinion there is something wrong with saving the interactions. I don't know how to fix the source. If you inspect interactions right after loading the simulation (before running it) you will notice that many parameters have 'nan' values. After running the simulation you can observe that some spheres just disappear. You can prevent this by resetting the interactions 'O.interactions.clear()', but I guess that is not a solution that you are looking for...

Cheers,
Karol

Revision history for this message
Roxana Saghafian Larijani (r-saghafian) said : 		#4

Hi Karol,

Thanks for your response! I changed the interaction physics to ViscEl instead of ViscElCap (Law2_ScGeom_ViscElPhys_Basic()), and surprisingly the O.save/O.load process is working correctly! So maybe the problem exists only in the model that I am using, and something related to liquid bridges!

Regards,
Roxana Saghafian

Revision history for this message
Karol Brzezinski (kbrzezinski) said : 		#5

Hi Roxana,

I wanted to try to store the interaction state separately and load it in the second script. But there is something wrong with the serialization of interaction physics. For example, if ii is an interaction, calling 'ii.phys.dict()' returns an error:
"No to_python (by-value) converter found for C++ type: yade::CapType"

If you can do your own modifications in the source code, you could try to investigate this line:
https://gitlab.com/yade-dev/trunk/-/blob/master/pkg/dem/ViscoelasticCapillarPM.hpp#L31

Best wishes,
Karol

Revision history for this message
Roxana Saghafian Larijani (r-saghafian) said : 		#6

Hi Karol,

Thanks for your suggestion! In fact, in line 47:
 https://gitlab.com/yade-dev/trunk/-/blob/master/pkg/dem/ViscoelasticCapillarPM.hpp#L47, if CapType gets changed to int, ii.phys.dict() will work fine.
However, in the second simulation with O.load, when I print ii.phys.Dict(), Fn,Fv, normalForce and ShearForce turn out to be Nan!

Regards,
Roxana

Revision history for this message
Best Karol Brzezinski (kbrzezinski) said : 		#7

Hi Roxana,

I don't have your full code, so I cannot run exactly the same example. But maybe you could try my previous idea to store only the parameters you need, and then load it in the second script. I prepared an example for storing Fv and Fn, but you can add the parameters you need.

###### put this at the end of the first script
import pandas as pd

intrState = pd.DataFrame(columns = ['id1','id2','Fn','Fv'])

for ii in O.interactions:
    iiState = pd.DataFrame({'id1':[ii.id1],'id2':[ii.id2],'Fn':[ii.phys.Fn],'Fv':[ii.phys.Fv]})
    intrState = intrState.append(iiState,ignore_index = True)

intrState.to_csv('tmpIntrState.csv')

### put this in the second script before run
import pandas as pd
intrState = pd.read_csv('tmpIntrState.csv')

for i in range(len(intrState)):
    iiSaved = intrState.iloc[i]
    ii = O.interactions[int(iiSaved.id1), int(iiSaved.id2)]
    ii.phys.Fn = iiSaved.Fn
    ii.phys.Fv = iiSaved.Fv

####

Cheers,
Karol

Revision history for this message
Bruno Chareyre (bruno-chareyre) said : 		#8

Hi,
I reported the bug [1]. It sounds like an initialization problem that could be solved by inspecting how the YADE_CLASS_... macro is used.

[1] https://gitlab.com/yade-dev/trunk/-/issues/292

Bruno

Revision history for this message
Roxana Saghafian Larijani (r-saghafian) said : 		#9

Thanks Karol Brzezinski, that solved my question.