Monday, December 14, 2009

Visualizing force chains out of ESyS-Particle simulations: fcconv

I report here my initial experience in visualizing force chains out of numerical simulations of sheared granular layers using ESyS-Particle.
In ESyS-Particle 2.0, there is a post-processing tool called fcconv that can help you with doing that. Its logic is similar to dump2vtk's one.

Let's start with the basics.

Where can you find fcconv ?
The program itself is in ESyS-Particle's binary directory. Its source code is in

<ESyS-Particle source tree>/Tools/ForceChains .

What's its command line ?

fcconv -i
< input filename > -o < ouput filename > -vvf -vtk ,

where <
input filename> is the input filename, <output filename> is the output filename, the option -vvf is used if you want an output file in a format readable with Itasca's PFC2D(3D) DEM software while -vtk it's used if you want an output file in VTK-XML format.
You can use either -vvf or -vtk or both.

The output file is an unstructured grid constituted by points and cells.
The points are centers of particles belonging to couples of interacting particles, where the interaction is of a specific type (see below). The cells are type=3 VTK cells, i.e., lines, connecting these couples of particles. Associated to the cells there is a DataArray attribute, as it is called in VTK language, specifically, a scalar. This scalar variable corresponds, in the simulation, to the norm of the interaction force acting on one of the two particles, for each couple of interacting particles.
If you wanted the output files in VTK-XML format, such that you could read them in with Paraview, it would be useful to define the file name suffix as .vtu, to remind Paraview that it is a VTK unstructured grid dataset file in XML format. Otherwise you have to specify to Paraview which type of dataset you are importing.

The input file is the output of a specific FieldSaver used inside your ESyS-Particle simulation script. This output refers to a specific time step of the simulation.
That means that, contrary to dump2vtk which manages multiple snapshots in time, you have to use fcconv inside a for loop to convert multiple snapshots in time into the corresponding VTK unstructured grid datasets.

Each input file is an ASCII file and contains the following information (for a specific time step) :
  • each line corresponds to a couple of particles interacting with each other according to a specific interaction type, whose name is specified as a member of the FieldSaver object; this interaction shall be either RotFriction or NRotFriction, so a pure frictional, repulsive, contact force;
  • each line contains 14 fields, each one separated by a blank space from the following one.
The fields (observables) are the following:
  1. particle #1's X coordinate
  2. particle #1's Y coordinate
  3. particle #1's Z coordinate
  4. particle #1's radius
  5. particle #2's X coordinate
  6. particle #2's Y coordinate
  7. particle #2's Z coordinate
  8. particle #2's radius
  9. contact point between particle #1 and particle #2, X coordinate
  10. contact point between particle #1 and particle #2, Y coordinate
  11. contact point between particle #1 and particle #2, Z coordinate
  12. contact force's X component
  13. contact force's Y component
  14. contact force's Z component
The FieldSaver that produces such files is the following one:

InteractionVectorFieldSaverPrms(
interactionName = "friction",
fieldName = "force",
fileName = "/pool2/grm118/ShearFaultGouge2DEx3_33/ContactForce",
fileFormat = "RAW2",
beginTimeStep = 0,
endTimeStep = nt,
timeStepIncr = 500
)

In this example, friction is the name of an interaction group of type RotFriction.
The fieldName indicates that we want to record the full interaction force for each couple of particles subject to this type of interaction. You can choose also fieldName = normal_force and only the projection of the interaction force along the direction passing through the two particles' centers will be recorded.
fileFormat has to be = RAW2 if you are using ESyS-Particle 2.0. I tried to use
fileFormat = RAW_SERIES, like for other FieldSavers, but it did not work. Maybe it's something related to ESyS-Particle's previous versions.
timeStepIncr indicates the time step in between two successive files. It's the sampling time step.
fileName is the FieldSaver's output filename prefix. The program actually add to that prefix a .n.dat suffix, where n is a increasing counter numbering the files (n = 0,1,2,...).
This means that each output filename does not track the simulation time step it refers to.

The VTK-XML unstructured grid dataset files (one for each recording time step) allows you to visualize the contact force network as a network of solid "beams", each one with the same thickness, but with different color scaled according to the scalar attribute, the norm of the contact force.
This contact force network visualization is not the best one: in the literature, you can find several simulation works showing force chains represented as cylindrical beams of different colors AND thicknesses, both scaled according to the value of the norm of the contact force.
However, it's still a useful pictorial representation of contact forces.

For more information, see threads #91547 and #92702 at ESyS-Particle launchpad Web site.

-------------------------------------------------------------------------------------------------
Disclaimer

This post is based on the author's personal experience in using the ESyS-Particle Discrete Element Method code.
This post should not be considered as an official document about the code itself.
This post does not come with any warranty about the correctness of its content.
If you find any error in the content, please write to the author.
-------------------------------------------------------------------------------------------------

