yade-users team mailing list archive

[Mithushan Soundaranathan <question699471@xxxxxxxxxxxxxxxxxxxxx>]

New question #699471 on Yade:
https://answers.launchpad.net/yade/+question/699471

Hi,

I am loading (o.load) results from  a saved simulation (o.save). I running a new simulation on loaded data. 
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
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
#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


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


##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='ParticleSwelling()',iterPeriod=10000)
  #DeformControl(label="DefControl")
]


def ParticleSwelling():
    time_current=(o.iter-initial_save[0])*o.dt
    Liq_pos=0.24*(time_current**0.574) #m
    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[1]-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')
    initial_save.append(O.iter)
    initial_save.append(utils.aabbExtrema()[1][2]+r1)
    
    






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