An adaptive, implicit, conservative 1D-2V multi-species Vlasov-Fokker-Planck multiscale solver in planar geometry

# An adaptive, implicit, conservative 1D-2V multi-species Vlasov-Fokker-Planck multiscale solver in planar geometry

W. T. Taitano L. Chacón A. N. Simakov Theoretical Division Los Alamos National Laboratory, Los Alamos, NM 87545 Theoretical Design Division, Los Alamos National Laboratory, Los Alamos, NM 87545
###### Abstract

We consider a 1D-2V Vlasov-Fokker-Planck multi-species ionic description coupled to fluid electrons. We address temporal stiffness with implicit time stepping, suitably preconditioned. To address temperature disparity in time and space, we extend the conservative adaptive velocity-space discretization scheme proposed in [Taitano et al., J. Comput. Phys., 318, 391–420, (2016)] to a spatially inhomogeneous system. In this approach, we normalize the velocity-space coordinate to a temporally and spatially varying local characteristic speed per species. We explicitly consider the resulting inertial terms in the Vlasov equation, and derive a discrete formulation that conserves mass, momentum, and energy up to a prescribed nonlinear tolerance upon convergence. Our conservation strategy employs nonlinear constraints to enforce these properties discretely for both the Vlasov operator and the Fokker-Planck collision operator. Numerical examples of varying degrees of complexity, including shock-wave propagation, demonstrate the favorable efficiency and accuracy properties of the scheme.

###### keywords:
Conservative discretization, thermal velocity based adaptive grid, 1D2V, Fokker-Planck, Rosenbluth potentials

## 1 Introduction

The Vlasov-Fokker-Planck (VFP) collisional kinetic description, coupled with Maxwell’s equations, is regarded as a first-principles physical model for describing weakly coupled plasmas in all collisionality regimes, and accordingly, has a wide range of applications in laboratory (e.g., magnetic and inertial thermonuclear fusion), space (e.g., Earth’s magnetosphere), and astrophysical (e.g., stellar mass ejections) plasmas. In the VFP system, collisions are modeled by the Fokker-Planck collision operator, which describes collisional relaxation of particle distribution functions in plasmas under the assumption of binary, grazing-angle collisions (rosenbluth, ; arsenev_1991conn_boltz_lfpe, ; desvillettes_1992_asymp_boltz_eqn_graze_col, ; degond_1992_fp_asymp_bolt_col_op_coul, ; goudon_1997_jsp, ; landau_1937_fokker_planck, ). Mathematically, the Fokker-Planck operator is integro-differential, non-local, and very difficult to invert.

The system of VFP equations for various plasma species supports disparate length and time scales, as well as arbitrary temperature disparity in time and space, which makes this system particularly challenging to solve with grid-based approaches. The challenges of temperature disparity are evident when one considers the thermal speed, , which provides a characteristic width of the plasma species distribution function and is a function of the plasma temperature, , and particle species mass, . In many practical applications of interest, variation for a given species can span several orders of magnitude in configuration space. In addition, mass differences result in strong disparities for different species. Since the velocity-space domain size is determined for a given species by the hottest region (large ), and the velocity-grid spacing must resolve the coldest region (small ), velocity-space discretizations with uniform Cartesian grids in such scenarios may lead to impractical grid size requirements.

Several studies recognized and tried to address these challenges by normalizing the velocity coordinate to the local thermal velocity (larroche_2007_lsse_jcp, ; larroche_EPJ_2003_icf_fuel_ion_implosion_sim, ; jarema_CPC_2015_block_structured_grid_adaptive_vpace, ; peigney_JCP_2014_fp_kinetic_modeling_of_alpha, ). In this fashion, the grid will expand as the plasma heats, and contract as it cools. Particularly relevant to this study is the work in Ref. (larroche_EPJ_2003_icf_fuel_ion_implosion_sim, ), where the velocity-space domain was adapted for multiple ion species based on a single local average (over the ion species) and hydrodynamic velocity of the plasma. This powerful strategy enabled the fully kinetic implosion simulations of inertial confinement fusion (ICF) capsules (larroche_EPJ_2003_icf_fuel_ion_implosion_sim, ; larroche_pop_2012_Dhe_3_icf_sim, ; inglebert_epl_2014_species_separation_kinetic_effect_neutron_diagnostics, ), but required intermittent remapping in both the physical and velocity space. None of these strategies conserve mass, momentum, and energy, and some of them (jarema_CPC_2015_block_structured_grid_adaptive_vpace, ) break the structured nature of the underlying computational mesh.

