Neural Stress Fields for Reduced-order Elastoplasticity and Fracture

We propose a hybrid neural network and physics framework for reduced-order modeling of elastoplasticity and fracture. State-of-the-art scientific computing models like the Material Point Method (MPM) faithfully simulate large-deformation elastoplasticity and fracture mechanics. However, their long runtime and large memory consumption render them unsuitable for applications constrained by computation time and memory usage, e.g., virtual reality. To overcome these barriers, we propose a reduced-order framework. Our key innovation is training a low-dimensional manifold for the Kirchhoff stress field via an implicit neural representation. This low-dimensional neural stress field (NSF) enables efficient evaluations of stress values and, correspondingly, internal forces at arbitrary spatial locations. In addition, we also train neural deformation and affine fields to build low-dimensional manifolds for the deformation and affine momentum fields. These neural stress, deformation, and affine fields share the same low-dimensional latent space, which uniquely embeds the high-dimensional simulation state. After training, we run new simulations by evolving in this single latent space, which drastically reduces the computation time and memory consumption. Our general continuum-mechanics-based reduced-order framework is applicable to any phenomena governed by the elastodynamics equation. To showcase the versatility of our framework, we simulate a wide range of material behaviors, including elastica, sand, metal, non-Newtonian fluids, fracture, contact, and collision. We demonstrate dimension reduction by up to 100,000X and time savings by up to 10X.


INTRODUCTION
Physical simulation plays a crucial role in computational mechanics, digital twins, computational design, robotics, animation, visual effects, and virtual reality.A crucial class of these physical simulations are those governed by the conservation of momentum equation [Gonzalez and Stuart 2008], where  is the first Piola-Kirchhoff stress,  is the deformation map,  is the initial density,  is the body force, and  ∈ Ω 0 is the reference position defined over domain Ω 0 .This partial differential equation (PDE) governs a wide range of elastoplastic behaviors.
To numerically solve this PDE, one has to spatially and temporally discretize it, e.g., via finite difference, finite element, or finite volume methods.A particularly flexible discretization framework is the material point method (MPM) [Jiang et al. 2016;Sulsky et al. 1995].MPM discretizes the spatial field via both Lagrangian particles and Eulerian grids.Thanks to this dual discretization paradigm, MPM thrives at handling large deformations, topology changes, and self-contact.
Nevertheless, MPM's versatility also comes at the cost of computation burden, in terms of both long runtime and excessive memory consumption.To obtain accurate results, MPM tracks a large number of state variables through the particles, often at the order of millions.Such a computation bottleneck significantly hinders the feasibility of deploying MPM in time-critical and memory-bound applications.Notably, MPM's high-dimension state variables also pose a challenge in applications where synchronization is required.For example, in virtual reality and cloud gaming, multiple users share the same simulated physical environment; each user's simulation state needs to be efficiently shared with others via internet streaming.Synchronizing millions of MPM particle data at frame rate is simply not possible.
We propose to solve these computational challenges via reducedorder modeling (ROM), also known as model reduction [Barbič and James 2005].ROM reduces the computation cost by training a low-dimensional latent embedding of the original high-dimensional simulation data.After training, instead of evolving the original highdimensional state variables over time, ROM only needs to time-step in the low-dimensional latent space, and synchronization between users only requires sharing the low-dimensional latent vector.The classic reduced-order, elasticity-only finite element method (FEM) [Sifakis and Barbic 2012] trains a low-dimensional embedding for the (discretized) deformation map  in eq.(1).However, the lowdimensional deformation embedding alone is not enough for MPM and elastoplasticity simulations in general.
History-dependent plasticity state variables.MPM simulations feature history-dependent effects, e.g., plastic deformations of sands or metals.The low-dimensional deformation embedding by itself is unable to determine the plasticity state variables that are crucial for MPM time-stepping.
Deformation gradients as independent state variables.MPM treats the deformation gradient as a separate state variable that evolves independently from deformation state variables.Again, the lowdimensional deformation embedding cannot capture these deformation gradients.
Our key observation is that the ultimate purpose of all these additional state variables is computing the stress field  in eq. ( 1).As such, we can bypass the need to capture these intermediate state variables by directly training a low-dimensional embedding for the stress field itself.The low-dimensional stress and deformation embedding together capture all the information necessary for MPM time-stepping.We construct the low-dimensional stress embeddings via implicit neural representations, also known as neural fields.Our neural stress field (NSF) approach enables stress evaluation and, in turn, force evaluation at arbitrary spatial locations.In a similar vein, we build low-dimensional neural deformation fields.To support MPM's affine particle-in-cell transferring scheme [Jiang et al. 2015], we also build low-dimensional neural affine fields for the affine momentum field.All these three neural fields share the same latent space.
After training, we solve new physical simulation problems via projection-based latent space dynamics [Benner et al. 2015;Carlberg et al. 2017].During this PDE-constrained latent space dynamics stage, we obtain computation savings by evaluating the neural fields only at a small spatial subset, similar to the idea of cubature [An et al. 2008].Our general, stress-based ROM approach works with any problem governed by the momentum equation eq.(1).To showcase the versatility of our approach, we validate NSF on a wide range of elastoplastic phenomena, including elastica, fracture, metal, sand, non-Newtonian fluids, contact, and collision.We demonstrate dimension reduction of 100,000× and computation savings of 10×.

