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 numpy array. In markers[ip, j] The row index ip refers to a specific particle, the column index j to its attributes. The columns are indexed 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_rows, 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#: Total umber of columns in markers array. The last 2 columns refer to box number and particle ID, respectively.

property n_rows#: Total number of rows in markers array.

property 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_rows 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 numpy array holding the marker information, including holes. The i-th row holds the i-th marker info.

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

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

property holes#: Array of booleans stating if an entry in the markers array is a hole.

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

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

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

property 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 holding the marker positions in logical space. The i-th row holds the i-th marker info.

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

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

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

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

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

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

property 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#: Collection of mandatory arguments for pusher 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.

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

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

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. Numpy’s histogramdd is used, following the algorithm outlined in Particle binning.

Parameters:

components (tuple[bool]) – List of length 3 + vdim; an entry is True if the direction in phase space is to be binned.
bin_edges (tuple[array]) – List of bin edges (resolution) having the length of True entries in components.
output_quantity (BinningOutput) – String literal used to determine weights in binning and the type of output
divide_by_jac (bool) – Whether to divide the weights by the Jacobian determinant for binning (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.

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 (rows 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 row index.

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

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

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 eval_kernels_sph.sph_mean_velocity_coeffs (via a Pyccel kernel) to assemble mean-velocity coefficients into the markers array, then calls eval_sph() for each velocity component.

eval_div_viscosity(eta1, eta2, eta3, h1, h2, h3, kernel_type='gaussian_1d', mu: float = 1.0, fast=True) → tuple[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). It then evaluates the necessary derivatives via eval_sph() and sums contributions to produce the three 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)

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': 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]#: 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.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.

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

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

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

Particel-to-grid accumulation#

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

class struphy.pic.accumulation.particles_to_grid.Accumulator(particles: Particles, space_id: str, kernel: PyccelKernel, mass_ops: WeightedMassOperators, args_domain: DomainArguments, *, add_vector: bool = False, symmetry: str | None = None, filter_params: FilterParameters | None = None)[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 (pyccelized function) – The accumulation kernel.
derham (Derham) – Discrete FE spaces object.
args_domain (DomainArguments) – Mapping infos.
add_vector (bool) – True if, additionally to a matrix, a vector in the same space is to be accumulated. Default=False.
symmetry (str) – In case of space_id=Hcurl/Hdiv, the symmetry property of the block matrix: diag, asym, symm, pressure or None (=full matrix, default)
filter_params (dict) – Params for the accumulation filter: use_filter(string, either `three_point or `fourier), repeat(int), alpha(float) and modes(list with int).

Note

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

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

__call__(*optional_args, **args_control)[source]#

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

Parameters:

particles (Particles) – Particles object holding the markers information in format particles.markers.shape == (n_markers, :).
optional_args (any) – Additional arguments to be passed to the accumulator kernel, besides the mandatory arguments which are prepared automatically (spline bases info, mapping info, data arrays). Examples would be parameters for a background kinetic distribution or spline coefficients of a background magnetic field. Entries must be pyccel-conform types.
args_control (any) – Keyword arguments for an analytical control variate correction in the accumulation step. Possible keywords are ‘control_vec’ for a vector correction or ‘control_mat’ for a matrix correction. Values are a 1d (vector) or 2d (matrix) list with callables or xp.ndarrays used for the correction.

property particles#: Particle object.

property kernel: PyccelKernel#: The accumulation kernel.

property derham#: Discrete Derham complex on the logical unit cube.

property args_domain#: Mapping info for evaluating metric coefficients.

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

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

property symmetry#: Symmetry of the accumulation matrix (diagonal, symmetric, asymmetric, etc.).

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

property vectors#: List of Stencil-/Block-/PolarVectors of the accumulator.

property accfilter#: Callable filters

init_control_variate(mass_ops)[source]#: Set up the use of noise reduction by control variate.

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: PyccelKernel, mass_ops: WeightedMassOperators, args_domain: DomainArguments, filter_params: FilterParameters | None = None)[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 (pyccelized function) – The accumulation kernel.
derham (Derham) – Discrete FE spaces object.
args_domain (DomainArguments) – Mapping infos.

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

__call__(*optional_args, **args_control)[source]#

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

Parameters:

optional_args (any) – Additional arguments to be passed to the accumulator kernel, besides the mandatory arguments which are prepared automatically (spline bases info, mapping info, data arrays). Examples would be parameters for a background kinetic distribution or spline coefficients of a background magnetic field. Entries must be pyccel-conform types.
args_control (any) – Keyword arguments for an analytical control variate correction in the accumulation step. Possible keywords are ‘control_vec’ for a vector correction or ‘control_mat’ for a matrix correction. Values are a 1d (vector) or 2d (matrix) list with callables or xp.ndarrays used for the correction.

property particles#: Particle object.

property kernel: PyccelKernel#: The accumulation kernel.

property derham#: Discrete Derham complex on the logical unit cube.

property args_domain#: Mapping arguments.

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

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

property vectors#: List of Stencil-/Block-/PolarVectors of the accumulator.

property accfilter#: Callable filters

init_control_variate(mass_ops)[source]#: Set up the use of noise reduction by control variate.

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: 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 (PyccelKernel) – Pyccelized accumulation kernel matching accum_space, for example PyccelKernel(accum_kernels.charge_density_0form).

Examples

>>> from struphy.pic.accumulation import accum_kernels
>>> from struphy.pic.accumulation.particles_to_grid import ParticlesToGrid
>>> from struphy.propagators.poisson_solve import PoissonSolve
>>> from cunumpy import PyccelKernel
>>> rho = ParticlesToGrid(
...     kinetic_ions.var,
...     "H1",
...     PyccelKernel(accum_kernels.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: PyccelKernel | None = None) → None#

__repr__()#: Return repr(self).

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

Accumulation kernels#

struphy.pic.accumulation.accum_kernels.cc_lin_mhd_6d_1()#

Accumulates into V1 with the filling functions

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

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

Parameters:

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

Note

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

struphy.pic.accumulation.accum_kernels.cc_lin_mhd_6d_2()#

Accumulates into V1 with the filling functions

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

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

Parameters:

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

Note

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

struphy.pic.accumulation.accum_kernels.charge_density_0form()#

Kernel for AccumulatorVector into V0 with filling function

\[B_p = w_p \,,\]

where $w_p$ is the marker weight.

struphy.pic.accumulation.accum_kernels.div_u_weak_1form()#

Kernel for AccumulatorVector into V1 with filling function

\[\mathbf{B}_p = \frac{w_p}{n_p} \mathbf{v}_p \,,\]

where $w_p$ is the marker weight, $\mathbf{v}_p$ the marker velocity and $n_p$ the (previously accumulated) density evaluated at the marker position, stored in column args_markers.first_free_idx of the markers array.

struphy.pic.accumulation.accum_kernels.linear_vlasov_ampere()#

Accumulates into V1 with the filling functions

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

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

Note

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

struphy.pic.accumulation.accum_kernels.pc_lin_mhd_6d()#: Accumulates into V1 with the filling functions

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

Note

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

struphy.pic.accumulation.accum_kernels.pc_lin_mhd_6d_full()#: Accumulates into V1 with the filling functions

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

Note

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

struphy.pic.accumulation.accum_kernels.vlasov_maxwell()#: Accumulates into V1 with the filling functions

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

Note

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

struphy.pic.accumulation.accum_kernels_gc.cc_lin_mhd_5d_D()#

Accumulation kernel for the propagator CurrentCoupling5DDensity.

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

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

Parameters:

epsilon (float) – scaling factor.
b2_1 (array[float]) – FE coefficients c_ijk of the magnetic field as a 2-form.
b2_2 (array[float]) – FE coefficients c_ijk of the magnetic field as a 2-form.
b2_3 (array[float]) – FE coefficients c_ijk of the magnetic field as a 2-form.
norm_b11 (array[float]) – FE coefficients c_ijk of the unit magnetic field as a 1-form.
norm_b12 (array[float]) – FE coefficients c_ijk of the unit magnetic field as a 1-form.
norm_b12 – FE coefficients c_ijk of the unit magnetic field as a 1-form.
curl_norm_b1 (array[float]) – FE coefficients c_ijk of the curl of the unit magnetic field as a 2-form.
curl_norm_b2 (array[float]) – FE coefficients c_ijk of the curl of the unit magnetic field as a 2-form.
curl_norm_b3 (array[float]) – FE coefficients c_ijk of the curl of the unit magnetic field as a 2-form.

Note

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

struphy.pic.accumulation.accum_kernels_gc.cc_lin_mhd_5d_M()#

Accumulation kernel for the propagator ShearAlfvenCurrentCoupling5D and MagnetosonicCurrentCoupling5D.

Accumulates 2-form vector with the filling functions:

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

Parameters:

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

Note

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

struphy.pic.accumulation.accum_kernels_gc.cc_lin_mhd_5d_curlb()#

Accumulation kernel for the propagator CurrentCoupling5DCurlb.

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

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

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

struphy.pic.accumulation.accum_kernels_gc.cc_lin_mhd_5d_gradB()#

Accumulation kernel for the propagator CurrentCoupling5DGradB.

Accumulates math:alpha -form vector with the filling functions

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

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

Parameters:

b1 (array[float]) – FE coefficients c_ijk of the magnetic field as a 2-form.
b2 (array[float]) – FE coefficients c_ijk of the magnetic field as a 2-form.
b3 (array[float]) – FE coefficients c_ijk of the magnetic field as a 2-form.
norm_b11 (array[float]) – FE coefficients c_ijk of the normalized magnetic field as a 1-form.
norm_b12 (array[float]) – FE coefficients c_ijk of the normalized magnetic field as a 1-form.
norm_b13 (array[float]) – FE coefficients c_ijk of the normalized magnetic field as a 1-form.
curl_norm_b1 (array[float]) – FE coefficients c_ijk of the curl of normalized magnetic field as a 2-form.
curl_norm_b2 (array[float]) – FE coefficients c_ijk of the curl of normalized magnetic field as a 2-form.
curl_norm_b3 (array[float]) – FE coefficients c_ijk of the curl of normalized magnetic field as a 2-form.
grad_PB1 (array[float]) – FE coefficients c_ijk of gradient of parallel magnetic field as a 1-form.
grad_PB2 (array[float]) – FE coefficients c_ijk of gradient of parallel magnetic field as a 1-form.
grad_PB3 (array[float]) – FE coefficients c_ijk of gradient of parallel magnetic field as a 1-form.
grad_PBeq1 (array[float]) – FE coefficients c_ijk of gradient of equilibrium parallel magnetic field as a 1-form; added to grad_PB for all u-spaces.
grad_PBeq2 (array[float]) – FE coefficients c_ijk of gradient of equilibrium parallel magnetic field as a 1-form; added to grad_PB for all u-spaces.
grad_PBeq3 (array[float]) – FE coefficients c_ijk of gradient of equilibrium parallel magnetic field as a 1-form; added to grad_PB for all u-spaces.

Note

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

struphy.pic.accumulation.accum_kernels_gc.cc_lin_mhd_5d_gradB_dg()#: TODO

struphy.pic.accumulation.accum_kernels_gc.cc_lin_mhd_5d_gradB_dg_init()#: TODO

struphy.pic.accumulation.accum_kernels_gc.gc_density_0form()#: Kernel for AccumulatorVector into V0 with the filling

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

struphy.pic.accumulation.accum_kernels_gc.gc_mag_density_0form()#: Kernel for AccumulatorVector into V0 with the filling

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

Pusher class#

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

self.add_init_kernel(
    KernelSetup(
        kernel=eval_kernels_gc.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: 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.

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, 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 accelerated pusher_kernels or pusher_kernels_gc for guiding-center models.

Notes

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

Spline/geometry evaluations at $\boldsymbol \eta^n_p$ can be be done via init_kernels.
Pusher kernel and eval_kernels can perform evaluations at arbitrary weighted averages $\eta_{p,i} = \alpha_i \eta_{p,i}^{n+1,k} + (1 - \alpha_i) \eta_{p,i}^n$, for $i=1,2,3$.
MPI sorting is done automatically before kernel calls according to the specified values $\alpha_i$ for each kernel.

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 (pyccelized function) – The pusher kernel.
args_kernel (tuple) – Optional arguments passed to the kernel.
args_domain (DomainArguments) – Mapping infos.
alpha_in_kernel (float | int | tuple | list) – For i=0,1,2, the spline/geometry evaluations in kernel are at alpha[i]*markers[:, i] + (1 - alpha[i])*markers[:, buffer_idx + i]. If float or int or then alpha = (alpha, alpha, alpha). alpha must be between 0 and 1. alpha[i]=0 means that evaluation is at the initial positions (time n), stored at markers[:, buffer_idx + i].
init_kernels (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#: The pyccelized pusher kernel.

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.

Table of Contents

Previous topic

Next topic

PIC#

Particle base class#

Particle subclasses#

Particel-to-grid accumulation#

Accumulation kernels#

Pusher class#

Marker evaluation setup#