PIC#

Particle base class#

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

Bases: object

Base class for particle species.

The marker information is stored in a 2D array of shape (n_cols, n_markers_max), one column per marker (see struphy.pic.marker_layout). In markers[j, ip] the index ip refers to a specific particle, the index j to its attributes; n_cols is the number of attributes and n_markers_max the number of markers (including holes), and the attribute indices are called “columns” below and in the rest of the code. They are as follows:

0:3: position in the logical unit cube ($\boldsymbol \eta_p \in [0, 1]^3$)
3:3 + vdim: velocities
3 + vdim: (time-dependent) weight $w_k(t)$
4 + vdim: PDF $s^0 = s^3/\sqrt g$ at particle position
5 + vdim: initial weight $w_0$
6 + vdim <= j < -2: buffer columns, laid out in consecutive blocks (each block’s starting column and width are given by a pair of attributes/properties):
- first_diagnostics_idx (width n_cols_diagnostics): free columns for model-specific diagnostics (e.g. canonical momentum, magnetic moment, …).
- first_pusher_idx (width n_cols_pusher $= 3 + \mathrm{vdim}$): scratch space for a Pusher call, used to hold the phase space coordinates at the start of a push (or of a sub-stage, for multi-stage pushers).
- first_shift_idx (width n_cols_shift $= 3$): accumulated shifts in $\eta$-space due to boundary conditions (e.g. periodic wrap-around), added back onto the pusher’s initial positions when reconstructing a marker’s unwrapped trajectory.
- residual_idx (width 1): residual of the current iteration, for pushers that solve a nonlinear/implicit equation iteratively.
- first_free_idx (width n_cols_aux): general-purpose scratch columns available to any routine that needs temporary per-marker storage (e.g. field evaluations).
The total number of columns is given by n_cols, i.e. first_free_idx + n_cols_aux + 2 (the + 2 accounts for the last two columns below).
-2: number of the sorting box the particle is in
-1: particle ID

Direct indexing into markers is rarely needed outside of the Particles class itself. Instead, the most commonly used columns are exposed as convenience properties, each returning (or setting) a 2D array of shape (n_mks_loc, ...) restricted to the valid markers on the current process (i.e. markers[..., self.valid_mks]):

positions (columns 0:3): marker positions $\boldsymbol\eta_p$.
velocities (columns 3:3 + vdim): marker velocities.
phasespace_coords (columns 0:3 + vdim): positions and velocities combined.
weights (column 3 + vdim): current weights $w_k(t)$.
sampling_density_values (column 4 + vdim): PDF $s^0$ at the particle position.
weights0 (column 5 + vdim): initial weights $w_0$.
marker_ids (column -1): unique particle IDs.

Each of these properties has a matching setter (e.g. self.positions = new_positions) that validates the shape of the assigned array and writes it back into self._markers at the corresponding columns, for the valid markers only.

Parameters:

comm_world (Intracomm) – World MPI communicator.
clone_config (CloneConfig) – Manages the configuration for clone-based (copied grids) parallel processing using MPI.
domain_decomp (tuple) – The first entry is a domain_array (see domain_array) and the second entry is the number of MPI processes in each direction.
loading_params (LoadingParameters) – Parameterts for particle loading.
weights_params (WeightsParameters) – Parameters for particle weights.
boundary_params (BoundaryParameters) – Parameters for particle boundary conditions.
sorting_params (SortingParameters) – Parameters for particle sorting.
saving_params (SavingParameters) – Parameters for particle saving.
bufsize (float) – Size of buffer (as multiple of total size, default=.25) in markers array.
domain (Domain) – Struphy domain object.
equil (FluidEquilibrium) – Struphy fluid equilibrium object.
projected_equil (ProjectedFluidEquilibrium) – Struphy fluid equilibrium projected into a discrete Derham complex.
background (KineticBackground) – Kinetic background.
initial_condition (KineticBackground) – Kinetic initial condition.
n_as_volume_form (bool) – Whether the number density n is given as a volume form or scalar function (=default).
perturbations (Perturbation | list) – Kinetic perturbation parameters.
equation_params (dict) – Normalization parameters (epsilon, alpha, …)
dry_run (bool) – If True, only compute the sizing of the marker array (n_markers_max, n_cols, …) and return early, without allocating any of the (potentially large) marker/sorting/buffer arrays. Used by nbytes_local to estimate the memory footprint before actually allocating the particles, see estimate_mem().

abstract property vdim#: Dimension of the velocity space.

abstract property coordinate_labels: tuple[str]#: Labels for the coordinates in the phase space. Length must be 3 + vdim, where the first 3 are the spatial coordinates and the last vdim are the velocity coordinates.

orbit_quantities: tuple[tuple[int, str, str, str], ...] = ()#: Marker columns saved in the orbits post-processing product, as (column, name, long_name, description). The first three are ORBIT_POSITIONS; the marker index becomes the marker coordinate of the product.

abstract property mu_idx#: Index of the column in the marker array where the magnetic moment is stored.

abstract property default_background#: The default background (of type Maxwellian).

abstract property default_n_cols#

12} for default number of columns.

Type:: Dictionary of the form {‘diagnostics’
Type:: 3, ‘aux’

abstract __post_init__()[source]#: Can be used for checks on the constructor arguments and for setting additional attributes in subclasses.

abstract property sampling_density#: Marker sampling density function $s^ extrm{vol}$ as a volume form, see Monte-Carlo integrals. Must be normalized to 1. Its coordinates are the coordinates used in Monte-Carlo Integrals approximated by the particles.

abstract s0(eta1, eta2, eta3, *v, flat_eval=False, remove_holes=True)[source]#: 0-form corresponding to ~struphy.pic.base.Particles.sampling_density. This is the quantity stored in each marker’s s0 column (see the class docstring) and used to compute initial weights w0 = f_init / s0 / Np.

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#: Number of quantities (“columns”) per marker, the length of the first axis of the markers array. The last 2 columns refer to box number and particle ID, respectively.

property n_markers_max#

the most markers this process can hold.

Type:: Number of marker slots (columns of markers), holes included

property mean_velocity_index#: Index in marker array where mean velocity for noslip BC is stored.

property nbytes_local: int#: Estimated local (per-MPI-rank) memory footprint, in bytes, of all marker-related arrays (markers, sorting buffers, lost-marker container). Only depends on n_markers_max and n_cols, so it is valid whether or not the arrays were actually allocated (see the dry_run argument of __init__()).

property kinds#: Name of the class.

property index#: Dict holding the column indices referring to specific marker parameters (coordinates).

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 loading_params: LoadingParameters#: Parameters for particle loading.

property weights_params: WeightsParameters#: Parameters for particle weights.

property boundary_params: BoundaryParameters#: Parameters for marker loading.

property sorting_params: SortingParameters#: Parameters for marker sorting.

property saving_params: SavingParameters#: Parameters for marker/distribution function saving.

property loading: Literal['pseudo_random', 'sobol_standard', 'sobol_antithetic', 'external', 'restart', 'tesselation']#: Type of particle loading.

property spatial#: Drawing particles uniformly on the unit cube(‘uniform’) or on the disc(‘disc’)

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 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 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 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 n_mks_load#: Array of number of markers on each process at loading stage

property markers#

2D array of shape (n_cols, n_markers_max) holding the marker information, including holes. The i-th column holds the i-th marker info (see struphy.pic.marker_layout).

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 row 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 column 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 column holds the i-th marker info.

property lost_markers#: Array containing the last infos of removed markers

property n_lost_markers#: Number of removed particles.

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 of shape (3, n_mks_loc) holding the marker positions in logical space. The i-th column holds the i-th marker info.

property velocities#: Array of shape (vdim, n_mks_loc) holding the marker velocities in logical space. The i-th column holds the i-th marker info.

property phasespace_coords#: Array of shape (3 + vdim, n_mks_loc) holding the marker positions and velocities in logical space. The i-th column holds the i-th marker info.

property weights#: Array holding the current marker weights. The i-th entry holds the i-th marker info.

property sampling_density_values#: Array holding the current marker 0form sampling density s0. The i-th entry holds the i-th marker info.

property weights0#: Array holding the initial marker weights. The i-th entry holds the i-th marker info.

property marker_ids#: Array holding the marker id’s on the current process.

property f_coords#: Coordinates of the distribution function.

property f_jacobian_coords#: Coordinates of the velocity jacobian determinant of the distribution fuction.

property args_markers: MarkerArguments | CudaMarkerArguments#: Arguments for marker kernels, selected when this particle object is created.

kernels() → tuple[source]#

The kernels these particles call when they apply their boundary conditions and sort into boxes.

reflect() if an axis has reflecting boundary conditions, and the box kernels of struphy.pic.kernels if there are sorting boxes. Used by compile_cuda_kernels().

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 reject_weights#: Whether to reect weights below threshold.

property threshold#: Threshold for rejecting weights.

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 f0: Maxwellian#: Callable background distribution function, used as the control variate in update_weights().

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 control_variate#: Boolean for whether to use the Control variate method during time stepping.

property boxes_per_dim#: Tuple, number of sorting boxes per dimension.

property sorting_boxes#: The SortingBoxes instance holding the sorting-box data structure used by put_particles_in_boxes() and do_sort().

property tesselation#: Tesselation of the current process domain.

property equation_params#: Parameters appearing in model equation due to Struphy normalization.

property domain: Domain#: From struphy.geometry.domains.

property equil: FluidEquilibrium#: From struphy.fields_background.equils.

property projected_equil: ProjectedFluidEquilibrium#: MHD equilibrium projected on 3d Derham sequence with commuting projectors.

draw_markers(sort: bool = True)[source]#

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:

Uniform distribution on the unit cube: $n^3(\eta) = 1$
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 Particles5Dvperp, 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.

On several processes, the markers are always sent to the process of their domain after drawing (mpi_sort_markers()). They are also sorted into boxes if sort is True (default) and sorting boxes were set up.

Parameters:: sort (Bool) – Whether to sort the markers into boxes after initial drawing (only if sorting boxes were set up).

initialize_weights(*, bckgr_params: dict | None = None, pert_params: dict | None = None)[source]#

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()[source]#: 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.

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)[source]#

