Coarse graining LAMMPS output

In this example, we demonstrate how to use the CoarseGraining class to perform coarse graining of LAMMPS output: a 3D granular assembly of spherical particles, read from LAMMPS dump files.

The files required for this example can be found in the repository:

• The DEM simulation data files in examples/data_examples/lammps_reading/DEM_data/ (dump_particles/ and dump_forcechain/)

• The configuration file config.ini in examples/LAMMPS_reading/

The same workflow can be run from the compute_CG.ipynb notebook in examples/LAMMPS_reading/. This example focuses on reading LAMMPS data; for the visualisation of the coarse-grained fields, see the MFiX examples.

Import packages

from pysammos.utils.config_loader import load_config
from pysammos.coarse_graining import CoarseGraining

Hello from pysammos
Loading macroscopic_fields package...
Loading data_read package...
Loading particle_phase package...
Loading spatial_weights package...
Loading neighbour_search package...
Loading grid_generation package...
Loading data_handle package...
Loading data_write package...

from pathlib import Path
import os

# The paths inside config.ini (DEM data, output folder) are relative to the example folder,
# so we work from examples/LAMMPS_reading/ and read its config.ini.
if "EXAMPLE_DIR" not in globals():          # resolve once, so the cell can be re-run
    EXAMPLE_DIR = (Path.cwd() / "../../examples/LAMMPS_reading").resolve()
os.chdir(EXAMPLE_DIR)
CONFIG_FILE = EXAMPLE_DIR / "config.ini"
print("Working directory:", Path.cwd())

Working directory: /exports/csce/datastore/geos/users/s1857688/Coarse_Graining/pysammos/examples/LAMMPS_reading

The configuration file

The configuration file config.ini contains the parameters for the coarse-graining process. You can modify the parameters in this file to suit your needs. For a description of every section, see the bedload transport example. The file used here is:

[paths]
particles_path = ../data_examples/lammps_reading/DEM_data/dump_particles/particles.
contacts_path = ../data_examples/lammps_reading/DEM_data/dump_forcechain/forcechain.
output_path = ./PysammosCG/

[timesteps]
t0 = 11991
tf = 11991
dt = 1
t_pad = 9

[smoothing_function]
type = Lucy

[flags]
partialignore = True

[key_mapping]
Global_ID = id
Particle_Position = [x, y, z]
Particle_Velocity = [vx, vy, vz]
Particle_Diameter = diameter
Particle_Density = v_pdens
Particle_Volume = v_pvol
Particle_Mass = mass
Particle_Radius = radius
Coordination_Number = c_Z
Particle_i_ID = c_pl[1]
Particle_j_ID = c_pl[2]
Force_ij = [c_pr[1], c_pr[2], c_pr[3]]
Contact_ij = None

[DEM_software]
name = lammps

[grid_info]
grid_dimension = 3
grid_axes = xyz
automatic_grid = True
x_min = 0.005
x_max = 0.10
y_min = 0.005
y_max = 0.44
z_min = 0.005
z_max = 0.05
x_transect = None
y_transect = None
z_transect = None
x_axis_periodic = True
y_axis_periodic = False
z_axis_periodic = True

[fields_to_export]
volume_fraction = True
density_particle = False
density_mixture = True
momentum_density = False
velocity = True
velocity_gradient = False
kinetic_tensor = False
contact_tensor = False
total_stress_tensor = True
pressure = True
granular_temperature = True
granular_temperature_slices = True
fabric_tensor = False
inertial_number = True
coordination_number = False
d43 = False
d32 = False
frictional_coefficient = True
shear_rate_tensor = False

[output_options]
vkthdf_output = True
h5_output = True

Mapping the LAMMPS dump columns

With DEM_software = lammps, the key_mapping section gives the column names of the dump files. The headers of the two dump files show the available columns:

particles_dump = sorted(Path("../data_examples/lammps_reading/DEM_data/dump_particles").glob("particles.*"))[0]
contacts_dump = sorted(Path("../data_examples/lammps_reading/DEM_data/dump_forcechain").glob("forcechain.*"))[0]
for dump in (particles_dump, contacts_dump):
    print(f"--- {dump.name}")
    with open(dump) as f:
        print("".join(next(f) for _ in range(9)))

--- particles.000011991.dump
ITEM: TIMESTEP
11991
ITEM: NUMBER OF ATOMS
35648
ITEM: BOX BOUNDS pp pp fm
-1.0000000000000001e-01 1.0000000000000001e-01
-1.0000000000000001e-01 1.0000000000000001e-01
0.0000000000000000e+00 5.0000000000000000e-01
ITEM: ATOMS id type diameter radius mass v_pdens x y z xu yu zu vx vy vz v_pvol c_Z

--- forcechain.000011991.dump
ITEM: TIMESTEP
11991
ITEM: NUMBER OF ENTRIES
105362
ITEM: BOX BOUNDS pp pp fm
-1.0000000000000001e-01 1.0000000000000001e-01
-1.0000000000000001e-01 1.0000000000000001e-01
0.0000000000000000e+00 5.0000000000000000e-01
ITEM: ENTRIES index c_pl[1] c_pl[2] c_pr[1] c_pr[2] c_pr[3] c_pr[4] c_pr[5] c_pr[6] c_pr[7]

• Vector quantities are given as a list of columns, e.g. Particle_Position = [x, y, z] and Particle_Velocity = [vx, vy, vz].

• Per-atom variables and computes keep their LAMMPS names, e.g. Particle_Density = v_pdens, Particle_Volume = v_pvol, Coordination_Number = c_Z.

• In the contact (local) dump, Particle_i_ID = c_pl[1] and Particle_j_ID = c_pl[2] are the IDs of the two particles in contact, and Force_ij = [c_pr[1], c_pr[2], c_pr[3]] are the force components.

• Contact_ij = None: the contact points are not in the dump, so Pysammos calculates them from the particle positions and diameters.

In grid_info, the three *_axis_periodic flags should match the boundary flags of the simulation box, given in the dump header (ITEM: BOX BOUNDS ..., where p means periodic). They are used to correct the branch vectors of contacts that cross a periodic boundary. With automatic_grid = True, the grid range is set from the particle bounds and the x/y/z_min/max values are not used.

Initialise the CoarseGraining class

load_config returns the arguments of CoarseGraining as a dictionary, so they can be passed directly with **. You can also give each argument explicitly, as shown in the bedload transport example.

cfg = load_config(CONFIG_FILE)
CG = CoarseGraining(**cfg)

Output path created: ./PysammosCG/

CG.DEM_software, CG.time_steps

('lammps', array([11991]))

Run the coarse graining

The whole workflow can be run with a single command:

CG.run()

Here we run it step by step instead, to show what each step does.