Abstract:

Most numerical simulations contain several steps, for example
initialization of model state, defining boundary conditions, stepping
through simulation time, statistical analysis on results, and
visualiation.  In this talk I demonstrate how UNIX principles of
process modularity and text streams allow for unparalleled flexibility
for scientific applications.

## brcon2021 - 2021-06-14

       title: Unix principles for science simulations

      author: Anders Damsgaard (adc)

     contact: [email protected]
              gopher://adamsgaard.dk
              https://adamsgaard.dk


## About me
* 34 y/o Dane
* #bitreich-en since 2019-12-16, member since 2020-08-28

present:

* Aarhus University (DK)

previous:

* Stanford University (US)
* Princeton University (US)
* National Oceanic and Atmospheric Administration (NOAA, US)
* Danish Environmental Protection Agency (DK)
* Scripps Institution of Oceanography, University of California, San Diego (US)

#pause
academic interests:

* ice sheets, glaciers, and climate
* earthquake and landslide physics
* modeling of fluid flows and granular materials


## Granular materials are everywhere
   "The handling and transport of granular materials consumes
    roughly 10% of the world energy consumption."

   - J. Duran, "Sands, Powders, and Grains", doi://10.1007/978-1-4612-0499-2

#pause
Granular materials can be simulated with the Discrete Element Method (DEM).

#pause
Basics of the DEM:
* represent each grain/particle individually (Lagrangian framework)
#pause
* resolve interactions between grains and external forcings
#pause
* update grain kinematics based on interactions and go to next time step

#pause
DEM usage examples:

* sediment mechanics
* beach stability during tsunamis
* landslides
* earthquake initiation
* planetary formation
* meteor impacts
* rock failure
* iceberg or sea ice drift
* planetary rover design
* industrial processes
* pharmaceutical manufacturing

## Configuring a DEM simulation

1. Specify grain positions and properties
#pause

2. Specify boundary conditions
#pause

3. Run simulation through model time in small increments
#pause
   a. Identify contacts between grains
#pause
   b. Resolve contact mechanics, find force on grain for each grain-grain contact
#pause
   c. Add any external forces
#pause
   d. Update acceleration based on sum of forces for each grain
#pause
   e. Use acceleration to determine new velocity and position (linear and rotational)
#pause
   f. Advance model time (t += dt), optionally save state to disk
#pause
   g. End simulation if t >= t_end
#pause

4. Perform statistics on output and/or visualize grains through time

## Configuring a simulation
Many implementations of DEM software, typically with a custom markup for simulation setup.

#pause
Example from DEM software https://github.com/lammmps/lammmps (examples/ASPHERE/box/in.box)

   atom_style      sphere
   atom_modify map array
   dimension   3
   boundary    p p f
   newton      off
   comm_modify vel yes
   units       si

   read_restart    restart30000000.bishear

   region      reg block 0 0.06 0 0.06 -1 0.046 units box
   neighbor    0.001 bin
   neigh_modify    delay 0

   # define pair style
   pair_style  gran/hertz/history 37e9 45e9 1000 1000 0.5  1
   pair_coeff  * *

   # add walls
   fix             tbwall nve_group wall/gran hertz/history 37e9 45e9 1000 1000 0.0 1 ...

   ...

## Problems

- Each program redesigns syntax, custom file format (or Python module)
- Doesn't integrate well with other Unix tools
- No easy user extensibility

## Solution?
Unix-like design:
* Split funcionality into smaller programs
#pause
* Do one thing, do it well
#pause
* Communicate in plaintext in common format
#pause
* Use pipes (or files) as inter-process communication
#pause
* Document use with man pages
#pause
* Modular design with intercompatibility in mind

## granular
   DEM software covered in this talk:

       git://src.adamsgaard.dk/granular
           - granular dynamics simulation
           - soft-body discrete element method
           - C99, make, shell script, POSIX compatible

#pause
   intended as successor to:

       git://src.adamsgaard.dk/sphere
           - GPU-based 3D discrete element method with optional fluid coupling

       git://src.adamsgaard.dk/Granular.jl
           - Julia package for granular dynamics simulation

       git://src.adamsgaard.dk/simple_DEM
           - Basic DEM for educational purposes