Computes full-f and delta-f distribution functions via marker binning in logical space. xp.histogramdd (NumPy or CuPy, on the device of the markers) 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 (default: True).

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: list[bool], bin_edges: list[ndarray], do_plot=False)[source]#

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 3+vdim giving the directions in phase space in which to bin. Up to two entries can be True, the rest must be False. The True entries correspond to the axes of the binning.
bin_edges (list[np.ndarray]) – List of bin edges (resolution) having the length of True entries in components.
do_plot (bool) – Whether to show the plot (default: False).

Returns:

err – Maximum relative error between the binned distribution function and the analytic initial condition.

Return type:

float

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

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:

apply_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.

Notes

On the CuPy backend, the markers are exchanged as device buffers with CUDA-aware MPI (staged through host memory otherwise), see _sendrecv_markers_device(). The per-rank counts are small host arrays on every backend (_sendrecv_all_to_all()).

apply_kinetic_bc(newton=False)[source]#

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.

finish_kernel_bc(newton=False)[source]#

Bookkeeping after a pusher kernel that applied the kinetic boundary conditions per marker (see apply_kinetic_bc_marker()). Refilling is not done in the kernel and is applied here by apply_kinetic_bc(). Markers removed in the kernel have become holes; they are counted as lost markers.

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

update_holes()[source]#: Recompute the holes mask (markers with markers[0] == -1) and, from it, refresh valid_mks. Must be called after any operation that creates, removes or moves markers (e.g. sorting, boundary conditions, refilling), since holes are tracked per marker index. The markers (columns) from n_markers_used on are holes and are not read.

property n_markers_used: int#

Number of marker slots (columns of markers, out of n_markers_max) that can hold markers: every slot from this one on is a hole.

Kernels loop over these markers only (args_markers.n_markers), and so do the row operations of the pushers and the boundary conditions (markers[j, :n_markers_used]). It is an upper bound kept on the host, so that it never waits for the device: markers always fill the lowest holes, so k incoming markers end below (occupied slots) + k (see _track_markers()). Loading and box sorting count it exactly (count_markers_used()), and so must code that writes markers into holes itself, outside the methods of this class.

count_markers_used()[source]#

Set n_markers_used to one past the last column that is not a hole, and count the occupied columns.

Called after loading and after box sorting; call it, before update_holes(), after writing markers into holes outside the methods of this class. Reads two numbers from the device on the CuPy backend.

set_velocities_comp(velocity, comp)[source]#

Set one or several velocity components to the same constant value, for all valid markers.

Parameters:

velocity (float) – The constant value to assign to the selected velocity components.
comp (iterable[int]) – Velocity components to set (0-based, e.g. 0 for v1, 1 for v2, …).

put_particles_in_boxes()[source]#: 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 neighbouring processes are also communicated (as ghost particles).

do_sort(use_numpy_argsort=False)[source]#

Assign the particles to their sorting boxes and reorder the markers array accordingly, so that markers in the same box occupy contiguous columns.

The box array boxes is rebuilt for the sorted markers, so it stays valid after the sort.

Parameters:: use_numpy_argsort (bool) – If True, sort via numpy.argsort() on the box row; if False (default), use the kernel sort_and_fill_boxes().

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

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[source]#

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 the Pyccel kernel sph_mean_velocity_coeffs 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[source]#

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 in Piola form). It then evaluates the necessary logical derivatives via eval_sph(), sums contributions and divides by the Jacobian determinant to produce the three Cartesian divergence components.

classmethod ker_dct()[source]#: Dict mapping the name of each available SPH smoothing kernel (e.g. "gaussian_1d") to its integer kernel ID, used in the Pyccel kernels. Kernel IDs must have three digits, see smoothing_kernel().

gather_scalar_in_subcomm_array(scalar: int, out: ndarray | None = None)[source]#

Return an array of length sub_comm.size, where the i-th entry corresponds to the value of the scalar on process i.

Parameters:

scalar (int) – The scalar value on each process.
out (xp.ndarray) – The returned array (optional).

gather_scalar_in_intercomm_array(scalar: int, out: ndarray | None = None)[source]#

Return an array of length inter_comm.size, where the i-th entry corresponds to the value of the scalar on clone i.

Parameters:

scalar (int) – The scalar value on each clone.
out (xp.ndarray) – The returned array (optional).

__weakref__#: list of weak references to the object (if defined)

Marker array layout#

Memory layout of the marker buffer.

Particles.markers has shape (n_cols, n_markers_max), one column per marker, and is stored in C order (struct-of-arrays): every quantity (row j) is one contiguous block of n_markers_max values. markers[j, ip] is the j-th entry of marker ip; markers[j] is quantity j of all markers, and markers[:, ip] all quantities of marker ip. Operations on quantities (markers[j], markers[a:b] = ..., a mask on one row) then touch only the rows involved, and on the GPU neighbouring threads (markers) read neighbouring addresses.

The same convention, one column per point, holds for every array of markers or points passed to an evaluation (Domain, SplineFunction, the equilibria): (k, n) with the logical coordinates in rows 0 to 2.

The kernels take C-ordered arrays of this shape. Selecting markers (markers[:, valid_mks]) gives a Fortran-ordered copy (NumPy keeps the selected axis innermost): fine to read or to write back with markers[:, idx] = ..., but convert it with xp.ascontiguousarray before passing it to a kernel or to MPI.

See struphy-hub/struphy#707.

struphy.pic.marker_layout.MARKER_ALIGNMENT = 32#: The number of marker slots n_markers_max is a multiple of this, so that every row of the buffer starts on a 256-byte boundary (32 doubles).

struphy.pic.marker_layout.padded_n_markers(n_markers_max: int) → int[source]#: Smallest multiple of MARKER_ALIGNMENT that is at least n_markers_max.

struphy.pic.marker_layout.marker_array(n_cols: int, n_markers_max: int)[source]#: Zero marker array of shape (n_cols, n_markers_max) on the active backend.

Particle subclasses#

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

Bases: Particles

Particles in the full 6D phase space $(\boldsymbol \eta, \mathbf v) \in [0, 1]^3 \times \mathbb R^3$, as used e.g. in full-orbit (Vlasov) kinetic models.

Each marker carries a logical (curvilinear) position $\boldsymbol \eta_p$ together with a velocity $\mathbf v_p$ expressed in the Cartesian velocity space attached to that position (i.e. velocities are not transformed by the curvilinear map, unlike positions).

