the Creative Commons Attribution 4.0 License.
the Creative Commons Attribution 4.0 License.
| 01 Apr 2025
Matter (v1): An open-source MPM solver for granular matter
Abstract. Simulating the mechanics and flow of granular media requires numerical methods that can handle extreme deformations, along with accurate constitutive models capable of capturing phenomena such as elasticity, shear banding, viscous behavior, compressibility, intergranular attractive forces and rate-dependent friction. In pursuit of this, this article introduces Matter, a Material Point Method (MPM) solver equipped with a range of models to describe dry and cohesive granular media. Rooted in a finite strain elastoplastic framework, this software features Drucker-Prager models, overstress models, critical state mechanics models as well as the µ(I)-rheology. This includes the recently proposed "critical state µ(I)-rheology model" for cohesive and compressible flows. Moreover, Matter provides a simple way of dealing with complex terrains and introduces a novel material-induced frictional boundary condition. Implemented in C++ with few required dependencies and parallelized on shared memory, it represents a lightweight yet computationally efficient option for laptops and desktops.
- Preprint
(3338 KB) - Metadata XML
- BibTeX
- EndNote
Status: final response (author comments only)
-
RC1: 'Comment on egusphere-2025-1157', Anonymous Referee #1, 09 May 2025
The authors present Matter, a code package written in C++ to simulate the flow and mechanics of granular media using the material point method (MPM). In general, the article is well-written with good English and describes a code package which has clear use cases (for example avalanche research and modelling) and clearly benefits the scientific community. The theoretical foundation is sound, and the chosen computational method suits well the intended uses.
The code itself seems to be well thought out in terms of organization and adheres to good coding practices. The package was easy to set up in Windows and VSCode thanks to the clear instructions, and it was simple to run the included example on granular collapse. I modified the example script for a few different cases and animated the results with the help of the included Python post-processing script. Everything ran smoothly without issues.
For the above reasons I am happy to recommend the paper to be published. I have written down a few minor comments below that I’d like the authors to consider first, however, but feel free to only make those changes which you also feel will improve manuscript.
- Introduction: Positioning Matter better in the context of existing open source MPM solvers for granular media.
- The authors describe some history of modelling granular matter and a large variety of existing open source MPM solvers, but at the moment it is not entirely clear from the Introduction how Matter differs from the other codes (is this the first code that implements the critical state \mu(I)-rheology model, or for example is the material-induced frictional boundary the main novelty?). The authors state that the existing codes “display a limited availability of constitutive laws and rheologies”, which gives a perfect opportunity to mention how Matter improves on the situation.
- The first two paragraphs of the Introduction, which go through the history of different methods for modelling granular media, could be condensed.
- The authors state that Matter is “lightweight yet computationally efficient”. I don’t doubt the claim, but it should be backed up with reasons for why it is so, or at least what the authors mean by computational efficiency. The current code is implemented using the Eigen library with dense grids which the authors state as a possible limitation for large sparse problems, but apart from using Eigen, which is known to be fast, and parallelization using OpenMP I didn’t find in the article reasons for computational efficiency.
Discussing more specifically the advantages Matter provides compared to existing (open source) codes would, in my opinion, strengthen the paper and the impact of Matter further. Highlighting the differences to other codes would also help the reader to understand in which case would it be beneficial to use Matter over, for example, the sparse-memory-encoding GPU framework described by Chen et al. (2025) that the authors cite.
- Code validation: The authors have done a set of validation tests against analytical solutions which show expected behaviour. Furthermore, the effect of various particle-grid transfer schemes is shown. These provide confidence in the code, but if there exists an “industry standard” code which implements some relevant constitutive laws and rheologies, a comparison in some simple case between the results of that code and Matter would be useful for further validation (I am not familiar enough with the field to say what the code and example should be). However, the already presented validation tests are sufficient.
- Paper organization: Section 7 on elasto-viscoplastic models currently seems out of place with regards to the rest of the paper. In my opinion it would be clearer if these models were described with the theory and other models available in the code, before validation and examples.
- Examples: I would like to see estimates of the compute times for the first two examples as well. There are compute times reported in Table 1 for a case similar to the first example (in Figure 9), but they don’t seem to be the same.
- The source code: Readability of the code is on a good level, although could always be improved for newcomers to the software, for example by being even more explicit about the variable names (as an example, the names of the different plasticity models in “plasticity.cpp”, such as DPSoft, DPVisc, DPMui, could be given more descriptive names or comments in the code on what they are, even though they are stated in the readme.md). But this is just a minor point, and it probably doesn’t make sense to start changing variable names now.
- Typos and small changes:
- Eq. (3): exponent of first m should be n+1 instead of n?
- Line 172: Check grammar (“time step is adapted for each time step”)
- Line 178: Check grammar
- Eq. (16): Function Ψ seems to not be defined anywhere?
- Eq. (29): Double check if the exponent of the last εp in the right hand side should be n instead of n+1.
- Line 293: Double citation. (After checking the References it seems they should be 2008a and 2008b)
- Line 301: The word “conservative”. Do you mean “conserving”?
- Figure 8: Grid lines would be helpful in assessing the results.
Citation: https://doi.org/10.5194/egusphere-2025-1157-RC1 -
AC1: 'Reply on RC1', Lars Blatny, 26 Aug 2025
We sincerely thank the Reviewer for taking the time to carefully read our manuscript and for providing thoughtful and constructive feedback. We greatly appreciate their effort and consideration in reviewing our work. Thanks to their input, we believe the manuscript has been significantly improved. Below, we present an overview of the changes made, along with general comments, following the order of the Reviewer’s remarks.
Introduction: The introduction in the revised manuscript has been clarified to highlight how Matter differs from the other available codes, most importantly related to its applicability for granular media, with various constitutive models (especially the ones based on the µ(I) rheology with/without cohesion and with/without compressibility), all implemented in a common finite strain elastoplastic framework. The MIBF boundary condition is also highlighted.
Following the Reviewer’s suggestion, the two first paragraphs of the introduction have been condensed into one shorter paragraph in the revised manuscript.
We have taken the opportunity to further clarify why we claim Matter is “lightweight yet computationally efficient” in Section 6 "Implementation, dependencies and performance" in the revised manuscript. Matter only has few required dependencies, a fairly straightforward code structure, and we dare to say computationally efficient due to its implementation, which also leverages OpenMP and OpenVDB for complex terrain handling. This claim is supported by the timing (and speedup) reported for various simulation examples.
Code validation: We are glad to hear that the Reviewer appreciates our efforts to verify our code through various tests. In particular, we have demonstrated convergence under discretization refinement, quantified numerical dissipation through various transfer schemes and verified the different boundary conditions. We have also added a new verification test of the frictional boundary conditions. The various constitutive laws are validated further in their respective papers. We are not aware of an “industry standard” code for the simulating large deformation in granular media as continua using the advanced constitutive laws as proposed here, but we would welcome any suggestions from the community on potential candidates for such a comparison that would enhance the code verification already provided.
Paper organization: Following the Reviewer’s suggestion, the original Section 7 has now been moved and merged with Section 2 in the revised manuscript.
Examples: Timing information has now been added to all examples in the revised manuscript.
The source code: In the README markdown file (the user manual) on the GitHub repository, there is a table listing the names used for all constitutive models alongside all parameter names used in the code. As the variable names might change (as the Reviewer suggests), we would propose keeping and maintaining this list in the user manual which can be updated accordingly at any time.
Typos and small changes:
- Eq. (3): This is actually correct as is, see e.g. the derivation in Pradhana (2017), chapter 2.5.2
- Line 172: Rewritten in the revised manuscript as “The size of the time step is adapted following the Courant–Friedrichs–Lewy (CFL) condition”
- Line 178: Rewritten in the revised manuscript as “The mass and velocity on each grid node are determined by interpolating the masses and velocities of all particles within the support of the interpolation functions. The grid mass is calculated as…”
- Eq. (16): Ψ has now been defined in Equation (4) as well as in the algorithm in the revised manuscript.
- Eq. (29): This is correct assuming an implicit return mapping. This has now been clarified in the revised manuscript.
- Line 293: The double citation has now been corrected in the revised manuscript.
- Line 301: Rewritten as “conserving the most energy” in the revised manuscript.
- Figure 8: Grid lines have been added in the revised manuscript
Citation: https://doi.org/10.5194/egusphere-2025-1157-AC1
-
RC2: 'Comment on egusphere-2025-1157', Anonymous Referee #2, 30 May 2025
This manuscript introduces Matter, an open-source Material Point Method (MPM) solver designed for modeling granular materials undergoing large deformations. The solver includes support for multiple constitutive models, interpolation functions, and particle-grid transfer schemes, and it proposes a novel Material-Induced Boundary Friction (MIBF) approach for boundary interaction. The accessibility, modularity, and broad applicability of the platform make it a valuable potential contribution to the computational mechanics and geomechanics communities.
While the technical scope is strong, the manuscript would benefit from several improvements in theoretical clarity, structure, and usability. I outline below a series of suggestions that I believe would help the authors strengthen the scientific quality and utility of the paper.
Major Comments
Clarify the finite-strain framework (Lines 5, 118). The solver is described as using a finite-strain elastoplastic formulation, but it does not state explicitly that it adopts the multiplicative decomposition of the deformation gradient. This decomposition is an important modeling choice and should be clearly mentioned, with a reference to Lee (1969) or similar foundational work.
Balance the comparison of MPM vs. SPH (Line 20). The introduction highlights the limitations of Smoothed Particle Hydrodynamics (SPH), including computational cost and difficulty enforcing boundary conditions. These challenges also apply to MPM, especially due to the need for particle-grid mapping and frictional contact handling. A more balanced comparison—either through literature references or brief performance benchmarks—would help contextualize the claimed advantages of MPM.
Discuss the implications of mass lumping (Line 114). The manuscript notes that mass lumping introduces energy dissipation, but this issue deserves further attention. In dynamic granular simulations, this dissipation can significantly impact the overall solution quality. It would be helpful to comment on when this effect becomes significant and whether the solver supports alternatives for better energy conservation.
Check the terminology for time integration schemes (Lines 151, Algorithm 1). The manuscript refers to the integration method as “symplectic Euler,” but the scheme shown is a basic explicit Euler update. Since symplectic integrators are typically semi-implicit and conserve Hamiltonian structure, this terminology seems incorrect. Please revise the description or provide justification if the method does retain symplectic properties.
Elaborate on viscoplastic modeling (Section 2.3). The viscoplastic formulation based on overstress regularization is introduced with minimal context. I recommend providing a bit more background on the motivation for using overstress and how it compares to other rate-dependent formulations. It would also help to clarify why this particular approach was selected.
Reorganize Sections 2.3 and 7 for better flow. Section 7 discusses rheology model selection, which is conceptually tied to Section 2.3. To improve readability, I suggest moving Section 7 earlier in the manuscript or merging it with Section 2.3 so that users get the full modeling picture up front.
Correct the formulation of stress invariant (Equation 34). Equation (34) implies a simplified 2D form of the first stress invariant p, but this formulation is misleading. In plane strain, for example, the out-of-plane stress component is non-zero and needs to be included in the calculation of p. Please revise the equation to reflect this (i.e. just write its definition for 3D stress tensor)
Acknowledge von Mises use in undrained soils (Section 7.1). The manuscript states that von Mises is rarely used for soil modeling. While this is true in general, it is worth noting that von Mises is often used for total stress analysis of undrained clays, especially in simplified or conceptual models. Adding this nuance would give a more accurate picture.
Provide guidance on model selection Although the paper supports various types of rheologies, it lacks sufficient guidance on when practitioners should use each one. A short summary of their strengths and appropriate use cases would be very helpful for users.
Comment on interpolation basis functions. Users can choose between quadratic and cubic B-spline shape functions, but no commentary is provided on the trade-offs. It would be useful to briefly explain this aspect.
Elaborate on velocity transfer schemes. The solver supports PIC, FLIP, APIC, and AFLIP schemes, but the manuscript doesn’t explain their relative advantages at the outset. In particular, AFLIP appears to offer good performance for granular flows. If the older schemes are retained for specific reasons (e.g., simplicity, stability), please clarify when users might choose one over the other.
Verify the MIBF scheme (Section 4.3). The MIBF scheme is a potentially valuable contribution, but there’s no verification of it. I recommend including a simple test case that demonstrates its behavior and compares it against standard boundary treatments.
Use correct terminology for verification and validation (Line 70, Section 5). Most of the results in Section 5 are verification exercises (e.g., convergence and numerical consistency), but the term “validation” is used throughout. Unless there are comparisons with experimental or field data, this should be corrected to “verification.”
Clarify how grain diameter in \mu(I) rheology is obtained (Line 415). The \mu(I) model uses an intrinsic grain diameter d_g, but the manuscript doesn’t say how this parameter is estimated or calibrated. Is it based on experimental data or empirical fitting? A sentence or two would help users set this parameter appropriately.
Convergence behavior. No refinement study is presented. I encourage the authors to include a basic convergence test to build confidence in the numerical results.
Minor Comments
- Line 30: Replace “a FEM” with “an FEM.”
- Use “computation time” instead of “compute time” throughout.
- Check the consistency of the vector/tensor notation.
Citation: https://doi.org/10.5194/egusphere-2025-1157-RC2 -
AC2: 'Reply on RC2', Lars Blatny, 26 Aug 2025
We are grateful to the Reviewer for dedicating time to carefully examine our manuscript and for offering insightful and constructive feedback. We truly value their effort and the thoughtful consideration given to our work. Their comments have helped us to substantially improve the manuscript. In the following, we provide a summary of the changes made and general responses, organized according to the order of the Reviewer’s points.
Clarify the finite-strain framework (Lines 5, 118): This has now been clarified in the abstract as well as in Section 2.3 using the suggested reference Lee (1969).
Balance the comparison of MPM vs. SPH (Line 20): Following the suggestion of Reviewer 1, the first part of the introduction has now been significantly condensed. Nevertheless, we have included a sentence later on in the introduction where MPM is more clearly and fairly contrasted to SPH: “Due to its hybrid nature, MPM comes with a larger memory footprint compared to both FEM and meshless methods such as the Smoothed-Particle Hydrodynamics (SPH). However, it does not necessarily demonstrate a poorer performance compared to SPH as such a meshless scheme requires keeping track of the neighboring discretization points”.
Discuss the implications of mass lumping (Line 114): This has now been specified more clearly: “Commonly applied in explicit MPM, mass lumping is a simplification also commonly adopted in explicit FEM schemes. This avoids having to invert the full mass matrix, also bringing about improved stability, however, it may result in increased energy dissipation”.
Check the terminology for time integration schemes (Lines 151, Algorithm 1): We have revised the description accordingly to describe the time integration as explicit. It is interesting to note that the MPM community has occasionally been using the term symplectic to describe this kind of time integration (see e.g. Berzins (2018, 2020)).
Elaborate on viscoplastic modeling (Section 2.3): This has been clarified further in the revised manuscript: “Viscoplasticity, i.e., rate-dependent plasticity, may be introduced with various approaches. Here, the general overstress model of Perzyna (1963, 1966) is considered. In this approach, a rate-independent yield surface y = 0 defines the onset of viscoplastic deformation which occurs for y > 0, and a separate expression for the plastic strain rate is imposed. This is in contrast to, e.g., the consistency model (Wang et al., 1997) where the yield surface itself evolves with the rate of deformation and this model can therefore obey all three Kuhn–Tucker relations. Interestingly, both approaches are in the end to a certain extent identical (Heeres et al., 2002).”
Reorganize Sections 2.3 and 7 for better flow: We agree with the Reviewer’s suggestion and have reorganized the section accordingly.
Correct the formulation of stress invariant (Equation 34): We thank the Reviewer for emphasizing that this needs further clarification, which we have now provided in the revised manuscript. The presented definitions of p and q are the ones used to define yield and plastic flow for all models in Matter. They follow from the common choice of assuming Poisson’s ratio ν = 0.5 for yield and plastic flow, see for example Tian (2009) Meccanica 44, 105-107. In general, under plane strain assumptions, the out-of-plane stress component τ3 = ν(τ1+τ2) which could be inserted into the 3D definitions of p and q. Interestingly, this would result in yield and plastic flow heavily dependent on Poisson’s ratio. For example, simulating a granular collapse with the simplest Drucker-Prager model would result in significantly different results by even a small change in Poisson’s ratio. The definitions of p and q used in Matter also ensure that the μ(I)-rheology models are consistent, providing the expected Bagnold solutions under steady-state flow. This has been clarified in the revised manuscript.
Acknowledge von Mises use in undrained soils (Section 7.1): This has now been included in the revised manuscript: “Often used to describe ductile metals, it has also been used to model soils under undrained conditions (Yu, 2006). However, due to its pressure-independency, its applicability for granular media in general is limited.”
Provide guidance on model selection: Following the Reviewer’s advice, we have now added a new sub-section giving a summary of the models and a brief overview of their key characteristics and for which materials they should be used (Section 2.10 in the revised manuscript). Hopefully, this will help users in choosing appropriate constitutive models for their problems.
Comment on interpolation basis functions: In the revised manuscript, we now write that “Higher-order B-splines are also considered, in particular cubic B-splines, which require an increased local support. Although cubic B-splines provide greater accuracy and smoothness within the interior domain and help mitigate numerical fracture, their performance near boundaries can degrade due to the increased local support required. Consequently, for simulations where near-boundary behavior is critical, quadratic B-splines may offer a more robust and reliable alternative.”
Elaborate on velocity transfer schemes: We have in the revised manuscript added a recommendation for using APIC or AFLIP. The latter is particularly used for situations where reducing numerical dissipation is of importance. Keeping the older schemes, PIC and FLIP, allows for assessing the sensitivity of the problem to numerical dissipation.
Verify the MIBF scheme (Section 4.3): This is now further verified in the new Figure 8 in the revised manuscript. In this case, we constructed an artificial material property depending only on time. With this material property used in the MIBF framework, there is an analytical solution to the sliding block problem. With a spatially varying friction, we are not aware of any analytical solutions to which we can compare. We remain open to any further suggestions how this can be verified.
Use correct terminology for verification and validation (Line 70, Section 5): We fully agree with the Reviewer and the term verification is now used in the revised manuscript.
Clarify how grain diameter in μ(I) rheology is obtained (Line 415): Following the μ(I)-rheology, this is simply the actual diameter of the grains (that make up the continuum) to be simulated. It is not a tunable parameter in that sense (unless the grain size is unknown).
Convergence behavior: A verification of the convergence under refinement of the resolution is provided in Figure 3.
Minor Comments:
The Reviewer’s minor grammatical suggestions have been fixed. Regarding notation, we have been consistent with using bold symbols vectors and tensors. However, in the journal’s manuscript template, it is visually difficult to distinguish bold fonts from normal fonts in the specific case of greek symbols. We will address this with the production team of the journal.Citation: https://doi.org/10.5194/egusphere-2025-1157-AC2
Viewed
| HTML | XML | Total | BibTeX | EndNote | |
|---|---|---|---|---|---|
| 923 | 192 | 17 | 1,132 | 19 | 40 |
- HTML: 923
- PDF: 192
- XML: 17
- Total: 1,132
- BibTeX: 19
- EndNote: 40
Viewed (geographical distribution)
| Country | # | Views | % |
|---|
| Total: | 0 |
| HTML: | 0 |
| PDF: | 0 |
| XML: | 0 |
- 1