dynamics_variational_integrators Module

Provides a maximal-coordinate variational integrator for constrained rigid-body systems. The implementation follows the discrete Euler-Lagrange equations and graph factorization described by Brudigam, Sosnowski, Manchester, and Hirche (2023).

Each body contributes three world-frame translational velocities and three body-frame angular velocities to the Newton system. Holonomic equality constraints contribute scalar Lagrange multipliers. The dense solver is useful as a reference implementation; the graph-factorized solver eliminates six-variable body nodes and scalar constraint nodes according to the sparsity graph of the Newton matrix.


Uses


Contents

Variables

Abstract Interfaces

Derived Types

Subroutines


Variables

Type Visibility Attributes Name Initial
integer(kind=int32), public, parameter :: VI_DENSE_SOLVER = 1

Selects conventional dense LU factorization of the Newton matrix.
integer(kind=int32), public, parameter :: VI_FORCE_IMPLICIT_ENDPOINT = 2

Evaluates applied loads at the trial next state inside Newton.
integer(kind=int32), public, parameter :: VI_FORCE_LEFT_ENDPOINT = 1

Evaluates applied loads at the current state, explicitly.
integer(kind=int32), public, parameter :: VI_FORCE_MIDPOINT = 3

Evaluates applied loads at the interpolated midpoint state inside Newton. This selects force quadrature, not a different conservative state update or formal integrator order.
integer(kind=int32), public, parameter :: VI_GRAPH_FACTORIZED_SOLVER = 2

Selects graph-ordered block LDU factorization of the Newton matrix, including graph fill generated by Schur-complement updates.

Abstract Interfaces

abstract interface

  • public subroutine variational_constraint(state, value, args)

    Evaluates the holonomic equality constraints for a supplied maximal-coordinate state.

    Arguments

    Type IntentOptional Attributes Name
    type(variational_state), intent(in) :: state

    The state at which to evaluate the constraints.
    real(kind=real64), intent(out), dimension(:) :: value

    The constraint residual vector, which is zero for a constraint-compatible state.
    class(*), intent(inout), optional :: args

    Optional user-supplied data.

abstract interface

  • public subroutine variational_constraint_jacobian(state, jacobian, args)

    Computes the reduced maximal-coordinate constraint Jacobian. Translational columns are ordinary position derivatives. Rotational columns use the local quaternion variation

    Arguments

    Type IntentOptional Attributes Name
    type(variational_state), intent(in) :: state

    The state at which to evaluate the Jacobian.
    real(kind=real64), intent(out), dimension(:,:) :: jacobian

    The nconstraint-by-(6*nbody) Jacobian. Columns are ordered [dx, quaternion-vector variation] for each body.
    class(*), intent(inout), optional :: args

    Optional user-supplied data.

abstract interface

  • public subroutine variational_force(t, state, force, torque, args)

    Computes applied forces and body-frame torques at the selected force-evaluation state. Depending on the integrator settings, this is the current state, trial next state, or interpolated midpoint state. Velocity-dependent forces are reevaluated as Newton's trial state changes for the latter two modes.

    Arguments

    Type IntentOptional Attributes Name
    real(kind=real64), intent(in) :: t

    The current simulation time.
    type(variational_state), intent(in) :: state

    The current maximal-coordinate state.
    real(kind=real64), intent(out), dimension(:,:) :: force

    The 3-by-nbody world-frame force array.
    real(kind=real64), intent(out), dimension(:,:) :: torque

    The 3-by-nbody body-frame torque array.
    class(*), intent(inout), optional :: args

    Optional user-supplied data.

Derived Types

type, public ::  variational_integrator

Implements the maximal-coordinate variational integrator in Brudigam et al. (2023), equations (18)-(20). The nonlinear system is solved with damped Newton iterations and either dense LU or the paper's graph-structured block factorization.

Components

Type Visibility Attributes Name Initial
type(variational_integrator_settings), public :: settings

Numerical settings used by the integrator.

Type-Bound Procedures

procedure , public :: solve => vi_solve Function

Computes the solution.
procedure , public :: step => vi_step Subroutine

Advances a maximal-coordinate rigid-body state by one step.

type, public ::  variational_integrator_info

Reports convergence diagnostics for a step or trajectory solve.

Components

Type Visibility Attributes Name Initial
logical, public :: converged = .false.

True when the requested solve completed successfully.
integer(kind=int32), public :: iterations = 0

Newton iterations used; aggregated across steps by solve.
logical, public :: jacobian_singular = .false.

True if any attempted Newton Jacobian was detected as singular, including one later regularized successfully.

type, public ::  variational_integrator_settings

Defines numerical settings for the nonlinear variational step.

Components

Type Visibility Attributes Name Initial
real(kind=real64), public :: constraint_rotation_scale = 1.0d0

The characteristic dimensionless quaternion-tangent magnitude used when perturbing rotational components in a numerical constraint Jacobian.
real(kind=real64), public :: constraint_translation_scale = 1.0d0

The characteristic translation magnitude used as the absolute floor when perturbing position components in a numerical constraint Jacobian. This value has the same length units as the model positions.
real(kind=real64), public :: finite_difference_step = 1.0d-7

The relative forward-difference step used for numerical Jacobians.
integer(kind=int32), public :: force_evaluation = VI_FORCE_LEFT_ENDPOINT

The applied-load evaluation point, used for every force and torque callback, including springs, dampers, and user loads. VI_FORCE_LEFT_ENDPOINT uses the current state explicitly; VI_FORCE_IMPLICIT_ENDPOINT uses the trial next state inside Newton; VI_FORCE_MIDPOINT uses the interpolated midpoint state inside Newton. The latter two modes reevaluate loads as the Newton trial changes. These modes change load quadrature, not the conservative state update or its formal order.
integer(kind=int32), public :: linear_solver = VI_DENSE_SOLVER

The linear solver used for each Newton correction. Valid values are VI_DENSE_SOLVER and VI_GRAPH_FACTORIZED_SOLVER.
integer(kind=int32), public :: maximum_iterations = 50

The maximum number of Newton iterations allowed per time step.
integer(kind=int32), public :: maximum_line_search_iterations = 12

The maximum number of step halvings allowed by the residual line search.
real(kind=real64), public :: tolerance = 1.0d-10

The Euclidean residual tolerance used to terminate Newton's method.

type, public ::  variational_state

Defines the maximal-coordinate state of a collection of rigid bodies. Angular velocities and inertia tensors are expressed in each body's frame. Positions, translational velocities, forces, and orientations use the world frame.

Components

Type Visibility Attributes Name Initial
real(kind=real64), public, allocatable, dimension(:,:) :: angular_velocity

Body-frame angular velocities, dimensioned 3-by-nbody.
type(quaternion), public, allocatable, dimension(:) :: orientation

Unit body-to-world orientation quaternions.
real(kind=real64), public, allocatable, dimension(:,:) :: position

Body center-of-mass positions, dimensioned 3-by-nbody.
real(kind=real64), public :: time = 0.0d0

The time associated with the state.
real(kind=real64), public, allocatable, dimension(:,:) :: velocity

Center-of-mass velocities, dimensioned 3-by-nbody.

Subroutines

public subroutine initialize_variational_state(state, nbody)

Allocates and initializes a maximal-coordinate state. Positions and velocities are set to zero, and every orientation is set to the identity quaternion.

Arguments

Type IntentOptional Attributes Name
type(variational_state), intent(out) :: state

The state to initialize.
integer(kind=int32), intent(in) :: nbody

The number of rigid bodies represented by the state.