See Particles for the structure of the numpy marker array and the meaning of its columns.

vdim = 3#: Dimension of the (Cartesian) velocity space, here 3.

coordinate_labels = ('$\\eta_1$', '$\\eta_2$', '$\\eta_3$', '$v_x$', '$v_y$', '$v_z$')#: Labels for the coordinates in the phase space. Length is 6, with the first 3 being the spatial coordinates and the last 3 being the velocity coordinates.

orbit_quantities: tuple[tuple[int, str, str, str], ...] = ((0, 'x', '$x$', 'physical position x'), (1, 'y', '$y$', 'physical position y'), (2, 'z', '$z$', 'physical position z'), (3, 'v1', '$v_x$', 'Cartesian velocity x'), (4, 'v2', '$v_y$', 'Cartesian velocity y'), (5, 'v3', '$v_z$', 'Cartesian velocity z'), (6, 'weight', '$w$', 'marker weight'))#: Marker columns saved as orbits, see orbit_quantities.

default_background = Maxwellian3D( n=(1.0, None), u1=(0.0, None), u2=(0.0, None), u3=(0.0, None), vth1=(1.0, None), vth2=(1.0, None), vth3=(1.0, None), uniform_on_disc=False, )#: Default kinetic background is a 3D Cartesian Maxwellian.

default_n_cols = {'aux': 6, 'diagnostics': 0}#: Default number of buffer columns reserved for diagnostics and auxiliary (pusher/free) use (6 for the symmetric 3x3 particle-matrix fillings stored per marker by the matrix-free accumulation, see MarkerOperator).

property mu_idx#: Index of the column in the marker array where the magnetic moment is stored.

__post_init__()[source]#: If the background is a CanonicalMaxwellian, set up the discrete magnetic field (needed to evaluate canonical invariants) from the projected equilibrium.

property sampling_density#

Sampling density function as volume form, used to draw markers via inverse transform/rejection sampling and to compute their initial weights (see draw_markers()).

This is a Maxwellian3D in the Cartesian velocities vx, vy, vz, parametrized by the mean velocities and thermal velocities in loading_params, with density normalized to 1 (i.e. uniform in eta1, eta2, eta3), further multiplied by the Jacobian factor 2 * eta1 if spatial is "disc" (to sample uniformly in physical space on a disc, where eta1 plays the role of a normalized radius).

s0(eta1, eta2, eta3, vx, vy, vz, flat_eval=False, remove_holes=True)[source]#

Sampling density function as 0 form. This is the quantity stored in each marker’s s0 column (see the class docstring) and used to compute initial weights w0 = f_init / s0 / Np.

Parameters:

eta1 (array_like) – Logical evaluation points.
eta2 (array_like) – Logical evaluation points.
eta3 (array_like) – Logical evaluation points.
vx (array_like) – Cartesian velocity evaluation points.
vy (array_like) – Cartesian velocity evaluation points.
vz (array_like) – Cartesian velocity evaluation points.
flat_eval (bool) – If true, perform flat (marker) evaluation (etas must be same size 1D).
remove_holes (bool) – If True, holes are removed from the returned array. If False, holes are evaluated to -1.

Returns:

out (array-like) – The 0-form sampling density.
——-

save_constants_of_motion()[source]#

Calculate each marker’s guiding-center constants of motion (only the equilibrium magnetic field is considered) and assign them into the diagnostics columns of the marker array:

0:3: guiding-center position (logical $\boldsymbol \eta$)
3: energy
4: magnetic moment
5: canonical toroidal momentum
6: parallel velocity

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

Bases: Particles6D

A class for kinetic species in full 6D phase space that solve for delta_f = f - f0.

See Particles6D for more information.

__post_init__()[source]#: Force the control-variate weight update off, since delta-f weights already evolve the perturbation directly (there is no separate background contribution to subtract).

set_n_to_zero(background: Maxwellian | SumKineticBackground)[source]#

Recursively set the density moment n of background (and, if it is a SumKineticBackground, of both its summands) to zero, keeping any perturbation attached to it.

Parameters:: background (Maxwellian | SumKineticBackground) – The kinetic background whose density is to be zeroed.

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

Bases: Particles6D

Classical Pauli particles (Xiao and Qin, Comput. Phys. Commun. 265, 107981 (2021)): full-orbit (6D) markers with an additional, constant magnetic moment $\mu_p$, stored in the first diagnostics column (mu_idx).

The markers are loaded on the slow manifold of the classical Pauli particle, whose dynamics is that of the guiding center: after drawing the markers from the initial condition as for Particles6D (and computing their initial weights and sampling densities), each marker gets the magnetic moment of its perpendicular velocity and its velocity is reduced to the parallel component,

\[\mu_p = \frac{|\mathbf v_p \times \mathbf b_0(\mathbf x_p)|^2}{2 |\mathbf B_0(\mathbf x_p)|}\,,\qquad \mathbf v_p \to (\mathbf v_p \cdot \mathbf b_0)\, \mathbf b_0(\mathbf x_p)\,,\]

see reduce_to_slow_manifold(). The marker position is from then on the guiding-center position, and the markers represent the gyro-average of the initial condition. The magnetic moment of the remaining perpendicular velocity, $\mu'_p$ (mu_prime()), stays close to zero as long as the dynamics remains on the slow manifold.

With the control variate, the distribution functions (background and initial condition) are evaluated at the gyrotropic velocity $v_\parallel \mathbf b_0 + \sqrt{2 \mu_p |\mathbf B_0|}\, \mathbf e_\perp$ with some unit vector $\mathbf e_\perp \perp \mathbf b_0$ (see f_coords), i.e. the background is the guiding-center distribution $f_0(\mathbf x, v_\parallel, \mu)$. It must therefore be gyrotropic: Maxwellians with isotropic thermal velocities and a mean velocity parallel to $\mathbf B_0$ (no diamagnetic drift). The polarization drift on the slow manifold is compressible, so the phase-space density of the markers is not constant along their orbits, and subtracting $f_0$ along the markers removes the polarization density from the deposited charge: it must then be added in the field equation (as for guiding centers, see AdiabaticPoissonSolve with polarization=True).

See Particles for the structure of the numpy marker array and the meaning of its columns.

orbit_quantities: tuple[tuple[int, str, str, str], ...] = ((0, 'x', '$x$', 'physical position x'), (1, 'y', '$y$', 'physical position y'), (2, 'z', '$z$', 'physical position z'), (3, 'v1', '$v_x$', 'Cartesian velocity x'), (4, 'v2', '$v_y$', 'Cartesian velocity y'), (5, 'v3', '$v_z$', 'Cartesian velocity z'), (6, 'weight', '$w$', 'marker weight'), (9, 'mu', '$\\mu$', 'magnetic moment'))#

those of Particles6D and the magnetic moment.

Type:: Marker columns saved as orbits

default_n_cols = {'aux': 6, 'diagnostics': 1}#: One diagnostics column holds the magnetic moment.

property mu_idx#: Column of the magnetic moment (the first diagnostics column).

__post_init__()[source]#: Check that there is a magnetic equilibrium and, with the control variate, that the kinetic backgrounds are gyrotropic Maxwellians (isotropic thermal velocities).

property magnetic_moments#: Array of shape (n_mks_loc,) holding the magnetic moments of the markers.

initialize_weights(**kwargs)[source]#: Compute the initial weights and sampling densities as for Particles6D with the drawn velocities, then reduce the markers to the slow manifold (reduce_to_slow_manifold()). Markers loaded from a file (loading="restart" or "external") are already on the slow manifold and carry their magnetic moments (column mu_idx), so they are not reduced again.

reduce_to_slow_manifold()[source]#: Set the magnetic moments $\mu_p = |\mathbf v_{\perp,p}|^2 / (2 |\mathbf B_0|)$ and remove the perpendicular velocities, $\mathbf v_p \to (\mathbf v_p \cdot \mathbf b_0) \mathbf b_0$.

mu_prime()[source]#: Magnetic moment of the perpendicular marker velocities, $\mu'_p = |\mathbf v_p \times \mathbf b_0|^2 / (2 |\mathbf B_0|)$ (Eq. (7) of Xiao and Qin), array of shape (n_mks_loc,). It is close to zero on the slow manifold.

property f_coords#