RELATED WORK
The Material Point Method.Sulsky et al. [1995] introduced MPM by combining Lagrangian and Eulerian techniques for solid mechanics, drawing upon the earlier works by Brackbill and Ruppel [1986]; Harlow [1962] on PIC/FLIP.Since its introduction to the graphics community [Hegemann et al. 2013;Stomakhin et al. 2013], MPM has garnered considerable attention.Its primary advantage in modeling elastoplastic materials lies in its capability to handle extreme deformation and topological changes.MPM has been successfully applied to simulate various phenomena, including granular media [Chen et al. 2021;Daviet and Bertails-Descoubes 2016;Klár et al. 2016;Yue et al. 2018], non-Newtonian fluids [Fei et al. 2019;Yue et al. 2015], viscoelasticity [Fang et al. 2019], fracture [Wang et al. 2019;Wolper et al. 2020Wolper et al. , 2019]], and thermomechanics [Ding et al. 2019].Efforts have been made to speed up MPM simulations through GPU [Fei et al. 2021;Gao et al. 2018;Wang et al. 2020b], multi-node [Qiu et al. 2023], and multigrid [Wang et al. 2020a] accelerations, as well as compiler optimization [Hu et al. 2019].However, the substantial computational cost and memory consumption of MPM still present challenges that need to be addressed.
Reduced-order Modeling.Classic reduced-order modeling methods employ linear subspaces [Barbič and James 2005;Sifakis and Barbic 2012].These subspaces are often constructed via principal component analysis and, equivalently, proper orthogonal decomposition [Berkooz et al. 1993;Holmes et al. 2012].These linear subspaces have been successively applied to solids [An et al. 2008;Barbič and Zhao 2011;Kim and James 2009;Xu et al. 2015;Yang et al. 2015] and fluids [Kim et al. 2019;Kim and Delaney 2013;Treuille et al. 2006;Wiewel et al. 2019].Recently ROM methods have been exploring nonlinear low-dimensional manifolds, often leveraging autoencoder neural networks [Lee and Carlberg 2020].These nonlinear approaches enable smaller latent space dimensions in comparison with the classic linear approaches [Fulton et al. 2019;Shen et al. 2021].Our technique also falls into this nonlinear model reduction category.
Relatedly, there has been lots of progress in data-driven latent space dynamics [Lusch et al. 2018], and the entire latent space evolution is strictly learned via another neural network, e.g., recurrent neural networks [Wiewel et al. 2019].By contrast, our method follows the classic, invasive ROM literature and evolves the latent space using the numerical methods and PDEs that were used to generate the training data.In our method, the latent space dynamics are entirely PDE-based without any data-driven component.
Neural Fields.A neural field [Xie et al. 2021] parameterizes a spatially dependent vector field via a neural network.The pioneering works by Chen and Zhang [2019]; Mescheder et al. [2019]; Park et al. [2019] employ this representation for signed distance fields, where different latent space vector corresponds to different geometries.Since then, it has been widely adopted for neural rendering [Mildenhall et al. 2020], topology optimization [Zehnder et al. 2021], geometry processing [Aigerman et al. 2022;Yang et al. 2021], and various PDE problems [Chen et al. 2022;Raissi et al. 2019].Recently, Chen et al. [2023a,b]; Pan et al. [2023] have leveraged neural fields for ROM.Notably, Chen et al. [2023a] build a neuralfield-based, reduced-order framework for MPM.Their approach constructs a low-dimensional embedding only for the deformation field.Consequently, their method is unable to handle historydependent plasticity, and the deformation gradient computed from differentiating the learned deformation field is too inaccurate for large deformation phenomena such as a fracture.As a major point of departure, we train a low-dimensional manifold directly for the stress field and can therefore handle both plasticity and fracture.Furthermore, we achieve angular momentum conservation by training a low-dimensional neural affine field while [Chen et al. 2023a]'s formulation suffers from excessive dissipation.

