Advection methods

The advection/ folder gathers the backward semi lagrangian scheme classes. There are two main classes, AdvectionSpatial and AdvectionVelocity. Implementing these operators separately makes sense, since we are using time splitting.

The feet of the characteristic curves are used to interpolate the updated distribution function on mesh points. It uses batched interpolators so the interpolation step is done over the whole distribution function.

Spatial advection

Here the purpose is the advection along a direction on the physical space dimension of the phase space. The dynamics of the motion on the spatial dimension are governed by the following equation.

\frac{d f_{s}}{d t} = \sqrt{\frac{m_{e}}{m_{s}}} v \frac{\partial f_{s}}{\partial x}

Velocity advection

Here the purpose is the advection along a direction on the velocity space dimension of the phase space. The dynamics of the motion on the velocity dimension are governed by the following equation, where E is the electric field.

\frac{d f_{s}}{d t} = q_{s} \sqrt{\frac{m_{e}}{m_{s}}} E \frac{\partial f_{s}}{\partial v}

1D advection with a given advection field

The purpose of the BslAdvection1D operator is an advection along a given direction of the phase space. The advection field is given as input. The dynamics of the motion are governed by the following equation.

\partial_{t} f (t, x) + A \partial_{x_{i}} (x^{'}) f (t, x) = 0, x \in Ω, x^{'} \in Ω^{'},

with

$f$ the function to advect. It is defined on $Ω$ domain and the time dimension;
$A$ the advection field. It could be defined on a subdomain $Ω^{'} \subset Ω$ ;
and $x_{i}$ a given direction. The advection field domain has to be defined on this dimension for the time integration method.

Example of use

Here are some examples of equation types the BslAdvection1D operator can solve:

Equation on a 1D domain (and time dimension):

\partial_{t} f (t, x) + A_{x} (x) \partial_{x} f (t, x) = 0, x \in Ω,

Here $Ω^{'} = Ω \in R$ . In the code, it would correspond to

DFieldX f(idx_range_x);
DFieldX A(idx_range_x);
using IDimInterest = IDimX;

Equation on a 2D domain (and time dimension):

\partial_{t} f (t, x, y) + A_{x} (x, y) \partial_{x} f (t, x, y) = 0,

Here $Ω^{'} = Ω \in R^{2}$ . In the code, it would correspond to

DFieldXY f(idx_range_xy);
DFieldXY A(idx_range_xy);
using IDimInterest = IDimX;

Equation on a 2Dx2V domain (and time dimension):

\partial_{t} f (t, x, y, v_{x}, v_{y}) + A_{x} (x, y) \partial_{x} f (t, x, y, v_{x}, v_{y}) = 0,

Here $Ω^{'} \in R^{2}$ and $Ω \in R^{4}$ . In the code, it would correspond to

DFieldXYVxVy f(idx_range_xyvxvy);
DFieldXY A(idx_range_xy);
using IDimInterest = IDimX;

Equation on a 1Dx1V domain (and time dimension):

\begin{aligned} \partial_{t} f (t, x, v_{x}) + A_{x} (x, v_{x}) \partial_{x} f (t, x, v_{x}) = 0, with for instance, A_{x} (x, v_{x}) = v_{x}, \\ or \partial_{t} f (t, x, v_{x}) + A_{v_{x}} (x, v_{x}) \partial_{v_{x}} f (t, x, v_{x}) = 0, with for instance, A_{v_{x}} (x, v_{x}) = E (x), \end{aligned}

Here $Ω^{'} = Ω \in R^{2}$ . In the code, it would correspond to

DFieldXVx f(idx_range_xvx);
DFieldXVx A(idx_range_xvx);
using IDimInterest = IDimX;

or

using IDimInterest = IDimVx;

Equation on a 1Dx1V domain with species dimension (and time dimension):

\begin{aligned} \partial_{t} f_{s} (t, x, v_{x}) + A_{s, x} (x) \partial_{x} f_{s} (t, x, v_{x}) = 0, \end{aligned}

Here $Ω^{'} = Ω \in R^{2}$ . In the code, it would correspond to

DFieldSpXVx f(idx_range_sp_xvx);
DFieldSpX A(idx_range_sp_x);
using IDimInterest = IDimX;

Parameters

The operator takes as templated parameters:

IDimInterest: a dimension of interest (or advection dimension) which refers to the dimension along which we advect the function;
AdvectionDomain: an advection domain, which refers to the domain where the advection field is defined. It has to contain the dimension of interest for the interpolation of the advection field in the time integration method;
FunctionDomain: the full domain where the function we want to advect is defined;
AdvectionFieldBuilder: a spline builder type for the advection field.
AdvectionFieldEvaluator: a spline evaluator type for the advection field.
TimeStepper: a time integration method (see timestepper) to solve the characteristic equation. It has to be defined on the advection field domain. The feet have to be a Field of coordinates of the dimension of interest defined on the advection field domain.

Remark/Warning: the BslAdvection1D operator is built with builder and evaluator for the advection field and interpolator for the function we want to advect. Theses operators have to be defined on the same domain as the advection field and function. For instance, if the advection field and/or the function are defined on the species dimension, then the interpolators have to contain the species dimension in its batched dimensions (see tests in the tests/advection/ folder).

Remark/Warning: The advection field need to use interpolation on B-splines. So we cannot use other type of interpolator for the advection field. However there is no constraint on the interpolator of the advected function.

PolarFootFinder

These methods are designed to calculate the foot of the characteristic on the polar plane for a 2D transport equation of the type:

\partial_{t} f (t, x (r, θ), y (r, θ)) + A (t, x (r, θ), y (r, θ)) \cdot \nabla f (t, x (r, θ), y (r, θ)) = 0,

with $f (0, x, y) = f_{0} (x, y)$ and A the advection field.

The characteristics are the solutions $X$ and $Y$ of the equations:

\begin{array}{r} \partial_{t} X (t; s, x, y) = A_{x} (t, X (t; s, x, y), Y (t; s, x, y)), \\ \partial_{t} Y (t; s, x, y) = A_{y} (t, X (t; s, x, y), Y (t; s, x, y)), \\ X (s; s, x, y) = x, \\ Y (s; s, x, y) = y . \end{array}

The characteristic $X$ represents the trajectory on the $x$ -dimension of the solution $f$ (idem for $Y$ on the $y$ -dimension). The parametrisation of the trajectory are given after the ";". $(X (t; s, x, y), Y (t; s, x, y))$ indicates that at the time $s$ , the trajectory passes by the point $(x, y)$ . In the backward semi-Lagrangian method, we solve the equation of the characteristics to identify the trajectory of the solution passing by each mesh point $(x, y)_{i j}$ at the time $s = t + d t$ . We are interested in its position at the previous time step $t$ . The conservation property along the characteristics informs us that the value of the function at this position $(X (t; s = t + d t, x, y), Y (t; s = t + d t, x, y))$ at $t$ is the same as the value of the function at the mesh point $(x, y)_{i j}$ at the time $s = t + d t$ .

The characteristic feet are calculated using a time integration method. For details of available methods see timestepper.

Advection domain

There are two advection domains to consider:

the physical domain;
the pseudo-Cartesian domain.

The logical domain is not used as it would be hard to calculate the feet close to the O-point in this domain.

It seems logical to use the physical domain, where the studied equation is given, as the advection domain.

However, we want to solve this equation on a polar grid. So before advecting, we have to compute the mesh points in the physical domain using a mapping function $F$ :

F : (r, θ)_{i, j} \mapsto (x, y)_{i, j} .

This adds some steps to the advection operator, we now have to compute

the mesh points in the physical domain using $F$ ;
the characteristic feet in the physical domain;
the characteristic feet in the logical domain (polar grid) using $F^{- 1}$ ;
then interpolate the advection function at the characteristic feet in the logical domain.

The third step can be difficult especially if the mapping function $F$ is not analytically invertible. It is not impossible but the computations can be costly.

That is why, we introduce the pseudo-Cartesian domain. We use another mapping function $G$ such that:

G : (r, θ)_{i, j} \mapsto (x, y)_{i, j} = (r \cos (θ), r \sin (θ))_{i, j} .

Then the four previous steps become

calculate the mesh points in the pseudo-Cartesian domain using $G$ ;
calculate the advection field $A$ in the pseudo-Cartesian domain using the Jacobian matrix of $(F \circ G^{- 1})^{- 1}$ ;
calculate the characteristic feet in the pseudo-Cartesian domain;
calculate the characteristic feet in the logical domain (polar grid) using $G^{- 1}$ ;

Here, $G$ is analytically invertible (we can fix $G^{- 1} (x = x_{0}, y = y_{0}) = (r = 0, θ = 0)$ ) and $(J_{F} J_{G}^{- 1})^{- 1}$ is well-defined. The details are given in Edoardo Zoni's article [1].

Remark 1: if $F$ is the circular mapping function, then the physical domain and the pseudo-Cartesian domain are the same.

Remark 2: if the mapping function is analytically invertible, it is less costly to advect in the physical domain.

References

[1] Edoardo Zoni, Yaman Güçlü. "Solving hyperbolic-elliptic problems on singular mapped disk-like domains with the method of characteristics and spline finite elements". (https://doi.org/10.1016/j.jcp.2019.108889.) Journal of Computational Physics (2019).