the positions and, on the slow manifold, the gyrotropic velocities $v_\parallel \mathbf b_0 + \sqrt{2 \mu_p |\mathbf B_0|}\, \mathbf e_\perp$.

Type:: Coordinates of the distribution function

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

Bases: Particles

Particles in the 5D guiding-center, drift-kinetic or gyro-kinetic phase space $(\boldsymbol \eta, v_\parallel, \mu) \in [0, 1]^3 \times \mathbb R \times \mathbb R_{\geq 0}$.

Each marker carries a logical (curvilinear) position $\boldsymbol \eta_p$ together with the velocity coordinates

\[v_{\parallel, p} = \mathbf v_p \cdot \mathbf b_0(\boldsymbol \eta_p) \,, \qquad \mu_p = \frac{1}{2 |\mathbf B_0|} m |\mathbf v_p|^2 - v_{\parallel, p}^2 \,,\]

defined with respect to the equilibrium magnetic field $\mathbf B_0$ and its unit vector $\mathbf b_0 = \mathbf B_0 / |\mathbf B_0|$ (unlike Particles6D, velocities are thus not Cartesian but expressed in a field-aligned basis that itself depends on $\boldsymbol \eta_p$).

By default, two diagnostics columns are reserved (default_n_cols["diagnostics"] = 2), holding each marker’s perpendicular energy and canonical toroidal momentum (see save_constants_of_motion()).

See Particles for the structure of the numpy marker array and the meaning of its columns.

vdim = 2#: Dimension of the velocity space, here 2 ($v_\parallel, \mu$).

coordinate_labels = ('$\\eta_1$', '$\\eta_2$', '$\\eta_3$', '$v_\\parallel$', '$\\mu$')#: Labels for the coordinates in the phase space. Length is 5, with the first 3 being the spatial coordinates and the last 2 being the velocity coordinates.

orbit_quantities: tuple[tuple[int, str, str, str], ...] = ((0, 'x', '$x$', 'physical position x'), (1, 'y', '$y$', 'physical position y'), (2, 'z', '$z$', 'physical position z'), (3, 'v_par', '$v_\\parallel$', 'parallel velocity'), (4, 'mu', '$\\mu$', 'magnetic moment'), (5, 'weight', '$w$', 'marker weight'), (9, 'p_phi', '$p_\\phi$', 'canonical toroidal momentum (set by save_constants_of_motion)'))#: Marker columns saved as orbits, see orbit_quantities.

mu_idx = 4#: Column index of particle magnetic moment.

default_background = GyroMaxwellian2D( n=(1.0, None), u_para=(0.0, None), u_perp=(0.0, None), vth_para=(1.0, None), vth_perp=(1.0, None), volume_form=True, B0=2.0, uniform_on_disc=False, )#: Default kinetic background is a gyrotropic Maxwellian in $(v_\parallel, \mu)$.

default_n_cols = {'aux': 12, 'diagnostics': 2}#: Default number of buffer columns is 2 diagnostics (perpendicular energy, canonical toroidal momentum, see save_constants_of_motion()) and 12 auxiliary columns.

__post_init__()[source]#: Retrieve the discrete equilibrium magnetic-field quantities ($|B_0|$, unit 1-form $\mathbf b_0$, Derham complex) needed to project marker velocities onto $v_\parallel, \mu$ and to evaluate diagnostics, and allocate the temporary FE coefficient vectors used for that.

property magn_bckgr#: Equilibrium fluid background carrying the magnetic field $\mathbf B_0$ with respect to which $v_\parallel, \mu$ are defined.

property absB0_h#: Discrete 0-form coefficients of $|B_0|$.

property unit_b1_h#: Discrete 1-form coefficients of the equilibrium field-aligned unit vector $\mathbf b_0 = \mathbf B_0/|B_0|$.

property epsilon#: Normalization parameter $\epsilon$ (from equation_params) entering the guiding-center equations of motion, e.g. the canonical toroidal momentum evaluation.

property derham#: Discrete Derham complex of the projected equilibrium.

property sampling_density#

Sampling density function as volume form, used to draw markers via inverse transform/rejection sampling and to compute their initial weights (see draw_markers()).

This is a GyroMaxwellian2D in $(v_\parallel, \mu)$, parametrized by the mean/thermal parallel velocity and by the equilibrium magnetic field in loading_params. It is normalized to 1 in logical space (i.e. uniform in eta1, eta2, eta3) and already includes the polar-coordinate Jacobian factor $|\mathbf B_0|$ (volume_form=True), further multiplied by 2 * eta1 if spatial is "disc".

s3(eta1, eta2, eta3, v_para, mu)[source]#

Sampling density function as 3-form, i.e. sampling_density() with the velocity-space ($B_0$) Jacobian factor divided back out, leaving a density that is a volume form in $\boldsymbol \eta$ only.

Parameters:

eta1 (array_like) – Logical evaluation points.
eta2 (array_like) – Logical evaluation points.
eta3 (array_like) – Logical evaluation points.
v_para (array_like) – Parallel velocity and magnetic moment evaluation points.
mu (array_like) – Parallel velocity and magnetic moment evaluation points.

Returns:

out (array-like) – The 3-form sampling density.
——-

s0(eta1, eta2, eta3, v_para, mu, flat_eval=False, remove_holes=True)[source]#

Sampling density function as 0-form, i.e. s3() pushed forward to a pointwise density by dividing out the spatial metric Jacobian determinant. This is the quantity stored in each marker’s s0 column and used to compute initial weights w0 = f_init / s0 / Np.

Parameters:

eta1 (array_like) – Logical evaluation points.
eta2 (array_like) – Logical evaluation points.
eta3 (array_like) – Logical evaluation points.
v_para (array_like) – Parallel velocity and magnetic moment evaluation points.
mu (array_like) – Parallel velocity and magnetic moment evaluation points.
flat_eval (bool) – If true, perform flat (marker) evaluation (etas must be same size 1D).
remove_holes (bool) – If True, holes are removed from the returned array. If False, holes are evaluated to -1.

Returns:

out (array-like) – The 0-form sampling density.
——-

save_constants_of_motion()[source]#

Calculate each marker’s guiding-center energy and canonical toroidal momentum (only the equilibrium magnetic field is considered) and assign them into the diagnostics columns of the marker array:

first_diagnostics_idx + 0: energy
first_diagnostics_idx + 1: canonical toroidal momentum

The magnetic moment itself is not a diagnostics column here (unlike in Particles5Dvperp) since it is already a phase-space coordinate, see mu_idx.

save_magnetic_energy(PBb)[source]#

Calculate the (time-dependent) magnetic field energy at each marker’s position and assign it into the energy diagnostics column (self.first_diagnostics_idx).

Parameters:: PBb (BlockVector) – Finite element coefficients of the time-dependent magnetic field, projected onto V0.

save_magnetic_background_energy()[source]#: Evaluate the equilibrium magnetic-moment energy $\mu_p |B_0(\boldsymbol \eta_p)|$ for each marker. The result is stored in the energy diagnostics column (self.first_diagnostics_idx).

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

Bases: Particles

Particles in the 5D guiding-center, drift-kinetic or gyro-kinetic phase space $(\boldsymbol \eta, v_\parallel, v_\perp) \in [0, 1]^3 \times \mathbb R \times \mathbb R_{\geq 0}$.

Each marker carries a logical (curvilinear) position $\boldsymbol \eta_p$ together with the parallel and perpendicular velocity coordinates

\[v_{\parallel, p} = \mathbf v_p \cdot \mathbf b_0(\boldsymbol \eta_p) \,, \qquad v_{\perp, p} = \left| \mathbf v_p - v_{\parallel, p} \, \mathbf b_0(\boldsymbol \eta_p) \right| \,,\]

defined with respect to the equilibrium magnetic field $\mathbf B_0$ and its unit vector $\mathbf b_0 = \mathbf B_0 / |\mathbf B_0|$ (unlike Particles6D, velocities are thus not Cartesian but expressed in a field-aligned basis that itself depends on $\boldsymbol \eta_p$).

By default, three diagnostics columns are reserved (default_n_cols["diagnostics"] = 3), holding each marker’s guiding-center energy, magnetic moment and canonical toroidal momentum (see save_constants_of_motion()).

See Particles for the structure of the numpy marker array and the meaning of its columns.

vdim = 2#: Dimension of the velocity space, here 2 ($v_\parallel, v_\perp$).