BACKGROUND: FULL-ORDER MPM
This section will briefly recap the essential ingredients of the fullorder MPM model.Sections 4 and 5 will introduce the corresponding reduced-order model.We refer to Jiang et al. [2016]; Sulsky et al. [1995] for additional MPM details.

Finite strain elasticity and elastoplasticity
Let Ω 0 ⊂ R 3 denote the material space and Ω  the world space at time  .We are interested in the dynamics of a continuum in time  ∈ [0, ].The deformation map  :=  (, ) maps  ∈ Ω 0 to world space coordinate  ∈ Ω  .From the Lagrangian view, the dynamics of a continuum can be described by a density field (, ) : Ω 0 × [0, ] → R and a velocity field  (, ) =  (, )  : Ω 0 × [0, ] → R 3 .They are governed by the conservation of mass (, ) (, ) = (, 0), (2) and the conservation of momentum Here  = det( ),  =   (, ) is the deformation gradient,  is the first Piola-Kirchhoff stress, and  is the gravity term. can be related to the Kirchhoff stress  as  =   − .
For a hyperelastic solid, the Kirchhoff stress can be computed as where  is the energy density function of the chosen constitutive model.For an elastoplastic continuum, the deformation gradient is multiplicatively decomposed into  =     , with the former being the elastic deformation that supplies elastic force, and the latter being the permanent plastic deformation gradient.The decomposition requires that  (  ) lies within an admissible region defined by some yield condition  () < 0. Given  ,   evolves from  , following some plastic flow until the yield condition is satisfied.The procedure is often called return mapping.

MPM discretization
MPM discretizes a continuum bulk into a set of Lagrangian particles , and discretizes time  into a sequence of timesteps  0 = 0,  1 ,  2 , ...Here we take a fixed stepsize Δ, so   = Δ .The advection is performed on particles so eq.( 2) is naturally satisfied.If we approximate   by 1 Δ ( +1 −   ), and assume no gravity and free surface for clarity, for an arbitrary test function , the weak form of eq. ( 3) is then given by Pushing forward the integral from Ω 0 to Ω  = Ω   , we obtain where ,   ,  +1 and  are the Eulerian counterparts of ,   ,  +1 and , respectively [Jiang et al. 2016].
MPM adopts B-Spline-based interpolations and uses material particles  as quadratures to approximate the integration eq. ( 5).Let   denote the mass of particle  with initial position   .Denote its position and velocity at time   by    and    .Let   and   denote the mass and velocity on background grid node  at position   .Let  () denote the weight function, and    =     −   .Employing mass lumping, we can express the force equilibrium as thus providing a way to update the next stage grid velocities  +1  .Here  0  and    are the initial volume and Kirchhoff stress at time   of material particle .

MPM algorithm
At each step, particle mass and momentum are transferred to grid nodes.Grid velocities are updated and then transferred back to particles for advection.Let    denote the affine momentum of particle  at time   .The explicit MPM algorithm can therefore be summarized as the following: (1) P2G.Transfer mass and momentum from particles to grid as ).
Here  is the B-spline degree, and Δ is the Eulerian grid spacing.If additional damping is desired, RPIC can be added in the computation of   as in [Fang et al. 2019].

REDUCED-ORDER MODEL: KINEMATICS
To reduce the full-order MPM model, we will construct a nonlinear approximation to the solution of eq. ( 3) over a low-dimensional manifold.A schematic illustration is shown in fig. 2.

