In a similar way, the acceleration update profits from the structure of the electromagnetic forces acting on the particles. In the nonrelativistic Vlasov equation, the action of the Lorentz force transforms the phase space distribution inside one simulation time step like a solid rotator in velocity space: the magnetic field causes a gyration of the distribution function by the Larmor motion
while the electric field causes (1) acceleration and (2) drift motion perpendicular to the magnetic field direction.
Hence, the overall velocity space transformation is a rotation around the local magnetic field direction, with the gyration center given by the drift velocity. This transformation can be decomposed into three successive shear motions along the coordinate axes, thus repeating the same fundamental algorithmic structure as the spatial translation update. The acceleration update affects the velocity space completely locally within each spatial simulation cell.
In Vlasiator 5, the default solver uses a 5th-order interpolation for the acceleration of w.r.t. , a 3rd-order interpolation for the translation of w.r.t. . The Strang splitting scheme in time is 2nd order. (Actually this part is not hard to implement. Strangely from my advection test results, Vlasiator is only 1st-order in space for translation and 1st-order in time?)
Interestingly, the same idea appears in theoretical plasma physics for studying the phase space evolution. An extremely hard to solve problem can become quite easy and straightforward by moving in the phase space and shift your calculation completely to another spatial-temporal location.
Let us write it down in a more rigorous way. Let be the spatial translation operator for advection
and be the acceleration operator including rotation
The splitting is performed using a leapfrog scheme:
can generally be described by an offset 3D rotation matrix (due to gyration around ). As every offset rotation matrix can be decomposed into three shear matrices , each performing an axis-parallel shear into one spatial dimension, the numerically efficient semi-Lagrangian acceleration update using the SLICE-3D approach is possible: before each shear transformation, the velocity space is rearranged into a single-cell column format parallel to the shear direction in memory, each of which requires only a 1D remapping with a high reconstruction order. This update method comes with a maximum rotation angle limit due to the shear decomposition of about 22∘ (WHY???), which imposes a further time step limit. For larger rotational angles per time step (caused by stronger magnetic fields), the acceleration can be subcycled.
can generally be described by a translation matrix with no rotation, and the same SLICE-3D approach lends itself to it in a similar vein as for velocity space. The main difference is the typical use of 3rd order reconstruction in order to keep the stencil width at two.
The magnetic field is propagated using the algorithm by Londrillo and del Zanna (2004), which is a 2nd-order accurate (both in time and space) upwind constrained transport method. All reconstructions between volume-, face- and edge-averaged values follow Balsara+ (2009). The face-averaged values of are propagated forward in time using the integral form of Faraday's law
where the integral on the LHS is taken over the cell face and the line integral is evaluated along the contour of that face. Once has been computed based on Ohm's law, it is easy to propagate using the discretized form of the above. When computing each component of on an edge, the solver computes the candidate values on the four neighboring cells of the edge. In the supermagnetosonic case, when the plasma velocity exceeds the speed of the fast magnetosonic wave mode, the upwinded value from one of the cells is used. In the submagnetosonic case the value is computed as a weighted average of the electric field on the four cells and a diffusive flux is added to stabilize the scheme.
For the time integration a order Runge–Kutta method is used. To propagate the field from to we need the and values at both and . With the leap-frog algorithm there are no real values for the distribution function at these times. A order accurate interpolation is used to compute the required values:
The field solver also contributes to the dynamic computation of the timestep. With the simplest form of Ohm's law the fastest characteristic speed is the speed of the fast magnetosonic wave mode and we use that speed to compute the maximum timestep allowed by the field solver. For the field solver Courant numbers is used, as higher values cause numerical instability of the scheme.
In Vlasiator the magnetic field has been split into a perturbed field updated during the simulations, and a static background field. The electric field is computed based on the total magnetic field and all changes to the magnetic field are only added to the perturbed part of the magnetic field. The background field must be curl-free and thus the Hall term can be computed based on the perturbed part only. This avoids numerical integration errors arising from strong background field gradients. In magnetospheric simulations the background field consists of the Earth's dipole, as well as a constant IMF in all cells. As can be seen in ldz_main.cpp, to handle the potential stiffness of the generalized Ohm's law, the subcycling technique (typically within 50 cycles) is applied for all terms on the RHS. Compared this to the MHD treatment in BATSRUS: BATSRUS solves the advection term and electron pressure gradient term are marched with explicit scheme, while the Hall term is marched with implicit GMRES scheme (typically ~ 10 steps). This lies in the fact that the Hall term is mathematically stiff, while other terms are nonstiff. This makes me wondering if it is possible to use different treatment for different terms in Vlasiator as well.
When adapting to AMR, Vlasiator developers made the decision of first keeping the same field solver on the finest level only. This makes the efficient implementation easier, at the cost of memory usage.
A typical ion population is localized in velocity space, for example a Maxwellian distribution, and a large portion of velocity space is effectively empty. A key technique in saving memory is the sparse velocity grid storage. The velocity grid is divided into velocity blocks comprising 4x4x4 velocity cells. The sparse representation is done at this level: either a block exists with all of its 64 velocity cells, or it does not exist at all. In the sparse representation we define that a block has content if any of its 64 velocity cells has a density above a user-specified threshold value. A velocity block exists if it has content, or if it is a neighbor to a block with content in any of the six dimensions. In the velocity space all 26 nearest neighbors are included, while in ordinary space blocks within a distance of 2 in the 6 normal directions and a distance of 1 in nodal directions are included.
There are two sources of loss when propagating the distribution function:
fluxes that flow out of the velocity grid;
distributions that go below the storage threshold.
Usually the term dominates.
|||Eulerian and Lagrangian descriptions of field appear most notably in fluid mechanics.|
|||This can be easily observed with a small velocity space, where the moments calculated from the distribution functions deviates from the analytical values.|