How to write a Propagator#

A Propagator is where the magic of your numerical method happens. It is where you implement the numerical update rules for variables of your model.

Basics#

Each Propagator must have at least the two magic methods __init__(self, ...) and __call__(self, dt). The former defines the variables to be updated and instantiates all necessary objects and data structures used during the simulation. The latter takes the time step dt as an argument and updates the variables by one time step.

The __init__(self, args, *, kwargs) constructor can take arguments args and keyword arguments kwargs, separated by a star *. The args are the variables to be updated, whereas the kwargs are parameters or options needed in the update step. The latter usually come from the input parameter file, as described in 5. Add Propagators. As an example, let us look at struphy.propagators.propagators_coupling.VlasovAmpere:

def __init__(self,
    e: BlockVector,
    particles: Particles6D,
    *,
    c1: float = 1.,
    c2: float = 1.,
    solver=options(default=True)['solver'],
    ):

Here, a variable e of type BlockVector and a variable particles of type Particles6D are updated by the Propagator. Check available Struphy data structures for possible types of variables. The parameters of this Propagator are the two floats c1 and c2, as well as some solver parameters whose default is taken from the options attribute discussed below.

Accessing Struphy objects#

Within a propagator, and within the __init__(self, ...) constructor in particular, you have access to several features like spline spaces, geometry, mass matrices, Hodge operators and projected MHD equilibria:

self.derham returns the simulation instance of Derham
self.domain returns the simulation instance of Domain
self.mass_ops returns the simulation instance of WeightedMassOperators
self.basis_ops returns the simulation instance of BasisProjectionOperators
self.projected_mhd_equil returns the simulation instance of ProjectedMHDequilibrium

You can check the Struphy API for more details on how to use these objects. As an example, here is a small snippet of struphy.propagators.propagators_fields.Maxwell:

# Define block matrix [[A B], [C I]] (without time step size dt in the diagonals)
_A = self.mass_ops.M1

self._B = -1/2 * self.derham.curl.T @ self.mass_ops.M2
self._C = 1/2 * self.derham.curl

Particle push and deposition#

Due to the large number \(N \gg 1\) of particles in a PIC simulation, all particle loops in Struphy are accelerated with pyccel, and can be distributed with MPI.

In order to facilitate the integration of particle loops in Propagators, Struphy provides the following “wrapper classes”:

Pusher for solving particle ODEs (explicit, implicit, iterative)
Accumulator for particle deposition into a matrix, and optionally a vector
AccumulatorVector for particle deposition into a vector only

These wrapper classes take care of MPI communication, nonlinear iteration, multi-stage time stepping, among other technical issues. All the developer has to provide is the particle loop in the same way as for a single process, a single iteration, or a single stage of the chosen algorithm.

The single-process/single-stage/single-iteration particle loop must be written into a so-called “kernel file”. The above wrapper classes then take the kernel as an input (plus some other arguments discussed below) and take care of the mentioned abstractions. As an example, let us look at the Propagator struphy.propagators.propagators_markers.PushEta for position update in logical space,

\[\frac{\textnormal d \boldsymbol \eta_p(t)}{\textnormal d t} = DF^{-1}(\boldsymbol \eta_p(t)) \,\mathbf v_p\,,\]

by means of an explicit Runge-Kutta scheme obtained from a ButcherTableau:

def __init__(self,
             particles: Particles,
             *,
             algo: str = options(default=True)['algo'],
             ):

    # base class constructor call
    super().__init__(particles)

    # get kernel
    kernel = pusher_kernels.push_eta_stage

    # define algorithm
    butcher = ButcherTableau(algo)
    args_kernel = (butcher.a,
                   butcher.b,
                   butcher.c)

    # instantiate Pusher
    self._pusher = Pusher(particles,
                          kernel,
                          args_kernel,
                          self.derham.args_derham,
                          self.domain.args_domain,
                          alpha_in_kernel=1.,
                          n_stages=butcher.n_stages,
                          mpi_sort='last')

def __call__(self, dt):
    # push markers
    self._pusher(dt)

Some remarks to this code:

This Propagator updates one variable particles of type Particles
The Propagator base class has to be instantiated with a call to super().__init__()
The kernel (for a single RK stage particle loop) is push_eta_stage from the kernel file pusher_kernels.py
The coefficients of the RK algorithm are retrieved from the ButcherTableau; the algo is a keyword argument of the Propagator
The Pusher is instantiated as self._pusher with its necessary arguments, in particular the particles object and the kernel
The object self._pusher is callable and called in the __call__() method of the Propagator with the time step dt as the single argument.

Particle depositions have a similar structure. As an example, let us look at the Propagator struphy.propagators.propagators_coupling.VlasovAmpere for the charge deposition in the Vlasov-Ampère system:

def __init__(self,
             e: BlockVector,
             particles: Particles6D,
             *,
             c1: float = 1.,
             c2: float = 1.,
             solver=options(default=True)['solver']):

    super().__init__(e, particles)

    # get accumulation kernel
    accum_kernel = accum_kernels.vlasov_maxwell

    # Initialize Accumulator object
    self._accum = Accumulator(particles,
                              'Hcurl',
                              accum_kernel,
                              self.mass_ops,
                              self.domain.args_domain,
                              add_vector=True,
                              symmetry='symm')

    ...

