Save and restart simulation for frictional contacts

Asked by Luc Sibille

Dear all,

I would like to save a simulation and be able to restart it latter as described and explained in this question:
https://answers.launchpad.net/esys-particle/+question/132497

However, in this post only bond interactions are considered. I would like to make a similar thing for frictional contacts with rotational particles.
However, I imagine that for frictional contacts, the tangential (or shear) contact forces are computed incrementally. This would mean that to restart a new simulation from a geo file, this file should include the tangential contact forces computed at the last iteration of the previous simulation.
Besides, I understand from this post:
https://answers.launchpad.net/esys-particle/+faq/877
that a geo file does not include such information about the tangential contact force, thus tangential contact forces cannot be saved and reloaded.

Thus, could you please tell me if I am right? Does it mean there is no way to restart a simulation with frictional contacts? Have you suggestions to make a restart possible?

Best regards
Luc

Question information

Language:
English Edit question
Status:
Solved
For:
ESyS-Particle Edit question
Assignee:
No assignee Edit question
Solved by:
Luc Sibille
Solved:
Last query:
Last reply:
Revision history for this message
Luc Sibille (luc-sibille) said :
#1

Dear all,

I try to revive my question with more information.
In the CRotFrictionInteraction class there are two functions:
saveRestartData
and
loadRestartData

Are these functions currently used in ESyS and for what? Do you think there would be a way to use this functions to save and restart a simulation?

Thank you in advance for your replies.

Best regards,
Luc

Revision history for this message
Luc Sibille (luc-sibille) said :
#2

Hi again,

I tried the sim.createRestartCheckPointer, it seems to work fine for RotFriction interactions by saving all the necessary data (among others the friction force). Thus, save a simulation for RotFriction interactions does not seem a problem.
Then I tried the sim.loadCheckPoint with the following file header:

V 2
10000
0.0001
0
LSMGeometry 0
BoundingBox 0 0 0 0 0 0
PeriodicBoundaries 0 0 0
Dimension 3D

