Semi-Lagrangian
Contents
Semi-Lagrangian¶
Meshes¶
Description¶
The mesh types aim at storing array data, but including additional information needed to put this data in context. There are several types of meshes:
- cartesian mesh
an N-D fully regular mesh defined by an N-D data array and additional parameters that define the domain (i.e.: xmin, xmax) and other parameters like the number of cells in each dimension and the spacing of the data (i.e.: delta). The services provided should be like (with corresponding name changes, indicating the dimensionality of the data):
allocation(new) and initialization functions,
deletion,
a copy constructor,
computation of data interpolants (cubic splines for example),
get_value(mesh,x1,x2,...)
: wherex1
,x2
,… belong to the corresponding interval(x1_min, x1_max)
,(x2_min, x2_max)
, etc., and which returns the interpolated value at the desired point. This operation launches spline interpolations under the hood.get_node_value(mesh,i,j,...)
: analogous toget_value()
but does not need to launch any interpolation as the indices are integers and the requested value falls on a mesh node. This is implemented by a macro.set_node_value(mesh,i,j,...,val)
: sets the node datum at i,j,… with value. Implemented with a macro.
There are some pitfalls with the suggested interface:
the naming convention could get a little complicated depending on the type of data stored in the array (scalar, multiple-valued, etc.), as we had originally intended with the field/vec naming convention.
The interface may need to directly expose its underlying data, or pointers to sections of it for certain operations, like FFT. The whole data field of a mesh could need to be set to a whole data array in one step, during the initialization, such as after a remap operation. At least in these cases, the access to the data might be safer given the nature of the operations.
- structured mesh
a structured mesh is a mapped cartesian mesh in which the coordinates are decoupled. An ND mesh data is stored as an ND array. In addition, however, we need N 1D arrays to store the actual coordinates of the node locations in each dimension. As an alternative, we could use N functions
, , … , , to represent each transformation. We follow the convention that , , etc. are values in the . The offered services ought to be:allocation, initialization, deletion and copy.
get_node_value(mesh, i, j, ... )
: which reads directly from the data array.set_node_value(mesh, i, j, ... )
: which writes directly to the data array.In an analogous way to the
get_value()
functions described above, we could provide something likeget_value(mesh,
,
)
. This would also trigger an interpolation step using uniform splines generated with the ND data. This would require the user to be thinking in terms of the logical variables .Similar functions could be provided such that the user could also request values at points
. For this, we could launch non-uniform splines, with the spacing determined by the arrays and the ND data.This type may also need to grant direct access to the data array for use in operations like FFT and similar.
- tensor product mesh
This is a mapped cartesian mesh in which the coordinates are coupled. That is, the mappings have the form:
, , … , . The specification of this type of mesh requires an ND data array, N ND arrays that specify the transformations numerically, or N functions of arity N to specify the transformations. The services provided should be:allocation, initialization, deletion and copy.
get_node_value(mesh, i, j, ... )
: which reads directly from the data array.set_node_value(mesh, i, j, ... )
: which writes directly to the data array.get_value(mesh,
,
, ... )
which would also trigger uniform spline interpolations, andget_value(mesh, x1, x2, ... )
, which would trigger nonuniform spline interpolations.
We may need to include other operations here which would entail inverse mappings (maybe also with splines), or NURBS or something else. Need to fill in the details more here.
Scalar Fields¶
Description¶
The physical quantities of interest are normally defined in a physical
space. For instance, the electric potential
(place picture of an example coordinate transformation here.)
Consider a scalar field and the services that it should offer to permit us to use it in our logical grid:
The scalar field is a derived type.
The type contains a simple array to store the data, i.e.: the values of the fields on a collection of points. Both, logical mesh values or physical mesh values can be stored in the same array.
The type contains an object that fully specifies the associated coordinate transformation (a mapped mesh). The transformation should provide all the needed services to permit moving the representation of the data from the physical space to the logical space and vice-versa:
something
something else
Quasi-Neutral Equation Solver¶
Description¶
Here we present a simplified but hopefully reasonably clear step-by-step derivation of the quasi-neutral equation (the gyrokinetic Poisson equation) which this module solves. The intent is to make clear some of the assumptions that are built into this model. We follow the general argument given by Krommes (cite reference here for “Nonlinear gyrokinetics: a powerful tool for the description of microturbulence in magnetized plasma”).
The model is based on the idea that the fast gyrations of particles around the magnetic field lines can be averaged away while still preserving the most important long-term physics (hence the name of this type of approach: gyrokinetics). The simplest model that we will look into first deals directly in the distribution function of the gyrocenters, instead of the particles’.
We start with the Poisson equation:
where
where the indices
Here
Here,
Here we have introduced the assumption of
we can recast equation (3) into
where
is called the dielectric constant of the gyrokinetic vacuum. With equation (4), and neglecting the polarization of the electrons, the modified Poisson equation (1) can be written as:
In fusion applications,
The charge density of the electrons also receives a special treatment. The basic assumption here is that the electrons can move very quickly along a magnetic field line and thus are able to rapidly react to electric potential variations through changes in the electron particle density. Thus, it is assumed that along a field line, the electrons obey the Boltzmann relation:
In equation (5),
At the risk of being too sloppy, we will equate the electron particle density with the electron gyrocenter density. A more careful step would involve gyroaveraging directly the original electron distribution function. With this assumption, our modified poisson becomes:
Rearranging terms and using the relation:
we arrive at:
But for a normalization of the variables, equation
(8) is virtually the same as the one stated in
the report by Latu (list here reference for “Scalable Quasineutral
solver for gyrokinetic simulation”). Equation
(8) is not yet ready to solve due to the
presence of the
By making two final assumptions: that the quantities involved in the
calculation of
The average ion charge density can be computed from
Once we calculate
Exposed Interface¶
Fundamental type: None. It is a function that operates on other top-level types. Function:
sll_solve_quasi_neutral_equation( electron_T_profile_2D,
initial_rho_ion_profile_2D,
charge_density,
phi )
Usage¶
Advection Field¶
Description¶
Exposed Interface¶
Fundamental type:
sll_advection_field_3D_t
This implies that one of the options is to have multiple representations, for 3D, 2D, 1D.
Usage¶
Advection¶
Description¶
Exposed Interface¶
Fundamental type: None. This is a function that operates on multiple top-level types. Function:
sll_advect( distribution_function,
advection_field,
dt,
space_mesh
scheme )
Above, scheme
is the functional parameterization of the various
methods in use (PSM, BSL, …) and for which we need a standardized
interface. The above assumes that we can devise a standard functional
interface.