Recently, a novel strategy that deals with strong temperature disparity, avoids remapping, and works on structured meshes was proposed in Ref. (Taitano_2016_rfp_0d2v_implicit, ) for the 0D-2V multispecies Fokker-Planck equation. The strategy employs a multiple-grid approach by normalizing each species’ velocity to its thermal speed. The Fokker-Planck equation was transformed analytically, and then discretized on a mesh. The transformed equations exposed the continuum conservation symmetries, which were then enforced in the discrete via nonlinear constraints. This strategy ensures that the species’ distribution function is always well resolved regardless of temperature or mass disparity.

In this study, we extend the conservative, multiple-dynamic velocity-space adaptive strategy developed in Ref. (Taitano_2016_rfp_0d2v_implicit, ) to a spatially inhomogeneous, 1D Cartesian system. We consider a quasi-neutral plasma with multiple kinetic ion species and fluid electrons. As before, ionic species are evolved on a velocity-space grid normalized to a temporally and spatially varying characteristic speed, (a function of their ), and the VFP equation is analytically transformed. This transformation introduces additional inertial terms, which are carefully discretized to ensure simultaneous conservation of mass, momentum, and energy.

The rest of the paper is organized as follows. Section 2 introduces the VFP and fluid-electron equations and discusses their conservation properties. In Sec. 3, we introduce the normalized Vlasov equation, and provide a detailed discussion on the implementation of the proposed schemes in the following order: 1) a discretization of the Vlasov-Fokker-Planck equation with the additional inertial terms, 2) a discretization of the fluid electron equation, 3) our discrete conservation strategy for fluid electrons and kinetic ions, 4) a discrete conservation strategy for the Vlasov component with the added inertial terms, and 5) temporal and spatial evolution of . The numerical performance of the scheme is demonstrated with various multi-species tests of varying degrees of complexity in Sec. 4. Finally, we conclude in Sec. 5.

## 2 The multi-species Vlasov-Rosenbluth-Fokker-Planck equation with fluid electrons

A dynamic evolution of weakly-coupled collisional plasmas is described by the Vlasov-Fokker-Planck equation for the particle distribution function (PDF), in configuration space, , velocity space, , and time, :

 ∂tfα+∇x⋅(vfα)+qαmα∇v⋅[(E+v×B)fα]=Ns∑β=1Cαβ, (2.1)

where is the electric field, is the magnetic field, is the total number of plasma species in the system, and is the Fokker-Planck collision operator for species colliding with species :

 Cαβ=Γαβ∇v⋅[\lx@stackrel↔\mathsfslDβ⋅∇vfα−mαmβAβfα]. (2.2)

Here, and are the tensor-diffusion and friction coefficients for species , and are the masses of species and , respectively, is the ionization state of species , is the proton charge, and is the Coulomb logarithm ( is assumed for simplicity in this study for all species).

The Rosenbluth formulation of the Fokker-Planck collision operator (rosenbluth, ) computes the velocity-space-transport coefficients from the so-called Rosenbluth potentials , as:

 \lx@stackrel↔\mathsfslDβ = ∇v∇vGβ, (2.3) Aβ = ∇vHβ, (2.4)

which, in turn, are computed from the distribution function of species as:

 ∇2vHβ=−8πfβ, (2.5)
 ∇2vGβ=Hβ. (2.6)

The Rosenbluth form is completely equivalent to the integral Landau form (landau_1937_fokker_planck, ), but more advantageous algorithmically because it can be inverted with complexity (with the number of degrees of freedom in velocity space) (Taitano_2015_rfp_0d2v_implicit, ).

The collision operator, Eq. (2.2), preserves the positivity of , and conserves mass, momentum, and energy. The conservation properties stem from the following symmetries (braginskii, ):

 ⟨1,Cαβ⟩v = 0, (2.7) mα⟨v,Cαβ⟩v = −mβ⟨v,Cβα⟩v, (2.8) mα⟨v22,Cαβ⟩v = −mβ⟨v22,Cβα⟩v, (2.9)

