How to restart a simulation from a saved point with restoring geometry movement?

Asked by gaoxuesong on 2020-02-14

Hi. I want to try some parameters that come into play at a late time point, like 1.0s. So i save the simulation at 0.9s and then restart it by O.load() method, followed by trying new parameters.

I have a roller constructed by facets doing translation and rotation controlled by RotationEngine and TranslationEngine, respectively. The parameter of "ids" in those engines is used to bond the engine to the roller (facets).

I use the pickle module to save the roller factets ids into a file during the first run and load these ids again in the restart run. So in the new run, the engines enable to find those roller facets to control. I have checked those ids of facets and they are right.

The problem is that the roller can not rotate or translate in the restart run.

My MWE is as the following, ( the scenario)
https://drive.google.com/file/d/1kQN9ZHgJzktLZQ57nWBMbXDa2joxDzUz/view?usp=sharing )
the first run saves every 0.05s and the restart run loads data at time 0f 0.3s. Because the rotation is active till 0.4s, the roller in the restart run should rotate at the beginning. It's stationary yet.

########## the first run to save at a time point ######## ##

#!/usr/bin/python
#-*- coding: utf-8 -*-
from yade import pack,geom,utils
import numpy as np
from yade import qt

## material
matSph = CohFrictMat(density=8e11, young=193e9, poisson=0.3)
Mat = O.materials.append(matSph)

## roller
radiusSweeper = 0.3e-3
lengthSweeper = 1.0e-3
numSweeperParts = 20
Sweeper=[]
for i in np.linspace(0, 2*pi, num=numSweeperParts, endpoint=True):
 Sweeper.append(Vector3(radiusSweeper*cos(-i), 0.0, radiusSweeper*sin(-i)+1.3e-3))
SweeperP=[Sweeper, [p+Vector3(0.0,lengthSweeper,0.0) for p in Sweeper]]
SweeperPoly = pack.sweptPolylines2gtsSurface(SweeperP, threshold=1e-7)
SweeperIDs = O.bodies.append(pack.gtsSurface2Facets(SweeperPoly, wire=False, material=Mat))
## save roller facets' id for restart ###
import pickle
fd = open('rollerids.txt', 'wb')
pickle.dump(SweeperIDs, fd)
fd.close()

## platform
O.bodies.append(geom.facetBox((2.0e-3,0.5e-3,0.8e-3),(2.0e-3,0.5e-3,0.2e-3), wallMask=63, wire=False, material=Mat))

### process parameters
velsw = 2.5e-3 ## translation vel
angvelsw = -50.0 ## rotation vel
t1 = 0.4 ## 0~0.2s, rotation + translation
t2 = 0.8 ## 0.2~0.4, translation
def change_motion():
    if O.time < t1:
        rotaEngineSw.dead = False
        transEngineSw.dead = False
        rotaEngineSw.zeroPoint = (O.time*velsw, 0, 1.30e-3)
    elif O.time < t2:
        rotaEngineSw.angularVelocity = 0.0
        transEngineSw.velocity = velsw
    else:
        transEngineSw.velocity = 0.0

## save per 0.05s
svdt = 0.05
saveTime = 0
def svfun():
    global svdt, saveTime
    if O.time> saveTime:
        O.save('time-'+str(O.time)[:4])
        saveTime += svdt
    if O.time>1.2:
        O.pause()

O.engines=[
 ForceResetter(),
 InsertionSortCollider([Bo1_Sphere_Aabb(),Bo1_Facet_Aabb()]),
 InteractionLoop(
         [Ig2_Sphere_Sphere_ScGeom6D(),Ig2_Facet_Sphere_ScGeom6D()],
         [Ip2_CohFrictMat_CohFrictMat_CohFrictPhys()],
          [Law2_ScGeom6D_CohFrictPhys_CohesionMoment()]
    ),
    PyRunner(iterPeriod=100,command='change_motion()'),
    RotationEngine(dead=True, label='rotaEngineSw', rotateAroundZero=True, zeroPoint=(0,0,1.30e-3), angularVelocity=angvelsw, rotationAxis=[0,1,0], ids=SweeperIDs)+ \
    TranslationEngine(dead=True, label='transEngineSw', translationAxis=[1,0,0], velocity=velsw, ids=SweeperIDs),
    NewtonIntegrator(damping=0.75, exactAsphericalRot=True, gravity=(0,0,-9.81)),
    PyRunner(iterPeriod=100,command='svfun()'),
]