Low-dimensional Manifold Construction
Let the continuous field  (, ; ) : Ω 0 × [0, ] → R  denote any relevant state variable in the solution to eq. ( 3) for  ∈ Ω 0 at time .Example state variables include the deformation map, stress, etc.
Here,  is the generalized problem parameter, including but not limited to material parameters, initial conditions, and boundary conditions.Choice of  for each experiment will be detailed in section 6.We seek a continuous field f (•; x) defined over Ω 0 and parameterized by x ∈ L, a low-dimensional latent space, such that The dimension  of L ⊂ R  is taken to be a small number so that the dynamics of a continuum becomes the evolution of the latent space vector x in a low-dimensional latent space L. For notational simplicity, we will omit explicit dependence on .To computationally construct any of these low-dimensional manifolds, we will employ a neural field, also known as implicit neural representation [Xie et al. 2021].Next, we will discuss specific MPM state variables for which we will build neural fields.

Neural Deformation Fields
Similar to classic elastic-only FEM, one must build a low-dimensional manifold for the deformation field  (, ) [Barbič and James 2005].
We achieve this by constructing a manifold (, x) [Chen et al. 2023a] such that

Neural Stress Fields
Unlike elasticity-only FEM, MPM features additional history-based plastic effects and state variables.Moreover, the deformation gradient  is treated as an evolving state variable independent of .To address these various state variables, we observe that representing the stress field is a neat yet effective approach.Since the eventual goal of all these state variables is computing the stress tensor, by directly building a low-dimensional manifold for the stress tensor, we avoid cumbersome treatment of numerous plasticity state variables as well as inaccurate calculation of deformation gradient.We approximate the Kirchhoff stress field  (, ) by a manifold (, x) such that The right hand side of eq. ( 5) can thus be approximated as which naturally fits within the spatial discretization of MPM. Remarks.
(1) An alternative approach is to use the deformation gradient to compute the stress.The deformation gradient can be computed by differentiating the neural deformation field [Chen et al. 2023a].However, the numerical deformation gradient  MPM in the full-order MPM is not computed from   , but rather numerically integrated.Consequently, this approach will cease to provide accurate grid forces when   does not resemble  MPM , e.g., in numerical fracture.A well-trained neural stress field, on the other hand, directly supplies the correct grid forces for MPM grid update.(2) Since stress is computed from the elastic part of the deformation gradient   = returnMap( trial ), the plastic flow is implicitly stored.Evaluation of the return map can be avoided in deployment, thus reducing the computational cost.Overall, our neural-stress-field approach is a general approach that allows for reduced-order solutions for all the standard plasticity models.

Network training
(2) Denote the latent space vectors obtained from the encoder above as x =    (  ).Here  * denotes network weights.Additional training details and network architecture are listed in the supplementary material.
If the problem parameter  contains information about return mapping, we can make the stress decoder explicitly depend on , i.e., (, x, ).
Together with the three neural networks, we have equipped ourselves with all ingredients needed to perform one step of MPM algorithm.

REDUCED-ORDER MODEL: DYNAMICS
After training, we can run new simulations by time-stepping in the latent space L, from x to x+1 .For this, we follow the projectionbased ROM approach by Chen et al. [2023a].Our projection-based ROM approach takes three steps: (1) network inference, (2) MPM time-stepping, and (3) network inversion.The pipeline is shown in fig. 3.As we will see, since the dimension of the manifold  is much much smaller than that of the full order problem |P |, only a small subset S ⊂ P of particles, which are named sample particles, are needed to determine x dynamics.Nevertheless, due to the non-local nature of MPM, time integration of this subset will involve a larger subset N ⊂ P, which we refer to as integration particles.Note that S ⊂ N and  ≤ 3|S| < 3|N | ≪ |P |.These sample and integration particles bear similarities to the cubature points often employed in reduced-order FEM [An et al. 2008].Their exact choice will be deferred to section 5.4.

Network inference
At timestep   , given x , the states for all initial location  ∈ Ω 0 , and in particular for the integration particles  ∈ N with initial position   can be obtained by inferencing the neural networks ).Note that the particle velocity here is obtained by backward differencing the position field, consistent with the explicit MPM framework.

MPM time-stepping
One step of the MPM algorithm (section 3.3) is performed on the integration particles N to advance to  +1 .Integrating all the particles belonging to N guarantees that the states on sample particles  +1  |  ∈ S are the same as if we perform the full-order MPM on all particles  ∈ P.There is no approximation in this step.

