compute continuum/chunk command

Syntax

compute ID group-ID continuum/chunk chunkID cutoff width field1 field2 ... keyword args ...
  • ID, group-ID are documented in compute command

  • continuum/chunk = style name of this compute command

  • chunkID = ID of compute chunk/atom command

  • cutoff = cutoff for the truncated Gaussian kernel

  • width = standard deviation of the corresponding Gaussian kernel without truncation

  • one or more input fields can be listed

  • field = natoms, density, volume/fraction, momentum/a, velocity/a, momentum/grad/ab, velocity/grad/ab, strain/rate/ab, stress/ab, stress/ke/ab, stress/contacts/ab, boundary/force/a, fabric/ab, temperature

    natoms = number of atoms in the kernel's support
    density = density field
    volume/fraction = volume fraction field
    momentum/a = a-component of the momentum field
    velocity/a = a-component of the velocity field
    momentum/grad/ab = ab-component of the momentum gradient field
    velocity/grad/ab = ab-component of the velocity gradient field
    strain/rate/ab = ab-component of the strain-rate field
    stress/ab = ab-component of the total stress field
    stress/ke/ab = ab-component of the kinetic stress field
    stress/contacts/ab = ab-component of the contact stress field
    boundary/force/a = a-component of the boundary force density
    fabric/ab = ab-component of the fabric tensor field
    temperature = granular temperature field
  • zero or more keyword/arg pairs may be appended

  • keyword = boundary/atom or boundary/fix

    boundary/atom arg = groupbound
      groupbound = group-ID for atoms that make up a boundary
    boundary/fix arg = none
      enables boundary corrections from fix wall/gran

Examples

