=Paper=
{{Paper
|id=Vol-2964/article_195
|storemode=property
|title=Model Reduction for the Material Point Method on Nonlinear Manifolds Using Deep Learning
|pdfUrl=https://ceur-ws.org/Vol-2964/article_195.pdf
|volume=Vol-2964
|authors=Peter Yichen Chen,Maurizio Chiaramonte,Eitan Grinspun,Kevin Carlberg
|dblpUrl=https://dblp.org/rec/conf/aaaiss/ChenCGC21
}}
==Model Reduction for the Material Point Method on Nonlinear Manifolds Using Deep Learning==
https://ceur-ws.org/Vol-2964/article_195.pdf
Model Reduction for the Material Point Method on Nonlinear Manifolds Using
Deep Learning
Peter Yichen Chen1,2 , Maurizio Chiaramonte1 , Eitan Grinspun2,3 , Kevin Carlberg1
1
Facebook Reality Labs Research
2
Columbia University
3
Toronto University
cyc@cs.columbia.edu, mchiaram@fb.com, eitan@cs.toronto.edu, carlberg@fb.com
Abstract Reduced-order model
We introduce a strategy for constructing a projection-based
This work proposes a projection-based model-reduction ap- reduced-order model (ROM) that introduces a kinematic
proach for the material point method (MPM), a popular hy-
brid Eulerian–Lagrangian method used in engineering and
approximation comprising constructing a low-dimensional
computer graphics for simulating solids, fluids, and multi- nonlinear manifold approximation to the deformation map
phase phenomena. The proposed technique employs a kine- and generates dynamics by first computing the velocity on
matic approximation of the deformation map by enforcing grid nodes of the Eulerian mesh, and then performing least-
deformation trajectories to reside on a low-dimensional man- squares projection of the computed velocity onto the tangent
ifold expressed in extrinsic form via a parameterization func- space of the manifold.
tion corresponding to a decoder. By explicitly approximating
the deformation map, the deformation gradient can be com- Kinematics: low-dimensional manifold
puted simply by differentiating the associated parameteriza- In analogue to constructing low-dimensional nonlinear man-
tion function. Further, this mapping can be inverted in order
to generate new material points as needed allowing for adap-
ifolds for finite-dimensional dynamical-system state spaces
tive refinement. The method generates dynamics via projec- (Lee and Carlberg 2020), one can restrict any element of the
tion, i.e., at every time step, the method performs two steps: reference domain Ωref ⊆ Rd to evolve on a low-dimensional
(1) compute the velocity on grid nodes of the Eulerian mesh, manifold. We first denote the approximated deformation
and (2) perform least-squares projection of the computed ve- map as φ̃ : Ωref × T × D → Rd with φ̃(·; ·, µ) ∈ S(µ),
locity onto the tangent space to this low-dimensional man- ∀µ ∈ D and
ifold. The approach supports hyper-reduction, which is es-
sential for achieving computational complexity independent φ̃(·; t, µ) : X 7→ x̃(t, µ) (1)
of the original number of material points. Numerical exam-
ples on large-scale problems illustrate the method’s ability
: Ωref → Ω̃(t, µ), ∀t ∈ T , µ ∈ D, (2)
to generate orders-of-magnitude computational-cost savings where S(µ) denotes the space of admissible solutions satis-
with negligible error. fying initial and essential boundary counditions, Ω̃(t, µ) ⊆
Rd denotes the deformed domain corresponding to the ap-
proximated deformation map at time t ∈ T and parameter
Introduction instance µ ∈ D, and enforce the kinematic constraint
The Material Point Method (MPM) (Sulsky, Zhou, and φ̃(X; ·, ·) ∈ M(X) := {g(X; ŷ) | ŷ ∈ Rr } ⊆ Rd , (3)
Schreyer 1995; Jiang et al. 2016) is a hybrid Eulerian–
Lagrangian discretization method widely employed in the ∀X ∈ Ωref , where g : Ωref × Rr → Rd denotes the parame-
computer-graphics community for solid, fluid, and multi- terization of an r-dimensional manifold.
phase simulations. Due to its dual Eulerian and Lagrangian The kinematic restriction (3) implies that there exists gen-
representations, MPM offers several advantages over the fi- eralized coordinates x̂ : T × D → Rr such that
nite element method: large deformation, fracture, as well as φ̃(X; t, µ) = g(X; x̂(t, µ)),
contact and collision are handled relatively easily. (4)
∀X ∈ Ωref , ∀t ∈ T , µ ∈ D.
However, MPM’s dual representations of the material
make it especially computationally expensive. This pre- Assuming the parameterization function is continuously dif-
cludes the ability of MPM to be used in time-critical and ferentiable in its first argument, Eq. (4) implies that the de-
real-time settings such as fast-turnaround design, interactive formation gradient of the approximate deformation map can
applications, and control. be calculated by differentiating the parameterization func-
tion as
Copyright © 2021 for this paper by its authors. Use permitted under F̃ : (X, t, µ) 7→ ∇φ̃(X; t, µ) ≡ ∇g(X; x̂(t, µ))
Creative Commons License Attribution 4.0 International (CC BY (5)
4.0) : Ωref × T × D → Rd×d .
, where ∇ denotes differentiation with respect to the first
argument. Further, assuming the parameterization function
is continuously differentiable in its second argument, Eq. (4)
implies that the velocity of the approximate solution can be
calculated via the chain rule as
˙ ∂g ˙ µ).
φ̃(X; t, µ) ≡ (X; x̂(t, µ))x̂(t, (6)
∂ x̂ Algorithm 1: ROM dynamics approach #1
Remark (Manifold construction via deep learning). While
the manifold-parameterization function g : Ωref × Rr → Rd Input: Generalized coordinates x̂(tn ) and velocities
˙ n ).
x̂(t
can be constructed in a variety of ways, we employ a fully-
connected, 5-layer, deep-learning architecture (i.e., a multi- Output: Generalized coordinates x̂(tn+1 ) and
layer perceptron) with ELU activation functions due to their ˙ n+1 ).
velocities x̂(t
b
continuous differentiability. ELU can be viewed as a smooth 1 Identify the grid nodes I ⊆ {1, . . . , n } needed to
extension of the more popular ReLU activation function. We compute dynamics for the sample material points,
set g ← gθ? where θ? , the network weights, are the (approx- i.e., I = {i | Ni (xpn ) 6= 0 for any p ∈ M}, where
imate) solutions to the problem Ni is the Eulerian basis function.
p p
minimize 2 Compute the deformation gradient Fn , velocity vn ,
θ, {x̂(t,µ)}t∈Ttrain , µ∈Dtrain p
and position xn for each sample and neighbor
X material point p ∈ M ∪ N at time instance tn by
(kgθ (X; x̂(t, µ)) − φ(X; t, µ)k22
evaluating (4)–(6) for X = X p , p ∈ M ∪ N and
µ∈Dtrain , t∈Ttrain , X∈Ωref,train
t = tn .
+ λ1 k∇gθ (X; x̂(t, µ)) − ∇φ(X; t, µ)k2F 3 Perform the ‘particle to grid’ transfer by computing
∂gθ for i ∈ I
+ λ2 k ˙ µ) − φ̇(X; t, µ)k2 ),
(X; x̂(t, µ))x̂(t,
∂ x̂ 2 X J(F p )
σ n
(7) fi,n =− σ(Fnp )∇x Ni (xpn ) mp
ρ0
p∈M∪N
where Dtrain ⊆ D, Ttrain ⊆ T , and Ωref,train ⊆ Ωref denote
parameter, time, and reference-domain instances at which e
X J(Fnp )
fi,n = b(xpn )Ni (xpn ) mp ,
the original (full-order) MPM has been solved and solutions ρ0
p∈M∪N
˙ µ) = x̂(t+∆t,µ)−x̂(t,µ) is computed via
are available, x̂(t, ∆t
explicit Euler with ∆t being the time step size of MPM, and where J := det(F ), σ denotes the Cauchy stress, b
λ1 , λ2 ∈ R+ denote penalty parameters for the deformation denotes body forces, mp is the material point mass,
gradient and velocity, respectively. and ρ0 is the initial material density.
4 Perform the update step by computing for i ∈ I
Dynamics: velocity projection
To generate projection-based dynamics at each time in- 1
v̇i,n+1 = (f σ + fi,n
e
)
stance, the proposed ROM (1) computes the velocity on grid mi i,n
nodes of the Eulerian mesh, and (2) performs least-squares ∆vi,n+1 = v̇i,n+1 ∆tn
projection of the computed velocity onto the tangent space to vi,n+1 = vi,n + ∆vi,n+1 .
the low-dimensional manifold. However, to ensure the com-
putational complexity is independent of the number of origi- 5 Perform the ‘grid to particle’ transfer by computing
nal material points np , we must restrict these steps to execute for p ∈ M ∪ N
on only a subset of domain. Toward this end, we propose two p
X
separate hyper-reduction approaches. vn+1 = vi,n+1 Ni (xpn ).
i∈I
Hyper-reduction approach #1: material-point projection
The first approach we propose is material-point centric. We 6 Update the generalized coordinates
begin by identifying a priori a set M ⊆ {1, . . . , np } of ˙ n+1 ), where x̂(t
x̂(tn+1 ) = x̂(tn ) + ∆tn x̂(t ˙ n+1 )
sample material via stochastic sampling. We also consider satisfies the minimization problem
the neighbors of these material points N ⊆ {1, . . . , np }, de-
fined as the set of material points needed for the computation X ∂g p
min k (X p ; x̂(tn ))x̂(tn+1 ) − vn+1 k22 .
of the dynamics of the sample material points for all t ∈ T x̂(tn+1 ) ∂ x̂
p∈M
and µ ∈ D; note that M ∩ N = ∅. With these material
points identified, Algorithm 1 provides the associated steps
taken by the resulting (online) reduced-order-model simula-
tion.
While this approach can incur an operation count inde-
pendent of the original number of material points np and Eu-
lerian basis functions nb , it requires computing (and track-
ing) the set of neighboring material points, which is difficult
to do in practice. The next method mitigates this issue by
adopting a sample-node-centric perspective.
Hyper-reduction approach #2: sample-node projection
This section presents an alternative that leverages the Algorithm 2: ROM dynamics approach #2
fact that—if constructed to be injective—the manifold- Input: Generalized coordinates x̂(tn ) and velocities
parameterization function g can be inverted for positions in ˙ n ).
x̂(t
the reference domain given the position in the deformed con- Output: Generalized coordinates x̂(tn+1 ) and
figuration and values of the generalized coordinates. Thus, ˙ n+1 ).
velocities x̂(t
the key elements of this approach are: (1) no need to track b
1 Select (any) sample nodes I ⊆ {1, . . . , n } of
individual material points, (2) the ability to generate effec-
interest that satisfy Ni (x) 6= 0 for some x ∈ Ω.
tive material points as needed, and (3) performing classical
2 Define a quadrature rule comprising quadrature
finite-element quadrature on the Eulerian grid to assemble
points and weights xqn ∈ Ω, wnq ∈ R+ ,
the governing equations.
q = 1, . . . , nq used to assemble the governing
Regarding quadrature, we first note that
equations at the sample nodes I.
ne Z ne Xnq q
3 Compute the deformation gradient Fn and velocity
Z X X
f ρdV = f ρdV ≈ f (xq )ρ(xeq )wq . vnq for each quadrature point q = 1, . . . , nq at time
Ω e=1 Ωe e=1 q=1 instance tn by evaluating (4)–(6) for X = Xnq ,
where ne is the number of elements, nq is the number of q = 1, . . . , nq and t = tn , where Xnq satisfy
quadrature points, and wq is the quadrature weight. Consider xqn = g(Xnq ; x̂(tn )), q = 1, . . . , nq .
4 Perform the ‘particle to grid’ transfer by computing
approximating the integrals from the weak-form of the equa-
tion of motion using classical techniques as for i ∈ I
q
n
nq
ne X
σ
X J(F q )
n
fi,n =− σ(Fnq )∇x Ni (xqn ) wq
X X
( ρv̇j Nj Ni )|xq wq ρ0
q=1
e=1 q=1 j
nq nq
ne X
e
X J(Fnq )
(8) fi,n = b(xqn )Ni (xqn ) wq .
X
=− (σ(F )∇Ni − bNi )|xq wq ρ
q=1 0
e=1 q=1
Z
5 Perform the update step by computing for i ∈ I
+ tNi , i = 1, . . . , nb .
∂Ω
1
In contrast to the first approach, this approach (Algorithm v̇i,n+1 = (f σ + fi,n
e
)
mi i,n
2) does not require tracking any fixed material points nor
∆vi,n+1 = v̇i,n+1 ∆tn
their neighbors. Rather, by selecting in step 1 a set of sample
nodes I and leveraging invertibility of the deformation map, vi,n+1 = vi,n + ∆vi,n+1 .
the approach can incur an operation count independent of 6 Select (any) sample material points xsn ,
the original number of material points np and Eulerian basis s = 1, . . . , ns that satisfy the condition: if
functions nb without any special tracking. Ni (xsn ) 6= 0, then i ∈ I.
7 If needed, perform the ‘grid to particle’ transfer by
Numerical experiments computing for s = 1, . . . , ns
Unit tests X
s
For unit testing, we shear the top of a cylinder made of elas- vn+1 = vi,n+1 Ni (xsn ),
tic material for one timestep. i∈I
Training with gradient penalty In Table 1, we compare where vns can be computed by evaluating Eq. (6) for
different training setups. With both λ1 and λ2 set to be zero, X = Xns , s = 1, . . . , ns , where Xns satisfies
the network is trained only with the position information xsn = x̂(Xns ; x̂(tn )), s = 1, . . . , ns .
without first-order gradient information. While the position 8 Update the generalized coordinates
error and deformation gradient error of the reduced-order ˙ n+1 ), where x̂(t
x̂(tn+1 ) = x̂(tn ) + ∆tn x̂(t ˙ n+1 )
simulations are relatively small, a large velocity error is ob- satisfies the minimization problem
served. Training with a positive λ2 while keeping λ1 zero, n s
the network is now trained with additional temporal gradi- X ∂g ˙ n+1 ) − v s k2 .
ent information (velocity). Consequently, the velocity error min k (X s ; x̂(tn ))x̂(t n+1 2
˙ n+1 )
x̂(t s=1
∂ x̂
is significantly improved, which also leads to significant im-
provement in the position error. However, a slight increase in
the deformation gradient is observed. We thereby train with
both the spatial (deformation gradient) and the temporal (ve-
locity) gradient information by setting both λ1 and λ2 set to
Training Position Velocity Deformation
Parameters Error Error Gradient
Error
λ1 = 0, λ2 = 0 0.19% 13% 0.2%
λ1 = 0, λ2 = 0.01 0.004% 0.3% 0.3%
λ1 > 10000, λ2 > 0.01 0.004% 0.2% 0.2%
Table 1: Training parameters’ influence on the accuracy of
the reduced-order simulations. Training with both spatial Figure 1: The material is poked at the top by different forces,
(λ1 ) and temporal (λ2 ) gradient penalty yields the best re- resulting in different deformed state.
sult. The material is discretized with 3 material points.
Specimen information Value Unit
Projection Position Velocity Deformation
Scheme Error Error Gradient Radius 10 cm
Error Height 40 cm
E 12500 Pa
Material-Point ν 0.3
0.01% 0.06% 0.11%
-Centric
Sample-Node Table 3: Poking test specimen information
0.01% 0.06% 0.11%
-Centric
Material-Point
-Centric, θ ∈ [0◦ , 9◦ ], and φ ∈ [10◦ , 15◦ ]. A total number of 27 simu-
0.01% 0.16% 0.11%
Hyper-reduction lations are generated via full factorial sampling, 22 of which
(6) for training while 5 for testing. Each simulation consists of
Sample-Node 36 timesteps or 0.25s. Therefore, a total of 945 simulation
-Centric, snapshots is used for training and testing.
0.01% 0.13% 0.11%
Hyper-reduction There are np = 1, 368 material points in the full-order
(6) simulation. The dimension of the generalized coordinate is
chosen to be 5, effectively reducing the dimension of the
Table 2: Projection approaches and hyper-reduction. Both simulation by a factor of 821. The number of material points
projection approaches achieve the same level of accuracy for hyper-reduction is 15; these points are chosen from the
with or without hyper-reduction. The material is discretized portion where the poking force is applied, the bottom where
with 17 material points. the cylinder is fixed, as well as the middle section where the
material is allowed to deform freely.
The reduced-order simulation is able to reproduce the
be positive. As expected, optimal errors across position, ve- training cases with a position error of 0.50% and is able to
locity, and deformation gradient are observed. predict the testing cases with a 0.55% error.
Material-point-centric projection vs. sample-node-
centric projection In Table 2, we compare the reduced- Future work
order simulation’s accuracy when using different projection In the future, we would like to extend our work to support
methods. All the simulations use the same network with the plasticity, fracture, contact, and collision.
same weights. Both the material-point-centric approach and
the sample-node-centric approach achieve relatively small References
error. Jiang, C.; Schroeder, C.; Teran, J.; Stomakhin, A.; and Selle,
Hyper-reduction Table 2 also demonstrates the effective- A. 2016. The material point method for simulating contin-
ness of hyper-reduction where we use only a small subset uum materials. In ACM SIGGRAPH 2016 Courses, 1–52.
of the original material points for projection, which leads to Lee, K.; and Carlberg, K. T. 2020. Model reduction of dy-
cheaper computational cost. While the original number of namical systems on nonlinear manifolds using deep convo-
material points is np = 17, using 6 for projection yields the lutional autoencoders. Journal of Computational Physics
same level of position error. 404: 108973.
Sulsky, D.; Zhou, S.-J.; and Schreyer, H. L. 1995. Applica-
Poking test tion of a particle-in-cell method to solid mechanics. Com-
An elastic cylinder, whose information is listed in Table 3, puter physics communications 87(1-2): 236–252.
is poked at the top (Figure 1). The poking force is char-
acterized by its spherical coordinate, f (ρ, θ, φ). We gener-
ate training and testing data by varying the three parame-
ters characterizing the poking force, where ρ ∈ [1.0, 1.2],