4.5. Numerical Scheme

4.5.1. Velocity Verlet Scheme

The time integration in exaDEM is performed using the velocity form of the Störmer-Verlet time integration algorithm, well-known as the velocity-Verlet algorithm. It is advantageous for DEM simulations due to its simplicity, stability, and accuracy. It provides a straightforward method for integrating Newton’s equations of motion, efficiently calculating both positions and velocities. The velocity-verlet algorithm is integrated using the following scheme at each time step:

Calculate the position vector at full timestep:

\[\mathbf{x} \left( t + \Delta t \right) = \mathbf{x} \left( t \right) + \mathbf{v} \left( t \right) \Delta t + \mathbf{a} \left(t\right)\frac{\Delta t}{2}\]

Calculate the velocity vector at half timestep:

\[\mathbf{v} \left( t + \frac{\Delta t}{2} \right) = \mathbf{v} \left( t \right) + \mathbf{a} \left( t \right) \frac{\Delta t}{2}\]

Compute the acceleration vector at full time-step \( \mathbf{a} \left( t + \Delta t\right) \) from the interatomic potential using the position at full timestep \( \mathbf{x} \left( t + \Delta t\right) \)
Finally, calculate the velocity vector at full timestep:

\[\mathbf{v} \left( t + \Delta t \right) = \mathbf{v} \left( t + \frac{\Delta t}{2} \right) + \frac{1}{2} \mathbf{a} \left( t + \Delta t\right) \Delta t\]

In exaDEM, the numerical scheme definition can be found in exaDEM/data/config/config_numerical_schemes.msp and the YAML block for the Velocity-Verlet scheme reads

numerical_scheme: numerical_scheme_verlet

numerical_scheme_verlet:
  name: scheme
  body:
   - push_f_v_r: { dt_scale: 1.0 }
   - push_f_v: { dt_scale: 0.5 }
   - push_to_quaternion: { dt_scale: 1.0 }
   - check_and_update_particles
   - reset_force_moment
   - compute_force_op
   - force_to_accel
   - update_angular_acceleration
   - update_angular_velocity: { dt_scale: 1 }
   - push_f_v: { dt_scale: 0.5 }

Note

Some operators are combined in operators to optimize code and avoid loading the same fields several times: combined_compute_prolog and combined_compute_epilog.

The exaNBody code provides a generic operator for 1st order time-integration purposes. For example, the file exaNBody/src/exanb/push_vec3_1st_order_xform.cpp provides 3 different variants:

push_v_r : for updating positions from velocities
push_f_v : for updating velocities from forces (i.e. acceleration)
push_f_r : for updating positions from forces (i.e. acceleration)

In addition, the exaNBody code also provides a generic operator for 2nd order time-integration purposes. For example, the file exaNBody/src/exanb/push_vec3_2nd_order_xform.cpp provides the following variant:

push_f_v_r : for updating positions from both velocities and forces (i.e. accelerations)

Since in exaDEM positions are expressed in a reduced frame, the argument xform_mode: INV_XFORM is mandatory when using any operator that updates the particles positions. The operators are described in detail in the following section.

4.5.2. Operators

In this section, we’ll explain operators and how to use them. The idea is to be able to easily reorder / add / remove operators to easily build another time integration scheme than the Verlet Vitesse scheme. For performance reasons, some operators are merged into a single operator, such as combined_compute_epilog and combined_compute_prolog.

4.5.2.1. Reset Forces and Moments

Operator Name: reset_force_moment
Description: This operator resets two grid fields: moments and forces.

Here is a YAML example:

- reset_force_moment

4.5.2.2. Update Particle Acceleration

Operator Name: force_to_accel
Description: This operator computes particle accelerations from forces and mass.

Formula:

\[f = a.m;\]

with f the sum of forces applied to the particle, a the acceleration of the particle, and m its mass.

Here is a YAML example:

- force_to_accel

4.5.2.3. Update Particle Orientation

Operator Name: push_to_quaternion
Description: This operator computes particle orientations from angular velocities and angular accelerations.
Parameter:
- dt_scale: Coefficient applied to the increment time (\(\Delta_t\))

Formula:

\[Q = Q+Q.av.\Delta_t\]

\[Q = \frac{Q}{||Q||}\]

\[av = av + aa.\frac{\Delta_t^2}{2}\]

with aa the angular acceleration, av the angular velocity, and Q the particle orientation.

Here is a YAML example:

- push_to_quaternion: { dt_scale: 1.0 }

4.5.2.4. Update Angular Velocity

Operator Name: push_to_angular_velocity
Description: This operator computes particle angular velocity values from angular velocities and angular accelerations.
Parameter:
- dt_scale: Coefficient applied to the increment time (\(\Delta_t\))

Formula:

\[av = av + aa.\frac{\Delta_t^2}{2}\]

with aa the angular acceleration, av the angular velocity, and Q the particle orientation.

Here is a YAML example:

- push_to_angular_velocity: { dt_scale: 1.0 }

Note

This operator is not (directly) used, it has been merged in the operator combined_compute_epilog

4.5.2.5. Update Angular Acceleration

Operator Name: push_to_angular_acceleration
Description: This operator computes angular accelerations.

Formula:

\[\omega = \bar{Q}.av\]

\[aa = Q.\dot{\omega}\]

\[\]

with aa the angular acceleration, av the angular velocity, I the particle inertia, and Q the particle orientation (and \(\bar{Q}\) its conjugate). To compute \(\dot{\omega}\), we need the particle moment and the particle inertia values.

Here is a YAML example:

- push_to_angular_acceleration

Note

This operator is not (directly) used, it has been merged in the operator combined_compute_epilog

4.5.2.6. Combined Prolog

Operator Name: combined_compute_prolog
Description: This is an operator that combined 3 operators:
- push_f_v_r
- push_f_v
- push_to_quaternion
Parameter:
- dt_scale: Coefficient applied to the increment time (\(\Delta_t\))

Here is a YAML example:

- combined_compute_prolog

4.5.2.7. Combined Epilog

Operator Name: combined_compute_epilog
Description: This is an operator that combined 3 operators:
- push_to_angular_accelaration
- push_angular_velocity
- push_f_v

Here is a YAML example:

- combined_compute_epilog