Thursday, November 5, 2009

SSH login welcome message

Here there is a very simple tip for changing the initial message appearing on the shell terminal of a user connecting via SSH to the front-end node of a Linux cluster.

  • Modify with whatever text editor you prefer the /etc/motd file

Tuesday, October 27, 2009

Simple ASCII text files: from MS Windows to Linux

Here there is something so trivial that it'll be considered too much trivial for a real Linux geek.
However, since I'm not a Linux hyper-geek, I like to report the following. I think it'll turn out to be useful for not wasting time around for such a small detail.

State of the problem:
  • I edit a simple ASCII text file under MS Windows, for example a PBS script to run a parallel job on a compute cluster.
  • I move the file to the cluster, which, of course, runs Linux
  • I submit the script, via the qsub command. Bash is the shell interpreter in this case. It gives me the following error: ^m command not found
What the hell is going on ?
You may say: "of course, buddy, you're trying to mess up editing files in one environment and running them in another". Well, when your institution makes you use MS Windows for simple desktop computing, there is nothing you can.

Actually, you can do the following:

when the file is on the Linux machine (the front-end node of the cluster, for example), before running it pass it through the dos2unix command. It's a very nice program for converting a text file formatted under DOS/MAC into a text file format for the Unix environment.

Wednesday, August 19, 2009

ESyS-Particle: when you want to use periodic boundary conditions

In ESyS-Particle, when you wanna use periodic (circular) boundary conditions along one (or more) direction(s), you MUST choose to divide the spatial domain of the model into AT LEAST 2 sub-domains along that/those) direction(s). Meaning, you MUST run the simulation in parallel with AT LEAST 2 sub-domains along that(those) direction(s).

That's a really important feature of ESyS-Particle cited in the following post on ESyS-Particlelaunchpad User's Forum .

I was distracted several times by strange results by the simulations and did not run them properly. I looked for the problem and suddenly realized that I forgot to run it in parallel, at least with 3 threads (i.e., processors), one as master and the other two as slaves (the real computational threads).

Do not loose your time in making stupid errors ...

If you want to run a simulation with circular boundary condition, always check about this detail first.

------------------------------------------------------------------------------
Disclaimer

This post is based on the author's personal experience in using the ESyS-Particle Discrete Element Method code.
This post should not be considered as an official document about the code itself.
This post does not come with any warranty about the correctness of its content.
If you find any error in the content, please write to the author.
------------------------------------------------------------------------------

Monday, August 17, 2009

ESyS-Particle: RotBondPrms interaction group 's parameters

Since version 2.0, for 3D simulations, RotBondPrms interaction group's parameters identifying elastic spring constants have been redefined in order to take into account the effect of the particle size distribution.
The parameters, members of the class RotBondPrms, are the following ones:

  • normalK
  • shearK
  • torsionK
  • bendingK
  • normalBrkForce
  • shearBrkForce
  • torsionBrkForce
  • bendingBrkForce
In input, inside a ESyS-Particle Python script, they must be assigned as elastic constants, i.e. moduli. For example, if the basic units used for mass, length and time are Kg, m and s respectively, then these parameters must be assigned in Pascal.
Then, for example, normalK must be expressed as a Young modulus and the corresponding breaking threshold normalBrkForce as a yield stress.

Again, this is valid only for 3D simulations and for ESyS-Particle version 2.0 on.
The same principle applies for the corresponding members of two other classes implementing other different interaction groups: RotFrictionPrms and NRotElasticWallPrms.

-------------------------------------------------------------------------------------------------
Disclaimer

This post is based on the author's personal experience in using the ESyS-Particle Discrete Element Method code.
This post should not be considered as an official document about the code itself.
This post does not come with any warranty about the correctness of its content.
If you find any error in the content, please write to the author.
-------------------------------------------------------------------------------------------------

Thursday, July 30, 2009

dump2vtk: tips & tricks

