Skip to content

Tutorial 1. Granular column collapse

In this tutorial, we'll simulate a very simple basic granular column collapse with the modified Cam-Clay. The simulation involves two parts: (1) gravity pack and (2) collapse.

Pre-requisites:

Before starting, please make sure that you have completed the following items:

  1. Installed the HydraxMPM.
  2. Comfortable with JAX - basics, and JAX - The Sharp Bits
  3. Read through an overview of the code structure
  4. Have Paraview

It is recommended that you follow along with the tutorial. Stuck or in a rush? then see code available in GitHub under /tutorials/t1_granular_column.py.

Learning objectives

By the end of this tutorial you will be able to

  • Understand the overall code style of a HydraxMPM script

Step 1: create the project file and folders

Create a project directory containing the driver script and two sub-folders with as follows:

# driver script
/tutorial/t1_granular_column.py  

# output of gravity pack
/tutorial/output/t1_pack/ 

# output of collapse
/tutorial/output/t2_collapse
The output is stored as vtk files in the /tutorial/output/t1_pack/ and /tutorial/output/t2_collapse/ folders, respectively.

Step 2: import modules

Import HydraxMPM and supporting the JAX dependencies. 

import jax
import jax.numpy as jnp
import hydraxmpm as hdx

HydraxMPM is top-level package. It contains all the bits and pieces under the same namespace, e.g., hdx.Particles or hdx.Gravity.

Step 3: create points representing initial granular column

Create rectangular column a rectangular column of material points. Particles are spaced evenly given a cell size and the number of particles per cell (in one direction). We pad so that material points do not touch the boundary. 

column_width = 0.3  # [m]
column_height = 0.4  # [m]
ppc = 2


cell_size = 0.025  # [m]
sep = cell_size / ppc
padding = 4.5 * sep

# create padded meshgrid of positions 
x = jnp.arange(0, column_width + sep, sep) + padding - sep 
y = jnp.arange(0, column_height + sep, sep) + padding - sep
xv, yv = jnp.meshgrid(x, y)

position_stack = jnp.array(list(zip(xv.flatten(), yv.flatten())))
Tip

There are several ways of initializing material points. See the how-to initialize material points .

Tip

Matplotlib or Pyvista may be used visualize initial positions. To visualize the particles copy and paste the following (note plt.show works only when GUI is active). 

    import matplotlib.pyplot as plt

    plt.scatter(*position_stack.T ) # unpack x and y and plot
    plt.show() # or plt.savefig("packing.png")

Step 4: create the simulation config

This is a juicy sandwich of all common general simulation parameters.

config = hdx.MPMConfig(
    origin=jnp.array([0.0, 0.0]),
    end=jnp.array([column_width +  padding+2*sep, 0.5]),
    project="t1_pack",
    ppc=ppc,
    cell_size=cell_size,
    num_points=len(position_stack),
    shapefunction="cubic",
    num_steps=60000,
    store_every=500,
    default_gpu_id=0,
    dt=3 * 10**-5,  # time step[s]
    file=__file__,  # store current location of file for output
)
Ok... it seems like a lot, but actually pretty straight forward:

  • hdx.MPMConfig is the config dataclass used for all MPM simulations.
  • origin contains start x and y coordinates of the domain boundary
  • end contains end x and y coordinates of the domain
  • project is the name of the project - mainly used for output purposes
  • ppc is the number of particles per cell (in one direction). this is only used in hdx.discretize function (more on that later)
  • cell_size is background grid cell size
  • num_points is number of material points
  • shapefunction the interpolation type from particle to grid. Two shapefunctions are supported, either linear or cubic. Cubic is better in most cases.
  • num_steps the total iteration count
  • store_every output every nth step
  • default_gpu_id if you are working on a shared GPU workstation, this is the parameter you change to avoid making the other person(s) angry! Run nvidia-smi, in your terminal to find the id of an empty GPU.
  • dt is a constant time step
  • file=__file__ this records the path of your driver script, which is important to save relative output in the correct folder.
Tip

origin and end will be of len 2 for 2D plane strain (x,y), and len 3 for 3D (x,y,z).

why a master config?

This is an easy way to ensure, same size arrays are set throughout all dataclasses.

Lets see the summary of the config 