Network Inversion
With the new particle positions at  +1 in hand, we are able to find the corresponding x+1 by inverting the neural deformation field, In this optimization problem, both the unknown x+1 and the number of summands |S| are significantly reduced.As the latent space trajectory generally evolves smoothly, with x as an initial guess, eq. ( 10) can be rapidly solved via the Gauss-Newton method [Nocedal and Wright 1999], converging in 2-3 iterations.We can optionally further speed up this nonlinear solver via a first-order Taylor approximation [Chen et al. 2023a].

Construction of Sample and Integration Particles
The least-squares problem is well-posed provided |S| ≥  /3.The projection will be more accurate if a decent number of sample particles can reflect the deformation of the geometry.For example, there should not be a group of sample particles that stand still

EXPERIMENTS
We validate the proposed reduced-order framework on a wide range of elastoplastic examples.The choice of the problem parameter  is stated in each experiment.The Experimental statistics are summarized in table 1.
In addition to visual results, we will also report the total relative deformation error across space and time,

Fracture
One remarkable feature of our neural stress field is its ability to capture fracture.We first simulate the tearing of a piece of bread with |P | = 4 × 10 4 particles governed by pure elasticity under different Young's moduli.The problem parameter  is the Young's modulus of the material.Weak elements are inserted in the middle region to seed the fracture.Figure 5(a) shows that our method is able to accurately generate the fracture pattern under various unseen Young's moduli.In MPM, numerical fracture happens when two (or more) sets of particles cannot see each other via the grid, after which the deformation gradient  for the two (or more) fractured pieces evolve independently.Thus, in this scenario,  computed from  other hand, explicitly equips the reduced-order model with the bona fide stress  that is used in the ground truth MPM simulation.
Here S consists of 450 randomly chosen particles  ∈ P, and the total relative deformation error is  = 1.2%, as is defined in eq. ( 11).The error for the baseline method is  = 6.5%.
Our neural stress field is also applicable to fracture with plastic models, such as von Mises plasticity, as is shown in cake cutting in fig. 1.Here we adopt the plasticity model in [Wang 2020].The cake is simulated with |P | = 2 × 10 5 particles.A spatula is slicing the cake at different angles, represented by the problem parameter .We select 700 particles clustered toward the middle and then reduce the sample size to 400 after  2 .The number of integration particles is 1.35 × 10 4 on average.The full-order and reduced-order MPM simulators are both implemented in WARP [Macklin 2022] under double precision.The neural networks are implemented in PyTorch.The total wall clock time of the full-order simulation is 14.495, while the wall time of our reduced method is 1.417.We achieve an overall speedup of 10.23× with an error of 1.3%.In general, since the dynamics are constrained to the low-dimensional manifold, we are also able to take a larger time step (1.5Δ .) at deployment time.In both fracture examples we choose  = 6, and  = 2 × 10 4 and 1 × 10 6 , respectively.With our reduced model, we are also able to achieve considerable memory saving.In this scenario, the average memory consumption of the full-order MPM model is 1.61G, while ours is 0.79G, including both latent space physics and neural networks.The computing setup is detailed in the supplementary material.

Sand plasticity
MPM is particularly suitable for simulating granular media.We simulate a column of sand falling onto the ground under gravity.
Here  represents different friction angles.Our neural stress field can perfectly capture such a noisy stress distribution and yields excellent results on D test , with an average error of  = 0.4%.(See fig.7) The ground truth is simulated with |P | = 72, 000 particles, while we set  = 6,  = 3.6 × 10 4 and |S| = 150.The memory consumption of full-order MPM for this scenario is 0.91G, and that of our reduced model is 0.65G.
Once trained using a low-resolution simulation, our approach can arbitrarily boost the resolution with no cost by simply evaluating the neural deformation fields at more   ∈ Ω 0 .In fig.6, we boost the resolution by 100× when running latent space dynamics.

