Improved Implementation of Approximate Full Mass Matrix Inverse Methods into Material Point Method Simulations

Abstract: Approximate full mass matrix methods for the material point method, known as FMPM(k) of order k, can improve the calculation of grid velocities from grid momentum. It can be implemented in any MPM code by inserting a new calculation task whenever grid velocities are needed. The implementation recommended in this paper only needs these calculations once per time step just before when updating particle positions and velocities. FMPM implementation issues arise, however, when its methods are mixed with other MPM feature that rely on lumped mass calculations. Some common lumped-mass MPM features are grid-based, velocity boundary condition, multimaterial contact calculations, crack contact calculations, and imperfect interfaces. This paper first derives a revised FMPM(k) implementation that both simplifies and clarifies the "FMPM Loop" that can be added to MPM codes. Next, that loop is modified to allow FMPM(k) to work well even in simulations that need other MPM features that previously caused conflicts. Two other FMPM(k) issues are apparent loss of stability at very higher order k and inherent computational cost. These issues are discussed in an analysis of temporal stability as a function of order k and in consideration of options to improve efficiency.

• The paper introduces an incremental FMPM(k) approach that reduces noise and improves boundary condition enforcement in MPM simulations.
• The revised method recasts contact handling using 'Evolving' and 'Net' algorithms, ensuring artifact-free multimaterial interface responses.
• The work implements dynamic order selection and energy stabilization, balancing high-order accuracy with computational efficiency.

Improved Implementation of Approximate Full Mass Matrix Inverse Methods in MPM

Introduction and Context

This work addresses several core challenges in integrating high-order full mass matrix inverse approximations (FMPM($k$)) into the Material Point Method (MPM) framework. Building on prior developments [Nairn:198], FMPM($k$) aims to mitigate null-space noise and enhance stability and accuracy in MPM by introducing $k$-order extrapolation for grid velocities. Despite improvements documented in prior works, FMPM($k$) faced implementation obstacles when combined with standard MPM features such as lumped-mass-based boundary conditions and contact handling, encountered additional stability constraints at higher orders, and incurred significant computational cost as $k$ increased. This paper presents a revised FMPM($k$) algorithmic formulation, supporting modifications to boundary condition and contact enforcement, convergence-based dynamic order selection, and practical stability-enhancing variants.

Revised FMPM($k$) Algorithm

The standard MPM update involves mapping particle data to the grid, computing grid momenta and velocities, applying forces and contact, updating particle states, and then propagating back. Traditional methods employ a diagonal (lumped) mass matrix for efficient but approximate grid velocity computation. FMPM($k$) leverages a Taylor expansion of the full mass matrix inverse, performing $k-1$ recursively defined extrapolation iterations:

• Each iteration computes an incremental velocity correction $\Delta \mathbf{v}(\ell)$, with the total grid velocity $k$0.
• The recursion is parameterized such that each update is independent of $k$1, which is vital for modular code and facilitates new options like dynamic $k$2.

A schematic of the revised flow is shown in Figure 1.

Figure 1: Flow chart detailing a single MPM time step, highlighting integration points of FMPM($k$3) calculations under various update strategies.

This incremental approach not only clarifies the role of each FMPM iteration but is central to the improved integration of boundary conditions and contact as detailed below.

Grid Velocity Boundary Conditions

Previous conflicts between high-order FMPM and conventional velocity boundary condition (BC) imposition, arising from the inconsistent application of lumped versus full mass matrices, led to error amplification at higher $k$4. The revised algorithm applies BC corrections incrementally after each FMPM extrapolation step—enforcing that the prescribed BCs are satisfied precisely after each correction, rather than solely at the end. This enforces correct propagation of BC constraints throughout the FMPM chain.

A canonical moving wall test, using CPDI and B2CPDI shape functions, demonstrates the superior accuracy and robustness of the revised BC enforcement (Figure 2, Figure 3):

Figure 2: Particle domain under a right-edge moving wall, with BC nodes circled and adjustable wall depth $k$5.

Figure 3: Velocity error in FMPM($k$6) as a function of $k$7 (A) and as a function of wall BC depth $k$8 (B) for both new and prior BC implementations.

Results confirm that, unlike previous algorithms which degrade or stagnate at high $k$9, the incremental BC application sustains low errors and high precision for all tested orders.

Contact and Multimaterial Interface Handling