where the inner product is defined as (for the cylindrically symmetric coordinate system in the velocity space employed herein). These conservation symmetries can be enforced in the discrete following the general procedures discussed in Refs. (Taitano_2015_rfp_0d2v_implicit, ; Taitano_2016_rfp_0d2v_implicit, ).

In this study, we consider a 1D planar geometry in the configuration space without a magnetic field. Without loss of generality, we consider a 2V cylindrically symmetric coordinate system in the velocity space. We adopt a fluid-electron model with a reduced ion-electron collision operator. We obtain the following simplified system of equations comprised of the ion-Vlasov-Fokker-Planck equation (per species ),

 ∂tfα+∂x(v||fα)+qαmαE||∂v||fα=Ns∑βCαβ+Cαe (2.10)

and the electron-temperature equation,

 32∂∂t[neTe]+52∂x[u||,eneTe]+∂xQ||,e−qeneu||,eE||=Ns∑αWeα. (2.11)

Here, is the parallel electron-heat flux, is the electron density, is the electron temperature, is the parallel electron fluid velocity, is the electron pressure, describes the electron-ion energy exchange,

 Weα=F||,αeu||,α+3νeαmemαne(Tα−Te), (2.12)

is the friction force between the -ion species and electrons, and

 νeα=4√2πnαq2αe4Λeα3√meT3/2e (2.13)

is the electron-ion collision frequency. The frictional force between the -ion species and electrons is given by,

 (2.14)

and the electron heat flux by

 Qe=β0neTe(ue−⟨uα⟩)−κe∇Te. (2.15)

Definitions of coefficients , , and the collision-frequency-averaged ion velocity, , can be found in App. A and in Ref. (simakov_PoP_2014_e_transp_wh_multi_ion, ) .

The parallel (to the x-axis) component of the ambipolar-electric field, , in Eq. (2.11) is found from the inertialess electron momentum equation ,

 E||=∑NsαF||,αe+∂xPeqene. (2.16)

Finally, we assume quasi-neutrality,

 ne=Ns∑αZαnα, (2.17)

and ambipolarity,

 u||,e=∑NsαZαnu||,αne, (2.18)

to close the system.

The electron-ion collision operator in Eq. (2.10) is given by:

 Cαe=Γαe∇v⋅[\lx@stackrel↔\mathsfslDαe⋅∇vfα−mαmeAαefα], (2.19)

where we adopt the reduced ion-electron potentials (hazeltine_1991_plas_confin, ):

 (2.20)
 (2.21)

Here, , , is the electron thermal speed, and the transport coefficients are:

 \lx@stackrel↔\mathsfslDαe=∇v∇vGαe=43ne√π\lx@stackrel↔\mathsfslIvth,e, (2.22)
 Aαe=Aαe,u+Aαe,F, (2.23)

 Aαe,u=∇vHαe,u=−83ne√πv3th,e(v−uα) (2.24)

and

 Aαe,F=∇vHαe,F=83ne√πv3th,eFαemeneνeα. (2.25)

### 2.1 Conservation properties of the kinetic-ion/fluid-electron system

The coupled kinetic-ion and fluid-electron system possesses continuum conservation properties. We remark that, although these properties are well known and have been discussed by others in the past (hazeltine, ), we reproduce them here to explicitly expose the continuum symmetries that are required to ensure the conservation properties. The goal is to develop a strategy that can ensure these properties discretely (Sec. 3.4).

Mass conservation follows trivially from the ion mass conservation and quasi-neutrality. Momentum conservation follows from taking the moment of the ion Vlasov equation summed over all ions:

 Ns∑α[mα∂tI||,α+mα∂xS2,||,α−qαnαE||]=Ns∑αmα⎡⎣Ns∑β⟨v||,C(fα,fβ)⟩v+⟨v||,Cαe⟩v⎤⎦, (2.26)

where is the parallel specific momentum density flux, is the parallel-parallel component of the pressure tensor, and is the parallel electrostatic force. The first term on the right-hand side (RHS) of Eq. (2.26) vanishes due to momentum conservation across ion species (Taitano_2015_rfp_0d2v_implicit, ). The second term on the RHS uses the reduced-collision operator between ion species and electrons, and it can be expanded as follows:

 (2.27)