Metal plasticity
Our neural stress field can also handle history-based plasticity models, such as the effect of hardening [Wang et al. 2019].The squeezing and bouncing back of a metal frame is simulated with von Mises return mapping under different hardening coefficients  =   .In the ground truth simulation, the yield condition is  () < 0, and thus the return mapping is constantly updated to account for hardening, de facto making the yield condition another path-based state  (•) =  (•, ).Since our neural stress field directly approximates the stress computed after the return mapping, such complexity is circumvented.In other words, the hardening state is implicitly learned by our neural stress field .In fig.8, we compare our deployment results and ground truth under different hardening coefficients.A sampling of |S| = 50 particles out of |P | = 49, 978 yields a remarkably small error of  = 0.2% averaging over all testing data, where we choose  = 5 and  = 29.987.While endto-end ML frameworks [Sanchez-Gonzalez et al. 2020] can only predict particle positions at rollout time, our PDE-based reducedorder model captures various physical quantities beyond positions.Indeed, our neural stress field can also accurately predict the stress distribution, as is shown in fig. 9 Furthermore, we can sample even fewer points to still obtain reasonably good results.Sampling only 20 particles results in an error of  = 0.5%, while sampling merely 30 particles results in an error of  = 0.3%, and the results are almost indistinguishable visually compared with the ground truth.In addition, with a randomly chosen 30 sample particles, and with the timestep in deployment set to 1.5Δ, we are able to speed up the total wall clock time from 7.83 in the full-order MPM to 1.55 in the reduced model, achieving a speedup more than 5 × .In this setup, the full-order MPM memory consumption is 0.82G, while ours is 0.51G.
We follow the Herschel-Bulkley model in [Yue et al. 2015].With Figure 12: A jelly cube hits onto a jelly wall.Our approach accurately reflects the rotation of the cube.The baseline approach [Chen et al. 2023b] is much more dissipative since it does not support angular momentum.The error  for our approach is 0.20%, and for the baseline approach is 16.6%.The problem parameter  represents different initial velocities of the jelly cube.
just 50 sampling points, we can predict the dynamics of toothpaste with an averaging total relative deformation error  = 0.6%.Further, the sample size can be even reduced without too much discount on the overall visual quality.As is shown in fig.11, with only 30 points, the deployment result still looks reasonably good, with an error of  = 0.9%.

Rotation and Collision
We simulate a collision scenario that yields salient rotation (fig.12) with |P | = 10 4 particles.With a manifold dimension of  = 5 and Compared with the baseline [Chen et al. 2023b], ours accurately captures both the self-contact and the rotation.The error  for our approach is 0.19%, and for the baseline approach is 4.7%. represents different inclinations of the plane.
|S| = 50 sample particles, our approach is able to accurately capture the rotational dynamics.The baseline approach, nevertheless, suffers from noticeable artifacts due to its flawed representation of stress and affine fields.Our approach can also phonograph complex contact scenarios (fig.13).We simulate an elastic squishy ball falling onto an inclined plane with |P | = 10 5 particles.The manifold dimension is set to  = 6.A randomly selected set of |S| = 300 sample particles suffices to delineate the contact of tentacles.The baseline method performs poorly as the affine momentum is missing, and the representation of stress is inaccurate in extreme contact.Notice that we do not need to sample all tentacles to capture their motion; rather, a small S is used to determine x in the latent space so that our neural deformation and neural stress field can generate their motion.The error  for either of the above experiments is less than 0.2%.The dimension reduction ratios are 6, 000 and 5 × 10 4 .