## particles
sp=pack.SpherePack()
sp.makeCloud((0.5e-3,0,1.0e-3),(2.0e-3,1e-3,1.4e-3),rMean=40e-6,rRelFuzz=0.5)
sp.toSimulation(material=Mat)

O.dt = 0.85*utils.PWaveTimeStep()

#O.run()

######## the restart run ########
##########################

#!/usr/bin/python
#-*- coding: utf-8 -*-
from yade import pack,geom,utils
import numpy as np
from yade import qt

## load a saved data
O.load('time-0.30')

import pickle
fd = open('rollerids.txt', 'rb')
SweeperIDs = pickle.load(fd)
fd.close()

### process parameters
velsw = 2.5e-3 ## translation vel
angvelsw = -50.0 ## rotation vel
t1 = 0.4 ## 0~0.2s, rotation + translation
t2 = 0.8 ## 0.2~0.4, translation
def change_motion():
    if O.time < t1:
        rotaEngineSw.dead = False
        transEngineSw.dead = False
        rotaEngineSw.zeroPoint = (O.time*velsw, 0, 1.30e-3)
    elif O.time < t2:
        rotaEngineSw.angularVelocity = 0.0
        transEngineSw.velocity = velsw
    else:
        transEngineSw.velocity = 0.0

## save per 0.05s
svdt = 0.05
saveTime = 0
def svfun():
    global svdt, saveTime
    if O.time>saveTime:
        O.save('time-r-'+str(O.time)[:4])
        saveTime += svdt
    if O.time>1.2:
        O.pause()
O.engines=[
 ForceResetter(),
 InsertionSortCollider([Bo1_Sphere_Aabb(),Bo1_Facet_Aabb()]),
 InteractionLoop(
         [Ig2_Sphere_Sphere_ScGeom6D(),Ig2_Facet_Sphere_ScGeom6D()],
         [Ip2_CohFrictMat_CohFrictMat_CohFrictPhys()],
          [Law2_ScGeom6D_CohFrictPhys_CohesionMoment()]
    ),
    PyRunner(iterPeriod=100,command='change_motion()'),
    RotationEngine(dead=True, label='rotaEngineSw', rotateAroundZero=True, zeroPoint=(0,0,1.30e-3), angularVelocity=angvelsw, rotationAxis=[0,1,0], ids=SweeperIDs)+ \
    TranslationEngine(dead=True, label='transEngineSw', translationAxis=[1,0,0], velocity=velsw, ids=SweeperIDs),
    NewtonIntegrator(damping=0.75, exactAsphericalRot=True, gravity=(0,0,-9.81)),
    PyRunner(iterPeriod=100,command='svfun()'),
]

O.dt = 0.85*utils.PWaveTimeStep()

#O.run()

#### code finishing ###

Thanks!

Question information

Language:
English Edit question
Status:
Solved
For:
Yade Edit question
Assignee:
No assignee Edit question
Solved by:
gaoxuesong
Solved:
2020-02-18
Last query:
2020-02-18
Last reply:
2020-02-16
Jan Stránský (honzik) said : #1

Hello,

you O.load a simulation and newly set exactly the same list of engines, why?
Why simply not let the loaded engines?
With original O.engines, the cylinder does rotate and move..

With newly set not, I guess it is a gotcha, but before more detailed investigation, please consider using originally O.loaded engines

cheers
Jan

gaoxuesong (260582472-9) said : #2

Hi Jan,

Thank you very much. Yes, you are right. The simulation will proceed if i use the original O.engines.

I've got some tips:
1. There's no need to save the ids of the facets binded to a movement engine. Even though i can not get the ids in the restart simulation, the engine knows those facets.
2. The functions invoked by PyRunner should be included in the restart python script and they can be modified.
3. The arguments in Rotation and Translation engine should be assigned explicitly when updated.

Thanks.