An Eulerian Vlasov-Fokker-Planck Algorithm for Spherical Implosion Simulations of Inertial Confinement Fusion Capsules
William T. Taitano, Brett D. Keenan, Luis Chacon, Steven E. Anderson, Hans R. Hammer, Andrei N. Simakov
AAn Eulerian Vlasov-Fokker-Planck Algorithm for Spherical ImplosionSimulations of Inertial Confinement Fusion Capsules
W. T. Taitano a, ∗ , B. D. Keenan b , L. Chacón a , S. E. Anderson a , H. R. Hammer d , A. N. Simakov c a T-5 Applied Mathematics and Plasma Physics Group, Theoretical Division Los Alamos National Laboratory, Los Alamos, NM87545 b XCP-6 Plasma Theory and Appellations Group, Computational Physics Division, Los Alamos National Laboratory, Los Alamos,NM 87545 c Theoretical Design Division, Los Alamos National Laboratory, Los Alamos, NM 87545 d T-3 Fluid Dynamics and Solid Mechanics Group, Theoretical Division Los Alamos National Laboratory, Los Alamos, NM 87545
Abstract
We present a numerical algorithm that enables a phase-space adaptive Eulerian Vlasov-Fokker-Planck (VFP)simulation of an inertial confinement fusion (ICF) capsule implosion. The approach relies on extending arecent mass, momentum, and energy conserving phase-space moving-mesh adaptivity strategy to sphericalgeometry. In configuration space, we employ a mesh motion partial differential equation (MMPDE) strategywhile, in velocity space, the mesh is expanded/contracted and shifted with the plasma’s evolving tempera-ture and drift velocity. The mesh motion is dealt with by transforming the underlying VFP equations intoa computational (logical) coordinate, with the resulting inertial terms carefully discretized to ensure con-servation. To deal with the spatial and temporally varying dynamics in a spherically imploding system, wehave developed a novel nonlinear stabilization strategy for MMPDE in the configuration space. The strategyrelies on a nonlinear optimization procedure that optimizes between mesh quality and the volumetric ratechange of the mesh to ensure both accuracy and stability of the solution. Implosions of ICF capsules aredriven by several boundary conditions: 1) an elastic moving wall boundary; 2) a time-dependent MaxwellianDirichlet boundary; and 3) a pressure-driven Lagrangian boundary. Of these, the pressure-driven Lagrangianboundary driver is new to our knowledge. The implementation of our strategy is verified through a set oftest problems, including the Guderley and Van-Dyke implosion problems –the first-ever reported using aVlasov-Fokker-Planck model.
Keywords:
ICF, spherical implosion, nonlinearly-stabilized MMPDE, 1D2V, Moving Grid,Vlasov-Fokker-Planck
PACS:
1. Introduction
In many high-energy-density (HED) laboratory experiments –including inertial confinement fusion (ICF)ones– experimental design employs single-fluid radiation hydrodynamic models (rad-hydro for short), wherethe constituent ions are modeled as a single-average fluid with an additional electron energy equation (fora comprehensive description overview, see Refs. [1, 2]). Contrary to the rad-hydro predictions, however,the National Ignition Facility (NIF) has not achieved ignition of ICF targets to date. Various authors haveproposed the importance of kinetic effects [3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20,21, 22, 23, 24, 25, 26] and their potential role in altering the rad-hydro predictions. Some of these effects ∗ Corresponding author
Email address: [email protected] (W. T. Taitano)
Preprint submitted to Elsevier September 18, 2020 a r X i v : . [ phy s i c s . c o m p - ph ] S e p nclude, but are not limited to: 1) plasma ion stratification due to differential motion of multiple ions; 2)long-mean-free path modification to transport physics (e.g., heat flux and viscosity); and 3) laser-plasmainteractions.The Vlasov-Fokker-Planck (VFP) kinetic equation, coupled with the Maxwell equations is regardedas the first-principles model for describing the dynamic evolution of weakly collisional (i.e., classical)plasma particle distribution function in 6-dimensional phase-space (3-dimensional in configuration and 3-dimensional in velocity). Consequently, there is wide applicability of the VFP-Maxwell equation set tolaboratory (e.g., fusion experiments, space propulsion) and astrophysical systems. However, VFP simula-tions of ICF capsules are challenged by the disparate length and velocity scales supported, which must beefficiently and accurately dealt with. In configuration space, there exist over 9 orders of magnitude of sep-aration between the Debye length ( ∼ − m), and capsule size (10 − m). In velocity space, there existsmore than 3 orders of magnitude separation between the thermal speed, v th = (cid:112) T / m (where T is the tem-perature, and m is the particle mass) of the initially preheated fuel, the compressed heated fuel, cold shell,electrons, and fusion-produced α -particles. Even if one considers fluid electrons (thus analytically removingkinetic electron length and velocity scales), more than 6 orders of magnitude of separation remains in spacebetween the ion mean free path and the capsule size.As far as we are aware, in the literature there are only three codes capable of performing Eulerian andSemi-Lagrangian VFP implosion simulation of ICF capsules: 1) The FPion and FUSE code by the Com-missariat à l’ènergie atomique et aux ènergies alternatives (CEA) group [27, 28, 29, 30, 31, 32, 33, 34, 35,36, 37, 38], which includes α -particle transport and burn physics capabilities; 2) the HIMICO code by theInstitute for Laser Energetics (ILE) group [39, 40, 41, 42], which in addition to burn physics, includes thetreatment of multi-group radiation transport; and 3) the iFP code at the Los Alamos National Laboratory(LANL) [43, 13, 44, 15, 14, 45, 46, 47, 24, 48, 49, 50, 51]. The approach taken by the CEA group relieson a 1D-2V semi-Lagrangian formalism with a periodic remapping of solution in the phase-space to trackheating/cooling of plasmas and implosion of the capsule, while that by the ILE group relies, effectively, ona spherically symmetric (1D), multi-energy group P1 approximation, where the energy space is discretizedin a finite difference fashion while a first order Legendre polynomial expansion is used in the pitch-anglecoordinate. The CEA approach can handle arbitrary distribution functions in phase-space and can thus inves-tigate capsules that are effectively collisionless such as exploding pushers. However, the formalism (whichrelies on remapping) is not conservative. On the other hand, the ILE approach relies on an Eulerian dis-cretization –on a static energy grid– hence conservation symmetries related to the Vlasov spatial advectionoperator are trivially satisfied. However, their approach assumes linear deviations from isotropic plasmas,and is thus limited in the distributions that it can describe. Further, the ILE approach relies on a Lagrangianformalism for the configuration space grid, where the grid is evolved according to the single-fluid hydro-dynamic velocity, and is thus neither able to selectively adapt the grid near regions where strong kineticfeatures may evolve (i.e., the shock front), nor handle multiple species.In this study, we present a multiscale 1D2V phase-space moving-grid kinetic-ion/fluid-electron hybridalgorithm, used in iFP to simulate spherical implosions of ICF capsules. The algorithm addresses the lim-itations of the previous approaches by FPION/FUSE and HIMICO by allowing the description of arbitrarydistribution functions for an arbitrary number of plasma species while enforcing conservation of mass andenergy (note that, global linear momentum conservation is enforced automatically by spherical symmetry),and moving the grid to resolve evolving sharp gradient structures. In the proposed algorithm, we extendto spherical geometry a recent moving phase-space grid strategy developed for Cartesian geometry [50].The hybrid equations are formulated in Cartesian logical coordinates in phase-space, resulting in severaladditional inertial terms –representing mesh motion– that are discretized so as to enforce the underlyingcontinuum conservation symmetries discretely. In velocity space, the grid is normalized and shifted basedon the local and instantaneous plasma thermal speed and drift velocity, while in configuration space a mesh-motion partial differential equation (MMPDE) formalism [52, 53, 54] is used to track evolving macroscopic2oment structures. The MMPDE strategy employed in our previous study, Ref. [50], defines a fixed gridrelaxation time scale to provide a temporal smoothing mechanism for grid evolution. However, in realisticproblems where different parts of the domain evolve at different rates, a constant relaxation time scale oftenleads to disparate cell volume evolution rate that can cause numerical brittleness and accuracy degradation[55]. To address these limitations, we develop the nonlinearly stabilized (NS)-MMPDE algorithm, whichnonlinearly optimizes the MMPDE grid with the local rate of change of volume. The algorithm’s advertisedcapabilities are demonstrated with known Guderley and Van-Dyke asymptotic hydrodynamic solutions aswell as published implosion results. For reproducibility purposes, we also provide details on the numeri-cal implementations of several implosion driving boundary conditions, namely: 1) elastic moving wall; 2)time-dependent Maxwellian Dirichlet; and 3) pressure-driven Lagrangian boundaries.The rest of paper is organized as follows. Section 2 introduces hybrid ion-VFP and fluid-electron equa-tions in a 1D2V spherical geometry. In Sec. 3, we introduce the transformed VFP ion and the fluid-electronequations in terms of the logical phase-space coordinates. In Sec. 4, we discuss our NS-MMPDE algorithm.In Sec. 5, we discuss in detail the numerical implementation of the proposed scheme in the following order:1) the conservative finite differencing discretization of the ion VFP equation; 2) the discretization of thefluid-electron energy equation; 3) the discretization of the NS-MMPDE algorithm ; 4) the various implosiondriving boundary conditions; 5) our adaptive time-stepping strategy; and 6) a brief description of the solverand overall algorithm combining the various components. The numerical performances and properties ofthe new algorithm are demonstrated with various implosion tests of varying degrees of complexity in Sec.6. Finally, we conclude in Sec. 7.
2. Hybrid Vlasov-Fokker-Planck Ion and Fluid Electron Formulation
In this study, we employ an electrostatic hybrid formulation with a Vlasov-Fokker-Planck kinetic de-scription for the ions and a fluid description for electrons. The Vlasov-Fokker-Planck (VFP) equation de-scribes the particle distribution function for ion species α , f α = f α ( x , v , t ) , interacting with the electric field, E , and colliding with β species (which includes both ions and electrons), given as: ∂ t f α + v · ∇ x f α + q α m α E · ∇ v f α = N s + ∑ β C αβ (cid:0) f β , f α (cid:1) . (2.1)Here, m α is the mass, q α is the charge, v is the particle velocity, ∇ x = { ∂ x , ∂ y , ∂ z } (cid:0) ∇ v = (cid:8) ∂ v x , ∂ v y , ∂ v z (cid:9)(cid:1) isthe configuration- (velocity-) space gradient operator, C αβ is the bilinear Fokker-Planck collision operatorof species α colliding with species β , C αβ = Γ αβ ∇ v · (cid:104) ↔ D β · ∇ v f α − A β f α (cid:105) , (2.2) N s is the total number of ion species, Γ αβ = q α q β Λ αβ πε m α , Λ αβ is the Coulomb logarithm, ε is the permittivityof vacuum, ↔ D β and A β are the collisional velocity space tensor diffusion and friction coefficients of species β , and the detailed numerical treatments are discussed in Refs. [56, 44, 47]. Electrons are described byquasi-neutrality, n e = − ∑ N s α q α n α q e , (2.3)ambipolarity, u e = − ∑ N s α q α n α u α q e n e , (2.4)3n conjunction with a temperature equation:32 ∂ t ( n e T e ) + ∇ x · ( u e n e T e ) + ∇ x · Q e − q e n e u e · E = N s ∑ α W e α . (2.5)Here, n e ( n α ) is the electron ( α -ion species) number density, u e ( u α ) is the electron ( α -ion species) driftvelocity, T e is the electron temperature, Q e is the electron heat flux, W e α is the electron-ion energy exchange,and E is computed from the generalized Ohm’s law which includes the friction contributions with ions. Forcompleteness, we provide the expressions for Q e , W e α , and E in (Appendix A). For detailed descriptions,we refer the readers to Ref. [57], and to Ref. [47] for numerical implementations in the context of our hybridVFP ions/fluid electrons.For spherical implosion problems, we adopt a spherically symmetric radial coordinate system in theconfiguration space and an azimuthally symmetric cylindrical velocity coordinate system. The transformedion VFP equation reads: ∂ t ( J S f α ) + ∂ r (cid:0) J S v || f α (cid:1) + q α m α E || ∂ v || ( J S f α ) + r ∂ v · ( (cid:101) a f α ) = J S (cid:34) N s ∑ β C (cid:0) f β , f α (cid:1) + C α e (cid:0) n e , u || , e T e , f α (cid:1)(cid:35) (2.6)and the fluid electron equation reads:32 ∂ t ( J S n e T e ) + ∂ r (cid:0) J S u || , e n e T e (cid:1) + ∂ r (cid:0) J S Q || , e (cid:1) − J S q e n e u || , e E || = J S N s ∑ α W e α . (2.7)Here, J S = r is the Jacobian of transformation going from Cartesian to a spherical coordinate system; v || and v ⊥ are the parallel (to the radial coordinate) and perpendicular particle velocity; r is the radial coordinate; u || , e is the parallel electron bulk velocity; T e ← k B T e is the electron temperature in units of energy and k B isthe Boltzmann constant; Q || , e is the parallel electron heat flux; E || is the parallel electric field; and C α e is thereduced ion-electron collision operator (for details, refer to Ref. [47]). The fourth term on the left-hand-sideof Eq. (2.6) is the inertial term that arises from the coordinate transformation, with (cid:101) a = (cid:2) v ⊥ , − v || v ⊥ (cid:3) T theassociated pseudo-acceleration vector and ∂ v · A = (cid:104) ∂ v || A || , v − ⊥ ∂ v ⊥ ( v ⊥ A ⊥ ) (cid:105) T the velocity-space divergenceoperator on a vector A = (cid:2) A || , A ⊥ (cid:3) T . The details of the transformation are provided in Appendix B.
3. Transformation of the Governing Equations with Moving Phase-Space Grid
To accommodate the disparate spatio-temporal variation in the thermal speed and drift velocity encoun-tered during the course of ICF capsule implosions, we adopt a moving phase-space grid strategy similar tothat proposed in Ref. [50]. For this, we transform the governing equations on a static-computational grid, (cid:0) ξ , (cid:101) v || , (cid:101) v ⊥ (cid:1) , where ξ ∈ [ , ] is the logical coordinate in configuration space with a mapping to the physicalcoordinate given as r ( ξ , t ) = r min + (cid:90) ξ J ξ (cid:0) ξ (cid:48) , t (cid:1) d ξ (cid:48) , (3.1)and (cid:101) v || and (cid:101) v ⊥ are the transformed velocity-space coordinates related to the physical coordinates as: v || = v ∗ α (cid:101) v || + u ∗|| , α (3.2)and v ⊥ = v ∗ α (cid:101) v ⊥ . (3.3)4ere, r min is the inner-radial domain in the configuration space ( r min = J ξ = ∂ ξ r is theJacobian of transformation from physical to logical coordinates in the configuration space, v ∗ α is the nor-malization speed and u ∗|| , α is the transformation velocity. Accordingly, the VFP equation in the transformedcomputational-coordinate system is given as (details shown in Appendix C): ∂ t (cid:16) J S ξ (cid:101) f α (cid:17) + ∂ ξ (cid:104)(cid:0) J S v || − ˙ J r (cid:1) (cid:101) f α (cid:105) − v ∗ − α ∂ (cid:101) v · (cid:110)(cid:2) J S ξ ∂ t v + (cid:0) J r v || − ˙ J r (cid:1) ∂ ξ v (cid:3) (cid:101) f α (cid:111) + rv ∗ − α J ξ ∂ (cid:101) v · (cid:16)(cid:101) a (cid:101) f α (cid:17) − q α m α E || J S ξ v ∗ − α ∂ (cid:101) v || (cid:101) f α = J S ξ (cid:34) N s ∑ β (cid:101) C αβ + (cid:101) C α e (cid:35) . (3.4)Here, (cid:101) f α = ( v ∗ α ) f α is the normalized distribution function, J S ξ = J S J ξ is the composite Jacobian, ˙ J r = ∂ t r is the local volumetric rate of change in configuration space; ∂ (cid:101) v · A = (cid:104) ∂ (cid:101) v || A || , (cid:101) v − ⊥ ∂ (cid:101) v ⊥ ( (cid:101) v ⊥ A ⊥ ) (cid:105) is thevelocity-space divergence operator on a vector, A = (cid:2) A || , A ⊥ (cid:3) T , in the transformed coordinate system; and v = v ∗ α (cid:101) v + u ∗|| , α e || , where e || = [ , ] T is the unit vector along the parallel velocity. For convenience, we definecommonly used velocity moments of the normalized ion distribution function such as the number density, n α = (cid:68) , (cid:101) f α (cid:69) v , (3.5)the drift velocity u α = u || , α e || = (cid:68) v || , (cid:101) f α (cid:69) v (cid:68) , (cid:101) f α (cid:69) v e || , (3.6)and the temperature, T α = m α (cid:68) ( v − u α ) , (cid:101) f α (cid:69) v (cid:68) , (cid:101) f α (cid:69) v , (3.7)where (cid:104) A ( v ) , B ( v ) (cid:105) v = π (cid:82) ∞ d (cid:101) v ⊥ (cid:82) ∞ − ∞ (cid:101) v ⊥ A ( (cid:101) v ) B ( (cid:101) v ) d (cid:101) v || denotes the transformed velocity-space inner prod-uct operation between A and B , and e || = [ , ] is the unit vector along the radial axis. Similarly, the fluidelectron equation is transformed as:32 ∂ t (cid:0) J S ξ n e T e (cid:1) + ∂ ξ (cid:20)(cid:18) J S u || , e −
32 ˙ J r (cid:19) n e T e (cid:21) + ∂ ξ (cid:0) J S Q || , e (cid:1) − J S ξ q e n e u || , e E || = J S ξ N s ∑ α W e α . (3.8)
4. Nonlinearly-Stabilized Moving Mesh Partial Differential Equation (NS-MMPDE)
We evolve the configuration space grid by a predictor-corrector strategy. In the predictor stage, wenumerically solve the MMPDE [52, 53]: ∂ t r ∗ = τ − r ∂ ξ (cid:2) ω∂ ξ r ∗ (cid:3) , (4.1)where r ∗ is the predictor grid location (to be discussed shortly) in configuration space with boundary condi-tions, r ∗ | ξ = = r min = r ∗ | ξ = = R , τ r is the grid-equilibration time-scale. The grid evolves according toa given monitor function, ω = ω ( ξ , t ) , which depends on inverse gradient-length scales of number density, n , drift velocity, u || , thermal speed, v th , and electron temperature, T e , as: ω ∗ = (cid:115) N s ∑ α (cid:20)(cid:0) L − n , α (cid:1) + (cid:0) L − v th , α (cid:1) + (cid:16) L − u || , α (cid:17) (cid:21) + (cid:0) L − T e (cid:1) , (4.2)5here the ( ∗ ) superscript denotes the pre-smoothed (to be discussed shortly) monitor function. The inverse-gradient-length scales, L − , are computed as: L − n , α = (cid:12)(cid:12) δ ξ ln n α (cid:12)(cid:12) + δ min , L − v th , α = (cid:12)(cid:12) δ ξ ln v th , α (cid:12)(cid:12) + δ min , (4.3) L − u || , α = (cid:12)(cid:12)(cid:12)(cid:12) δ ξ u || , α v th , α (cid:12)(cid:12)(cid:12)(cid:12) + δ min , L − T e = (cid:12)(cid:12) δ ξ ln T e (cid:12)(cid:12) + δ min , where δ min is the user-specified floor for the relative variation ( δ min = .
01 unless otherwise specified). Inaddition, to avoid an arbitrarily fine mesh, we limit the minimum-to-maximum ratio of the monitor function, ω ∗ min / ω ∗ max , to a cutoff value, η , by modifying the monitor function by a constant, a , ω ∗ ← ω ∗ + a , (4.4)where a = ( ηω ∗ max − ω ∗ min ) / ( − η ) . Note that this limiting is performed only when ω ∗ min / ω ∗ max ≤ η . Unlessotherwise stated, η = .
01 is used in this study. We also employ a Winslow grid smoothing operation [58]for the monitor function to ensure the regularity of the mesh: (cid:104) − λ ω ∂ ξ ξ (cid:105) ω = ω ∗ . (4.5)Here, ω is the smoothed monitor function, and λ ω is an empirically chosen smoothing coefficient (5 × − in this study).In problems where the boundary deforms (i.e., imploding ICF capsules) we prevent the grid from tan-gling at the boundary by evolving the MMPDE equation in terms of the normalized radial coordinate, ∂ t (cid:98) r ∗ = τ − r ∂ ξ (cid:2) ω ( ξ , t ) ∂ ξ (cid:98) r ∗ (cid:3) , (4.6)where (cid:98) r ∗ ( ξ , t ) = r ∗ ( ξ , t ) / R ( t ) is the normalized radial coordinate with boundary conditions, (cid:98) r ∗ | ξ = = (cid:98) r ∗ | ξ = =
1, and R ( t ) is the instantaneous radial domain size (determined by the implosion boundaryconditions – to be discussed in Sec. 5.4). Since (cid:98) r ∗ ∈ [ , ] at all times, the regularity of the renormalizedgrid, r ∗ i + / = (cid:98) r ∗ i + / R , is trivially ensured.In ICF capsule implosions, where a variety of time-scales exist, it is often not possible to choose a singlevalue for τ r that simultaneously ensures 1) a stable grid variation that does not excite numerical instabilities,and 2) a fast enough grid velocity that tracks important physical features (e.g., shocks). To address theselimitations with the classic MMPDE scheme, we apply a corrector step where the predicted grid obtainedfrom solving Eq. (4.6) is optimized by minimizing the following cost function: F = (cid:90) d ξ ( ln ζ ) + ψ vol − J ( p ) S ξ ζ J ∗ S ξ + − ζ J ∗ S ξ J ( p ) S ξ − λ (cid:20) (cid:90) d ξ ζ J ∗ S ξ − (cid:90) d ξ J ∗ S ξ (cid:21) , (4.7)for the modifier function, ζ ( ξ ) , and the Lagrange multiplier, λ , for the total volume preservation constraint(both discussed shortly). This is achieved by solving the following system of nonlinear equations, (cid:20) ∂ ζ F ∂ λ F (cid:21) = , (4.8)as described later in this study (Sec. 5.3). Here, the function ζ ∈ ( , ∞ ) is the modification function forthe Jacobian of the predicted grid (i.e., J S ξ = ζ J ∗ S ξ ), ψ vol = .
05 is the empirically chosen penalty factor for6he cost function based on the volumetric-rate-of-change of the grid with respect to the previous time step, p , and J ∗ r ξ = ( r ∗ ) ∂ ξ r ∗ is the composite Jacobian for the predicted grid . We remark that the form of thecost function includes ln ζ to ensure the positivity of ζ , and the last term in the Eq. (4.7) is a constraint thatensures that the optimized grid preserves the total volume of the domain determined by the predicted grid(inspired by similar considerations by Ref. [59]), which is dictated by the boundary conditions for r ∗ . Theoptimized grid is finally computed from, r ( ξ ) = r min + (cid:90) ξ ζ (cid:0) ξ (cid:48) (cid:1) J ∗ ξ (cid:0) ξ (cid:48) (cid:1) d ξ . (cid:48) (4.9)
5. Numerical Implementation
We employ a conservative finite-difference discretization in the 1D2V phase space, where the phase-space cell volume is defined on the grid (cid:0) ξ i , (cid:101) v || , j , (cid:101) v ⊥ , k (cid:1) as: ∆ V i , j , k = π J S ξ , i (cid:101) v ⊥ , k ∆ (cid:101) v || ∆ (cid:101) v ⊥ , (5.1)with, ∆ (cid:101) v || = (cid:101) L || N || (cid:16) ∆ v ⊥ = (cid:101) L ⊥ N ⊥ (cid:17) the parallel- (perpendicular-) velocity cell size, (cid:101) L || ( (cid:101) L ⊥ ) is the transformedparallel- (perpendicular-) velocity domain size, N || ( N ⊥ ) is the number of cells in the parallel (perpendic-ular) velocity space, J S ξ , i = J S , i J ξ , i , J S , i = r i , and J ξ , i = ∆ r i ∆ ξ . Here, ∆ r i = r i + / − r i − / , is the discreteconfiguration-space cell size, and r i = (cid:118)(cid:117)(cid:117)(cid:116) r i + / − r i − / r i + / − r i − / (5.2)is the cell-center location defined such that the discrete-volume integral preserves the exact volume of asphere, i.e., r i ∆ r i = (cid:82) r i + / r i − / r dr = r i + / − r i − / . We define the discrete velocity-space moment of function B as, (cid:104) A ( v ) , B ( v ) (cid:105) v ≈ (cid:104) A ( v ) , B ( v ) (cid:105) δ v = π N || ∑ j = ∆ (cid:101) v || N ⊥ ∑ k = (cid:101) v ⊥ , k ∆ (cid:101) v ⊥ A j , k B j , k . (5.3)Additionally, with the distribution function defined on cell centers and fluxes defined on cell faces, thediscretized VFP equation in the transformed coordinate system reads: (cid:104) δ t (cid:16) J r ξ (cid:101) f α (cid:17)(cid:105) ( p + ) i , j , k + (cid:2) δ ξ F r , α (cid:3) ( p + ) i , j , k (cid:124) (cid:123)(cid:122) (cid:125) a (cid:13) + (cid:2) δ ξ F ˙ r , α (cid:3) ( p + ) i , j , k (cid:124) (cid:123)(cid:122) (cid:125) b (cid:13) + (cid:104) δ (cid:101) v || J acc , α (cid:105) ( p + ) i , j , k (cid:124) (cid:123)(cid:122) (cid:125) c (cid:13) + [ δ (cid:101) v · ( γ t , α J t , α )] ( p + ) i , j , k (cid:124) (cid:123)(cid:122) (cid:125) d (cid:13) + (cid:2) δ (cid:101) v · (cid:0) γ r , α , i + / J − r , α + γ r , α , i − / J + r , α (cid:1)(cid:3) ( p + ) i , j , k (cid:124) (cid:123)(cid:122) (cid:125) e (cid:13) +[ δ (cid:101) v · ( γ S , α J S , α )] ( p + ) i , j , k (cid:124) (cid:123)(cid:122) (cid:125) f (cid:13) = J ( p + ) r ξ , i (cid:34) N s ∑ β (cid:101) C αβ + (cid:101) C α e (cid:35) ( p + ) i , j , k . (5.4)Here, we introduce a shorthand notation for discrete-time derivatives,7 δ t (cid:16) J S ξ (cid:101) f α (cid:17)(cid:105) ( p + ) i , j , k = c ( p + ) J ( p + ) S ξ , i (cid:101) f ( p + ) α , i , j , k + c ( p ) J ( p ) S ξ , i (cid:101) f ( p ) α , i , j , k + c ( p − ) J ( p − ) S ξ , i (cid:101) f ( p − ) α , i , j , k ∆ t ( p ) , (5.5)where c ( p + ) , c ( p ) , and c ( p − ) are the coefficients for the second-order backward-differencing scheme (BDF2)[60] and superscripts ( p ) indicate the discrete time index. The discrete spatial divergence operator is definedas, (cid:2) δ ξ F (cid:3) ( p + ) i , j , k = F ( p + ) i + / , j , k − F ( p + ) i − / , j , k ∆ ξ , (5.6)and velocity-space divergence operator as, (cid:104) δ (cid:101) v · J ( p + ) (cid:105) i , j , k = (cid:104) δ (cid:101) v || J || + (cid:101) v − ⊥ δ (cid:101) v ⊥ ( (cid:101) v ⊥ J ⊥ ) (cid:105) ( p + ) i , j , k = J ( p + ) || , i , j + / , k − J ( p + ) || , i , j − / , k ∆ (cid:101) v || + (cid:101) v ⊥ , k + / J ( p + ) ⊥ , i , j , k + / − (cid:101) v ⊥ , k − / J ( p + ) ⊥ , i , j , k − / (cid:101) v ⊥ , k ∆ (cid:101) v ⊥ . (5.7)Here, the subscripts i ± / j ± / k + / (cid:13) cor-responds to the discrete representation of the spatial streaming term, b (cid:13) corresponds to the inertial termin the configuration space (arising from the moving grid), c (cid:13) corresponds to the electrostatic-accelerationterm, d (cid:13) corresponds to the velocity-space inertial terms due to temporal variation of the velocity-spacemetrics (i.e., v ∗ α and (cid:98) u ∗|| , α ), e (cid:13) corresponds to the velocity-space inertial terms due to the spatial variationof the metrics, and f (cid:13) corresponds to the inertial term arising from transformation from Cartesian phase-space to a spherical-configuration and cylindrical-velocity-space coordinate system. Additionally, γ t ( ξ , (cid:101) v , t ) , γ r ( ξ , (cid:101) v , t ) , and γ S ( ξ , (cid:101) v , t ) in terms d (cid:13) , e (cid:13) , and f (cid:13) are discrete-nonlinear-constraint functions, which enforcethe underlying integral-discrete-conservation symmetries, and therefore the discrete conservation of energyof the system. The functional forms of γ t , γ r , and γ S are discussed in Appendix D. For reproducibil-ity purposes, the detailed treatment of fluxes for individual terms are discussed in Appendix E. Finally,the right-hand-side of Eq. (5.4) corresponds to the discrete collision operator, for which we follow Refs.[43, 44, 61]. The electron temperature equation in the transformed coordinate system, Eq. (3.8), is discretized usinga finite-difference scheme in space and BDF2 in time:32 (cid:2) δ t (cid:0) J S ξ n e T e (cid:1)(cid:3) ( p + ) i + (cid:104) δ ξ (cid:16) J S u || , e (cid:100) n e T e (cid:17)(cid:105) ( p + ) i − (cid:104) δ ξ (cid:16) ˙ J r (cid:100) n e T e (cid:17)(cid:105) ( p + ) i + (cid:2) δ ξ (cid:0) J S Q || , e (cid:1)(cid:3) ( p + ) i − J ( p + ) S ξ , i q e (cid:2) n e u || , e E || (cid:3) ( p + ) i = J ( p + ) S ξ , i N s ∑ α W ( p + ) e α , i . (5.8)Here, the hat denotes cell-face interpolated quantities, for which we use a SMART discretization [62].Following Ref. [47], we discretize the Joule heating term (last term on the left hand side) to ensure exactenergy balance with the ion electrostatic acceleration term. The modified MMPDE equation, Eq. (4.6), is discretized on a staggered grid as: c ( p + ) (cid:98) r ∗ i + / + c ( p ) (cid:98) r ( p ) i + / + c ( p − ) (cid:98) r ( p − ) i + / ∆ t ( p ) − τ r ω ( p ) i + (cid:16)(cid:98) r ∗ i + / − (cid:98) r ∗ i + / (cid:17) − ω ( p ) i (cid:16)(cid:98) r ∗ i + / − (cid:98) r ∗ i − / (cid:17) ∆ ξ = . (5.9)8o reduce the nonlinearity between the solution and the grid for a given time step, the monitor function isevaluated from the previous time solution. The cost function for the MMPDE optimization in Eq. (4.7) isdiscretized as: F = N ξ ∑ i = ∆ ξ ( ln ζ i ) + ψ vol − J ( p ) S ξ , i ζ i J ∗ S ξ , i + − ζ i J ∗ S ξ , i J ( p ) S ξ , i − λ (cid:34) N ξ ∑ i = ∆ ξ ζ i J ∗ S ξ , i − N ξ ∑ i = ∆ ξ J ∗ S ξ , i (cid:35) . (5.10)Here, J ∗ r ξ , i = ( r ∗ i ) (cid:104) r ∗ i + / − r ∗ i − / (cid:105) ∆ ξ is the discrete composite Jacobian for the predicted grid and the cell-centergrid location is defined by Eq. (5.2). The discrete cost function is minimized with respect to ζ i and λ bysolving : ∂ ζ F ... ∂ ζ N ξ F ∂ λ F = ∆ ξ (cid:26) ln ζ ζ + ψ vol (cid:20)(cid:18) − J ( p ) S ξ , ζ J ∗ S ξ , (cid:19) J ( p ) S ξ , ζ J ∗ S ξ , − (cid:18) − ζ J ∗ S ξ , J ( p ) S ξ , (cid:19) J ∗ S ξ , J ( p ) S ξ , (cid:21) − λ J ∗ S ξ , (cid:27) ... ∆ ξ (cid:40) ln ζ N ξ ζ N ξ + ψ vol (cid:34)(cid:32) − J ( p ) S ξ , N ξ ζ N ξ J ∗ S ξ , N ξ (cid:33) J ( p ) S ξ , N ξ ζ N ξ J ∗ S ξ , N ξ − (cid:32) − ζ N ξ J ∗ S ξ , N ξ J ( p ) S ξ , N ξ (cid:33) J ∗ S ξ , N ξ J ( p ) S ξ , N ξ (cid:35) − λ J ∗ S ξ , N ξ (cid:41) − ∑ N ξ i = ∆ ξ ζ i J ∗ S ξ , i + ∑ N ξ i = ∆ ξ J ∗ S ξ , i = , (5.11)which is a nonlinearly coupled system of equations and must be solved iteratively. We employ a Newtonsolver, where the solution is iterated as: x ( k + ) = x ( k ) + δ x ( k ) , (5.12)where, the superscript, k , denotes the Newton iteration index, x = (cid:104) ζζζ T , λ (cid:105) T is the solution vector, δ x ( k ) = − J ( k ) , − R ( k ) (5.13)is the Newton update, J = ∂ R ∂ x is the Jacobian matrix, and R = (cid:2) ∂ ζζζ F T , ∂ λ F (cid:3) T is the nonlinear residualvector. Due to the trivial 1D nature of the problem Eq. (5.13) is solved using a direct inversion of J ( k ) . Theiteration is continued until | R ( k ) | ≤ ε grid | R ( ) | , where ε grid = − is the nonlinear convergence tolerance.Finally, the new-time grid is updated as: r ( p + ) i + / = r ( p + ) i − / + i ∑ j = ζ j J ∗ ξ , j . (5.14)The procedure for the NS-MMPDE algorithm is given in Alg. 1. We note that, relative to the overallnonlinear solve of the coupled hybrid equations (to be discussed in Sec. 5.6), the cost of the nonlinear solvefor Eq. (5.11) is negligible.Once the new-time grid is obtained, the volumetric rate of change of the cell ˙ J r (which is related to thegrid speed as ∂ t r = ˙ r = ˙ J r / J S ), is obtained in such a way as to discretely ensure the geometric conservationlaw, ∂ t J S ξ = ∂ ξ ˙ J r . (5.15)Discretizing and using Eq. (5.2), we find: 9 lgorithm 1 Grid update procedure with the NS-MMPDE algorithm.1. Compute the non-smoothed monitor function, ω ∗ , ( p ) , based on the previous time solution.2. Limit the monitor function via ω ∗ ← ω ∗ + ( ηω ∗ max − ω ∗ min ) / ( − η ) if ω ∗ min / ω ∗ max ≤ η .3. Smooth the monitor function by solving (cid:104) − λ ω ∂ ξ ξ (cid:105) ω = ω ∗ where λ ω = × − in this study.4. Solve the normalized predictor grid, (cid:98) r ∗ i + / , from Eq. (5.9).5. Compute the renormalized predictor grid, r ∗ i + / = (cid:98) r ∗ i + / R ( p + ) , where R ( p + ) is the instantaneous domainsize, determined from implosion-driving boundary conditions (to be discussed in Sec 5.4).6. Compute the predictor Jacobian, J ∗ S ξ , i = ( r ∗ i ) (cid:104) r ∗ i + / − r ∗ i − / (cid:105) ∆ ξ .7. Solve for ζ and λ from Eqs. (5.11)-(5.13).8. Update the new time grid from Eq. (5.14). c ( p + ) (cid:20)(cid:16) r ( p + ) i + / (cid:17) − (cid:16) r ( p + ) i − / (cid:17) (cid:21) + c ( p ) (cid:20)(cid:16) r ( p ) i + / (cid:17) − (cid:16) r ( p ) i − / (cid:17) (cid:21) + c ( p − ) (cid:20)(cid:16) r ( p − ) i + / (cid:17) − (cid:16) r ( p − ) i − / (cid:17) (cid:21) ∆ t ( p ) ∆ ξ − ˙ J ( p + ) r , i + / − ˙ J ( p + ) r , i − / ∆ ξ = . (5.16)Equating the terms with common cell indices gives the definition of ˙ J r as:˙ J ( p + ) r , i + / = c ( p + ) (cid:16) r ( p + ) i + / (cid:17) + c ( p ) (cid:16) r ( p ) i + / (cid:17) + c ( p − ) (cid:16) r ( p − ) i + / (cid:17) ∆ t ( p ) . (5.17) We discuss next the various ICF implosion-driving boundary conditions. Unless otherwise specified, forthe spatial advection terms, we employ an upwind discretization at the configuration-space boundary face tocompute the flux, i.e., F B , j , k = v || , j (cid:40) (cid:101) f B − / , j , k if v || , j ≥ (cid:101) f B + / , j , k otherwise (cid:41) . Here, the subscript B denotes the boundary surface and B ± / To model an infinitely massive piston compressing a plasma, we consider the elastic scattering of acollection of particles off a moving wall. Consider a particle with velocity, v , reflecting off a moving wallwith velocity, u W . One can derive the exact kinematic relation between the pre- and post-scattered particlevelocities, ( v , v s ) by transforming the velocity in the reference of the wall as: v (cid:48) = v − u W , (5.18)and v (cid:48) s = v s − u W . (5.19)In the frame of reference of the moving wall, an elastic scatter will reflect the particle’s parallel (to wallnormal vector) velocity component, v (cid:48)|| = v (cid:48) · n W = v || − u W , || , (5.20)10n the equal and opposite direction while the perpendicular/tangent component, v (cid:48)⊥ = v (cid:48) − v (cid:48)|| , (5.21)is unchanged, i.e., v (cid:48)|| , s = − v (cid:48)|| , (5.22) v (cid:48)⊥ , s = v (cid:48)⊥ . (5.23)Here, n W is the wall normal vector (in 1D spherically symmetric geometry, simply -1/1 at inner/outer radialboundary). Finally, substituting Eq. (5.22) into Eq. (5.19) and solving for v s (i.e., in the laboratory frame),we obtain for the scattered-particle parallel velocity as v || , s = − v || + u W , || . Given a distribution function of particles, this kinematic relationship can be expressed equivalently in themoving frame as: (cid:101) f (cid:16) v (cid:48)|| (cid:17)(cid:12)(cid:12)(cid:12) B = (cid:101) f (cid:16) − v (cid:48)|| (cid:17)(cid:12)(cid:12)(cid:12) B if v (cid:48)|| n W , || ≤ (cid:101) f (cid:16) v (cid:48)|| (cid:17)(cid:12)(cid:12)(cid:12) B else . Numerically, the boundary condition in the laboratory frame is ensured accordingly for the distributionfunction: (cid:101) f (cid:0) ξ B , v || (cid:1) = (cid:40) C (cid:101) f (cid:0) ξ B ± / , − v || + u W , || (cid:1) if (cid:0) v || − u W , || (cid:1) n W , || ≤ (cid:101) f (cid:0) ξ B ± / , v || (cid:1) else , (5.24)where C is the scaling factor that ensures discrete mass conservation at the boundary (discussed shortly).We note that, in general, it is not possible to reconstruct numerically an exact one-to-one map between thepre- and post-scattered particle distribution function due to the discrete nature of velocity-space grid (i.e.,there may be no grid point at the specified reflection velocity). For numerical purposes, we reconstruct thescattered distribution function by splining the incoming velocity distribution - in this study, we consider asecond-order spline to ensure monotonicity of the reflected distribution function. Further, because of thediscrete violation of symmetry in the reconstruction procedure, the mass is not automatically conserved fora finite-wall velocity. In order to discretely conserve mass, we rescale the reflected distribution function toensure the following constraint: (cid:68) , (cid:16) J S v || (cid:101) f − ˙ J r (cid:101) f (cid:17) B (cid:69) δ v = C (cid:68) , (cid:16) J S v || (cid:101) f − ˙ J r (cid:101) f (cid:17) B (cid:69) v || ≥ u W , || δ v + (cid:68) , (cid:16) J S v || (cid:101) f − ˙ J r (cid:101) f (cid:17) B (cid:69) v || < u W , || δ v , (5.25)where J r v || (cid:101) f − ˙ J r (cid:101) f is the effective configuration-space flux including the boundary mesh motion, and: (cid:104) , ( · · · ) (cid:105) v || ≥ u W , || δ v = π N ⊥ ∑ k = N || ∑ j = j (cid:48) (cid:101) v ⊥ , k ∆ ˜ v ⊥ ∆ ˜ v || ( · · · ) j , k , (cid:104) , ( · · · ) (cid:105) v || < u W , || δ v = π N ⊥ ∑ k = j (cid:48) − ∑ j = (cid:101) v ⊥ , k ∆ ˜ v ⊥ ∆ ˜ v || ( · · · ) j , k , j (cid:48) corresponds to the discrete parallel velocity grid point nearest to the 2 u W , || n W boundary, j (cid:48) = floor (cid:16) u W , || − u ∗|| − v ∗ (cid:101) v || , min v ∗ ∆ (cid:101) v || (cid:17) if 1 ≤ floor (cid:16) u W , || − u ∗|| − v ∗ (cid:101) v || , min v ∗ ∆ (cid:101) v || (cid:17) ≤ N || (cid:16) u W , || − u ∗|| − v ∗ (cid:101) v || , min v ∗ ∆ (cid:101) v || (cid:17) < N || otherwise . We note that the symmetry boundary condition –identical to the specular reflection boundary condition considered in Ref. [63]– is a special case of the treatment above in the limit of vanishing wall velocity andis used for the r = u W B
In a full ICF capsule implosion, kinetic plasma physics is not sufficient for fidelity. For a predictivecapability, additional physics are required such as radiation transport, α particle generation and transport,laser ray tracing, as well as non-ideal equations of state for cryogenic and high-Z materials. However, incertain experiments (at the Omega facility for example), these effects become important only in regionsaway from the fuel-pusher interface. Thus, in these instances, it is possible to model only part of the capsuleby decoupling the outer pusher dynamics from the fuel region. In the literature, Ref. [33] was amongst thefirst to explore this coupling with rad-hydro codes such as XRAGE [64], HYDRA [65], and LILAC [66].In their approach, the authors used the rad-hydro simulations to drive the kinetic simulation by sourcing the12oundary condition with a time-dependent Maxwellian distribution function, in which the state variables( n , u , T ) were provided by rad-hydro simulations as a function of time. The particle distribution function atthis VFP-rad-hydro interface is defined as, (cid:101) f B = (cid:40) (cid:101) f M , H (cid:0) n H , B ( t ) , u || , H ( t ) , T H ( t ) , v (cid:1) for v · n B ≤ (cid:101) f B − / ( v ) otherwise , (5.26)where, the subscript H denotes rad-hydro simulation quantities.In a Lagrangian formulation of the hydrodynamic equations, grid velocities are chosen to be the center-of-mass velocity and consequently, lead to a local mass conservation theorem of the associated computa-tional cell. In order to ensure that the boundary conserves total mass, we evolve the boundary velocity toensure a zero total mass flux. A total mass conservation theorem can be derived by taking the m α v momentof Eq. (3.4), integrating over ξ , and summing over all species, N s ∑ α (cid:90) d ξ (cid:2) ∂ t (cid:0) J S ξ ρ α + ∂ ξ (cid:2) J S ρ u || , α − ˙ J r ρ α (cid:3)(cid:1)(cid:3) = , (5.27)where ρ α = m α (cid:68) , (cid:101) f α (cid:69) δ v and ρ u || , α = m α (cid:68) v || , (cid:101) f α (cid:69) δ v . We can reduce the above equation to: ∂ t M + (cid:2) J S ρ u || − ˙ J r ρ (cid:3) ξ = ξ = = , (5.28)where ∂ t M = M = ∑ N s α (cid:82) d ξ J S ξ ρ α , ρ = ∑ N s α ρ α , and ρ u || = ∑ N s α ρ u || , α . We note that, with r ( ξ = , t ) = (cid:2) J S ρ u || − ˙ J r ρ (cid:3) ξ = = J r ρ (cid:12)(cid:12) ξ = = J S ρ u || (cid:12)(cid:12) ξ = . By expanding the relationship and semi-discretizing, we obtain the following dis-crete constraint: δ t r ξ = (cid:28) , (cid:101) f α (cid:12)(cid:12)(cid:12) ξ = (cid:29) δ v = J S , ξ = (cid:28) , v || (cid:101) f α (cid:12)(cid:12)(cid:12) ξ = (cid:29) δ v . (5.29)Since δ t r ξ = = (cid:20) c ( p + ) (cid:16) r ( p + ) ξ = (cid:17) + c ( p ) (cid:16) r ( p ) ξ = (cid:17) + c ( p − ) (cid:16) r ( p − ) ξ = (cid:17) (cid:21) ∆ t ( p ) , this provides us with a relationship to solve the new-time radius, r ( p + ) ξ = that ensures total mass conservation in the system. Arguably, the more physically consistent approach to couple the hydrodynamic simulations with the ki-netic approach is to drive the implosion solely with the external pressure provided by rad-hydro simulationsand self-consistently account for internal (back) pressure to evolve the boundary location. To achieve thisgoal, we introduce an auxiliary equation that evolves the boundary velocity, du || , B dt = − ρ B (cid:2) ∂ r P B − r − B ∂ r (cid:0) r B τ |||| , B (cid:1)(cid:3) . (5.30)Discretizing the above equation, we obtain: δ t u ( p + ) || , B = − ρ ( p + ) B (cid:34) P ( p + ) B + / − P ( p + ) B − / ∆ r ( p + ) B − (cid:16) r ( p + ) B (cid:17) (cid:16) r ( p + ) B , i + / (cid:17) τ ( p + ) |||| , B + / − (cid:16) r ( p + ) B , i − / (cid:17) τ ( p + ) |||| , B − / ∆ r ( p + ) B (cid:35) , (5.31)13here the boundary total-mass density, ρ B , is obtained as an average of the inner and ghost-cell boundaryvalues, ρ ( p + ) B = (cid:104) ρ ( p + ) B + / + ρ ( p + ) B − / (cid:105) , (5.32)where ρ ( p + ) = N s ∑ α (cid:68) m α , (cid:101) f ( p + ) α (cid:69) δ v . (5.33)The total hydrodynamic pressure, P , is computed as the sum of the ion and electron pressures, P ( p + ) = P ( p + ) e + N s ∑ α m α (cid:28)(cid:16) v − u ( p + ) α (cid:17) , (cid:101) f ( p + ) α (cid:29) δ v (5.34)and the total radial viscous stress, τ |||| , is computed as: τ |||| = N s ∑ α ( P ⊥ , α − P || , α ) , (5.35)where the parallel and perpendicular pressures are defined, respectively, as: P || , α = m α (cid:68)(cid:0) v || − u || , α (cid:1) , (cid:101) f α (cid:69) δ v , (5.36)and P ⊥ , α = m α (cid:68) v ⊥ , (cid:101) f α (cid:69) δ v . (5.37)Once u || , B is evaluated, we source a Maxwellian distribution at the boundary ghost cell with the prescribednumber density, temperature, and u || , B . To ensure total mass conservation, we follow the procedure in Eqs.(5.27)-(5.29) to evolve the boundary grid such that the total mass flux at the boundary vanishes. Time-scales in ICF implosions range from the stiff collision time-scale, particle advection and accelera-tion time-scale, dynamical (hydrodynamic) time-scale (in which slow physical features evolve), all the wayto the system-time scale (e.g., implosion time). In order to efficiently and accurately integrate through thevarious time-scales, we employ an adaptive time-stepping strategy where the step size for a given time step, p , is computed as: ∆ t ( p ) = min (cid:110) ∆ t ( ∗ ) , ∆ t max , τ ( p ) v , τ ( p ) ˙ r , τ ( p ) vol , τ ( p ) M (cid:111) , (5.38)where ∆ t ( ∗ ) = γ ∆ t ( p − ) (5.39)is the predictor time-step size based on the previous time-step size, with γ = .
15, and ∆ t max is the maximumallowable time-step size. The advection time scale is given as τ ( p ) v = . (cid:110) τ ( p ) v , n α , τ ( p ) v , u α , τ ( p ) v , T α (cid:111) for α = { , · · · , N s } , (5.40)with τ ( p ) v , M α = L ( p ) M α (cid:12)(cid:12)(cid:12) u ( p ) α (cid:12)(cid:12)(cid:12) + φ v ( p ) th , α , (5.41)14nd L ( p ) M α = (cid:12)(cid:12)(cid:12) ∂ r ln M ( p ) α (cid:12)(cid:12)(cid:12) − , (5.42)the gradient length-scale for moment quantity, M , and φ = . τ ( p ) ˙ r = . (cid:110) τ ( p ) ˙ r , n α , τ ( p ) ˙ r , u α , τ ( p ) ˙ r , T α (cid:111) for α = { , · · · , N s } , (5.43)where τ ( p ) ˙ r , M α = L ( p ) M α (cid:12)(cid:12) ˙ r ( p ) (cid:12)(cid:12) + − . (5.44)We remark that, even with the configuration-space grid nonlinear stabilization strategy employed here, thereare certain instances in realistic implosions where the local volume rate of change of the grid remains large.Time-step-size limiting is used to address these situations. The time scale in which the cell volume varies isgiven as, τ ( p ) vol = . (cid:12)(cid:12)(cid:12) ∂ t ln J ( p ) r ξ (cid:12)(cid:12)(cid:12) − . (5.45)Finally, the exponentiation time-scale for moments is given as, τ ( p ) M = . (cid:110) τ ( p ) n α , τ ( p ) u α , τ ( p ) T α (cid:111) for α = { , · · · , N s } , (5.46) τ ( p ) M α = (cid:12)(cid:12)(cid:12) ∂ t ln M ( p ) α (cid:12)(cid:12)(cid:12) − . (5.47)In this manner, the time-step size will adapt to the dynamically varying features of the physics and the grid. The discretized set of ion-VFP and fluid-electron equations lead to a system of nonlinear equations,which we solve implicitly using an Anderson acceleration scheme [67, 68]. Consider a fixed point map, x l + = G (cid:16) x l (cid:17) , (5.48)which maps the solution from one iteration to another. Here x l = (cid:104)(cid:101) f l , T le (cid:105) T (5.49)is the solution vector at iteration l and and (cid:101) f l − = (cid:104) (cid:101) f l − , · · · , (cid:101) f l − N s (cid:105) . The Anderson acceleration algorithmaccelerates the iterative convergence of the solution by employing the history of previous iterations as: x l + = m l ∑ i = θ li G (cid:16) x l − m l + i (cid:17) , (5.50)where m l = θ li are found min-imizing (cid:13)(cid:13) ∑ m l i = θ li (cid:0) G (cid:0) x l − m l + i (cid:1) − x l − m l + i (cid:1)(cid:13)(cid:13) , subject to the constraint ∑ m l i = θ li =
1. For G (cid:0) x l (cid:1) , we consider a15icard linearized solver where we: 1) solve the fluid electron equation with given ion distribution functions;and 2) solve the ion VFP equations with given fluid electron temperature and the collisional transport co-efficients from the l − T l , s + e = G T e (cid:16) T l , se (cid:17) = T l , se − (cid:16) P l , sT e (cid:17) − R l , sT e . (5.51)Here, P l , sT e ◦ = ∂ t (cid:16) J l − S ξ n l − e ◦ (cid:17) + ∂ ξ (cid:20)(cid:18) J l − S u l − || , e −
32 ˙ J r (cid:19) n l − e ◦ (cid:21) + ∂ ξ (cid:104) J l − S Q l − e ◦ (cid:105) − J l − S ξ q e n l − e u l − || , e ∂ ξ (cid:16) n l − e ◦ (cid:17) − J l − S ξ N s ∑ α = W l − e α ◦ , (5.52)is the preconditioning operator (where ◦ simply denotes the quantity that is being operated on) for theelectron temperature equation subsystem, and R l , sT e ( (cid:101) f l − , T l , se ) = ∂ t (cid:16) J l − S ξ n l − e T l , se (cid:17) + ∂ ξ (cid:20)(cid:18) J l − S u l − || , e −
32 ˙ J r (cid:19) n l − e T l , se (cid:21) + ∂ ξ (cid:104) J l − S Q || , e (cid:16)(cid:101) f l − , T l , se (cid:17)(cid:105) − J l − S ξ q e n l − e u l − || , e E || (cid:16)(cid:101) f l − , T l , se (cid:17) − J l − S ξ N s ∑ α = W e α (cid:16)(cid:101) f l − , T l , se (cid:17) (5.53)is the nonlinear residual we seek the root of. We remind the readers that n l − e = − ∑ N s α (cid:68) , (cid:101) f l − α (cid:69) δ v / q e and u l − || , e = − ∑ N s α q α (cid:68) v || , (cid:101) f l − α (cid:69) δ v / q e n l − e are the electron number density and parallel drift velocity, respec-tively, computed from quasineutrality and ambipolarity conditions; Q l − e and W l − e α in Eq. (5.52) are thelinearized operators for the electron heat flux and the ion-electron energy exchange, respectively, and forthe fluid electron subsystem we use a nonlinear history length of m s =
5. The Quasi-Newton system in Eq.(5.51) is solved using a direct inversion of a resulting tridiagonal system for P l , sT e . The nonlinear iteration iscontinued until (cid:13)(cid:13)(cid:13) R l , sT e (cid:13)(cid:13)(cid:13) ≤ − (cid:13)(cid:13)(cid:13) R l , T e (cid:13)(cid:13)(cid:13) .The ion VFP equation inner solve is similarly performed with an Anderson acceleration with a Quasi-Newton fixed point map, (cid:101) f l , s + α = G (cid:101) f α (cid:16) (cid:101) f l , s α (cid:17) = (cid:101) f l , s α − (cid:16) P l − V FP , α (cid:17) − R l , s α , (5.54)where (cid:16) P l − V FP , α (cid:17) − ◦ = (cid:16) P l − t + P l − ξ , α (cid:17) − (cid:20) I − P l − (cid:101) v (cid:16) P l − t + P l − (cid:101) v (cid:17) − (cid:21) ◦ (5.55)is the operator-split preconditioner, I is the identity operator, P l − t ◦ = ∂ t (cid:16) J l − S ξ ◦ (cid:17) (5.56)is the temporal derivative, P l − ξ , α ◦ = ∂ ξ (cid:104)(cid:16) J l − S v || − ˙ J l − r (cid:17) ◦ (cid:105) (5.57)is the configuration space operator, P l − (cid:101) v , α ◦ = − v ∗ − α ∂ (cid:101) v · (cid:110)(cid:104) J l − S ξ ∂ t v + ∂ ξ (cid:16)(cid:104) J l − S v || − ˙ J l − r (cid:105) v (cid:17)(cid:105) ◦ (cid:111) + r l − v ∗ − α J l − ξ ∂ (cid:101) v · ( (cid:101) a ◦ ) − q α m α E || (cid:16) n l − e , T le (cid:17) J l − S ξ v ∗ − α ∂ (cid:101) v || ◦ − J l − S ξ (cid:34) N s ∑ β (cid:101) C αβ (cid:16) (cid:101) f l − β , (cid:101) f l − α (cid:17) + (cid:101) C α e (cid:16) T le , n l − e , u l − || , e (cid:17)(cid:35) ◦ (5.58)16 lgorithm 2 Evaluation of the outer fixed point map, G (cid:0) x l − (cid:1) , in Eq. (5.48).1. Compute the new boundary location, r l − ξ = , and grid locations, r l − ( (cid:101) f l − , T l − e ), and derived quantities( J l − S , J l − ξ , J l − S ξ )2. Initialize s = T l , e = T l − e
3. Invert T e system for T le
4. For α = N s Invert α -ion species’ VFP system for (cid:101) f l α end5. Set x l = (cid:104)(cid:101) f l , T le (cid:105) T is the velocity-space operator, and R l , s α ( (cid:101) f l , s , T le ) = ∂ t (cid:16) J l − S ξ (cid:101) f l , s α (cid:17) + ∂ ξ (cid:104)(cid:16) J l − S v || − ˙ J l − r (cid:17) (cid:101) f l , s α (cid:105) − v ∗ − α ∂ (cid:101) v · (cid:110)(cid:104) γ l − t , α J l − S ξ ∂ t v + γ l − x , α ∂ ξ (cid:16)(cid:104) J l − S v || − ˙ J l − r (cid:105) v (cid:17)(cid:105) (cid:101) f l , s α (cid:111) + r l − J l − ξ v ∗ − α ∂ (cid:101) v · (cid:16) γ l − S , α (cid:101) a (cid:101) f l , s α (cid:17) − q α m α E l , s || J l − S ξ v ∗ − α ∂ (cid:101) v || (cid:101) f l , s α − J l − S ξ (cid:34) N s ∑ β (cid:101) C l , s αβ + (cid:101) C l , s α e (cid:35) (5.59)is the nonlinear-residual of the ion VFP equation for the α species. Here, (cid:101) C αβ and (cid:101) C α e are the Picardlinearized (to the previous outer nonlinear iteration) α − β ion, and ion-electron collision operators, respec-tively, and we use a nonlinear history length of m s = P t + P ξ and P t + P (cid:101) v , where we employ one V cycle of classical geometricmultigrid, smoothed with 3 passes of damped Jacobi and employ agglomeration for restriction and secondorder prolongation. The nonlinear iteration is continued until (cid:13)(cid:13)(cid:13) R l , s α (cid:13)(cid:13)(cid:13) ≤ − (cid:13)(cid:13)(cid:13) R l , α (cid:13)(cid:13)(cid:13) . We note that the nonlinearity in the domain size –determined by the implosion boundary conditionsin Sec. 5.4– is not stiff and therefore the configuration-space grid locations and other derived quantitiesare Picard linearized with respect to the electron temperature and the ion distribution functions inside theevaluation of the outer fixed point map. All advection operators in the preconditioner, both in the fluidelectron and ion VFP systems, are discretized using a first-order upwind discretization. The evaluationof the outer fixed-point map in Eq. (5.48) is summarized in Alg. 2. The outer fixed point iteration iscontinued until (cid:13)(cid:13) R l α (cid:13)(cid:13) ≤ ε r (cid:13)(cid:13) R l = α (cid:13)(cid:13) , where ε r is the relative nonlinear convergence tolerance ( ε r = − unless otherwise specified).The full algorithm containing the NS-MMPDE for the configuration-space-grid adaptivity, velocity-space-grid adaptivity (for completeness, briefly summarized in Appendix F), time-step adaptivity, and non-linear solvers is described in Algorithm 3. Algorithm 3
Integrated algorithm of adaptive time-stepping, NS-MMPDE, velocity-space grid adaptivity,and the AA solver for the hybrid ion Vlasov-Fokker-Planck and fluid electron equation.1. Adapt time step-size according to Eq. (5.38).2. Adapt velocity-space grid according to Alg. 4.3. Adapt configuration space grid according to Alg. 1.4. Solve coupled ion-VFP and fluid electron system according to Eq. (5.48) and Alg. 2.17 . Numerical Results
We test the proposed integrated algorithm on a set of benchmark problems that demonstrate its advertisedcapabilities. In this study, the VFP and fluid electron temperature equations are non-dimensionalized interms of the proton mass, m p = . × − [ kg ], proton charge, q p = . × − [ C ], normalizationnumber density, n ∗ = . × [particles / m ], normalization temperature, T ∗ = .
654 [ eV ] = . × − [ J ]; and the resulting proton-proton collision time, τ ∗ = √ m p ε ( π T ∗ ) / n ∗ q p = . × − [ s ], speed, u ∗ = (cid:112) k B T ∗ / m p = . × [ m / s ], and length, L ∗ = u ∗ τ ∗ = . × − [ m ]. The v ∗ normalizeddistribution function is initialized based on the Maxwellian distribution function, (cid:101) f M = n M π / (cid:18) v ∗ v th , M (cid:19) exp (cid:34) − ( v − u M ) v th , M (cid:35) , (6.1)where the numerical Maxwellian moments, n M , u M , and v th , M = (cid:112) T M / m are computed to ensure thatthe numerical integrals of the Maxwellian distribution function return the prescribed density, drift velocity,and temperature [45]. Unless otherwise mentioned, we use a density- and temperature-dependent Coulomblogarithm, defined in Ref. [69]. Also, unless otherwise stated, an initial Poisson grid generation and op-timization (Appendix G) is employed to avoid initial numerical pollution of simulation due to unresolvedgradients. We consider an initially perturbed plasma in a spherical cavity, confined in a reflective system (i.e.,symmetry boundary conditions from Sec. 5.4.1 at both r | ξ = = r | ξ = = R ) to test the discretizationconvergence rate and conservation properties of the proposed algorithm. We consider a domain with anormalized radius of R = ξ × (cid:101) v || × (cid:101) v ⊥ ∈ [ , ] × [ − , ] × [ , ] . A deuterium-tritium plasma is used with masses, m D = m T =
3; charges, q D = q T =
1; number densities, n D = n T = + . ( k r r ) ; drift velocities, u D = u T =
0; and temperatures, T D = T T = T e = + . ( k r r ) , where k r = π R ; and the prescribed analytical grid is given as: r ( ξ i , t ) = r , i + . ∆ r sin (cid:0) ω grid t (cid:1) cos ( k r r , i ) . (6.2)Here, ω grid = π / ∆ r = R / N ξ is the average grid size in the physicalspace, and r = ξ R is the initial grid. We also use a Coulomb log of 10 for all interactions.For this section, the time-step ramping discussed in Sec. 5.5 is not used to allow for a time-convergencestudy, and the first time-step size is taken to be ∆ t max . We demonstrate the discrete conservation theoremfor mass and energy with a grid of N ξ × N v || × N v ⊥ = × ×
32, a time-step size of ∆ t max = − , andvarying nonlinear convergence tolerance of ε rel = − , 10 − , and 10 − ; refer to Fig. 6.1. As can be seen,discrete conservation error in mass, ∆ M ( p ) M ( ) , and energy, ∆ U ( p ) U ( ) , reduces commensurately with the nonlinearconvergence tolerance. Here, M ( p ) = N s ∑ α m α N ξ ∑ i = ∆ ξ J ( ) S ξ , i (cid:68) , (cid:101) f ( p ) α , i (cid:69) δ v , (6.3) ∆ M ( p ) = (cid:12)(cid:12)(cid:12) M ( p ) − M ( ) (cid:12)(cid:12)(cid:12) , (6.4) U ( p ) = N s ∑ α m α N ξ ∑ i = ∆ ξ J ( p ) S ξ , i (cid:28) v , (cid:101) f ( p ) α , i (cid:29) δ v + N ξ ∑ i = ∆ ξ J ( p ) S ξ , i n ( p ) e , i T ( p ) e , i , (6.5)and ∆ U ( p ) = (cid:12)(cid:12)(cid:12) U ( p ) − U ( ) (cid:12)(cid:12)(cid:12) . (6.6)18 igure 6.1: Doubly reflective boundary problem with prescribed grid motion: Discrete conservation of mass and energy for varyingnonlinear convergence tolerance. Next, we demonstrate the discretization order of convergence of the proposed algorithm. The temporal,configuration-space, and velocity-space convergence study is performed by measuring the L -norm of therelative difference of the temperature, E ∆ tT = (cid:118)(cid:117)(cid:117)(cid:116) N ξ ∑ i = ∆ ξ (cid:34) N sp ∑ α = (cid:16) T ∆ t α , i − T ∆ t ref α , i (cid:17) + (cid:16) T ∆ te , i − T ∆ t ref e , i (cid:17) (cid:35) , (6.7) E ∆ ξ T = (cid:118)(cid:117)(cid:117)(cid:116) N ξ ∑ i = ∆ ξ (cid:34) N sp ∑ α = (cid:16) T ∆ ξα , i − T ∆ ξ ref α , i (cid:17) + (cid:16) T ∆ ξ e , i − T ∆ ξ ref e , i (cid:17) (cid:35) , (6.8) E ∆ (cid:101) vT = (cid:118)(cid:117)(cid:117)(cid:116) N ξ ∑ i = ∆ ξ (cid:34) N sp ∑ α = (cid:16) T ∆ (cid:101) v α , i − T ∆ (cid:101) v ref α , i (cid:17) + (cid:16) T ∆ (cid:101) ve , i − T ∆ (cid:101) v ref e , i (cid:17) (cid:35) . (6.9)Here, the superscript ∆ t denotes the solution obtained with a larger time-step size while ∆ t re f denotes thereference solution, ∆ ξ denotes the solution obtained with a coarse logical space grid size and ∆ ξ re f denotesthe reference solution, and ∆ (cid:101) v denotes the solution obtained with a coarse velocity-space grid and ∆ (cid:101) v re f denotes the reference grid.The temporal discretization relies on a BDF2 scheme and is second order accurate (refer to Fig. 6.2-left).We compute the L -norm of the relative difference of the temperatures. The reference solution uses a gridof N ξ = N v || = N v ⊥ =
16, and a time-step size of ∆ t re f = − , and is time integrated to t max = . . The convergence study is performed by fixing the grids while varying the time-step size.The configuration-space discretization relies on a SMART discretization for the flux interpolation, whichasymptotically is third order accurate. However, we see in Fig. 6.2-center that we recover only second-orderaccuracy. This is due to the second-order accurate operations employed elsewhere (e.g., linear interpolationemployed for the configuration-space cell-face quantities [e.g., J r , J ξ , v ∗ , u ∗|| ] as well as cell-face evaluationof the gradients [second order centered differencing] of velocity-space metrics [e.g., ∂ ξ v ∗ and ∂ ξ u ∗|| ]). Wenote that the solution in the physical space is compared by interpolating the reference solution to the coarsegrid. The reference solution uses a grid of N ξ = N v || = N v ⊥ =
16, and a time-step size of ∆ t = . t max = . . The convergence study is performed by fixing the time-step size andvelocity-space resolution while varying the number of configuration-space grid points.Similarly to the configuration space, the velocity-space discretization also relies on the SMART schemefor the flux interpolation. However, due to the linear interpolation used elsewhere (e.g., interpolation of19dvection coefficients to the cell face), the overall discretization is second order as is seen in Fig. 6.2-right. The reference solution uses a grid of N ξ = N v || = N v ⊥ = ∆ t = . t max = . . The convergence study is performed by fixing the time-step size andconfiguration-space grid while uniformly refining the transformed velocity-space grid in both (cid:101) v || and (cid:101) v ⊥ . Figure 6.2: Doubly reflective boundary problem with prescribed grid motion: Temporal (left), configuration-space (center), andvelocity-space (right) convergence study. As can be seen, second order convergence rate is observed in all discretization parameters.
We consider the Guderley problem [70, 71, 72, 73] of a converging/diverging shock to test the proposedtime-dependent boundary conditions in spherical geometry and to demonstrate the ability of the algorithmto capture the hydrodynamic asymptotic limit for short mean-free-paths. Semi-analytic solutions to theEuler equations exist for the Guderley problem (see Appendix H for details), which we will appeal to forverification. Suppose that a strong shock is spherically converging into a uniform, Deuterium-Tritium, (cid:104) DT (cid:105) ,plasma with constant mass density of ρ = m (cid:104) DT (cid:105) n (cid:104) DT (cid:105) = .
173 [ g / cc ] and temperature T =
10 [ eV ] (i.e.,conditions relevant to double-shell and revolver ICF implosions[74, 75, 76, 77, 78, 79, 80, 81, 82, 83]). Wemodel the (cid:104) DT (cid:105) plasma as an equimolar, averaged single ion, fully ionized, and initialize the distributionfunction with density, drift velocity, and pressure/temperature profiles obtained from a Guderley solution.The Guderley profiles are obtained 20 [ ps ] after the shock was initialized from a radius of R =
200 [ µ m ],and the Mach number there is M = .
23. Similarly, the hydrodynamic boundary conditions for the sphericaldomain are sourced from the same Guderley solution. The logical grid sizes for these calculations are chosenas N ξ =
192 and 384, N v ⊥ =
64 and N v || =
128 with domain limits ξ × (cid:101) v || × (cid:101) v ⊥ ∈ [ , ] × [ − , ] × [ , ] . Theprofiles for the initial and the boundary conditions for the moments are shown in Fig. 6.3.We start our verification study by testing on a static uniform mesh. For the verification metric, we choosethe shock speed, u s = (cid:12)(cid:12)(cid:12)(cid:12) dr s dt (cid:12)(cid:12)(cid:12)(cid:12) , (6.10)where r s is the shock position. The iFP shock trajectory is defined as either: 1) the location of the maximumion temperature, or 2) the location of the maximum ion viscous (energy) dissipation rate. The first definitionis the most reliable (i.e., it has less dependence on small-scale fluctuations in the moment quantities), butthe maximum temperature only unambiguously identifies the shock location in the converging phase. Forour simulations, we choose definition 1 for the converging shock phase, while 2 for the diverging phase. Wedefine the viscous dissipation rate as: (cid:18) dT i dt (cid:19) visc = − (cid:18) τ |||| n i (cid:19) ∂ r (cid:0) r u || , i (cid:1) r , (6.11)20
200 400 600t [ps]2.02.53.03.54.04.5 n [ / m ] u || [ μ m / n s ] μ [ e V ] Figure 6.3: Guderley problem: The initial condition (top) and the time-dependent Maxwellian boundary conditions (bottom) for thenumber density (left), radial drift velocity (center), and temperature (right), obtained from semi-analytically solving the Guderleyproblem in Ref. [71]. For the initial condition, the black markers denotes the grid density (i.e., coarse grid regions are lighter incolor than finer regions). The boundary temperature is assumed to be equilibrated across ions and electrons at all time. where, T i is the ion temperature, u || , i is the radial component of the ion drift velocity, n i is the total ionnumber density, and τ |||| is the radial component of the ion viscosity tensor [Eq. (5.35)]. Fig. 6.4-leftshows a comparison of the shock speed versus time for the two different spatial-grid resolutions. The twocurves show strong deviations from each other near the shock collapse time of approximately 500 [ps]. Thenoticeable fluctuations in the curves suggest that the simulations are poorly resolving the shock structure,where the spatial gradients can change dramatically in the course of the shock’s collapse. We see a dramatic S h o c k s p ee d [ c m / µ s ] static mesh (N ξ = 192)static mesh (N ξ = 384) 0 100 200 300 400 500Time [ps]20304050607080 S h o c k s p ee d [ c m / µ s ] NS-MMPDE (N ξ = 192)NS-MMPDE (N ξ = 384) Figure 6.4: Guderley problem: The shock speed versus time with a static and uniform grid (left) and a moving mesh with anoptimized initial mesh generation (right). improvement in the quality of the solution when we employ an initial grid optimization in conjunction withthe NS-MMPDE strategy with δ min = × − and λ ω = − . Indeed, as shown in Fig. 6.4-right, thedifferences between N ξ =
192 and N ξ =
384 disappear, and the fluctuations seen in Fig. 6.4-left are nearlyeliminated for both cases. Moreover, comparing the static and the moving mesh cases for N ξ =
384 reveals21he origin of the deviations in the shock trajectory near the collapse time: As seen in Fig. 6.5, the shocktrajectory with the uniform mesh veers away from the correct trajectory immediately after the simulationbegins. The initial shock structure is not well resolved in the static mesh simulation, and this seems to be the S h o c k s p ee d [ c m / µ s ] static mesh (N ξ = 384)NS-MMPDE (N ξ = 384) Figure 6.5: Guderley problem: The shock speed versus time comparing the two simulations for N ξ = cause of the strong deviation. The drifting Maxwellian for the temperature, density, and drift velocity arenot the true distribution functions for the shock, and this leads to the production of transients in the solutionwhich, if not properly resolved, can cause the solution to strongly deviate from initial condition.We expect that such a high (cid:104) DT (cid:105) density (0 .
173 [ g / cc ]) justifies the use of boundary conditions sourcedfrom a dissipation-less fluid theory (i.e., Euler equations), and that a meaningful comparison to this theoryis actually possible. To this end, we compare our N ξ =
192 NS-MMPDE results to the correspondingGuderley solution. Fig. 6.6-left shows the shock trajectory as a function of time from the simulation and thesemi-analytic Guderley prediction. We see that the shock collapse times agree very well. Remarkably, theGuderley shock collapse time and the simulation agree to within 0 . S h o c k p o s i t i o n [ µ m ] iFPGuderley 500 550 600 650 700Time [ps]0102030405060 S h o c k p o s i t i o n [ µ m ] iFPGuderley Figure 6.6: Guderley problem: The shock trajectory vs. time. The blue (“solid”) curve is obtained from our algorithm for N ξ =
192 with the proposed NS-MMPDE and initial grid optimization algorithms, and the red (“dashed”) curve is the semi-analyticalGuderley prediction over the entire simulation (left) and the isolated time at the vicinity of shock collapse (right). in the neighborhood of the collapse, as depicted in the “zoomed in” trajectory plot shown in Fig. 6.6-right,we see that the Guderley and simulated shock trajectories differ in the early post-collapse phase. Insight intothis discrepancy may be gained by directly comparing the simulation and Guderley hydrodynamic profiles.Fig. 6.7-top shows pre-collapse profiles for density (left), drift velocity (center), and total hydrodynamicpressure, P = n i T i + n e T e (right). We see that there is some separation between the simulation and Guderleydensities at an early time (t = 100 [ ps ] after the simulation is initialized). Moreover, the transients created22n the simulation initialization have yet to be fully eliminated. By t = 400 [ ps ], however, the simulation andGuderley profiles largely agree, although we see some remaining differences near the shock collapse time(t = 552.5 [ ps ]), which persist into the diverging phase – as depicted in Fig. 6.7-bottom. Nevertheless, thequalitative features are nearly identical between simulation and Guderley. The strong deviations in the hydro N u m b e r d e n s i t y [ c m − ] t = 100.0 pst = 400.0 pst = 552.5 ps D r i f v e l o c i y [ c m / µ s ] iFPGuderley = 100.0 ps = 400.0 ps = 552.5 ps P r e ss u r e [ e V / cc ] t = 100.0 pst = 400.0 pst = 552.5 ps N u m b e r d e n s i t y [ c m − ] t = 600.0 ps t = 640.0 ps t = 700.0 ps D r i f v e l o c i y [ c m / µ s ] iFPGuderley = 600.0 ps = 640.0 ps = 700.0 ps P r e ss u r e [ e V / cc ] t = 600.0 ps t = 640.0 ps t = 700.0 ps Figure 6.7: Guderley problem: Pre-collapse hydro profiles from simulation (“solid”) and Guderley (“dashed”). profiles seen in Fig. 6.7-bottom hint at the origin of the trajectory discrepancy in the post-collapse phase –as shown in Fig. 6.6-right. The Guderley solution predicts a zero density and (infinite temperature) at theorigin at the shock collapse time. This behavior is unphysical, and it results from the infinite collisionalityassumption in the Euler equations [72]. The VFP equation models the full dissipation physics, and thuscaptures the correct behavior near shock collapse (where ion heat conduction, viscosity, and kinetic effectsbecome important).In Fig. 6.8 we show, as a function of time, the max ( ∆ r ) / min ( ∆ r ) and (cid:104) ∆ r (cid:105) / min ( ∆ r ) (where (cid:104) ∆ r (cid:105) = ∑ N ξ i = ∆ r i / N ξ ) to demonstrate the computational complexity reduction (relative to a uniform grid) affordedby the proposed NS-MMPDE scheme. As can be seen, both ratios exhibit large values and a reduction time [ps]0 200 400 600 m a x ( ∆ r) / m i n ( ∆ r) < ∆ r > / m i n ( ∆ r) Figure 6.8: Guderley problem: Ratios of maximum to minimum (left) and average to minimum (right) grid size as a function oftime. of approximately 20 in the computational complexity is achieved by the mesh motion for the Guderleyproblem. In Fig. 6.9, a Lagrangian radius-time (RT) diagram of 1 / ∆ r is shown to demonstrate how the grid23racks the shock front as a function of time. As can be seen, the grid tightly concentrates near the sharp Figure 6.9: Guderley problem: RT diagram of 1 / ∆ r for N ξ = converging shock front, and diffuses away after shock collapse due to kinetic effects smearing the gradients.Finally, to demonstrate the computational complexity reduction afforded by the velocity-space adaptivityscheme, we present the parallel velocity-space marginal distribution function, f || ,< DT > = π (cid:90) ∞ dv ⊥ v ⊥ f (cid:104) DT (cid:105) , (6.12)in the r − v || plane in Fig. 6.10. For the calculation, max ( v th ) / min ( v th ) ≈ , thus a total computational Figure 6.10: Guderley problem: Parallel velocity-space marginal distribution functions at t ≈ , , ,
700 [ps]. complexity reduction of max (cid:16) (cid:104) ∆ r (cid:105) min ( ∆ r ) (cid:17) (cid:16) max ( v th ) min ( v th ) (cid:17) ≈ Consider a spherical piston surrounding a gas at constant pressure and density. At t =
0, the piston beginsmoving at a constant velocity inward. A shock is created in the gas, and converges to the center. We simulatethis setup directly, using the moving elastic wall boundary condition discussed in Sec. 5.4.1. Once more,we consider a fully ionized (cid:104) DT (cid:105) gas with an initial mass density of ρ (cid:104) DT (cid:105) = m (cid:104) DT (cid:105) n (cid:104) DT (cid:105) = .
173 [ g / cc ] andtemperature of T (cid:104) DT (cid:105) =
10 [ eV ]. The piston velocity was chosen to be u W = −
19 [ cm / µ s ] inwards, and theinitial radius was R =
326 [ µ m ]. To avoid any violent, abrupt changes in the ion distribution function atthe wall, we accelerate the wall to the desired piston velocity (-19 [ cm / µ s]) at a constant rate. Consideringtwo different acceleration times (2.5 and 5 [ ps ]), we found negligible differences in the final solution. Ourconvergence study begins with a static mesh ( N ξ = N v ⊥ = N v || = ps ] later than the Van Dykeprediction [84] of 908 [ ps ]. Nonetheless, this is only about a 3% difference. The collapse time furtherimproves when we use the NS-MMPDE algorithm, going from 938 →
928 [ ps ]. This figure is largelyinsensitive to the grid resolution, and the non-linear convergence tolerance; changing only to 927 [ ps ] for N ξ =
384 (keeping all other parameters constant), and 930 [ ps ] when N ξ =
192 and the relative tolerance ischanged from ε rel = − to ε rel = − . The largest gain was attained by enhancing the grid resolution nearthe piston using an initial grid optimization. We once more used N ξ = ps ], putting the agreement with the Van-Dyke solution within ≈ . S h o c k p o s i t i o n [ µ m ] iFPVan Dyke 200 400 600 800 1000Time [ps]050100150200250300350 S h o c k p o s i t i o n [ m i c r o n s ] iFPGuderley Figure 6.11: Van Dyke problem: The shock trajectory versus time. The blue (“solid”) curve is obtained from iFP using the elasticmoving wall boundary condition, and the red (“dashed”) curve is the Van Dyke (left) and Guderley (right) predictions.
Finally, the Van Dyke solution asymptotically approaches the Guderley solution as the shock collapses.Fig. 6.11-right shows decent agreement between iFP and the equivalent Guderley-shock trajectories. Asbefore, we can also compare the numerical and Guderley hydro profiles. Figs. 6.12-top shows pre-collapseprofiles for density (left), drift velocity (center), and total hydrodynamic pressure (right). Similarly to thetime-dependent boundary condition simulations in Section 6.2, the simulation and Guderley profiles agreeexceedingly well. The shock collapse times agree within 0 .
50 100 150 200 250 300Radius [µ]0.51.01.52.02.53.03.54.04.5 N u m b e r d e n s i t y [ c m − ] t = 100.0 pst = 550.0 pst = 903.0 ps D r i f t v e l o c i t y [ c m / µ ] iFPGuderley t = 100.0 p t = 550.0 p t = 903.0 p P r e ss u r e [ e V / cc ] t = 100.0 pst = 550.0 pst = 903.0 ps N u m b e r d e n s i t y [ c m − ] t = 920.0 ps t = 1000.0 ps t = 1100.0 ps D r i f v e l o c i y [ c m / µ s ] iFPGuderley = 920.0 ps = 1000.0 ps = 1100.0 ps P r e ss u r e [ e V / cc ] t = 920.0 ps t = 1000.0 ps t = 1100.0 ps Figure 6.12: Van Dyke problem: Pre-collapse (top) and post-collapse (bottom) hydro profiles for density (left), drift velocity(center), and total hydrodynamic pressure (right) – given the moving wall boundary condition – from simulation (“solid”) andGuderley (“dashed”). the reflected shock. Further, since the Guderley problem is strictly only applicable for an infinite domain,the solutions continue to deviate as the wall comes closer to the diverging shock front.Similarly to the Guderley problem, we demonstrate the grid savings in the configuration space byshowing the maximum- and average-to-minimum grid size ratios in Fig. 6.13. As can be seen, a com-
Figure 6.13: Van Dyke problem: Ratios of maximum to minimum (left) and average to minimum (right) grid size as a function oftime. plexity reduction in configuration space of roughly 50 is achieved for the Van Dyke problem. Similarlyto the Guderley problem, max ( v th ) / min ( v th ) ≈ , thus a total computational complexity reduction ofmax (cid:16) (cid:104) ∆ r (cid:105) min ( ∆ r ) (cid:17) (cid:16) max ( v th ) min ( v th ) (cid:17) ≈ / ∆ r to demonstrate how the grid tracks the shock front as a function of time.As can be seen, the grid tracks several shock reflections from the center and the wall, until gradients broadendue to finite transport effects and the grid relaxes. 26 igure 6.14: Van Dyke problem: RT diagram of 1 / ∆ r for N ξ = We reproduce the results reported in Ref. [46] in which a fuel-only simulation was performed for acompression yield target to study the so-called
Rygg effect [7, 8, 46]. We consider this problem for threereasons: 1) to test the implementation of the pressure-driven Lagrangian boundary conditions on a realisticimplosion experiment; 2) to demonstrate the performance of our nonlinear solver strategy applied to arealistic implosion scenario; and 3) to demonstrate the algorithmic advantage of our proposed NS-MMPDEscheme over the classic MMPDE scheme.The Rygg effect refers to the anomalous (non-hydrodynamic) scaling of experimental yield observed byRef. [7]. In the reference, experiments were performed to test the validity of the hydrodynamic assumptionsin a plastic (CH) capsule with D- He filled capsules at the Omega laser ICF facility. In the experiments, thefuel species concentrations were varied from shot to shot such that the initial mass density, ρ = ∑ N s α m α n α ,and total hydrodynamic pressure, P = ∑ N s α ( + Z α ) n α T α , was kept fixed, thus ensuring hydro-equivalence .Here, m D = m p , m He = m p and Z α is the ionization state of ion-species α . Accordingly, if hydrodynamicpredictions are valid, the nuclear yield of D-D reaction, Y DD ∝ n D , should follow a simple f D scaling law,where f D is the deuterium number fraction, defined as: f D = n D n D + n He . (6.13)The hydrodynamic variables were obtained from the LILAC Lagrangian rad-hydro simulation for Omegafacility experimental shot 22862 at the last fuel region. The simulation was initialized at t = . ps ] atthe moment in which the shock breaks out of the fuel-pusher interface; refer to Fig. 6.15. We note that atthis point the fuel species’, Deuterium and Helium-3, are assumed to be fully ionized ( Z D = Z He = N ξ = N v || = N v ⊥ =
64, and a transformed velocity-domain size of (cid:101) v || ∈ [ − , ] and (cid:101) v ⊥ ∈ [ , ] with grid parameters, λ v ∗ = λ u ∗ = − , δ min = . λ ω = − , τ r = . ps ], and ∆ t max = . ps ]. In Fig. 6.16, the simulation results are shown for the case of f D = . ps ]). The separations are supported due to the differencein v th = (cid:112) T / m amongst the ion species and the long D- He mean-free-paths –which scales as λ D − He ∝ T / (cid:0) Z D Z He (cid:1) – supported within the shock; refer to Fig. 6.17. As can be seen, a distinct beam feature formsin the Deuterium species, allowing a significant population to run ahead of the Helium-3. These structuresare inherently kinetic and far from a linear perturbation in pitch angle, necessitating a general description27 igure 6.15: Fuel-only compression yield simulation: The initial condition (top) and the time-dependent Maxwellian boundaryconditions (bottom) for the number density, n , radial drift velocity, u , and temperature, T , obtained from the LILAC rad-hydrocalculation for the last Lagrangian fuel zone at the time in which the first shock breaks out of the fuel-pusher interface (~500.5[ ps ]). For the initial condition for number densities, the blue (red) curve denotes Deuterium (Helium-3); for the temperature, blue(green) denotes ions (electrons); and the black markers denotes the grid density (i.e., coarse grid regions are lighter in color thanfiner regions). We remind the readers that the boundary position and velocity are evolved from the mass conservation constraintand the Lagrangian equation of motion [Eqs. (5.29) and 5.31], respectively, and are not explicitly provided to drive the simulation.Further the boundary temperatures and drift velocities are assumed to be equilibrated across species at all times. for the distribution function. In fact, one observes that the plasma near the shock front is appreciably kineticthroughout the entire course of the implosion process as can be seen in the Knudsen number RT diagramin Fig. 6.18. As seen, although the Knudsen number does reduce to Kn ≤ − at peak-compression time,hence equilibrating the distribution function roughly to a Maxwellian, the early separation in the speciesdensities persists due to the early kinetic imprinting in the solution. As a consequence, Deuterium speciesare depleted from the hot core and transported to the periphery of the domain, leading to a degradation inthe D-D yield.In Fig. 6.19 the normalized (to the f D = DD -neutron nuclear yield, (cid:101) Y f D DD = Y f D DD Y f D = DD (cid:32) n f D = D n f D D (cid:33) , (6.14)is shown, where Y DD = π (cid:90) t max t dt (cid:90) r dr (cid:90) d v f D ( r , v ) (cid:90) d v (cid:48) f D (cid:0) r , v (cid:48) (cid:1) (cid:12)(cid:12) v − v (cid:48) (cid:12)(cid:12) g (cid:0)(cid:12)(cid:12) v − v (cid:48) (cid:12)(cid:12)(cid:1) , (6.15) g ( | v − v (cid:48) | ) is the DD − n cross section [85], and n f D D is the initial (un-shocked) Deuterium number densityfor a given f D . As can be seen, with decreasing f D (increasing f He ), the scaled relative yield decreases.In a single-fluid hydrodynamic formalism, the total mass and pressure would be independent of speciesconcentrations, leading to no variation in (cid:101) Y f DDD . In kinetic theory however, differential-ion motion, viscousheating, and other kinetic physics break hydro-equivalence and lead to dependence of (cid:101) Y f D DD with Deuterium28 igure 6.16: Fuel-only compression yield simulation: The density (top), drift velocity (center), and temperature (bottom) for the f D = . ps ]), upon shock reflection (~1650 [ ps ]) and at near peak compression andstagnation (~2050 [ ps ]). The black markers denotes the grid density. concentration. For further insights and discussion into the yield variation physics, we refer the reader to Ref.[46].In Fig. 6.20, we demonstrate the performance of our solver strategy by showing, for the f D = . ≤ t ≤ ps ] time interval, where nontrivial structures are evolvingin the simulation (e.g., second shock launched into the system, interaction and merging of the two shocks,and the subsequent fast evolution of physical structures), leading to occasional spikes in the number ofnonlinear iterations. However, as can be seen from a single pass of binomial smoothing, the times in whichthe nonlinear iteration exceeds 10 is rare.Finally, once again for the f D = . ∆ t afforded by the NS-MMPDE scheme(relative to the standard MMPDE) by comparing the instantaneous ∆ t as a function of time; refer to Fig.6.21. As can be seen, the NS-MMPDE scheme permits a generally larger ∆ t with less erratic behaviorthan the standard MMPDE approach. The erratic behavior of MMPDE is attributed to the periodic jitter in the grid caused by the traveling shock front, which occasionally causes the grid near the shock frontto evolve rapidly [and consequently a smaller ∆ t is taken according to Eq. (5.38)]. In fact, the erraticgrid evolution eventually leads to very small ∆ t , rendering the simulation impractical for engineering-scalesimulations. The only stable integration for the classic MMPDE scheme is obtained for a τ r =
10 [ ps ] case,where the grid evolution time-scale is set to be longer than the dynamical time-scales in the system. Even29 igure 6.17: Fuel-only compression yield simulation: The temperature (top) of the various plasma species, and the correspondingion marginal parallel-velocity distribution function (bottom) for Deuterium (left) and Helium-3 (right). then, occasional erratic behaviors are seen at around 1725 [ ps ]. In contrast, the proposed NS-MMPDEscheme detects regions where the grid is evolving rapidly in the predictor stage and attempts to minimizeit through the optimization procedure, leading to a more robust behavior in the time-integration throughoutthe simulation. Thus, the proposed NS-MMPDE scheme greatly relaxes the process of having to empirically tune the τ r .
7. Conclusions and Future Work
We have demonstrated a conservative (mass, energy) phase-space moving-grid strategy in a sphericallyimploding system with a variety of boundary conditions - relevant to modeling ICF capsule implosions. Thephase-space grid-adaptivity algorithm leverages the work from Ref. [50] by adapting the velocity-spacegrid based on the instantaneous thermal speed and bulk flow of the individual plasma species (i.e., multiplevelocity-space grid). In the configuration space, the grid is evolved based on the MMPDE formalism. Todeal with strong shocks and complex features encountered in implosion problems, we have developed anonlinear stabilization strategy for the classic MMPDE algorithm (NS-MMPDE). The approach nonlinearlystabilizes MMPDE against numerical instabilities that can arise based on a fixed grid relaxation time-scale, τ r . The strategy is based on splitting the grid evolution process in two stages: 1) predict the grid based onMMPDE; and 2) pass the predicted grid through a nonlinear constrained optimization stage which optimizesthe grid against the predicted value and the volumetric rate change. Our fully implicit solver is based on anested Anderson accelerated fixed point iteration strategy that efficiently deals with the nonlinear couplingbetween the kinetic ion and fluid electrons. We have demonstrated the capability of the NS-MMPDE strat-egy (coupled with the previous velocity-space adaptivity capability and the adaptive time-stepping strategy)to robustly simulate ICF implosion problems including the Guderley/Van-Dyke problem –as far as we areaware, the first to do so in the literature with a Vlasov-Fokker-Planck approach. We have also demonstrated30 igure 6.18: Fuel-only compression yield simulation: The Deuterium Knudsen number RT diagram is shown. Here, the Knudsennumber is defined as Kn = v th , D τ e f f ∂ r ln n ion , where n ion = n D + n he , v th , D = (cid:112) T D / m D , τ e f f = (cid:16) τ − DD + τ − D He + τ − De (cid:17) − , and τ D β ∝ T / β / (cid:16) Z D Z β n β m / β (cid:17) . As can be seen, the Knudsen number is relatively large near the shock trajectory while at peakcompression, it reduces to below 10 − . with a realistic Omega implosion experiment the integrated capability of the proposed algorithm in effec-tively dealing with multi-length and -velocity scale problems. Future work will include the extension of theproposed capability to: 1) a fully kinetic system (e.g., kinetic electron and Ampère equations) by advancingthe work performed by Ref. [51] to spherical geometry; 2) to multiple spatial dimensions (e.g., 2D3V) ina general curvilinear coordinate system; 3) couple with radiation transport physics; and 4) couple with ahigh-order low-order (HOLO) algorithm [86, 87, 88, 89, 90] to step over the stiff collision time-scales whendense and cold ablators are included in the simulation. We note that all works have already been performedand will be documented in follow-on publications. Acknowledgments
We greatly appreciate Dr. Brian J. Albright from LANL for the many helpful discussions and particu-larly suggestions regarding the physics applications for the proposed capability. We appreciate Dr. OlivierLarroche from CEA for providing the LILAC rad-hydro initial and boundary conditions that were used forthe Rygg problem and the many helpful discussions regarding clarifications on the FPION and FUSE ca-pabilities. We also thank Dr. Hong Sio from the Lawrence Livermore National Laboratory for discussionsin simulating exploding pusher capsules, which ultimately lead to the need for improving our mesh mo-tion strategy (and the eventual development of the NS-MMPDE scheme). This work was sponsored by theMetropolis Postdoctoral Fellowship for W.T.T. between the years 2015-2017, the Thermonuclear Burn Ini-tiative of the Advanced Simulation and Computing Program between years 2018-2020, the revolver LDRDproject for B.D.K for year 2019, and the Institutional Computing program at the Los Alamos National Lab-oratory. This work was performed under the auspices of the National Nuclear Security Administration ofthe U.S. Department of Energy at Los Alamos National Laboratory, managed by Triad National Security,LLC under contract 89233218CNA000001.
References [1] Y. B. Zel’dovich and Y. P. Raizer,
Physics of shock waves and high-temperature hydrodynamic phe-nomena . Dover Publications, INC., 2 ed., 1967.31 igure 6.19: Fuel-only compression yield simulation: (cid:101) Y f D DD for varying f D . The blue markers denote values obtained in this study,while the reference values obtained from Ref. [46] are shown in red. We note that the observed small disagreement with thereference result is expected due to the significant algorithmic improvements we have made to iFP (e.g., improved treatment for thephase-space grid adaptivity, dynamic time-step adaptivity, etc.)Figure 6.20: Fuel-only compression yield simulation: Instantaneous number of outer fixed point iteration, N , versus time (blue) andthe same quantity with a single pass of binomial smoothing (red), N ∗ , ( p ) = N ( p + ) + N ( p ) + N ( p − ) , where the superscript (*) denotesthe smoothed value, and (p) denotes the time-index. [2] R. P. Drake, High-energy-density physics . Springer, 2 ed., 2018.[3] L. Welser-Sherrill, D. A. Haynes, R. C. Mancini, J. H. Cooley, R. Tommasini, I. E. Golovkin, M. E.Sherrill, and S. W. Haan, “Inference of icf implosion and core mix using experimental data and theo-retical mix modeling,”
High Energy Dens. Phys. , vol. 5, pp. 249–257, 2009.[4] J. A. Baumgaertel, P. A. Bradley, S. C. Hsu, J. A. Cobble, P. Hakel, I. Tregillis, N. S. Krasheninnikova,T. J. Murphy, M. J. Schmitt, R. C. Shah, K. D. Obrey, S. Batha, H. Johns, T. Joshi, D. Mayes, R. C.Mancini, and T. Nagayama, “Observation of early shell-dopant mix in omega direct-drive implosionsand comparisons with radiation-hydrodynamic simulations,”
Phys. Plasmas , vol. 21, p. 052706, 2014.[5] H. G. Rinderknecht, H. Sio, C. K. Li, A. B. Zylstra, M. J. Rosenber, P. Amendt, J. Delettrez, C. Bellei,J. A. Frenje, M. G. Johnson, and et al., “First observations of nonhydrodynamic mix at the fuel-shellinterface in shock-driven inertial confinement implosions,”
Phys. Rev. Lett. , vol. 112, p. 135001, 2014.[6] T. Ma, P. K. Patel, N. Izumi, P. T. Springer, M. H. Key, L. J. Atherton, L. R. Benedetti, D. K. Bradley,32 nitial Δ𝑡 ramp-up phase Adiabatic compression phase. Domain getting very small, Δ𝑡 is reduced, then increased after entering expansion phaseShock collapseat r=0 center, sharp gradients evolving. Δ𝑡 is reduced Figure 6.21: Fuel-only compression yield simulation: Instantaneous ∆ t versus time for NS-MMPDE and MMPDE for various gridrelaxation time-scales. D. A. Callahan, P. M. Celliers, and et al., “Onset of hydrodynamic mix in high-velocity, highly com-pressed inertial confinement fusion implosions,”
Phys. Rev. Lett. , vol. 111, p. 085004, 2013.[7] J. R. Rygg, J. A. Frenje, C. K. Li, F. H. Seguin, and R. D. Petrasso, “Tests of the hydrodynamicequivalence of direct-drive implosions with different d and he mixtures,” Phys. Plasmas , vol. 13,p. 052702, 2006.[8] H. W. Hermann, J. R. Langenbrunner, J. M. Mack, J. H. Cooley, D. C. Wilson, S. C. Evans, T. J. Sedillo,G. A. Kyrala, S. E. Caldwell, and et al., “Anomalous yield reduction in direct-drive deuterium/tritiumimplosions due to he addition,” Phys. Plasmas , vol. 16, p. 056312, 2009.[9] A. Thomas, M. Tzoufras, A. Robinson, R. Kingham, C. Ridgers, M. Sherlock, and A. Bell, “A re-view of Vlasov-Fokker-Planck numerical modeling of inertial confinement fusion plasma,”
J. Comput.Phys. , vol. 231, pp. 1051–1079, 2012.[10] M. J. Rosenberg, H. G. Rinderknecht, N. M. Hoffman, P. A. Amendt, S. Atzeni, A. B. Zylstra, C. K.Li, F. H. Séguin, H. Sio, M. G. Johnson, J. A. Frenje, R. D. Petrasso, V. Y. Glebov, C. Stoeckl, W. Seka,F. J. Marshall, J. A. Delettrez, T. C. Sangster, R. Betti, V. N. Goncharov, D. D. Meyerhofer, S. Skupsky,C. Bellei, J. Pino, S. C. Wilks, G. Kagan, K. Molvig, and A. Nikroo, “Exploration of the transitionfrom the hydrodynamic like to the strongly kinetic regime in shock-driven implosions,”
Phys. Rev. ,vol. 112, p. 185001, 2014.[11] C. Bellei, P. A. Amendt, S. C. Wilks, M. G. Haines, D. T. Casey, C. K. Li, R. Petrasso, and D. R.Welch, “Species separation in inertial confinement fusion fuels,”
Phys. Plasmas , vol. 20, p. 012702,2013.[12] P. A. Amendt, S. C. Wilks, C. Bellei, C. K. Li, and R. D. Petrasso, “The potential role of electricfields and plasma barodiffusion on the inertial confinement fusion database,”
Phys. Plasmas , vol. 18,p. 056308, 2011.[13] L. Yin, B. J. Albright, W. Taitano, E. L. Vold, L. Chacón, and A. N. Simakov, “Plasma kinetic effectson interfacial mix,”
Phys. Plasmas , vol. 23, p. 112302, 2016.3314] B. D. Keenan, A. N. Simakov, W. T. Taitano, and L. Chacón, “Ion species stratification within strongshocks in two-ion plasmas,”
Phys. Plasmas , vol. 25, p. 032103, 2018.[15] B. D. Keenan, A. N. Simakov, L. Chacón, and W. T. Taitano, “Deciphering the kinetic structure ofmulti-ion plasma shocks,”
Phys. Rev. E , vol. 96, p. 053203, 2017.[16] G. Kagan and X. Tang, “Thermo-diffusion in inertially confined plasmas,”
Phys. Lett. A , vol. 378,pp. 1531–1535, 2014.[17] A. N. Simakov and K. Molvig, “Hydrodynamic description of an unmagnetized plasma with multipleion species. ii. two and three ion species plasmas,”
Phys. Plasmas , vol. 23, p. 032116, 2016.[18] A. N. Simakov, B. D. Keenan, W. T. Taitano, and L. Chacón, “Plasma ion stratification by weak planarshocks,”
Phys. Plasmas , vol. 24, p. 092702, 2017.[19] A. Le, T. J. T. Kwan, M. J. Schmitt, H. W. Herrmann, and S. H. Batha, “Simulation and assessment ofion kinetic effects in a direct-drive capsule implosion experiment,”
Phys. Plasmas , vol. 23, p. 102705,2016.[20] H. G. Rinderknecht, M. Rosenberg, C. Li, N. Hoffman, G. Kagan, A. Zylstra, H. Sio, J. A. Frenje,M. G. Johnson, F. Séguin, R. Petrasso, and et al., “Ion thermal decoupling and species separation inshock-driven implosions,”
Phys. Rev. Lett. , vol. 114, p. 025001, 2015.[21] H. Sio, J. Frenje, A. Le, S. Atzeni, T. Kwan, M. G. Johnson, G. Kagan, C. Stoeckl, C. Li, C. Parker,C. Forrest, and et al., “Observations of multiple nuclear reaction histories and fuel-ion species dy-namics in shock-driven inertial confinement fusion implosions,”
Phys. Rev. Lett. , vol. 122, p. 035001,2019.[22] H. Sio, J. Frenje, J. Katz, C. Stoeckl, D. Weiner, M. Bedzyk, V. Glebov, C. Sorce, M. G. Johnson,H. Rinderknecht, A. B. Zylstra, and et al., “A particle X-ray temporal diagnostic (PXTD) for studiesof kinetic, multi-ion effects, and ion-electron equilibration rates in inertial confinement fusion plasmasat OMEGA,”
Rev. Sci. Instrum. , vol. 87, p. 11D701, 2016.[23] J. D. Sadler, Y. Lu, B. Spiers, M. W. Mayr, A. Savin, R. Wang, R. Aboushelbaya, K. Glize, R. Bingham,H. Li, K. A. Flippo, and P. A. Norreys, “Kinetic simulation of fusion ignition with hot-spot ablatormix,”
Phys. Rev. E , vol. 100, p. 033206, 2019.[24] E. Vold, L. Yin, W. Taitano, K. Molvig, and B. J. Albright, “Diffusion-driven fluid dynamics in idealgases and plasmas,”
Phys. Plasmas , vol. 25, p. 062102, 2018.[25] E. Vold, R. Rauenzahn, and A. Simakov, “Multi-species plasma transport in 1D direct-drive ICF sim-ulations,”
Phys. Plasmas , vol. 26, p. 032706, 2019.[26] V. T. Tikhonchuk, “Physics of laser plasma interaction and particle transport in the context of inertialconfinement fusion,”
Nuclear Fusion , vol. 59, p. 032001, 2019.[27] M. Casanova, O. Larroche, and J. Matte, “Kinetic simulation of a collisional shock wave in a plasma,”
Phys. Rev. , vol. 67, no. 16, pp. 2143–2146, 1991.[28] O. Larroche, “Kinetic simulation of a plasma collision experiment,”
Phys. Fluids B , vol. 5, pp. 2816–2840, 1993.[29] F. Vidal, J. P. Matte, M. Casanova, and O. Larroche, “Ion kinetic simulations of the formation andpropagation of a planar collisional shock wave in a plasma,”
Phys. Plasmas , vol. 5, p. 3182, 1993.3430] F. Vidal, J. Matte, M. Casanova, and O. Larroche, “Spherical ion kinetic simulations of dt implosions,”
Phys. Rev. E , vol. 52, no. 4, p. 4568, 1995.[31] F. Vidal, J. P. Matte, M. Casanova, and O. Larroche, “Modeling and effects of nonlocal electron heatflow in planar shock waves,”
Phys. Plasmas , vol. 2, p. 1412, 1995.[32] O. Larroche, “Kinetic simulations of fuel ion transport in ICF target implosions,”
Eur. Phys. J. D ,vol. 27, pp. 131–146, 2003.[33] O. Larroche, “Kinetic simulation of fuel ion transport in ICF target implosions,”
European PhysicalJournal D , vol. 27, pp. 131–146, 2003.[34] O. Larroche, “An efficient explicit numerical scheme for diffusion-type equations with a highly inho-mogeneous and highly anisotropic diffusion tensor,”
J. Comput. Phys. , vol. 223, pp. 436–450, 2007.[35] O. Larroche, “Ion Fokker-Planck simulation of D-3He gas target implosions,”
Phys. Plasmas , vol. 19,p. 122706, 2012.[36] B. E. Peigney, O. Larroche, and V. Tikhonchuk, “Fokker-Planck kinetic modeling of supra thermal α -particles in a fusion plasma,” J. Comput. Phys. , vol. 278, pp. 416–444, 2014.[37] A. Inglebert, B. Canaud, and O. Larroche, “Species separation and modification of neutron diagnosticsin inertial-confinement fusion,”
Euro. Phys. Lett. , vol. 107, p. 65003, 2014.[38] O. Larroche, “Nuclear yield reduction in inertial confinement fusion exploding-pusher targets ex-plained by fuel-pusher mixing through hybrid kinetic-fluid modeling,”
Phys. Rev. E , vol. 87, p. 031201,2018.[39] A. Nishiguchi, K. Mima, H. Azechi, N. Miyanaga, and S. Nakai, “Kinetic effects of electron thermalconduction on implosion hydrodynamics,”
Phys. Plasmas , vol. 4, no. 417, pp. 417–422, 1992.[40] M. Honda, A. Nishiguchi, H. Takabe, H. Azechi, and K. Mima, “Kinetic effects on the electron thermaltransport in ignition target design,”
Phys. Plasmas , vol. 3, p. 3420, 1996.[41] K. Mima, M. Honda, S. Miyamoto, and S. Kato, “Effects of nonlocal heat transport on laser implosion,”
AIP Conference Proceeding , vol. 369, no. 179, pp. 179–185, 1996.[42] M. Honda, A. Nishiguchi, H. Takabe, H. Azechi, and K. Mima, “Effects of non-local electron thermaltransport on ablative Rayleigh-Taylor instability,”
Phys. Plasmas , vol. 3, p. 3420, 1996.[43] W. T. Taitano, L. Chacón, A. N. Simakov, and K. Molvig, “A mass, momentum, and energy conserving,fully implicit, scalable algorithm for the multi-dimensional, multi-species Rosenbluth-Fokker-Planckequation,”
J. Comput. Phys. , vol. 297, pp. 357–380, 2015.[44] W. T. Taitano, L. Chacón, and A. N. Simakov, “An adaptive, conservative 0D-2V multispeciesRosenbluth-Fokker-Planck solver for arbitrarily disparate mass and temperature regimes,”
J. Comput.Phys. , vol. 318, pp. 391–420, 2016.[45] W. T. Taitano, L. Chacón, and A. N. Simakov, “An equilibrium-preserving discretization for the non-linear Fokker-Planck operator in arbitrary multi-dimensional geometry,”
J. Comput. Phys. , vol. 339,pp. 453–460, 2017. 3546] W. T. Taitano, A. N. Simakov, L. Chacón, and B. Keenan, “Yield degradation in inertial-confinement-fusion implosions due to shock-driven kinetic fuel-species stratification and viscous heating,”
Phys.Plasmas , vol. 25, p. 056310, 2018.[47] W. T. Taitano, L. Chacón, and A. N. Simakov, “An adaptive, implicit, conservative 1D-2V multi-species Vlasov-Fokker-Planck multiscale solver in planar geometry,”
J. Comput. Phys. , vol. 365,pp. 173–205, 2018.[48] B. D. Keenan, W. Taitano, and K. Molvig, “Physics of the implosion up until the time of ignition in arevolver (triple-shell) capsule,”
Phys. Plasmas , vol. 27, p. 042704, 2020.[49] B. D. Keenan, W. Taitano, A. N. Simakov, L. Chacón, and B. Albright, “Shock-driven kinetic anddiffusive mix in high-z pusher ICF designs,”
Phys. Plasmas , vol. 27, p. 022704, 2020.[50] W. Taitano, L. Chacón, A. Simakov, and S. Anderson, “A conservative phase-space moving-grid strat-egy for a 1d-2v vlasov-fokker-planck solver,”
Comp. Phys. Comm. , p. in review, 2020.[51] S. Anderson, W. Taitano, L. Chacón, and A. Simako, “An efficient, conservative, time-implicit solverfor the fully kinetic arbitrary-species 1D-2V Vlasov-ampère system,”
J. Comp. Phys. , vol. 419,p. 109686, 2020.[52] W. Huang, Y. Ren, and R. D. Russell, “Moving mesh methods based on moving mesh partial differen-tial equations,”
J. Comput. Phys. , vol. 113, no. 2, pp. 279–290, 1994.[53] W. Huang, Y. Ren, and R. D. Russell, “Moving mesh partial differential equations (MMPDEs) basedon the equidistribution principle,”
SIAM J. Numer. Anal. , vol. 31, no. 3, pp. 709–730, 1994.[54] C. J. Budd, W. Huang, and R. D. Russell, “Adaptivity with moving grids,”
Acta Numerica , vol. 18,pp. 111–241, 2009.[55] S. Li, L. Petzold, and Y. Ren, “Stability of moving mesh systems of partial differential equations,”
SIAM J. Sci. Comput. , vol. 20, no. 2, pp. 719–738, 1998.[56] W. T. Taitano, L. Chacón, A. N. Simakov, and K. Molvig, “A mass, momentum, and energy conserving,fully implicit, scalable algorithm for the multi-dimensional, multi-species rosenbluth-fokker-planckequation,”
J. Comput. Phys. , vol. 297, pp. 357–380, 2015.[57] A. N. Simakov and K. Molvig, “Electron transport in a collisional plasma with multiple ion species,”
Phys. Plasmas , vol. 21, p. 024503, 2014.[58] A. Winslow, “Numerical solution of the quasi-linear Poisson equation in a nonuniform triangle mesh,”
J. Comput. Phys. , vol. 1, p. 149, 1967.[59] L. Chacón, G. Delzanno, and J. Finn, “Robust, multidimensional mesh-motion based on monge-kantorovich equidistribution,”
J. Comp. Phys. , vol. 230, pp. 87–103, 2011.[60] G. D. Byrne and A. C. Hindmarsh, “A polyalgorithm for the numerical solution of ordinary differentialequations,”
ACM Transactions on Mathematical Software , vol. 1, no. 1, pp. 71–96, 1975.[61] E. J. D. Toit, M. R. O’Brien, and R. G. L. Vann, “Positivity-preserving scheme for two-dimensionaladvection-diffusion equations including mixed derivatives,”
Comp. Phys. Comm. , vol. 228, pp. 61–68,2018. 3662] P. H. Gaskell and A. K. C. Lau, “Curvature-compensated convective transport: SMART, a newboundedness-preserving transport algorithm,”
Intern. Journ. Num. Meth. in Fluids , vol. 8, pp. 617–641, 1988.[63] G. Vogman, U. Shumlak, and P. Colella, “Conservative fourth-order finite-volume Vlasov-Poissonsolver for axisymmetric plasmas in cylindrical ( r , v r , v θ ) phase space coordinates,” J. Comp. Phys. ,vol. 373, pp. 877–899, 2018.[64] M. Gittings, R. Weaver, M. Clover, T. Betlach, N. Byrne, R. Coker, E. Dendy, R. Hueckstaedt,K. New, W. Oakes, D. Ranta, and R. Stefan, “The RAGE radiation-hydrodynamic code,” eprint ,p. arXiv:1903.05467, 2008.[65] M. M. Marinak, R. E. Tiption, O. L. Landen, T. J. Murphy, P. Amendt, S. W. Haan, S. P. Hatchett, C. J.Keane, R. McEachem, and R. Wallace, “Three-dimensional simulations of Nova high growth factorcapsule implosion experiments,”
Phys. Plasmas , vol. 3, p. 2070, 1996.[66] J. Delettrez and E. Goldman, “Laboratory for laser energetics report no. 36,” 1976.[67] D. G. Anderson, “Iterative procedures for nonlinear integral equations,”
J. Assoc. Comput. Mach. ,vol. 12, pp. 547 – 560, 1965.[68] A. Toth and C. Kelley, “Convergence analysis for anderson acceleration,”
SIAM J. Numer. Anal. ,vol. 53, no. 2, pp. 805–819, 2015.[69] J. D. Huba,
NRL Plasma Formulary . NRL/PU/6790-98-358, Washington, DC: Naval Research Labo-ratory, 1998.[70] V. G. Guderley, “Strake kugelige und zylindrische verdichtungsstöe in der nähe des kugelmittelpunktesbzw. des zylinderachse,”
Luftfahrtforshung , vol. 19, pp. pp. 302–312, 1942.[71] R. B. Lazarus, “Self-similar solutions for converging shocks and collapsing cavities,”
SIAM J. Numer.Anal. , vol. 18, no. 2, pp. 316–371, 1981.[72] A. Vallet, X. Ribeyre, and V. Tikhonchuk, “Finite Mach number spherical shock wave, application toshock ignition,”
Phys. Plasmas , vol. 20, p. 082702, 2013.[73] S. D. Ramsey and J. F. Lilieholm, “Verification assessment of piston boundary conditions for la-grangian simulation of the guderley problem,”
J. Verif. Valid. and Uncert. Quant. , vol. 2, p. 031001,2017.[74] R. Kirkpatrick, C. C. Cremer, L. C. Madsen, H. H. Rogers, and R. S. Cooper, “Structured fusion targetdesigns,”
Plasma Phys. Control. Fusion , vol. 15, p. 333, 1975.[75] R. Kirkpatrick and J. A. Wheeler, “The physics of DT ignition in small fusion targets,”
Plasma Phys.Control. Fusion , vol. 21, p. 333, 1981.[76] S. A. Colgate, “Minimum fusion ignition conditions,”
Los Alamos Laboratory Report No. LAUR-88-1268 , 1988.[77] S. A. Colgate, A. G. Petschek, and R. C. Kirkpatrick, “Minimum fusion ignition conditions,”
LosAlamos Laboratory Report No. LAUR-92-2599 , 1992.[78] K. Lackner, S. Colgate, N. Johnson, R. Kirkpatrick, R. Menikoff, and A. Petschek, “Equilibriumignition for icf capsules,”
AIP Conf. Proc. , vol. 318, p. 356, 1994.3779] P. Amendt, J. D. Colvin, R. E. Tipton, D. E. Hinkel, M. J. Edwards, O. L. Landen, J. D. Ramshaw,L. J. Suter, W. S. Varnum, and R. G. Watt, “Indirect-drive noncryogenic double-shell ignition targetsfor the national ignition facility: Design and analysis,”
Phys. Plasmas , vol. 9, p. 2221, 2002.[80] P. Amendt, C. Cerjan, A. Hamza, D. Hinkel, J. Milovich, and H. Robey, “Assessing the prospectsfor achieving double-shell ignition on the national ignition facility using vacuum hohlraums,”
Phys.Plasmas , vol. 14, p. 056312, 2007.[81] D. S. Montgomery, W. S. Daughton, B. J. Albright, A. N. Simakov, D. C. Wilson, E. S. Dodd, R. C.Kirkpatrick, R. G. Watt, M. A. Gunderson, E. N. Loomis, E. C. Merritt, and et al., “Design considera-tions for indirectly driven double shell capsules,”
Phys. Plasmas , vol. 25, p. 092706, 2018.[82] K. Molvig, M. J. Schmitt, B. Albright, E. Dodd, N. Hoffman, G. McCall, and S. Ramsey, “Low fuelconvergence path to direct-drive fusion ignition,”
Phys. Rev. Lett. , vol. 116, p. 255003, 2016.[83] K. Molvig, M. Schmitt, R. Betti, E. Campbell, and P. McKenty, “Stable and confined burn in a revolverignition capsule,”
Phys. Plasmas , vol. 25, p. 082708, 2018.[84] M. V. Dyke and A. J. Guttmann, “The converging shock wave from a spherical or cylindrical piston,”
Phys. J. Fluid Mech. , vol. 120, pp. pp. 451–462, 1982.[85] H. S. Bosch and G. M. Hale, “Improved formulas for fusion cross sections and thermal reactivities,”
Nucl. Fusion , vol. 32, no. 4, pp. 611–631, 1992.[86] K. S. Smith and J. D. R. III, “Full-core, 2-D LWR core calculations with CSMO-4E,” in
Proceedingsof International Conference on New Frontiers of Nuclear Technology: Reactor Physics, Safety andHigh-Performance Computing (PHYSOR) , (Seoul, Korea), 2002.[87] H. Park, D. A. Knoll, R. M. Rauenzahn, C. K. Newman, J. D. Densmore, and A. B. Wollaber, “An ef-ficient and time accurate, moment-based scale-bridging algorithm for thermal radiative transfer prob-lems,”
SIAM Journal of Scientific Computing , vol. 35, no. 5, pp. S18–S41, 2013.[88] E. Aristova and D. Baydin, “Implementation of the quasi diffusion method for calculating the criticalparameters of a fast neutron reactor in 3D hexagonal geometry,”
Math. Models. Comput. Simul. , vol. 5,pp. 145–155, 2013.[89] W. T. Taitano and L. Chacón, “Charge-and-energy conserving moment-based accelerator for a multi-species Vlasov-Fokker-Planck-Ampère system, part I: Collisionless aspects,”
J. Comput. Phys. ,vol. 284, pp. 718–736, 2015.[90] L. Chacón, G. chen, D. Knoll, C. Newman, H. Park, W. Taitano, J. Willert, and G. Womeldorff,“Multiscale high-order/low-order (HOLO) algorithms and applications,”
J. Comp. Phys. , vol. 330,pp. 21–45, 2017.[91] M. Murakami, J. Sanz, and Y. Imamoto, “Stability of spherical converging shock wave,”
Phys. Plas-mas , vol. 22, p. 072703, 2015.
Appendix A. Details on the Fluid Electron Model
The friction between the α -ion species and electrons is modeled as: F α e = − m e n e ν e α ( u α − (cid:104) u α (cid:105) ) + α m e n e ν e α ( u e − (cid:104) u α (cid:105) ) + β n e ν e α ∇ x T e ∑ N s α ν e α , (A.1)38here ν e α = n e e Λ e α ε m / e ( π T e ) / , (A.2)is the collision frequency between electrons and the α -ion species (cid:104) u α (cid:105) = ∑ N s α ν e α u α ∑ N s α ν e α , (A.3)is the collision frequency averaged drift velocity, α = (cid:16) Z e f f + √ Z e f f + (cid:17) Z e f f + √ Z e f f + , (A.4) β = Z e f f (cid:16) Z e f f + √ (cid:17) Z e f f + √ Z e f f + , (A.5)and the effective charge is defined as: Z e f f = − ∑ N s α q α n α q e n e . (A.6)The electron heat flux is given as: Q e = β n e T e ( u e − (cid:104) u α (cid:105) ) − κ e ∇ x T e , (A.7)where the electron-thermal conductivity is given as: κ e = γ n e T e m e ∑ N s α ν e α , (A.8)with γ = Z e f f (cid:16) Z e f f + √ (cid:17) (cid:16) Z e f f + √ Z e f f + (cid:17) . (A.9)The generalized Ohm’s law for an electrostatic plasma is given as: E = ∑ N s α F α e + ∇ x P e q e n e , (A.10)and the electron-ion energy exchange as: W e α = − F α e · u α + ν e α m e m α n e ( T e − T α ) . (A.11) Appendix B. Cartesian-Cartesian to Spherical-Cylindrical Phase-Space Coordinate Transformation
We derive the Vlasov equation in a 1D spherically symmetric configuration space, ( r ) , and a 2V cylin-drically symmetric velocity space, (cid:0) v || , v ⊥ (cid:1) , from the 3D3V Cartesian coordinate system. An illustration ofthe spherical-cylindrical hybrid coordinate system is illustrated in Fig. B.1. We define a map to transformfrom a Cartesian to spherical coordinate system in the configuration space as:39
The contravariant component of the advection velocity is derived for individual components, (cid:8) ˙ t , ˙ r , ˙ θ , ˙ φ , ˙ v || , ˙ v ⊥ , ˙ φ v (cid:9) .The ˙ t component is trivially evaluated to unity, as t depends neither on space nor velocity.The ˙ r is also trivially evaluated as: ˙ r = v · ∂ r ∂ r = v · (cid:98) r = v || . (B.30)The ˙ θ is evaluated starting from: v · ∂ θ∂ r = v · (cid:98) θθθ | ∂ r θ | . (B.31)Here, v · (cid:98) θθθ = v (cid:16) cos θ v (cid:98) r + sin θ v cos φ v (cid:98) θθθ + sin θ v sin φ v (cid:98) φφφ (cid:17) · (cid:98) θθθ = v sin θ v cos φ v = v ⊥ cos φ v , (B.32)and (cid:12)(cid:12)(cid:12)(cid:12) ∂ θ∂ r (cid:12)(cid:12)(cid:12)(cid:12) = r . (B.33)Thus, putting the two terms together, we obtain˙ θ = v · ∂ θ∂ r = v ⊥ cos φ v r . (B.34)The ˙ φ is evaluated similarly from, ˙ φ = v · ∂ φ∂ r = v · ˆ φφφ | ∂ r φ | . (B.35)Here, v · (cid:98) φφφ = v (cid:16) cos θ v (cid:98) r + sin θ v cos φ v (cid:98) θθθ + sin θ v sin φ v (cid:98) φφφ (cid:17) · (cid:98) φφφ = v sin θ v sin φ v = v ⊥ sin φ v , (B.36)and (cid:12)(cid:12)(cid:12)(cid:12) ∂ φ∂ r (cid:12)(cid:12)(cid:12)(cid:12) = r sin θ . (B.37)Thus, putting the two terms together, we obtain˙ φ = v ⊥ r sin φ v sin θ . (B.38)The ˙ v || is evaluated from, ˙ v || = v · ∂ v || ∂ r + a · ∂ v || ∂ v . (B.39)Since v || = v · (cid:98) r = v · r √ x + y + z , v · ∂ v || ∂ r = v · ∂∂ r (cid:16) v · r r (cid:17) = v · (cid:40) v · (cid:34) r ↔ I − r ⊗ r r (cid:35)(cid:41) . (B.40)43ith r = { x , y , z } , r = (cid:112) x + y + z , and completing the square, we obtain v · (cid:40) v · (cid:34) r ↔ I − r ⊗ r r (cid:35)(cid:41) = v − ( v · (cid:98) r ) r = v − v || r . (B.41)Since v ⊥ = v − v || , ˙ v || = v ⊥ r + a · ∂ v || ∂ v . (B.42)The ˙ v ⊥ is evaluated from, ˙ v = v · ∂ v ⊥ ∂ r + a · ∂ v ⊥ ∂ v . (B.43)Since v ⊥ = (cid:113) v − v || , v · ∂ v ⊥ ∂ r = v · ∂∂ r (cid:113) v − v || = v · − v || ∂ r v || (cid:113) v − v || = − v || v ⊥ v · (cid:40) v · (cid:34) r ↔ I − r ⊗ r r (cid:35)(cid:41) . (B.44)Thus, − v || v ⊥ v · (cid:40) v · (cid:34) r ↔ I − r ⊗ r r (cid:35)(cid:41) = − v || v ⊥ r (B.45)and finally ˙ v ⊥ = − v || v ⊥ r + a · ∂ v ⊥ ∂ v . (B.46)The ˙ φ v is obtained from: ˙ φ v = v · ∂ φ v ∂ r + a · ∂ φ v ∂ v , (B.47)where given φ v = tan − ( v / v ) , we obtain ∂ φ v ∂ r = v + v (cid:18) v ∂ v ∂ r − v ∂ v ∂ r (cid:19) . (B.48)From, v = v · (cid:98) θθθ , and v = v · (cid:98) φφφ , (B.49) ∂ v ∂ r = v · ∂ (cid:98) θθθ∂ r , and ∂ v ∂ r = v · ∂ (cid:98) φφφ∂ r , (B.50)and v = ( v sin θ cos φ + v cos θ cos φ − v sin φ ) (cid:98) x + ( v sin θ sin φ + v cos θ sin φ + v cos φ ) (cid:98) y + ( v cos θ − v sin θ ) (cid:98) z , (B.51)one can show, v · ∂ (cid:98) θθθ∂ r = − r sin θ (cid:2)(cid:0) v cos φ sin θ + v cos θ sin φ (cid:1)(cid:98) x ++ (cid:18) − v cos θ cos φ + v cos θ sin φ + v (cid:0) + sin θ sin φ (cid:1)(cid:19)(cid:98) y + v sin θ cos θ (cid:98) z (cid:21) (B.52)44nd v · ∂ (cid:98) φφφ∂ r = − v r sin θ ( cos φ (cid:98) x + sin φ (cid:98) y ) . (B.53)We can then trivially obtain: v · (cid:34) v · ∂ (cid:98) φφφ∂ r (cid:35) = − r (cid:104) v v + v v tan θ (cid:105) , and v · (cid:34) v · ∂ (cid:98) θθθ∂ r (cid:35) = − r (cid:20) v v − v tan θ (cid:21) . (B.54)Assembling the terms, v · ∂ φ v ∂ r = v + v (cid:34) v v · v · ∂ (cid:98) φφφ∂ r − v v · v · ∂ (cid:98) θθθ∂ r (cid:35) = − r (cid:0) v + v (cid:1) (cid:20) v (cid:16) v v + v v tan θ (cid:17) − v (cid:18) v v − v tan θ (cid:19)(cid:21) = − v r θ . (B.55)Finally, since v = v ⊥ sin φ v , ˙ φ v = − v ⊥ r sin φ v tan θ + a · ∂ φ v ∂ v . (B.56) Appendix C. Physical to Logical Phase-Space Coordinate Transformation
Starting from our 1D-2V spherically symmetric in configuration space and azimuthally symmetric cylin-drical velocity-space coordinate system, T = (cid:8) t , r , v || , v ⊥ (cid:9) , we transform our Vlasov equation, Eq. (B.28), ∂ t ( J Sv ⊥ f ) + ∂ r (cid:0) J Sv ⊥ v || f (cid:1) + qm E || ∂ v || ( J Sv ⊥ f ) + r ∂ v · ( J Sv ⊥ (cid:101) a f ) = , (C.1)(where (cid:101) a = (cid:2) v ⊥ , − v || v ⊥ (cid:3) T is the acceleration vector associated with the inertial term due to the coordinatetransformation) to logical in configuration and normalized in velocity coordinate system, τ = (cid:8) t , ξ , (cid:101) v || , (cid:101) v ⊥ (cid:9) .We begin by rewriting the equation as: ∇ T · (cid:0) ˙T ¯ f (cid:1) ≡ ∂ t (cid:0) ¯ f (cid:1) + ∂ r (cid:0) v || ¯ f (cid:1) + ∂ v || (cid:20) qm E || + v ⊥ r (cid:21)(cid:124) (cid:123)(cid:122) (cid:125) ¯ a || ¯ f + ∂ v ⊥ − v || v ⊥ r (cid:124) (cid:123)(cid:122) (cid:125) ¯ a ⊥ ¯ f = , (C.2)where ¯ f = J Sv ⊥ f and ˙T = (cid:8) , v || , ¯ a || , ¯ a ⊥ (cid:9) . From steps similar to those in Appendix B, we obtain: ∇ T · (cid:0) ˙T ¯ f (cid:1) = J ξ (cid:98) v ∂∂ τ i (cid:0) J ξ (cid:98) v ˙ τ i ¯ f (cid:1) = , (C.3)where, J ξ (cid:98) v = ∂ r ∂ξ v ∗ is the Jacobian of transformation for our system. By multiplying by J ξ (cid:98) v , we obtain: ∂∂ τ i (cid:0) J ξ (cid:98) v ˙ τ i ¯ f (cid:1) = . (C.4)The advection velocity in the transformed coordinate system is given by:˙ τ i = ˙T · g i = ∂ τ i ∂ t + v || ∂ τ i ∂ r + ¯ a || ∂ τ i ∂ v || + ¯ a ⊥ ∂ τ i ∂ v ⊥ , (C.5)45nd the components are given as: ˙ τττ = ∂ξ∂ t + v || ∂ξ∂ r ∂ (cid:101) v || ∂ t + v || ∂ (cid:101) v || ∂ r + ¯ a || ∂ (cid:101) v || ∂ v || ∂ (cid:101) v ⊥ ∂ t + v || ∂ (cid:101) v ⊥ ∂ r + ¯ a ⊥ ∂ (cid:101) v ⊥ ∂ v ⊥ . (C.6)We derive the expression for each terms individually, beginning with the temporal derivatives. Recallthat ∂∂ t ≡ (cid:18) ∂∂ t (cid:19) r , v || , v ⊥ (C.7)and that (cid:18) ∂ τ i (cid:54) = t ∂ t (cid:19) T j (cid:54) = t = − (cid:18) ∂ T ∂ t (cid:19) τ k (cid:54) = t · g i (cid:54) = t . (C.8)We obtain, ( ∂ t t ) r , v || , v ⊥ ( ∂ t ξ ) r , v || , v ⊥ (cid:0) ∂ t (cid:101) v || (cid:1) r , v || , v ⊥ ( ∂ t (cid:101) v ⊥ ) r , v || , v ⊥ = − (cid:16) ∂ r ∂ t (cid:17) τττ k (cid:54) = t (cid:16) ∂ r ∂ξ (cid:17) − τττ k (cid:54) = ξ − (cid:20)(cid:16) ∂ v || ∂ t (cid:17) τττ k (cid:54) = t − (cid:16) ∂ r ∂ t (cid:17) τττ k (cid:54) = t (cid:16) ∂ r ∂ξ (cid:17) − τττ k (cid:54) = ξ (cid:16) ∂ v || ∂ξ (cid:17) τττ k (cid:54) = ξ (cid:21) (cid:16) ∂ v || ∂ (cid:101) v || (cid:17) − τττ k (cid:54) = (cid:101) v || − (cid:20)(cid:16) ∂ v ⊥ ∂ t (cid:17) τττ k (cid:54) = t − (cid:16) ∂ r ∂ t (cid:17) τττ k (cid:54) = t (cid:16) ∂ r ∂ξ (cid:17) − τττ k (cid:54) = ξ (cid:16) ∂ v ⊥ ∂ξ (cid:17) τττ k (cid:54) = ξ (cid:21) (cid:16) ∂ v ⊥ ∂ (cid:101) v ⊥ (cid:17) − τττ k (cid:54) = (cid:101) v ⊥ . (C.9)Here, recall that v || = (cid:101) v || v ∗ + u ∗|| and v ⊥ = (cid:101) v ⊥ v ∗ , (C.10)therefore (cid:18) ∂ v || ∂ t (cid:19) τττ k (cid:54) = t = (cid:101) v || ∂ v ∗ ∂ t + ∂ u ∗ ∂ t , (C.11) (cid:18) ∂ v ⊥ ∂ t (cid:19) τττ k (cid:54) = t = (cid:101) v ⊥ ∂ v ∗ ∂ t , (C.12) (cid:18) ∂ v || ∂ (cid:101) v || (cid:19) τττ k (cid:54) = (cid:101) v || = (cid:18) ∂ v ⊥ ∂ (cid:101) v ⊥ (cid:19) τττ k (cid:54) = (cid:101) v ⊥ = v ∗ . (C.13)The rest of the terms are trivially evaluated as: v || ∂ ξ∂ r = v ∗ (cid:16)(cid:101) v || + (cid:98) u ∗|| (cid:17) ∂ ξ∂ r , (C.14) ∂ (cid:101) v || ∂ r = ∂ (cid:101) v || ∂ v ∗ ∂ v ∗ ∂ ξ ∂ ξ∂ r + ∂ (cid:101) v || ∂ u ∗|| ∂ u ∗|| ∂ ξ ∂ ξ∂ r , (C.15) ∂ (cid:101) v ⊥ ∂ r = ∂ (cid:101) v || ∂ v ∗ ∂ v ∗ ∂ ξ ∂ ξ∂ r , (C.16) ∂ (cid:101) v || ∂ v || = v ∗ , and ∂ (cid:101) v ⊥ ∂ v ⊥ = v ∗ , (C.17)46 (cid:101) v || ∂ v ∗ = − v || − u ∗|| v ∗ = − (cid:101) v || − (cid:98) u ∗|| v ∗ , ∂ (cid:101) v || ∂ u ∗|| = v ∗ , (C.18)and ∂ (cid:101) v ⊥ ∂ v ∗ = − v ⊥ v ∗ = − (cid:98) v ⊥ v ∗ . (C.19)Since ¯ f = J Sv ⊥ f , (cid:101) f = v ∗ f , by defining J S (cid:101) v ⊥ = r (cid:101) v ⊥ , J ξ = ∂ r ∂ξ ; and noting that (cid:101) v ∂ v ∗ ∂ t + ∂ u ∗ ∂ t = ∂∂ t (cid:0) v ∗ (cid:101) v + u ∗ e || (cid:1) = ∂ v ∂ t , (cid:101) v ∂ v ∗ ∂ξ + ∂ u ∗ ∂ξ = ∂∂ξ (cid:0) v ∗ (cid:101) v + u ∗ e || (cid:1) = ∂ v ∂ t , where e || is the unit vector in the parallel-velocity direction; v || = v ∗ (cid:101) v || + u ∗|| , v ⊥ = v ∗ (cid:101) v ⊥ ; and substituting the above results into the Eq. C.3, the following transformed equationfor (cid:101) f ( ξ , (cid:101) v , t ) is obtained: ∂ (cid:16) J S (cid:101) v ⊥ J ξ (cid:101) f (cid:17) ∂ t + ∂∂ ξ (cid:104) J S (cid:101) v ⊥ (cid:0) v || − ∂ t r (cid:1) (cid:101) f (cid:105) − v ∗ ∂∂ (cid:101) v · (cid:20) ∂ v ∂ t J S (cid:101) v ⊥ J ξ (cid:101) f (cid:21) − v ∗ ∂∂ (cid:101) v · (cid:20) ∂ v ∂ ξ (cid:0) v || − ∂ t r (cid:1) J S (cid:101) v ⊥ (cid:101) f (cid:21) + v ∗ ∂∂ (cid:101) v || (cid:16) qm E || J S (cid:101) v ⊥ J ξ (cid:101) f (cid:17) − v ∗ ∂∂ (cid:101) v · (cid:20) J S (cid:101) v ⊥ J ξ (cid:101) a r (cid:101) f (cid:21) . (C.20)Finally by defining ˙ J r = ∂ r ∂ t = J S ∂ r ∂ t , J S ξ = J S J ξ , and dividing by (cid:101) v ⊥ , we obtain: ∂ (cid:16) J S ξ (cid:101) f (cid:17) ∂ t + ∂∂ ξ (cid:104)(cid:0) J S v || − ˙ J r (cid:1) (cid:101) f (cid:105) − v ∗ ∂∂ (cid:101) v · (cid:20) ∂ v ∂ t J S ξ (cid:101) f (cid:21) − v ∗ ∂∂ (cid:101) v · (cid:20) ∂ v ∂ ξ (cid:0) J S v || − ˙ J r (cid:1) (cid:101) f α (cid:21) + v ∗ ∂∂ (cid:101) v || (cid:16) J S ξ qm E || (cid:101) f (cid:17) + v ∗ ∂∂ (cid:101) v · (cid:20) J S ξ (cid:101) a r (cid:101) f (cid:21) = . (C.21)Here, ∂ (cid:101) v · A = (cid:104) ∂ (cid:101) v || (cid:16) A v || (cid:17) , (cid:101) v − ⊥ ∂ (cid:101) v ⊥ ( (cid:101) v ⊥ A v ⊥ ) (cid:105) T is the velocity-space divergence operator acting on a vector A = (cid:104) A v || , A v ⊥ (cid:105) T . Appendix D. Vlasov Conservation Symmetries for the Transformed Coordinate System
Following similar procedures outlined in Ref. [50], we derive the conservation symmetries for theVlasov equation in a 1D spherically symmetric configuration space and 2V cylindrically symmetric velocity-space coordinate system. We note that the conservation properties associated with the temporal and spatialvariation of velocity-space metrics, as well as the inertial terms introduced by the spherical configurationspace transformation can all be shown independently. We begin by developing discretization for exact mass,linear momentum, and energy conservation in a spatially homogeneous system (0D), and then a spatiallyinhomogeneous case (1D) in a periodic spatial domain without any background field (conservation with theelectric field can be shown separately, as demonstrated in Ref. [47]).
Appendix D.0.1. Temporal variation of v ∗ and (cid:98) u ∗|| For simplicity, we drop the subscript denoting ion species. By consider only the terms associated withthe temporally varying velocity-space metrics in Eq. (3.4), we obtain the following simplified form of theVlasov equation: ∂ t (cid:16) J r ξ (cid:101) f (cid:17) − J r ξ v ∗ ∂∂ (cid:101) v · (cid:16) (cid:101) f ∂ t v (cid:17) = , (D.1)47here ∂ t v = ∂ t (cid:104) v ∗ (cid:16)(cid:101) v + (cid:98) u ∗|| e || (cid:17)(cid:105) . In the continuum, mass conservation is defined as: (cid:90) d ξ (cid:68) m , ∂ t (cid:16) J r ξ (cid:101) f (cid:17)(cid:69) (cid:101) v = . (D.2)This can be shown trivially due to the divergence form of the inertial terms.Momentum and energy conservation is defined by following Ref. [50] and using the chain rule andtaking the (cid:10) mv || , ( · ) (cid:11) (cid:101) v moment of Eq. (D.1) as: (cid:42) , m ∂ t (cid:16) v || J r ξ (cid:101) f (cid:17) − (cid:20) J r ξ (cid:101) f ∂ t v || + J r ξ v || v ∗ ∂∂ (cid:101) v · (cid:16) ∂ t v (cid:101) f (cid:17)(cid:21)(cid:124) (cid:123)(cid:122) (cid:125) a (cid:13) (cid:43) (cid:101) v = , and (cid:42) , m ∂ t (cid:18) v J r ξ (cid:101) f (cid:19) − (cid:34) J r ξ (cid:101) f ∂ t v + J r ξ v v ∗ ∂∂ (cid:101) v · (cid:16) (cid:101) f ∂ t v (cid:17)(cid:35)(cid:124) (cid:123)(cid:122) (cid:125) b (cid:13) (cid:43) (cid:101) v = . It can be shown by using integration by parts that, (cid:104) m , a (cid:13)(cid:105) (cid:101) v = (cid:104) m , b (cid:13)(cid:105) (cid:101) v =
0, and the key is to ensurethese symmetries discretely. These symmetries are ensure by multiplying the inertial term by a velocity-space dependent function, γ t ( v ) , such that, at a time-step p , (cid:42) v ( p ) || , i , δ (cid:101) v · γ ( p + ) t , i J ( p + ) t , i (cid:124) (cid:123)(cid:122) (cid:125) c (cid:13) (cid:43) δ (cid:101) v − (cid:20)(cid:28) v ( p ) || , i , δ t (cid:16) J r ξ (cid:101) f (cid:17) ( p + ) i (cid:29) δ (cid:101) v − (cid:28) , δ t (cid:16) v || J r ξ (cid:101) f (cid:17) ( p + ) i (cid:29) δ (cid:101) v (cid:21) = (cid:42) (cid:104) v ( p ) i (cid:105) , δ (cid:101) v · (cid:16) γ ( p + ) t , i J ( p + ) t , i (cid:17)(cid:43) δ (cid:101) v − (cid:42) (cid:104) v ( p ) i (cid:105) , δ t (cid:16) J r ξ (cid:101) f (cid:17) ( p + ) i (cid:43) δ (cid:101) v − (cid:42) , δ t (cid:18) v J r ξ (cid:101) f (cid:19) ( p + ) i (cid:43) δ (cid:101) v = . (D.4)Here, i is the spatial cell index (velocity-space indices are dropped for brevity), δ (cid:101) v · is the discrete velocity-space divergence operator, δ t is the discrete temporal derivative, γ t ( (cid:101) v ) = + O (cid:16) ∆ β v , ∆ ζ t (cid:17) is the discrete-nonlinear-constraint function [47], where β and ζ are the velocity-space and temporal discretization trun-cation order, the term c (cid:13) is the discrete flux for the inertial term in Eq. 5.4 , due to the temporal variation ofthe velocity-space transformation metrics [details shown in Eqs. (E.4)-(E.7)], and v ( p ) i = v ∗ , ( p ) i (cid:104)(cid:101) v + (cid:98) u ∗ , ( p ) || e || (cid:105) , t v ( p + ) i = c ( p + ) v ( p ) i + c ( p ) v ( p − ) i + c ( p − ) v ( p − ) i ∆ t ( p ) , δ t (cid:16) J S ξ (cid:101) f (cid:17) ( p + ) i = c ( p + ) J ( p + ) S ξ , i (cid:101) f ( p + ) i + c ( p ) J ( p ) S ξ , i (cid:101) f ( p ) i + c ( p − ) J ( p − ) S ξ , i (cid:101) f ( p − ) i ∆ t ( p ) , δ t (cid:16) v || J S ξ (cid:101) f (cid:17) ( p + ) i = c ( p + ) v ( p ) || , i J ( p + ) S ξ , i (cid:101) f ( p + ) i + c ( p ) v ( p − ) || , i J ( p ) S ξ , i (cid:101) f ( p ) i + c ( p − ) v ( p − ) || , i J ( p − ) S ξ , i (cid:101) f ( p − ) i ∆ t ( p ) , δ t (cid:18) v J S ξ (cid:101) f (cid:19) ( p + ) i = c ( p + ) (cid:104) v ( p ) i (cid:105) J ( p + ) S ξ , i (cid:101) f ( p + ) i + c ( p ) (cid:104) v ( p − ) i (cid:105) J ( p ) S ξ , i (cid:101) f ( p ) i + c ( p − ) (cid:104) v ( p − ) i (cid:105) J ( p − ) S ξ , i (cid:101) f ( p − ) i ∆ t ( p ) . To evaluate γ t , we follow Refs. [47, 50] and begin by assuming a velocity space dependent local (in config-uration space) functional representation: γ t (cid:0)(cid:101) v || , (cid:101) v ⊥ (cid:1) = + P ∑ l = C l B l (cid:0)(cid:101) v || , (cid:101) v ⊥ (cid:1) , (D.5)where P ∑ l = C l B l (cid:0)(cid:101) v || , (cid:101) v ⊥ (cid:1) = P || ∑ r = P ⊥ ∑ s = C rs B || , r (cid:0)(cid:101) v || (cid:1) B ⊥ , s ( (cid:101) v ⊥ ) , B || , r is a r th representation in the parallel velocity component, B ⊥ , s is a similar quantity in the perpendicularvelocity component, and C rs is the coefficient corresponding to the respective functions. In this study, wechose a Fourier representation where: B || , r = r = (cid:104) rk || (cid:16)(cid:101) v || + (cid:98) u ∗|| − u || / v ∗ (cid:17)(cid:105) if mod ( r , ) = (cid:104) ( r − ) k || (cid:16)(cid:101) v || + (cid:98) u ∗|| − u || / v ∗ (cid:17)(cid:105) if mod ( r , ) = , B ⊥ , s = s = [ sk ⊥ (cid:101) v ⊥ ] if mod ( s , ) = [( s − ) k ⊥ (cid:101) v ⊥ ] if mod ( s , ) = , and k || = π / L || , k ⊥ = π / (cid:101) L ⊥ are the wave vectors. We also choose r = s = ( , , ) . The coefficients, C l ,are obtained by minimizing their amplitude while satisfying the discrete symmetry constraints as given byEqs. (D.3) and (D.4). This is done by solving a constrained-minimization problem for the following costfunction: F ( C , λλλ ) = P ∑ l = C l − λλλ T · M . (D.6)Here, C = [ C , C , · · · , C P ] T , λλλ is a vector of Lagrange multipliers and M is the vector of constraints [Eqs.(D.3) and (D.4)]; and C is obtained from the linear system: (cid:20) ∂ CCC F ∂ λλλ F (cid:21) = . (D.7)We note that, similarly to previous studies employing discrete nonlinear constraints [56, 89, 43, 44, 47],since γ t is an implicit function of the solution, for nonlinearly implicit system such as ours the quality ofdiscrete conservation properties depends on the prescribed nonlinear convergence tolerance of our solver (asdemonstrated in Sec. 6.1). 49 ppendix D.0.2. Spatial variation of v ∗ and (cid:98) u ∗|| Similarly to the temporal variation, conservation symmetries for the case of spatial variation of thevelocity-space metrics can be shown independently. Consider only the spatial gradient terms in the Vlasovequation, Eq. (3.4), to obtain the following expression: ∂ ξ (cid:16) v || , e f f (cid:101) f (cid:17) (cid:101) v , t − v ∗ ∂∂ (cid:101) v · (cid:16) ∂ ξ v (cid:12)(cid:12) (cid:101) v , t v || , e f f (cid:101) f (cid:17) . (D.8)Here, v || , e f f = J S v ∗ (cid:16)(cid:101) v || + (cid:98) u ∗|| (cid:17) − ˙ J r and ∂ ξ v (cid:12)(cid:12) (cid:101) v , t = ∂ ξ (cid:104) v ∗ (cid:16) v + (cid:98) u ∗|| e || (cid:17)(cid:105) , and the mass conservation theorem isrevealed by taking the mv moment, assuming a periodic boundary condition, and integrating in ξ to find (cid:90) (cid:68) , ∂ ξ (cid:16) mv || , e f f (cid:101) f (cid:17)(cid:69) (cid:101) v d ξ = . (D.9)Note that the inertial term is in a divergence form in the velocity space, and therefore its mv momenttrivially vanishes both in the continuum and discretely.In the transformed coordinate, the momentum and energy conservations are defined as: (cid:90) d ξ (cid:28) v || , ∂ ξ (cid:16) v || , e f f (cid:101) f (cid:17) (cid:101) v , t (cid:29) (cid:101) v − (cid:42) v || , ∂∂ (cid:101) v · (cid:32) ∂ ξ v (cid:12)(cid:12) (cid:101) v , t v ∗ v || , e f f (cid:101) f (cid:33) ξ , t (cid:43) (cid:101) v = . (D.10)and (cid:90) d ξ (cid:26)(cid:28) v , ∂ ξ (cid:16) v || , e f f (cid:101) f (cid:17) (cid:101) v , t (cid:29) (cid:101) v − (cid:28) v , v ∗ ∂∂ (cid:101) v · (cid:16) ∂ ξ v (cid:12)(cid:12) (cid:101) v , t v || , e f f (cid:101) f (cid:17) (cid:101) v , t (cid:29) (cid:101) v (cid:27) = . (D.11)These symmetries are ensured in the discrete by modifying the inertial term by a velocity-space dependentfunction, γ r ( (cid:101) v ) , such that the following relationships are satisfied: N ξ ∑ i = ∆ ξ (cid:110)(cid:68) v ( p ) || , i , δ ξ ( F r + F ˙ r ) ( p + ) i (cid:69) δ (cid:101) v − (cid:42) v ( p ) || , i , δ (cid:101) v · γ ( p + ) r , i + / J − , ( p + ) r (cid:124) (cid:123)(cid:122) (cid:125) a (cid:13) + γ ( p + ) r , i − / J + , ( p + ) r (cid:124) (cid:123)(cid:122) (cid:125) b (cid:13) i (cid:43) δ (cid:101) v = N ξ ∑ i = ∆ ξ (cid:42) (cid:104) v ( p ) i (cid:105) , δ ξ ( F r + F ˙ r ) ( p + ) i (cid:43) δ (cid:101) v − (cid:42) (cid:104) v ( p ) i (cid:105) , δ (cid:101) v · (cid:16) γ ( p + ) r , i + / J − , ( p + ) r + γ ( p + ) r , i − / J + , ( p + ) r (cid:17) i (cid:43) δ (cid:101) v = . (D.13)Here, v ( p ) i = v ∗ , ( p ) i (cid:104)(cid:101) v + (cid:98) u ∗ , ( p ) || , i e || (cid:105) , γ x ( (cid:101) v ) = + O (cid:16) ∆ β v , ∆ η x (cid:17) is the spatial discrete-nonlinear-constraint func-tion, with the functional form defined similarly to γ t [Eq. (D.5)], and terms a (cid:13) and b (cid:13) are the discrete fluxesfor the inertial term in Eq. (5.4) [details shown in Eqs. (E.8)-(E.15)]. Performing a discrete integrationby parts (i.e., telescoping the summation) on Eqs. (D.12) and (D.13), we obtain the following constraints50hat relate γ r , i + / , the discrete configuration-space flux, and the velocity-space inertial terms for momentumconservation: (cid:68) v ( p ) || , i − v ( p ) || , i + , ( F r − F ˙ r ) ( p + ) i + / (cid:69) δ (cid:101) v − ∆ ξ (cid:20)(cid:68) v ( p ) || , i , δ (cid:101) v · (cid:16) γ ( p + ) r , i + / J − , ( p + ) r (cid:17) i (cid:69) δ (cid:101) v + (cid:68) v ( p ) || , i + , δ (cid:101) v · (cid:16) γ ( p + ) r , i + / J + , ( p + ) r (cid:17) i + (cid:69) δ (cid:101) v (cid:21) = , (D.14)and for energy conservation: (cid:42) (cid:104) v ( p ) i (cid:105) − (cid:104) v ( p ) i + (cid:105) , ( F r − F ˙ r ) ( p + ) i + / (cid:43) δ (cid:101) v − ∆ ξ (cid:42) (cid:104) v ( p ) i (cid:105) , δ (cid:101) v · (cid:16) γ ( p + ) r , i + / J − , ( p + ) r (cid:17) i (cid:43) δ (cid:101) v + (cid:42) (cid:104) v ( p ) i + (cid:105) , δ (cid:101) v · (cid:16) γ ( p + ) r , i + / J + , ( p + ) r (cid:17) i + (cid:43) δ (cid:101) v = . (D.15)The vector of coefficients, C , for γ r , i + / is evaluated by solving a constrained minimization problem as inEq. (D.7) with the vector of vanishing constraints, M , being Eqs. (D.14) and (D.15). We end by notingthat, at boundaries we set γ r , i + / = ∂ ξ v ∗ = ∂ ξ (cid:98) u || = ∂ ξ v =
0) as the boundary condition violates thecontinuum conservation principle.
Appendix D.0.3. Spherical inertial term
Consider the inertial term due to the transformation from Cartesian to spherical geometry: v ∗ − r − ∂ (cid:101) v · ( J S (cid:101) a f ) . The mass conservation is trivially ensured due to the divergence form of the operator while conservationof linear momentum is automatically satisfied via spherical symmetry. The energy conservation theorem isshown by, r − v ∗ − (cid:28) v ∂ (cid:101) v · ( J S (cid:101) a f ) (cid:29) v = − J S r − v ∗ − (cid:68) v || , v ⊥ (cid:101) f (cid:69) (cid:101) v (cid:124) (cid:123)(cid:122) (cid:125) a (cid:13) − (cid:68) v ⊥ , v || v ⊥ (cid:101) f (cid:69) (cid:101) v (cid:124) (cid:123)(cid:122) (cid:125) b (cid:13) = . (D.16)In the discrete, we ensure the cancellation between terms a (cid:13) and b (cid:13) by modifying the parallel-velocitycomponent of spherical geometry inertial term by a constant, γ S , such that: (cid:42) (cid:104) v ( p ) i (cid:105) , γ S , i δ (cid:101) v || (cid:16) J ( p + ) S , || (cid:17) i + (cid:101) v − ⊥ δ (cid:101) v ⊥ (cid:16)(cid:101) v ⊥ J ( p + ) S , ⊥ (cid:17) i (cid:43) δ (cid:101) v = . (D.17)Here, a (cid:13) and b (cid:13) are the discrete fluxes for inertial term in Eq. (5.4) [details shown in Eqs. (E.16) and (E.17)],and γ S is trivially evaluated as: γ S , i = − (cid:42) (cid:104) v ( p ) i (cid:105) , δ (cid:101) v ⊥ (cid:16) J ( p + ) S , ⊥ (cid:17) i (cid:43) δ (cid:101) v (cid:42) (cid:104) v ( p ) i (cid:105) , δ (cid:101) v || (cid:16) J ( p + ) S , || (cid:17) i (cid:43) δ (cid:101) v . (D.18)51 ppendix E. Discretization of the Vlasov Equation Phase-Space Fluxes The individual terms in Eq. (5.4) are elaborated. The term a (cid:13) corresponds to the discrete representationof the spatial streaming term, with F ( p + ) r , α , i + / , j , k = J S , i + / v ∗ , ( p ) α , i + / (cid:16)(cid:101) v || , j + (cid:98) u ∗ , ( p ) || , α , i + / (cid:17) SMART (cid:16)(cid:101) v || , j + (cid:98) u ∗ , ( p ) || , α , i + / , (cid:101) f ( p + ) α (cid:17) i + / , j , k , (E.1) v ∗ , ( p ) α , i + / = v ∗ , ( p ) α , i + + v ∗ , ( p ) α , i , and (cid:98) u ∗ , ( p ) || , α , i + / = (cid:98) u ∗ , ( p ) || , α , i + (cid:98) u ∗ , ( p ) || , α , i + , where SMART ( a , φ ) f ace denotes a SMART cell-face interpolation operation [62] of a scalar φ at a cell facewith a given velocity a , SMART ( a , φ ) f ace = N ∑ i (cid:48) = c f ace , i (cid:48) ( a , φ ) φ i (cid:48) . Here, coefficients c f ace , i (cid:48) are the interpolation weights for cell index i (cid:48) surrounding the cell face.The term b (cid:13) corresponds to the inertial term in the configuration space, arising from the moving grid,with F ( p + ) ˙ r , α , i + / , j , k = − ˙ J ( p + ) r , i + / SMART (cid:16) − ˙ J ( p + ) r , i + / , (cid:101) f ( p + ) α (cid:17) i + / , j , k , (E.2)where the definition of ˙ J r , i + / is given in Sec. 5.3.The term c (cid:13) corresponds to the electrostatic-acceleration term with J ( p + ) acc , α , i , j + / , k = J ( p + ) S ξ , i q α m α E ( p + ) || , i v ∗ , p α , i SMART (cid:16) q α E ( p + ) || , i , (cid:101) f ( p + ) α (cid:17) i , j + / , k . (E.3)The term d (cid:13) corresponds to the inertial terms due to temporal variation of the velocity-space metrics(i.e., v ∗ α and (cid:98) u ∗|| , α ) with J ( p + ) t , || , α , i , j + / , k = I ( p ) || , t , α , i , j + / , k SMART (cid:16) I ( p + ) || , t , α , i + / , j + / , k , (cid:101) f ( p + ) α (cid:17) i , j + / , k (E.4)and J ( p + ) t , ⊥ , α , i , j , k + / = I ( p + ) ⊥ , t , α , i , j , k + / SMART (cid:16) I ( p + ) ⊥ , t , α , i , j , k + / , (cid:101) f ( p + ) α (cid:17) i , j , k + / , (E.5) I ( p + ) t , || , α , i , j + / , k = − J ( p + ) S ξ , i v ∗ , ( p ) α , i δ t (cid:16) v ∗ α (cid:104)(cid:101) v || , j + / + (cid:98) u ∗|| , α (cid:105)(cid:17) ( p + ) i , (E.6)and I ( p + ) t , ⊥ , α , i , j , k + / = − J ( p + ) S ξ , i v ∗ , ( p ) α , i (cid:104) δ t ( v ∗ α (cid:101) v ⊥ ) ( p + ) i , k + / (cid:105) . (E.7)As also described in Ref. [44], we lag the time level between the BDF2 coefficients and the normalizationspeed (and similarly, the shift velocity) to avoid over-constraining the nonlinear residual (we refer the readersto the reference for further detail). 52he term e (cid:13) corresponds to the inertial terms due to the spatial variation of the metrics with J ( p + ) , − r , || , α , i , j + / , k = I ( p + ) , − r , || , α , i , j + / , k SMART (cid:16) I ( p + ) , − r , || , α , i , j + / , k , (cid:101) f ( p + ) α (cid:17) i , j + / , k (E.8) J ( p + ) , + r , || , α , i , j + / , k = I ( p + ) , + r , || , α , i , j + / , k SMART (cid:16) I ( p + ) , + r , || , α , i , j + / , k , (cid:101) f ( p + ) α (cid:17) i , j + / , k (E.9)and J ( p + ) , − r , ⊥ , α , i , j , k + / = I ( p + ) , − r , ⊥ , α , i , j , k + / SMART (cid:16) I ( p + ) , − r , ⊥ , α , i , j , k + / , (cid:101) f ( p + ) α (cid:17) i , j , k + / , (E.10) J ( p + ) , + r , ⊥ , α , i , j , k + / = I ( p + ) , + r , ⊥ , α , i , j , k + / SMART (cid:16) I ( p + ) , + r , ⊥ , α , i , j , k + / , (cid:101) f ( p + ) α (cid:17) i , j , k + / , (E.11)where I ( p + ) , − r , || , α , i , j + / , k = − (cid:26)(cid:104) J S v ∗ α (cid:16)(cid:101) v || + (cid:98) u ∗|| , α (cid:17) − ˙ J r (cid:105) ∂∂ ξ (cid:104) v ∗ α (cid:16)(cid:101) v || + (cid:98) u ∗|| , α (cid:17)(cid:105)(cid:27) ( p + ) i + / , j + / , k ≈ − (cid:104) J S , i + / v ∗ , ( p ) α , i + / (cid:16)(cid:101) v || , j + / + (cid:98) u ∗ , ( p ) || , α , i + / (cid:17) − ˙ J ( p + ) r , i + / (cid:105) × v ∗ , ( p ) α , i + (cid:16)(cid:101) v || , j + / + (cid:98) u ∗ , ( p ) || , α , i + (cid:17) − v ∗ , ( p ) α , i (cid:16)(cid:101) v || , j + / + (cid:98) u ∗ , ( p ) || , α , i (cid:17) ∆ ξ , (E.12) I ( p + ) , + r , || , α , i , j + / , k = − (cid:26)(cid:104) J S v ∗ α (cid:16)(cid:101) v || + (cid:98) u ∗|| , α (cid:17) − ˙ J r (cid:105) ∂∂ ξ (cid:104) v ∗ α (cid:16)(cid:101) v || + (cid:98) u ∗|| , α (cid:17)(cid:105)(cid:27) ( p + ) i − / , j + / , k ≈ − (cid:104) J S , i − / v ∗ , ( p ) α , i + / (cid:16)(cid:101) v || , j + / + (cid:98) u ∗ , ( p ) || , α , i − / (cid:17) − ˙ J ( p + ) r , i − / (cid:105) × v ∗ , ( p ) α , i (cid:16)(cid:101) v || , j + / + (cid:98) u ∗ , ( p ) || , α , i (cid:17) − v ∗ , ( p ) α , i − (cid:16)(cid:101) v || , j + / + (cid:98) u ∗ , ( p ) || , α , i − (cid:17) ∆ ξ , (E.13) I ( p ) , − r , ⊥ , α , i , j , k + / = − (cid:26) J S (cid:104) v ∗ α (cid:16)(cid:101) v || + (cid:98) u ∗|| , α (cid:17) − ˙ J r (cid:105) ∂∂ ξ ( v ∗ α (cid:101) v ⊥ ) (cid:27) ( p + ) i + / , j , k + / ≈ − (cid:104) J S , i + / v ∗ , ( p ) α , i + / (cid:16)(cid:101) v || , j + (cid:98) u ∗ , ( p ) || , α , i + / (cid:17) − ˙ J ( p + ) r , i + / (cid:105) v ∗ , ( p ) α , i + (cid:101) v ⊥ , k + / − v ∗ , ( p ) α , i (cid:101) v ⊥ , k + / ∆ ξ , (E.14)and I ( p ) , + r , ⊥ , α , i , j , k + / = − (cid:26) J S (cid:104) v ∗ α (cid:16)(cid:101) v || + (cid:98) u ∗|| , α (cid:17) − ˙ J r (cid:105) ∂∂ ξ ( v ∗ α (cid:101) v ⊥ ) (cid:27) ( p + ) i − / , j , k + / ≈ − (cid:104) J S , i − / v ∗ , ( p ) α , i − / (cid:16)(cid:101) v || , j + (cid:98) u ∗ , ( p ) || , α , i − / (cid:17) − ˙ J ( p + ) r , i − / (cid:105) v ∗ , ( p ) α , i (cid:101) v ⊥ , k + / − v ∗ , ( p ) α , i − (cid:101) v ⊥ , k + / ∆ ξ . (E.15)The term f (cid:13) corresponds to the velocity-space inertial term due to the transformation of VFP equationfrom Cartesian to spherical geometry with J ( p + ) S , || , α , i , j + / , k = r ( p + ) i (cid:104) v ∗ , ( p ) α , i (cid:101) v ⊥ , k (cid:105) SMART (cid:16)(cid:101) v ⊥ , k , (cid:101) f ( p + ) α (cid:17) i , j + / , k (E.16)53nd J ( p + ) S , ⊥ , α , i , j , k + / = r ( p + ) i (cid:104) v ∗ , ( p ) α , i (cid:105) (cid:101) v ⊥ , k + / (cid:16)(cid:101) v || , j + (cid:98) u ∗ , ( p ) || , i (cid:17) SMART (cid:16)(cid:101) v ⊥ , k + / (cid:16)(cid:101) v || , j + (cid:98) u ∗ , ( p ) || , i (cid:17) , (cid:101) f ( p + ) α (cid:17) i , j , k + / . (E.17) Appendix F. Velocity-Space Grid Adaptivity
For completeness, we briefly review the velocity-space grid adaptivity strategy from Ref. [50]. Thevelocity-space grid metrics, are updated as: v ∗ , ( p + ) α = v ∗ , ( p ) α + ∆ t ( p ) ˙ v ∗ , ( p + ) α if ∆ t ( p ) (cid:12)(cid:12)(cid:12) ˙ v ∗ , ( p + ) α (cid:12)(cid:12)(cid:12) v ∗ , ( p ) α ≤ . v ∗ , ( p ) α (cid:104) + . ∆ t sign (cid:16) ˙ v ∗ , ( p + ) α (cid:17)(cid:105) otherwise , (F.1)and u ∗ , ( p + ) || , α = u ∗ , ( p ) || , α + ∆ t ( p ) ˙ u ∗ , ( p + ) || , α if ∆ t ( p ) (cid:12)(cid:12)(cid:12)(cid:12) ˙ u ∗ , ( p + ) || , α v ∗ , ( p ) α (cid:12)(cid:12)(cid:12)(cid:12) ≤ . u ∗ , ( p ) || , α (cid:104) + . ∆ t sign (cid:16) ˙ u ∗ , ( p + ) || , α (cid:17)(cid:105) otherwise , (F.2)where ˙ v ∗ , ( p + ) α = (cid:104) v ∗ , ( p + ) α (cid:105) † − v ∗ , ( p ) α ∆ t ( p ) , (F.3)˙ u ∗ , ( p + ) || , α = (cid:104) u ∗ , ( p + ) || , α (cid:105) † − u ∗ , ( p ) || , α ∆ t ( p ) , (F.4) (cid:104) v ∗ , ( p + ) α (cid:105) † = (cid:115) T ( p + ) α m α , (F.5) (cid:104) u ∗ , ( p + ) || , α (cid:105) † = u ( p + ) || , α + ∆ w ( p + ) || , α , (F.6)and ∆ w ( p + ) || , α = (cid:28)(cid:16) v ( p ) || − u ( p + ) || , α (cid:17) (cid:16) v ( p ) − u ( p + ) α (cid:17) , f ( p + ) α (cid:29) v (cid:28)(cid:16) v ( p ) − u ( p + ) α (cid:17) , f ( p + ) α (cid:29) v . (F.7)To ensure that the profiles of v ∗ and u ∗|| are smooth in space, we employ a Winslow smoother. The gridadaptivity algorithm is summarized in Alg. 4. 54 lgorithm 4 Update procedure for the velocity-space grid metrics.1. Compute the velocity-space grid metrics, v ∗ and u ∗|| from Eqs. (F.1)-(F.7).2. Winslow smooth the v ∗ and u ∗|| by solving for (cid:104) − λ ω ∂ ξ ξ (cid:105) v ∗ = v ∗∗ and (cid:104) − λ ω ∂ ξ ξ (cid:105) u ∗|| = u ∗∗|| where λ ω = − in this study and the superscript, ∗∗ , denotes the non-smoothed quantities. Appendix G. Poisson Grid Generation and Initial Optimization
It is critical to ensure an initial high quality grid. Severe numerical pollution of the solution can oc-cur if one fails at this task and solely relies on MMPDE to transiently resolve the initial gradients as isdemonstrated in Sec. (6.2). To address this issue, we employ an initial Poisson mesh generation and gridoptimization to resolve initial gradients.Consider a reference solution, F ∗ ( r ∗ ) ∈ R N ∗ ξ and N ∗ ξ (cid:29) N ξ , defined on a highly resolved reference grid, r ∗ ∈ R N ∗ ξ . We begin by interpolating this solution onto the coarse, uniform computational grid, r ( ξ ) ∈ R N ξ , such that F = M F ∗ , where F ∈ R N ξ is the reference solution interpolated onto the initial coarse-computational grid and M = M (cid:0) r , F ∗ (cid:1) ∈ R N ξ × N ∗ ξ is the interpolation map. Note that, to ensure themonotonicity of the interpolated quantity, we restrict ourselves to linear interpolation in this study. The gridoptimization strategy is shown in Algorithm 5. In Fig. G.1, we illustrate the iterative procedure for the grid Algorithm 5
Grid optimization procedure.1. Initialize iteration index, k = k = k + ω (cid:0) F k − (cid:1) .4. Solve the Poisson equation, ∂∂ξ (cid:104) ω (cid:0) F k − (cid:1) ∂ r † ∂ξ (cid:105) = , for r † .5. Update the grid location: r k = r k − + w relax (cid:0) r † − r k − (cid:1) . Here w relax = .
75 is the under-relaxation factorused in this study.6. Interpolate reference solution onto the new grid: F k = M k F ∗ .
7. If (cid:114) N ξ + ∑ N ξ i = (cid:16) r ki − r k − i r k − i + − (cid:17) ≤ − exit iteration, else return to 2.optimization. Appendix H. Guderley and Van-Dyke Problem
The Guderley solution is the self-similar solution to the Euler equations for a strong, spherically con-verging/diverging shock. In 1D spherical geometry, we may write the Euler equations as [91]: ∂ t ρ + r ∂ r (cid:0) r ρ u r (cid:1) = , (H.1) ( ∂ t + u r ∂ r ) u r = − ρ ∂ r P , (H.2) ( ∂ t + u r ∂ r ) ln P ρ γ = , (H.3)where P is the gas pressure, u r is the radial drift velocity, and ρ is the mass density. Guderley was the firstto show that Eqs. (H.1-H.3) admits a self-similar solution [70, 72, 91, 71, 73] with the shock position: R s ( t ) ∝ | t | α , (H.4)55 [ A . U .] Initial grid F ∗ F r [A.U.]0 0.5 1 ω [ A . U .] Early iteration r [A.U.]0 0.5 100.511.5 0 0.5 100.51
Intermediate iteration r [A.U.]0 0.5 100.511.5 0 0.5 100.51
Converged grid r [A.U.]0 0.5 100.511.5
Figure G.1: Illustration of the iterative procedure for the initial grid generation. The top row shows the reference solution, F ∗ (blue)and the interpolated solution, F (red), while the bottom row shows the corresponding monitor function. The black markers denotethe grid points. As can be seen, the grid adapts to resolve the initial gradients. where α is a constant which differs between the converging and diverging phases. This constant must beobtained by numerical integration of the Euler equations, with the coordinate transformation: r → ξ ≡ r / R s ( t ) . In this self-similar variable, ξ , the fluid quantities follow as [91]: u r = rt V ( ξ ) , (H.5) ρ = ρ G ( ξ ) , (H.6) c s = (cid:16) rt (cid:17) Z ( ξ ) , (H.7)where c s ≡ γ P / ρ is the sound speed, γ is the adiabatic index (which is 5 / ρ is theinitial (upstream) density, and the functions V, G, and Z are determined via numerical integration.Whereas the Guderley solution is self-similar (i.e., it assumes a converging shock that has “forgotten”its initial conditions), Van Dyke & Guttmann [84] consider the creation of a shock from the surface of aspherically converging piston (moving at a constant velocity). The key assumption is that the shock, afteremerging from the piston face, is approximately planar. The Euler equations are then solved with a Taylorseries expansion of the shock position: R s ( t ) = ∞ ∑ n = R n t n , (H.8)where each subsequent term represents a spherical correction to the planar trajectory (which is given by theRankine-Hugoniot conditions). Likewise, Taylor series expansions were employed for the pressure, density,and drift velocities: u r = ∞ ∑ n = U n ( ξ ) t n − , (H.9) ρ r = ∞ ∑ n = R n ( ξ ) t n − , (H.10) P = ∞ ∑ n = P n ( ξ ) t n − , (H.11)where the similarity variable is ξ ≡ γ − (cid:16) ru p t − (cid:17) , and u pp