How to save/restore/continue a simulation

Asked by Feng Chen

Hi, All:

I am going to divide my simulation to two steps:

Step 1, packing of particles under gravity
Step 2, compress the particles using results from Step 1

It seems to me that the most convenient way of doing this is have complete two independent runs, Step 2 will read input from Step 1, e.g. I will have something like:

mySim.readGeometry(step1Output.geo),

in the Step 2 script

however how can we save the particle positions in a .geo text file at end of Step 1?
Is there any function like mySim.SaveGeometry(step1Output.geo) that I can put at the end of Step1.py? PS: I understand the velocity info (translational and angular) will be lost, but it is fine for my problem)

I have read the previous post:
https://answers.launchpad.net/esys-particle/+question/94168
Does the function "sim.createRestartCheckPointer" work properly?

Thanks a lot!

Feng

Question information

Language:
English Edit question
Status:
Solved
For:
ESyS-Particle Edit question
Assignee:
No assignee Edit question
Solved by:
Dion Weatherley
Solved:
Last query:
Last reply:
Revision history for this message
Will (will-hancock) said :
#1

Hi Feng,

I think you are on the right track with this. The easiest way would be to write a GEO file using the end state of the simulation and read it back in using sim.readGeometry(). Provided you settle for long enough and kinetic energy goes to near zero, I view the loss of velocities as a bonus as you are helping the settling come to a complete rest without any vibration through it.

There is currently no SaveGeometry function within ESyS-Particle, it's relatively easy to write a simple geo writer using runnables and getParticleList(), you just have to make sure the formatting is in the correct order. What i would do though if i had my time again.. and i think this would be an excellent contribution to the community is to write a dump2geo.py tool which reads in from the checkpoint files and outputs a GEO file.

I think the second option is actually much less effort with the only trick being to read from the master slave structure of the checkpoint files. Let me know if you need any help with this one.

Regards,
Will

Revision history for this message
Feng Chen (fchen3-gmail) said :
#2

Hi, Will:

Thanks for your prompt reply, I will try both way and see which one is actually easier and more useful.

Feng

Revision history for this message
Feng Chen (fchen3-gmail) said :
#3

Hi, Will:

I am now trying to implement the easy way, using the mySim.getParticleList() at the end of simulation, however I also need to store the bond information between particles, it seems there is no function like mySim.getBondList()??:-), is there a good way to get this list?

Thanks for any suggestions!

Feng

Revision history for this message
Feng Chen (fchen3-gmail) said :
#4

I read the tutorial again, can we achieve this by using a FieldSaver described below?
https://twiki.esscc.uq.edu.au/bin/view/ESSCC/DocumentationAndPresentations#Tables_of_ESyS_Particle_Interact

Feng

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

Hi Feng,

The easy trick here provided your bonds don't break is to copy the list from your old GEO file. I never needed to go beyond that as i ensured that after settling the number of bonds in my simulation remained the same.

You are able to track the number of bonds in you simulation using the following field saver:
sim.createFieldSaver (
   InteractionScalarFieldSaverPrms(
      interactionName="bondGroup",
      fieldName="count",
      fileName="NumBonds.dat",
      fileFormat="SUM",
      beginTimeStep=0,
      endTimeStep=SettlingTimesteps,
      timeStepIncr=1
   )
)

Regards,
Will

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

Hi Feng,

The Checkpoint files contain information on remaining bonds (IDs of particles bonded together and their bond tag) so Will's second suggestion (dump2geo) is probably a better option for you. The Checkpoint file format is not terribly difficult to read from a python script so I'd recommend you try that instead of a combination of sim.getParticleList() and reading a FieldSaver output file.

Cheers,

Dion.

Revision history for this message
Feng Chen (fchen3-gmail) said :
#7

Thanks Dion and Will for detailed suggestions, since I am considering the breakage of the bonds, simply copying the initial bond info from the original geo file might not be an accurate solution for the problem, I am writing a script something like dump2geo to solve this issue.

Best Regards,

Feng

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

Hi Feng,

I just noticed your question about the sim.createRestartCheckPointer(). This feature remains in "beta" but typically it can be convinced to work with a bit of experimentation. Implementing a bullet-proof checkpoint-restart facility for a mature code like ESyS-Particle is very challenging. There are numerous special cases to consider that are not immediately obvious. If you are certain that the loss of information on particle velocities etc. is irrelevant for your purposes, then dump2geo is definitely your best option. If you resort to the RestartCheckPointer, post any errors you encounter and we'll endeavour to help you resolve them.