But then I do not manage to go further because of inconsistent data... (I imagine ESyS does not appraciate for instance "BoundingBox 0 0 0 0 0 0" and I even don't know if the loadCheckPoint is able to read file with RotFriction interactions)
However in revision 825 it is indicated:
" Restarting
        - added save/load for walls
        - fixed bug in RotParticle save/load
        - tested on triaxial compression"
including modifications in Model/RotFricInteraction.cpp. Thus it seems that someone managed to restart a simulation with RotFriction interactions.
Have you got more information about that? Do you know if loadCheckPoint works?

Regards,
Luc

Revision history for this message
Dion Weatherley (d-weatherley) said :
#3

Hi Luc,

The RestartCheckPointer is the correct one to use if you wish to restart simulations from a previously saved state (regardless of whether you have frictional interactions or not). This feature is still a bit experimental though, so you will need to be patient. There are some hints and tips here:
https://answers.launchpad.net/esys-particle/+question/94168

Often you will need to explicitly set the bounding box (using LsmMpi.setSpatialDomain) in the script used to restart simulations. You may also need to manually change the LSMGeometry version to 1.2.

Cheers,

Dion.

Revision history for this message
Luc Sibille (luc-sibille) said :
#4

Dear Dion,

Thank you for your help and advice.
I finally managed to restart a simulation for rotational particles with frictional contacts. As you said, I needed to change manually the file here named "restart_t=5000_0.txt" written by sim.createRestartCheckPointer. Initially the file contains:

V 2
10000
0.0001
5000
LSMGeometry 0
BoundingBox 0 0 0 0 0 0
PeriodicBoundaries 0 0 0
restart_t=5000_1.txt

I changed it into:

V 2
10000
0.0001
5000
LSMGeometry 1.2
BoundingBox -7 -1 -7 7 11 7
PeriodicBoundaries 0 0 0
Dimension 3D
vrai_restart_t=5000_1.txt

Apparently the version number 1.2 for LSMGeometry is not an obstacle and it seems to work fine even with version number 0.
For the BoundingBox, I used the same definition as the one in the python script defining the simulation (but once again I am not sure that is necessary)
The line about the Dimension was missing, thus I added it.

In the python script I use the command line:
sim.loadCheckPoint("restart_t=5000_0.txt")
as you explained after the definition of the spatial domain, the neighbour search, and the creation of the interaction group between particle;
but before the creation of the interaction groups between particles and walls and for gravity and damping.

I made a very short comparison of the time evolution of the kinetic energy, for a granular assembly falling under gravity into a box, between a non stopping simulation and a simulation saved and reloaded at the half total simulation time. Kinetic energy is exactly the same, so this is encouraging!

Thank you again,
Luc

Revision history for this message
Dion Weatherley (d-weatherley) said :
#5

Hi Luc,

Glad to hear you can now restart simulations. Thanks for sharing details on how you achieved this. I'm sure others will find that helpful in future.

Cheers,

Dion.

Revision history for this message
AndyRathbun (andy-p-rathbun) said :
#6

Hi All,

I am trying to restart my shearbox simulations also and am having some problems.

I create my file with:
 compress_simulation.createRestartCheckPointer(
  CheckPointPrms(
  fileNamePrefix="SAVE",
  beginTimeStep=1000,
  endTimeStep=1000,
  timeStepIncr=1000
  )
 )

And SAVE_t=1000_0.txt contains:
V 2
1000
0.04
1000
LSMGeometry 1.2
BoundingBox 0 0 0 40 30 17.5
PeriodicBoundaries 1 0 0
Dimension 3D
SAVE_t=1000_1.txt SAVE_t=1000_2.txt SAVE_t=1000_3.txt SAVE_t=1000_4.txt SAVE_t=1000_5.txt SAVE_t=1000_6.txt

It seems like the wall loader does not run properly at the restart
compress_simulation.loadCheckPoint("SAVE_t=1000_0.txt")

My complete code is below (or these files can be veiwed at:
http://www3.geosc.psu.edu/~apr137/restart/

I have tried the .loadCheckPoint command in several places, all are marked.

Andy

from esys.lsm import*
from esys.lsm.util import*
from time import*

class Loading(Runnable):
 def __init__(self,sim,dy,t0,fn):
  Runnable.__init__(self)
  self.sim=sim
  self.dy=dy
  self.counter=0
  self.t0=t0
  self.fn=fn
  self.start_time=time()

 def run(self):
  if(self.counter>self.t0):
   self.sim.moveWallBy("lowerWall",Vec3(self.dy,0.0,0.0))
   self.sim.moveWallBy("upperWall",Vec3(-1.0*self.dy,0.0,0.0))
  else:
   dyt=(float(self.counter)/float(self.t0))*self.dy
   self.sim.moveWallBy("lowerWall",Vec3(dyt,0.0,0.0))
   self.sim.moveWallBy("upperWall",Vec3(-1.0*dyt,0.0,0.0))
  self.sim.applyForceToWall("lowerWallInteraction",Vec3(0.0,self.fn,0.0))
  self.sim.applyForceToWall("upperWallInteraction",Vec3(0.0,-1.0*self.fn,0.0))

  self.counter=self.counter+1
  if((self.counter%100)==0):
   runtime=time()-self.start_time
   print "timestep: ",self.counter,self.counter," at ", runtime

def runSimulation():
 cpux=3
 cpuy=2
 cpuz=1
 geometry_name="/home/rathbuna/data/geometry/geometries/poly_8_0.4_200_3000_seis.geo"
 t0=600
 nt=2000
 dt=0.04
 xmax=40.0
 ymax=30.0
 zmax=17.5
 dy=0.0004
 fn=0.0005*xmax*zmax

 compress_simulation=LsmMpi(cpux*cpuy*cpuz,[cpux,cpuy,cpuz])

 compress_simulation.initNeighbourSearch("RotSphere",2.5,0.1)
 compress_simulation.setNumTimeSteps(nt)
 compress_simulation.setTimeStepSize(dt)
 compress_simulation.readGeometry(geometry_name)

 # particle interactions
 fric_coeff=0.6
 stiffness=1.0
 fric_param=RotFrictionPrms("friction",stiffness,fric_coeff,fric_coeff,stiffness)
 compress_simulation.createInteractionGroup(fric_param)

 k=1.0
 kbreak=1000.0
 bond_param=RotBondPrms("bonded",k,k*0.2,k*0.08,k*0.035,kbreak*0.00025,
 kbreak*0.00125,
 kbreak*0.000125,
 kbreak*0.000125,
 0)

 #compress_simulation.loadCheckPoint("SAVE_t=1000_0.txt")
 #errors
#expected Walls , got
#number of bonded interaction groups differ between snapshot and script!
#number of dynamic interaction groups differ between snapshot and script!

 compress_simulation.createInteractionGroup(bond_param)
 compress_simulation.createExclusion("bonded","friction")

 #compress_simulation.loadCheckPoint("SAVE_t=1000_0.txt")
 #bombs the calculations after the old model

 #create upper and lower walls to move and compress the model
 compress_simulation.createWall("lowerWall",Vec3(0.0,0.0,0.0),Vec3(0.0,1.0,0.0))
 compress_simulation.createWall("upperWall",Vec3(0.0,ymax,0.0),Vec3(0.0,-1.0,0.0))

 #side walls
 compress_simulation.createWall("frontWall",Vec3(0.0,0.0,0.0),Vec3(0.0,0.0,1.0))
 compress_simulation.createWall("backWall",Vec3(xmax,ymax,zmax),=Vec3(0.0,0.0,-1.0))

 #compress_simulation.loadCheckPoint("SAVE_t=1000_0.txt")
 #bombs like above

 wp1=NRotElasticWallPrms("frontWallInteraction","frontWall",1.0)
 wp2=NRotElasticWallPrms("backWallInteraction","backWall",1.0)
 wp3=NRotBondedWallPrms("upperWallInteraction","upperWall",1.0,4)
 wp4=NRotBondedWallPrms("lowerWallInteraction","lowerWall",1.0,8)

 #compress_simulation.loadCheckPoint("SAVE_t=1000_0.txt")
 #bombs, like above.

 #Create particle wall interactions
 compress_simulation.createInteractionGroup(wp1)
 compress_simulation.createInteractionGroup(wp2)
 compress_simulation.createInteractionGroup(wp3)
 compress_simulation.createInteractionGroup(wp4)

 #compress_simulation.loadCheckPoint("SAVE_t=1000_0.txt")
 #done makes the model run twice from initial conditions,
 #pins new data at the end of old data

 compress_simulation.createInteractionGroup(
  LinDampingPrms(
  name="damping1",
  viscosity=0.01,
  )
 )
 compress_simulation.createInteractionGroup(
  RotDampingPrms(
  name="damping2",
  viscosity=0.01,
  )
 )

# compress_simulation.createRestartCheckPointer(
# CheckPointPrms(
# fileNamePrefix="SAVE",
# beginTimeStep=1000,
# endTimeStep=1000,
# timeStepIncr=1000
# )
# )

 load=Loading(compress_simulation,dy,t0,fn)
 compress_simulation.addPreTimeStepRunnable(load)
 compress_simulation.run()

if (__name__ == "__main__"):
 runSimulation()

Revision history for this message
Will (will-hancock) said :
#7

Hi Andy,

Here's a few changes which will hopefully get you on your way:

1. When loading from a checkpointer, the line "compress_simulation.readGeometry(geometry_name)" becomes unnecessary as all the particle location information is contained within your checkpointer file. So you can go ahead and comment that out.

2. As we are not using readGeometry() anymore, we have to manually set the domain size and tell it we have periodic boundary conditions. We can do this by adding the following lines of code after we call initNeighbourSearch():
        domain = BoundingBox(Vec3(0.0,0.0,0.0), Vec3(40.0,30.0,17.5))
        compress_simulation.setSpatialDomain (domain,[True,False,False])

3. As you have discovered the location of loadCheckPoint() can be a bit tricky. Steffen might be better placed to comment on this. The general rule i have been using is to place it after you create the particle interactions and exclusions and before anything else.

4. You should also remove the compress_simulation.createWall() function for all 4 walls as they will also be read in automatically.

By the way, nice setup.. i like the tagging. There is one odd thing i notice though is that the front wall doesn't seem to sit flush with the sample. After restart, the gouge material seems to bulge slightly on the front side. Anyway, let us know how you get along with those changes.

Regards,
Will

Revision history for this message
AndyRathbun (andy-p-rathbun) said :
#8

Thanks a lot Will. That got it.

Also thanks for noticing that problem with how I make the shear cell. That was a remnant of how I overlap the ridges, it adds a little bit to the box length but I didn't add that part to the gouge.

Andy

Revision history for this message
Behrooz Ferdowsi (bferdowsi) said :
#9

Hi all,

I would like to remark here that the development version of the code has a bit of change for restarting a simulation that enables arbitrary precision in the restart checkpoints (rev. 1066).

Copy of the changes from a communication with Steffen Abe:

the parameter class for createRestartCheckPointer is now called RestartCheckPointPrms instead of CheckPointPrms and it takes an (optional) 6th argument called "Precision" - which gives the number of output digits wanted.
For example:

----
# Setup restart checkpointers for restarting simulation at a later time
 mySim.createRestartCheckPointer(
  RestartCheckPointPrms(
  fileNamePrefix="snap"
  beginTimeStep = snap_t0,
  endTimeStep = snap_tend,
  timeStepIncr =snap_dt,
                Precision = 30
  )
 )