coordinate_labels = ('$\\eta_1$', '$\\eta_2$', '$\\eta_3$', '$v_\\parallel$', '$v_\\perp$')#: Labels for the coordinates in the phase space. Length is 5, with the first 3 being the spatial coordinates and the last 2 being the velocity coordinates.

orbit_quantities: tuple[tuple[int, str, str, str], ...] = ((0, 'x', '$x$', 'physical position x'), (1, 'y', '$y$', 'physical position y'), (2, 'z', '$z$', 'physical position z'), (3, 'v_par', '$v_\\parallel$', 'parallel velocity'), (4, 'v_perp', '$v_\\perp$', 'perpendicular velocity'), (5, 'weight', '$w$', 'marker weight'), (10, 'p_phi', '$p_\\phi$', 'canonical toroidal momentum (set by save_constants_of_motion)'))#: Marker columns saved as orbits, see orbit_quantities.

default_background = GyroMaxwellian2Dvperp( n=(1.0, None), u_para=(0.0, None), u_perp=(0.0, None), vth_para=(1.0, None), vth_perp=(1.0, None), equil=None, volume_form=True, uniform_on_disc=False, )#: Default kinetic background is a gyrotropic Maxwellian in $(v_\parallel, v_\perp)$.

default_n_cols = {'aux': 12, 'diagnostics': 3}#: Default number of buffer columns is 3 diagnostics (energy, magnetic moment, canonical toroidal momentum, see save_constants_of_motion()) and 12 auxiliary columns.

property mu_idx#: Index of the column in the marker array where the magnetic moment is stored.

__post_init__()[source]#: Retrieve the discrete equilibrium magnetic-field quantities ($|B_0|$, unit 1-form $\mathbf b_0$, Derham complex) needed to project marker velocities onto $v_\parallel, v_\perp$ and to evaluate diagnostics, and allocate the temporary FE coefficient vectors used for that.

property magn_bckgr#: Equilibrium fluid background carrying the magnetic field $\mathbf B_0$ with respect to which $v_\parallel, v_\perp$ are defined.

property absB0_h#: Discrete 0-form coefficients of $|B_0|$.

property unit_b1_h#: Discrete 1-form coefficients of the equilibrium field-aligned unit vector $\mathbf b_0 = \mathbf B_0/|B_0|$.

property epsilon#: Normalization parameter $\epsilon$ (from equation_params) entering the guiding-center equations of motion, e.g. the canonical toroidal momentum evaluation.

property derham#: Discrete Derham complex of the projected equilibrium.

property sampling_density#

Sampling density function as volume form, used to draw markers via inverse transform/rejection sampling and to compute their initial weights (see draw_markers()).

This is a GyroMaxwellian2Dvperp in $(v_\parallel, v_\perp)$, parametrized by the mean/thermal parallel and perpendicular velocities in loading_params. It is normalized to 1 in logical space (i.e. uniform in eta1, eta2, eta3) and already includes the polar-coordinate Jacobian factor $|v_\perp|$ (volume_form=True), further multiplied by 2 * eta1 if spatial is "disc".

Returns:

out (array-like) – The volume-form sampling density.
——-

s3(eta1, eta2, eta3, v_para, v_perp)[source]#

Sampling density function as 3-form, i.e. sampling_density() with the velocity-space ($|v_\perp|$) Jacobian factor divided back out, leaving a density that is a volume form in $\boldsymbol \eta$ only.

Parameters:

eta1 (array_like) – Logical evaluation points.
eta2 (array_like) – Logical evaluation points.
eta3 (array_like) – Logical evaluation points.
v_para (array_like) – Parallel and perpendicular velocity evaluation points.
v_perp (array_like) – Parallel and perpendicular velocity evaluation points.

Returns:

out (array-like) – The 3-form sampling density.
——-

s0(eta1, eta2, eta3, v_para, v_perp, flat_eval=False, remove_holes=True)[source]#

Sampling density function as 0-form, i.e. s3() pushed forward to a pointwise density by dividing out the spatial metric Jacobian determinant. This is the quantity stored in each marker’s s0 column and used to compute initial weights w0 = f_init / s0 / Np.

Parameters:

eta1 (array_like) – Logical evaluation points.
eta2 (array_like) – Logical evaluation points.
eta3 (array_like) – Logical evaluation points.
v_para (array_like) – Parallel and perpendicular velocity evaluation points.
v_perp (array_like) – Parallel and perpendicular velocity evaluation points.
flat_eval (bool) – If true, perform flat (marker) evaluation (etas must be same size 1D).
remove_holes (bool) – If True, holes are removed from the returned array. If False, holes are evaluated to -1.

Returns:

out (array-like) – The 0-form sampling density.
——-

draw_markers(sort: bool = True)[source]#

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:

Uniform distribution on the unit cube: $n^3(\eta) = 1$
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 Particles5Dvperp, 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.

On several processes, the markers are always sent to the process of their domain after drawing (mpi_sort_markers()). They are also sorted into boxes if sort is True (default) and sorting boxes were set up.

Parameters:: sort (Bool) – Whether to sort the markers into boxes after initial drawing (only if sorting boxes were set up).

save_constants_of_motion()[source]#

Calculate each marker’s guiding-center energy and canonical toroidal momentum (only the equilibrium magnetic field is considered) and assign them into the diagnostics columns of the marker array:

first_diagnostics_idx + 0: energy
first_diagnostics_idx + 1: magnetic moment (set once in draw_markers(), unchanged here since it is an adiabatic invariant)
first_diagnostics_idx + 2: canonical toroidal momentum

save_magnetic_energy(PBb)[source]#

Calculate the (time-dependent) magnetic field energy at each marker’s position and assign it into the energy diagnostics column (self.first_diagnostics_idx).

Parameters:: PBb (BlockVector) – Finite element coefficients of the time-dependent magnetic field, projected onto V0.

save_magnetic_background_energy()[source]#: Evaluate the equilibrium magnetic-moment energy $\mu_p |B_0(\boldsymbol \eta_p)|$ for each marker. The result is stored in the energy diagnostics column (self.first_diagnostics_idx).

save_magnetic_moment()[source]#: Calculate the magnetic moment of each marker and assign it into the magnetic-moment diagnostics column (self.first_diagnostics_idx + 1).

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

Bases: Particles

Particles in pure 3D configuration space $\boldsymbol \eta \in [0, 1]^3$, with no velocity space attached (vdim = 0) — each marker only carries a logical (curvilinear) position, used e.g. to represent a (massless) tracer or cold-plasma fluid density.

See Particles for the structure of the numpy marker array and the meaning of its columns.

Parameters:

name (str) – Name of particle species.
Np (int) – Number of particles.
bc (list) – Either ‘remove’, ‘reflect’, ‘periodic’ or ‘refill’ in each direction.
loading (str) – Drawing of markers; either ‘pseudo_random’, ‘sobol_standard’, ‘sobol_antithetic’, ‘external’ or ‘restart’.
**kwargs (dict) – Parameters for markers, see Particles.

vdim = 0#: Dimension of the velocity space, here 0 (no velocity coordinates).

coordinate_labels = ('$\\eta_1$', '$\\eta_2$', '$\\eta_3$')#: Labels for the coordinates in the phase space. Length is 3, with all being spatial coordinates.

orbit_quantities: tuple[tuple[int, str, str, str], ...] = ((0, 'x', '$x$', 'physical position x'), (1, 'y', '$y$', 'physical position y'), (2, 'z', '$z$', 'physical position z'), (3, 'weight', '$w$', 'marker weight'))#: Marker columns saved as orbits, see orbit_quantities.

default_background = ColdPlasma( n=(1.0, None), u1=(0.0, None), u2=(0.0, None), u3=(0.0, None), equil=None, uniform_on_disc=False, vth1=(0.0, None), vth2=(0.0, None), vth3=(0.0, None), )#: Default kinetic background is a cold-plasma (velocity-independent) density.

default_n_cols = {'aux': 5, 'diagnostics': 0}#: Default number of buffer columns reserved for diagnostics and auxiliary (pusher/free) use.

property mu_idx#: Index of the column in the marker array where the magnetic moment is stored.

__post_init__()[source]#: No additional setup is required for this class.

property sampling_density#: Sampling density function as volume form.

s0(eta1, eta2, eta3, flat_eval=False, remove_holes=True)[source]#

