A mass, momentum, and energy conservative dynamical low-rank scheme for the Vlasov equation
AA mass, momentum, and energy conservative dynamical low-rank schemefor the Vlasov equation
Lukas Einkemmer a, ∗ , Ilon Joseph b a Department of Mathematics, University of Innsbruck, Austria b Physics Division, Lawrence Livermore National Laboratory, United States of America
Abstract
The primary challenge in solving kinetic equations, such as the Vlasov equation, is the high-dimensionalphase space. In this context, dynamical low-rank approximations have emerged as a promising way toreduce the high computational cost imposed by such problems. However, a major disadvantage of thisapproach is that the physical structure of the underlying problem is not preserved. In this paper, we proposea dynamical low-rank algorithm that conserves mass, momentum, and energy as well as the correspondingcontinuity equations. We also show how this approach can be combined with a conservative time and spacediscretization.
Keywords: dynamical low-rank approximation, conservative numerical methods, complexity reduction,Vlasov equation, kinetic equation
1. Introduction
Solving kinetic equations efficiently is important in applications ranging from plasma physics to radiativetransfer. The main challenge in this context is the up to six-dimensional phase space and the associatedunfavorable scaling of computational cost and memory requirements, usually referred to as the curse ofdimensionality. Because of this, particle methods, such as the particle in cell (PIC) scheme, have been andare still widely used. However, it also well known that particle methods miss or do not resolve certain physicalphenomena (see, e.g., [2]) or require an immense number of particles thereby negating their advantage. Inaddition, complexity reduction techniques such as sparse grids have been investigated. While they canprovide some advantage compared to a full grid simulation the gain is usually modest [20].More recently, using dynamical low-rank approximations to solve kinetic problems has received considerableinterest. Such methods have been developed for both the Vlasov equation [19, 11, 12, 13] and radiationtransport problems [8, 27, 10]. Dynamical low-rank integrators approximate a six-dimensional Vlasov equa-tion by a set of only three dimensional advection problems. Moreover, they have a range of properties thatmakes them well suited for performing kinetic simulation. In particular, such methods are able to resolvefilamentation and are almost exact if the dynamics can be well represented by a linearized equation [13].In addition, dynamical low-rank schemes can capture the limiting fluid or diffusive regime [9, 8, 10, 28].Therefore, dynamical low-rank integrators can drastically decrease the numerically effort that is required tosolve a number of kinetic problems. This enables the simulation of such problems on desktop computers orsmall clusters that are otherwise either unfeasible or require large supercomputersA major disadvantage of dynamical low-rank integrators, however, is that they do not respect the physicalstructure of the underlying equations. In particular, mass, momentum, and energy is not preserved by the ∗ Corresponding author
Email address: [email protected] (Lukas Einkemmer)
Preprint submitted to Elsevier February 1, 2021 a r X i v : . [ m a t h . NA ] J a n ow-rank approximation. This is in stark contrast to Eulerian and semi-Lagrangian Vlasov solvers, whereat least mass and momentum is conserved [26] and methods with good long-time behavior with respect toenergy have been obtained [7], and particle methods, where usually mass and either momentum or energycan be conserved [31].Some approaches to improve this deficiency have been proposed. In [27] the solution is simply rescaled suchthat the mass is preserved. Such an approach, however, can not be extended to simultaneously conservemomentum or energy. It also does not respect the underlying continuity equation for the mass density, whichcould actually be more important for the long time behavior of the integrator [12]. In [12] a method basedon Lagrange multipliers is proposed. This approach succeeds in improving the conservative properties ofthe low-rank approach, but it does not allow us to simultaneously conserve both the continuity equationsfor density and momentum as well as the corresponding invariants. It also does not change the dynamicallow-rank approximation. Instead, this methods corrects each step of the low-rank integrator after the fact.In [28] the fluid moments, i.e. the continuity equations associated to the invariants, are integrated explicitly.These moments are then coupled to a low-rank approximation that resolves the kinetic dynamics. Since themoments are integrated explicitly this yields a conservative scheme. However, coupling the fluid momentsto the dynamical low-rank approximation without having to explicitly form the full solution, and thusencountering the curse of dimensionality, is currently only possible for linear equations. This severely limitsthe widespread applicability of this approach. In particular, it is not possible to use such a scheme for theVlasov equation considered here.Neither of the conservative or quasi-conservative methods developed in the literature solve the fundamentalproblem. Namely, that the classic dynamical low-rank approximation does not take the structure of theequations that we are trying to solve into account. In this paper we introduce such an approach thatconserves both mass, momentum, and energy and ensures that the corresponding continuity equations aresatisfied as well. The proposed method is based on the observation that if certain functions of velocitybelong to the approximation space, then the desired conservation follows. To accomplish this two steps arenecessary. First, we have to formulate the dynamical low-rank scheme in a different function space thanthe L space that is usually used (we will use an appropriately weighted L space instead). Second, wehave to constrain the low-rank factors such that the desired functions belong to the approximation spaceat all times. This is done by fixing certain basis functions and using a modified Petrov–Galerkin conditionin such a way that this is compatible with the remainder of the dynamical low-rank approximation. Thisresults in equations of motions that are somewhat different from the ones used in the classic dynamicallow-rank algorithm introduced in [18] (and that subsequently was used heavily in the literature; see, e.g,[25, 24, 11, 12, 27, 13, 10]). We also introduce a time and space discretization of the resulting equations ofmotions that conserves mass and momentum as well as a discrete version of the corresponding continuityequations up to machine precision.The remainder of the paper is structured as follows. In section 2 we introduce the Vlasov–Poisson equa-tion and discuss the invariants and corresponding physical structure of this model. This is followed by arecollection of the classic dynamical low-rank integrator in section 3. The proposed conservative dynamicallow-rank integrator is then introduced in section 4 and its properties are discussed in section 5. The fullydiscretized integrator, i.e. time and space discretization aspects, are discussed in section 6. Finally, a numberof numerical simulation are presented in section 7.
2. Vlasov–Poisson equations and conservation
In this work we consider the Vlasov–Poisson equations ∂ t f ( t, x, v ) + v · ∇ x f ( t, x, v ) − E ( f )( t, x ) · ∇ v f ( t, x, v ) = 0 , (1) ∇ x · E ( f )( t, x ) = 1 − Z f ( t, x, v ) d v, ∇ × E ( f )( t, x ) = 0 (2)2n ( x, v ) ∈ Ω x × Ω v with Ω x ⊂ R d (physical space) and Ω v ⊂ R d (velocity space), d ≤
3. The sought afterquantity is the particle-density function f . This equation models the dynamics of electrons that are subjectto a homogeneous ion background distribution. We consider the single species Vlasov–Poisson equation here,but it should be emphasized that the described algorithm can be easily extended to the multi-particle caseas well as to more complicated models, such as the Vlasov–Maxwell equations.The dynamics of the Vlasov–Poisson equations conserves a number of physically relevant invariants. Inparticular, mass M , momentum JM ( t ) = Z Ω f ( t, x, v ) d ( x, v ) , J ( t ) = Z Ω vf ( t, x, v ) d ( x, v )and energy E E ( t ) = 12 Z Ω v f ( t, x, v ) d ( x, v ) + 12 Z Ω x E ( t, x ) dx are invariants of the analytic solution. We note that in the single species setting mass conservation isequivalent to charge conservation (since charge and mass are proportional) and momentum conservation isequivalent to current conservation.For each of these invariants, a continuity or moment equation, that is posed in physical space only, is satisfiedby an associated density. In the case of mass, the mass density ρ ( t, x ) = Z Ω x f ( t, x, v ) dv satisfies ∂ t ρ ( t, x ) + ∇ x · j ( t, x ) = 0 . Conservation of mass can be easily derived from the continuity equation by integrating in x . It is in factthese continuity equations that our dynamical low-rank integrator satisfies and from which conservation ofthe global invariants is then followed. We note that this is a much stronger result than simply conservingthe invariants, as the present approach also preserves the underlying physical structure of the equation.The associated density for momentum conservation is the momentum density j ( t, x ) = Z Ω v vf ( t, x, v ) dv which satisfies the following continuity equation ∂ t j ( t, x ) + ∇ x · σ ( t, x ) = − E ( t, x ) ρ ( t, x ) , σ ( t, x ) = Z Ω v ( v ⊗ v ) f ( t, x, v ) dv. We can follow conservation of momentum by recognizing that E (1 − ρ ) = ∇ · ( E ⊗ E − E ) and integratingin physical space. Note that due to the normalization of the particle-density function we have R E dx = 0.The energy density e ( t, x ) = 12 Z Ω v v f ( t, x, v ) dv + 12 E ( t, x )satisfies the following continuity equation ∂ t e ( t, x ) + ∇ x · Q ( t, x ) = E ( t, x ) · ( ∂ t E ( t, x ) − j ( t, x )) , Q ( t, x ) = 12 Z Ω v vv f ( t, x, v ) dv. We follow global conservation, i.e. conservation of energy, by ∂ t E ( t, x ) = ∂ t Z e ( t, x ) dx = Z E ( t, x ) · ( ∂ t E ( t, x ) − j ( t, x )) dv = 0 . ∂ t E ( t, x ) = j ( t, x ) is just the electrostatic version of Ampere’s law. A moredetailed discussion can be found in [29]. Looking at the conserved quantities in light of their correspondingcontinuity equations also relates the kinetic equation considered here to their corresponding fluid model. Formore details we refer the reader to [16, 10, 9].
3. The classical dynamical low-rank scheme
For a dynamical low-rank scheme we seek an (approximate) solution of the form f ( t, x, v ) = r X i,j =1 X i ( t, x ) S ij ( t ) V j ( t, v ) , (3)where r is the rank of the approximation. The dynamics is represented by the low-rank factors X i and V j ,that only depend on physical space x and velocity space v , respectively, and the low-rank factors S ij , whichcarries no spatial or velocity dependence.From a physical point of view, the decomposition into physical space and velocity space is natural. If thedynamics is integrable then there are conserved action variables j ( x, v ) and conjugate angle variables θ ( x, v ).In this case, the natural basis functions are Θ k ( θ, t ) = e ikθ and J k ( j, t ) = δ ( j − j k ). Note, however, that thereis a nonlinear coordinate transformation that must be determined in order for this representation to hold.The separation into physical space and velocity used here does hold for linearized equations, such as the Case-Van Kampen eigenfunctions [30] for the linear Vlasov–Poisson equation. For a mathematical perspectivein the context of the dynamical low-rank approximation, see [13]. In the collisional case, the nonlinearcollision operator generates scattering in v and is considered local in x . This simplification would not hold ingeneralized phase space coordinates { θ, j } . Hence, when collisions are dominant, the decomposition betweenphysical space and velocity space is still the natural one. Due to the fact that sharp structures, such as shocksand boundary layers, can form in the x direction, good basis functions are localized in x . For the v direction,good basis functions include functions orthogonal over a Maxwellian distribution and eigenfunctions of thecollision operator. The low-rank approximation can accommodate this without any difficulty, as has beendemonstrated in [9, 8, 10].From now on, we will, for reasons of simplicity, write f for f ( t, x, v ), X i for X i ( t, x ), etc. The equations ofmotions for the low-rank factors X i , S ij , and V j are then determined by imposing a Galerkin condition, see[18]. Following this argument for the PDE case we obtain [11] ∂ t X i X i S ij ! = ( V j D [ f ]) v , (4) ∂ t S ij = ( X i V j D [ f ]) xv , (5) ∂ t X j S ij V j = ( X i D [ f ]) x , (6)where ( f ) x = Z Ω x f dx, ( f ) v = Z Ω v f dv, ( f ) xv = Z Ω f d ( x, v ) . For the Vlasov–Poisson equation, D [ f ] = − v · ∇ x f + E · ∇ v f with f in the form of equation (3). We notethat the low-rank factors X i and V j are orthonormal; that is, they satisfy ( X i X j ) x = δ ij and ( V i V j ) v = δ ij .The main utility of the low-rank approximation is that instead of a problem in dimension 2 d , as in equation(1), we only have to solve 2 r equations of dimension d and r ordinary differential equations. This is thereason why the dynamical low-rank approximation drastically reduces the memory and computational effortrequired to solve such problems, assuming that r is not too large, which is often true in practice. For a moredetailed discussion we refer the reader to [11, 13] and the subsequent discussion in section 4.4quations (4)-(6), however, do not preserve the invariants of the original system. That is, even if the systemis solved exactly (i.e. no time or space error is introduced) conservation of mass, momentum, and energy islost. This is not very surprising as the dynamically low-rank approximation described above does not takethe physical structure of the equations into account. Thus, there is no mechanism that prevents, e.g., massto be removed from the system due to the truncation performed by the low-rank approximation.
4. The conservative dynamical low-rank scheme
In this section, we propose a novel dynamical low-rank integrator with equations of motions, that will takethe place of (4)-(6), that ensure mass, momentum, and energy conservation.For the classical dynamical low-rank integrator we have ∂ t f = X ij ∂ t ( X i S ij ) V j − X ij X i ( ∂ t S ij ) V j + X ij X i ∂ t ( S ij V j )and thus by integrating over v∂ t ρ = X j ( V j D [ f ]) v (1 V j ) v − X ij X i ( X i V j D [ f ]) xv (1 V j ) v + X i X i ( X i D [ f ]) xv , (7)where we have used equations (4)-(6). Now, if we could ensure that 1 ∈ V , where V = span { V j } , then wecould find coefficients c j such that 1 = P j c j V j . Plugging this into equation (7) we get ∂ t ρ = (cid:16) D [ f ] X j c j V j (cid:17) v − X i X i (cid:16) X i D [ f ] X j c j V j (cid:17) xv + X i X i ( X i D [ f ]) xv and thus ∂ t ρ = ( D [ f ]) v − X i X i ( X i D [ f ]) xv + X i X i ( X i D [ f ]) xv = ( D [ f ]) v . By using D [ f ] = − v · ∇ x f + E · ∇ v f we at once obtain the continuity equation ∂ t ρ + ∇ · j = 0 , which would imply mass conservation. However, 1 ∈ V is clearly not true as constant functions do not liein L ( R d ). Thus, the continuity equation is not satisfied and the classic dynamical low-rank integrator doesnot conserve mass. A very similar argument can be made for the momentum, with v ∈ V , and energy,with v ∈ V . One might object at this point that in a numerical simulation we necessarily use a truncateddomain and thus the functions 1, v, and v lie in L (Ω v ). However, on a finite domain the approximationspaces have to be equipped with appropriate boundary conditions (usually either periodic or homogeneousDirichlet conditions are imposed) which essentially leads to the same problem.To obtain a conservative dynamical low-rank integrator we proceed as follows. First, we will use a weightedfunction space that includes the possibility of representing constant functions in V .
However, this on itsown is not yet sufficient. In fact, the dynamic low-rank algorithm automatically chooses appropriate basisfunctions in order to minimize the overall error, according to some Galerkin condition, in the particle-densityfunction f . There is no guarantee that such a choice satisfies 1 ∈ span { V } . Thus, we have to constrain theapproximation space in such a way that constant functions are always part of the basis. This, in particular,requires us to modify the Galerkin condition that is used in the classic integrator. Those two ideas combinedallows us to formulate a conservative dynamical low-rank algorithm. The details of this procedure will bethe topic of the remainder of this section. 5or the conservative dynamical low-rank scheme the solution is approximated by a function in the followingform f ( t, x, v ) = f ( t, x, v ) r X i,j =1 X i ( t, x ) S ij ( t ) V j ( t, v ) , (8)where r is the rank of the approximation and f ( x, v ) = f x ( t, x ) f v ( t, v ) is a, yet to be chosen, weightfunction. The low-rank factors X i ( t, · ) and V j ( t, · ) are assumed to lie in the L spaces weighted by f x and f v , respectively. That is, X i ( t, · ) ∈ L (Ω x , f x ) and V j ( t, · ) ∈ L (Ω v , f v ). The corresponding weightedinner products are denoted by h X i ( t, · ) , X j ( t, · ) i x = Z Ω x f x ( t, x ) X i ( t, x ) X j ( t, x ) dx and h V i ( t, · ) , V j ( t, · ) i v = Z Ω v f v ( t, x ) V i ( t, x ) V j ( t, v ) dv, respectively.In the following we will choose f x ( t, x ) = 1 and f v ( t, v ) = f v ( v ), i.e. f v is time independent. The choice of f v must guarantee that 1, v , and v lie in L (Ω v , f v ) in order to obtain conservation of mass, momentum,and energy, respectively. It is also important to choose the temperature large enough such that the relevantpart of the phase space is captured. For many problems that model plasma instabilities f v ( v ) = exp( − v / v lie in theapproximation space. However, other than the constraints outlined above the choice of f v is arbitrary.The second crucial ingredient of the algorithm is that some of the functions V j are held fixed as the systemevolves in time. We write U a ( v ) = V a ( t, v ) , ≤ a ≤ m and W p ( t, v ) = V p ( t, v ) , m < p ≤ r, where the U a ( v ) are fixed, i.e. they are not changed by the dynamical low-rank integrator, and the W p ( t, v )are allowed to vary in time according to the low-rank algorithm. In the following, indices i , j , k , l span1 , . . . , r , indices p , q span m + 1 , . . . , r , and the indices a , b span 1 , . . . , m . For example, to obtain mass,momentum, and energy conservation in 1+1 dimension and for f v ( v ) = exp( − v /
2) we choose m = 3 with U ( v ) = 1 / k k , U ( v ) = v/ k v k , and U ( v ) = ( v − / k v − k . It is clear that then the orthogonalityconstraint h U a , U b i v = δ ab is satisfied. As usual, we also impose the orthogonality conditions h W p , V j i = δ pj .Note that this implies that the W p are orthogonal both with respect to each other and with respect to the U a .The equations of motion for the low-rank factors are derived by considering the low-rank manifold of allfunctions of a given rank r . In our case this manifold M can be written as M = (cid:26) f ∈ L (Ω , f ) : f ( x, v ) = f v ( v ) X ij X i ( x ) S ij V j ( v ) with invertible S = ( S ij ) ∈ R r × r ,X i ∈ L (Ω x ) , V j ∈ L (Ω v , f v ) with h X i , X j i x = δ ij , h V i , V j i v = δ ij (cid:27) . As usual we will impose the gauge conditions h ∂ t X i , X j i x = 0 and h ∂ t W p , W q i v = 0 in order to make surethat the dynamics of the low-rank factors is uniquely determined. The tangent space of the manifold is then T f M = (cid:26) ˙ f ∈ L (Ω , f ) : ˙ f = f v X ij (cid:0) ˙ X i S ij V j + X i ˙ S ij V j (cid:1) + f v X ip X i S ip ˙ W p , with ˙ S ∈ R r × r , ˙ X i ∈ L (Ω x ) , ˙ V j ∈ L (Ω v , f v ) , and h X i , ˙ X j i x = 0 , h V i , ˙ W q i v = 0 (cid:27) , ∂ t U a = 0. While this construction is similar to the classic dynamicallow-rank integrator, it is crucial that we are cognizant for which low-rank factors the gauge conditions areimposed and that we make sure that the entire set of the V j are orthogonal to each other (and not only the W p ).We now derive the equations of motion for the low-rank factors X i , S ij , and V ij that satisfy the followingPetrov–Galerkin condition (cid:18) vf v , ( ∂ t f − D [ f ]) (cid:19) xv = 0 ∀ v ∈ T f M , (9)where ( f, g ) xv = ( f g ) xv and ∂ t f ∈ T f M . We note that the Petrov–Galerkin condition is different from theclassic dynamical low-rank scheme; an additional factor of 1 /f v has been introduced. The reason for this isthat in the following our goal is to choose an appropriate v ∈ T f M such that the equation of motion for aspecific low-rank factor is isolated. This change in the Petrov–Galerkin condition is what makes this possiblein the weighted approximation spaces that we consider here. Equations for X i : In this case there is little difference to the classic algorithm. We consider a family oftest functions ν k = f v χ ( x ) V k , where χ is an arbitrary function of x . Since we can write ν = f v P ij ˙ X i S ij V j with ˙ X i = χ ( x ) S − ki it clearly holds that ν k ∈ T f M . The Petrov–Galerkin condition (9) then becomes (cid:18) V k χ ( x ) , f v X ij (cid:0) ˙ X i S ij V j + X i ˙ S ij V j (cid:1) + f v X ip X i S ip ˙ W p (cid:19) xv = ( V k χ ( x ) , D [ f ]) xv . We can rewrite this as (cid:28) V k χ ( x ) , X ij (cid:0) ˙ X i S ij V j + X i ˙ S ij V j (cid:1) + X ip X i S ip ˙ W p (cid:29) xv = ( V k χ ( x ) , D [ f ]) xv and thus by using the orthogonality and gauge conditions as well as the fact that χ is arbitrary we obtain X i ˙ X i S ik = ( V k , D [ f ]) v − X i X i ˙ S ik , (10)which is precisely the first equation of motion that is obtained for the classic algorithm (see, e.g., [11]).Now, we can plug the right hand-side of the Vlasov–Poisson equation into right-hand side of equation (10).This yields ( V k , D [ f ]) v = ( V k , − v · ∇ x f + E · ∇ v f ) v = − X ij ( V k , f v v · ( ∇ x X i ) S ij V j ) v + X ij ( V k , X i S ij E · ∇ v ( f v V j )) v = − X ij h V k , vV j i v ∇ x X i S ij + X ij E · ( V k , ∇ v ( f v V j )) v X i S ij = − X ij c kj ∇ x X i S ij + X ij ( c kj · E ) X i S ij , with c kj = h V k , vV j i v , c kj = ( V k , ∇ v ( f v V j )) v . We note that due to the weight function f v the coefficients c kl and c kl are changed compared to the classicalgorithm in [11]. Equations for W p : In this case we consider the family of test functions ν q = f v ζ ( v ) P i X i S iq , where ζ isan arbitrary function in v . Those ν q lie in the tangent space as we can write ν q = f v P ip X i S ip ˙ W p with˙ W p = δ pq ζ ( v ). The Petrov–Galerkin condition (9) then becomes X i (cid:18) ζ ( v ) X i S iq , f v X kl (cid:0) ˙ X k S kl V l + X k ˙ S kl V l (cid:1) + f v X kp X k S kp ˙ W p (cid:19) xv = X i ( ζ ( v ) X i S iq , D [ f ]) xv ζ is arbitrary and using the gauge conditions) X ip S iq S ip ˙ W p = 1 f v X i S iq ( X i , D [ f ]) x − X il S iq ˙ S il V l . (11)On the left-hand side we have the matrix T qp = P i S iq S ip . Since S has full rank the same is true for T andthus we can invert it in order to obtain the equations of motion for the W p .For the Vlasov–Poisson equation we have1 f v ( X i , D [ f ]) x = 1 f v ( X i , − v · ∇ x f + E · ∇ v f ) x = − f v X kl f v h X i , ∇ x X k i x · vS kl V l + 1 f v X kl S kl ∇ v ( f v V l ) · h X i , EX k i x = − X kl ( v · d ik ) S kl V l + 1 f v X kl d ik [ E ] · ∇ v ( f v S kl V l )= − X kl ( v · d ik ) S kl V l + X kl d ik [ E ] · [ ∇ v ( S kl V l ) + ∇ v (log f v ) S kl V l ] , where d ik [ E ] = h X i , EX k i x , d ik = h X i , ∇ x X k i x . Equation for S ij : In this case we consider ν kl = f v X k V l , which lies in the tangent space as ν kl = f v P ij X i ˙ S ij V j with ˙ S ij = δ ki δ jl . The Petrov–Galerkin condition (9) becomes (cid:18) X k V l , f v X ij (cid:0) ˙ X i S ij V j + X i ˙ S ij V j (cid:1) + f v X ip X i S ip ˙ W p (cid:19) xv = ( X k V l , D [ f ]) xv and thus ˙ S kl = ( X k V l , D [ f ]) xv . (12)For the Vlasov–Poisson equation( X k V l , D [ f ]) xv = ( X k V l , − v · ∇ x f + E · ∇ v f ) xv = − X ij h X k , ∇ x X i i x · h V l , vV j i v S ij + X ij h X k , EX i i x · ( V l , ∇ v ( f v V j )) v S ij = − X ij ( d ki · c lj ) S ij + X ij ( d ki · c lj ) S ij . Together equations (10)-(12) are the equations of motions for the proposed conservative dynamical low-rankintegrator. They take the place of equations (4)-(6) that have been used in the classic algorithm. Theprimary difference lies in the equation for the W p , where it has to be ensured that the W p are updated insuch a way that they remain orthogonal to not only the other W p but also the U a . In addition, due to the useof the weighted L spaces all coefficients are changed as well and f v and its derivative make an appearancein the equations of motion.
5. Mass, momentum, and energy conservation for the proposed low-rank approximation
In this section we will show that the numerical scheme derived in section 4 is indeed mass, momentum, andenergy conservative. 8o ensure mass conservation we choose, as is discussed in the previous section, U ∝ K j = P i X i S ij , the density is given by ρ = 1 U K . This can be easily seen as ρ = Z Ω v f dv = 1 U X ij X i S ij h U , V j i v . Using the orthogonality condition h U , V j i = δ j we obtain the desired relation.Deriving the corresponding continuity equation is now straightforward. We have ∂ t ρ = 1 U ∂ t K . From equation (10) we follow ∂ t K = ( V , D [ f ]) v and thus ∂ t ρ = 1 U ( U , D [ f ]) v = Z Ω v D [ f ] dv. This is precisely the relation that is obtained for the continuous evolution, i.e. the solution of the Vlasov–Poisson equation without any low-rank approximation present. Since Z D [ f ] dv = −∇ x · Z Ω v vf dv + E · Z Ω v ∇ v f dv = −∇ x · j we have ∂ t ρ + ∇ x · j = 0 . From the continuity equation conservation of mass follows at once (by integrating in space).For conservation of momentum we proceed in precisely the same manner. For simplicity we only considerthe 1+1 dimensional case here. By choosing U such that v = k v k U , i.e. U ∝ v , momentum is given by j = k v k K . We thus have ∂ t j = k v k ∂ t K = Z Ω v vD [ f ] dv = −∇ x · σ − Eρ, which is the desired continuity equation that implies conservation of momentum. The extension to themulti-dimensional case is straightforward. We have to use one fixed basis function for each of the directionsin which we want to conserve momentum (e.g. U ∝ v , U ∝ v , and U ∝ v ).For energy conservation we can not choose U ∝ v since v is not orthogonal to 1. Thus, we use U ∝ v − . In this setting we have h U , U i v = h U , U i v = 0, as desired. We then can represent v as follows v = k v − k U + k k U . Using this relation we have Z v f dv = X ij X i S ij (cid:0) k v − kh U , V j i v + k kh U , V j i v (cid:1) = k v − k K + k k K . e = k v − k K + k k K + E . Our goal is once again to derive the corresponding continuity equation. We have ∂ t e = k v − k ∂ t K + k k ∂ t K + E · ( ∂ t E )= k v − k ( U , D [ f ]) v + k k ( U , D [ f ]) v + E · ( ∂ t E )= 12 Z v D [ f ] dv + E · ( ∂ t E ) . Evaluating the integral and using the relations for the electric field derived in section 2 we get the desiredcontinuity relation ∂ t e + ∇ x · Q = E · ( ∂ t E − j )which implies conservation of energy.We have here outlined an approach with m = 3, U ∝ U ∝ v , and U ∝ v − U ∝ v and theargument above is even simpler as we do not have to take the orthogonality between U and U into account.We can clearly see that for our algorithm the argument with respect to conservation made for the Vlasov–Poisson equation in section 2 carries over to the low-rank approximation. This is by design and highlightsthe ability of our approach to preserve the physical structure of the Vlasov–Poisson equations.
6. Discretization
In the previous section we have established that the proposed dynamical low-rank integrator conservesmass, momentum, and energy. This is done completely within a continuous formulation. Thus, the onlyapproximation made is due to the fact that the proposed scheme uses a low-rank representation of thesolution. However, to actually implement the conservative dynamical low-rank integrator on a computer, wehave to also introduce a time and space discretization. To devise a numerical method that discretizes theequations of motions for the proposed dynamical low-rank integrator, while maintaining conservation, is themain purpose of this section.In section 6.1 we will introduce an explicit integrator for equations (10)-(12) that preserves mass and mo-mentum up to machine precision. For dynamical low-rank approximations integrators that are robust to thepresence of small singular values have been developed recently [23]. Unfortunately, it turns out that the thisprojector splitting integrator can not be easily adapted to the present situation. In section 6.2 we proposea robust dynamical low-rank integrator that fits within the framework considered in this paper.
Simply applying a classic Runge–Kutta scheme to the evolution equations (10)-(12) does destroy the con-servative properties of the method. As an example, let us consider the classic explicit Euler scheme S n +1 kl = S nkl + τ ( X nk V nl , D [ f n ]) xv ,X n +1 i = X ni + τ X k ( S n ) − ik (cid:20) ( V nk , D [ f n ]) v − X l X nl ( X nl V nk , D [ f n ]) xv (cid:21) ,W n +1 p = W np + τ X q (( S n ) T S n ) − pq (cid:20) f v X i S niq ( X ni , D [ f n ]) x − X il S niq ( X ni V nl , D [ f n ]) xv V nl (cid:21) , τ is the time step size and the upper indices denote the value of the corresponding quantity at thediscrete time t n . The reason why this scheme fails to be conservative is that we use the inverse of S n , i.e. S at time t n , to compute X n +1 . This inconsistency means that there is no well defined K ni and K n +1 i andthus the argument in section 5 can not be applied.However, if we write equation (10) in the following conservative form ∂ t (cid:18)X i X i S ik (cid:19) = ( V k , D [ f ]) v (13)and discretize it using the explicit Euler scheme X i X n +1 i S n +1 ik = X i X ni S nik + τ ( V nk , D [ f n ]) . No change is made to the other two equations of motions. Putting this together we obtain S n +1 kl = S nkl + τ ( X nk V nl , D [ f n ]) xv , (14) X n +1 i = X k ( S n +1 ) − ik (cid:20)X j X nj S njk + τ ( V nk , D [ f n ]) v (cid:21) , (15) W n +1 p = W np + τ X qi (( S n ) T S n ) − pq S niq (cid:20) f v ( X ni , D [ f n ]) x − X l ( X ni V nl , D [ f n ]) xv V nl (cid:21) . (16)In the following we will call this method the conservative Euler scheme. This is still a fully explicit methodas we can first compute S n +1 using equation (14), which is then used in the computation of X n +1 .The scheme satisfies a discrete version of the continuity equation ρ n +1 − ρ n τ = 1 U K n +11 − K n τ = Z Ω v D [ f n ] dv = −∇ x · j n as can be easily seen from equation (15). Integrating this equation in x we get M n +1 = M n and thusconservation of mass.For the momentum we proceed in the same way. We obtain the following discrete continuity equation j n +1 − j n τ = k v k K n +12 − K n τ = Z Ω v vD [ f n ] dv = −∇ x · σ n − E n ρ n which implies conservation of momentum.For energy we have e n +1 − e n τ = k v − k K n +13 − K n τ + k k K n +11 − K n τ + ( E n +1 ) − ( E n ) τ = 12 Z Ω v v D [ f n ] dv + E n E n +1 − E n τ + ( E n +1 − E n ) τ = ∇ x · Q n + E n · (cid:18) E n +1 − E n τ − j n (cid:19) + ( E n +1 − E n ) τ . The reason why we have rewritten the electric field in that particular way should become clear shortly.Integrating the above relation in space yields E n +1 − E n = τ Z Ω x E n · (cid:18) E n +1 − E n τ − j n (cid:19) dx + 12 Z Ω x ( E n +1 − E n ) dx. (17)11e again write the electric field using its potential, i.e. E n = −∇ φ n , to get Z Ω x E n · (cid:18) E n +1 − E n τ − j n (cid:19) dx = Z Ω x φ n (cid:18) ∇ · E n +1 − ∇ · E n τ − ∇ · j n (cid:19) dx = − Z Ω x φ n (cid:18) ρ n +1 − ρ n τ + ∇ · j n (cid:19) dx = 0 . Thus, the first in term on the right-hand side of (17) vanishes. That leaves the second term Z Ω x ( E n +1 − E n ) dx = O ( τ ) , which can clearly not be zero (except for trivial solutions that satisfy E n = const for all n ). Thus, theconservative Euler scheme commits a first order error in the energy. The reason for this is that in thediscrete setting the relation ∂ t ( E ) = 2 E∂ t E does not hold true. This introduces the second term in thecalculations above and thus destroys the conservation of energy in the discrete setting.Fundamentally, the issue is that the Euler scheme is not symmetric under time reversal. If one allows forimplicit methods this deficiency can be remedied. For example, using E n +1 / = ( E n +1 + E n ) / e n +1 − e n τ = ∇ x · Q n − E n +1 / · j n + ( E n +1 − E n )( E n +1 + E n )2 τ = ∇ x · Q n + E n +1 / · (cid:18) E n +1 − E n τ − j n (cid:19) . Integrating in physical space, as above, then shows conservation of energy. However, since the electric fielddepends on the particle-density function this approach results in a fully implicit scheme that has to be solvedup to machine precision, if energy conservation to the same level of accuracy is desired. We consider theconstruction of efficient energy conservative time integrators a subject of future research.Let us now turn our attention to the discretization of space. Fortunately, this is relatively straightforward. Aslong as the discrete approximation of the derivatives and the quadrature rule chosen to define the invariantsallows us to perform integration by parts without introducing any error, the calculations made in thissection carry over to the fully discretized case. This is true for a number of space discretization strategies.For example, using fast Fourier techniques (FFT) or the standard second-order centered finite differencescheme to compute the derivatives in combination with the trapezoidal rule to evaluate the integrals satisfiesthis property. The former will be used in the numerical results in section 7. Moreover, for a number of moreadvanced and higher-order finite difference, finite volume, and discontinuous Galerkin schemes this propertyis satisfied as well (see, e.g., [1, 14, 32]).There is one additional concern in an actual implementation that, by necessity, operates in finite precisionarithmetic (i.e. with doubles or floats on a computer). To show conservation we employed the orthogonalitycondition h W p , U a i = 0 . This is true for the dynamical low-rank integrator by construction (see section 4).The conservative Euler discretization also preserves this property. This can be seen explicitly from equation(16) by taking the inner product with U a h W n +1 p , U a i v = h W np , U a i v + τ X qi (( S n ) T S n ) − pq S niq (cid:28) f v ( X ni , D [ f n ]) x − X l ( X nk V nl , D [ f n ]) xv V nl , U a (cid:29) v . Since h W p , U a i v = 0 is true for the initial value we can assume that h W np , U a i v = 0 and thus h W n +1 p , U a i v = τ X qi (( S n ) T S n ) − pq S niq [( X ni U a , D [ f n ]) xv − ( X nk U a , D [ f n ]) xv ] = 0 . W p fail to be exactly orthogonal to the U a . In the worst case, this can cause a small linear driftin the error of mass and momentum. To avoid this we perform an orthonormalization of the W p at theend of each step. This does not negatively impact the accuracy of the numerical method and the incurredcomputational cost is negligible. The integrator described in the previous section is not robust if the matrix S has small singular values. Thereason for this is that inverting the matrix S is then numerically ill-conditioned. Small singular values occurcommonly if the rank of the solution is lower than the rank chosen to conduct the numerical simulation.We refer to [17] for a more detailed discussion, but note that robustness is generally considered a desirableproperty especially if the rank is adaptively changed during the simulation.For the classic dynamical low-rank approximation a projector splitting integrator has been proposed byLubich & Oseledets [23] that remedies this deficiency. The method has also been extended to a variety oftensor formats [21, 22, 11, 5, 4]. The main utility of the projector splitting integrator is that by using aQR decomposition it avoids the inversion of S . Unfortunately, this projector splitting integrator can not beused in the present situation, because if we solve equation (11) to obtain P ip S iq S ip W p , while holding the X i constant, it is not possible to use a QR decomposition to extract the low-rank factors S and W .Instead of the projector splitting integrator, we consider an approach based on the the recently developedunconventional integrator by Ceruti & Lubich [3]. The main idea is that the dynamical low-rank approxima-tion has two parts. On the one hand, the low-rank factors X i and V j determine the subspaces X = span { X i } and V = span { V j } in which an approximation is sought. However, it does not matter how exactly the X i are chosen as long as they span the same approximation space X . On the other hand, the low-rank factor S contains the coefficients that combine the basis functions in an appropriate way in order to obtain a goodapproximation. Inverting S is only required in equations (10) and (11). That is, it is only required to obtainthe low-rank factors X i and V j .We thus proceed as follows. First, we discretize equation (13) to compute K n +1 k = K nk + τ ( V nk , D [ f n ]) v with K k = X i X i S ik . and discretize equation (11) to compute L n +1 q = L nq + τf v X i S niq ( X ni , D [ f n ]) x − τ X il S niq ( X nk V nl , D [ f n ]) V nl with L n +1 q = X ip S niq S nip W n +1 p and L nq = X ip S niq S nip W np . This does not give us W n +1 p (not even using a QR decomposition) nor X n +1 i . However, the approximationspaces V n +1 and X n +1 are uniquely defined by L n +1 q and K n +1 k , respectively. Thus, we perform a QRdecomposition L n +1 q = X p W n +1 p R pq , K n +1 k = X i X n +1 i R ik to obtain W n +1 p and X n +1 i . The R parts of the QR decomposition, i.e. R pq and R ik , are simply discarded.We then determine S n +1 as follows S n +1 kl = X ij M ki S nij N Tjl + τ (cid:0) X n +1 k V n +1 l , D [ f ( X n +1 , M S n N T , V n +1 )] (cid:1) xv , M ki = h X n +1 k X ni i x , N Tjl = h V nj V n +1 l i v . This procedure is a discretization of equation (12), where we have used that X ij X ni S nij V nj ≈ X kijl X n +1 k M ki S nij N Tjl V n +1 l . That is, we transform the coefficient matrix S nij to the new basis spanned by the previously determined V n +1 j and X n +1 i and then solve equation (12) with the basis held fixed over one time step. This gives us arobust dynamical low-rank integrator as the matrix S needs not be inverted. The downside of this integrator,however, is the the mass and momentum conservation up to machine precision is lost, as we will also see insection 7.
7. Numerical experiments
In this section we illustrate the conservative dynamical low-rank scheme using a number of numerical ex-amples. The numerical algorithm outlined in sections 4 and 6 will be employed. In order to compute thederivatives in space we employ techniques based on the fast Fourier transform (FFT).
As the first example we consider a variant of Landau damping on the domain ( x, v ) ∈ [0 , π ] × [ − , f (0 , x, v ) = (cid:0) α cos( x ) (cid:1) e − v / √ π + X k =1 (cid:15) k cos( k +12 ) v k e − v / , (18)where α = 10 − , (cid:15) k = 10 − for k ∈ { , , } and (cid:15) k = 10 − for k ∈ { , } . This is the classic Landau dampingproblem with an added perturbation that ensures that the initial value is rank 6.The results of a numerical simulation with the conservative Euler dynamical low-rank integrator and rank r = 6 is shown in Figure 1. For this problem the linearized decay rate can be determined analytically andmatches the observed numerical results well. To study the conservation properties we run simulations using m = 0 (no conservation), m = 1 ( U ∝ m = 2 ( U ∝ U ∝ v and thus massand momentum conservation), and m = 3 ( U ∝ U ∝ v , U ∝ v − m = 3 the error in energy is significantly reducedcompared to all the other configurations. Thus, using the energy conservative dynamical low-rank integratoris clearly beneficial, even if the time integration error is taken into account. Reducing the time step sizefurther also improves energy conservation, indicating that the energy error in the m = 3 configuration islimited by the time integration error, as expected. We also note that for the m = 3 configuration the fidelityof the simulation is somewhat worse. The reason for this that the dynamical low-rank integrator is only ableto choose 3 basis functions (all others are already determined by imposing conservation) and has thus lessfreedom to decrease the error in the particle density function (see also the discussion in section 7.3).14 e l e c t r i c e n e r g y m=0m=1m=2m=3analytic decay rate 0 5 10 15 20 25 30time10 m a ss e rr o r m=0m=1m=2m=30 5 10 15 20 25 30time10 m o m e n t u m e rr o r m=0m=1m=2m=3 0 5 10 15 20 25 30time10 e n e r g y e rr o r m=0m=1m=2m=3m=3 ( =1e-5)m=3 ( =1e-6) Figure 1: Time evolution of the electric energy (top-left), error in mass (top-right), momentum (bottom-left), and energy(bottom-right) for the linear Landau damping problem with initial value (18) is shown. For mass and energy the relative andfor momentum the absolute error is reported. The simulation is conducted using the conservative Euler dynamical low-rankintegrator with step size τ = 10 − (except if otherwise indicated), rank r = 6, and 128 grid points in both the spatial andvelocity direction are used. The configurations considered are m = 0, m = 1 ( U ∝ m = 2 ( U ∝ U ∝ v ), and m = 3( U ∝ U ∝ v , U ∝ v − As our second example, we consider a Maxwellian particle density f (0 , x, v ) = ρ √ π exp (cid:0) − ( v − u ( x )) / (cid:1) (19)with a position-dependent velocity u ( x ) = α cos( x ) on the domain ( x, v ) ∈ [0 , π ] × [ − , ρ = 1 + (cid:15) cos( x ). In the simulation the parameters are chosen as α = 0 . (cid:15) = 10 − . Note thatthe initial value given in (19) is not low-rank. However, for small to moderate u , as is the case here, we canuse the following expansion (see, e.g., [6, 9])exp (cid:0) − ( v − u ( x )) / (cid:1) = exp( − v / (cid:0) uv + u ( v −
1) + u ( v − v ) + u (3 − v + v )+ u (15 v − v + v ) (cid:1) + O ( u ) . Thus, for the initial value of the dynamical low-rank integrator we use the following rank 6 approximation f (0 , x, v ) = ρ √ π exp( − v / (cid:0) uv + u ( v −
1) + u ( v − v ) + u (3 − v + v )+ u (15 v − v + v ) (cid:1) (20)The results of a numerical simulation with the conservative Euler dynamical low-rank integrator and rank r = 6 is shown in Figure 1. We again observe excellent agreement between theory and the numerical results.15n particular, mass and momentum are conserved up to machine precision (and are between 9 and 14 ordersof magnitude smaller than for the m = 0 configuration) and the error in energy is reduced and shown to bedominated by the time integration error. e rr o r i n t h e e l e c t r i c e n e r g y m=0m=1m=2m=3 0 2 4 6 8 10time10 m a ss e rr o r m=0m=1m=2m=30 2 4 6 8 10time10 m o m e n t u m e rr o r m=0m=1m=2m=3 0 2 4 6 8 10time10 e n e r g y e rr o r m=0m=1m=2m=3m=3 ( =1e-5)m=3 ( =1e-6) Figure 2: Time evolution of the error in electric energy (top-left), error in mass (top-right), error in momentum (bottom-left),and error in energy (bottom-right) for the Maxwellian with position-dependent velocity (20) is shown. To determine the errorin the electric energy a reference solution is computed using a full grid semi-Lagrangian scheme that uses 1024 degrees offreedom in both the x and v direction. For mass and energy the relative and for momentum the absolute error is reported. Thesimulation is conducted using the conservative Euler dynamical low-rank integrator with step size τ = 10 − (except if otherwiseindicated), rank r = 6, and 128 grid points in both the spatial and velocity direction are used. The configurations consideredare m = 0, m = 1 ( U ∝ m = 2 ( U ∝ U ∝ v ), and m = 3 ( U ∝ U ∝ v , U ∝ v − The conservative Euler dynamical low-rank integrator can not be used if the rank with which the simulationis run is larger than the rank of the solution (see the discussion in section 6.2). This is often the case for anumber of commonly considered plasma instabilities, where the initial value has rank 1. However, for theunconventional dynamical low-rank integrator (see section 6.2) this is not an issue. To demonstrate this wewill consider the following linear Landau damping problem f (0 , x, v ) = (cid:0) α cos( x ) (cid:1) e − v / √ π , (21)which is rank 1. For the simulation we use α = 10 − on the domain ( x, v ) ∈ [0 , π ] × [ − , r = 10 is shown in Figure 3. We observe that the analytic decayrate is reproduced accurately by all configurations. In particular, no reduction in fidelity is observed for the m = 3 configuration (as is the case in section 7.1). Thus, having a certain number of basis functions that thealgorithm can choose freely is clearly beneficial. As explained in section 6.2, the unconventional integratoris not conservative up to machine precision. The error in mass, momentum, and energy is dominated by thetime integration error and thus reducing the time step size reduces the error in the conserved quantities, asis illustrated for mass in Figure 3. 16 e l e c t r i c e n e r g y m=0m=1m=2m=3analytic decay rate 0 5 10 15 20 25 30 35time10 e rr o r mass errormomentum errorenergy errormass error ( =1e-4) Figure 3: Time evolution of the electric energy (left) and error in mass, momentum, and energy (right) for the linear Landaudamping problem with initial value (21) is shown. For mass and energy the relative and for momentum the absolute error isreported. The simulation is conducted using the unconventional dynamical low-rank integrator with step size τ = 10 − , rank r = 10, and 128 grid points in both the spatial and velocity direction are used. The configurations considered are m = 0, m = 1( U ∝ m = 2 ( U ∝ U ∝ v ), and m = 3 ( U ∝ U ∝ v , U ∝ v − Finally, we consider the two-stream instability f (0 , x, v ) = √ π (cid:16) e − ( v − ¯ v ) / + e − ( v +¯ v ) / (cid:17) (1 + α cos( kx )) (22)with α = 10 − , k = 0 .
2, and ¯ v = 2 . x, v ) ∈ [0 , π ] × [ − , e l e c t r i c e n e r g y m=0m=1m=2m=3full grid solution 0 5 10 15 20 25 30 35 40time10 e rr o r mass errormomentum errorenergy errormass error ( =1e-4)mass error ( =1e-5) Figure 4: Time evolution of the electric energy (left) and error in mass, momentum, and energy (right) for the two-streaminstability (22) is shown. For mass and energy the relative and for momentum the absolute error is reported. The simulation isconducted using the unconventional dynamical low-rank integrator with step size τ = 10 − , rank r = 10, and 128 grid pointsin both the spatial and velocity direction are used. The configurations considered are m = 0, m = 1 ( U ∝ m = 2 ( U ∝ U ∝ v ), and m = 3 ( U ∝ U ∝ v , U ∝ v − x and v direction. cknowledgements Work by I. Joseph was performed under the auspices of the U.S. DOE by LLNL under Contract DE-AC52-07NA27344.
References [1] A. Arakawa. Computational design for long-term numerical integration of the equations of fluid motion-two-dimensional incompressible flow.
J. Comput. Phys. , 135(2), 1997.[2] E. Camporeale, G.L. Delzanno, B.K. Bergen, and J.D. Moulton. On the velocity space discretizationfor the Vlasov–Poisson system: Comparison between implicit Hermite spectral and Particle-in-Cellmethods.
Comput. Phys. Commun. , 198:47–58, 2016.[3] G. Ceruti and C. Lubich. An unconventional robust integrator for dynamical low-rank approximation. arXiv:2010.02022 , 2020.[4] G. Ceruti and C. Lubich. Time integration of symmetric and anti-symmetric low-rank matrices andTucker tensors.
BIT Numer. Math. , pages 1–24, 2020.[5] G. Ceruti, C. Lubich, and H. Walach. Time integration of tree tensor networks. arXiv:2002.11392 ,2020.[6] S. Chen and G.D. Doolen. Lattice Boltzmann method for fluid flows.
Annu. Rev. Fluid Mech. , 30(1):329–364, 1998.[7] N. Crouseilles, L. Einkemmer, and E. Faou. Hamiltonian splitting for the Vlasov–Maxwell equations.
J. Comput. Phys. , 283:224–240, 2015.[8] Z. Ding, L. Einkemmer, and Q. Li. Error analysis of an asymptotic preserving dynamical low-rankintegrator for the multi-scale radiative transfer equation. arXiv:1907.04247 , 2019.[9] L. Einkemmer. A low-rank algorithm for weakly compressible flow.
SIAM J. Sci. Comput. , 41(5):A2795–A2814, 2019.[10] L. Einkemmer, J. Hu, and Y. Wang. An asymptotic-preserving dynamical low-rank method for themulti-scale multi-dimensional linear transport equation. arXiv:2005.06571 , 2020.[11] L. Einkemmer and C. Lubich. A low-rank projector-splitting integrator for the Vlasov–Poisson equation.
SIAM J. Sci. Comput. , 40:B1330–B1360, 2018.[12] L. Einkemmer and C. Lubich. A quasi-conservative dynamical low-rank algorithm for the Vlasov equa-tion.
SIAM J. Sci. Comput. , 41(5):B1061–B1081, 2019.[13] L. Einkemmer, A. Ostermann, and C. Piazzola. A low-rank projector-splitting integrator for the Vlasov–Maxwell equations with divergence correction.
J. Comput. Phys. , 403:109063, 2020.[14] L. Einkemmer and M. Wiesenberger. A conservative discontinuous Galerkin scheme for the 2D incom-pressible Navier–Stokes equations.
Comput. Phys. Commun. , 185(11):2865–2873, 2014.[15] E. Hairer, C. Lubich, and G. Wanner.
Geometric numerical integration: structure-preserving algorithmsfor ordinary differential equations , volume 31. Springer Science, 2006.[16] J. Hu, S. Jin, and Q. Li. Asymptotic-Preserving Schemes for Multiscale Hyperbolic and Kinetic Equa-tions. In
Handbook of Numerical Analysis , pages 103–129. 2017.1817] E. Kieri, C. Lubich, and H. Walach. Discretized dynamical low-rank approximation in the presence ofsmall singular values.
SIAM J. Numer. Anal. , 54(2):1020–1038, 2016.[18] O. Koch and C. Lubich. Dynamical tensor approximation.
SIAM J. Matrix Anal. Appl. , 31:2360–2375,2010.[19] K. Kormann. A semi-Lagrangian Vlasov solver in tensor train format.
SIAM J. Sci. Comput. , 37:613–632, 2015.[20] K. Kormann and E. Sonnendrücker. Sparse grids for the Vlasov–Poisson equation.
Sparse Grids Appl. ,pages 163–190, 2014.[21] C. Lubich, I. V. Oseledets, and B. Vandereycken. Time integration of tensor trains.
SIAM J. Numer.Anal. , 53:917–941, 2015.[22] C. Lubich, B. Vandereycken, and H. Walach. Time integration of rank-constrained Tucker tensors.
SIAM J. Numer. Anal. , 56:1273–1290, 2018.[23] Christian Lubich and Ivan V. Oseledets. A projector-splitting integrator for dynamical low-rank ap-proximation.
BIT , 54:171–188, 2014.[24] E. Musharbash, F. Nobile, and T. Zhou. Error analysis of the Dynamically Orthogonal approximationof time dependent random PDEs.
SIAM J. Sci. Comput. , 37(2):A776–A810, 2015.[25] A. Nonnenmacher and C. Lubich. Dynamical low-rank approximation: applications and numericalexperiments.
Math. Comput. Simul. , 79(4):1346–1357, 2008.[26] M. Palmroth, U. Ganse, Y. Pfau-Kempf, M. Battarbee, L. Turc, T. Brito, M. Grandin, S. Hoilijoki,A. Sandroos, and S. von Alfthan. Vlasov methods in space physics and astrophysics.
Living Reviews inComputational Astrophysics , 4(1):1, 2018.[27] Z. Peng, R. McClarren, and M. Frank. A low-rank method for two-dimensional time-dependent radiationtransport calculations. arXiv, 1912.07522 , 2019.[28] Z. Peng and R.G. McClarren. A high-order/low-order (HOLO) algorithm for preserving conservationin time-dependent low-rank transport calculations. arXiv:2011.06072 , 2020.[29] E. Sonnendrücker. Numerical methods for the Vlasov equations. , 2013.[30] N.G. Van Kampen. On the theory of stationary waves in plasmas.
Physica , 21(6-10):949–963, 1955.[31] J. P. Verboncoeur. Particle simulation of plasmas: review and advances.
Plasma Phys. Control. Fusion ,47(5A):A231, 2005.[32] M. Wiesenberger, L. Einkemmer, M. Held, A. Gutierrez-Milla, X. Saez, and R. Iakymchuk. Repro-ducibility, accuracy and performance of the Feltor code and library on parallel computer architectures.