Here, the terms and vanish independently:

 mα⟨v||,Γαe∇v⋅[\lx@stackrel↔\mathsfslDαe⋅∇vfα]⟩v=mαΓαe43nevth,e√π⟨v||,∇v⋅[\lx@stackrel↔\mathsfslI⋅∇vfα]⟩v=0 (2.28)

and

 mαΓαe⟨v||,∇v⋅[mαmeAαe,ufα]⟩v=−mαΓαemαme83nev3th,e√π⟨v||,∇v⋅[(v−uα)fα]⟩v=0. (2.29)

The term becomes [by Eq. (2.25)]:

 mαΓαe⟨v||,∇v⋅[mαmeAαe,Ffα]⟩v=−F||,αe. (2.30)

When combined with quasi-neutrality, Eq. (2.17), and the ion electrostatic acceleration term:

 Ns∑αqαnαE||=−(∂xPe+Ns∑αF||,αe), (2.31)

Eq. 2.26 yields the total plasma momentum equation

 Ns∑αmα[∂tIα+∂xS2,||,α]+∂xPe=0, (2.32)

which is in a conservative form.

To show energy conservation, we take the second moment () of the ion Vlasov equation as:

 Ns∑α[mα∂tUα+mα∂xS3,||,α−qαnαu||,αE||]=Ns∑αmα⎡⎣Ns∑β⟨v22,C(fα,fβ)⟩v+⟨v22,Cαe⟩v⎤⎦. (2.33)

Here, is the specific-energy density, is the parallel component of the specific-energy density flux, and . The first term on the RHS vanishes owing to energy conservation across ion species. The second term on the RHS can be expanded as follows:

 mα⟨v22,Cαe⟩v=mαΓαe⎡⎢ ⎢ ⎢ ⎢ ⎢⎣⟨v22,∇v⋅[\lx@stackrel↔\mathsfslDαe⋅∇fα]⟩v\textcircleda−mαme⟨v22,∇v⋅⎡⎢ ⎢⎣⎛⎜ ⎜⎝Aαe,u\textcircledb+Aαe,F\textcircledc⎞⎟ ⎟⎠fα⎤⎥ ⎥⎦⟩v⎤⎥ ⎥ ⎥ ⎥ ⎥⎦. (2.34)

Here, terms , , and independently yield [using Eq. (2.22) and (2.23)]:

 Γαemα⟨v22,∇v⋅[\lx@stackrel↔\mathsfslDαe⋅∇vfα]⟩v=mαΓαe43nevth,e√π⟨v22,∇v⋅[\lx@stackrel↔\mathsfslI⋅∇vfα]⟩v=3νeαmemαneTe, (2.35)
 mαΓαe⟨v22,∇v⋅[mαmeAαe,ufα]⟩v=3νeαmemαneTα, (2.36)

and

 mαΓαe⟨v22,∇v⋅[mαmeAαe,Ffα]⟩v=−u||,αF||,αe. (2.37)

Gathering terms, the energy moment of the ion-electron collision operator yields:

 mα⟨v22,Cαe⟩v=3νeαmemαne(Te−Tα)+u||,αF||,αe=−Weα. (2.38)

Using the fluid electron temperature equation and ambipolarity, we finally obtain:

 ∂t(32neTe+Ns∑αmαUα)+∂x(52u||,eneTe+Ns∑αmαS3,α,||+Q||,e)=0, (2.39)

which is a conservative form of the evolution equation for the total plasma energy density .

## 3 Numerical implementation

### 3.1 VFP equation in normalized velocity variables

We consider the normalization of all velocity-space quantities for a species to a reference speed, , related to their thermal speed, as follows:

 ˆvα=vv∗α,∂ˆv||=v∗α∂ˆv||,ˆfα=(v∗α)3fα.

Here, the hat denotes quantities normalized to . As an example, the density, drift, and temperature moments are defined as:

 nα=⟨1,ˆfα⟩ˆv,u||,α=v∗αˆu||,α=v∗α⟨ˆv||,ˆfα⟩ˆv,Tα=(v∗α)2ˆTα=mα(v∗α)2⟨(ˆv−ˆuα)2,ˆfα⟩ˆv3⟨1,ˆfα⟩ˆv

where .The normalization of other relevant quantities and the collision operator are discussed in Ref. (Taitano_2016_rfp_0d2v_implicit, ). We note that, as is a function of local for a given plasma species (elaborated in Sec. 3.6), the grid will expand as the plasma heats, and contract as it cools; refer to Fig. 3.1.