compute cc1 all chunk/atom bin/2d x 0.0 1.0 y 0.0 1.0 units box
compute prop1 all property/chunk cc1 coord1 coord2 count
compute cont1 all continuum/chunk cc1 1.0 0.3 density velocity/* stress/*
fix 1 all ave/time 10 10 100 c_prop1[*] c_cont1[*] file continuum.dat mode vector

compute cc2 all chunk/atom bin/1d z lower 0.5 units reduced
compute cont2 flow continuum/chunk cc2 2.0 1.0 volume/fraction temperature boundary/force/* boundary/fix

Description

Added in version TBD.

Define a computation that calculates coarse-grained continuum fields for chunks of atoms on demand using the construction in (Goldhirsch). This formulation is designed to construct rigorous macroscopic continuum fields from discrete particle data, commonly used in discrete element method simulations. The fields are evaluated at the chunk centers with a truncated Gaussian kernel. This compute does no time averaging and writes no output files itself. To time average or write out the per-chunk data, use fix ave/time exactly as with other chunk-based compute styles.

In LAMMPS, chunks are collections of atoms defined by a compute chunk/atom command, which assigns each atom to a single chunk (or no chunk). The ID of that command is specified as chunkID. For example, a chunk may be a molecule or a spatial bin. See compute chunk/atom and Howto chunk for details.

This compute is only compatible with the binning styles of compute chunk/atom (bin/1d, bin/2d, or bin/3d). If binning is not performed along one of the box dimensions, outputs are normalized by the box length in that dimension when appropriate. For example, if a 3d system is binned only along z, the reported fields are normalized by Lx and Ly. Normalization is not performed on velocities, gradients, strain rates, and the number of atoms.

The available fields include scalar, vector, and tensor quantities. A vector field such as the velocity may be requested as an individual component such as velocity/x or with a wildcard such as velocity/*. This expands to all components (two in 2d, three in 3d). Tensor fields such as the stress use component names like stress/xy. Wildcards such as stress/* expand to all components (four in 2d, nine in 3d). Tensors are not assumed to be symmetric such that the xy component may differ from the yx component.

Only atoms in the specified group contribute to the calculations. In pair sums, both atoms must be in the group. The referenced compute chunk/atom command defines its own group and optional region. Atoms with chunk ID = 0 are not assigned to any chunk and do not contribute.

One may want to use this compute in conjunction with compute property/chunk to additionally obtain the positions of each chunk’s center (where kernels are centered) and fix ave/time to provide optional time averaging and data output.


The kernel used in the coarse-graining is a truncated Gaussian function. The standard deviation of the function is defined by the width parameter. The maximum cutoff of the Gaussian function is defined by the cutoff parameter. These values are not restricted to the size of a chunk and can extend beyond its boundaries. The width has effectively no impact on performance, but a larger cutoff does increase the number of atoms that add contributions to each coarse-graining site.


The natoms field is simply a count of the number of atoms within the cutoff distance of the kernel. Note that this is not equivalent to the number of atoms in the chunk (as calculated by compute property/chunk) as the kernel may not span an entire chunk or it may span beyond a single chunk.

The density field is

\[\sum_i m_i W(\vec{r}_\mathrm{chunk} - \vec{r}_i)\]

where the summation is across all atoms \(i\) in the group and within the cutoff of the kernel from the center of the chunk. \(m_i\) is the atom mass, \(\vec{r}_\mathrm{chunk}\) is the chunk center, \(\vec{r}_i\) is the atom position, and \(W\) is the kernel.

The volume/fraction field is

\[\sum_i V_i W(\vec{r}_\mathrm{chunk} - \vec{r}_i)\]

where \(V_i\) is the finite sized particle volume in 3d or area in 2d.

The momentum field is

\[\sum_i m_i v_{i,a} W(\vec{r}_\mathrm{chunk} - \vec{r}_i)\]

where \(v_{i,a}\) is the \(a\) component of the atom velocity. The momentum/grad field is then obtained with centered finite differences between neighboring chunks. Thus, if a box dimension is represented by a single chunk, gradients along that dimension are zero. Gradient values are calculated using a one sided finite difference for chunks adjacent to a non-periodic boundary. The first index is the direction of the momentum while the second index is the direction of the derivative, such that xy is calculating the derivatives of \(v_x\) in the \(y\) direction.

The velocity field is the ratio of the momentum and density fields. The velocity/grad field is then obtained with centered finite differences. Thus, if a box dimension is represented by a single chunk, gradients along that dimension are zero. Gradient values are calculated using a one sided finite difference for chunks adjacent to a non-periodic boundary. The first index is the direction of the momentum while the second index is the direction of the derivative, such that xy is calculating the derivatives of \(v_x\) in the \(y\) direction.

The boundary/force field is the interaction force density of the particles on the boundaries as defined in (Weinhart):

\[- \sum_i \sum_k f_{ik,a} W(\vec{r}_\mathrm{chunk} - \vec{r}_{\mathrm{contact},ik})\]

where the sum over \(k\) is over boundary elements and \(\vec{r}_{\mathrm{contact},ik}\) is the contact point between atom \(i\) and boundary element \(k\). At least one of the boundary/atom or boundary/fix keywords must be used to request this quantity.

The stress/ke field is the kinetic contribution to the stress:

\[-\sum_i m_i (v_{i,a} - v_{\mathrm{chunk},a}) (v_{i,b} - v_{\mathrm{chunk},b}) W(\vec{r}_\mathrm{chunk} - \vec{r}_i)\]

where \(v_{i,a}\) is the \(a\)-th component of the velocity of atom \(i\) and \(v_{\mathrm{chunk},a}\) is the \(a\)-th component of the average velocity of the chunk defined by the velocity option above. Note that this average velocity is evaluated at the center of the binned chunk which contains the atom.

The stress/contacts field is the contact contribution to the stress:

\[-\frac{1}{2}\sum_{i,j} f_{ij,a} r_{ij,b} \int_0^1 ds\, W(\vec{r}_\mathrm{chunk} - \vec{r}_i + s \vec{r}_{ij})\]

where \(f_{ij,a}\) is the force on atom \(i\) from atom \(j\) and \(\vec{r}_{ij}\) is the displacement between the two atoms. Here, the double summation is over all pairs of atoms \(i\) and \(j\) where each pair can be counted twice if both atoms are in the kernel’s cutoff.

The stress field is the sum of the kinetic and contact contributions.

The fabric field is

\[\frac{1}{2}\sum_{i,j} V_i n_{ij,a} n_{ij,b} \int_0^1 ds\, W(\vec{r}_\mathrm{chunk} - \vec{r}_i + s \vec{r}_{ij})\]

where \(V_i\) is the volume of the atom in 3D and area in 2D and \(\vec{n}_ij\) is the normal vector pointing from \(i\) to \(j\). See compute fabric for more general discussion on fabric tensors.

The strain/rate field is

\[\frac{1}{2} \left( \nabla_{ab} v + \nabla_{ba} v \right)\]

where \(\nabla_{ab} v\) is the \(ab\) component of the velocity/grad field.

The temperature field is a local granular temperature defined as:

\[\frac{1}{2} \sum_i m_i (v_i - v_\mathrm{chunk})^2 W(\vec{r}_\mathrm{chunk} - \vec{r}_i)\]

This is the kinetic energy density (energy per volume or area) and does not include any factors of density or dimension which may be included in a statistical mechanical temperature calculation.


The optional boundary/atom and boundary/fix keywords turn on the boundary corrections for stress and stress/contacts described in (Weinhart). The boundary/atom keyword designates a group of atoms as a boundary. Those atoms are removed from the atom and pair sums above, except for the boundary/force contribution. The boundary/fix keyword applies the analogous correction for boundaries created by fix wall/gran.

Output info

This compute calculates a global array where the number of rows is the number of chunks \(N_\text{chunk}\) defined by the referenced compute chunk/atom command. The number of columns is the number of requested values after wildcard expansion. The columns appear in the same order as specified in the command. Internal intermediate quantities needed to evaluate derived fields are not exposed as output columns. An array is used even if the number of requested column is one.

These values can be accessed by any command that uses global arrays from a compute as input. See Howto output for an overview of output options. The array values are intensive. Units depend on the requested fields.

Restrictions

This compute is part of the GRANULAR package. It is only enabled if LAMMPS was built with that package. See the Build package page for more information.

Only bin/1d, bin/2d, and bin/3d styles of compute chunk/atom are supported. Furthermore, the compress, limit, and discard no options of the chunk/atom compute are not supported If a boundary is periodic, results from bins on the boundary will be incorrect unless the simulation box in that dimension is evenly divisible by the bin width such that no bin extends beyond the simulation boundaries. A warning will be issued if this condition is not met. It is recommended that the user pick reasonable settings for chunks as not all permutations of settings may be covered. Tests are recommended.

Triclinic boxes are not supported.

The volume/fraction, stress, stress/contacts, boundary/force, and fabric values require particles with a radius attribute.

Pair-dependent quantities require a pair style that can compute the force for a single pair of atoms. The boundary/fix keyword requires at least one fix wall/gran instance. All such instances must have the contacts keyword enabled.

Default

No boundary corrections are applied unless boundary/atom or boundary/fix is specified.


(Goldhirsch) Goldhirsch, Granular Matter, 12, 3, 239-252 (2010).

(Weinhart) Weinhart, Thornton, Luding, Bokhove, Granular Matter, 14, 2, 289-294 (2012).