CGHandler: Live coarse-graining during a discrete particle simulation

The old fstatistics tool was based on the templated StatisticsVector class. It added functionality to DPMBase that allowed the evaluation of continuum fields, either during a simulation or as a post-processing step. For example, to add global statistical output to your simulation, you had to derive your Driver class from both Mercury3D and StatisticsVector<O>:

class DriverClass : public Mercury3D, public StatisticsVector<O> {
// .. class definition
}
This adds on the hierarchical grid code for 3D problems.
Definition: Mercury3D.h:37
This class is used to extract statistical data from MD simulations.
Definition: StatisticsVector.h:62

This setup was inflexible, however, as it allowed only the evaluation of one type of coarse-graining. This was replaced in MercuryCG by the cgHandler, which can contain several CG objects if desired:

class DriverClass : public Mercury3D {
// .. class definition
}
int main () {
DriverClass dpm; //declare your DPM class
CG<CGCoordinates::O> cg; //declare a cg object
cg.statFile.setSaveCount(50); //set parameters such as the output frequency
dpm.cgHandler.copyAndAddObject(cg); // add the CG object to the cgHandler
dpm.solve(); //run your DPM simulation
}
int main(int argc, char *argv[])
Definition: T_protectiveWall.cpp:215
File statFile
File class to handle the output into a .stat file.
Definition: BaseCG.h:376
Evaluates time-resolved continuum fields and writes the data into a stat file.
Definition: CG.h:76
void setSaveCount(unsigned int saveCount)
Sets File::saveCount_.
Definition: File.cc:273

Output to stat-files

If a cg object is added to the cgHandler, it will produce analytical data while the simulation is running. This data is stored in a file named $name.$n.stat, where $name is the problem name specified by setName, and $n is the index of the cg object in the cgHandler. A basic example is shown in TutorialCG0.cpp. It produces spatially averaged statistics at every 50-th time step, resulting in the following output:

::::::::::::::
TutorialCG0.0.stat
::::::::::::::
CG n 1 1 1 width 1 timeMin 0.1601 min 0 0 0 max 9.92738 9.92738 9.92738 Lucy cutoff 1
time 2:volumeFraction 3:density 4-6:momentum 7-12:momentumFlux 13-21:contactStress 22-24:interactionForceDensity
0.1601 0.535174 1.02211 3.11529e-17 6.62723e-17 -1.38841e-17 0.000348305 9.52371e-06 -6.18274e-06 0.000314944 -1.0499
8e-05 0.000332432 583.523 -93.4013 -191.383 -56.5639 1100.79 -52.4474 -192.268 -48.4432 853.035 0 0 0
0.1651 0.535174 1.02211 1.19044e-17 4.79952e-16 1.43027e-15 75.5124 -29.6753 -15.9208 212.387 10.134 93.6674 1019.24
-2.79241 5.15266 1.58597 1525.6 9.29084 12.7538 -0.135959 956.726 0 0 0
...
const unsigned n
Definition: CG3DPackingUnitTest.cpp:32

The first header line gives the user information about the type of coarse-graining applied; in this case, the coordinate type is O, which denotes globally-averaged statistics.
The second header line denotes the column numbers of the coordinates (time) and the fields (volume fraction, etc) evaluated at the coordinates.
All other lines contain coordinate and field values; e.g. the third line tells you that the volume fraction at time 0 was 0.6.
The simple file format allow easy plotting of the data in gnuplot, e.g.

> set xlabel 'volume fraction'
> set ylabel 'confining pressure'
> p 'TutorialCG0.0.stat' u 2:13
void gnuplot(std::string command)
Plots to a gnuplot window.
Definition: MiscHelpers.cc:38

One can also plot the same data in Matlab using loadStatistics.m provided in the Matlab folder

>> cd ~/MercuryDPM/MercuryBuild/Drivers/Tutorials/
>> addpath("~/MercuryDPM/MercurySource/Matlab")
>> data = loadStatistics("TutorialCG0.0.stat");
>> plot(data.VolumeFraction,data.ContactStressXX);
>> xlabel('volume fraction'), ylabel('confining pressure');

CGCoordinates

The example above produced globally averaged statistics. However, usually coarse-graining is used to extract spatially dependent continuum fields. This can be done by changing the CGCoordinate type:

int main () {
DriverClass dpm;
dpm.setName("Test")
dpm.cgHandler.copyAndAddObject(cgO); // produces globally averaged statistics
cgXYZ.setWidth(1); //set cg width
cgXYZ.setN(20); //set number of mesh points in each spatial direction
dpm.cgHandler.copyAndAddObject(cgXYZ); // produces spatially resolved statistics
dpm.solve();
}
void setN(std::size_t n)
Sets nX_, nY_, nZ_, the number of spatial mesh points in each cartesian direction.
Definition: BaseCG.cc:132
void setWidth(Mdouble width) override
Definition: CG.h:220

The code above produces two sets of output: Test.0.stat contains globally averaged data, Test.1.stat contains spatially resolved data. In the latter case, the number of evaluation points and the coarse-graining width has to be specified – in this case, a Gaussian of width (standard deviation) 1 and a 20x20x20 mesh of evaluation points will be used. Note that the latter statistics will produce huge amounts of data (8000 times more than the globally averaged statistics), even though we have not chosen a particularly high resolution. The data will further be highly fluctuating, as the continuum fields only depend on the small amount of data in the sphere around the evaluation point. Thus, it is advisable to average in space when ever possible, i.e. if the data is homogeneous in space. Often the data is not homogeneous in all spatial directions, but in some. In that case, spatial averaging can be applied in specific directions only. This is done by using partially resolved CGCoordinates: X, Y, Z, XY, XZ, or YZ. For example, CG<CGCoordinates::Z> averages over the yz-plane, but resolves the fields in z-direction. This is shown exemplary in TutorialCG1.1.stat:

::::::::::::::
TutorialCG1.1.stat
::::::::::::::
...

A simple matlab plot reveals that the confining stress increases in the z-direction.

>> cd ~/MercuryDPM/MercuryBuild/Drivers/Tutorials/
>> addpath("~/MercuryDPM/MercurySource/Matlab")
>> data = loadStatistics("TutorialCG1.1.stat");
>> plot(data.z,data.ContactStressZZ);
>> xlabel('z'), ylabel('confining pressure');

  • Coarse-graining function: The default cg function is a cutoff Gaussian, but it can be changed to other functions, such as polynomial radial basis functions (a Heaviside, Linear, or a Lucy polynomial). ...

Temporal averaging

Coarse-graining can also be applied in the time dimension. This can be done in MercuryCG by using a TimeAveragedCG or TimeSmoothedCG object:

int main () {
DriverClass dpm;
dpm.setName("Test")
dpm.cgHandler.copyAndAddObject(cgO); // produces time-resolved statistics
dpm.cgHandler.copyAndAddObject(cgTS); // produces time-smoothed statistics
dpm.cgHandler.copyAndAddObject(cgTA); // produces fully time-averaged statistics
dpm.solve();
}
Evaluates time-averaged continuum fields and writes the data into a stat file.
Definition: TimeAveragedCG.h:59
Evaluates time-smoothed continuum fields and writes the data into a stat file.
Definition: TimeSmoothedCG.h:60

This can best be studied by looking at a time-dependent problem, such as a slowly sheared Lees-Edwards-style simulation: The friction in this system initially builds up until it reaches a maximum, then decreases towards a smaller, critical value.

If your DPM problem is steady, it is not necessary to resolve the continuum fields in time. This can be done in Mercury by using time-averaging cg, e.g.

TimeAveragedCG<CGFunctions::GaussZ> cg;
cg.setNZ(10);
cg.setWidth(0.15);
cg.statFile.setSaveCount(20000);
problem.cgHandler.copyAndAddObject(cg);

To limit the time interval over which the continuum fields are evaluated to e.g. \([10,20]\), use

cg.setTimeMin(10);
cg.setTimeMax(20);

Post-processing

Processing restart files

Documentation of MercuryCG

The list below links to the documentation of the different classes that comprise MercuryCG.

  • class CGHandler:
    Handler that stores all CG objects.
  • classes CG, TimeSmoothedCG, TimeAveragedCG:
    The three available CG classes, used to obtain time-resolved, time smoothed, and time-averaged continuum fields. Contains a vector of CGPoints.
  • templated class CGPoint:
    Derived from a CGFields and a CGFunctions class. Does the main work of evaluating the continuum fields at a given position.
  • namespace CGFields:
    Contains the classes StandardFields etc.
  • namespace CGFunctions:
    Contains the templated classes Gauss, Heaviside, Linear, Lucy (the latter three are derived from Polynomial). Contains the definition of the cg function used. Derived from a CGCoordinates class.
  • namespace CGCoordinates:
    Contains the classes O, X, Y, Z, XY, XZ, YZ, XYZ: Contains the position of the CGPoint.