DISCUSSIONS AND FUTURE WORK
We proposed Neural Stress Fields (NSF), a novel, reduced-order framework for elastoplastic and fracture simulations.NSF significantly alleviates the computational burden of simulating complex elastoplasticity and fracture effects by training a unified, lowdimensional latent space for the neural deformation, stress, and affine fields.Following the training phase, we efficiently conserve computational resources by leveraging these low-dimensional latent variables for evolution.Our approach sets a compelling precedent for multiple potential research trajectories.
Generalization.Our work supports both interpolation and extrapolation of the training data (see experiments on sand friction angles and bread weak elements).Nevertheless, our approach cannot handle extremely out-of-distribution extrapolation.We trade aggressive generalizability for massive compression and speedup.Future work may consider exploring alternative balancing between generalizability and performance.In addition, for each experiment, we train a network using data from this particular scenario [Sifakis and Barbic 2012].An exciting future direction is training on one scenario but generalizing to multiple materials and objects.
Training time.Currently, training time is long, between 2hrs and 20hrs.Our target applications are cases where the model would be re-used multiple times.For example, after training, our model can be deployed in VR and gaming applications, where millions of users will interact with it.In these cases, training time is not the main bottleneck.That said, improving training time will help capture larger scenes and accelerate development cycles.
High-frequency neural fields.MPM simulations often involve stress fields with high-frequency details and large spatial variations.In practice, we find it challenging to train neural fields that correctly capture these distributions, preventing us from capturing larger scenes.Future work may consider developing more advanced neural architectures [Sitzmann et al. 2020;Tancik et al. 2022] to improve performance when high-frequency details are presented.
Path-dependent plasticity.Our latent space vector x is only determined by the position   .Nevertheless, since plasticity is pathdependent [Borja 2013], the same position field does not imply the same stress field.A potential fix to this issue would be to, instead of training two distinct networks, concatenate  and  and train (  , ([  ,   ])) ≈ [   ,    ].In addition, to more explicitly enforce history dependency, future work may consider evolving the latent space according to both stress updates and deformation updates.
Data-free training.Sharp et al. [2023] introduces a data-free reduced-order modeling framework by incorporating a physicsinformed loss term.Extending it to include MPM's plasticity and fracture phenomena is another exciting direction.Figure 18: Pie-chart breakdown for the runtime of each component in our reduced-order pipeline for two examples: the cake and the metal.The pie chart for the metal example is plotted smaller to show that this is a smaller problem and thus has a shorter runtime.Overall, the main bottleneck of the reduced-order algorithm remains to be MPM timestepping while the network operations introduce little overhead.
In the example of cake, the runtime of full-order MPM is 4.53/100 frames.The runtime of the reduced-order MPM is 0.469/100 frames.The cost in network inference, projection, as well as communication between WARP and PyTorch is 0.195/100 frames.
In the example of metal, the runtime of full-order MPM is 1.305/100 frames.The runtime of the reduced-order MPM is 0.209/100 frames.The cost in network inference, projection, as well as communication between WARP and PyTorch is 0.178/100 frames.
We also provide two pie charts for the breakdown of the runtime of each component in our reduced-order scheme for the two examples.The runtime for network-related operations grows relatively mildly when problem size increases, whereas the runtime for MPM timestepping in theory grows linearly in problem size (see next paragraph for a discussion).Thus, problems with larger scales tend to enjoy more time savings from our algorithm.
For a paralleled explicit MPM simulator, one can usually observe that the runtime is roughly linearly proportional to the total number of particles.The two majority of costs come from the SVD computation for stress and the atomic addition in P2G.In our experiments, reduced MPM with about 9% of original particles costs about 16% of original runtime, and with 6.7% of original particles costs about 10.3% of original runtime.In theory, this reduced MPM runtime should be even smaller.In addition, SVD is not required in the reduced MPM.We reason that this is because, during deployment, PyTorch has already occupied some memory.Thus, the computing power allocated for the reduced MPM solver in WARP would be less than that allocated for a full-order MPM solver alone.

D TRAINING DETAILS
The training data is generated from an MPM solver written in WARP [Macklin 2022] under double precision.The Adam optimizer [Kingma and Ba 2014] for stochastic gradient descent is used for training.The Xavier initialization is used for the ELU layers.We fix the learning rates to be (10 −3 , 5 × 10 −4 , 2 × 10 −4 , 10 −4 , 5 × 10 −5 ).For the Neural Deformation Field, 300 epochs are trained for the learning rate above.For the Neural Stress Field and Neural Affine Field, 600 epochs are trained for the learning rate above.The batch size for the three manifolds is  • |P | where we choose  to be the largest value such that the training data fits within memory or 32, whichever is smaller.The training data is normalized to have mean zero and variance one.The whole training pipeline is implemented in Pytorch Lightning [Falcon 2019].

E NETWORK ARCHITECTURE
The (input, output) dimension pairs for ,  and  are (3 + , 3), (3 + , 6), and (3 +, 9), respectively.The Kirchhoff stress  has 6 degrees of freedom since it is symmetric.Each MLP network contains 5 hidden layers, each of which has a width of  • , where  = 3 for ,  = 6 for , and  = 9 for . is a hyperparameter, the exact value of which for each experiment is listed in Table 1 in the main text.We adopt the ELU activation function [Clevert et al. 2015].The encoder network is devised as the following: several 1D convolution layers of kernel size 6, stride size 4, and output channel size 3 are applied until the length of the 1D output vector reaches or below 12.The vector is then reshaped to 1 channel.One MLP layer transforms its dimension to 32, followed by the last MLP layer that outputs a vector with dimension  .