Cheers,

Dion.

Revision history for this message
Feng Chen (fchen3-gmail) said :
#9

Hi, Dr Dion and all other ESyS users:

I made a very lazy version for the dump2geo with minimum change from dump2vtk, I posted my solution here:

https://sites.google.com/site/fchen3/discreteelementmethod/codesamples

Basically my solution was using dump2vtk's particle list and bond list, then write to a geo file.

I know this is a very crappy way, however ...

To use:

First copy the dump2vtk folder to dump2geo in the ESyS-Particle src/Tools

download my version of frame_vtk.cpp at the website

overwrite the old frame_vtk.cpp

then: in the dump2geo folder, type "make"

if there is no error, there should be a version of dump2vtk executable, change the name to dump2geo, then copy to ESyS-particle's bin folder, then you can use the dump2geo just the same as dump2vtk, you can test the dump2geo using:

dump2geo -i rod_data -o rodvtk_ -t 0 8 1000 -rot

with my sample data posted at the webpage:rod_data.tar.gz (you need to first extract it), note this is only for Rotational, Bonded Particles.

Hope this will help some of the ESyS users. I might be able to write some .py file, but sorry I really really have no time...

Feng

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

Hi Feng,

Many thanks for sharing your dump2geo code. With your permission, I would like to add this to the ESyS-Particle development version and include it in future stable releases. Your contribution will be acknowledged in the source code and I'll add you to the authors/contributors list. Please let me know if it is OK to do this.

Thanks again.

Dion.

Revision history for this message
Feng Chen (fchen3-gmail) said :
#11

Hi, Dr Dion:

I will be very glad to be listed as authors/contributors, and I welcome anyone that find this code useful to change/modify it. This piece of dump2geo still need a lot of cleanup, at least would need a global change in the code from dump2vtk to dump2geo, and handle different cases for Rotational/Non-Rotational particles. One might use "kompare" to easily find what I changed to the dump2vtk's "frame_vtk.cpp"

Another possible problem is that, when we use checkerpointer files, I noticed we might lose one to two significant digits (during the C++ iostream ">>" or "<<" operation), I do not know how this is going to affect simulation results, however I think we might need to control the output of the checkerpoint file format to have a better control..., just a thought:-)

Thanks again for all your suggestions!

Feng

Revision history for this message
Michele Griffa (michele-griffa) said :
#12

Hello Feng and everyone else

Thanks a lot for the fantastic tool dump2geo.

As a user of ESyS-Particle, I agree that this is a very useful pre-/post-processing tool for saving & restarting a simulation.

I just would like to address the following question:

can it be used for a parallel simulation with decomposed domain, thus multiple checkPointer output files (<checkPointer output file name>_1.txt, <checkPointer output file name>_2.txt, ....) ?
dump2vtk is able to collect the information from the multiple checkPointer output files <checkPointer output file name>_?.txt for the different subdomains and condense together the system's snapshot into a signle VTK file.
Does dump2geo do the same ?

Another comment:

I found out that within the thread #98314, titled "write a geometry file", Steffen proposed a Python script called in the same way, dump2geo, for producing a .geo file out of a checkPointer output file. Maybe the two dump2geo do similar things but not completely, I have not checked in details yet.

Thanks a lot again Feng for your contribution.

It will certainly turn out to be very useful for my daily use of ESyS-Particle

Best regards

Michele

Revision history for this message
Feng Chen (fchen3-gmail) said :
#13

Hi, Michele:

This version of dump2geo should work exactly the same way as the dump2vtk with multiple chekerpointer files since I modified the source code based on dump2vtk, and I have tested with multiple checkpointer files (rotational particles) without problems (you can test with the file "rod_data.tar.gz: on the link above). I have not yet tested the python file in thread #98314, although it seems very similar.

Note that at the time I worked on this code was for a conference deadline so I just implemented my particle type "rotational particle", therefore it is not surprising if you have trouble with non-rotational particles, once I get a chance I would do a detailed clean up for this tool.

Let me know if you find any problems or bugs (certainly there is:-)).

Feng

Revision history for this message
Michele Griffa (michele-griffa) said :
#14

Hello Feng and everyone else who's reading this thread.

