Accurately Solving Physical Systems with Graph Learning
Han Shao, Tassilo Kugelstadt, Torsten Hädrich, Wojciech Pałubicki, Jan Bender, Sören Pirk, Dominik L. Michels
AAccurately Solving Physical Systemswith Graph Learning
Han Shao
KAUST [email protected]
Tassilo Kugelstadt
RWTH Aachen [email protected]
Wojciech Pałubicki
Jan Bender
RWTH Aachen [email protected]
S¨oren Pirk
Google Brain [email protected]
Dominik L. Michels
KAUST [email protected]
Abstract
Iterative solvers are widely used to accurately simulate physical systems. Thesesolvers require initial guesses to generate a sequence of improving approximatesolutions. In this contribution, we introduce a novel method to accelerate iter-ative solvers for physical systems with graph networks (GNs) by predicting theinitial guesses to reduce the number of iterations. Unlike existing methods thataim to learn physical systems in an end-to-end manner, our approach guaranteeslong-term stability and therefore leads to more accurate solutions. Furthermore,our method improves the run time performance of traditional iterative solvers. Toexplore our method we make use of position-based dynamics (PBD) as a com-mon solver for physical systems and evaluate it by simulating the dynamics ofelastic rods. Our approach is able to generalize across different initial conditions,discretizations, and realistic material properties. Finally, we demonstrate that ourmethod also performs well when taking discontinuous effects into account suchas collisions between individual rods.
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 thecontext of physics-based systems, the mathematical model is usually based on first principles anddepending on the properties of the simulated system, the numerical integration of a complex sys-tem can be very resource demanding [37], e.g., hindering interactive applications. Enabled by thesuccess of deep neural networks to serve as effective function approximators, researchers recentlystarted investigating the applicability of neural networks for simulating dynamic systems. Whilemany physical phenomena can well be described within fixed spatial domains (e.g., in fluid dynam-ics) that can be learned with convolutional neural network (CNN) architectures [9, 17, 49, 57], alarge range of physical systems can more naturally be represented as graphs. Examples include sys-tems based on connected particles [36], coupled oscillators [30, 33], or finite elements [37]. Existingmethods enable learning these systems often in an end-to-end manner and with a focus on replac-ing the entire integration procedure. A number of methods show initial success in approximatingphysical systems; however, they often fail to reliably simulate the state of a system over longer timehorizons if significant disadvantages are not accepted, such as the use of large data sets containinglong-term simulations and the employment of specific memory structures [44].In this paper, we aim to improve the performance of iterative solvers for physical systems with graphnetworks (GN). An iterative solver requires an initial guess, and based on it generates a sequenceof improving approximate solutions. The initial guess can be computed by simply using valuesobtained in the previous iteration or by solving a few steps of a similar but simpler physical system.The performance of an iterative solver is significantly influenced by the calculation of the initialguess, which we aim to replace with the prediction of a GN. To demonstrate our approach, weuse a position-based dynamics (PBD) solver that approximates physical phenomena by using setsof connected vertices [4, 5, 28, 36]. To simulate a physical system, PBD first computes updatedlocations of vertices using symplectic Euler integration to then correct the initial position estimatesso as to satisfy a set of predefined constraints. The correction step is known as constraint projection a r X i v : . [ phy s i c s . c o m p - ph ] J un igure 1: Renderings taken from real-time simulations of the elastic deformation of a helix fallingdown on the ground plane (left) and two rods colliding with each other (right).and is commonly solved iteratively. The explicit forward integration for predicting the system’supdated state has negligible cost, whereas the projection step is computationally expensive. Ourgoal 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 a classic PBD solver, while providingbetter run-time performance.To showcase the capabilities of our combined PBD solver, we aim to simulate the physically plau-sible mechanics of elastic rods. Rods play an important role for a variety of application domains,ranging from surgical simulation of sutures [14], catheters, and tendons [38], to human hair [32]and vegetation [39] in animated movies. Furthermore, approaches exist to realistically simulate rodsas sets of connected vertices accurately capturing their mechanical properties [7, 23, 32, 38]. Ourapproach is able to generalize across different initial conditions, rod discretizations, and realisticmaterial parameters such as Young’s modulus and torsional modulus [10]. Moreover, we demon-strate that our approach can handle discontinuous collisions between individual rods. Figure 1 showsexamples of elastic rod deformations of a helix falling down (left) and two colliding rods (right). Fi-nally, we show that the data-driven prediction of the initial guesses of the constraint projection leadsto a decreased number of required iterations, which – in turn – results in a significant increase ofperformance compared to canonical initial guesses.In summary, our contributions are: (1) we show how to accelerate iterative solvers with GNs; (2) weshow that our network-enabled solver ensures long-term stability required for simulating physicalsystems; (3) we showcase the effectiveness of our method by realistically simulating elastic rods;(4) we demonstrate accuracy and generalizability of our approach by simulating different scenariosand various mechanical properties of rods including collisions. In the following we provide an overview of the related work, which spans from data-driven physicssimulations and graph learning to position-based dynamics and elastic rods.
Data-driven Physics Simulations.
It has been recognized that neural networks can be used as ef-fective function approximators for physical and dynamic systems. To this end, early approachesfocus on emulating the dynamics of physics through learned controllers [16] or by designing sub-space integrators [1]. Today, a range of approaches exist that enable learning ordinary and partialdifferential equations [25, 40, 41], for example, to transform them into optimization problems [11],to accelerate their computation [34, 48], or to solve for advection and diffusion in complex ge-ometries [6]. Other methods focus on specific data-driven solutions for non-linear elasticity [19],for approximating Maxwell’s equation in photonic simulations [50], or for animating cloth [54].More recently, research on data-driven approaches for modeling the intricacies of fluid dynamicshas gained momentum [24, 53]. Due to fixed-size spatial representation of Eulerian fluid solvers,a number of approaches rely on CNN-type architectures [9, 17, 49, 57]. Furthermore, it has beenshown that data-driven approaches can even be used to approximate the temporal evolution of fluidflows [56], to compute liquid splashing [51], artistic style-transfer [21], or to derive fluid dynamicsfrom reduced sets of parameters [20].
Graph-based Learning.
Graphs have proven to be a powerful representation for learning a widerange of tasks [13, 46]. In particular, it has been shown that graphs enable learning knowledgerepresentations [22], message passing [15], or to encode long-range dependencies, e.g., as foundin video processing [55]. A variety of methods uses graph-based representations to learn proper-ties of dynamic physical systems, e.g. for climate prediction [47], with an emphasis on individualobjects [8] and their relations [45], for partially observable systems [27], the prevalent interactionswithin physical systems [22], hierarchically-organized particle systems [35], or – more generally –physical simulation [43, 44]. While many of the existing approaches learn the time integration of2 iscretized Rod Forward Integration Constraint Projection
Figure 2: Illustration of the discretization of a single rod using several rod segments arranged alongits centerline (left). Each rod segment is described by its position and orientation within the gener-alized coordinates p i . The Lagrange multipliers λ i represent the interaction between rod segments.The forward integration path is illustrated in red (middle) and constraint projection in green (right).physical systems in an end-to-end manner, we use a graph network to predict the outcome of a PBDsolver for rod dynamics to enable more efficient computations. Position-based Dynamics and Elastic Rods.
PBD is a robust and efficient approach for simulatingposition changes of connected sets of vertices [4, 5, 28, 36]. Compared to forced-based methodsthat compute the force directly, the interaction between different vertices in PBD is realized by aconstraint projection step in an iterative manner. To avoid the dependency of the system’s stiffnesson the number of iterations and the time step size, an extended position-based dynamics approachwas introduced (XPBD) [28]. A number of methods exist that model the dynamic properties ofrods that can even simulate more complicated rod mechanics [38]. Moreover, particle systemswere employed to simulate the dynamics of rods [31] and, in particular, for the physically accuratesimulation of thin fibers [32] such as present in human hair or textiles. On a different trajectory,it has been recognized that rods can be simulated based on PBD [52]. The initial formulation wasimproved [23] by including the orientation of rod segments in the system’s state to account fortorsion effects. Later, the XPBD framework was utilized [10] to address the non-physical influenceof iteration numbers and steps sizes, which enables the more accurate simulation of elastic rods.
We propose a novel approach to simulate the temporal evolution of a dynamic system which con-sists of elastic rods. Each rod is discretized using several rod segments arranged along its centerline(Figure 2). For each rod segment, its state is described by its position, orientation, velocity and an-gular velocity. The state of the system is given as the set of the individual states of all rod segments.The simulation is carried out by employing PBD [36] directly manipulating the system’s state. Ori-entations are represented as quaternions allowing for a convenient implementation of bending andtwisting effects [23]. Extended PBD (i.e. XPBD) [28] is implemented to avoid that the rods’ stiff-nesses depend on the time step size and the number of iterations [10].The generalized coordinates of a rod segment i at time t is given by p i,t ∈ R × H , which includesits position x i,t ∈ R given in Cartesian coordinates and its orientation described by a quaternion q i,t ∈ H . Correspondingly, υ i,t ∈ R refers to the generalized velocity of the rod segment, whichincludes velocity and angular velocity. The system is continuously updated during the simulationby applying corrections ∆ p i = (∆ x i , ∆ φ i ) T ∈ R with position shifts ∆ x i ∈ R and orientationshifts ∆ φ i ∈ R representing the integration of the angular velocity. A single time integration step is presented in Algorithm 1. In the beginning (lines 1 to 4), general-ized velocity and generalized coordinates are updated by employing a symplectic Euler integrationstep. In this regard, a ext denotes the generalized acceleration due to the external net force, e.g., givenby gravity. XPBD [28] employs the Lagrange multiplier λ which is initialized as zero (line 5) alongwith the integrated generalized coordinates p ∗ . They are then used for the rod constraint projection(line 9). Collision detection results are stored in Coll r-r and
Coll r-p (line 6), where
Coll r-r includesall the pairs of two rod segments that potentially collide with each other and
Coll r-p includes in-formation of all rod segments that potentially collide with another object such as a plane. Withinseveral solver iterations, we alternate between rod constraint projection and the collision constrainprojection (lines 7 to 15). The rod constraints include shear-stretch and bend-twist constraints repre-senting the corresponding elastic energy. The Lagrange multiplier represents the interaction betweenrod segments. Figure 2 illustrated the discretization for a single rod into several interacting segments.The correction values ∆ p and ∆ λ in line 9 are computed by constraint projection [10, 23]. The gen- Please note, that ∆ q i = G ( q )∆ φ i ∈ R , in which the matrix G ( q ) ∈ R × describes the relationshipbetween quaternion velocity and angular velocity [3]. lgorithm 1 Numerical integration procedure updating p i,t (cid:55)→ p i,t +∆ t and υ i,t (cid:55)→ υ i,t +∆ t . for all rod segments do υ ∗ i ← υ i,t + ∆ t a ext p ∗ i ← p i,t + ∆ t H ( q i,t ) υ ∗ i with H ( q i,t ) := [ × , × ; × , G ( q i,t )] end for λ ← , p ← p ∗ ( Coll r-r , Coll r-p ) ← generateCollisionConstraints ( p , p ∗ ) for j ← to number of required solver iterations do for all rods do (∆ p , ∆ λ ) ← rodConstraintProjection ( p j , λ j ) λ j +1 ← λ j + ∆ λ p j +1 ← p j + ∆ p end for p j +1 ← updateCollisionConstraintProjection ( p j +1 , Coll r-r , Coll r-p ) j ← j + 1 end for for all rod segments do p i,t +∆ t ← p ji υ i,t +∆ t ← H T ( q i,t )( p i,t +∆ t − p i,t ) / ∆ t end for D ec od e r GN GN GN M ... E n c od e r rodConstrainedProjectioncorrection Guess Output
Figure 3: Illustration of our approach incorporating a network which consists of M graph networks(GN-blocks) into the position-based dynamics framework.eralized coordinates and Lagrange multipliers are updated for each rod (lines 8 to 12), and rod-rodand rod-plane collisions are addressed to update the generalized coordinates. For details about thecollision projection procedure, we refer to Macklin et al. [29]. For the non-collision case, the stepswithin line 6 and 13 are not needed, and less solver iterations are necessary in line 7.The most expensive part within Algorithm 1 involves the computation of the corrections of gener-alized coordinates and Lagrange multipliers (line 9). This projection step requires the solution of alinear system which is a linearization of a non-linear one so that the matrix depends on the system’sstate making it impossible to precompute its inverse. Instead, a new system at every point in timeis solved iteratively using the conjugated gradient (CG) solver. Such iterative solvers are widelyused in the context of physical simulations and regularly described as the de facto standard [2, 42]since they often show superior performance and usually scale well allowing for exploiting parallelhardware. However, we would like to point out that also highly efficient direct solvers can be foundin the literature [10].Instead of fully replacing the projection step in an end-to-end learning manner, we follow the strategyof accelerating it by first computing a guess (∆ p , ∆ λ ) ← correctionGuess ( p j ) , (1)for the iterative procedure (line 9) (∆ p , ∆ λ ) ← rodConstraintProjection ( p j , λ j , ∆ p , ∆ λ ) . (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 convergenceof the CG solver compared to the canonical initialization with zeros. We obtain our final frameworkby replacing line 9 in Algorithm 1 with Eq. (1) and Eq. (2) as illustrated in Figure 3. We name4he data-driven part COPINGNet (“COnstraint Projection INitial Guess Network”) which learns tocompute the correction guess.
COPINGNet is a graph network based architecture which we apply in order to compute initialguesses for ∆ p and ∆ λ . In this regard, we need to incorporate the rods’ state into a graph de-scription [13]. A graph G = ( V, E, U ) usually consists of nodes (or vertices) V , edges E as well asglobal features U . However, in our framework, global features U are not used. For example, gravitycould be a potentially meaningful global feature, but it also can be easily included as an externalacceleration. Hence, U = ∅ and the graph can just be represented as G = G ( V, E ) . In our case, therods’ segments within the scene are represented by the graph’s nodes while the interactions betweenthe rods’ segments are represented by the edges.COPINGNet provides a graph-to-graph mapping: G in → G out , from an input graph G in ∈ G in to anoutput graph G out ∈ G out . Nodes and edges of both graphs are equipped with specific features. Inthe case of the input graph, the node features describe the state of the rods’ segments, i.e. v in ,i = ( x i , q i , r i , ρ i , (cid:96) i , α i , f i , f i , f i ) T ∈ V in ⊆ R , in which the positions are denoted with x i ∈ R , the orientations with q i ∈ H , the radii with r i ∈ R > , the densities with ρ i ∈ R > , and the segment lengths with (cid:96) i ∈ R > . Moreover, asegment position indicator α i ∈ [0 , ⊂ R is included corresponding to a parameterization byarc length. Binary segment flag f i ∈ { , } , “left” end flag f i ∈ { , } and “right” end flag f i ∈ { , } are set to zero if the specific segment respectively the left or right segment of the rodis fixed and to one otherwise. The nodes of G in are given as the set of V in = ∪ ni =1 { v in ,i } , in which n = | V in | denotes the number of segments in the scene. The nodes of G out contain the correctionvalues of the generalized coordinates, i.e. v out ,i = ∆ p i ∈ V out ⊆ R , and V out = ∪ ni =1 { v out ,i } .While rod segments are represented as node features, it is natural to represent constraints betweenrod segments as edge features: e in ,i = ( ω i , Y i , T i ) T ∈ E in ⊆ R , in which the (rest) Darboux vector ω ∈ R describes the static angle of two rod segments, andYoung’s modulus Y ∈ R > and torsion modulus T ∈ R > are corresponding to extension,bending, and twist constraint parameters. The set of edges of the input graph is then given by E in = ∪ mi =1 { e in ,i } , in which m = | E in | denotes the number of interactions between different seg-ments. The correction of the Lagrange multiplier ∆ λ i ∈ R is stored in the output edges: e out ,i = ∆ λ i ∈ E out ⊆ R . The set of output edges is then given by E out = ∪ mi =1 { e out ,i } .The connectivity information C of each graph is stored in two vectors c sd and c rv containing thesender node index and the receiver node index of each corresponding edge. This concludes thespecification of the input graph G in = G ( V in , E in ) and the output graph G out = G ( V out , E out ) withconnectivity information C = ( c sd , c rv ) . In the following, we describe the structure of COPINGNet after we formalized its input and output inthe previous section. As illustrated in Figure 3, the network consists of an encoder network, multiplestacks of GN-blocks, and a decoder network. The graph network from Battaglia et al. [13] is usedas a basic element and denoted as a GN-block. Residual connection between blocks could improveperformance of neural networks in both CNN [18], and graph neural network [26]. As in relatedwork [44], we employ the residual connection between the GN-blocks, but our encoder/decodernetwork directly deals with the graph. The encoder network performs a mapping: G in → G latent and is implemented using two multi-layer perceptrons (MLPs), MLP edge : E in → E latent ⊆ R l and For a single rod in the scene which consists of N segments of equal lengths, for the i -th segment, we obtain α i = ( i − / ( N − for i ∈ { , , . . . , n } . rain/Val N Young’s Modulus Y Initial Angle φ Rod Length (cid:96) U d (10 ,
55) 10 a Pa, a ∼ U (4 . , . U (0 ◦ , . ◦ ) U (1 . m , . m ) Train/Val N Torsion Modulus G Helix Radius / Height Winding Number U d (45 , a Pa, a ∼ U (4 . , . U (0 . m , . m ) / U (0 . m , . m ) U (2 . , . Table 1: Specification of training and validation data sets for the two scenarios of an initially straightbending rod (top) and an elastic helix (bottom). The data sets are comprised of a number of datapoints (left) each describing the rod’s dynamics within t ∈ [0 s , t ] discretized with a time stepsize of ∆ t = 0 . s. The number of nodes N is sampled from a discrete uniform distribution U d while the remaining parameters are sampled from a continuous uniform distribution U . MLP node : V in → V latent ⊆ R l , in which l denotes the latent size. They work separately and thus E en = MLP edge ( E in ) and V en = MLP node ( V in ) . Edge features E in from the input are constant fora rod during the simulation and this results in constant encoded edge features E en , which could berecorded after the first run and used afterwards during inference. The edge feature e in ,i containsthe material parameters which could vary by different orders of magnitude. Hence, we normalizeYoung’s modulus and torsion modulus before feeding them into the network. After encoding, thegraph G en ( V en , E en ) ∈ G latent is passed to several GN-blocks with residual connections. Each GN-block also contains two MLPs. However, they use message passing taking advantage of neighbour-hood nodes’/edges’ information [13]. A number of M GN-blocks enable the use of neighbourhoodinformation with distances smaller or equal to M . The graph network performs a mapping withinthe latent space: G latent → G latent , and after M GN-blocks, we obtain G (cid:48) en ( V (cid:48) en , E (cid:48) en ) ∈ G latent .The decoder network performs a mapping: G latent → G out , which has a similar structure as theencoder network. Two MLPs MLP edge : E latent → E out and MLP node : V latent → V out compute E out = MLP edge ( E (cid:48) en ) and V out = MLP node ( V (cid:48) en ) . A tanh -function at the end of MLP edge and
MLP node is used restricting the output to the interval [ − , ⊂ R . COPINGNet learns the relativecorrection values. The generated data set is normalized, and the maximum correction value of thegeneralized coordinates and the Lagrange multipliers will be recorded as norms. The final correctionvalue is achieved by multiplying the relative correction values and the norms. This normalizationprocess damps the noise caused by the network and leads to a more stable performance. For sim-plicity, all the MLPs in different blocks have the same number of layers, and the same width aslatent size l within the latent layers. The input and output sizes of each MLP are consistent with thecorresponding node/edge feature dimensions. We generate training and validation data sets based on two scenarios: an initially straight bending rodand an elastic helix each fixed at one end oscillating under the influence of gravity. The specificationof these data sets is provided in Table 1. The PBD code is written in C++ [10] while the COPINGNetis implemented in PyTorch. The training is performed on an NVIDIA R (cid:13) Tesla R (cid:13) V100 GPU. Thetraining time varies from 8 to 30 hours for different architecture parameters. A constant learning rateof η = 0 . was used and a mean square loss function was employed. Our approach generalizesacross different initial conditions, geometries, discretizations, and material parameters. It is able torobustly generate various dynamical results as shown in the supplementary material (Figure 8). Architecture.
The evaluation of COPINGNet’s architecture, introduced in Section 3.2, is presentedin Figure 4. The performance is studied for different numbers of GN-blocks and MLP layers, anddifferent MLP widths (latent sizes). The measurements demonstrate that a larger number of GN-blocks usually increases the performance while the performance improvement is not longer signifi-cant for more than four GN-blocks. This is plausible and can be considered analogously to the sizeof the stencil of a numerical integrator: a larger number of GN-blocks means that a specific node cantake advantage of information gained from neighbourhoods further away. Correspondingly, a largerstencil can do so as well. However, increasing the size of the stencil usually does not result in moreaccurate results once the critical size is reached. This can be observed here as well. In contrast, it canbe clearly observed that deeper MLP networks do not improve the performance. This is consistentwith other research on graph networks [44]. However, increasing the MLP width can increase thenetwork’s ability to generalize. It is shown that the highest performance is achieved with the largestMLP width, especially in the case of the helix scenarios. Since increasing the number of GN-blocksand MLP width will lead to longer inference time, we compromise this by choose medium numbers.For further evaluations, we employ three GN-blocks, two MLP layers, and a MLP width of .6 C G I t e r a t i o n R a t i o ( B e n d i n g R o d ) C G I t e r a t i o n R a t i o ( H e li x ) Figure 4: Illustration of the evaluation of COPINGNet’s architecture (blue: benchmark architecture).The result is averaged from test simulations each running for time steps.
20 40 60 80 100 120 140Number of Nodes0.81.01.21.41.6 Sp ee d u p
25 50 75 100 125 150 175 200Number of Nodes 0.00.10.20.30.4 t i n f e r / t C G Figure 5: Illustration of the ratio of COPINGNet’s inference time t infer and the vanilla CG solver’srun time t CG (purple curves; right vertical axis) for the initially straight bending rod (left) and theelastic helix (right) simulations. Moreover, the black curves show the CG iteration number ratiowhile the red curves show the total speedup of the constraint projection when taking into accountCOPINGNet’s inference time (left vertical axis). The pink curves show the total speedup of theentire simulations. The orange and the green dashed lines indicate the lower and upper boundariesof the nodes number’s sampling range. The result is averaged from test simulations each runningfor time steps. Discretization.
We evaluate the influence of the system’s complexity defined by the number ofnodes discretizing a rod. Our approach specifically addresses the acceleration of the most expensivepart within PBD by providing an accurate initial guess of the constraint projection. In this regard,the costs coming with the computations within the COPINGNet part have to be taken into account aswell. Figure 5 (purple curves) shows the ratio of COPINGNet’s inference run time compared to therun time of the vanilla CG solver (initialized without COPINGNet) for different numbers of nodes.We can clearly observe that the additional run time added by COPINGNet is getting more trivialwhen the number of nodes is already relatively large. The measured ratio is used to indicate howmuch speedup the constraint projection solver can gain by making use of COPINGNet. Taking intoaccount COPINGNet’s inference run time, we can analyse the total (net) speedup of the constraintprojection. This is illustrated in Figure 5 (black and red curves) showing that for a small number ofnodes, the guess computed by COPINGNet does not reduce the computational costs that much sincethe inference costs can not yet be compensated. Once the number of nodes is increasing, a significantspeedup can be obtained of up to about for the constraint projection. More precisely, it keepson increasing until about the middle of the data set sampling range. This is natural since neuralnetworks perform better when getting closer to their training data set distribution. However, ourapproach also performs well when going far beyond the sampling range successfully demonstratingits ability to generalize. Since the constraint projection is the most time-consuming part of the entiresimulation, specifically for a large number of nodes, the speedup of the whole simulation becomessimilar to the one of the constraint projection as shown in Figure 5 (pink curves). For a small numberof nodes, assembling the linear system takes more time, which results in a smaller speedup.
Temporal Evolution.
We analyze the required number of CG iterations for the vanilla constraintprojection compared to the one accelerated using COPINGNet over time as shown in Figure 6. Weclearly observe that we obtain a significant speedup in total. As stated in Table 1, our training datacontains dynamical simulations over time steps. In this range, we naturally observe the highestspeedup. However, we obtain a significant speedup even after hundreds of time steps demonstratingCOPINGNet’s ability to generalize. Long-term Stability.
The previously described result shown in Figure 6 can also be con-sidered as a demonstration for long-term stability. This is a consequence of PBD’s stability7
200 400 600 800 1000Time Frame Number1.01.52.0 C G I t e r a t i o n R a t i o Figure 6: The ratio of the required iteration numbers within the CG solver carrying out the constraintprojection. For the initially straight bending rod (orange curve) simulation, the parameters φ = 0 ◦ , N = 30 , (cid:96) = 3 . m, and Y = 1 . · Pa are used. For the elastic helix (green curve) simulation,the parameters HR = 0 . m, HH = 0 . m, HW = 2 . , G = 1 . · Pa, and N = 60 are used. l / l Y =1.0 10 Pa Y =2.0 10 Pa Y =5.0 10 Pa Y =1.0 10 Pa Y =2.0 10 Pa Y =5.0 10 Pa Y =1.0 10 Pa Figure 7: Illustration of the relativechange of the total rod length.which is accurately preserved within our approach. Thisis due to the fact that COPINGNet’s output serves as theinitial guess of the iteratively enforced constraint projec-tion. In contrast, if the constraint projection is completelyreplaced by COPINGNet, significant stability issues canbe observed as error accumulation takes place. This is il-lustrated in Figure 7 showing the temporal evolution ofthe relative change of the total rod length simulated (a)using COPINGNet supplying the initial guess for the con-straint projection, and (b) using COPINGNet fully replac-ing the conventional constraint projection step. An ini-tially straight rod bending under the influence of gravity is simulated using the parameters φ = 0 ◦ , N = 30 , (cid:96) = 4 . m, and varying Young’s modulus Y corresponding to dynamic sequences in whichthe rod lengths should be approximately constant. While this is clearly the case in setup (a) even forthousand frames as shown in Figure 6, the length can not be preserved at all in setup (b) as shownin Figure 7 (colored curves correspond to (b); black visually overlapping curves to (a) for varyingYoung’s modulus). After the time steps covered within the training data set are exceeded, weobtain a strongly implausible result demonstrating the superiority of setup (a). Collisions.
We further showcase that our approach can handle discontinuous events such as col-lisions between individual rods and other objects. This is illustrated in Figure 1 showcasing thecollision of an elastic helix with the ground plane of the scene and the collision between a rodfalling on another rod which is fixed at both ends. Since collisions are handled “outside” the neu-ral network part within our approach, they can be treated as described in the literature. Collisiondetection is efficiently implemented using the hierarchical spatial hashing scheme according to Eitzand Lixu [12]. Rod-rod collisions are then treated with line-to-line distance constraints, and col-lisions with the ground plane using half-space constraints. Moreover, frictional effects are imple-mented according to Macklin et al. [29]. For the scene comprising the collision of the elastic helix( HR = HH = 0 . m, HW = 2 . , G = 1 . · Pa, N = 50 ) with the ground plane, we mea-sure a total speedup of approx. . In the case of two colliding rods ( φ = 0 ◦ , N ∈ { , } , (cid:96) ∈ { . m , . m } , Y = 1 . · Pa), a speedup of approx. is measured. We presented a novel approach for accelerating iterative solvers with GNs. Unlike existing end-to-end learned approaches, our network-enabled solver ensures long-term stability inherited fromtraditional solvers for physical systems. We evaluated the architecture of our GN and demonstratedthat it is able to generalize across different initial conditions, rod discretizations, and material pa-rameters, and is even able to handle discontinuous effects such as collisions. Although, it is apparentthat our method extends to other physical systems modeled with PBD, we did not conclusively ex-plore different network topologies, i.e. vertices with multiple neighbors, or various segment lengths(e.g. by simulating vegetation such as trees). Our approach to accelerate iterative solvers with GNsopens multiple avenues for future work. For one, it would be interesting to explore mechanical sys-tems describing general deformable (e.g. textiles) or volumetric objects, which have been simulatedwith PBD. Second, our approach can be applied to other iterative methods, such as in finite ele-ments analysis or in the context of linear complementarity problems (LCP). This would allow us toaccelerate physical simulations, when iterative solvers are applied, without compromising stability.8 roader Impact
Apart from being useful for accelerating the simulation of mechanical systems, e.g., using itera-tive methods such as in position-based dynamics, finite elements analysis or in the context of linearcomplementarity problems, we see a broader impact of our work in contributing to the further de-velopment of synergy effects of machine learning (e.g. graph learning) and physical simulations.
Acknowledgements
This work has been supported by KAUST under individual baseline and center partnership funding.The authors are grateful to Torsten H¨adrich, Libo Huang, and Jonathan Klein for helpful discussions.
References [1] Jernej Barbiˇc and Doug L. James. Real-time subspace integration for st. venant-kirchhoffdeformable models.
ACM Trans. Graph. , 24(3):982990, July 2005.[2] R. Barrett, M. Berry, T. F. Chan, J. Demmel, J. Donato, J. Dongarra, V. Eijkhout, R. Pozo,C. Romine, and H. Van der Vorst.
Templates for the Solution of Linear Systems: BuildingBlocks for Iterative Methods, 2nd Edition . SIAM, Philadelphia, PA, 1994.[3] Jan Bender, Kenny Erleben, and Jeff Trinkle. Interactive Simulation of Rigid Body Dynamicsin Computer Graphics.
Computer Graphics Forum , 33(1):246–270, 2014.[4] Jan Bender, Matthias M¨uller, and Miles Macklin. Position-based simulation methods in com-puter graphics. In
EUROGRAPHICS 2017 Tutorials . Eurographics Association, 2017.[5] Jan Bender, Matthias Mller, Miguel A. Otaduy, Matthias Teschner, and Miles Macklin. A sur-vey on position-based simulation methods in computer graphics.
Computer Graphics Forum ,33(6):228–251, 2014.[6] Jens Berg and Kaj Nystrm. A unified deep artificial neural network approach to partial differ-ential equations in complex geometries.
Neurocomputing , 317:28 – 41, 2018.[7] Mikl´os Bergou, Max Wardetzky, Stephen Robinson, Basile Audoly, and Eitan Grinspun. Dis-crete elastic rods.
ACM Transactions on Graphics , 27(3):63:1–63:12, 2008.[8] Michael B Chang, Tomer Ullman, Antonio Torralba, and Joshua B Tenenbaum. A composi-tional object-based approach to learning physical dynamics. arXiv preprint arXiv:1612.00341 ,2016.[9] Mengyu Chu and Nils Thuerey. Data-Driven Synthesis of Smoke Flows with CNN-BasedFeature Descriptors.
ACM Trans. Graph. , 36(4), July 2017.[10] Crispin Deul, Tassilo Kugelstadt, Marcel Weiler, and Jan Bender. Direct position-based solverfor stiff rods.
Computer Graphics Forum , 37(6):313–324, 2018.[11] M. W. M. G. Dissanayake and N. Phan-Thien. Neural-network-based approximations for solv-ing partial differential equations.
Communications in Numerical Methods in Engineering ,10(3):195–201, 1994.[12] Mathias Eitz and Lixu Gu. Hierarchical spatial hashing for real-time collision detection.
IEEEInternational Conference on Shape Modeling and Applications 2007 (SMI ’07) , pages 61–70,2007.[13] Peter W. Battaglia et al. Relational inductive biases, deep learning, and graph networks.
CoRR ,abs/1806.01261, 2018.[14] Stefan Feess, Kathrin Kurfiss, and Dominik L. Michels. Accurate Simulation of Wound Heal-ing and Skin Deformation. In
Proceedings of the ACM SIGGRAPH/Eurographics Symposiumon Computer Animation , SCA 16, page 129137, Goslar, DEU, 2016. Eurographics Associa-tion.[15] Justin Gilmer, Samuel S. Schoenholz, Patrick F. Riley, Oriol Vinyals, and George E. Dahl.Neural message passing for quantum chemistry.
CoRR , abs/1704.01212, 2017.916] Radek Grzeszczuk, Demetri Terzopoulos, and Geoffrey Hinton. Neuroanimator: Fast neuralnetwork emulation and control of physics-based models. In
Proceedings of the 25th AnnualConference on Computer Graphics and Interactive Techniques , SIGGRAPH 98, page 920,New York, NY, USA, 1998. Association for Computing Machinery.[17] 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 Knowl-edge Discovery and Data Mining , KDD ’16, pages 481–490, New York, NY, USA, 2016.ACM.[18] Kaiming He, Xiangyu Zhang, Shaoqing Ren, and Jian Sun. Deep residual learning for imagerecognition.
CoRR , abs/1512.03385, 2015.[19] Ruben Iba˜nez, Domenico Borzacchiello, Jose Vicente Aguado, Emmanuelle Abisset-Chavanne, Elias Cueto, Pierre Ladeveze, and Francisco Chinesta. Data-driven non-linearelasticity: Constitutive manifold construction and problem discretization.
Comput. Mech. ,60(5):813826, November 2017.[20] B. Kim, V. C. Azevedo, N. Thuerey, T. Kim, M. Gross, and B. Solenthaler. Deep fluids: Agenerative network for parameterized fluid simulations.
CGF , 38(2):59–70, 2019.[21] B. Kim, V.C. Azevedo, M. Gross, and B. Solenthaler. Lagrangian neural style transfer forfluids.
ACM Transaction on Graphics (SIGGRAPH) , 2020.[22] Thomas Kipf, Ethan Fetaya, Kuan-Chieh Wang, Max Welling, and Richard Zemel. Neuralrelational inference for interacting systems. arXiv preprint arXiv:1802.04687 , 2018.[23] Tassilo Kugelstadt and Elmar Schoemer. Position and orientation based cosserat rods. In
Proceedings of the 2016 ACM SIGGRAPH/Eurographics Symposium on Computer Animation .Eurographics Association, 2016.[24] Lubor Ladick´y, SoHyeon Jeong, Barbara Solenthaler, Marc Pollefeys, and Markus Gross.Data-driven fluid simulations using regression forests.
ACM Trans. Graph. , 34(6), 2015.[25] I. E. Lagaris, A. Likas, and D. I. Fotiadis. Artificial neural networks for solving ordinary andpartial differential equations.
IEEE Transactions on Neural Networks , 9(5):987–1000, 1998.[26] Guohao Li, Matthias Mller, Ali Thabet, and Bernard Ghanem. DeepGCNs: Can GCNs Go asDeep as CNNs? In
The IEEE International Conference on Computer Vision (ICCV) , 2019.[27] Yunzhu Li, Jiajun Wu, Jun-Yan Zhu, Joshua B. Tenenbaum, Antonio Torralba, and RussTedrake. Propagation networks for model-based control under partial observation, 2018.[28] Miles Macklin, Matthias M¨uller, and Nuttapong Chentanez. Xpbd: Position-based simulationof compliant constrained dynamics. In
Proceedings of the 9th International Conference onMotion in Games , MIG 16, page 4954, New York, NY, USA, 2016. ACM.[29] Miles Macklin, Matthias M¨uller, Nuttapong Chentanez, and Tae-Yong Kim. Unified particlephysics for real-time applications.
ACM Trans. Graph. , 33(4), July 2014.[30] Dominik L. Michels and Mathieu Desbrun. A Semi-analytical Approach to Molecular Dynam-ics.
Journal of Computational Physics , 303:336 – 354, 2015.[31] Dominik L. Michels, Vu Thai Luan, and Mayya Tokman. A Stiffly Accurate Integrator forElastodynamic Problems.
ACM Trans. Graph. , 36(4), July 2017.[32] Dominik L. Michels, J. Paul T. Mueller, and Gerrit A. Sobottka. A Physically Based Ap-proach to the Accurate Simulation of Stiff Fibers and Stiff Fiber Meshes.
Comput. Graph. ,53(PB):136146, December 2015.[33] Dominik L. Michels, Gerrit A. Sobottka, and Andreas G. Weber. Exponential Integrators forStiff Elastodynamic Problems.
ACM Transactions on Graphics , 33(1):7:1–7:20, 2014.[34] Siddhartha Mishra. A machine learning framework for data driven acceleration of computa-tions of differential equations, 2018.[35] Damian Mrowca, Chengxu Zhuang, Elias Wang, Nick Haber, Li Fei-Fei, Joshua B. Tenen-baum, and Daniel L. K. Yamins. Flexible neural representation for physics prediction, 2018.[36] Matthias Mller, Bruno Heidelberger, Marcus Hennix, and John Ratcliff. Position based dy-namics.
Journal of Visual Communication and Image Representation , 18(2):109 – 118, 2007.1037] Andrew Nealen, Matthias Mller, Richard Keiser, Eddy Boxerman, and Mark Carlson. Physi-cally based deformable models in computer graphics.
Computer Graphics Forum , 25(4):809–836, 2006.[38] Dinesh K. Pai. Strands: Interactive simulation of thin solids using cosserat models.
ComputerGraphics Forum , 21(3):347–352, 2002.[39] S¨oren Pirk, Michal Jarzabek, Torsten H¨adrich, Dominik L. Michels, and Wojciech Palubicki.Interactive wood combustion for botanical tree models.
ACM Trans. Graph. , 36(6):197:1–197:12, November 2017.[40] M. Raissi, P. Perdikaris, and G.E. Karniadakis. Physics-informed neural networks: A deeplearning framework for solving forward and inverse problems involving nonlinear partial dif-ferential equations.
Journal of Computational Physics , 378:686 – 707, 2019.[41] Maziar Raissi and George Em Karniadakis. Hidden physics models: Machine learning ofnonlinear partial differential equations.
Journal of Computational Physics , 357:125 – 141,2018.[42] Y. Saad.
Iterative Methods for Sparse Linear Systems . Society for Industrial and AppliedMathematics, USA, 2nd edition, 2003.[43] Alvaro Sanchez-Gonzalez, Victor Bapst, Kyle Cranmer, and Peter Battaglia. Hamiltoniangraph networks with ode integrators, 2019.[44] Alvaro Sanchez-Gonzalez, Jonathan Godwin, Tobias Pfaff, Rex Ying, Jure Leskovec, and Pe-ter W. Battaglia. Learning to simulate complex physics with graph networks, 2020.[45] Alvaro Sanchez-Gonzalez, Nicolas Heess, Jost Tobias Springenberg, Josh Merel, Martin Ried-miller, Raia Hadsell, and Peter Battaglia. Graph networks as learnable physics engines forinference and control, 2018.[46] F. Scarselli, M. Gori, A. C. Tsoi, M. Hagenbuchner, and G. Monfardini. The graph neuralnetwork model.
IEEE Transactions on Neural Networks , 20(1):61–80, Jan 2009.[47] Sungyong Seo and Yan Liu. Differentiable physics-informed graph networks.
CoRR ,abs/1902.02950, 2019.[48] Justin Sirignano and Konstantinos Spiliopoulos. Dgm: A deep learning algorithm for solvingpartial differential equations.
Journal of Computational Physics , 375:1339 – 1364, 2018.[49] Jonathan Tompson, Kristofer Schlachter, Pablo Sprechmann, and Ken Perlin. Acceleratingeulerian fluid simulation with convolutional networks.
CoRR , abs/1607.03597, 2016.[50] Rahul Trivedi, Logan Su, Jesse Lu, Martin F Schubert, and Jelena Vuckovic. Data-drivenacceleration of photonic simulations, 2019.[51] Kiwon Um, Xiangyu Hu, and Nils Thuerey. Liquid splash modeling with neural networks,2017.[52] Nobuyuki Umetani, Ryan Schmidt, and Jos Stam. Position-based elastic rods. In
Proceedingsof the ACM SIGGRAPH/Eurographics Symposium on Computer Animation , SCA ’14, pages21–30, Aire-la-Ville, Switzerland, Switzerland, 2014. Eurographics Association.[53] B. Ummenhofer, L. Prantl, N. Thuerey, and V. Koltun. Lagrangian fluid simulation with con-tinuous convolutions. In
ICLR , 2020.[54] Huamin Wang, James F. O’Brien, and Ravi Ramamoorthi. Data-driven elastic models forcloth: Modeling and measurement. In
ACM SIGGRAPH 2011 Papers , SIGGRAPH ’11, pages71:1–71:12, New York, NY, USA, 2011. ACM.[55] Xiaolong Wang, Ross B. Girshick, Abhinav Gupta, and Kaiming He. Non-local neural net-works.
CoRR , abs/1711.07971, 2017.[56] Steffen Wiewel, Moritz Becher, and Nils Thuerey. Latent-space physics: Towards learning thetemporal evolution of fluid flow.
CoRR , abs/1802.10123, 2018.[57] X. Xiao, Y. Zhou, H. Wang, and X. Yang. A novel cnn-based poisson solver for fluid simula-tion.
IEEE Transactions on Visualization and Computer Graphics , 26(3):1454–1465, 2020.11 upplementary Material
Dynamic Results
Figure 8 illustrates the temporal evolution of different scenarios (a) to (d) for the initially straightbending rod, and (e) to (f) for the elastic helix. The default parameters of the initially straightbending rod are φ = 0 ◦ , N = 30 , (cid:96) = 4 . m, and Y = 1 . · Pa. In (a), wemodify φ ∈ { . ◦ , . ◦ , . ◦ , . ◦ } . In (b), we modify N ∈ { , , , } . In (c),we modify (cid:96) ∈ { . m , . m , . m , . m , . m , . m } . In (d), we modify Y ∈ { . · Pa , . · Pa , . · Pa , . · Pa , . · Pa , . · Pa } . The default parametersof the elastic helix are HR = 0 . m (helix radius), HH = 0 . m (helix height), HW = 2 . (winding number), G = 1 . · Pa, and N = 60 . In (e), we modify ( HR , HH , HW ) ∈{ (0 . m , . m , . m ) , (0 . m , . m , . m ) , (0 . m , . m , . m ) , (0 . m , . m , . m ) } . In (f),we modify N ∈ { , , , } . In (g), we modify G ∈ { . · Pa , . · Pa , . · Pa , . · Pa } . For each experiment, the rod colors indicate the corresponding parameters in the followingorder: blue , orange , green , red , purple , brown. These results are obtained by employing three GNblocks, two MLP layers, and a MLP width of . a)b)c)d)e)f)g)a)b)c)d)e)f)g)