Sampling density function as 0 form, i.e. sampling_density() pushed forward to a pointwise (non-volume-form) density by dividing out the metric Jacobian determinant.

Parameters:

eta1 (array_like) – Logical evaluation points.
eta2 (array_like) – Logical evaluation points.
eta3 (array_like) – Logical evaluation points.
flat_eval (bool) – If true, perform flat (marker) evaluation (etas must be same size 1D).
remove_holes (bool) – If True, holes are removed from the returned array. If False, holes are evaluated to -1.

Returns:

out (array-like) – The 0-form sampling density.
——-

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

Bases: Particles

Particles for Smoothed Particle Hydrodynamics (SPH) models. The particle distribution itself lives in pure 3D configuration space $\boldsymbol \eta \in [0, 1]^3$, exactly as for Particles3D (sampling_density() and s0() depend only on $\boldsymbol \eta_p$).

Each marker additionally carries a Cartesian velocity $\mathbf v_p$ in its marker-array columns, but this is a per-particle helper quantity (e.g. the SPH velocity-field sample used by pushers and kernel-based reconstructions) rather than a coordinate of a sampled phase-space density.

See Particles for the structure of the numpy marker array and the meaning of its columns.

Parameters:

name (str) – Name of the particle species.
**params (dict) – Parameters for markers, see Particles.

vdim = 3#: Dimension of the per-marker Cartesian velocity attribute, here 3 (not a sampled coordinate, see class docstring).

coordinate_labels = ('$\\eta_1$', '$\\eta_2$', '$\\eta_3$', '$v_x$', '$v_y$', '$v_z$')#: Labels for the coordinates in the phase space. Length is 6, with the first 3 being the spatial coordinates and the last 3 being the velocity coordinates.

orbit_quantities: tuple[tuple[int, str, str, str], ...] = ((0, 'x', '$x$', 'physical position x'), (1, 'y', '$y$', 'physical position y'), (2, 'z', '$z$', 'physical position z'), (3, 'v1', '$v_x$', 'Cartesian velocity x'), (4, 'v2', '$v_y$', 'Cartesian velocity y'), (5, 'v3', '$v_z$', 'Cartesian velocity z'), (6, 'weight', '$w$', 'marker weight'))#: Marker columns saved as orbits, see orbit_quantities.

default_background = ConstantVelocity( ux=0.0, uy=0.0, uz=0.0, velocity_step_function_in_y=None, n=1.0, n1=0.0, density_profile=constant, upper_x=None, lower_x=None, upper_y=None, lower_y=None, p0=1.0, )#: Default fluid background is a spatially constant velocity field.

default_n_cols = {'aux': 24, 'diagnostics': 0}#: Default number of buffer columns reserved for diagnostics and auxiliary (pusher/free) use.

property mu_idx#: Index of the column in the marker array where the magnetic moment is stored.

__post_init__()[source]#: Attach the domain to the background (needed to evaluate it at marker positions). SPH does not support clone-based (tile-copied) MPI parallelization.

property sampling_density#

Sampling density function as volume form, used to draw markers via inverse transform/rejection sampling and to compute their initial weights (see draw_markers()).

This density is purely spatial: uniform (normalized to 1) if spatial is "uniform", or multiplied by the Jacobian factor 2 * eta1 if spatial is "disc".

s0(eta1, eta2, eta3, *v, flat_eval=False, remove_holes=True)[source]#

Sampling density function as 0 form, i.e. sampling_density() pushed forward to a pointwise (non-volume-form) density by dividing out the metric Jacobian determinant.

Parameters:

eta1 (array_like) – Logical evaluation points.
eta2 (array_like) – Logical evaluation points.
eta3 (array_like) – Logical evaluation points.
*v (array_like) – Accepted for a call signature compatible with generic phase-space evaluation, but unused (see sampling_density()).
flat_eval (bool) – If true, perform flat (marker) evaluation (etas must be same size 1D).
remove_holes (bool) – If True, holes are removed from the returned array. If False, holes are evaluated to -1.

Returns:

out (array-like) – The 0-form sampling density.
——-

Particle-to-grid accumulation#

Three classes deposit marker data on the FE spaces:

Accumulator assembles particle matrices (and optionally vectors), $\sum_p \Lambda^\alpha_{\mu, ijk}(\boldsymbol \eta_p) F_{p, \mu\nu} \Lambda^\alpha_{\nu, mno}(\boldsymbol \eta_p)$; operators[0] is the assembled matrix as a linear operator;
AccumulatorVector accumulates vectors only;
MarkerOperator is the matrix-free counterpart of Accumulator.operators[0]: a linear operator that applies the particle matrix to coefficients $s$ by evaluating their spline function $S_h$ (see SplineFunction) at the markers and depositing $F_p\, S_h(\boldsymbol \eta_p)$, in each iteration of a linear solver. The fillings $F_p$ are stored in the auxiliary marker columns by the accumulation kernel of the propagator (*_mf); the generic kernels marker_operator_symm_dot, marker_operator_asym_dot and marker_operator_rank1_dot apply symmetric (6 entries), antisymmetric (3 entries) and rank-one (4 entries) fillings, and marker_operator_pressure_dot the 3x3 block matrix of the pressure coupling (n_blocks=3). Assembling costs the products of all pairs of basis functions per marker, $\mathcal O((p+1)^6)$, once per step; the matrix-free application costs one spline evaluation and one deposition per marker, $\mathcal O((p+1)^3)$, per solver iteration. It is faster for 3D grids with higher degrees, expensive fillings or few solver iterations. All propagators that assemble particle matrices (except the disabled FaradayExtended) have the option matrix_free: CurrentCoupling5DCurlb, CurrentCoupling5DDensity, VlasovAmpereCoupling, EfieldWeightsCoupling, CurrentCoupling6DDensity, CurrentCoupling6DCurrent and PressureCoupling6D. It is the default in the 5D current couplings and in the pressure coupling, whose assembly is the most expensive (nine blocks).

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

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

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 (Kernel | PyccelKernel) – The accumulation kernel, e.g. charge_density_0form from struphy.pic.accumulation.kernels.charge_density_0form. On the CuPy backend, a kernel without CUDA version raises NotImplementedError.
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 live in struphy/pic/accumulation/kernels/, one folder per kernel (see struphy.pic.accumulation.kernels).

__init__(particles: Particles, space_id: str, kernel: Kernel | PyccelKernel, mass_ops: WeightedMassOperators, args_domain: DomainArguments, *, add_vector: bool = False, symmetry: str | None = None, filter_params: FilterParameters | None = None)[source]#

__call__(*optional_args)[source]#

Performs the accumulation into the matrix/vector by calling the chosen accumulation kernel.

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.

property particles#: Particle object.

property kernel: Kernel#: The accumulation kernel.

kernels() → tuple[Kernel, ...][source]#: The kernels this accumulator calls (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

show_accumulated_spline_field(mass_ops: WeightedMassOperators, eta_direction=0, component=0)[source]#

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: Kernel | PyccelKernel, mass_ops: WeightedMassOperators, args_domain: DomainArguments, filter_params: FilterParameters | None = None, n_vectors: int = 1)[source]#

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 (Kernel | PyccelKernel) – The accumulation kernel, e.g. charge_density_0form from struphy.pic.accumulation.kernels.charge_density_0form. On the CuPy backend, a kernel without CUDA version raises NotImplementedError.
derham (Derham) – Discrete FE spaces object.
args_domain (DomainArguments) – Mapping infos.
filter_params (FilterParameters) – Params for the accumulation filter.
n_vectors (int) – Number of vectors accumulated by one kernel call (in the order of the kernel arguments), e.g. 3 for the pressure coupling. Default=1.

__init__(particles: Particles, space_id: str, kernel: Kernel | PyccelKernel, mass_ops: WeightedMassOperators, args_domain: DomainArguments, filter_params: FilterParameters | None = None, n_vectors: int = 1)[source]#

__call__(*optional_args)[source]#

Performs the accumulation into the vector by calling the chosen accumulation kernel.

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.

property particles#: Particle object.

property kernel: Kernel#: The accumulation kernel.

kernels() → tuple[Kernel, ...][source]#: The kernels this accumulator calls (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

show_accumulated_spline_field(mass_ops, eta_direction=(True, False, False), save_L2=False)[source]#

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)