I had a look at the frame_vtk.cpp code and I noticed that it was automatically considering only LSMGenGeo models without periodic boundary conditions and 3D systems, meaning that, whatever written in the checkPointer output file with index "_0" (the "header" file), 3D and non-periodic boundary conditions were defined for the .geo output file.

So, I change frame_vtk.cpp a little bit in a very C-oriented way in order to keep track of the actual dimensionality of the system and of the flag for the periodic boundary conditions (only along the X-axis for the time being, as implemented till now in ESyS-Particle 2.0, see ) as read from the checkPointer output file.

Here is the result of running the diff command between the original version of frame_vtk.cpp (the one for dump2vtk) and the one for dump2geo:

[....]$ diff -u frame_vtk.cpp ../dump2vtk/frame_vtk.cpp
--- frame_vtk.cpp 2011-01-18 11:44:17.684051184 +0100
+++ ../dump2vtk/frame_vtk.cpp 2009-07-09 04:59:17.000000000 +0200
@@ -53,7 +53,6 @@
   Vec3 circ_shift;
   double rad;
   double mass;
- int id;
   int tag;
   double q1,q2,q3,q4;
   Vec3 angvel;
@@ -147,7 +146,7 @@
   return version;
 }

-vector<string> get_filenames(const string& infilename, int version, double *bdbx=NULL,int *PBCXAxisFlagTemp=NULL, int *DimensionalityTemp=NULL)
+vector<string> get_filenames(const string& infilename, int version)
 {
   cout << "infilename : " << infilename << endl;
   ifstream headerfile(infilename.c_str());
@@ -168,45 +167,19 @@
   // get bounding box
   headerfile >> dummystring;
   headerfile >> xmin >> ymin >> zmin >> xmax >> ymax >> zmax ;
-
- //double boundingBox[6]={xmin,ymin,zmin,xmax,ymax,zmax};
- if (bdbx!=NULL)
- {
- bdbx[0]=xmin,bdbx[1]=ymin,bdbx[2]=zmin;
- bdbx[3]=xmax,bdbx[4]=ymax,bdbx[5]=zmax;
- }
- //std::cerr << bdbx[0] << " " << bdbx[1]<< " " << bdbx[2]<< " " << bdbx[3]<< " " << bdbx[4]<< " " << bdbx[5] << endl;
-
- // Read the flag for periodic boundary conditions along the X-axis (which is the only axis along which periodic BC can be used up to ESyS-Particle 2.0)
- int dummyPBCXAxis;
- headerfile >> dummystring >> dummyPBCXAxis >> dummy >> dummy;
- *PBCXAxisFlagTemp = dummyPBCXAxis;
+
+ // ignore periodic bdry
+ headerfile >> dummystring >> dummy >> dummy >> dummy;

   // v. 1.1 geometry files didn't have dimension info
   if(geo_version>1.15){
     headerfile >> dummystring >> dummystring;
- if (dummystring == "3D")
- {
- *DimensionalityTemp = 3;
- }
- else if (dummystring == "2D")
- {
- *DimensionalityTemp = 2;
- }
- else
- {
- cerr << "Wrong definition of the dimensionality: " << dummystring << endl;
- }
   }

- // Get the file names of the file where the checkPointer recoreded the actual particle system's information
+ // get file names
   copy(istream_iterator<string>(headerfile),istream_iterator<string>(),back_inserter(filenames));

   headerfile.close();
-
- //==================================
-
- //==================================

   cout << "nr. of filenames: " << filenames.size() << endl;
   return filenames;
@@ -357,12 +330,7 @@
 void do_single_frame_vtk_r(const string& infilename,const string& outfilename,int iframe,bool with_list,const string& listfilename,bool remove_xbonds,double bond_remove_dist,bool with_brklist,const string& brklistname)
 {
   int version=get_version(infilename);
- double bdbx[6];
- int *PBCXAxisFlag, *Dimensionality;
- vector<string> filenames=get_filenames(infilename,version,bdbx,PBCXAxisFlag,Dimensionality);
-
- std::cerr << bdbx[0] << " " << bdbx[1]<< " " << bdbx[2]<< " " << bdbx[3]<< " " << bdbx[4]<< " " << bdbx[5] << endl;
-
+ vector<string> filenames=get_filenames(infilename,version);
   map<int,r_part> datamap;
   map<int,int> ng_id_map;
   map<int,float> ng_mass_map;
@@ -478,36 +446,6 @@
     brklistfile.close();
   }

- //std::cerr << "I am being called!!!!!!!!!!!!!!!!!!" << std::endl;
-
- ostringstream geofilename;
- geofilename << outfilename << iframe << ".geo";
- ofstream geofile(geofilename.str().c_str());
- geofile << "LSMGeometry 1.2" << endl;
- geofile << "BoundingBox " << bdbx[0] << " " << bdbx[1] << " " << bdbx[2] << " " << bdbx[3] << " " << bdbx[4] << " " << bdbx[5] << endl;
- geofile << "PeriodicBoundaries " << *PBCXAxisFlag << " 0 0" << endl;
- geofile << "Dimension " << *Dimensionality << "D" << endl;
- geofile << "BeginParticles" << endl;
- geofile << "Simple" << endl;
- geofile << datamap.size() << endl;
- int geopid = 0;
- for(map<int,r_part>::iterator iter=datamap.begin();
- iter!=datamap.end();
- iter++){
- geofile << (iter->second).pos << " "<< (iter->second).rad << " "<< geopid++ << " " << (iter->second).tag << endl;
- }
- geofile << "EndParticles" << endl;
- //write bond info
- geofile << "BeginConnect" << endl;
- geofile << bonds.size() << endl;
- for(vector<bond>::iterator iter=bonds.begin();
- iter!=bonds.end();
- iter++){
- geofile << iter->id1 << " " << iter->id2 << " " << iter->tag << endl;
- }
- geofile << "EndConnect" << endl;
- geofile.close();
- /* The following part, commented, is the original output part of dump2vtk which is not used here anymore
   // open output file
   ostringstream vtkfilename;
   vtkfilename << outfilename << iframe << ".vtu";
@@ -641,7 +579,6 @@

   // close file
   vtkfile.close();
- */
 }