## granular(5)
* Adheres to Unix principles
* Design heavily inspired by sfeed (git://git.codemadness.org/sfeed)  #sfeed-fantasies


## One standard format for inter-process communication
GRANULAR(5)                   File Formats Manual                  GRANULAR(5)

NAME
    granular – format for granular dynamics simulations

DESCRIPTION
    granular(1) writes individual grain data in a TAB-separated format to
    stdout.

TAB-SEPARATED FORMAT FIELDS
    Each line represents one grain in a TSV-like format.

    The order, content, and units of the fields are:

    1. diameter [m]
    2. position, x [m]
    3. position, y [m]
    4. position, z [m]
    5. velocity, x [m]
    6. velocity, y [m]
    7. velocity, z [m]
    8. acceleration, x [m]
    9. acceleration, y [m]
    10. acceleration, z [m]
    ...
    54. thermal energy [J]
    55. color [-]



## A simulation with git://src.adamsgaard.dk/granular
Consistent behavior across programs:
   - Communication in granular(5) format
   - Reuse code between programs
   - Optional argument naming always consistent
#pause


1. Specify grain positions and properties
   - granulargrain(1)
   - granularpacking(1)

#pause
2. Specify boundary conditions
   - granular(1)

#pause
3. Run simulation through model time in small increments
   a. Identify contacts between grains
   b. Resolve contact mechanics, find force on grain for each grain-grain contact
   c. Add any external forces
   d. Update acceleration based on sum of forces for each grain
   e. Use acceleration to determine new velocity and position (linear and rotational)
   f. Advance model time (t += dt), optionally save state to disk
   g. End simulation if t >= t_end
   - granular(1)

#pause
4. Perform statistics on output and/or visualize grains through time
   - granularenergy(1)
   - granular2img(1)
   - granular2vtu(1)

## granular.c: Core functionality
In main():

   66     sim_read_grains(&sim, stdin);
   67     sim_add_acceleration_scalar(&sim, grav, 1);
   68     sim_set_timestep(&sim);
   69     if (sim.t < sim.t_end)
   70         sim_run_time_loop(&sim);
   71     sim_print_grains(stdout, &sim);
   72     sim_free(&sim);

#pause
Need something else?  Simple to rewrite or incorporate in another program.

   66     sim_read_grains(&sim, stdin);
   67     sim_add_acceleration_scalar(&sim, grav, 1);
   68     sim_set_timestep(&sim);
   69     if (sim.t < sim.t_end) {
   70         sim_step_time(&sim);
   71         /* add some extra physics here */
   72     }
   73     sim_print_grains(stdout, &sim);
   74     sim_free(&sim);

## Basic examples
Generate a packing, simulate behavior over time, show final state visually:
   $ granularpacking | granular | granular2img | zathura -

#pause
Generate two grains and collide:
   $ (granulargrain -R -u 0.1; granulargrain -f -x 1.2) | granular -e 4.0 -I 0.1 collision

#pause
Determine mean grain size:
   $ granularpacking | cut -f1 | mean

#pause
Plot frequency of grain sizes in generated power-law distribution
   $ granularpacking -P | cut -f1 | histpdf > grain-size-distribution.pdf

#pause
mean(1) and histpdf(1) from git://src.adamsgaard.dk/mathtools

## More thorough example: two-grain-collision.sh
   #!/bin/sh
   id=two-grain-collision
   rm -f ${id}.grains.*.{tsv,png} "${id}.mp4"
#pause

   (granulargrain -R -u 0.1; granulargrain -f -x 1.2) | granular -e 4.0 -I 0.1 "$id"
#pause

   for f in ${id}.grains.*.tsv; do
       granular2img -f '$14' -l 'force_x [N]' -t png < "$f" > "${f%.tsv}.png"
   done
#pause

   ffmpeg -y -framerate 5 -i ${id}.grains.%05d.png \
       -c:v libx264 -r 30 -pix_fmt yuv420p "${id}.mp4"
#pause

   > "${id}.energy.tsv"
   for f in ${id}.grains.*.tsv; do
       granularenergy < "$f" >> "${id}.energy.tsv"
   done
#pause

   gnuplot -e "set term png;\
               set xlabel 'time step';\
               set ylabel 'Energy [J]';\
               plot '-' u 0:1 w lp t 'Total energy'" \
               < "${id}.energy.tsv" > "${id}.energy.png"

## TODO
* Learn additional useful ways of Unix program design
* Add tool to generate packings from sattelite images
* Implement more efficient contact search
* Refine contact mechanics
* Add MPI parallelization

It's just for fun! (for now)

## Thanks for listening!

   Thanks to everyone sharing well-designed software!

#pause
   Thanks to 20h for arranging BRCON2021!