Accurately Solving Rod Dynamics with Graph Learning
Authors: Han Shao
Tassilo Kugelstadt
Torsten Hädrich
Wojciech Pałubicki
Jan Bender
Soren Pirk
Dominik Michels
Venue:  (2021)
Keywords: Dynamical Systems
Representation Learning

Abstract: Iterative solvers are widely used to accurately simulate physical systems. These solvers require initial guesses to generate a sequence of improving approximate solutions. In this contribution, we introduce a novel method to accelerate iterative solvers for rod dynamics with graph networks (GNs) by predicting the initial guesses to reduce the number of iterations. Unlike existing methods that aim to learn physical systems in an end-to-end manner, our approach guarantees long-term stability and therefore leads to more accurate solutions. Furthermore, our method improves the run time performance of traditional iterative solvers for rod dynamics. To explore our method we make use of position-based dynamics (PBD) as a common solver for physical systems and evaluate it by simulating the dynamics of elastic rods. Our approach is able to generalize across different initial conditions, discretizations, and realistic material properties. We demonstrate that it also performs well when taking discontinuous effects into account such as collisions between individual rods. Finally, to illustrate the scalability of our approach, we simulate complex 3D tree models composed of over a thousand individual branch segments swaying in wind fields.

Under review as a conference paper at ICLR 2021
ACCURATELY SOLVING ROD DYNAMICS WITH GRAPH LEARNING
Anonymous authors Paper under double-blind review
ABSTRACT
Iterative solvers are widely used accurately simulate physical systems. These solvers require initial guesses to generate a sequence of improving approximate solutions. In this contribution, we introduce a novel method to accelerate itera- tive solvers for rod dynamics with graph networks (GNs) by predicting the initial guesses to reduce the number of iterations. Unlike existing methods that aim to learn physical systems in an end-to-end manner, our approach guarantees long- term stability and therefore leads to more accurate solutions. Furthermore, our method improves the run time performance of traditional iterative solvers for rod dynamics. To explore our method we make use of position- based dynamics (PBD) common solver for physical systems and evaluate by simulating the dynam- ics of clastic rods. Our approach is able to generalize across different initial condi- tions, discretizations, and realistic material properties. We demonstrate that it also performs well when taking discontinuous effects into account such as collisions between individual rods. Finally, to illustrate the scalability of our approach, we simulate complex 3D tree models composed of over a thousand individual branch segments swaying wind fields.
INTRODUCTION
The numeric simulation of a dynamic system commonly comprises the derivation of the mathemati- cal model given by the underlying differential equations and their integration forward in time. In the context of physics-based systems, the mathematical model is usually based first principles and depending on the properties of the simulated system, the numerical integration of a complex sys- tem be very resource demanding (Nealen al. 2006), e.g., hindering interactive applications. Enabled by the success of deep neural networks to serve as effective function approximators, re- searchers recently started investigating the applicability of neural networks for simulating dynamic systems. While many physical phenomena can well be described within fixed spatial domains (e.g.. in fluid dynamicc) that can be learned with convolutional neural network (CNN) architectures (Chu & Thuereyy 2017] Guo et al. [0166] Tompson et al. 2016 ct al. 22200) a large range of physical systems can more naturally
represented as graphs. Examples include systems based on connected particles (Müller et al., 2007), coupled oscillators (Michels & Desbrun 2015) Michels et al.) 2014), or finite elements (Nealen et al., (2066). Existing methods enable learning these systems often in an end-to-end manner and with focus on replacing the entire part of the integration procedure. A number methods show initial success in approximating physical systems; however, they often fail to reliably simulate the state of system over longer time horizons if significant disadvantages not accepted, such as the use of large datasets containing long-term simulations and the employment of specific memory structures (Sanchez-Gonzalez al. 2020)-
In this paper, we aim to improve the performance of iterative solvers for physical systems with graph networks (GN). An iterative solver requires an initial guess, and based on it generates sequence of improving approximate solutions. The initial guess can be computed by simply using values obtained in the previous iteration or by solving a few steps of similar but simpler physical system. The performance of an iterative solver is significantly influenced by the calculation of the initial guess, which we aim to replace with the prediction of a GN. To demonstrate our approach, we use position-basee dynamics (PBD) solver that approximates physical phenomena by using sets of connected vertices (Bender et al. 2017 2014b; Macklin et al.) 2016 Müller al., 2007). To simulate a physical system, PBD first computes updated locations of vertices using symplectic Euler integration to then correct the initial position estimates so as to satisfy a set of predefined constraints. The correction step is known as constraint projection and is commonly solved iteratively. The explicit forward integration for predicting the system's 'ss updated state has negligible cost, whereas the projection step is computationally expensive. Our goal is to employ a GN to predict the outcome of the constraint projection step as an initial guess. This way, our approach inherits the long-term stability of classic PBD solver, while providing better run-time performance.

---

