Small Steps In Physics Simulation - GitHub Pages
Small Steps in Physics SimulationMiles MacklinKier StoreyMichelle LuNVIDIAUniversity of comNVIDIAmichellel@nvidia.comPierre TerdimanNuttapong ChentanezStefan nvidia.comNVIDIAsjeschke@nvidia.comMatthias MüllerNVIDIAmatthiasm@nvidia.comABSTRACTIn this paper we re-examine the idea that implicit integrators withlarge time steps offer the best stability/performance trade-off forstiff systems. We make the surprising observation that performinga single large time step with n constraint solver iterations is lesseffective than computing n smaller time steps, each with a singleconstraint solver iteration. Based on this observation, our approachis to split every visual time step into n substeps of length t/n andto perform a single iteration of extended position-based dynamics(XPBD) in each such substep. When compared to a traditionalimplicit integrator with large time steps we find constraint error anddamping are significantly reduced. When compared to an explicitintegrator we find that our method is more stable and robust for awider range of stiffness parameters. This result holds even whencompared against more sophisticated implicit solvers based onKrylov methods. Our method is straightforward to implement, andis not sensitive to matrix conditioning nor is it to overconstrainedproblems.CCS CONCEPTS Computing methodologies Simulation by animation; Interactive simulation.Figure 1: High resolution cloth consisting of 150k particles,and 896k spring constraints draped over a Stanford bunny.With 1 substep and 30 XPBD iterations the simulation takes12.4ms/frame but shows visible stretching (left). With 30substeps, each performing only 1 XPBD iteration the simulation takes 13.5ms/frame but shows significantly less stretching and higher material stiffness (right). In both cases wehave performed collision detection once per-frame.KEYWORDSPhysics-based animation, real-time simulationACM Reference Format:Miles Macklin, Kier Storey, Michelle Lu, Pierre Terdiman, Nuttapong Chentanez, Stefan Jeschke, and Matthias Müller. 2019. Small Steps in Physics Simulation. In SCA ’19:The ACM SIGGRAPH / Eurographics Symposium on Computer Animation (SCA ’19), July 26–28, 2019, Los Angeles, CA, USA. ACM, NewYork, NY, USA, Article 39, 7 pages. https://doi.org/10.1145/3309486.3340247Permission to make digital or hard copies of all or part of this work for personal orclassroom use is granted without fee provided that copies are not made or distributedfor profit or commercial advantage and that copies bear this notice and the full citationon the first page. Copyrights for components of this work owned by others than theauthor(s) must be honored. Abstracting with credit is permitted. To copy otherwise, orrepublish, to post on servers or to redistribute to lists, requires prior specific permissionand/or a fee. Request permissions from permissions@acm.org.SCA ’19, July 26–28, 2019, Los Angeles, CA, USA 2019 Copyright held by the owner/author(s). Publication rights licensed to ACM.ACM ISBN 978-1-4503-6677-9/19/07. . . DUCTIONThe simulation of physical systems is a challenging problem ininteractive computer graphics. Ideal algorithms should not onlyreproduce the underlying physical model, they should be robustand efficient enough for real-time applications and user interaction.Many methods exist to evolve a physical system forward in time,however they can broadly be split into two categories: explicit andimplicit methods.Explicit time-integration is rarely used in physical simulationsbecause it is only conditionally stable. This means that the simulation can diverge any time if certain conditions are not met. For thisreason, implicit integrators are often preferred due to their stability.However, in contrast to explicit integration, where the state of thenext time step can be computed directly from the current state, forimplicit integration a system of non-linear equations must be solved
SCA ’19, July 26–28, 2019, Los Angeles, CA, USAM. Macklin, K. Storey, M. Lu, P. Terdiman, N. Chentanez, S. Jeschke, and M. Müllerat every time step. Typically this is done using a Newton methodthat repeatedly linearizes the non-linear system. Each linear system is then solved by a global solver such as Conjugate Gradients(CG) or direct methods, and the solution is used to get closer to thenon-linear one. First order implicit methods introduce numericaldamping, and are often more computationally expensive than explicit methods. In addition, they may not guarantee a time-boundon convergence, which makes them less attractive for interactiveapplications.For real-time applications, local iterative relaxation methodssuch as Position-based Dynamics (PBD) [Müller et al. 2007] arepopular. Unlike global solvers which treat the system of equationsas a whole, local solvers such as the Projected Gauss-Seidel (PGS)iteration used in PBD handle each equation individually, one afterthe other. These methods are known to be very robust on typicalproblems. This is due in part to the fact that during a global solve thesystem matrix is frozen, meaning all constraint gradients are heldfixed. This can cause the solution to move far from the constraintmanifold. In contrast, non-linear PGS methods work on the systemof non-linear equations directly. Each constraint projection usesthe current gradients, which minimizes overshooting. In addition,relaxation methods are robust for overconstrained systems that posea challenge for global linear solvers, which may fail to converge atall in such scenarios. Furthermore, solving each equation separatelyallows PGS to trivially handle unilateral (inequality) constraints.Here, only constraints which violate the inequality conditions (theactive set) are projected. This set can change in every iteration andeven after each constraint projection, which prevents sticking. Tosimulate the two-way coupling of different features like liquids andrigid bodies, their list of constraints are simply concatenated orinterleaved, allowing fine-grained coupling [Macklin et al. 2014;Stam 2009].Despite their advantages, due to their local nature, relaxationmethods propagate error corrections more slowly than if the equations are solved simultaneously, making them less suitable for stiffproblems. They also suffer from numerical dissipation like traditional implicit methods. The aim of this project was to increase theconvergence rate and energy preservation of local solvers whilekeeping all their advantages. We derived our solution from a rathersurprising observation. There are two ways to increase the accuracyof a simulation: either decrease the time step size, or increase thenumber of solver iterations. Both increase the amount of computational work that has to be done. Consequently, if we want to keepthe work constant we have to change both. If we decide to increasethe number of solver iterations or the complexity of the solver,we have to increase the time step size to stay within the computational budget. On the other hand, if we decide to take smaller timesteps, then we have to decrease the number of solver iterations.The question we asked was which direction is more effective. Is itmore effective to (a) solve one difficult problem accurately, or (b)many simpler problems approximately. Since the work of Baraffet al. [1998], the commonly accepted knowledge in the computergraphics community has been to prefer large steps and implicitmethods for stiff problems.However, in our studies we found that (b) is significantly moreeffective than (a). In fact, for PGS the optimum lies at the far end of(b), i.e. taking as many small steps as the time budget allows whileTable 1: Summary of the relative strengths and weaknessesof each method.MethodSemi-Implicit EulerImplicit EulerXPBDSmall StepsStability Efficiency Simplicity Energy performing a single iteration in each step. Based on this observationwe split a time step of size t into n substeps of size t/n andperform one iteration of Extended Position Based Dynamics (XPBD)in each substep. The effect of replacing iterations by substeps is sosubstantial that one-step XPBD is competitive with sophisticatedglobal matrix solvers in terms of convergence. Intuitively, a singlepass over the constraints seems hardly enough to yield meaningfulinformation about forces and torques, but our results show thatthis is actually not the case. Indeed, as was shown by Macklin etal [2016], the first iteration of XPBD is equivalent to the first stepof a Newton solver operating on the backward Euler equations.Thus, while our single iteration method has a close computationalresemblance to explicit integration, it comes with all the stabilityproperties of an implicit method.It is common in computer graphics literature to compare explicit integration with small time steps against implicit integrationwith large time steps. But we are not aware of works that explorethe middle-ground: approximate implicit integration with smalltime steps. The contributions of this work are a new approach tosimulation, but more importantly, a study on the effectiveness of decreasing the time step size and the investigation of single iterationimplicit integration.In Table 1 we broadly summarize these trade-offs. Explicit methods are simple, and can be efficient for moderately stiff problems.Traditional implicit methods such as backwards Euler are stable, buthave poor energy conservation, and are generally more expensive.Iterative implicit methods like XPBD are stable, and efficient formoderate stiffness, however they struggle to achieve high stiffness,and also suffer from numerical damping. Our proposed method improves the convergence and energy preservation of XPBD througha simple modification of the underlying algorithm.2RELATED WORKThe use of implicit integration for forward dynamics in computergraphics dates back to a seminal work by Terzopoulos et al. [1988;1987]. They utilize the alternating-direction implicit (ADI) method[Peaceman and Rachford 1955] to solve some of the forces implicitly. Since that time, the use of explicit integration became popularuntil Baraff and Witkin [1998] proposed a implicit backward Eulerscheme for handling all the forces implicitly, including dampingin cloth simulation. This provides a method that is stable even forlarge time step sizes. Desbrun et al. [1999] sped up the computationby using a predictor-corrector approach to compute an approximate solution to implicit integration. These approaches, however,suffer from artificial numerical damping. A second-order backwarddifference formula (BDF) was used for cloth simulation in [Choi andKo 2002] to reduce this numerical damping. Bridson et al. [2003]
Small Steps in Physics SimulationSCA ’19, July 26–28, 2019, Los Angeles, CA, USAdemonstrated the use of mixed implicit/explicit integration (IMEX)for cloth simulation. They used explicit integration for treating theelastic force and implicit integration for damping forces, yieldinga central Newmark scheme [Newmark 1959]. Along with the useof strain limiting, the method is stable for moderate time step sizeand does not suffer much from numerical damping. An IMEX typemethod was also used for a particle system simulation in [Eberhardt et al. 2000]. Fierz et al. [2011] combined the use of explicitand implicit integrators using element-wise criteria. In all theseworks, the use of an implicit integrator that requires a global linearsolver is commonplace.Variational approaches can also be used to derive integrators[Kane et al. 2000; Kharevych et al. 2006; Marsden and West 2001]with excellent energy preservation and controlled damping. Depending on the quadrature rule used for converting the continuousLangrangian to a discrete version, one can arrive at explicit or implicit integration of varying orders of accuracy. While energy isstable, the resulting animation can still oscillate and produce unnatural high frequency vibration. The trade-offs between explicitand implicit integration still apply. Another interesting alternativeis an exponential integrator [Michels et al. 2014] that combines ananalytic solution of the linear part of the ODE with a numericalmethod for the nonlinear parts. However, as with all integratorsstability is only guaranteed in the linear regime.Projective Dynamics [Bouaziz et al. 2014] produces impressivereal time simulation results by combining local and global solves.Recently Dinev et al. [2018] proposed an energy control strategyfor Projective Dynamics by mixing implicit midpoint with eitherforward Euler or Backward Euler. Nonetheless, the global solve stepneeds to pre-factorize the system matrix to be fast, which preventssimulated meshes from changing topology at runtime, for example.Asynchronous integrators [Lew et al. 2003; Thomaszewski et al.2008] and contact handling [Harmon et al. 2009; Zhao et al. 2016] allow for varying the time step over the simulation domain. However,the required computation can fluctuate greatly over time, which isnot desirable in real-time simulations.The usefulness of a single Newton solver step over explicit integration was also reported in [Gast et al. 2015], but they did notexplore utilizing this observation further. The idea of using thepredicted state in the next time step for force computation similarto an implicit integrator is also used in the context of PD controllerin [Tan et al. 2011], resulting in a more stable simulation.3TIME INTEGRATIONWe write our equations of motion using an implicit position-leveltime discretization as follows:M(xn 1 x̃) C(xn 1 )T λn 1 0(1)C(xn 1 ) α̃ λn 1 0.(2)Here M is the system mass matrix and xn 1 is the system stateat the end of the nth time step. C is a vector of constraint functions, C it’s gradient with respect to system coordinates, and λn 1 theassociated Lagrange multipliers. Constraints are regularized witha compliance matrix α̃ that results from factorizing a quadraticenergy potential [Macklin et al. 2016]. The predicted or inertialposition x̃ is obtained by an explicit integration of external forces:x̃ xn tvn t 2 M 1 fex t (xn ).(3)Examining equation (3), we make the observation that the effectof external forces on positions is proportional to t 2 . This is dueto the discretization of a second order differential equation, and ithas a strong influence on the error committed in a single step. Forexample, halving the time step will result in a quarter of the positionerror, and so on. This simple fact is what motivates our method,and makes smaller time steps so effective at reducing positionalerror, as we show in Section 6.Algorithm 1 Substep XPBD simulation loop1:2:3:4:5:6:7:8:9:10:11:12:13:perform collision detection using xn , vn tf ts ns t epswhile n nst eps dopredict position x̃ xn ts vn ts2 M 1 fex t (xn )for all constraints docompute λ using Eq (7)compute x using Eq (4)update λn 1 λ (optional)update xn 1 x x̃end for update velocities vn 1 t1s xn 1 xnn n 1end while14:4CONSTRAINT SOLVETo enforce the constraints on the system coordinates we use XPBD[Macklin et al. 2016] which performs a position projection for aconstraint with index i as follows,(4) x M 1 Ci (x)T λi .Where the associated Lagrange multiplier increment is given by, Ci (x) α̃ i λi.(5) Ci M 1 CiT α̃ iPosition-based dynamics would typically repeat this update multiple times per-constraint in a Gauss-Seidel or Jacobi fashion. Ouridea is to instead divide the whole frame’s time step t f into nsubsteps, t f ts ,(6)nst eps λi and then perform a single constraint iteration for that substep.This can be thought of as an approximate, or inexact implicit solvefor each substep. However, by dividing the time-step we benefitfrom the dependence of position error on t 2 in the discrete equations of motion.Since we perform a single iteration per-substep and the initialLagrange multiplier is always zero, the numerator in (5), can besimplified as follows: Ci (x) λi .(7) Ci M 1 CiT α̃ i
SCA ’19, July 26–28, 2019, Los Angeles, CA, USAM. Macklin, K. Storey, M. Lu, P. Terdiman, N. Chentanez, S. Jeschke, and M. 5-30010 iters10 substeps100 iters100 002003004005006007008009001000Figure 3: Left: The stretch in a 1D chain of particles connected by distance constraints hanging under gravity. Substeps (red) are significantly more effective at reducingstretch than the equivalent number of iterations (blue).Right: The same test with a large mass (105 kg) attached tothe bottom of the chain, emphasizing the effect.Figure 2: A position-based fluid simulation consisting of936k particles. Using 10 iterations the fluid shows significant compression and highly damped motion. In contrast,substepping with 10 iterations results in a visibly stifferfluid with more dynamic motion.10 124022010 020010 -118010 -216010 -314010 -4IterationsSubstepsWhen using a single iteration per-step it is not required to store theLagrange multipliers, however they can be useful to provide forceestimates back to the user. In the results section we demonstrate thatthese force estimates remain accurate even when using only a singleiteration per substep (Figure 4). An overview of our simulation loopis given in Algorithm 1.4.1DampingReducing the time step reduces the amount of numerical dissipationin the integrator. For this reason it becomes important to explicitlyinclude damping in our constraint model. We do this using theXPBD formulation: λi Ci (x) γi Ci (x xn )(1 γi ) Ci M 1 CiT α̃ i.(8)Given β i ts2 βi , a time step scaled damping parameter forα̃ β i iconstraint i, we then define γi t. We refer the reader tospublications by Macklin et al. [2016], and Servin et al. [2006] forthe derivation.4.2Collision DetectionDecreasing time step size typically improves the accuracy of collision detection. However, performing collision detection every timestep adds a significant computational overhead. A key idea thatmakes our substepping approach feasible is to amortize this costby performing collision detection once per-frame and re-using thecontact set over multiple substeps. To do this, we first predict thestate of the system using the current velocity and the whole frame’stime step, t f . We use this trajectory to detect potential collisionsand generate contact constraints accordingly.10 -5020406080100120140160180100 Iterations100 Substeps1202000102030405060708090100Figure 4: Left: We plot the residual error at frame 1000 forvarying iteration and substep counts. We see approximatelytwo orders of magnitude lower error for the equivalent number of substeps. Right: We plot the force estimate (Lagrangemultiplier) for the constraint at the top of the chain overtime. The true value is 190N. Surprisingly, even performinga single iteration per-substep provides accurate force estimates.Performing one collision detection phase per-frame could resultin missed collisions due to changing trajectories. Missed collisionswould then be processed the next frame, however to minimize thiseffect we generate contacts between features that come within someuser-controlled margin distance. This could be further addressed byupdating the contact set in a manner similar to constraint manifoldrefinement (CMR) [Otaduy et al. 2009].4.3ContactWe treat contacts as inelastic and prevent interpenetration betweenbodies through inequality constraints of the formCn (x) nT [a(x) b(x)] 0.(9)Here a and b are points on a rigid or deformable body and n R3is the contact plane normal. The contact normal may be fixed inworld space, or may itself be a function of the system coordinates.One artifact of reducing the time step is that any initial penetration error between bodies will be converted to large separatingvelocities by our implicit solve. A solution to
The simulation of physical systems is a challenging problem in interactive computer graphics. Ideal algorithms should not only . tem is then solved by a global solver such as Conjugate Gradients (CG) or direct methods, and the solution is used to get closer to the . real time simulation results by combining local and global solves.
Physics 20 General College Physics (PHYS 104). Camosun College Physics 20 General Elementary Physics (PHYS 20). Medicine Hat College Physics 20 Physics (ASP 114). NAIT Physics 20 Radiology (Z-HO9 A408). Red River College Physics 20 Physics (PHYS 184). Saskatchewan Polytechnic (SIAST) Physics 20 Physics (PHYS 184). Physics (PHYS 182).
Advanced Placement Physics 1 and Physics 2 are offered at Fredericton High School in a unique configuration over three 90 h courses. (Previously Physics 111, Physics 121 and AP Physics B 120; will now be called Physics 111, Physics 121 and AP Physics 2 120). The content for AP Physics 1 is divided
General Physics: There are two versions of the introductory general physics sequence. Physics 145/146 is intended for students planning no further study in physics. Physics 155/156 is intended for students planning to take upper level physics courses, including physics majors, physics combined majors, 3-2 engineering majors and BBMB majors.
Physics SUMMER 2005 Daniel M. Noval BS, Physics/Engr Physics FALL 2005 Joshua A. Clements BS, Engr Physics WINTER 2006 Benjamin F. Burnett BS, Physics SPRING 2006 Timothy M. Anna BS, Physics Kyle C. Augustson BS, Physics/Computational Physics Attending graduate school at Univer-sity of Colorado, Astrophysics. Connelly S. Barnes HBS .
PHYSICS 249 A Modern Intro to Physics _PIC Physics 248 & Math 234, or consent of instructor; concurrent registration in Physics 307 required. Not open to students who have taken Physics 241; Open to Freshmen. Intended primarily for physics, AMEP, astronomy-physics majors PHYSICS 265 Intro-Medical Ph
strong Ph.D /strong . in Applied Physics strong Ph.D /strong . in Applied Physics with Emphasis on Medical Physics These programs encompass the research areas of Biophysics & Biomedical Physics, Atomic Molecular & Optical Physics, Solid State & Materials Physics, and Medical Physics, in
Modern Physics: Quantum Physics & Relativity. You can’t get to Modern Physics without doing Classical Physics! The fundamental laws and principles of Classical Physics are the basis Modern Physics
Ib physics hl ia. Ib physics hl data booklet. Ib physics hl notes. Ib physics hl topics. Ib physics hl textbook. Ib physics hl past papers. Ib physics hl grade boundaries. If you are watching this program, you are probably thinking of taking IB Economics or are currently enrolled in the