class struphy.pic.accumulation.particles_to_grid.ParticlesToGrid(pic_variable: PICVariable | SPHVariable | None = None, accum_space: Literal['H1', 'Hcurl', 'Hdiv', 'L2', 'H1vec'] | None = None, accum_kernel: Kernel | PyccelKernel | None = None)[source]#

Bases: object

Lightweight, serializable description of a particle-to-grid coupling (for example charge- or current deposition) into FEEC degrees of freedom.

A ParticlesToGrid does not perform any accumulation itself: it simply bundles the pieces needed to build an AccumulatorVector.

Parameters:

pic_variable (PICVariable | SPHVariable) – The kinetic variable whose markers (pic_variable.particles) are deposited on the grid.
accum_space ({"H1", "Hcurl", "Hdiv", "L2", "H1vec"}) – FEEC space identifier of the vector to accumulate into.
accum_kernel (Kernel | PyccelKernel) – Accumulation kernel matching accum_space, for example charge_density_0form from struphy.pic.accumulation.kernels.charge_density_0form.

Examples

>>> from struphy.pic.accumulation.kernels.charge_density_0form import charge_density_0form
>>> from struphy.pic.accumulation.particles_to_grid import ParticlesToGrid
>>> from struphy.propagators.poisson_solve import PoissonSolve
>>> rho = ParticlesToGrid(
...     kinetic_ions.var,
...     "H1",
...     charge_density_0form,
... )
>>> poisson = PoissonSolve(rho=rho, rho_coeffs=alpha**2 / epsilon)

__eq__(other)#: Return self==value.

__hash__ = None#

__init__(pic_variable: PICVariable | SPHVariable | None = None, accum_space: Literal['H1', 'Hcurl', 'Hdiv', 'L2', 'H1vec'] | None = None, accum_kernel: Kernel | PyccelKernel | None = None) → None#

__repr__()#: Return repr(self).

__weakref__#: list of weak references to the object (if defined)

class struphy.pic.accumulation.particles_to_grid.MarkerOperator(particles: Particles, space_id: str, kernel: Kernel | PyccelKernel, mass_ops: WeightedMassOperators, args_domain: DomainArguments, args: tuple = (), symmetry: str = 'symm', n_blocks: int = 1)[source]#

Bases: LinearOperator

Matrix-free particle matrix: the matrix-free counterpart of the matrices assembled by Accumulator.

For a vector-valued space $V^\alpha_h$ ($\alpha = 1, 2, v$) with the basis functions $\Lambda^\alpha_{\mu, ijk}$, and a 3x3 matrix $F_p$ per marker (e.g. of rank one, $F_p = c_p\, \mathbf e_p \mathbf e_p^\top$, or antisymmetric), the particle matrix is

\[A_{(\mu, ijk), (\nu, mno)} = \sum_p \Lambda^\alpha_{\mu, ijk}(\boldsymbol \eta_p)\, F_{p, \mu\nu}\, \Lambda^\alpha_{\nu, mno}(\boldsymbol \eta_p)\,.\]

Assembling $A$ costs, per marker, the products of all pairs of basis functions that do not vanish at $\boldsymbol \eta_p$ in each block. This operator is never assembled; with the spline function $S_h$ of the coefficients $s$ (see SplineFunction),

\[\sum_{\nu, mno} F_{p, \mu\nu}\, \Lambda^\alpha_{\nu, mno}(\boldsymbol \eta_p)\, s_{\nu, mno} = \big(F_p\, S_h(\boldsymbol \eta_p)\big)_\mu\,,\]

it is applied as

\[(A s)_{\mu, ijk} = \sum_p \Lambda^\alpha_{\mu, ijk}(\boldsymbol \eta_p)\, \big(F_p\, S_h(\boldsymbol \eta_p)\big)_\mu\,,\]

i.e. one evaluation of $S_h$ and one deposition per marker (by kernel, called through an AccumulatorVector without filter). The data of $F_p$ must be stored in the marker array beforehand, e.g. by the accumulation kernel of a vector. Like the assembled particle matrices, the operator acts on the coefficients with boundary conditions and polar extraction: with the extraction and boundary operators $\mathbb E$ and $\mathbb B$ of the space,

\[s \mapsto \mathbb B \mathbb E\, A\, \mathbb E^\top \mathbb B^\top s\,.\]

Parameters:

particles (Particles) – The markers.
space_id (str) – The FE space of the coefficients ("Hcurl", "Hdiv" or "H1vec").
kernel (Kernel) – The kernel (args_markers, args_derham, args_domain, vec1, vec2, vec3, s1, s2, s3, *args) that deposits $A s$ into vec, where s1, s2, s3 are the coefficient arrays of the three components of $s$.
mass_ops (WeightedMassOperators) – Weighted mass operators (for the de Rham sequence).
args_domain (DomainArguments) – Arguments of the mapping.
args (tuple) – Further arguments of the kernel.
symmetry (str) – "symm" ($F_p^\top = F_p$) or "asym" ($F_p^\top = -F_p$), for the transpose.
n_blocks (int) – Number of copies of the space in the domain (and codomain). For n_blocks > 1, the operator is the n_blocks x n_blocks block matrix $(A_{ab})_{a,b}$ on the product space, e.g. the pressure coupling with $A_{ab} = \sum_p v_{p,a} v_{p,b}\, \vec \Lambda_p F_p \vec \Lambda_p^\top$; then kernel takes the n_blocks output vectors followed by the n_blocks coefficient vectors (three arrays each), and symmetry refers to the whole block matrix.

__init__(particles: Particles, space_id: str, kernel: Kernel | PyccelKernel, mass_ops: WeightedMassOperators, args_domain: DomainArguments, args: tuple = (), symmetry: str = 'symm', n_blocks: int = 1)[source]#

property domain: VectorSpace#: The domain of the linear operator - an element of Vectorspace

property codomain: VectorSpace#: The codomain of the linear operator - an element of Vectorspace

property dtype#: The data type of the coefficients of the linear operator, upon convertion to matrix.

kernels() → tuple[Kernel, ...][source]#: The kernels this operator calls (the kernel that applies the matrix).

tosparse()[source]#

Assemble the global matrix of the linear operator as a scipy.sparse matrix.

Default implementation calling the generic LinearOperator.toarray with is_sparse=True; see there for cost and parallel behavior. Subclasses with an explicit matrix representation should override it.

Parameters:: format (str) – One of ‘csr’ (default), ‘csc’, ‘bsr’, ‘lil’, ‘dok’, ‘coo’ or ‘dia’.
Returns:: The global matrix of shape self.shape, identical on all ranks.
Return type:: scipy.sparse matrix

toarray()[source]#

Assemble the global matrix of the linear operator column by column.

Column j is computed as self.dot(e_j), where e_j is the j-th canonical basis vector of the domain, in the global numbering of Vector.toarray() (C-ordered within each StencilVector, blocks concatenated). Since only dot is used, this works for any linear operator, including matrix-free ones, whose domain is a StencilVectorSpace or a (possibly nested) BlockVectorSpace thereof.

The cost is one call to dot per global degree of freedom of the domain, hence this default is meant for testing and small problems. Subclasses with an explicit matrix representation should override it.

In parallel, all ranks call dot collectively, each rank fills the rows it owns, and every rank receives the full matrix.

Parameters:

out (numpy.ndarray, optional) – Dense array of shape self.shape into which the result is written in place. Must be None if is_sparse is True.
is_sparse (bool) – If True, return a scipy.sparse matrix, otherwise a dense numpy.ndarray.
format (str) – Sparse format, only used if is_sparse is True: one of ‘csr’ (default), ‘csc’, ‘bsr’, ‘lil’, ‘dok’, ‘coo’ or ‘dia’.

Returns:

The global matrix of shape self.shape, identical on all ranks.

Return type:

numpy.ndarray or scipy.sparse matrix

transpose(conjugate=False)[source]#

Transpose the LinearOperator .

If conjugate is True, return the Hermitian transpose.

dot(v: Vector, out: Vector | None = None) → Vector[source]#

Apply the LinearOperator self to the Vector v.

The result is written to the Vector out, if provided.

Parameters:

v (Vector) – The vector to which the linear operator (self) is applied. It must belong to the domain of self.
out (Vector) – The vector in which the result of the operation is stored. It must belong to the codomain of self. If out is None, a new vector is created and returned.