@@ -680,7 +617,7 @@
     datafile >> npart;
     if(version < 2){
       for(int i=0;i<npart;i++){
- datafile >> data.pos >> data.rad >> data.id >> data.tag >> data.mass >> initpos >> oldpos >> data.vel >> data.force
+ datafile >> data.pos >> data.rad >> id >> data.tag >> data.mass >> initpos >> oldpos >> data.vel >> data.force
                 >> data.q1 >> data.q2 >> data.q3 >> data.q4 >> data.angvel;
        if(data.tag==ptag){
          datamap[id]=data;
@@ -688,7 +625,7 @@
       }
     } else {
      for(int i=0;i<npart;i++){
- datafile >> data.pos >> data.rad >> data.id >> data.tag >> data.mass >> initpos >> oldpos >> data.vel >> data.force >> circ_shift
+ datafile >> data.pos >> data.rad >> id >> data.tag >> data.mass >> initpos >> oldpos >> data.vel >> data.force >> circ_shift
                 >> data.q1 >> data.q2 >> data.q3 >> data.q4 >> data.angvel;
        if(unwrap){
         data.pos-=circ_shift;

This new version of frame_vtk.cpp can be downloaded from the following URL:

http://www.calcolodistr.altervista.org/en/work/ESySParticle/ESySParticle_en.html

When I try to build it inside its folder, <ESyS-Particle 2.0 root>/Tools/dump2geo, after having copied within it the Makefile.in and Makefile from the directory <ESyS-Particle 2.0 root>, I have some problems at linking time with the libstdc+ library.
For the errors I get see another thread by myself, #141911.

This problem at build time may be just related to the fact that I did not create a new Makefile limited to dump2geo. I used the general one for all ESyS-Particle 2.0, without re-building all ESyS-Particle.

Best regards

Michele

Revision history for this message
Feng Chen (fchen3-gmail) said :
#15

Hi, Michele and other ESyS developers/users:

I am not sure if we can use subversion to further develop this tool, that would be evidently more efficient and helpful.

Feng

Revision history for this message
Vince Boros (v-boros) said :
#16

Hello Feng and Michele:

I have added dump2geo to the repository, making changes to dump2vtk based on your version, Feng. I had expected to check in the files by the middle of last week, except that Brisbane floods interrupted access to my workplace and the internet. I am currently testing for nonrotational features and two-dimensional geometries. I will check in any refinements as they are needed.

Thanks, Michele, for your suggestions; they are all included. If you need to make manual changes to the build files, look at the new versions of configure.ac and Tools/Makefile.am. You will of course also need Tools/dump2geo/Makefile.am.

Thank you both for your help and advice.

Vince