PIC

Particle base class

class struphy.pic.base.Particles(comm_world: Intracomm | None = None, clone_config: CloneConfig | None = None, domain_decomp: tuple | None = None, n_cols_diagnostics: int | None = None, n_cols_aux: int | None = None, type: str = 'full_f', name: str = 'some_name', loading_params: LoadingParameters | None = None, weights_params: WeightsParameters | None = None, boundary_params: BoundaryParameters | None = None, sorting_params: SortingParameters | None = None, saving_params: SavingParameters | None = None, bufsize: float = 0.25, domain: Domain | None = None, equil: FluidEquilibrium | None = None, projected_equil: ProjectedFluidEquilibrium | None = None, background: KineticBackground | FluidEquilibrium | None = None, initial_condition: KineticBackground | None = None, perturbations: dict[str, Perturbation] | None = None, n_as_volume_form: bool = False, equation_params: dict | None = None)

Bases: object

Base class for particle species.

abstract property type

‘full_f’, ‘delta_f’ or ‘sph’.

Type:

Particle type

abstract property vdim

Dimension of the velocity space.

abstract property default_background

The default background (of type Maxwellian).

property default_n_cols

12} for default number of columns.

Type:

Dictionary of the form {‘diagnostics’

Type:

3, ‘aux’

abstract svol(eta1, eta2, eta3, *v)

Marker sampling distribution function \(s^\textrm{vol}\) as a volume form, see Monte-Carlo integrals.

abstract s0(eta1, eta2, eta3, *v, flat_eval=False, remove_holes=True)

Marker sampling distribution function \(s^0\) as 0-form, see Monte-Carlo integrals.

property n_cols_diagnostics

Number of columns for storing diagnostics for each marker.

property n_cols_aux

Number of auxiliary columns for each marker (e.g. for storing evaluation data).

property first_diagnostics_idx

Starting index for diagnostics columns: after 3 positions, vdim velocities, weight, s0 and w0.

property first_pusher_idx

Starting index for storing initial conditions for a Pusher call.

property n_cols_pusher

Dimension of the phase space (for storing initial conditions for a Pusher call).

property first_shift_idx

First index for storing shifts due to boundary conditions in eta-space.

property n_cols_shift

Number of columns for storing shifts due to boundary conditions in eta-space.

property residual_idx

Column for storing the residual in iterative pushers.

property first_free_idx

First index for storing auxiliary quantities for each particle.

property n_cols

Total umber of columns in markers array. The last 2 columns refer to box number and particle ID, respectively.

property n_rows

Total number of rows in markers array.

property kinds

Name of the class.

property name

Name of the kinetic species in DATA container.

property loading: Literal['pseudo_random', 'sobol_standard', 'sobol_antithetic', 'external', 'restart', 'tesselation']

Type of particle loading.

property bc

List of particle boundary conditions in each direction.

property bc_refill

How to re-enter particles if bc is ‘refill’.

property bc_sph

List of boundary conditions for sph evaluation in each direction.

property mean_velocity_index

Index in marker array where mean velocity for noslip BC is stored.

property Np

Total number of markers/particles, from user input.

property Np_per_clone

Array where i-th entry corresponds to the number of loaded particles on clone i. (This is not necessarily the number of valid markers per clone, see self.n_mks_on_each_clone).

property ppc

Particles per cell (=Np if no grid is present).

property ppb

Particles per sorting box.

property bufsize

Relative size of buffer in markers array.

property mpi_comm

MPI communicator.

property mpi_size

Number of MPI processes.

property mpi_rank

Rank of current process.

property clone_config

Manages the configuration for clone-based (copied grids) parallel processing using MPI.

property num_clones

Total number of clones.

property clone_id

Clone id of current process.

property background: KineticBackground

Kinetic background.

property perturbations: dict[str, Perturbation]

Kinetic perturbations, keys are the names of moments of the distribution function (“n”, “u1”, etc.).

property boundary_params: BoundaryParameters

Parameters for marker loading.

property sorting_params: SortingParameters

Parameters for marker sorting.

property reject_weights

Whether to reect weights below threshold.

property threshold

Threshold for rejecting weights.

property boxes_per_dim

Tuple, number of sorting boxes per dimension.

property equation_params

Parameters appearing in model equation due to Struphy normalization.

property initial_condition: KineticBackground

Kinetic initial condition

property f_init

Callable initial condition (background + perturbation). For kinetic models this is a Maxwellian. For SPH models this is a FluidEquilibrium.

property u_init

Callable initial condition (background + perturbation) for the Cartesian velocity in SPH models.

property control_variate

Boolean for whether to use the Control variate method during time stepping.

property domain_array

A 2d array[float] of shape (comm.Get_size(), 9). The row index denotes the process number and for n=0,1,2:

• domain_array[i, 3*n + 0] holds the LEFT domain boundary of process i in direction eta_(n+1).

• domain_array[i, 3*n + 1] holds the RIGHT domain boundary of process i in direction eta_(n+1).

• domain_array[i, 3*n + 2] holds the number of cells of process i in direction eta_(n+1).

property mpi_dims_mask

3-list | tuple; True if the dimension is to be used in the domain decomposition (=default for each dimension). If mpi_dims_mask[i]=False, the i-th dimension will not be decomposed.

property nprocs

Number of MPI processes in each dimension.

property n_mks_load

Array of number of markers on each process at loading stage

property markers

2D numpy array holding the marker information, including holes. The i-th row holds the i-th marker info.

index	0 | 1 | 2 |	3 | … | 3+(vdim-1)|	3+vdim	4+vdim	5+vdim	>=6+vdim	…	-2	-1
value	position (eta)	velocities	weight	s0	w0	other	…	box	ID

The column indices referring to different attributes can be obtained from index.

property holes

Array of booleans stating if an entry in the markers array is a hole.

property ghost_particles

Array of booleans stating if an entry in the markers array is a ghost particle.

property markers_wo_holes

Array holding the marker information, excluding holes. The i-th row holds the i-th marker info.

property markers_wo_holes_and_ghost

Array holding the marker information, excluding holes and ghosts (only valid markers). The i-th row holds the i-th marker info.

property domain

From struphy.geometry.domains.

property equil

From struphy.fields_background.equils.

property projected_equil

MHD equilibrium projected on 3d Derham sequence with commuting projectors.

property lost_markers

Array containing the last infos of removed markers

property n_lost_markers

Number of removed particles.

property index

Dict holding the column indices referring to specific marker parameters (coordinates).

property valid_mks

Array of booleans stating if an entry in the markers array is a true local particle (not a hole or ghost).

property n_mks_loc

Number of valid markers on process (without holes and ghosts).

property n_mks_on_each_proc

Array where i-th entry corresponds to the number of valid markers on i-th process (without holes and ghosts).

property n_mks_on_clone

Number of valid markers on current clone (without holes and ghosts).

property n_mks_on_each_clone

Number of valid markers on current clone (without holes and ghosts).

property n_mks_global

Number of valid markers on current clone (without holes and ghosts).

property positions

Array holding the marker positions in logical space. The i-th row holds the i-th marker info.

property velocities

Array holding the marker velocities in logical space. The i-th row holds the i-th marker info.

property phasespace_coords

Array holding the marker positions and velocities in logical space. The i-th row holds the i-th marker info.

property weights

Array holding the current marker weights. The i-th row holds the i-th marker info.

property weights_at_t0

Array holding the initial marker weights. The i-th row holds the i-th marker info.

property sampling_density

Array holding the current marker 0form sampling density s0. The i-th row holds the i-th marker info.

property weights0

Array holding the initial marker weights. The i-th row holds the i-th marker info.

property marker_ids

Array holding the marker id’s on the current process.

property is_volume_form

True means volume-form, False means 0-form.

Type:

Tuple of size 2 for (position, velocity), defining the p-form representation of f_init

property spatial

Drawing particles uniformly on the unit cube(‘uniform’) or on the disc(‘disc’)

property f_coords_index

Dict holding the column indices referring to coords of the distribution fuction.

property f_jacobian_coords_index

Dict holding the column indices referring to coords of the velocity jacobian determinant of the distribution fuction.

property f_coords

Coordinates of the distribution function.

property args_markers

Collection of mandatory arguments for pusher kernels.

property f_jacobian_coords

Coordinates of the velocity jacobian determinant of the distribution fuction.

property tesselation

Tesselation of the current process domain.

classmethod ker_dct()

Available smoothing kernels, numbers must be multiplies of 100.

draw_markers(sort: bool = True)

Drawing markers

• for PIC: according to the volume density \(s^\textrm{vol}_{\textnormal{in}}\)

• for SPH: from unity/disc in space and according to the vector-field representation of the fluid velocity

In Struphy, the initial marker distribution \(s^\textrm{vol}_{\textnormal{in}}\) is always of the form

\[s^\textrm{vol}_{\textnormal{in}}(\eta,v) = n^3(\eta)\, \mathcal M(v)\,,\]

with \(\mathcal M(v)\) a multi-variate Gaussian:

\[\mathcal M(v) = \prod_{i=1}^{d_v} \frac{1}{\sqrt{2\pi}\,v_{\mathrm{th},i}} \exp\left[-\frac{(v_i-u_i)^2}{2 v_{\mathrm{th},i}^2}\right]\,,\]

where \(d_v\) stands for the dimension in velocity space, \(u_i\) are velocity constant shifts and \(v_{\mathrm{th},i}\) are constant thermal velocities (standard deviations). The function \(n^3:(0,1)^3 \to \mathbb R^+\) is a normalized 3-form on the unit cube,

\[\int_{(0,1)^3} n^3(\eta)\,\textnormal d \eta = 1\,.\]

The following choices are available in Struphy:

1. Uniform distribution on the unit cube: \(n^3(\eta) = 1\)

2. Uniform distribution on the disc: \(n^3(\eta) = 2\eta_1\) (radial coordinate = volume element of square-to-disc mapping)

Velocities are sampled via inverse transform sampling. In case of Particles6D, velocities are sampled as a Maxwellian in each 3 directions,

\[r_i = \int^{v_i}_{-\infty} \mathcal M(v^\prime_i) \textnormal{d} v^\prime_i = \frac{1}{2}\left[ 1 + \text{erf}\left(\frac{v_i - u_i}{\sqrt{2}v_{\mathrm{th},i}}\right)\right] \,,\]

where \(r_i \in \mathcal R(0,1)\) is a uniformly drawn random number in the unit interval. So then

\[v_i = \text{erfinv}(2r_i - 1)\sqrt{2}v_{\mathrm{th},i} + u_i \,.\]

In case of Particles5D, parallel velocity is sampled as a Maxwellian and perpendicular particle speed \(v_\perp = \sqrt{v_1^2 + v_2^2}\) is sampled as a 2D Maxwellian in polar coordinates,

\[\begin{split}\mathcal{M}(v_1, v_2) \, \textnormal{d} v_1 \textnormal{d} v_2 &= \prod_{i=1}^{2} \frac{1}{\sqrt{2\pi}}\frac{1}{v_{\mathrm{th},i}} \exp\left[-\frac{(v_i-u_i)^2}{2 v_{\mathrm{th},i}^2}\right] \textnormal{d} v_i\,, \\ &= \frac{1}{v_\mathrm{th}^2}v_\perp \exp\left[-\frac{(v_\perp-u)^2}{2 v_\mathrm{th}^2}\right] \textnormal{d} v_\perp\,, \\ &= \mathcal{M}^{\textnormal{pol}}(v_\perp) \, \textnormal{d} v_\perp \,.\end{split}\]

Then,

\[r = \int^{v_\perp}_0 \mathcal{M}^{\textnormal{pol}} \textnormal{d} v_\perp = 1 - \exp\left[-\frac{(v_\perp-u)^2}{2 v_\mathrm{th}^2}\right] \,.\]

So then,

\[v_\perp = \sqrt{- \ln(1-r)}\sqrt{2}v_\mathrm{th} + u \,.\]

All needed parameters can be set in the parameter file, see params_yml.

An initial sorting will be performed if sort is given as True (default) and sorting_params were given to the init.

Parameters:

sort (Bool) – Wether to sort the particules in boxes after initial drawing (only if sorting params were passed)

mpi_sort_markers(apply_bc: bool = True, alpha: tuple | list | int | float = 1.0, do_test: bool = False, remove_ghost: bool = True)

Sorts markers according to MPI domain decomposition.

Markers are sent to the process corresponding to the alpha-weighted position alpha*markers[:, 0:3] + (1 - alpha)*markers[:, first_pusher_idx:first_pusher_idx + 3].

Periodic boundary conditions are taken into account when computing the alpha-weighted position.

Parameters:

• appl_bc (bool) – Whether to apply kinetic boundary conditions before sorting.

• alpha (tuple | list | int | float) – For i=1,2,3 the sorting is according to alpha[i]*markers[:, i] + (1 - alpha[i])*markers[:, first_pusher_idx + i]. If int or float then alpha = (alpha, alpha, alpha). alpha must be between 0 and 1.

• do_test (bool) – Check if all markers are on the right process after sorting.

• remove_ghost (bool) – Remove ghost particles before send.

initialize_weights(*, bckgr_params: dict | None = None, pert_params: dict | None = None)

Computes the initial weights

\[w_{k0} := \frac{f^0(t, q_k(t)) }{s^0(t, q_k(t)) } = \frac{f^0(0, q_k(0)) }{s^0(0, q_k(0)) } = \frac{f^0_{\textnormal{in}}(q_{k0}) }{s^0_{\textnormal{in}}(q_{k0}) }\]

from the initial distribution function \(f^0_{\textnormal{in}}\) specified in the parmeter file and from the initial volume density \(s^n_{\textnormal{vol}}\) specified in draw_markers(). Moreover, it sets the corresponding columns for “w0”, “s0” and “weights” in the markers array. If control_variate is True, the background f0 is subtracted.

Parameters:

• bckgr_params (dict) – Kinetic background parameters.

• pert_params (dict) – Kinetic perturbation parameters for initial condition.

update_weights()

Applies the control variate method, i.e. updates the time-dependent marker weights according to the algorithm in Control variate method. The background f0 is used for this.

reset_marker_ids()

Reset the marker ids (last column in marker array) according to the current distribution of particles. The first marker on rank 0 gets the id ‘0’, the last marker on the last rank gets the id ‘n_mks_global - 1’.

binning(components: tuple[bool], bin_edges: tuple[ndarray], output_quantity: Literal['density', 'current_1', 'current_2', 'current_3', 'energy_tensor_11', 'energy_tensor_22', 'energy_tensor_33', 'energy_tensor_12', 'energy_tensor_13', 'energy_tensor_23', 'heat_flux_1', 'heat_flux_2', 'heat_flux_3'] = 'density', divide_by_jac: bool = True)

Computes full-f and delta-f distribution functions via marker binning in logical space. Numpy’s histogramdd is used, following the algorithm outlined in Particle binning.

Parameters:

• components (tuple[bool]) – List of length 3 + vdim; an entry is True if the direction in phase space is to be binned.

• bin_edges (tuple[array]) – List of bin edges (resolution) having the length of True entries in components.

• output_quantity (BinningOutput) – String literal used to determine weights in binning and the type of output

• divide_by_jac (bool) – Whether to divide the weights by the Jacobian determinant for binning.

Returns:

• f_slice (array-like) – The reconstructed full-f distribution function.

• df_slice (array-like) – The reconstructed delta-f distribution function.

show_distribution_function(components, bin_edges)

1D and 2D plots of slices of the distribution function via marker binning. This routine is mainly for de-bugging.

Parameters:

• components (list[bool]) – List of length 6 giving the directions in phase space in which to bin.

• bin_edges (list[array]) – List of bin edges (resolution) having the length of True entries in components.

apply_kinetic_bc(newton=False)

Apply boundary conditions to markers that are outside of the logical unit cube.

Parameters:

newton (bool) – Whether the shift due to boundary conditions should be computed for a Newton step or for a strandard (explicit or Picard) step.

particle_refilling()

When particles move outside of the domain, refills them. TODO: Currently only valid for HollowTorus geometry with AdhocTorus equilibrium.

In case of guiding-center orbit, refills particles at the opposite poloidal angle of the same magnetic flux surface.

\[\begin{split}\theta_\text{refill} &= - \theta_\text{loss} \\ \phi_\text{refill} &= -2 q(r_\text{loss}) \theta_\text{loss}\end{split}\]

In case of full orbit, refills particles at the same gyro orbit until their guiding-centers are also outside of the domain. When their guiding-centers also reach at the boundary, refills them as we did with guiding-center orbit.

gyro_transfer(outside_inds)

Refills particles at the same gyro orbit. Their perpendicular velocity directions are also changed accordingly:

First, refills the particles at the other side of the cross point (between gyro circle and domain boundary),

\[\theta_\text{refill} = \theta_\text{gc} - \left(\theta_\text{loss} - \theta_\text{gc} \right) \,.\]

Then changes the direction of the perpendicular velocity,

\[\vec{v}_{\perp, \text{refill}} = \frac{\vec{\rho}_g}{|\vec{\rho}_g|} \times \vec{b}_0 |\vec{v}_{\perp, \text{loss}}| \,,\]

where \(\vec{\rho}_g = \vec{x}_\text{refill} - \vec{X}_\text{gc}\) is the cartesian radial vector.

Parameters:

outside_inds (xp.array (int)) – An array of indices of particles which are outside of the domain.

Returns:

out – An array of indices of particles where its guiding centers are outside of the domain.

Return type:

xp.array (bool)

class SortingBoxes(markers_shape: tuple, is_sph: bool, *, nx: int = 1, ny: int = 1, nz: int = 1, bc_sph: list | None = None, is_domain_boundary: dict | None = None, comm: Intracomm | None = None, box_index: int = -2, box_bufsize: float = 2.0)

Bases: object

Boxes used for the sorting of the particles.

Boxes are represented as a 2D array of integers, where each line coresponds to one box, and all entries of line i that are not -1 correspond to a particles in the i-th box.

Parameters:

• markers_shape (tuple) – shape of 2D marker array.

• is_sph (bool) – True if particle type is “sph”.

• nx (int) – number of boxes in the x direction.

• ny (int) – number of boxes in the y direction.

• nz (int) – number of boxes in the z direction.

• bc_sph (list) – Boundary condition for sph density evaluation. Either ‘periodic’, ‘mirror’, ‘fixed’ or ‘noslip’ in each direction.

• is_domain_boundary (dict) – Has two booleans for each direction; True when the boundary of the MPI process is a domain boundary.

• comm (Intracomm) – MPI communicator or None.

• box_index (int) – Column index of the particles array to store the box number, counted from the end (e.g. -2 for the second-to-last).

• box_bufsize (float) – additional buffer space in the size of the boxes

property is_domain_boundary

Dict with two booleans for each direction (e.g. ‘x_m’ and ‘x_p’); True when the boundary of the MPI process is a domain boundary (0.0 or 1.0).

property bc_sph

List of boundary conditions for sph evaluation in each direction.

property bc_sph_index_shifts

Dictionary holding the index shifts of box number for ghost particles in each direction.

__weakref__

list of weak references to the object (if defined)

put_particles_in_boxes()

Assign the right box to the particles and the list of the particles to each box. If sorting_boxes was instantiated with an MPI comm, then the particles in the neighbouring boxes of neighbours processors or also communicated

check_and_assign_particles_to_boxes()

Check whether the box array has enough columns (detect load imbalance wrt to sorting boxes), and then assigne the particles to boxes.

do_sort(use_numpy_argsort=False)

Assign the particles to boxes and then sort them.

prepare_ghost_particles()

Markers for boundary conditions and MPI communication.

Does the following: 1. determine which markers belong to boxes that are at the boundary and put these markers in a new array (e.g. markers_x_m) 2. set their last index to -2 to indicate that they will be “ghost particles” after sending 3. set their new box number (boundary conditions enter here) 4. optional: mirror position for boundary conditions

determine_markers_in_box(list_boxes)

Determine the markers that belong to a certain box (list of boxes) and put them in an array

get_destinations_box()

Find the destination proc for the particles to communicate for the box structure.

self_communication_boxes()

Communicate the particles in case a process is it’s own neighbour (in case of periodicity with low number of procs/boxes)

sendrecv_all_to_all_boxes()

Distribute info on how many markers will be sent/received to/from each process via all-to-all for the communication of particles in boundary boxes.

sendrecv_markers_boxes()

Use non-blocking communication. In-place modification of markers for the communication of particles in boundary boxes.

eval_density(eta1, eta2, eta3, h1, h2, h3, kernel_type='gaussian_1d', derivative=0, fast=True)

Evaluate particle number density (0-form) using an SPH smoothing kernel.

Parameters:

• eta1 (array_like) – Logical evaluation points. Inputs may be 1-D arrays (flat evaluation) or broadcastable meshgrid arrays; the output will match the shape of eta1.

• eta2 (array_like) – Logical evaluation points. Inputs may be 1-D arrays (flat evaluation) or broadcastable meshgrid arrays; the output will match the shape of eta1.

• eta3 (array_like) – Logical evaluation points. Inputs may be 1-D arrays (flat evaluation) or broadcastable meshgrid arrays; the output will match the shape of eta1.

• h1 (float) – Support radius of the smoothing kernel in each logical dimension.

• h2 (float) – Support radius of the smoothing kernel in each logical dimension.

• h3 (float) – Support radius of the smoothing kernel in each logical dimension.

• kernel_type (str, optional) – Name of the smoothing kernel (must be a key in self.ker_dct()).

• derivative (int, optional) – Selects whether to evaluate the kernel derivative along a coordinate direction: 0 (default) returns the scalar density, 1/2/3 returns the corresponding component of the density gradient with respect to logical coordinates.

• fast (bool, optional) – If True, use the box-based neighbor search (faster for many particles); if False, use the naive all-pairs evaluation (simpler, slower).

Returns:

out – Estimated number density (or requested derivative component) at the provided evaluation points. The array uses the same shape as eta1 and is returned as a cunumpy (xp) array.

Return type:

xp.ndarray

Notes

This method is a thin wrapper around eval_sph() and internally evaluates the column given by self.index[‘weights’] (particle weights).

eval_velocity(eta1, eta2, eta3, h1, h2, h3, kernel_type='gaussian_1d', derivative=0, fast=True) → tuple

Estimate mean velocity components using SPH smoothing.

Parameters:

• eta1 (array_like) – Logical evaluation points. May be 1-D arrays or broadcastable meshgrid arrays; the returned component arrays match the shape of eta1.

• eta2 (array_like) – Logical evaluation points. May be 1-D arrays or broadcastable meshgrid arrays; the returned component arrays match the shape of eta1.

• eta3 (array_like) – Logical evaluation points. May be 1-D arrays or broadcastable meshgrid arrays; the returned component arrays match the shape of eta1.

• h1 (float) – Support radius of the smoothing kernel in each logical dimension.

• h2 (float) – Support radius of the smoothing kernel in each logical dimension.

• h3 (float) – Support radius of the smoothing kernel in each logical dimension.

• kernel_type (str, optional) – Name of the smoothing kernel (must be a key in self.ker_dct()).

• derivative (int, optional) – If 0 (default) evaluate the mean velocity; if 1/2/3 return the corresponding component of the spatial derivative of the velocity.

• fast (bool, optional) – If True use the box-based neighbor search (faster for many particles); if False use the naive all-pairs evaluation.

Returns:

(v1, v2, v3) – Three arrays containing the estimated velocity components at the provided evaluation points. Each array has the same shape as eta1.

Return type:

tuple of xp.ndarray

Notes

This method first computes SPH coefficients by calling eval_kernels_sph.sph_mean_velocity_coeffs (via a Pyccel kernel) to assemble mean-velocity coefficients into the markers array, then calls eval_sph() for each velocity component.

eval_div_viscosity(eta1, eta2, eta3, h1, h2, h3, kernel_type='gaussian_1d', mu: float = 1.0, fast=True) → tuple

Compute divergence of the viscous stress (mu * viscosity tensor).

Parameters:

• eta1 (array_like) – Logical evaluation points where the divergence is evaluated.

• eta2 (array_like) – Logical evaluation points where the divergence is evaluated.

• eta3 (array_like) – Logical evaluation points where the divergence is evaluated.

• h1 (float) – Support radius of the smoothing kernel in each logical dimension.

• h2 (float) – Support radius of the smoothing kernel in each logical dimension.

• h3 (float) – Support radius of the smoothing kernel in each logical dimension.

• kernel_type (str, optional) – Name of the smoothing kernel (must be a key in self.ker_dct()).

• mu (float, optional) – Dynamic viscosity coefficient used in the viscosity kernel.

• fast (bool, optional) – If True use the box-based neighbor search; if False use naive evaluation.

Returns:

(gamma_x, gamma_y, gamma_z) – Components of the divergence of the viscous stress evaluated at the provided points. Each array matches the shape of eta1.

Return type:

tuple of xp.ndarray

Notes

The routine populates intermediate marker columns using two Pyccel kernels: sph_mean_velocity_coeffs (mean velocity) and sph_viscosity_tensor (viscosity tensor components). It then evaluates the necessary derivatives via eval_sph() and sums contributions to produce the three divergence components.

eval_sph(eta1: ndarray, eta2: ndarray, eta3: ndarray, index: int, out: ndarray | None = None, fast: bool = True, kernel_type: str = 'gaussian_1d', derivative: int = 0, h1: float = 0.1, h2: float = 0.1, h3: float = 0.1)

Perform a (meshgrid) SPH evaluation of a function \(\rho: [0, 1]^3 \to \mathbb R\) in the following sense:

\[\rho(\boldsymbol \eta_i) = \sum_{j=0}^{N-1} \rho_j\, W_h(\boldsymbol \eta_i - \boldsymbol \eta_j)\,.\]

The coefficients \(\rho_j\) must be available in the marker array, stored at some index self.markers[j, index]. In case that derivative=k where k is not zero, the k-th component of the gradient of \(\rho\) is computed:

\[\textrm{derivative}=k:\qquad [\nabla \rho(\boldsymbol \eta_i)]_k = \sum_{j=0}^{N-1} \rho_j \frac{\partial W_h}{\partial \eta_k}(\boldsymbol \eta_i - \boldsymbol \eta_j)\,.\]

The possible choices for \(W_h\) are listed in Smoothing kernels and in ker_dct().

Parameters:

• eta1 (array_like) – Logical evaluation points.

• eta2 (array_like) – Logical evaluation points.

• eta3 (array_like) – Logical evaluation points.

• index (int) – At which index of the markers array are located the coefficients \(\rho_j\).

• out (array_like) – Output will be store in this array. A new array is created if not provided.

• fast (bool) – If true, uses an optimized evaluation algorithm taking advantage of the box structure. This assume that the boxes are bigger then the radius used for the smoothing kernel.

• kernel_type (str) – Name of the smoothing kernel, see sph_smoothing_kernels and ker_dct().

• derivative (int) – 0: no kernel derivative 1: first component of grad 2: second component of grad 3: third component of grad

• h1 (float) – Radius of the smoothing kernel in each dimension.

• h2 (float) – Radius of the smoothing kernel in each dimension.

• h3 (float) – Radius of the smoothing kernel in each dimension.

update_holes()

Compute new holes, new number of holes and markers on process

update_ghost_particles()

Compute new particles that belong to boundary processes needed for sph evaluation

sendrecv_determine_mtbs(alpha: list | tuple | ndarray = (1.0, 1.0, 1.0))

Determine which markers have to be sent from current process and put them in a new array. Corresponding rows in markers array become holes and are therefore set to -1. This can be done purely with numpy functions (fast, vectorized).

Parameters:

alpha (list | tuple) – For i=1,2,3 the sorting is according to alpha[i]*markers[:, i] + (1 - alpha[i])*markers[:, first_pusher_idx + i]. alpha[i] must be between 0 and 1.

Returns:

• hole_inds_after_send (array[int]) – Indices of empty columns in markers after send.

• sorting_etas (array[float]) – Eta-values of shape (n_send, :) according to which the sorting is performed.

sendrecv_get_destinations(send_inds)

Determine to which process particles have to be sent.

Parameters:

send_inds (array[int]) – Indices of particles which will be sent.

Returns:

send_info – Amount of particles sent to i-th process.

Return type:

array[int]

sendrecv_all_to_all(send_info)

Distribute info on how many markers will be sent/received to/from each process via all-to-all.

Parameters:

send_info (array[int]) – Amount of markers to be sent to i-th process.

Returns:

recv_info – Amount of marticles to be received from i-th process.

Return type:

array[int]

sendrecv_markers(recv_info, hole_inds_after_send)

Use non-blocking communication. In-place modification of markers

Parameters:

• recv_info (array[int]) – Amount of markers to be received from i-th process.

• hole_inds_after_send (array[int]) – Indices of empty rows in markers after send.

__weakref__

list of weak references to the object (if defined)

Particel-to-grid accumulation

Base classes for particle deposition (accumulation) on the grid.

class struphy.pic.accumulation.particles_to_grid.Accumulator(particles: Particles, space_id: str, kernel: Pyccelkernel, mass_ops: WeightedMassOperators, args_domain: DomainArguments, *, add_vector: bool = False, symmetry: str | None = None, filter_params: FilterParameters | None = None)

Bases: object

Approximates integrals of the form

\[\begin{split}I_A &= \int_\Omega \int_{\mathbb R^3} \Lambda^\mu_{ijk}(\boldsymbol \eta) \, A^{\mu, \nu}(\boldsymbol \eta, \mathbf v) \, \Lambda^\nu_{mno}(\boldsymbol \eta) \, f^{\textrm{vol}}(\boldsymbol \eta, \mathbf v)\,\mathrm d\mathbf v \textrm d \boldsymbol \eta\,, \\[2mm] I_B &= \int_\Omega \int_{\mathbb R^3} \Lambda^\mu_{ijk}(\boldsymbol \eta) \, B^\mu(\boldsymbol \eta, \mathbf v) \, f^{\textrm{vol}}(\boldsymbol \eta, \mathbf v)\,\mathrm d\mathbf v \textrm d \boldsymbol \eta\,,\end{split}\]

for given weight functions \(A^{\mu,\nu}\) and \(B^\mu\) by Monte-Carlo quadrature through the particle distribution function \(f^{\textrm{vol}}\):

\[f^{\textrm{vol}}(\boldsymbol \eta, \mathbf v) \approx \sum_{p=0}^{N-1} w_p \, \delta(\boldsymbol \eta - \boldsymbol \eta_p) \, \delta(\mathbf v - \mathbf v_p)\,.\]

This results in stencil (block) matrices and vectors

\[\begin{split}M &= (M^{\mu,\nu})_{\mu,\nu}\,,\qquad && M^{\mu,\nu} \in \mathbb R^{\mathbb N^\alpha_\mu \times \mathbb N^\alpha_\nu}\,, \\[2mm] V &= (V^\mu)_\mu\,,\qquad &&V^\mu \in \mathbb R^{\mathbb N^\alpha_\mu}\,,\end{split}\]

where \(N^\alpha_\mu\) denotes the dimension of the \(\mu\)-th component of the Derham space \(V_h^\alpha\) (\(\mu,\nu = 1,2,3\) for vector-valued spaces), with entries obtained by summing over all particles \(p\),

\[\begin{split}M^{\mu,\nu}_{ijk,mno} &= \sum_{p=0}^{N-1} w_p\, \Lambda^\mu_{ijk}(\boldsymbol \eta_p) \, A^{\mu,\nu}_p \, \Lambda^\nu_{mno}(\boldsymbol \eta_p) \,, \\[2mm] V^\mu_{ijk} &= \sum_{p=0}^{N-1} w_p\, \Lambda^\mu_{ijk}(\boldsymbol \eta_p) \, B^\mu_p \,.\end{split}\]

Here, \(\Lambda^\mu_{ijk}(\boldsymbol \eta_p)\) denotes the \(ijk\)-th basis function of the \(\mu\)-th component of a Derham space.

Parameters:

• particles (Particles) – Particles object holding the markers to accumulate.

• space_id (str) – Space identifier for the matrix/vector (H1, Hcurl, Hdiv, L2 or H1vec) to be accumulated into.

• kernel (pyccelized function) – The accumulation kernel.

• derham (Derham) – Discrete FE spaces object.

• args_domain (DomainArguments) – Mapping infos.

• add_vector (bool) – True if, additionally to a matrix, a vector in the same space is to be accumulated. Default=False.

• symmetry (str) – In case of space_id=Hcurl/Hdiv, the symmetry property of the block matrix: diag, asym, symm, pressure or None (=full matrix, default)

• filter_params (dict) – Params for the accumulation filter: use_filter(string, either `three_point or `fourier), repeat(int), alpha(float) and modes(list with int).

Note

Struphy accumulation kernels called by Accumulator objects must be added to struphy/pic/accumulation/accum_kernels.py (6D particles) or struphy/pic/accumulation/accum_kernels_gc.py (5D particles), see accum_kernels and accum_kernels_gc for details.

__init__(particles: Particles, space_id: str, kernel: Pyccelkernel, mass_ops: WeightedMassOperators, args_domain: DomainArguments, *, add_vector: bool = False, symmetry: str | None = None, filter_params: FilterParameters | None = None)

__call__(*optional_args, **args_control)

Performs the accumulation into the matrix/vector by calling the chosen accumulation kernel and additional analytical contributions (control variate, optional).

Parameters:

• particles (Particles) – Particles object holding the markers information in format particles.markers.shape == (n_markers, :).

• optional_args (any) – Additional arguments to be passed to the accumulator kernel, besides the mandatory arguments which are prepared automatically (spline bases info, mapping info, data arrays). Examples would be parameters for a background kinetic distribution or spline coefficients of a background magnetic field. Entries must be pyccel-conform types.

• args_control (any) – Keyword arguments for an analytical control variate correction in the accumulation step. Possible keywords are ‘control_vec’ for a vector correction or ‘control_mat’ for a matrix correction. Values are a 1d (vector) or 2d (matrix) list with callables or xp.ndarrays used for the correction.

property particles

Particle object.

property kernel: Pyccelkernel

The accumulation kernel.

property derham

Discrete Derham complex on the logical unit cube.

property args_domain

Mapping info for evaluating metric coefficients.

property space_id

Space identifier for the matrix/vector (H1, Hcurl, Hdiv, L2 or H1vec) to be accumulated into.

property form

p-form (“0”, “1”, “2”, “3” or “v”) to be accumulated into.

property symmetry

Symmetry of the accumulation matrix (diagonal, symmetric, asymmetric, etc.).

property operators

List of WeightedMassOperators of the accumulator. Matrices can be accessed e.g. with operators[0].matrix.

property vectors

List of Stencil-/Block-/PolarVectors of the accumulator.

property accfilter

Callable filters

init_control_variate(mass_ops)

Set up the use of noise reduction by control variate.

show_accumulated_spline_field(mass_ops: WeightedMassOperators, eta_direction=0, component=0)

1D plot of the spline field corresponding to the accumulated vector. The latter can be viewed as the rhs of an L2-projection:

\[\mathbb M \mathbf a = \sum_p \boldsymbol \Lambda(\boldsymbol \eta_p) * B_p\,.\]

The FE coefficients \(\mathbf a\) determine a FE SplineFunction.

__weakref__

list of weak references to the object (if defined)

class struphy.pic.accumulation.particles_to_grid.AccumulatorVector(particles: Particles, space_id: str, kernel: Pyccelkernel, mass_ops: WeightedMassOperators, args_domain: DomainArguments, filter_params: FilterParameters | None = None)

Bases: object

Approximates integrals of the form

\[I_B = \int_\Omega \int_{\mathbb R^3} \Lambda^\mu_{ijk}(\boldsymbol \eta) \, B^\mu(\boldsymbol \eta, \mathbf v) \, f^{\textrm{vol}}(\boldsymbol \eta, \mathbf v)\,\mathrm d\mathbf v \textrm d \boldsymbol \eta\,,\]

for a given weight function and \(B^\mu\) by Monte-Carlo quadrature through the particle distribution function \(f^{\textrm{vol}}\):

\[f^{\textrm{vol}}(\boldsymbol \eta, \mathbf v) \approx \sum_{p=0}^{N-1} w_p \, \delta(\boldsymbol \eta - \boldsymbol \eta_p) \, \delta(\mathbf v - \mathbf v_p)\,.\]

This results in a stencil (block) vector

\[V = (V^\mu)_\mu\,,\qquad V^\mu \in \mathbb R^{\mathbb N^\alpha_\mu}\,,\]

where \(N^\alpha_\mu\) denotes the dimension of the \(\mu\)-th component of the Derham space \(V_h^\alpha\) (\(\mu,\nu = 1,2,3\) for vector-valued spaces), with entries obtained by summing over all particles \(p\),

\[V^\mu_{ijk} = \sum_{p=0}^{N-1} w_p\, \Lambda^\mu_{ijk}(\boldsymbol \eta_p) \, B^\mu_p \,.\]

Here, \(\Lambda^\mu_{ijk}(\boldsymbol \eta_p)\) denotes the \(ijk\)-th basis function of the \(\mu\)-th component of a Derham space.

Similar to Accumulator but only for vectors \(V\).

Parameters:

• particles (Particles) – Particles object holding the markers to accumulate.

• space_id (str) – Space identifier for the matrix/vector (H1, Hcurl, Hdiv, L2 or H1vec) to be accumulated into.

• kernel (pyccelized function) – The accumulation kernel.

• derham (Derham) – Discrete FE spaces object.

• args_domain (DomainArguments) – Mapping infos.

__init__(particles: Particles, space_id: str, kernel: Pyccelkernel, mass_ops: WeightedMassOperators, args_domain: DomainArguments, filter_params: FilterParameters | None = None)

__call__(*optional_args, **args_control)

Performs the accumulation into the vector by calling the chosen accumulation kernel and additional analytical contributions (control variate, optional).

Parameters:

• optional_args (any) – Additional arguments to be passed to the accumulator kernel, besides the mandatory arguments which are prepared automatically (spline bases info, mapping info, data arrays). Examples would be parameters for a background kinetic distribution or spline coefficients of a background magnetic field. Entries must be pyccel-conform types.

• args_control (any) – Keyword arguments for an analytical control variate correction in the accumulation step. Possible keywords are ‘control_vec’ for a vector correction or ‘control_mat’ for a matrix correction. Values are a 1d (vector) or 2d (matrix) list with callables or xp.ndarrays used for the correction.

property particles

Particle object.

property kernel: Pyccelkernel

The accumulation kernel.

property derham

Discrete Derham complex on the logical unit cube.

property args_domain

Mapping arguments.

property space_id

Space identifier for the matrix/vector (H1, Hcurl, Hdiv, L2 or H1vec) to be accumulated into.

property form

p-form (“0”, “1”, “2”, “3” or “v”) to be accumulated into.

property vectors

List of Stencil-/Block-/PolarVectors of the accumulator.

property accfilter

Callable filters

init_control_variate(mass_ops)

Set up the use of noise reduction by control variate.

show_accumulated_spline_field(mass_ops, eta_direction=(True, False, False))

1 or 2D plot of the spline field corresponding to the accumulated vector. The latter can be viewed as the rhs of an L2-projection:

\[\mathbb M \mathbf a = \sum_p \boldsymbol \Lambda(\boldsymbol \eta_p) * B_p\,.\]

The FE coefficients \(\mathbf a\) determine a FE SplineFunction.

Parameters:

eta_direction – axes of eta to show accumulation (eta1, eta2, eta3).

__weakref__

list of weak references to the object (if defined)

Accumulation kernels

struphy.pic.accumulation.accum_kernels.cc_lin_mhd_6d_1()

Accumulates into V1 with the filling functions

\[A_p^{\mu, \nu} = w_p * [ G^{-1}(\eta_p) * B2_{\times}(\eta_p) * G^{-1}(\eta_p) ]_{\mu, \nu}\]

where \(B2_{\times} * a := B2 \times a\) for \(a \in \mathbb R^3\).

Parameters:

• b2_1 (array[float]) – FE coefficients c_ijk of the magnetic field as a 2-form.

• b2_2 (array[float]) – FE coefficients c_ijk of the magnetic field as a 2-form.

• b2_3 (array[float]) – FE coefficients c_ijk of the magnetic field as a 2-form.

Note

The above parameter list contains only the model specific input arguments.

struphy.pic.accumulation.accum_kernels.cc_lin_mhd_6d_2()

Accumulates into V1 with the filling functions

\[ \begin{align}\begin{aligned}A_p^{\mu, \nu} &= w_p * [ G^{-1}(\eta_p) * B2_{\times}(\eta_p) * G^{-1}(\eta_p) * B2_{\times}(\eta_p)^\top * G^{-1}(\eta_p) ]_{\mu, \nu}\\B_p^\mu &= w_p * [ G^{-1}(\eta_p) * B2_{\times}(\eta_p) * DF^{-1}(\eta_p) * v_p ]_\mu\end{aligned}\end{align} \]

where \(B2_{\times} * a := B2 \times a\) for \(a \in \mathbb R^3\).

Parameters:

• b2_1 (array[float]) – FE coefficients c_ijk of the magnetic field as a 2-form.

• b2_2 (array[float]) – FE coefficients c_ijk of the magnetic field as a 2-form.

• b2_3 (array[float]) – FE coefficients c_ijk of the magnetic field as a 2-form.

Note

The above parameter list contains only the model specific input arguments.

struphy.pic.accumulation.accum_kernels.charge_density_0form()

Kernel for AccumulatorVector into V0 with weight \(B^\mu = 1\).

struphy.pic.accumulation.accum_kernels.hybrid_fA_Arelated()

Accumulates into V1 with the filling functions

\[ \begin{align}\begin{aligned}A_p^{\mu, \nu} &= f_0(\eta_p, v_p) * [ DF^{-1}(\eta_p) * v_p ]_\mu * [ DF^{-1}(\eta_p) * v_p ]_\nu\\B_p^\mu &= \sqrt{f_0(\eta_p, v_p)} * w_p * [ DF^{-1}(\eta_p) * v_p ]_\mu\end{aligned}\end{align} \]

Note

The above parameter list contains only the model specific input arguments.

struphy.pic.accumulation.accum_kernels.hybrid_fA_density()

Accumulates the values of density at quadrature points with the filling functions

\[n = \sum_p w_p S(x - x_p)\]

Parameters:

do (To)

Note

The above parameter list contains only the model specific input arguments.

struphy.pic.accumulation.accum_kernels.linear_vlasov_ampere()

Accumulates into V1 with the filling functions

\[ \begin{align}\begin{aligned}A_p^{\mu, \nu} &= \frac{\alpha^2 \kappa^2}{v_{\text{th}}^2} \frac{1}{N\, s_0} f_0(\mathbf{\eta}_p, \mathbf{v}_p) [ DF^{-1}(\mathbf{\eta}_p) \mathbf{v}_p ]_\mu [ DF^{-1}(\mathbf{\eta}_p) \mathbf{v}_p ]_\nu \,,\\B_p^\mu &= \alpha^2 \kappa \sqrt{f_0(\mathbf{\eta}_p, \mathbf{v}_p)} w_p [ DF^{-1}(\mathbf{\eta}_p) \mathbf{v}_p ]_\mu \,.\end{aligned}\end{align} \]

Parameters:

array[float] (f0_values ;) – Value of f0 for each particle.

Note

The above parameter list contains only the model specific input arguments.

struphy.pic.accumulation.accum_kernels.pc_lin_mhd_6d()

Accumulates into V1 with the filling functions

\[ \begin{align}\begin{aligned}{V_{p,i}}_\perp A_p^{\mu, \nu} {V_{p,j}}_\perp &= w_p * [ DF^{-1}(\eta_p) DF^{-\top}(\eta_p) ]_{\mu, \nu} * {V_{p,i}}_\perp * {V_{p,j}}_\perp\\{V_{p,i}}_\perp B_p^\mu &= w_p * [DF^{-1}(\eta_p)]_\mu * {V_{p,i}}_\perp\end{aligned}\end{align} \]

Note

The above parameter list contains only the model specific input arguments.

struphy.pic.accumulation.accum_kernels.pc_lin_mhd_6d_full()

Accumulates into V1 with the filling functions

\[ \begin{align}\begin{aligned}V_{p,i} A_p^{\mu, \nu} V_{p,j} &= w_p * [ DF^{-1}(\eta_p) DF^{-\top}(\eta_p) ]_{\mu, \nu} * V_{p,i} * V_{p,j}\\V_{p,i} B_p^\mu &= w_p * [DF^{-1}(\eta_p)]_\mu * V_{p,i}\end{aligned}\end{align} \]

Note

The above parameter list contains only the model specific input arguments.

struphy.pic.accumulation.accum_kernels.vlasov_maxwell()

Accumulates into V1 with the filling functions

\[\begin{split}A_p^{\mu, \nu} &= w_p \, G^{-1}_{\mu, \nu}(\boldsymbol \eta_p) \,, \\[2mm] B_p^\mu &= w_p [DF^{-1}(\boldsymbol \eta_p) \cdot \mathbf{v}_p ]_\mu \,.\end{split}\]

Note

The above parameter list contains only the model specific input arguments.

struphy.pic.accumulation.accum_kernels_gc.cc_lin_mhd_5d_D()

Accumulation kernel for the propagator CurrentCoupling5DDensity.

Accumulates \(\alpha\)-form matrix with the filling functions (\(\alpha = 2\))

\[A_p^{\mu, \nu} = w_p \frac{1}{\epsilon} \left( 1-\frac{\hat B_\parallel}{\hat B^*_\parallel} \right) g^{-1} (\mathbf B^2_\times)_{\mu, \nu} \,.\]

Parameters:

• epsilon (float) – scaling factor.

• b2_1 (array[float]) – FE coefficients c_ijk of the magnetic field as a 2-form.

• b2_2 (array[float]) – FE coefficients c_ijk of the magnetic field as a 2-form.

• b2_3 (array[float]) – FE coefficients c_ijk of the magnetic field as a 2-form.

• norm_b11 (array[float]) – FE coefficients c_ijk of the unit magnetic field as a 1-form.

• norm_b12 (array[float]) – FE coefficients c_ijk of the unit magnetic field as a 1-form.

• norm_b12 – FE coefficients c_ijk of the unit magnetic field as a 1-form.

• curl_norm_b1 (array[float]) – FE coefficients c_ijk of the curl of the unit magnetic field as a 2-form.

• curl_norm_b2 (array[float]) – FE coefficients c_ijk of the curl of the unit magnetic field as a 2-form.

• curl_norm_b3 (array[float]) – FE coefficients c_ijk of the curl of the unit magnetic field as a 2-form.

Note

The above parameter list contains only the model specific input arguments.

struphy.pic.accumulation.accum_kernels_gc.cc_lin_mhd_5d_M()

Accumulation kernel for the propagator ShearAlfvenCurrentCoupling5D and MagnetosonicCurrentCoupling5D.

Accumulates 2-form vector with the filling functions:

\[B^\mu_p = \omega_p \mu_p\left(\sqrt{g}^{-1} \hat{\mathbf{b}}¹_0\right)_\mu \,.\]

Parameters:

• norm_b11 (array[float]) – FE coefficients c_ijk of the normalized magnetic field as a 1-form.

• norm_b12 (array[float]) – FE coefficients c_ijk of the normalized magnetic field as a 1-form.

• norm_b13 (array[float]) – FE coefficients c_ijk of the normalized magnetic field as a 1-form.

Note

The above parameter list contains only the model specific input arguments.

struphy.pic.accumulation.accum_kernels_gc.cc_lin_mhd_5d_curlb()

Accumulation kernel for the propagator CurrentCoupling5DCurlb.

Accumulates \(\alpha\)-form matrix and vector with the filling functions (\(\alpha = 2\))

\[ \begin{align}\begin{aligned}A_p^{\mu, \nu} &= w_p \left[\left( \frac{v_{\parallel,p}}{g\hat B^*_\parallel}\right)^2 \mathbf B^2_{\times} \left| \hat \nabla \times \hat{\mathbf b}^1_0 \right|^2 (\mathbf B^2_{\times})^\top \right]_{\mu, \nu}\,,\\B_p^\mu &= w_p \left( \frac{v^2_{\parallel,p}}{g\hat B^*_\parallel} \mathbf B^2_{\times} \right)_\mu \,,\end{aligned}\end{align} \]

where \(\mathbf B^2_{\times} \mathbf a := \hat{\mathbf B}^2 \times \mathbf a\) for \(a \in \mathbb R^3\).

struphy.pic.accumulation.accum_kernels_gc.cc_lin_mhd_5d_gradB()

Accumulation kernel for the propagator CurrentCoupling5DGradB.

Accumulates math:alpha -form vector with the filling functions

\[B_p^\mu &= \omega_p \left[\left(\frac{\mu_p}{\sqrt{g}\hat B^*_\parallel}\right) \mathbf B^2_{\times} G^{-1} \mathbf b^2_{0 \times} G^{-1} \nabla B_\parallel¹\right]_\mu \,,\]

where \(B2_{\times} * a := B2 \times a\) for \(a \in \mathbb R^3\).

Parameters:

• b1 (array[float]) – FE coefficients c_ijk of the magnetic field as a 2-form.

• b2 (array[float]) – FE coefficients c_ijk of the magnetic field as a 2-form.

• b3 (array[float]) – FE coefficients c_ijk of the magnetic field as a 2-form.

• norm_b11 (array[float]) – FE coefficients c_ijk of the normalized magnetic field as a 1-form.

• norm_b12 (array[float]) – FE coefficients c_ijk of the normalized magnetic field as a 1-form.

• norm_b13 (array[float]) – FE coefficients c_ijk of the normalized magnetic field as a 1-form.

• curl_norm_b1 (array[float]) – FE coefficients c_ijk of the curl of normalized magnetic field as a 2-form.

• curl_norm_b2 (array[float]) – FE coefficients c_ijk of the curl of normalized magnetic field as a 2-form.

• curl_norm_b3 (array[float]) – FE coefficients c_ijk of the curl of normalized magnetic field as a 2-form.

• grad_PB1 (array[float]) – FE coefficients c_ijk of gradient of parallel magnetic field as a 1-form.

• grad_PB2 (array[float]) – FE coefficients c_ijk of gradient of parallel magnetic field as a 1-form.

• grad_PB3 (array[float]) – FE coefficients c_ijk of gradient of parallel magnetic field as a 1-form.

Note

The above parameter list contains only the model specific input arguments.

struphy.pic.accumulation.accum_kernels_gc.cc_lin_mhd_5d_gradB_dg()

TODO

struphy.pic.accumulation.accum_kernels_gc.cc_lin_mhd_5d_gradB_dg_init()

TODO

struphy.pic.accumulation.accum_kernels_gc.gc_density_0form()

Kernel for AccumulatorVector into V0 with the filling

\[B_p^\mu = \frac{w_p}{N} \,.\]

struphy.pic.accumulation.accum_kernels_gc.gc_mag_density_0form()

Kernel for AccumulatorVector into V0 with the filling

\[B_p^\mu = \mu \frac{w_p}{N} \,.\]

Pusher class

class struphy.pic.pushing.pusher.Pusher(particles: Particles, kernel: Pyccelkernel, args_kernel: tuple, args_domain: DomainArguments, *, alpha_in_kernel: float | int | tuple | list, init_kernels: list = [], eval_kernels: list = [], n_stages: int = 1, maxiter: int = 1, tol: float = 1e-08, mpi_sort: str | None = None)

Bases: object

Class for solving particle ODEs

\[\dot{\mathbf Z}_p(t) = \mathbf U(t, \mathbf Z_p(t))\,,\]

for each marker \(p\) in Particles class, where \(\mathbf Z_p\) are the marker coordinates and the vector field \(\mathbf U\) can contain discrete Derham splines and metric coefficients from accelerated evaluation_kernels.

The solve is MPI distributed and can handle multi-stage Runge-Kutta methods for any ButcherTableau as well as iterative nonlinear methods.

The particle push is performed via accelerated pusher_kernels or pusher_kernels_gc for guiding-center models.

Notes

For iterative methods with iteration index \(k\), spline evaluations at positions \(\alpha_i \eta_{p,i}^{n+1,k} + (1 - \alpha_i) \eta_{p,i}^n\) for \(i=1, 2, 3\) and different \(\alpha_i \in [0,1]\) need particle MPI sorting in between. This requires calling dedicated eval_kernels during the iteration. Here are some rules to follow for iterative solvers:

• Spline/geometry evaluations at \(\boldsymbol \eta^n_p\) can be be done via init_kernels.

• Pusher kernel and eval_kernels can perform evaluations at arbitrary weighted averages \(\eta_{p,i} = \alpha_i \eta_{p,i}^{n+1,k} + (1 - \alpha_i) \eta_{p,i}^n\), for \(i=1,2,3\).

• MPI sorting is done automatically before kernel calls according to the specified values \(\alpha_i\) for each kernel.

Parameters:

• particles (Particles) – Particles object holding the markers to push.

• kernel (pyccelized function) – The pusher kernel.

• args_kernel (tuple) – Optional arguments passed to the kernel.

• args_domain (DomainArguments) – Mapping infos.

• alpha_in_kernel (float | int | tuple | list) – For i=0,1,2, the spline/geometry evaluations in kernel are at alpha[i]*markers[:, i] + (1 - alpha[i])*markers[:, buffer_idx + i]. If float or int or then alpha = (alpha, alpha, alpha). alpha must be between 0 and 1. alpha[i]=0 means that evaluation is at the initial positions (time n), stored at markers[:, buffer_idx + i].

• init_kernels (dict) – Keys: initialization kernels for spline/ SPH evaluations at time n (initial state). Values: optional arguments.

• eval_kernels (dict) – Keys: evaluation kernels for splines before the pusher kernel is called. Values: optional arguments and weighting parameters alpha for sorting (before evaluation), according to alpha[i]*markers[:, i] + (1 - alpha[i])*markers[:, buffer_idx + i] for i=0,1,2. alpha must be between 0 and 1, see mpi_sort_markers().

• n_stages (int) – Number of stages of the pusher (e.g. 4 for RK4)

• maxiter (int) – Maximum number of iterations (=1 for explicit pushers).

• tol (float) – Iteration terminates when residual<tol.

• mpi_sort (str) – When to do MPI sorting: * None : no sorting at all. * each : sort markers after each stage. * last : sort markers after last stage.

__call__(dt: float)

Applies the chosen pusher kernel by a time step dt, applies kinetic boundary conditions and performs MPI sorting.

property particles

Particle object.

property kernel

The pyccelized pusher kernel.

property init_kernels

A dict of kernels for initial spline evaluation before iteration.

property eval_kernels

A dict of kernels for spline evaluation before execution of kernel during iteration.

property args_kernel

Optional arguments for kernel.

property args_domain

Mandatory Domain arguments.

property n_stages

Number of stages of the pusher.

property maxiter

Maximum number of iterations (=1 for explicit pushers).

property tol

Iteration terminates when residual<tol.

__weakref__

list of weak references to the object (if defined)

property mpi_sort

When to do MPI sorting: * None : no sorting at all. * each : sort markers after each stage. * last : sort markers after last stage.