The temporal and spatial dependence of introduces inertial terms in the Vlasov equation, Eq. (2.10), which after the transformation reads (App. B):

 ∂tˆfα+∂x(v∗αˆv||ˆfα) + ∂ˆv||[(qαmαE||v∗α)ˆfα]−ˆ∇ˆv⋅[(∂tv∗αv∗α)ˆvˆfα]−ˆ∇ˆv⋅[(∂xv∗α)ˆvˆv||ˆfα] (3.1) = (v∗α)3⎛⎝Ns∑βCαβ+Cαe⎞⎠=Ns∑βˆCαβ+ˆCαe.

### 3.2 Discretization of the VFP equation with inertial terms

We discretize the VFP equation using finite volumes in a 1D planar configuration space () and 2V cylindrical-velocity space ( and ) with azimuthal symmetry. We compute the discrete volume for cell (i,j,k) as:

 ΔVi,j,k=2πΔxΔˆv||ˆv⊥,kΔˆv⊥,

where , , and are mesh spacings in the configuration space and the parallel- and perpendicular-velocity space, respectively. For a uniform mesh (assumed henceforth), we have:

 Δx=LxNx,Δˆv||=L||N||,Δˆv⊥=L⊥N⊥,

where , , and are the configuration space, and parallel and perpendicular velocity-space domain sizes, respectively, and , , and are the corresponding number of cells. The mesh is arranged such that cell faces map to the domain boundary (and therefore outermost cell centers are half a mesh-spacing away from the boundary). We define the distribution function and the Rosenbluth potentials , at cell centers.

Velocity-space inner products are approximated via a mid-point quadrature rule as

 ⟨A,B⟩ˆv≈2πN||∑j=1N⊥∑k=1ˆv⊥,kΔˆv||Δˆv⊥Aj,kBj,k (3.2)

for scalars and

 ⟨A,B⟩ˆv≈2π⎡⎢⎣N||∑j=0N⊥∑k=1ˆv⊥,kΔˆv||Δˆv⊥A||,j+1/2,kB||,j+1/2,k+N||∑j=1N⊥∑k=0ˆv⊥,k+1/2Δˆv||Δˆv⊥A⊥,j,k+1/2B⊥,j,k+1/2⎤⎥⎦

for vectors (with components defined at cell faces as denoted by the half-integer indices , ).

We discretize Eq. (3.1) in a conservative form as

 cp+1ˆfp+1α,i,j,k+cpˆfpα,i,j,k+cp−1ˆfp−1α,i,j,kΔt+Fp+1α,i+1/2,j,k−Fp+1α,i−1/2,j,kΔx\textcircleda+Jp+1α,acc,i,j+1/2,k−Jp+1α,acc,i,j−1/2,kΔˆv||\textcircledb +⎡⎢⎣Jp+1||,α,t,i,j+1/2,k−Jp+1||,α,t,i,j−1/2,kΔˆv||+ˆv⊥,k+1/2Jp+1⊥,α,t,i,j,k+1/2−ˆv⊥,k−1/2Jp+1⊥,α,t,i,j,k−1/2ˆv⊥,kΔˆv⊥⎤⎥⎦\textcircledc +⎡⎢⎣Jp+1||,α,x,i,j+1/2,k−Jp+1||,α,x,i,j−1/2,kΔˆv||+ˆv⊥,k+1/2Jp+1⊥,α,x,i,j,k+1/2−ˆv⊥,k−1/2Jp+1⊥,α,x,i,j,k−1/2ˆv⊥,kΔˆv⊥⎤⎥⎦\textcircledd =Ns∑βCp+1αβ∣∣i,j,k+Cp+1αe∣∣i,j,k. (3.3)

Here, , , and are the coefficients for the second-order backwards difference formula (BDF2) (Taitano_2015_rfp_0d2v_implicit, ) and is the discrete time index.

The term corresponds to the discretization of the spatial streaming term, with

 Fp+1α,i+1/2,j,k=v∗,pα,i+1/2ˆv||,j%Interp(ˆv||,j,ˆfp+1α)i+1/2,j,k,
 v∗α,i+1/2=v∗α,i+1+v∗α,i2,