def __call__(self, dt):

    # accumulate
    self._accum()

    ...

Some remarks to this code:

This Propagator updates two variables: e of type BlockVector and particles of type Particles6D
The Propagator base class has to be instantiated with a call to super().__init__()
The kernel is vlasov_maxwell from the kernel file accum_kernels.py
The Accumulator is instantiated as self._accum with its necessary arguments, in particular the particles object, the space 'Hcurl' of the deposition and the accum_kernel
The object self._accum is callable and called in the __call__() method of the Propagator to perform the accumulation.

Particle kernels#

A “kernel” is where the particle loops are written in Struphy. The following kernel files are available:

pic/pushing/pusher_kernels.py for general particle pushing
pic/pushing/pusher_kernels_gc.py for guiding-center pushing
pic/pushing/eval_kernels_gc.py for particle evaluation of specific functions
pic/accumulation/accum_kernels.py for general particle deposition
pic/accumulation/accum_kernels_gc.py for particle deposition in guiding-center models

These kernel files are compiled when the struphy compile command is executed from the console.

In a kernel there are usually some of the following tasks to perform:

evaluate metric coefficients at the particle position \(\boldsymbol \eta_p\)
evaluate FEEC spline fields at the particle position \(\boldsymbol \eta_p\)
compute a sum over nearest-neighbor particles.

Within a kernel the metric coefficients of the map \(F:[0, 1]^3 \to \Omega\) are available through the following module, imported at the top of the kernel files:

import struphy.geometry.evaluation_kernels as evaluation_kernels

This provides callables to all things mapping. Linear algebra operations are available through the module:

import struphy.linear_algebra.linalg_kernels as linalg_kernels

which provides products, transpose, inverse, etc.

The evaluation of FEEC spline fields is managed through the following functions, which are imported at the top of the kernel files as well:

get_spans
eval_0form_spline_mpi
eval_1form_spline_mpi
eval_2form_spline_mpi
eval_3form_spline_mpi
eval_vectorfield_spline_mpi

Here is an example of a Struphy particle loop from the kernel struphy.pic.pushing.pusher_kernels.push_v_with_efield(), for updating the Cartesian velocity

\[\frac{\text{d} \mathbf{v}_p}{\text{d} t} = \kappa \, DF^{-T} \mathbf{E}(\boldsymbol \eta_p) \,,\]

using the following kernel:

for ip in range(n_markers):

    # only do something if particle is a "true" particle (i.e. not a hole)
    if markers[ip, 0] == -1.:
        continue

    eta1 = markers[ip, 0]
    eta2 = markers[ip, 1]
    eta3 = markers[ip, 2]

    # evaluate Jacobian, result in dfm
    evaluation_kernels.df(eta1, eta2, eta3,
                          args_domain,
                          dfm)

    # metric coeffs
    det_df = linalg_kernels.det(dfm)
    linalg_kernels.matrix_inv_with_det(dfm, det_df, dfinv)
    linalg_kernels.transpose(dfinv, dfinvt)

    # spline evaluation
    span1, span2, span3 = get_spans(eta1, eta2, eta3, args_derham)

    # electric field: 1-form components
    eval_1form_spline_mpi(span1, span2, span3,
                          args_derham,
                          e1_1,
                          e1_2,
                          e1_3,
                          e_form)

    # electric field: Cartesian components
    linalg_kernels.matrix_vector(dfinvt, e_form, e_cart)

    # update velocities
    markers[ip, 3:6] += dt * const * e_cart

Some remarks to this code:

For each marker with index ip, the Jacobian of the mapping df is evaluated at the logical marker position eta
Other metric coefficients such as det_df, dfinv and dfinvt are computed from df using linalg_kernels
The spline evaluation happens in two steps:
1. Compute the knot span indices (see below) in the three directions at eta using get_spans
2. Evaluate the 1-form by calling eval_1form_spline_mpi and store the result in e_form
The marker velocities are updated in the last line in the markers array.

Knot span index: given a spline space of degree \(d\), at each particle position \(\boldsymbol \eta_p\) there are \(d+1\) non-zero B-spline basis functions. The \(d+1\) indices of these basis functions are given by \([s(\boldsymbol \eta_p)-d, \ldots, s(\boldsymbol \eta_p)]\), where \(s(\boldsymbol \eta_p) \in \mathbb N\) is the so-called knot span index at position \(\boldsymbol \eta_p\).

Helper classes#

ButcherTableau for choosing a Runge-Kutta method used in Pusher.
SchurSolver for solving 2x2 block systems arising from mid-point rule
SchurSolverFull for solving 2x2 block systems with general rhs

Table of Contents

How to write a Propagator#

Basics#

Accessing Struphy objects#

Particle push and deposition#

Particle kernels#

Helper classes#

This Page