Under review as a conference paper at ICLR 2021
Figure/ 1: Renderings taken from real-time simulations of the elastic deformation of a helix falling down on the ground plane (left) and two rods colliding with each other (right).
To showcase the capabilities of our combined PBD solver, we aim to simulate the physically plau- sible mechanics of elastic rods. Rods play important role for a variety of application domains, ranging from surgical simulation of sutures (Feess et al. 2016), catheters, and tendons (Pai 2002). to human hair (Michels et al. (015) and vegetation (Pirk al. 2017) in animated movies. Furthermore, approaches exist to realistically simulate rods as sets of connected vertices accurately capturing their mechanical properties (Bergou et al., 2008 Kugelstadt & Schoemer, 2016 Michels et al. 2015 Pai 2002). Our approach is able to generalize across different initial conditions, rod discretizations, and realistic material parameters such as Young's 's modulus and torsional modulus (Deul et al. 2018). Moreover, we demonstrate that our approach can handle discontinuous collisions between individ- ual rods. Figure shows examples of elastic rod deformations of a helix falling down (left) and two colliding rods (right). Finally, we show that the data-driven prediction of the initial guesses of the constraint projection leads to a decreased number of required iterations, which in turn - results in significant increase of performance compared to canonical initial guesses.
In summary, our contributions are: (1) we show how to accelerate iterative solvers with GNs; (2) we show that our network-enabled solver ensures long-term stability required for simulating physical systems; (3) we showcase the effectiveness of our method by realistically simulating elastic rods; (4) we demonstrate accuracy and generalizability of approach by simulating different scenarios and various mechanical properties of rods including collisions and complex topologies (dynamic tree simulations).
2 RELATED WORK
In the following we provide overview of the related work, which spans from data-driven physics simulations and graph learning to position-basee based dynamics and clastic rods.
Data-driven Physics Simulations. It has been recognized that neural networks can be used as effec- tive function approximators for physical and dynamic systems. To this end, early approaches focus on emulating the dynamics of physies through learned controllers (Grzeszczuk al. 1998) by de-
signing subspace integrators (Barbič & James 2005). Today, a range of approaches exist that enable
learning ordinary and partial differential equations (Lagaris et al. 1998 Raissi et al. 2019 Raissi
& Karniadakis 2018), for example, to transform them into optimization problems (Dissanayake
& Than-Thien 1994). to accelerate their computation (Mishra 2018 Sirignano & Spiliopoulos/
2018). or to solve for advection and diffusion in complex geometries (Berg & Nyström 2018).
Other methods focus on specific data-driven solutions for non-linear elasticity (Ibañez et al. 2017
for approximating Maxwell's equation in photonic simulations (Trivedi et 2019), or for ani- mating cloth (Wang et al.] (01T), partially focusing on interactive applications (Holden al., 2019). More recently, research on data-driven approaches for modeling the intricacies of fluid dynamics has gained momentum (Ladický et al.) 2015; Ummenhofer et al. 2020). Due to fixed-size spatial rep- resentation of Eulerian fluid solvers, number of approaches rely on CNN-type architectures (Chul & Thuerey 2017: Guo et al.] 2016 Tompson et al. 2016; Xiao et al.] 2020).. Furthermore, it has been shown that data-driven approaches
even be used to approximate the temporal evolution of fluid flows (Wiewel et al. 2018), to compute liquid splashing (Um et al. (017). artistic style- transfer (Kim et al. 2020). or derive fluid dynamics from reduced sets of parameters (Kim al., 2019).
Graph-based Learning. Graphs have proven to be a powerful representation for learning wide range of tasks (Battaglia et al., 2018 Scarselli et al., 2009). In particular, it has been shown that graphs enable learning knowledge representations (Kipf al. 2018). message passing (Gilmer et al.] 2017), or to encode long-range dependencies. e.g., as found in video processing (Wang let al. 2017). A variety of methods uses graph- representations to learn properties of dynamic physical systems. e.g. for climate prediction (Seo & Liu 2019). with an emphasis on individual objects (Chang et al. 2016) and their relations (Sanchez-Gonzalez et al. 2018), for partially observ- able systems (Li et al. 2018), the prevalent interactions within physical systems (Kipf et al. 2018). particle systems (Mrowca et al.) 2018). or more generally physical sim- ulation (Sanchez- Gonzalez et al. 2019 2020). While many of the existing approaches learn the time integration of physical systems in an end-to-end manner, we use a graph network to predict the outcome of PBD solver for rod dynamics to enable more efficient computations.
hierarchically-orgaanl
2

---

Under review as a conference paper at ICLR 2021
Pt-1
Ai
A
A1-1
Discretized Rod
Forward Integration
Constraint Projection
Figure 2: Illustration of the discretization of single rod using several rod segments arranged along its centerline (left). Each rod segment is described by its position and orientation within the gener- alized coordinates p,. The Lagrange multipliers X; represent the interaction between rod segments. The forward integration path is illustrated in red (middle) and constraint projection in green (right).
Position-based Dynamics and Elastic Rods. PBD is a robust and efficient approach for simulating position changes of connected sets of vertices (Bender et al. 2017 7; 201 4b Macklin et al.) 2016 Müller et al. 2007). Compared to forced-bases methods that compute the force directly, the inter- action between different vertices in PBD is realized by constraint projection step in iterative manner. To avoid the dependency of the system's 's stiffness on the number of iterations and the time step size, an extended position-based based dynamics approach was introduced (XPBD) (Macklin al.) 2016). A number of methods exist that model the dynamic properties of rods that can even simulate more complicated rod mechanics 2002). Moreover, particle systems were employed simulate the dynamics of rods (Michels et al. 2017) and, in particular, for the physically accurate simulation of thin fibers (Michels et al.j (2033) such as present in human hair or textiles. On a different tra- jectory, it has been recognized that rods can be simulated based on PBD (Umetani et al.) 2014). The initial formulation was improved (Kugelstadt & Schoemer 2016) by including the orientation of rod segments in the system's ' state to account for torsion effects. Later, the XPBD framework was utilized (Deul et al.) 2018) to address the non-physical influence of iteration numbers and steps sizes, which enables the more accurate simulation of elastic rods.
METHODOLOGY
We propose a novel approach to simulate the temporal evolution of a dynamic system which consists of elastic rods. Each rod is discretized using several rod segments arranged along its centerline (Fig- ure2). For each rod segment, its state is described by its position, orientation, velocity and angular velocity. The state of the system is given the set of the individual states of all rod segments. The simulation is carried out by employing PBD (uüller et al., 2007) directly manipulating the sys- represented quaternions allowing for convenient implementation of bending and twisting effects (Kugelstadt & Schoemer 2016). Extended PBD (XPBD) (Macklin) et al. 2016) is implemented avoid that the rods' stiffnesses depends on the time step size and the number of iterations (Deul et al. (0018).
tem's state. Orientations
The generalized coordinates of a rod segment at time t is given by Pa, E R x H. which includes its position x,, E R given in Cartesian coordinates and its orientation described by a quaternion qr.c H.. Correspondingly, v E RSS refers the generalized velocity of the rod segment, which includes velocity and angular velocity. The system is continuously updated during the simulation by applying corrections Ap; = (Ax,, A T E R with position shifts Ax; E R and orientation shifts Ad; E R3 representing the integration of the angular velocity
A single time integration step is presented in Algorithm In the beginning (lines generalized velocity and generalized coordinates are updated by employing symplectic Euler integration step. In this regard, aеx denotes the generalized acceleration due the external net force, e.g., given by gravity. XPBD (Macklin et al., 2016) employs the Lagrange multiplier À which is initialized as zero (line(5) along with the integrated generalized coordinates p" Collision detection results are stored in Coll,, and Collep (line G). where Coll,, includes all the pairs of two rod segments that poten- tially collide with each other and Coll.p includes information of all rod segments that potentially collide with another object such plane. Within several solver iterations, we alternate between rod constraint projection and the collision constrain projection (lines 7 to 15). The rod constraints include shear-stretch and bend-twist constraints representing the corresponding elastic energy. The Lagrange multiplier represents the interaction between rod segments. Figure
2 illustrated the dis- cretization for single rod into several interacting segments. The correction values Ap and AA in line arre computed by constraint projection (Deul et al., 22018) Kugelstadt & Schoemer, 2016). The generalized coordinates and Lagrange multipliers are updated for each rod (lines (8lnes|81tt to 12). and rod- rod and rod-plane collisions are addressed to update the generalized coordinates. For details about the collision projection procedure, we refer to Macklin et al. (2014). For the non-collision case, the steps within line 6 and 13 are not needed.
Pleasee note, that Aq. GØg)AØ, E RR,, in which the matrix G(g) RX3 describes the relationship hetween quaternion velocity and angular velocity (Bender al. 2014a).
3

---

Under review as a conference paper at ICLR 2021
Algorithm Numerical integration procedure updating pi, + Pi,++Ar and Wi,t H Ar.
for all segments do
2: vj - Vi.t + At Best
3:
p + Pi.t + At of with H(q.e) = [1x3, Oax3; 0xx3, G(q.)))
4: end for
5: A 0.000pp
(Collar, Coll.p) < gemeraeCCollssinnnntttrnnnnniii p*)
7:
for jorj <- 0 number of required solver iterations do
8:
for all rods do rodConsttrinnnProecttonn// X*)
9:
(Ap, AA)
10:
1+1 --
AA
11:
p end for
p'
+ Ap
12:
13:
p 1+1
uppaacCulliioCCotttaanPPriiiiiii Coll.e Coll.p)
14:
j j+i+1
15: end for
16: for all rod segments do
17:
PA + P
18:
end 01 H' (q:)) Pa.)/At
(Ap AA)
Dinput
X
rodConstrainndProoccion
p + p Ap
(AP",AA))
correction Guess
Decoder
GN'
GN2
GN
Figure 3: Illustration of our approach incorporating a network which consists of M graph networks (GN-blocks) into the position-bast based dynamics framework.
The most expensive part of Algorithm involves the computation of the corrections of generalized coordinates and Lagrange multipliers (line 9 This projection step requires the solution of linear system which is a linearization of non-linear one so that the matrix depends on the system's 's state making it impossible solved iteratively using the conjugated gradient (CG) solver. Such iterative solvers are widely used the context of physical simulations and regularly described as the facto standard (Barrett al.) 1994 (Sadd 2003) since they often show superior performance and usually scale well allowing for exploiting parallel hardware. However, we would like to point out that also highly efficient direct solvers can be found in the literature (Deul al. 2018).
precompute its inverse. Instead, new system every point in time is
Instead of fully replacing the projection step in an end-to-end learning manner, we follow the strategy of accelerating it by first computing a guess
AAp%. 0 A < correctionGues((p),
(1)
for the iterative procedure (line
9
(Ap, AA) < rodConstrainnPProecttonpp X "',App" AAA).
(2)
A neural network is employed to compute the initial guess in Eq. 1 for the constraint projection. The motivation for this approach is to reduce the number of iterations required for the convergence of the CG solver which solves the linear system in Eq. compared to the canonical initialization with zeros. We obtain our final framework by replacing line 9 in Algorithm | with Eq. ①I and Eq. ( as illustrated in Figure which inherently as accurate as the traditional PBD method. We name the data-driven part COPINGNer (e ("COnstraint Projection INitial Guess Network") which learns to compute the correction guess.
4

---

Under review as a conference paper at ICLR 2021
3.1 GRAPH ENCODING
COPINGNet is graph network based architecture which we apply in order to compute initial guesses for Ap and AA. In this regard, we need to incorporate the rods' state into a graph de- scription (Battaglia et al. 2018). A graph G = (V, E, U) usually consists of nodes (or vertices) V. edges E as well as global features U. However, used. For example, gravity could be a potentially meaningful global feature, but it also can be eas- ily included as an external acceleration. Hence, U = and the graph can just be represented as == GVV E). In our case, the rods' segments within the scene are represented by the graph's 'ss nodes while the interactions between the rods' segments are represented by the edges.
our framework, global features U are not
COPINGNet provides traph-to-graphe mapping: G - G from input graph Gin E G to an output graph Gout E Gout, Nodes and edges of both graphs are equipped with specific features. In the case of the input graph, the node features describe the state of the rods' segments, i.e.
Vin,, (x:; Qi Ti; Pi, la, Cis fi.: E vin C RR44
in which the positions are denoted with x; E IR", the orientations with qi E H. the radii with
rs R>0 the densities with P: R=0 and the segment lengths with l E R>00 Moreover,
segment position indicator a; E [0, 1] C R is included corresponding to a parameterization by
are length [ Binary segment flag So. [0, "left" end flag fi, {0, and "right" end flag
fz, E {0. f are set to zero if the specific segment respectively the left or right segment of the rod
fixed and one otherwise. The nodes of Gin are given as the set of Vnn = (Vin.V J. in which
= |Vin denotes the number of segments in the scene. The nodes of Gout contain the correction values of the generalized coordinates, i.e. Vout, Ap; your C R and Vout = ( Vout, While rod segments represented as node features, represent constraints between rod segments as edge features:
Gin. = (w), Yi TT) E EE' C TR$
in which the (rest) Darboux vector w E RR's describes the static angle of two rod segments, and Young's modulus Y E R>0 and torsion modulus T E R>0 are corresponding to extension, bending, and twist constraint parameters. The set of edges of the input graph is then given by Ein = L {Cim,i). in which m = |iil denotes the number of interactions between different seg- ments. The correction of the Lagrange multiplier AA; E R is stored the output edges:
Cout, Aλ, Eour C IR
The set of output edges is then given by Eout U 1 (- The connectivity information C' of each graph is stored in two vectors Cad and Cty containing the sender node index and the receiver node index of each corresponding edge. This concludes the specification of the input graph Gin = G(Vin. Ein) and the output graph Gout = G(Vout: Eout) with connectivity information C (Cad. C.v).
3.2 NETWORK STRUCTURE
In the following, we describe the structure of COPINGNet after we formalized its input and output in the previous section. As illustrated in Figure 3 the network consists of an encoder network, multiple stacks of GN-blocks and decoder network. The graph network from Battaglia al. (2018) used as a basic element and denoted as GN-block. Residual connection between blocks could improve performance of neural networks in both CNN (He et al., 2013). and graph neural network (Li et al. 2.19)]. As in related work (Sanalez-Gonnzlez et al. 2020). we employ the residual connection between the GN-blocks, but our encode//decoder network directly deals with the graph. The encoder network performs a mapping: Gn + Glatent and is implemented using two multi-layer perceptrons (MLPs), MLPedge : Enn - glatent C IR' and MLP node : Vn + ylatent C R in which denotes the latent size. They work separately and thus Een MLP MLPedge (Ei)) and Ven MLP mode (Vin). Edge features Ein from the input are constant for a rod during the simulation and this results constant encoded edge features Een, which could be recorded after the first run and used afterwards during inference. The edge feature ein.i contains the material parameters which could vary by different orders of magnitude. Hence, we normalize Young's 's modulus and torsion modulus before feeding them into the network. After encoding, the graph Gen Eee)) E G Glatent is passed to several GN-blocks with residual connections. Each GN-block also contains two MLPs. However, they use message passing taking advantage of neighbourhood nodes/edges information (Battaglia et al.)
For a single rod in the scene which consists of N segments of equal lengths, for the i-th segment, we obtain 0, (6-1)) 1)/(N- - 1) for E {1, 1,2....,n). n..
5

---

Under review as a conference paper at ICLR 2021
Train/Val
#Steps
#Nodes N
Young's Modulus Y
Initial Angle do 21(0%, 45.0%)
Rod Length é 2(10 m. 5.0 m)
256/100
50
2a(1 65)
10" Pa, a ro 24(000
Train/Val
#Steps 50
#Nodes
Torsion Modulus
Helix Radius/ / Height
Winding Number ((200,3.0)
256/100
(45, 105)
10 Pa, a (4.000 6.0) (0.4 m, 0.6 m// /24(0.4 m, 0.6 m)
Table I: Specification of training and validation datasets for the two scenarios of an initially straight bending rod (top) and an elastic helix (bottom). The datasets are comprised of a number of data points (left) each describing the rod's dynamics within t [O s. 50 discretized with a time step size of At = 0.02 s. The number of nodes N is sampled from a discrete uniform distribution Ld while the remaining parameters are sampled from continuous uniform distribution 2. We trained our network for 5 hours for the bending rod scenario and 6 hours for the helix case.
0.0
20
40
60
80
100
120
1.00
25
50
75
100
125
150
175
200
Number Nodes
Number of Nodes
Figure 4: Illustration of the ratio of COPINGNet's inference time finfer and the vanilla CG solver's run time tcg (purple curves; right vertical axis) for the initially straight bending rod (left) and the elastic helix (right) simulations. The black curves show the CG iteration number ratio while the red curves show the total speedup of the constraint projection when taking into account COPINGNet's inference time (left vertical axis). The pink curves show the total speedup of the entire simulations. The orange and green dashed line indicate the lower and upper boundaries of the total number of nodes used in the training data. We observe a speedup even for rods and helices with greater node number than the ones used in the the training dataset. The result is averaged from 50 simulations each running 100 steps.
2018). A number of M GN-blocks enable the use of neighbourhood information with distances smaller or equal to M. The graph network performs mapping within the latent space: Glatent - Glatent and after M GN-blocks, we obtain Gemn (V''. E'n)) Glatent The decoder network performs mapping: Glatent + Gout, which has similar structure as the encoder network. Two MLPs MLPadge : Iplatent + Rout and MLPpodee 3ylatent -> youut compute Foou = MMPPadee(F)) (FEn) and Vout = MLP (V'')) A tanh-function at the end of MLPedge and MLP is used restricting the output the interval [-1, C R.. COPINGNet learns the relative correction values. The generated dataset is normalized, and the maximum correction value of the generalized coordinates and the Lagrange multipliers will be recorded norms. The final correction value is achieved by multiplying the relative correction values and the norms. This normalization process damps the noise caused by the network and leads to a more stable performance. For simplicity, all the ML.Ps in different blocks have the same number of layers, and the same width as latent size within the latent layers. The input and output sizes of each MLP are consistent with the corresponding node/edge feature dimensions. The loss L is computed from both parts, nodes and edges,
L = MSE(V Vout) + MSE(Eout Eout),
in which (Vout: Eout) is the output graph's ground truth in contrast to COPINGNet's prediction (Vouts F.out) The mean squared error between x and X is denoted with MSE(X, x)
For evaluation purposes, we also incorporated the k-nearest neighbor (k-NN) algorithm as described in the supplementary material setting k- = 3.
4 EVALUATION
We generate training and validation datasets based on two scenarios: an initially straight bending rod and an clastic helix each fixed at one end oscillating under the influence of gravity. The specification of these datasets is provided in Table [ The PBD code written in C++, while the COPINGNet is implemented in PyTorch. The training is performed on an NVIDIA @ Tesla V100 GPU. The training time varies from 8 to 30 hours for different architecture parameters. A constant learning rate of 1 0.001 was used and a mean square loss function was employed. Our approach gen- eralizes across different initial conditions, geometries, discretizations, and material parameters. In the supplementary material we show that it is possible to robustly generate various dynamical re- sults (Figure 9). For a discussion on the network architecture please see Figure 11 (supplementary material).

---

Under review as a conference paper at ICLR 2021
Realistic biomechanical simulation of 3D tree model composed of over 1 nodes swaying in wind field. Our GN approach performs correctly even under a large number of rod segments while increasing performance of the original PBD method.
Figure 5:
2.5
Iterationn Ratio
2.0
2.0
1.5
1.0
200
400
600
800
1000
200 400
600
800
1000
Number
Number
Figure 6: The dashed horizontal lines show the benchmark CG solver's constant performance for an initially straight bending rod simulation (left) and an elastic helix simulation (right). The parameters du 09, N 30, 3.0 m, Y = 1.0 10 Pa (left), and HR = 0.5 m, HH 0.5 m, HW 2.5, G = 1.0 10th Pa, N = 60 (right) are applied. We are comparing our COPINGNet-asiitedd PBD approach (orange) to a simplified version in which the k-NN (green) is used to predict the initial guess. Moreover (left), we added performance measurement in which we restricted the solution to be in -wo-dimensional plane (black). The results are averages over 20 simulations and smoothed with window of 10 frames. The dashed vertical lines mark the number of 50 frames contained in the training data set. For the helix simulation we observe a speed up of approximately 50% even beyond the training data (50 time frames). For the bending rod simulation observe speedup that further decreases with increasing time frame number.
Discretization. Our approach addresses the acceleration of the most expensive part within PBD by providing an accurate initial guess of the constraint projection. We measure the system's complexity by the number of nodes in rod. Figure 4 shows the ratio of COPINGNet's inference run time (black and red eurves) compared to the run time of the vanilla CG solver (purple curves) for different numbers of nodes. The increasing black and red curves indicate that the speedup of COPINGNet is more significant with greater number nodes. With only few nodes the CG solver performs better due to the inference overhead. Once the number of nodes is increasing, a significant speedup can be obtained of up to 50% for the constraint projection. Surprisingly, our approach also performs well when going far beyond the sampling range (orange green dashed lines). Since the constraint projection is the most time-consuming part of the entire simulation, the speedup converges to that of the whole simulation (pink curves) with increasing number of nodes as shown in Figure Temporal Evolution. In addition to the complexity analysis, we also analyze the required number of CG iterations for the vanilla constraint projection compared to the one accelerated using COP- INGNet over time as shown in Figure 6 We obtain significant total speedup compared to the CG solver (dashed blue line). As stated in Table [ our training data contains dynamical simulations of 50 time steps. In this range we observe the highest speedup. As was the case for the complexity anal- ysis, we again obtain a significant speedup for simulations beyond 50 time steps. This performance gain is more present for the helix (green curve) compared the rod simulation (orange curve).
Long-term Stability. In case the constraint projection is completely replaced by COPINGNet (end- to-end approach), stability of the PBD method decreases as error accumulation takes place. This is illustrated in Figure / showing the temporal evolution of the relative change of the total rod length. An initially straight rod bending under the influence of gravity simulated using the parameters 6n = 09, N = 30, f 4.0 m, and varying Young's 's modulus Y. In this scenario the rod length is expected to stay constant during the mechanical simulation. In case COPINGNet is used in an end- to-end manner (colored lines), where the whole constraint projection step is replaced, we observe that the rod length changes incorrectly. In fact, the divergence of using COPINGNet to replace the constraint projection step increases exponentially beyond the range of the training data (50 steps). To the contrary. when COPINGNet is used to only estimate the initial guess of the constraint projection (black lines), rod deformation takes place even beyond the training data range.

---

Under review as a conference paper at ICLR 2021
Collisions. Figure I Iillustrates collision of an elastic helix with the ground plane and collision be- tween two rods. Collision detection is efficiently implemented using the hierarchical spatial hashing scheme according to Eitz & Gu (007). Rod-rod collisions are then treated with line-to-line distance constraints, and collisions with the ground plane using half-space constraints. Moreover, frictional effects are implemented according to Macklin et al. (2014). This approach allows us to handle discontinuous events such as collisions between individual rods and other objects. Moreover, for this experiment we use the neural net- works trained on the rod and helix simulations. Em- ploying COPINGNet trained on the helix data to simulate the collision of the ground plane (HR = HH = 0.5 m, HW = 2.5, G = 1.0 10 Pa, N = 50), we measure a total speedup of approx. 10%. In the case of two collid- ing rods (on = 0-, N E (20, 30), l E [4.0 m, 4.5 m), y= 1.0 10 Pa). we obtain speedup of approx. 6%. Complex Scenarios. Our method is also capable to deal with complex scenarios such as 3D tree models swaying in wind fields as shown in Figure 5 We represent trees using the extended Cosserat rod theory introduced in Pirk et al. (0177) This method allows us to simulate realistic biomechanics of thousands of rod segments at interactive rates. We generated 100 different tree topologies using 70 individual rods (average node number: 1056) and sim- ulate the swaying motions of these tree models with vanilla PBD to generate the training dataset. For the evaluation, 30 different tree topologies have been generated for each of the following four experiments using 10 rods (204 nodes), 20 rods (355 nodes). 40 rods (654 nodes), and 70 rods (1061 nodes). We were able to significantly improve the runtime performance of the method by 17.0% (10 rods), 15.8% (20 rods). 13.0% (40 rods), and 11.1% (70 rods). This takes into account the inference time introduced by the neural network.
10
Y-10 1.0 10ff Pa
Y- 2.0 10 Pb
5-0 -10f
1/0 10
5-50 10"
Y- 1.0 -10 Pb
0
20
100
Time Framm Number
Illustration of the relative change of the total rod length (l/l)) for different values of Young's moduli (Y). Colored lines show different re- sults for COPINGNet replacing the con- straint projection. The thick black line represent the result for COPINGNet re- placing only the initial guess. This indi- cates that the end-to-end approach does not generalize past the 50 time steps used in preparing the training data.
Figure 7:
Generalization. Figure 7 and Figure 10 (supplementary material) indicate that using COPINGNet in end-to-end manner does not generalize beyond the training data. Specifically, the end-to-end setup diverges in terms of rod geometry and segment position from the correct solution. This effect increases significantly beyond the training data range. Although, we only use COPINGNet benchmark for this evaluation, other GNs are expected to perform similarly. Common workarounds to increase the stability of dynamical systems with neural networks are temporal skip-connections, recurrent training, and data augmentation. However, these approaches focus on runtime speed and memory performance rather than stability (Holden et al. 2019). Interestingly, replacing just the ini- tial guess with a GN does not deform the rod or introduces discontinuities at any observed stage of the bending simulation. This means that employing GNs can result in performance increase with- out a loss of stability. Furthermore, a GN that only provides initial guesses seems to also generalize to other scenarios. We observed performance improvements for the collision and complex tree case, although the network was never trained on these different rod discretizations and topologies. This indicates that our method is capable to generalize within a specific physical scenario and to lesser extent to other scenarios.
CONCLUSION
We discovered that applying GNs for replacing the initial guess has fundamental advantages over end-to-end approaches. First, our network-enabled solver ensures long-term stability inherited from traditional solvers for physical systems, while improving runtime performance. Second, our ap- proach is able to generalize across different initial conditions, rod discretizations, and material parameters, and it handles discontinuous effects such as collisions. While end-to-end approaches offer more significant speedups, our method is superior in cases where stability is an essential re- quirement. Our approach to accelerate iterative solvers with GNs opens multiple avenues for future work. For one, it would be interesting to explore mechanical systems describing general deformable (e.g. textiles) or volumetric objects, which have been simulated with PBD. Second, our approach can be applied to other iterative methods, such in finite elements analysis or in the context of lin- ear complementarity problems (LCP). This would allow us to accelerate physical simulations, when iterative solvers applied, without compromising stability.
8

---

Under review as a conference paper at ICLR 2021
REFERENCES
Jernej Barbič and Doug I.. James. Real-time subspace integration for st. venant-kirchhoff de- formable models. ACM Transactions on Graphics, 24(3):982-99) -900 July 2005.
Richard Barrett, Michael Berry, Tony F. Chan, James Demmel, June Donato, Jack Dongarra, Victor Eijkhout, Roldan Pozo, Charles Romine, and Henk van der Vorst. Templates for the Solution of Linear Systems: Building Blocks for Iterative Methods, 2nd Edition. SIAM, Philadelphia, PA, 1994.
Peter W. Battaglia, Jessica B. Hamrick, Victor Bapst, Alvaro Sanchez-Gonzan Vinícius Flores Zambaldi, Mateusz Malinowski, Andrea Tacchetti, David Raposo, Adam Santoro, Ryan Faulkner, Çaglar Gülçehre, H. Francis Song, Andrew J. Ballard, Justin Gilmer, George E. Dahl, Ashish Vaswani, Kelsey R. Allen, Charles Nash, Victoria Langston, Chris Dyer, Nicolas Heess, Daan Wierstra, Pushmeet Kohli, Matthew Botvinick, Oriol Vinyals, Yujia Li, and Razvan Pascanu. Relational inductive biases, deep learning, and graph networks. CoRR,
2018.
URL
org/ abs/ 0ttp///886.0285000 01261
Jan Bender, Kenny Erleben, and Jeff Trinkle. Interactive Simulation of Rigid Body Dynamics in Computer Graphics. Computer Graphics Forum, 33(1):246 2014a.
Jan Bender, Matthias Müller, Miguel A. Otaduy, Matthias Teschner, and Miles Macklin. A survey on position-based simulation methods in computer graphics. Computer Graphics Forum, 33(6): 228-251, 2014b.
Jan Bender, Matthias Müller, and Miles Macklin. Position-based simulation methods in computer graphics. In EUROGRAPHICS 20/7 Tutorials. Eurographics Association, 2017.
Jens Berg and Kaj Nyström. A unified deep artificial neural network approach to partial differential equations in complex geometries. Neurocomputing, 317:28 41, 2018.
Miklós Bergou, Max Wardetzky, Stephen Robinson, Basile Audoly, and Eitan Grinspun. Discrete elastic rods. ACM Transactions on Graphics, (73))66:1-63:12, 2008.
Michael B. Chang, Tomer Ullman, Antonio Torralba, and Joshua B. Tenenbaum. A compositional object-based approach to learning physical dynamics. arXiv preprint arXiv:1612.00341 2016.
Mengyu Chu and Nils Thuerey, Data-Driven Synthesis of Smoke Flows with CNN-Based Feature Descriptors. ACM Transactions on Graphics, 36(4), July 2017.
Crispin Deul, Tassilo Kugelstadt, Marcel Weiler, and Jan Bender. Direct dusition-based solver for stiff rods. Computer Graphics Forum, 37(6):313-33 334324,,, 2018.
Gamini Dissanayake and Nhan Phan-Thien. netwol-n approximations for solving par- tial differential equations. Communications Numerical Methods in Engineering, 10(3):195- 201, 1994.
Mathias Eitz and Lixu Gu. Hierarchical spatial hashing for real-time collision detection. IEEE International Conference on Shape Modeling and Applications 2007 (SMI '07), pp. 61-70, 2007.
Stefan Feess, Kathrin Kurfiss, and Dominik L. Michels. Accurate Simulation of Wound Healing and Skin Deformation. In Proceedings of the ACM SIGGRAPt//Erroraphiis Symposium on Computer Animation, SCA '16, pp. 129-137, Goslar, DEU, 2016. Eurographics Association.
Justin Gilmer, Samuel S. Schoenholz, Patrick F. Riley, Oriol Vinyals, and George E. Dahl. Neural message passing for quantum chemistry. CoRR, Ab4/11044, 2017.
Radek Grzeszczuk, Demetri Terzopoulos, and Geoffrey Hinton. Neuroanimator: Fast neural net- work emulation and control of bhysics-based models. In Proceedings of the 25th Annual Confer- ence on Computer Graphics and Interactive Techniques, SIGGRAPH '98, pp. 9-20, New York, NY, USA, 1998. Association for Computing Machinery.
Xiaoxiao Guo, Wei Li, and Francesco Iorio. Convolutional neural networks for steady flow ap- proximation. In Proceedings of the 22Nd ACM SIGKDD International Conference on Knowledge Discovery and Data Mining, KDD '16, pp. 481-490, New York, NY, USA, 2016. ACM.

---

Under review as a conference paper at ICLR 2021
Kaiming He, Xiangyu Zhang, Shaoqing Ren, and Jian Sun. Deep residual learning for image recog- nition. CoRR, abs/1512.03338,, 2015.
Daniel Holden, Bang Chi Duong, Sayantan Datta, and Derek Nowrouzezahrai. Subspace neural physics: Fast data-driven interactive simulation. In Proceedings of the 18th Annual ACM SIG- GRAPH/Eurographicc Symposium on Computer Animation, SCA '19, New York, NY, USA, 2019.
Association for Computing Machinery.
Ruben Ibañez, Domenico Borzacchiello, Jose Vicente Aguado, Emmanuelle Abisset- Chavanne,
Elias Cueto, Pierre Ladeveze, and Francisco Chinesta. Data-driven non-linear elasticity: Con- stitutive manifold construction and problem discretization. Computational Mechanics, 60(5):
813-826, November 2017.
Byungsoo Kim, Vinicius C. Azevedo, Nils Thuerey, Theodore Kim, Markus Gross, and Barbara Solenthaler. Deep fluids: A generative network for parameterized fluid simulations. CGE, 38(2): 59-70, 2019.
Byungsoo Kim, Vinicius C. Azevedo, Markus Gross, and Barbara Solenthaler. Lagrangian neural style transfer for fluids. ACM Transaction on Graphics (SIGGRAPH), 2020.
Thomas Kipf, Ethan Fetaya, Kuan-Chieh Wang, Max Welling, and Richard Zemel. Neural relational inference for interacting systems. arXiv preprint arXiv:1802.0400 2018.
Tassilo Kugelstadt and Elmar Schoemer. Position and orientation based cosserat rods. In Pro- ceedings of the 2016 ACM SIGGRAPi/EErrgraahiis Symposium on Computer Animation. Euro- graphics Association, 2016.
L'ubor Ladicky, SoHyeon Jeong, Barbara Solenthaler, Mare Pollefeys, and Markus Gross. Data- driven fluid simulations using regression forests. ACM Transactions on Graphics, 34(6), 2015.
Isaac E. Lagaris, Aristidis Likas, and Dimitrios I. Fotiadis. Artificial neural networks for solving ordinary and partial differential equations. IEEE Transactions on Neural Networks, 9(5)9887- 1000, 1998.
Guohao Li, Matthias Müller, Ali Thabet, and Bernard Ghanem. DeepGCNs: Can GCNs Go Deep as CNNs? In The IEEE International Conference on Computer Vision (ICCV), 2019.
Yunzhu Li, Jiajun Wu, Jun-Yan Yan Zhu, Joshua B. Tenenbaum, Antonio Torralba, and Russ Tedrake. Propagation networks for model-based control under partial observation, 2018.
Miles Macklin, Matthias Müller, Nuttapong Chentanez, and Tae- Yong Kim. Unified particle physics for real-time applications. ACM Transactions on Graphics, 33(4), July 2014.
Miles Macklin, Matthias Müller, and Nuttapong Chentanez. Xpbd: Position-based simulation of compliant constrained dynamies. In Proceedings the 9th International Conference Motion in Games, MIG 11,, pp. 49-54, New York, NY, USA, 2016. ACM.
Dominik L. Michels and Mathieu Desbrun. A Semi-analytical Approach to Molecular Dynamics. Journal of Computational Physics, 303:336 354, 2015.
Dominik L. Michels, Gerrit A. Sobottka, and Andreas G. Weber. Exponential Integrators for Stiff Elastodynamic Problems. ACM Transactions on Graphics, 33(1):7:1-7:20, 2014.
Dominik L. Michels, J. Paul T. Mueller, and Gerrit A. Sobottka. Physically Based Approach to the Accurate Simulation of Stiff Fibers and Stiff Fiber Meshes. Comput. Graph., 53(PB): December 2015.
Dominik L. Michels, Vu Thai Luan, and Mayya Tokman. A Stiffly Accurate Integrator for Elasto- dynamic Problems. ACM Transactions on Graphics, 36(4), July 2017.
Siddhartha Mishra. A machine learning framework for data driven acceleration of computations of differential equations, 2018.
Damian Mrowca, Chengxu Zhuang, Elias Wang, Nick Haber, Li Fei-Fei. Joshua B. Tenenbaum, and Daniel I. K. Yamins. Flexible neural representation for physics prediction, 2018.
10

---

Under review as a conference paper at ICLR 2021
Matthias Müller, Bruno Heidelberger, Marcus Hennix, and John Ratcliff. Position based dynamics. Journal of Visual Communication and Image Representation, 18(2):109 118, 2007.
Andrew Nealen, Matthias Müller, Richard Keiser, Eddy Boxerman, and Mark Carlson. Physically based deformable models in computer graphies. Computer Gruphics Forum, 25(4):809-836, 2006.
Dinesh K. Pai. Strands: Interactive simulation of thin solids using cosserat models. Computer Graphics Forum, (1(3):347-352 2002.
Sören Pirk, Michal Jarzabek, Torsten Hädrich, Dominik L. Michels, and Wojeiech Palubicki. In- teractive wood combustion for botanical tree models. ACM Transactions on Graphics, 36(6):
197:1-197:12, November 2017.
Maziar Raissi and George Em Karniadakis. Hidden physics models: Machine learning of nonlinear partial differential equations. Journal of Computational Physics, 357:125 141, 2018.
Maziar Raissi, Paris Perdikaris, and George Em Karniadakis. Physics-informed neural networks: A deep learning framework for solving forward and inverse problems involving nonlinear partial differential equations. Journal of Computational Physics, 378:686 2019.
Yousef Saad. Iterative Methods for Sparse Linear Systems. Society for Industrial and Applied Mathematics, USA, 2nd edition, 2003.
Alvaro Sanchez-Gonzalez, Nicolas Heess, Jost Tobias Springenberg, Josh Merel, Martin Riedmiller, Raia Hadsell, and Peter Battaglia. Graph networks as learnable physics engines for inference and control, 2018.
Alvaro Sanchez-Gonzalez, Victor Bapst, Kyle Cranmer, and Peter Battaglia. Hamiltonian graph networks with ode integrators, 2019.
Alvaro Sanzhez-Gonzzlezz, Jonathan Godwin, Tobias Pfaff, Rex Ying, Jure Leskovec, and Peter W. Battaglia. Learning to simulate complex physics with graph networks, 2020.
Franco Scarselli, Marco Gori, Ah Chung Tsoi, Markus Hagenbuchner, and Gabriele Monfardini. The graph neural network model. IEEE Transactions on Neural Networks, 20(1):61-80, Jan 2009.
Sungyong Seo and Yan Liu. abs/1902.02250, 2019.
Differentiable physics-informed graph networks.
CoRR,
Justin Sirignano and Konstantinos Spiliopoulos. Dgm: deep learning algorithm for solving partial differential equations. Journal of Computational Physics, 375:1339 1364, 2018
Jonathan Tompson, Kristofer Schlachter, Pablo Sprechmann, and Ken Perlin. Accelerating eulerian fluid simulation with convolutional networks. CORR, abs/1607.00 0005970 2016.
Rahul Trivedi, Logan Su, Jesse Lu, Martin F Schubert, and Jelena Vuckovic. Data-driven accelera- tion of photonic simulations, 2019.
Kiwon Um, Xiangyu Hu, and Nils Thuerey. Liquid splash modeling with neural networks, 2017.
Nobuyuki Umetani, Ryan Schmidt, and Jos Stam. Position-based elastic rods. In Proceedings of the ACM SICCRAPH/Eurographics Symposium on Computer Animation, SCA '14, pp. 21-30, Pire-la-Ville Switzerland, Switzerland, 2014. Eurographics Association.
Benjamin Ummenhofer, Lukas Prantl, Nils Thuerey, and Vladlen Koltun. Lagrangian fluid simu- lation with continuous convolutions. In International Conference Learning Representations, 2020.
Huamin Wang. James F. O'Brien, and Ravi Ramamoorthi. Data-driven elastic models for cloth: Modeling and measurement. ACM SIGGRAPH 2011 Papers, SIGGRAPH 11, pp. 71:1-71:12, New York, NY, USA, 201 ACM.

---

Under review as a conference paper at ICLR 2021
Xiaolong Wang, Ross B. Girshick, Abhinav Gupta, and Kaiming He. Non-local neural networks. CoRR, abs/1711.07911, 2017.
Steffen Wiewel, Moritz Becher, and Nils Thuerey. Latent-space physics: Towards learning the temporal evolution of fluid flow. CoRR, 1tss/880..1123, 2018.
Xiangyun Xiao, Yanqing Zhou, Hui Wang, and Xubo Yang- A novel cnn-based poisson solver for fluid simulation. IEEE Transactions on Visualization and Computer Graphics, (663):14441444 2020.
12

---

Under review as a conference paper at ICLR 2021
SUPPLEMENTARY MATERIAL
DYNAMIC RESULTS
t= 0.0s
t- 1.0s
tm 2.05
t= 3.0s
t- 4.05
-2-
-2
-6
-6
-6
-6
-6
0.0
2.5
5.0
0.0
2,5
5.0
0.0
2.5
5.0
0.0
2.5
5.0
0.0
2.5
5.0
b)
0
0-
0
-4
-4
6
2.5
5.0
0.0
2.5
0.0
5.0
0.0
2.5
5.0
0.0
2.5
6)
0-
0-
4
4-
44
44
4
d)
-2-
-2
-4
-4
-4
5,0
2.5
-1
-2
-2
M
3
M
-1
7
-2
2
2
VMM
-2
4
-3
33
0
0
Figure 8: Illustration of the temporal evolution of scenarios (a) to (d) for the initially straight bending rod, and (e) to (f) for the elastic helix. The different cases represent various material property configurations.
13

---

Under review as a conference paper at ICLR 2021
In Figure show the temporal evolution of different scenarios (a) (d) for the initially straight bending rod, and (e) to (f) for the elastic helix. The default parameters of the ini- tially straight bending rod are on 09, .N- 30, L = 4.0 m. and Y 1.0 10y Pa. In In (b), we modify N E {10, 20, 40, 60). In 6.0m). In (d), we modify Y (2.0 10 Pa, 5.0 10 Pa, 1.0 - 100 Pa, 2.0 -100 Pa. 5.0 10 Pa, 1.0 -10 100 Pa). The default parameters 0.5 m (helix radius). HH 0.5 m (helix height), HW = 2.5
(a), we modify Φn E {0.0", 15.09, 13.00, 30.0%,45.0%)
(c), modify L (1.0 2.0 m.0.
of the elastic helix are HR (winding number), T f(0.4 m, 0.4 m, m). (0.4m. m. 0.4 1..4m,3.0mm 3)0 (0.5)).0.5mm m. 0.5 m.0m,..0. 2) (0.5 m, 0.5 m. 3.0m)).
1.0 - -00 Pa, and N 60. In (e). we modify (HR, HH, HW) E
In (f),
we modify N E {30, 60, 80, 1003. In (g), we modify G E (2.0-10* Pa, 5.0- 10 Pa, 1.0- 10 Pa, 1.0- 10 Pa]. For each experiment, the rod colors indicate the corresponding parameters in the following (ascending) order: blue, orange, green, red, purple, brown.
Dynamic results of our graph learning assisted simulations can be found in the video
GrapilLaninnghhssialyyttmemmmm (2:20 min, 20.6 MB, 1920 1080, H.264)
provided as additional supplementary material. These results are obtained by employing three GN blocks, two MLP layers, and a MLP width of 32. Parameters have been used which are are not covered within the training dataset.
Moreover, to further demonstrate general applicability, Figure | shows an additional simulation of tightening knot for which COPINGNet-based PBD shows total speedup of 10.4% compared to vanilla PBD. COPINGNet trained using a dataset containing 20 simulations of the same knot scenario with different discretizations over 100 frames. It can clearly be observed that COPINGNet generalizes beyond the training data.
Framedas
fouretos
Frmm e
Fowee 423
Frum 227
Frmends
Frame:
hore 64B
hram t55
Furte 962
Frems 268
Prame
Frame 082
Huree (96
Nave t34
Panee 103
Fromc 138
Pranel
Fame
Ionre1
Faym 137
Faaee 161
Foom
Frmme Wt
17
Pamee 174
Faree
Fanm
Frame230
Figure 9: Illustration of the temporal evolution of a knot scenario. The rod is fixed both ends and the knot is pulled tight. The predicted results by end-to-end COPINGNet-based learning are shown in red while the results computed with COPINGNet-assitedd PBD are shown in black. The visualization of the last frame contains close-up to demonstrate that the knot structure is still preserved. The rod parameters are N 110, = 0.01 m, l 11.0 m. and Y = 1.0 10 Pa.
14

---

COMPARISONS
Figure 10 illustrates the temporal evolution of a bending rod and elastic helix scenarios using dif-
08, N 30, 4.0 m, and Y = 1.0 109 Pa are used. the helix case, HR = 0.5 m (helix radius), HH = 0.5 m (helix height),
ferent approaches. For the bending rod case, the parameters
HW 2.5 (winding number), T 1.0 10 Pa, and N 60 were applied. The temporal evolution of the positions' normalized root mean square error (NRMSE) is shown in Figure
-03
-05 6.3
25 -65 3.0
0.5
-25 0.6 -05 60
15
-6.5 6.0
Figure 10: Illustration of the temporal evolution of bending rod (upper row) and elastic helix (lower row) scenarios simulated using our COPINGNat-asisted PBD approach (blue curves), vanilla PBD (orange dotted curves), k-NN-based end-to-end learning (purple curves) and GN-based end-to-end learning (green curves). While the cosiNGNet-asistedd approach and the vanilla PBD solver allow to simulate both cases in stable manner, both end-to-end learning approaches diverge with increasing time due to error accumulation.
Pesition NRMSE
1.5
Position NRASE
1.0
0.5
100
60
Time Frame Number
Time Frame Number
Figure 11: Illustration of the positions' normalized root mean square error (NRMSE) compared vanilla PBD using our COPINGNet-assiseed PBD approach (blue curves), k-NN-based end-to- end learning (purple curves), and GN-based end-to-end learning (green curves) for the bending rod (left) and the elastic helix (right) scenarios. While the COPINGNet-assised PBD approach converges almost correct solution, both end-to-end learning approaches show large error rates with increasing time frame number.
ARCHITECTURE
The evaluation of COPINGNet's architecture, introduced in Section 3.2 is presented in Figure 12 The performance studied for different numbers of GN-blocks and MLP layers, and different MLP widths (latent sizes). The measurements demonstrate that larger number of GN-blocks usually increases the performance while the performance improvement is not longer significant for more than four GN-blocks. This is plausible and can be considered analogously to the size of the stencil of a numerical integrator: a larger number of GN-blocks means that a specific node can take advantage of information gained from neighbourhoods further away. Correspondingly. larger stencil can do

---

Under review as a conference paper at ICLR 2021
well. However, increasing the size of the stencil usually does not result in more accurate results once the critical size is reached. This can be observed here as well. contrast, it can be clearly observed that deeper MLP networks do not improve the performance. This is consistent with other research on graph networks (Sanchez-Gionzalez et al.) 2020). However, increasing the MLP width can increase the network's ability generalize. It shown that the highest performance is achieved with the largest MLP width, especially in the case of the helix scenarios. Since increasing the number of GN-blocks and MLP width will lead to longer inference time, we compromise by choosing medium numbers. For further evaluations, we employ three GN-blocks, two MLP layers, and MLP width of 32.
1.6
1.0
16
Number of blocks
MLP Layers
MLP wioth
collwabee Raiio
1.6
1.4
1.2
1.0
Number Gh-brocks
wheth
Figure 12:
Illustration of the evaluation of COPINGNet's architecture (blue: benchmark architec- ture). The result is averaged from 50 test simulations each running for 100 time steps.
k-NEAREST NEIGHBOR ALGORITHM
For evaluation purposes, we also incorporated the k-nearest neighbor (k-NN) algorithm into our framework compute an initial guess for the constraint projection step.
Consider k-NN data point yi). Its components contain node and edge feature information:
X; (Vin,i, Cin, )T, Y; (Vout, Cout,i .T k Nu-
Simulation data is collected to form dataset ((x, y.) During inference, the nearest &s data points of x are picked and linear interpolations are performed obtain the predicted y
For each input node feature vin, obtain the corresponding output node feature Vout, Vout.i- The output edge feature is averaged from the information of the corresponding component y):
1 M
Σ oott.ja
out
in which M denotes the number of nodes.
The k--N algorithm is applied as an end-to-end method completely replacing constraint projection or to predict initial guess.
16


Reviews:
[{"review_id": "vxN_YJV9d5", "reviewer": "AnonReviewer2", "summary": "", "questions": "", "limitations": "", "rating": 6, "confidence": 4, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "main_review": "", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "", "review_text": "Summary: the authors present a Graph Network (GN) architecture to speed up the running time of rod physical simulations by predicting the initial guesses to reduce the number of iterations of an iterative position-based solver. The proposed approach offers speed-ups varying from 6-50%, depending on the problem analyzed. \n\nStrengths: \nThe paper is well written and motivated, methods are very clear. There are substantial experiments to test for the paper's main claims on inference efficiency and out of distribution generalization, although mostly on small scale systems (up to 200 nodes) and short time horizons (100 steps).\n\nWeaknesses: \n\n1) The biggest weakness to me is shown on Figure 6, when plotting the CG interaction ratio as a function of time, for the straight bending rod and the elastic helix, up to 1000 timesteps. Especially on the straight bending rod simulation (orange curve), the speed-up is not consistent across time. Further insight on why those oscillations are happening would be appreciated. In the elastic helix simulation (green), one also sees a huge drop right after the training regime range (50-100 steps). End-to-end approaches will generally offer more significant speed-ups (orders of magnitudes rather than the 6-50% provided by COPINGNet), but they might struggle to remain stable for longer time horizons. Therefore what one would like to see from an alternative approach like the proposed here, is a system that scales well with time. Figure 6 shows that's not the case here yet. \n\n2) It would have been interesting to see how the GN compares to a different neural network predicting the initial guesses (say a simple LSTM or -- more ambitiously -- a transformer). I understand that many of the successful neural approaches for physics simulations are GN based, so the choice is well motivated -- but those systems are generally trained end-to-end, which isn't the case in the work from this paper. Therefore, it's not clear to me that other non GN approaches wouldn't perform well on this task of producing better initial guesses for an external solver.\n\n3) In figure 4, even if small, you can see a decay on the speed-ups as you have a greater number of nodes (towards the right tail of the curve). Since the main interest in an approach that purely offers speed-ups (rather than from example, additional better generalization or accuracy) would be on large scale simulations, I see this as a fundamental weakness.\n\n4) When comparing to end-to-end approaches, authors cite a few times in the text instability for longer time sequence rollouts. Note that in Sanchez et al. 2020, the authors propose adding noise to the inputs and correcting the predictions as a stabilizing technique, they present stable rollouts for systems with thousands of steps (longer than the ones considered proposed here). This could be integrated in the end-to-end approach tested on the paper. \n\n\nAll in all, despite its weaknesses, I still think the paper takes a meaningful step towards the efficient use of neural networks for physical simulations. So I am willing to adjust my score after some of my concerns have been addressed. \n\nSmall typos:\nIn the \"Generalization\" paragraph from the Evaluation section, should it be Figures 9 and 10, rather than 7 and 10?", "extra_scores": {}, "extra_text": {}}, {"review_id": "p_8vbOFe4Gm", "reviewer": "AnonReviewer3", "summary": "", "questions": "", "limitations": "", "rating": 4, "confidence": 4, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "main_review": "", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "", "review_text": "\nUsing the same bar as NeurIPS, I continue to recommend rejecting this paper.\nSince the paper remained largely the same. My review remains largely the same.\nIt is also doubtful that these changes have satisfyingly address the other three\nreviews.\n\nThis paper proposes an iterative PBD solver which uses a neural network to guess\nthe constraint force and position updates before polishing with a conjugate\ngradient solver. The method utilizes a graph network architecture which makes it\nagnostic to the particular discretization allowing generalization (in theory)\nacross scenarios.\n\nSpeeding up simulations is an evergreen topic. A strength of this work is the\nproblem that it is tackling. The choice of methodology makes it mildly of\ninterest to the ICLR community. It is an application (without adaptation) of\ngraph networks. The result is mildly positive (though not earth shattering),\nindicating success of graph networks to some degree. I have not seen graph\nnetworks applied to rod simulation or elasticity simulation, although they have\nbeen applied to meshes in geometry processing (similar and more challenging\nscenario) and neural networks have been applied more generally to speeding up\nelasticity simulations (with what would be a trivial extension to rods). A\npossibly unique strength of this paper is the combination of rod simulation and\na convergent solver (i.e., just using the network as an initial guess in an\niterative solve). However, this is a really straight forward idea and the gains\nare again small. So, while this is a strength it is minor and possibly not\nexciting to the broader ML/ICLR community.\n\nMy main criticism of this method have to do with three aspects of scaling: 1)\nthe network is applied at the inner most loop and message passing occurs by\niterating over edges (which I believe are based on the original connections).\nThis means that global information passing requires M = O(n) iterations. 2) the\nmethod operates only on the fine scale details (the input resolution). This runs\ncounter to multigrid literature which would suggest operating on all the\nresidual at all frequency levels. 3) the method shows a small number of small\nexamples (low resolution rods with one or two rods in a scene). PBD is often\nemployed in scenarios with millions of rod segments from thousands of rods\n(e.g., hair on a human head will have 100,000 rods each with hundreds/thousands\nof segments). For a 1.5x to become impressive I would like to see this operating\nat scale.\n\nI worry that the network has simply memorized a mapping from global positions to\nguesses. Since the network has access to the raw positions, I would be surprised\nif it has learned rotation and translation invariance. Lack of this would be a\ngood indication that it's learning a spatial mapping. This would be a severe\nlimitation.\n\nThe exposition of main results is very confusing. If I understand this part of\nFigure 5 (and I doubt I do), then the pink curve below the red and black curves\nis indicating that using this method does not simply speed up the CG solve but\nsome how negatively affects the total runtime (I guess by confusing the outer\nloop down an incorrect path). So the gains by taking an aggressive initial guess\nare tempered a bit, though still resulting in a (quite modest) overall speed up\n<1.5x.\n\nThis modest speedup coupled with the small number of simple experimental\nscenarios worries me that these results will not generalize to a statistically\nsignificant speed up in general. \n\nCollisions appear as an afterthought. If the speedups were more significant, I\nwould happily excuse this as collisions can often be an extra systems effort.\nHowever, part of the \"glory\" of a PBD solver rather than an FEM+LCP type\nsimulation is that one can throw all sorts of constraints into the same system.\nSo, when this paper tacks on collisions outside the learning aspect of this work\nit calls the choice of PBD into question as perhaps needlessly inaccurate or an\n\"purely-out-of-convenience\" type choice rather than a scientifically motivated\none.\n\nIf I understand correctly, the paper always compares to a baseline of resetting\nthe position and multiplier updates to zero. Is this the best baseline? What\nabout using the previous iteration? Or other momentum based strategies? \n\nThe for-loop on line 8 of the pseudocode is misleading/confusing/incorrect. This\nimplies that updates for rods are conducted independently. But then line 9\nappears to accept as input and send as output coordinates and lagrange\nmultipliers for the entire system (rather than the rod i). What is the role of\nthe for loop of line 8 if line 9 does not depend on which rod is being\nconsidered? it would be better to write out the linear solve that is being\nconducted (e.g., with CG). The figure and text below confuses me further. Is the\ncorrection guess applied: a) once per time step, b) once per constraint\nlineariztaion (just before CG), or c) once per rod?\n\nThe revised paper (and perhaps rebuttal) should include a clarified pseudocode\nto understand what this method is actually proposing.\n\nThe effectiveness of the network depends on a parameter M which would appear to\nscale poorly with the resolution of the input shapes. Does M need to be adjusted\nfor higher resolution examples?\n\n\nWhat is the video showing? Training data created using the groundtruth\nsimulation? Results of this method at test time? If so what was the training\ndata used for each? This video did not really help supplement this submission. I\nwish that the video had included experiments that could be used to gain\nintuition about what's being learned and about the accuracy of the initial\nguess.\n\nFor example, the paper mentions many times that it is *not* replacing the time\nintegration/constraint projection with an end-to-end trained network for\nrobustness reasons. Let's see it fail then!\n\nIn a future revision I would like to understand\n  - under what circumstances does end-to-end learning fail, but this method\n    succeeds.\n  - under what circumstances does simply using the previous frame's constraints\n    as an initial guess out perform this method?\n  \nIt'd be great to see a video where we also see an evolving plot per-frame (like\nFigure 6) showing the performance gains of applying this method rather than\nnaive initialization methods. This would be great to get match up whether the\ngains happen far from the rest state, in collision heavy scenarios, or the\nopposite etc.\n\nOn line 180, I did not understand how the set of input edges to the graph\nnetwork is defined. If inputs nodes are assigned to each rod segment, then input\nedges should connect two segments. Are only neighboring segments connected with\nthese edges? (e.g., the dual graph of the polyline representing the central axis\nof the rod). Or are all possible pairs of segments generated? Where are Young's\nand Torsion moduli stored? It's natural to store Young's modulus at segments but\nthis would correspond to nodes of the graph not edges.  \n\nThe graph network description (which gets at the core contribution of this\npaper) is poorly described. I read this section many times and via cross\nreferencing with [13] could finally understand with low confidence what is being\ndone.\n\nFigure 4 is very confusing. Is this figure showing that none of the hyper\nparameters matter except for MLP width? This is surprising to the point of\nindicating a bug/overfitting.\n\nI do not understand Figure 5. What is the purple curve? the vanilla CG solver?\nThen wouldn't its speed up be 1x? Are these plots showing two different y-axis\nor two different examples (straight vs helix)? The caption is very confusing.\n\nThe paper does not accurately categorize past works when it writes \"Existing\nmethods enable learning these systems often in an end-to-end manner and with a\nfocus on replacing the entire integration procedure.\" Many fluids papers retain\npressure projection to ensure divergence free-ness and within elasticity\nsimulation, \n\nLatent-space Dynamics for Reduced Deformable Simulation\nLawson Fulton, Vismay Modi, David Duvenaud, David I.W. Levin, Alec Jacobson\n\ndoes not replace the time integration procedure.\n\nIn the future, I would appreciate a pdf in supplemental material with\nhighlighted changes. It is potentially rude to continuous reviewers to have to\nhunt for small changes (and then see that their reviews were largely ignored).\n", "extra_scores": {}, "extra_text": {}}, {"review_id": "SUjYe5z7fuz", "reviewer": "AnonReviewer4", "summary": "", "questions": "", "limitations": "", "rating": 6, "confidence": 4, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "main_review": "", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "", "review_text": "The paper proposes an algorithm to accelerate simulations of deformable rods based on position-based dynamics. Despite what the title and abstract suggests, the scope of the paper is limited to a single physical system, and there is no evidence that this approach will work on generic physical systems. The title and abstract should be tuned down and made more specific.\n\nThe key idea is novel: instead of replacing the entire time integrator with a neural network, which is typical of previous approaches, the authors propose to use the network to accelerate the constraint project step, and in particular to still rely on the standard projection used in PDB, but using the network to generate the initial guess for the nonlinear optimization. A better initial guess will reduce the number of iterations (an ideal one will cause termination after 1 iteration), thus reducing runtime. Errors in the prediction will still likely be corrected by the non-linear optimization making the approach stable. \n\nWhile I think the idea is great on paper, I am concerned by the results and choices made in the paper, in particular:\n1. The abstract makes claims of generality, while the approach is very specific to rod simulation with a very specific solver.\n2. The approach is based on Daul et al. 2018, it is likely that the authors started from their codebase. However the key contribution of that paper is “However, this solver requires many iterations to converge for complex models and if convergence is not reached, the material becomes too soft. In contrast we use Newton iterations in combination with our direct solver to solve the non-linear equations which significantly improves convergence by solving all constraints of an acyclic structure (tree), simultaneously. Our solver only requires a few Newton iterations to achieve high stiffness and inextensibility”. My understanding is that the authors of this submission are purposely comparing against the CG solver instead of the new faster solver proposed in Daul et al. 2018. Is this correct and if so why? The comparison should be done with the state of the art, which is not CG.\n3. I am having a hard time understanding the plot in figure 4. Is this the entire runtime or just the inference time for the projection part? It seems that the speedup is ~20%  which is minor and I believe switching to newton as suggested in Daul et al would likely give you way more. Please report complete runtimes of the entire simulation, and also a plot of the errors. A speedup of 20% only on a specific step of the algorithm is in my opinion not worth training a neural network and thus losing generality. Especially since this is applied to an algorithm that is not the current state of the art.\n4. How is CG initialized? Do you use the previous solution?\n5. Please add a comparison of your network against a simple nearest neighbor search over the training data, as a simple baseline.\n6. Report complete timings of how long it takes to train. How many simulations do you need to run before the training time is amortized by the gained acceleration?\n7. Only the simpler examples from Daul et al 2018 are shown in the paper, the rope knots and the teaser model have been omitted. Why? These are the more challenging cases where errors and penetrations would be immediately visible.\n\nWhile I agree with the author's conclusions that “We discovered that applying GNs for replacing the initial guess has fundamental advantages over end-to-end approaches.” this is by itself not a surprising result. Endtoend approaches are, currently, usually worse than traditional time integrators for integrating physical systems. Instead of replacing the entire system, replacing a smaller piece does less damage than replacing the entire integrator, but it is still likely worse if the comparison is properly done with the state of the art. I don’t think these results are publishable until a clear advantage over the state of the art is shown, or an argument that is not superior performance is made against standard PBD simulators.\n\nA minor comment: “is included corresponding to a parameterization by arc length” I do not understand this sentence.\n\nUpdate:\n\nThe paper has considerably improved:\n\nthe title is now accurately describing the paper\ncomplex scenes have been added (and the method works there too)\ncomparison with knn has been added and it shows a clear improvement for the proposed method\n\nI think the paper has improved, but I still find the comparison with direct solvers problematic. For small scenes like the one shown in this paper, a direct solver is fine, there is no need to use an iterative one. If the scenes are large enough to require an iterative solver (i.e. a direct one runs out of memory) then it should be shown that the proposed methods provide benefits in that specific setting. It could be that my bar for comparison is too high, as I usually publish in a different community where quantitative improvement against the closest baseline is always required.\n\nOverall, I am still mildly positive, but not willing to champion this paper given the many issues raised in the reviews.\n", "extra_scores": {}, "extra_text": {}}, {"review_id": "Jqrl6mkIb0c", "reviewer": "AnonReviewer1", "summary": "", "questions": "", "limitations": "", "rating": 6, "confidence": 4, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "main_review": "", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "", "review_text": "This paper proposes a graph network(GN) called COPINGNet (“COnstraint Projection INitial Guess Network”), which learns to compute a good initialization for the traditional PBD method. To simulate a physical system, PBD first computes updated locations of vertices then corrects the estimates of the initial position by constraint projection.  The projection step is computationally expensive, and that is where the proposed COPINGNet is applied to generate a good initial guess for the built-in linear system solver (e.g., CG).\n\nStrengths:  \n1. The proposed method is superior to end-to-end approaches in terms of long-term stability and out-of-distribution transfer (initial conditions, discretizations, material parameters, etc.).\n2. This paper applies the idea of message passing to design the graph network architecture, which enables propagations of node information (location, speed, etc.) and edge information (Young's module, torsion module, etc.) on the graph.\n\nWeaknesses: \n1.  According to Figure 4 and Figure 6,  the total speedup of the entire simulations is approximately 1.4 compared to vanilla CG solver.  This improvement seems marginal.\n2. It seems that the current model cannot generalize to different shapes. If a network needs to be retrained for each different shape, it may significantly limit its applicability in practice.\n\nQuestions:\n1. Can the initialization guess network accelerate other iterative solvers? \n2. According to Figure 6, helix simulation and bending rod simulation show different speedup rates (when extrapolating in time). Can this be explained? \n\nGenerally speaking, this article provides a new idea. By combining graph neural networks with traditional methods, the speed of solving physical systems is improved. Due to some concerns about the speedup ratio and generalization, the practicability of this method needs further investigation.", "extra_scores": {}, "extra_text": {}}]

Author rebuttal:

Meta-review: The paper proposes speeding up iterative simulations of complex dynamics systems based on connected rods. Traditionally, these systems alternate between forward integration of the dynamics and constraint projections. Instead of replacing this entirely with end-to-end trained ML, here ML is only added a single point in the method to speed up the iterative solver itself, more precisely by providing initial estimates for the constraint projection step. This is done with graph networks.

At initial evaluation, the paper had two slightly favorable reviews (6) and two unfavorable reviews (4) and was therefore on the fence leaning towards rejection.

Reviewers appreciated a well motivated method and in an interesting problem.

However, on the downside, issues raised where lukewarm performance; novelty (a direct application of graph networks); lack of generality of the approach; similarity to graph networks applied on mesh based physics simulations, and similarity to NN applied for speeding up elasticity simulations; application on the finest level only; memorization/lack of generalization; simplicity of baselines; simplicity of tasks (including the added more complex tree task).

There seemed to be some confusion on whether one of the reviewers had read the initial NeurIPS submission only (which he also had reviewed) or also the ICLR submission; the authors seemed to be upset up this possibility and made it clear in their responses, but the AC can assure them that the proper version has been read, reviewed and discussed; the author's responses in that respect were not helpful.

The authors provided responses to most of the raised issues, and several reviewers acknowledged that the paper had been improved, in particular by adding comparisons (e.g NN search).

However, in spite of these improvements, the discussion among reviewers and AC revealed that the paper still has serious issues, in particular minor novelty, lukewarm improvements, and some issues re: comparisons to baselines. While the reviewers acknowledged merits in the idea, the weaknesses hindered them to champion the paper for acceptance at this point, and the AC concurs, recommending rejection.
Decision: Reject