where is an advection interpolation operator of a scalar at a cell face based on a given velocity , which can be written in general as

 Interp(a,ϕ)face=N∑i′=1ωface,i′(a,ϕ)ϕi′. (3.4)

The coefficients are the interpolation weights for the spatial cells surrounding the cell face of interest (in this study, they are determined by the SMART discretization (smart, )).

The term corresponds to the electrostatic-acceleration term with

 Jp+1α,acc,i,j+1/2,k=Ap+1α,iInterp(Ap+1α,i,ˆfp+1α)i,j+1/2,k, (3.5)

where

 Ap+1α,i=qαmαEp+1||,iv∗,pα,i.

The term corresponds to the inertial term due to temporal variation of the normalization velocity, with

 Jp+1||,α,t,i,j+1/2,k=Ipα,t,iˆv||,j+1/2Interp(Ipα,t,iˆv||,j+1/2,ˆfp+1α)i,j+1/2,k

and

 Jp+1⊥,α,t,i,j,k+1/2=Ipα,t,iˆv⊥,k+1/2Interp(Ipα,t,iˆv⊥,k+1/2,ˆfp+1α)i,j,k+1/2,

where

 Ipα,t,i=−(∂tv∗,pα)iv∗,pα,i≈−cp+1v∗,pα,i+cpv∗,p−1α,i+cp−1v∗,p−2α,iv∗,pα,iΔt. (3.6)

We lag the time level between the BDF2 coefficients and the normalization velocity for well-posedness of the velocity-space grid motion (Taitano_2016_rfp_0d2v_implicit, ).

The term corresponds to the inertial term due to the spatial variation of the normalization velocity, , with

 Jp+1||,α,x,i,j+1/2,k=Ipα,x,iˆv2||,j+1/2Interp(Ipα,x,i,ˆfp+1α)i,j+1/2,k, (3.7)
 Jp+1⊥,α,x,i,j,k+1/2=Ipα,x,iˆv||,jˆv⊥,k+1/2Interp(Ipα,x,iˆv||,j,ˆfp+1α)i,j,k+1/2, (3.8)

where

 Ipα,x,i=−(∂xv∗,pα)i≈−v∗,pα,i+1/2−v∗,pα,i−1/2Δx. (3.9)

Finally, the right-hand-side of Eq. (3.3) corresponds to the Fokker-Planck-collision operator and its treatment is discussed in detail in Refs. (Taitano_2015_rfp_0d2v_implicit, ; Taitano_2016_rfp_0d2v_implicit, ). In this study, we use the mimetic differencing approach for the tensor diffusion operator in the collision term proposed in Ref. (lipnikov_rjnamm_2012, ).

### 3.3 Discretization of the electron temperature equation

The electron temperature equation, Eq. (2.11), is also discretized using a finite-volume scheme in space and BDF2 in time:

 32cp+1np+1e,iTp+1e,i+cpnpe,iTpe,i+cp−1np−1e,iTp−1e,iΔtk+52up+1||,e,i+1/2(˜neTe)p+1i+1/2−up+1||,e,i−1/2(˜neTe)p+1i−1/2Δx+ Qp+1||,e,i+1/2−Qp+1||,e,i−1/2Δx−qe[(neu||,e)EiE||,i]p+1=3νp+1eα,imemαnp+1e,i(Tp+1α,i−Tp+1e,i)+Fp+1||,αe,iup+1||,α,i. (3.10)

Here the tilde denotes a cell-face discretization for the advection quantities (SMART in this study). The quantity with a superscript in Eq. (3.10) are defined so as to enforce conservation properties, and will be discussed shortly. Other terms in Eq. (3.10) are defined as:

 up+1||,e,i+1/2=0.5(up+1||,e,i+1+up+1||,e,i),
 Qp+1||,e,i+1/2=[β0ne(u||,e−⟨u||,α⟩)]p+1i+1/2˜Tp+1e,i+1/2−κp+1||,e,i+1/2Tp+1e,i+1−Tp+1e,iΔx, (3.11)

where

 [β0ne(u||,e−⟨u||,α⟩)]p+1i+1/2=0.5{[β0ne(u||,e−⟨u||,α⟩)]p+1i+1+[β0ne(u||,e−⟨u||,α⟩)]p+1i}, (3.12)
 κp+1