F ELASTICITY AND PLASTICITY DETAILS
We first list all parameters that shall be needed in discussing the models below.In all plasticity models used in our work, the deformation gradient is multiplicatively decomposed into  =     following some yield stress condition.A hyperelastic constitutive model is applied to   to compute the Kirchhoff stress  .For a pure elastic continuum, one simply takes   =  .

F.1 Fixed corotated elasticity
The Kirchhoff stress  is defined as where  =    and   =    is the singular value decomposition of elastic deformation gradient.[Jiang et al. 2015] F

Figure 2 :
Figure 2: Latent space kinematics.Given a latent space vector x ∈ L, evaluating the neural deformation, stress, and affine fields at any reference position   ∈ Ω 0 (e.g., the black dot in Ω 0 ) results in the corresponding deformation, stress, and affine momentum at time  at the current position (the boxed dot in Ω  ).

Figure 3 :
Figure 3: Latent space dynamics.We time-step the latent space via three steps.Each step involves a small spatial subset S ⊂ N ⊂ P of the original full-order MPM particles.

Figure 4 :
Figure 4: (a) Construction of sample particles and integration particles; (b) Sample and integration particles across the domain.
.., , ∈ P ||(  , x ) −  (  ,   )|| 2 =1,2,.., , ∈ P || (  ,   )|| 2 .(11) Throughout this section, dataset D is always split as non-overlapping D train and D test .Neural fields are constructed with D train and validated on D test .Furthermore, we will report the dimension reduction ratio defined by  = 3|P |/ , i.e., the dimension of the full-order model divided by the latent space dimension.See the supplementary material for additional details regarding experiments, the training dataset, generalizability, extrapolation, and elastoplastic models.

Figure 5 :
Figure5: Tear a piece of bread.Our method accurately captures the tearing behavior at different elastic moduli.Due to a lack of accurate stress representation and the inaccurate deformation gradients computed from neural fields, the baseline approach by[Chen et al. 2023a] fails to capture the fracturing behavior.
would provide a misleading stress  that is not what is being used in the MPM setting.As is shown in fig.5(b), the baseline method by Chen et al. [2023a], which uses  ≈  (   ), fails to reconstruct a clean fracture.Our neural stress field, on the

Figure 6 :
Figure 6: After training on low-res simulation (left), our method can directly infer high-resolution results (right) by querying the continuous neural deformation field.No additional post-processing is needed.

Figure 7 :
Figure 7: Simulate column collapse for sand under varying friction angles.

Figure 8 :
Figure 8: Our neural stress field can capture the hardening effect under different hardening coefficients.

Figure 9 :
Figure9: Unlike end-to-end ML frameworks that can only predict particle positions, our first-principal-based reducedorder approach also matches stress quantitatively.

Figure 10 :
Figure 10: A ribbon of toothpaste is smeared onto a toothbrush held at different angles.The four subplots show our deployment results with 50 sample points.The corresponding ground truth is shown in the top right corner of each subplot.

Figure 11 :
Figure 11: Results for different numbers of sampling points are shown.For this example where the toothpaste is held at 7.1 • , sampling 15 randomly chosen particles yields an error of  = 1.8%, while sampling 30 yields an error of  = 0.9%.

Figure 13 :
Figure 13: An elastic squishy ball falls onto an inclined plane.Compared with the baseline[Chen et al. 2023b], ours accurately captures both the self-contact and the rotation.The error  for our approach is 0.19%, and for the baseline approach is 4.7%. represents different inclinations of the plane.
the APIC transfer scheme is adopted.If the conventional PIC scheme is adopted, the latter is simply replaced by       =          .(2) Grid update.Update grid velocities at next timestep by  +1  =    − Δ       ∇    0  + Δ.Collision and Dirichlet boundary conditions are also handled at this stage.(3) G2P.Transfer velocities back to particles and update particle states. +1

Table 1 :
Simulation and reduction statistics.