dump2vtk is a very useful additional tool in ESyS-Particle.
It can be found in the folder ESyS-Particle root folder>/Tools/dump2vtk (for ESyS-Particle version 2.0 and not compiled with a custom folder for binaries, otherwise it'll be in the directory of all the binaries).

It can be used for creating VTK data files out of the files produced by CheckPointer.
What it does is essentially the following.
  • Take as input the different files produced by the CheckPointer (one file for each recording time step); these files are usually called
<CheckPointer fileNamePrefix> _t=????_ID.txt

where ID = 0 for the format file and = 1,2,3, .... for the different processes the simulations is
divided into, ???? stands for the number of the recording time step (in a series).
  • For each recording time step, a corresponding VTK unstructured grid (.vtu) file is created. In case of parallel simulations, dump2vtk takes care of putting together the different CheckPointer files and respective datasets stored by the different processes.
The basic syntax for dump2vtk is the following:

dump2vtk -i -o _ -t tini numsnap deltat [-list] [-bkrlist] [-t0] [-sz] [-rxb] [-single_tag] [-rot] [-unwrap]

Options and their arguments:

-i : setup the CheckPointer fileNamePrefix, which should be equal to the fileNamePrefix defined in the CheckPointPrms class.

-o: setup the VTK output fileNamePrefix

-t: define the initial recording time step (tini), the total number of snapshots that you want to produce (numsnap) and the gap between two recording time steps, i.e., the interval (in time steps) between two CheckPointer files (deltat), such that dump2vtk knows which CheckPointer file to convert next.

-list: instead of constructing the list of input files starting from the CheckPointer fileNamePrefix and from the arguments of the option -t, it takes as input a list of files

-brklist: not yet analyzed

-t0: not yet analyzed

-sz: take only a slice Z = constant

-rxb: not yet analyzed

-single_tag: not yet analyzed

-rot: flag for indicating that the particles are rotational

-unwrap: not yet analyzed

N.B. #1: the VTK output files produced by dump2vtk are numbered according to an increasing sequence. This sequence does not have any relation to the sequence of time steps the different CheckPointer files refer to. Example: suppose that, during an ESyS-Particle simulation, you recorded CheckPointer files every deltat = 100 time steps. You want to convert into VTK snapshots the CheckPointer files from time tini = 1000 to tfin = 5000. That means that you want to convert a total number of numsnap = 51 CheckPointer datasets into corresponding VTK snapshots. You can do that running the command

dump2vtk -i <CheckPointer filename prefix> -o <VTK unstructure grid filename prefix>_ -t 1000 51 100 [...]

The correspondingly generated VTK file names will be indexed from 0 to 50, being their index not related to the actual simulation time step the snapshot refers to.

N.B. #2: there could be a problem when your ESyS-Particle script is not writing CheckPointer
output files in the same directory where it's located. In order to get that, you have to define the CheckPointer fileNamePrefix with the full path, as in the following example:

ckptr_for_ParaView = CheckPointPrms(
fileNamePrefix = "/pool2/grm118/ShearFaultGouge2DEx3_31/2Dfg_test3_31_ckpt",
beginTimeStep = 0,
endTimeStep = nt,
timeStepIncr = 100
)

In such a case, the format file for each CheckPointer (the one whose name terminates with index 0) will point to the CheckPointer sub-domains files via their full paths. Example: for a simulation domain divided in 2 subdomains and for the CheckPointer at time t = 0, the file /pool2/grm118/ShearFaultGouge2DEx3_31/2Dfg_test3_31_ckpt_t=0_0.txt will point to the files /pool2/grm118/ShearFaultGouge2DEx3_31/2Dfg_test3_31_ckpt_t=0_1.txt and /pool2/grm118/ShearFaultGouge2DEx3_31/2Dfg_test3_31_ckpt_t=0_2.txt.

If you move the the CheckPointer files into a different directory where they were created or onto a different machine and then try to run dump2vtk, it won't find the sub-domain Checkpointer files, unless the same path is recreated, in the new machine case.
Suggestion: if you are using full path CheckPointer fileNamePrefix in your ESyS-Particle script, leave the CheckPointer files in the same directory and machine before running dump2vtk. If you need to move the CheckPointer files on a different machine, write the CheckPointer files in the same directory the ESyS-Particle script runs in.
-------------------------------------------------------------------------------------------------
Disclaimer

This post is based on the author's personal experience in using the ESyS-Particle Discrete Element Method code.
This post should not be considered as an official document about the code itself.
This post does not come with any warranty about the correctness of its content.
If you find any error in the content, please write to the author.
-------------------------------------------------------------------------------------------------

Friday, May 29, 2009

Create animations with ImageMagick

A frequent task in post-processing scientific numerical simulations consists of creating a movie of simulation runs.

If your simulation codes saves on files graphical snapshots (plots, images, etc. ...) regarding your system, you can create a corresponding movie at the end of the simulation run.

A simple way to create a movie is to put together in a sequence the single snapshots.

ImageMagick
is a software suite providing lots of commands for image processing.
Using its command convert you can also stack together in a sequence your simulation snapshots.

Simple example of its use:

convert -adjoin -delay X file_snapshot1 file_snapshot2 ... file_snapshotN file_movie.format

While convert is routinely used for converting an image file from one format into another, it can be used also for stacking together image files into a movie file passing to it the option -adjoin.

-delay X indicates the (real) time delay (X) between one snapshot and the following one. The delay is expressed in a default unit: one 100th of second. This default value can be changed.

After delay put the list of image files being your snapshots. If the files are numbered and ordered incrementally, you can just use wildcards (''*'' character in Linux/Unix) for indicating the changing index. Finally you have to indicate the name of the movie file followed by the format extension. ImageMagick will use the right encoder according to the format extension it finds in the movie filename.