MPM frequently involves multimaterial contact assignments (e.g., interface shock propagation, frictional sliding, and crack closure). Conventional contact enforcement in MPM is predicated on lumped-mass fields. Embedding FMPM($k$0) in these contexts led to interface artifacts and loss of convergence to correct single-material results.

This work introduces two incremental algorithms for coupling contact laws to the FMPM loop: the "Evolving" and "Net" methods. Both recast the contact response as an update to the velocity increment at each extrapolation order. Benchmark shock simulations validate that the "Net" method, which always evaluates contact in terms of the uncorrected total velocity, both guarantees reversion to single-material behavior and produces accurate, artifact-free interface responses (Figure 4, Figure 5).

Figure 4: Pressure shock profile traversing a material interface under “stick” contact, showing FMPM($k$1) accuracy for $k$2 matching single-material solutions.

Figure 5: Frictional shock interface, comparing results under the “Net” and “Evolving” methods; only “Net” produces artifact-free, converged profiles.

These findings highlight that only incremental, order-independent correction preserves both the energy consistency and the contact law fidelity at arbitrary $k$3.

Numerical Stability Analysis and Enhancement

Explicit MPM stability is Courant-limited, generally requiring $k$4. FMPM($k$5)'s improved mass matrix approximation restricts the admissible $k$6 severely as $k$7 rises, with plateauing at $k$8 for large orders (Figure 12A).

Figure 6: (A) Maximum stable Courant number for various FMPM($k$9) blending strategies; (B) associated energy loss versus blending parameter $k$0.

Potential stability-enhancing variants include:

• Blending the full mass matrix with FMPM(2) rather than the lumped mass matrix, which increases stability thresholds while maintaining low dissipation.
• Periodic application of FMPM($k$1) updates, reducing per-step dissipation at low $k$2.
• Utilization of velocity-gradient extrapolations, suggesting possible stability at $k$3 for FMPM($k$4).

Crucially, the choice of energy-preserving or dissipative blend is task-dependent: insufficiently low $k$5 increases dissipation (see energy loss in two-disk impact, Figure 7), while overly aggressive blending with the lumped-mass matrix destroys the fundamental FMPM benefits.

Figure 7: Kinetic plus internal energy loss in two-disk offset impact simulations, showing rapid reduction of loss with increasing $k$6.

Dynamic FMPM and Computational Efficiency

The modular incremental structure of the new FMPM loop supports dynamic order adaptation—halting extrapolation once velocity corrections fall below a convergence threshold. Experiments on block impact and shock wave benchmarks demonstrate that, while dynamic $k$7 can moderately enhance efficiency, the practical improvement is fragile and sensitive to the convergence criterion. In many scenarios, using a fixed modest order (e.g., $k$8–5) is both robust and nearly optimal (Figure 8, Figure 9).

Figure 8: Helmholtz free energy conservation in off-center disk impact, demonstrating near-perfect conservation for both FLIP and FMPM(4).

Figure 9: Spatial and temporal convergence of energy dissipation for various update algorithms and Courant numbers.

Practical Implications and Future Directions

The methods presented ensure that approximate full mass matrix inversion techniques can be robustly integrated with advanced MPM features—without sacrificing stability, accuracy, or efficiency. The revised FMPM($k$9):

• Provides algorithmic routes for handling velocity BCs and multimaterial or crack contact at arbitrary extrapolation orders;
• Supplies practical guidance on temporal stability constraints and efficient stabilization strategies;
• Enables the practical deployment of adaptive (dynamic) high-order corrections, though the heuristic for dynamic adaptation warrants refinement.

Potential future advances include further exploration of velocity-gradient stabilization, refined dynamic convergence metrics that explicitly account for diminishing returns in high $k$0, and formal extension to more complex interface (e.g., imperfect adhesion) and transport-coupled problems.

Conclusion

This work recasts the FMPM($k$1) approach for MPM in a modular, incrementally consistent framework that systematically resolves incompatibilities with essential MPM features. The derived methods are empirically validated for high-order accuracy, stability, and physical consistency across challenging benchmark simulations. In practice, FMPM($k$2) of order 4–5 achieves substantial reduction in artificial dissipation and noise with moderate cost; implementing the revised algorithm requires minimal additional infrastructure beyond what is needed for FMPM(2). The techniques and analyses herein set a new standard for the robust, high-fidelity MPM simulation of complex boundary- and contact-rich problems.

Reference:

John A. Nairn. "Improved Implementation of Approximate Full Mass Matrix Inverse Methods into Material Point Method Simulations" (2604.07307).