A fully implicit, exactly conserving algorithm for multidimensional particle-in-cell kinetic simulations L. Chacón Applied Mathematics and Plasma Physics Group Theoretical Division Los Alamos National Laboratory P.O. Box 1663, Los Alamos, NM 87544 Other collaborators: G. Chen, W. Taitano, D. A. Knoll (LANL) D. C. Barnes (Coronado Consulting) Numerical Methods for Large-Scale Nonlinear Problems and Their Applications Aug 31-Sept 4, 2015 ICERM, Brown U. Providence, RI Work funded by DOE-SC ASCR, and the LANL and ORNL LDRD programs Luis Chacon, chacon@lanl.gov
Outline ➤ Particle-in-cell (PIC) methods for plasmas ➤ Explicit vs. implicit PIC ➤ Preconditioned JFNK ➤ 1D ES implicit PIC ➫ Multirate implicit time integrator (particle subcycling) ➫ Moment-based preconditioning ➤ Generalization to Curvilinear Multi-D EM PIC: Vlasov-Darwin model ➫ Review and motivation for Darwin model ➫ Conservation properties (energy, charge, and canonical momenta) ➫ Numerical benchmarks Luis Chacon, chacon@lanl.gov
Introduction Luis Chacon, chacon@lanl.gov
Kinetic Plasma Simulation ➤ A fully ionized collisionless plasma: ions, electrons, and electromagnetic fields ➤ Challenge : integrate electron-ion-field kinetic system on an ion time-scale and a system length scale while retaining electron kinetic effects accurately. (We are developing a new implicit algorithm for long-term, system-scale simulations. ) ➤ Problem features a hierarchical description : ➫ How to design a multiscale algorithm? ➫ How to respect conservation laws, and constraints? Luis Chacon, chacon@lanl.gov
Particle-in-cell (PIC) methods for kinetic plasma simulation ∂ t f + v · ∇ f + F m · ∇ v f = 0 ➤ Lagrangian solution by the method of characteristics : � � t � t � 0 dt v , v − 1 f ( x , v , t ) = f 0 x − ; x ( t = 0 ) = x 0 ; v ( t = 0 ) = v 0 0 dt F m ➤ PIC approach follows characteristics employing macroparticles (volumes in phase space) f ( x , v , t ) = ∑ p δ ( x − x p ) δ ( v − v p ) ∂ t B + ∇ × E = 0 = x p ˙ v p − µ 0 ǫ 0 ∂ t E + ∇ × B = µ 0 j q p ∇ · B = 0 = ( E + v × B ) v p ˙ m p ρ ∇ · E = ǫ 0 → S ( x − x p ) ; E p = ∑ E i S ( x i − x p ) ; j i = ∑ δ ( x − x p ) − j p S ( x i − x p ) p i Luis Chacon, chacon@lanl.gov
State-of-the-art classical PIC algorithm is explicit ➤ Classical explicit PIC: “leap-frogs” particle positions and velocities, field-solve at position update: ➤ Implementation is straightforward, but... ➤ Performance limitations: ➫ CFL-type instability: min( ω pe ∆ t < 1 , c ∆ t < ∆ x ). Minimum temporal resolution ➫ Finite-grid instability: ∆ x < λ Debye . Minimum spatial resolution ➫ Memory bound: challenging for efficient use of modern computer architectures. ➤ Accuracy limitations: ➫ Lack of energy conservation , problematic for long-time-scale simulations ➤ To remove the stability/accuracy constraints of explicit methods, we consider implicit methods. Luis Chacon, chacon@lanl.gov
Implicit PIC methods ➤ Exploration of implicit PIC started in the 1980s ➫ Implicit moment method 1 ➫ Direct implicit method 2 ➤ Early approaches used linearized, semi-implicit formulations: ➫ Lack of nonlinear convergence ➫ Particle orbit accuracy (particle and fields integrated in lock-step) ➫ Inconsistencies between particles and moments ➫ Inaccuracies! → Plasma self-heating/cooling 3 ➤ Our approach: nonlinear implicit PIC ➫ Enforcing nonlinear convergence; complete consistency between particles, moments, and fields. ➫ Allowing stable and robust integrations with large ∆ t and ∆ x ( 2nd order accurate). ➫ Ensuring exact global energy conservation and local charge conservation properties. ➫ Allowing adaptivity in both time and space without loss of the conservation properties. ➫ Allowing particle subcycling → high operational intensities ( compute bound ). ➫ Allowing fluid preconditioning to accelerate the iterative kinetic solver! 1 Mason, R. J. (1981), Brackbill, J. U., and Forslund, D. W. (1982) 2 Friedman, A., Langdon, A. B. and Cohen, B. I.(1981) 3 Cohen, B. I., Langdon, A. B., Hewett, D. W., and Procassini, R. J. (1989) Luis Chacon, chacon@lanl.gov
Fully implicit PIC: 1D electrostatic PIC Chen et al, JCP 2011, 2012, 2013; Taitano et al, SISC (2013) Luis Chacon, chacon@lanl.gov
Fully implicit PIC formulation (at first glance) ➤ A fully implicit formulation couples particles and fields non-trivially (integro-differential PDE): f n + 1 − f n + v · ∇ f n + 1 + f n m ∇ Φ n + 1 + Φ n f n + 1 + f n q − · ∇ v = 0 ∆ t 2 2 2 ∇ 2 Φ n + 1 − Φ n d v v f n + 1 + f n � � ∇ 2 Φ n + 1 = d v f n + 1 = ∇ · or ∆ t 2 ➤ In PIC, f n + 1 is sampled by a large collection of particles in phase space, { x , v } n + 1 . p ➫ There are N p particles, each particle requiring 2 × d equations ( d → dimensions), ➫ Field requires N g equations, one per grid point. ➤ If implemented naively, an impractically large algebraic system of equations results: G ( { x , v } n + 1 , { Φ n + 1 } g ) = 0 → dim ( G ) = 2 dN p + N g p ➫ No current computing mainframe can afford the memory requirements ➫ Algorithmic issues are showstoppers (e.g., how to precondition it?) ➤ An alternative strategy exists: nonlinear elimination ( particle enslavement ) Luis Chacon, chacon@lanl.gov
Particle enslavement (nonlinear elimination) ➤ Full residual G ( { x , v } p , { Φ } g ) = 0 is impractical to implement ➤ Alternative: nonlinearly eliminate particle quantities so that they are not dependent variables : ➫ Formally, particle equations of motion are functionals of the electrostatic potential: x n + 1 = x p [ Φ n + 1 ] ; v n + 1 = v p [ Φ n + 1 ] p p n + 1 , v p n + 1 , Φ n + 1 ) = G ( x [ Φ n + 1 ] , v [ Φ n + 1 ] , Φ n + 1 ) = ˜ G ( Φ n + 1 ) G ( x p Nonlinear residual can be unambiguously formulated in terms of electrostatic potential only! ➤ JFNK storage requirements are dramatically decreased, making it tractable: ➫ Nonlinear solver storage requirements ∝ N g , comparable to a fluid simulation ➫ Particle quantities ⇒ auxiliary variables: only a single copy of particle population needs to be maintained in memory throughout the nonlinear iteration Luis Chacon, chacon@lanl.gov
Energy-conserving (EC) Vlasov-Ampère discretization ➤ Fully implicit Crank-Nicolson time discretization: E n + 1 − E n + j n + 1 / 2 i i − � j � = 0; ε 0 i ∆ t x n + 1 − x n p p − v n + 1 / 2 = 0; p ∆ t v n + 1 − v n − q p p p E n + 1 / 2 S ( x i − x n + 1 / 2 m p ∑ )= 0; p i ∆ t i j n + 1 / 2 q p v n + 1 / 2 S ( x n + 1 / 2 = ∑ − x i ) . p p i p ➤ C-N enforces energy conservation to numerical round-off: E n + 1 E n + 1 − E n + E n m p 1 1 ⇒ ∑ p + ∑ 2 ( v n + 1 p )( v n + 1 ∑ + v n − v n p ) = − ∑ i i i i 2 m p v 2 2 ε 0 E 2 i = const ε 0 p p ∆ t 2 p p i i ➫ No CFL condition. ➫ Does not suffer from finite-grid instabilities: ∆ x ≮ λ D !! ➫ Requires that particles and fields are nonlinearly converged. Luis Chacon, chacon@lanl.gov
Jacobian-Free Newton-Krylov Methods � x n + 1 ) = � ➤ After spatial and temporal discretization ⇒ a large set of nonlinear equations: G ( � 0 ➤ Converging nonlinear couplings requires iteration: Newton-Raphson method: � ∂ � � G x k = − � � G ( � x k ) δ � � ∂ � x � k ➤ Jacobian linear systems result, which require a linear solver ⇒ Krylov subspace methods (GMRES) ➫ Only require matrix-vector products to proceed. ➫ Jacobian-vector product can be computed Jacobian-free ( CRITICAL : no need to form Jacobian matrix): � � ∂ � � y ) − � G ( � x k + ǫ � G ( � x k ) G y = J k � y = lim � ∂ � x ǫ ǫ → 0 k ➫ Krylov methods can be easily preconditioned: P − 1 ∼ J − 1 k k J k P − 1 � x = − G k k P k δ � We will explore suitable preconditioning strategies later in this talk. Luis Chacon, chacon@lanl.gov
Algorithmic implementation details ➤ The nonlinear residual formulation G ( E n + 1 ) based on Vlasov-Ampere formulation is as follows: 1. Input E (given by JFNK iterative method) 2. Move particles (i.e., find x p [ E ] , v p [ E ] by solving equations of motion) (a) Requires inner (local) nonlinear iteration: Picard (not stiff) (b) Can be as complicated as we desire (substepping, adaptivity, etc) 3. Compute moments (current) 4. Form Vlasov-Ampere equation residual 5. return Because dependent variables are fields, we can explore moment-based preconditioning! ➤ Because particle move is performed within function evaluation, we have much freedom. ➫ We can explore improvements in particle mover to ensure long-term accuracy! ◮ Particle substepping and orbit averaging (ensures orbit accuracy and preserves exact energy conservation) ◮ Exact charge conservation strategy (a new charge-conserving particle mover) ◮ Orbit adaptivity (to improve momentum conservation) Luis Chacon, chacon@lanl.gov
Recommend
More recommend