Returns:

The result of the operation. If out is None, a new vector is returned. Otherwise, the result is stored in out and out is returned.

Return type:

Vector

Accumulation kernels#

Accumulation kernels, one folder per kernel (see CUDA_STRATEGY.md).

Each folder <name> holds the pyccel kernel <name>_kernels.py and, once ported, its CUDA version <name>_cuda.cu. The folder’s __init__.py declares the kernel, a cunumpy.kernels.Kernel that runs the version of the active backend:

from struphy.pic.accumulation.kernels.charge_density_0form import charge_density_0form

struphy.pic.accumulation.kernels.vlasov_maxwell.vlasov_maxwell_kernels.vlasov_maxwell(args_markers: MarkerArguments, args_derham: DerhamArguments, args_domain: DomainArguments, mat11: float[:, :, :, :, :, :], mat12: float[:, :, :, :, :, :], mat13: float[:, :, :, :, :, :], mat22: float[:, :, :, :, :, :], mat23: float[:, :, :, :, :, :], mat33: float[:, :, :, :, :, :], vec1: float[:, :, :], vec2: float[:, :, :], vec3: float[:, :, :])[source]#: 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.

Pusher class#

Marker evaluation setup#

Propagator.init_kernels and Propagator.eval_kernels are tuples of KernelSetup instances. Each setup names the callable, its additional arguments, its destination marker columns, and the evaluation weights alpha. Initialization kernels run once at the start of each push. Evaluation kernels run before every stage/iteration, after sorting particles using their spatial alpha weights.

For example, register a vector evaluation with three explicit output columns:

from struphy.pic.pushing.kernel_setup import KernelSetup
from struphy.pic.pushing.kernels.unit_b_1form import unit_b_1form

self.add_init_kernel(
    KernelSetup(
        kernel=unit_b_1form,
        args=(self.derham.args_derham, b1, b2, b3),
        output_indices=(first_free_idx, first_free_idx + 1, first_free_idx + 2),
    )
)

Output indices are absolute marker columns in component order. Use None to skip a component, for example (20, None, 24). Scalars use a one-element tuple; tensors use nine destinations in row-major order. This replaces the previous column_nr and comps arguments. Kernels receive an integer array of output indices, with -1 representing a skipped component.

For an evaluation during iteration, pass a setup with the desired alpha to self.add_eval_kernel. A scalar weight applies to all six phase-space coordinates. A tuple supplies three to six weights; omitted velocity weights are zero. Initialization setups must use alpha=0 (the default).

class struphy.pic.pushing.kernel_setup.KernelSetup(*, kernel: Kernel | Callable, args: tuple = (), output_indices: tuple[int | None, ...], alpha: float | tuple[float, ...] = 0.0)[source]#

A marker evaluation kernel and all of its configuration.

output_indices gives the destination marker column for each result component: one index for a scalar, three for a vector, or nine in row-major order for a tensor. Use None to skip a component. For example, (20, None, 24) writes vector components 0 and 2 to columns 20 and 24.

alpha selects the evaluation state, coordinate by coordinate: alpha * current + (1 - alpha) * initial. A scalar applies to all six phase-space coordinates; a tuple supplies three to six weights explicitly, with omitted velocity weights set to zero. Initialization kernels use zero weights. The arrays needed by compiled kernels are prepared once; args retains references to mutable field data.

Kernels follow the signature kernel(alpha, output_indices, args_markers, args_domain, *args). At this boundary, a skipped output is represented by the integer -1.

kernels() → tuple[Kernel, ...][source]#: The kernel of this setup, as a one-element tuple.

property name: str#: Kernel name used in profiling and validation messages.

property sorting_alpha#: Spatial evaluation weights used for MPI sorting before execution.

validate_outputs(n_cols: int)[source]#: Check that every destination fits in the marker array.

evaluate(args_markers, args_domain)[source]#: Evaluate into the configured marker columns.

class struphy.pic.pushing.pusher.Pusher(particles: Particles, kernel: PyccelKernel | Kernel, args_kernel: tuple, args_domain: DomainArguments, pushes_eta: bool, *, alpha_in_kernel: float | int | tuple | list, init_kernels: tuple[KernelSetup, ...] = (), eval_kernels: tuple[KernelSetup, ...] = (), n_stages: int = 1, maxiter: int = 1, tol: float = 1e-08, mpi_sort: str | None = None, local_eval_only: bool = False)[source]#

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 an accelerated kernel from struphy.pic.pushing.kernels, e.g. push_eta_stage from struphy.pic.pushing.kernels.push_eta_stage.

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.

MPI sorting is skipped whenever it is known to be a no-op: markers are assumed to be sorted according to the domain decomposition on entry, and they stay sorted for any $\alpha$ until the pusher kernel has moved them. Pushers with pushes_eta=False hence never sort. Pushers with pushes_eta=True always sort at least once after the last stage, such that the above assumption holds for the next pusher.

Parameters:

particles (Particles) – Particles object holding the markers to push.
kernel (PyccelKernel | Kernel) – The pusher kernel. A Kernel also holds its CUDA version; on the CuPy backend, a kernel without CUDA version raises NotImplementedError.
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 (tuple[KernelSetup, ...]) – Evaluations at the initial state, executed once per push in tuple order. Each setup specifies the kernel, arguments, and output marker indices.
eval_kernels (tuple[KernelSetup, ...]) – Evaluations before each pusher stage/iteration. Each setup’s alpha weights determine the evaluation state and preceding MPI sort.
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 (only allowed for pushes_eta=False, becomes “last” otherwise). * each : sort markers after each stage. * last : sort markers after last stage.
pushes_eta (bool) – Whether the kernel updates the marker positions $\boldsymbol \eta_p$. If False, no MPI sorting is performed at all.
local_eval_only (bool) – Set to True if the kernel does not evaluate distributed splines, i.e. it only calls metric coefficients or equilibrium quantities, which are available on every process. Markers then need not be on the right process during the stages; they are sorted only once after the last stage (mpi_sort=”last”).

__call__(dt: float)[source]#: Applies the chosen pusher kernel by a time step dt, applies kinetic boundary conditions and performs MPI sorting.

property particles#: Particle object.

property kernel: Kernel#: The pusher kernel for the active backend (pyccel or CUDA).

kernels() → tuple[Kernel, ...][source]#: The kernels this pusher calls: the pusher kernel and those of its init and eval kernel setups.

property init_kernels: tuple[KernelSetup, ...]#: Ordered setups for evaluations at the initial state.

property eval_kernels: tuple[KernelSetup, ...]#: Ordered setups for evaluations before each pusher stage/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.

property pushes_eta#: Whether the kernel updates the marker positions.

property local_eval_only#: Whether the kernel needs no distributed spline evaluations (sorting only after last stage).

__weakref__#: list of weak references to the object (if defined)

property mpi_sort#: When to do MPI sorting: * None : no sorting at all (only for pushes_eta=False). * each : sort markers after each stage. * last : sort markers after last stage.

Pusher kernels#

struphy.pic.pushing.kernels.push_v_with_efield.push_v_with_efield_kernels.push_v_with_efield(dt: float, stage: int, args_markers: MarkerArguments, args_domain: DomainArguments, args_derham: DerhamArguments, e1_1: float[:, :, :], e1_2: float[:, :, :], e1_3: float[:, :, :], const_factor: float)[source]#

Updates particle velocities as

\[\frac{\mathbf v^{n+1} - \mathbf v^n}{\Delta t} = c \, \bar{DF}^{-\top} (\mathbb L^1)^\top \mathbf e\]

where $\mathbf e \in \mathbb R^{N_1}$ are given FE coefficients of the 1-form spline field and $c \in \mathbb R$ is some constant (const_factor).

Parameters:

e1_1 (ndarray[float]) – 3d array of FE coeffs of E-field as 1-form.
e1_2 (ndarray[float]) – 3d array of FE coeffs of E-field as 1-form.
e1_3 (ndarray[float]) – 3d array of FE coeffs of E-field as 1-form.
const_factor (float) – A constant (usuallly related to the charge-to-mass ratio).

Table of Contents

Previous topic

Next topic

PIC#

Particle base class#

Marker array layout#

Particle subclasses#

Particle-to-grid accumulation#

Accumulation kernels#

Pusher class#

Marker evaluation setup#

Pusher kernels#