save and load

Hello,

> I was wondering, what it the best way to save the simulation in the middle and then load it and continue with simulation.

I would say O.save / O.load is pretty standard (depending on the situation and reason as always, of course).

> I am getting error message saying plate is not defined and some other parameter.

"plate" is Python variable.
O.save by default only saves C++ "internals", omitting Python.
If you want also to save Python variables, have a look at [1] and links therein.

Cheers
Jan

[1] https://answers.launchpad.net/yade/+question/698152

Hi Jan,

Thank you for your help.
I tried the suggestion in [1]

Here is how I did it:

For saving:

o.save('Tablet_swelling.xml')
utils.saveVars('somethingtest',plate=plate)

for load:

from yade.params.somethingtest import *
o.load('Tablet_swelling.xml')
utils.loadVars('somethingtest')

I am getting error message saying: ModuleNotFoundError: No module named 'yade.params.somethingtest'; 'yade.params' is not a package.

When I am saving, I am getting following warning message:
/home/mithushan/myYade/trunk/build/install/lib/x86_64-linux-gnu/yade-2021.01a/py/yade/utils.py:67: UserWarning: Overwriting yade.params.somethingtest which already exists.
  if mark in yade.params.__dict__: warnings.warn('Overwriting yade.params.%s which already exists.'%mark)

Best regards,
Mithu
[1] https://answers.launchpad.net/yade/+question/698152

Here is my 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
from yade.params.somethingtest import *

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)
    fCheck.command = 'unloadPlate()'

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

def stopUnloading():
    if plate.state.pos[2] > Cyl_height/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)]
        if save==1:
            utils.saveVars("somethingtest",plate=plate)
            o.save('Tablet_swelling.xml')

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==0:
    o.load('Tablet_swelling.xml')
    utils.loadVars("somethingtest")

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