config.print_summary()
If all went well you should get the following output:

~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
project: tutorial1
dim: 2
num_points: 625
num_cells: 1400
num_interactions: 10000
domain origin: (0, 0)
domain end: (1.399999976158142, 0.6000000238418579)
dt: 3.0000000000000004e-05
total time: 12.000000000000002
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

Step 4: create the Material

Lets create a material with some initial bulk density and a particle density

# Initial bulk density material parameters
phi_0 = 0.8  # [-] initial solid volume fraction
rho_p = 1400  # [kg/m^3] particle (skeletan density)
rho_0 = rho_p * phi_0  # [kg/m^3] initial bulk density

Now the material dataclass

material = hdx.ModifiedCamClay(
    config=config,
    nu=0.3,
    M=0.7,
    R=1,
    lam=0.0186,
    kap=0.0010,
    rho_p=rho_p,
    p_t=1.0,
    ln_N=jnp.log(1.29),
    phi_ref_stack=phi_0 * jnp.ones(config.num_points),
)

Lets break this down:

  • hdx.ModifiedCamClay is the Modified Cam Clay material class
  • We always pass config along when creating HydraxMPM dataclasses.
  • The only non-parameter input is the reference solid volume fraction.

Tip

Materials in HydraxMPM are normally initialize either with a reference pressure or reference solid volume fraction. Given one, we can find the other.

Tip

Material may be used in MPM analysis or single integration analysis. Sometimes its good to do both

Step 5: create the Nodes and Particles data classes

We initialize particle stresses to match the known pressure state, given a predefined solid-volume fraction above.

Here we finally make use of JAX vmap, to get the stress tensor. 

def get_stress_ref(p_ref):
    return -p_ref * jnp.eye(3)

stress_stack = jax.vmap(get_stress_ref)(material.p_ref_stack)

Pass all positions and stresses to the Particles dataclass. 

particles = hdx.Particles(
    config=config, position_stack=position_stack, stress_stack=stress_stack
)
We can create background grid nodes via the config.

nodes = hdx.Nodes(config)

The discretize function determines initial particle volume by dividing the number of particles in a cell by the cell size. 

particles, nodes = hdx.discretize(
    config=config, particles=particles, nodes=nodes, density_ref=rho_0
)

Step 6: create the Gravity and Domain

Gravity is slowly ramped up 

stop_ramp_step = config.num_steps

increment = jnp.array([0.0, -9.8]) / stop_ramp_step

gravity = hdx.Gravity(config=config, increment=increment, stop_ramp_step=stop_ramp_step)

Creating the outside domain box.

box = hdx.NodeLevelSet(config, mu=0.0)

Step 7: create the Solver

The solver determines how the background grid and material points interact.

solver = hdx.USL_ASFLIP(config)
Tip

We recommend using the hdx.ASFLIP for granular materials

Step 8: Run gravity pack

print("Start gravity pack")
carry, accumulate = hdx.run_solver_io(
    config=config,
    solver=solver,
    particles=particles,
    nodes=nodes,
    material_stack=[material],
    forces_stack=[gravity, box],
    callbacks=(
        hdx.io_vtk_callback(
            config,
            particle_output=(
                "stress_stack",
                "phi_stack",
            ),
        ),
    ),
)

solver, particles, new_nodes, material_stack, forces_stack = carry
print("Gravity pack done")

Step 9: Modify script to incorporate collapse

config = config.replace(
    origin=jnp.array([0.0, 0.0]),
    end=jnp.array([1.2, 0.5]),
    project="t2_collapse",
)

solver = hdx.USL_ASFLIP(config)

nodes = hdx.Nodes(config)

gravity = hdx.Gravity(config=config, gravity=jnp.array([0.0, -9.8]))

box = hdx.NodeLevelSet(config, mu=0.7)

Step 10: Granular column collapse

carry, accumulate = hdx.run_solver_io(
    config=config,
    solver=solver,
    particles=particles,
    nodes=nodes,
    material_stack=material_stack,
    forces_stack=[gravity, box],
    callbacks=(
        hdx.io_vtk_callback(
            config,
            particle_output=(
                "pressure_stack",
                "phi_stack",
            ),
        ),
    ),
)

print("Collapse done")