Base modules

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

Bases: object

Base class for particle species.

The marker information is stored in a 2D numpy array, see Tutorial on PIC data structures.

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 indices; see attributes first_diagnostics_idx, first_pusher_idx and first_free_idx below

• -2: number of the sorting box the particle is in

• -1: particle ID

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.

• mpi_dims_mask (list | tuple of bool) – 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.

• boxes_per_dim (tuple) – Number of boxes in each logical direction (n_eta1, n_eta2, n_eta3).

• box_bufsize (float) – Between 0 and 1, relative buffer size for box array (default = 0.25).

• type (str) – Either ‘full_f’ (default), ‘delta_f’ or ‘sph’.

• name (str) – Name of particle species.

• loading_params (LoadingParameters) – Parameterts for particle loading.

• weights_params (WeightsParameters) – Parameters for particle weights.

• boundary_params (BoundaryParameters) – Parameters for particle boundary conditions.

• 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, …)

• verbose (bool) – Show some more Particle info.

abstract classmethod default_background()

The default background (of type Maxwellian).

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

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

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

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

abstract property vdim

Dimension of the velocity space.

abstract property n_cols_diagnostics

Number of columns for storing diagnostics for each marker.

abstract property n_cols_aux

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

property first_diagnostics_idx

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

property first_pusher_idx

Starting index for storing initial conditions for a Pusher call.

property n_cols_pusher

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

property first_shift_idx

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

property n_cols_shift

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

property residual_idx

Column for storing the residual in iterative pushers.

property first_free_idx

First index for storing auxiliary quantities for each particle.

property n_cols

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

property n_rows

Total number of rows in markers array.

property kinds

Name of the class.

property name

Name of the kinetic species in DATA container.

property type

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

Type:

Particle type

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

Type of particle loading.

property bc

List of particle boundary conditions in each direction.

property bc_refill

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

property bc_sph

List of boundary conditions for sph evaluation in each direction.

property Np

Total number of markers/particles, from user input.

property Np_per_clone

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

property ppc

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

property ppb

Particles per sorting box.

property bufsize

Relative size of buffer in markers array.

property mpi_comm

MPI communicator.

property mpi_size

Number of MPI processes.

property mpi_rank

Rank of current process.

property clone_config

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

property num_clones

Total number of clones.

property clone_id

Clone id of current process.

property background: KineticBackground

Kinetic background.

property perturbations: dict[str, Perturbation]

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

property loading_params: LoadingParameters

property weights_params: WeightsParameters

property boundary_params: BoundaryParameters

Parameters for marker loading.

property reject_weights

Whether to reect weights below threshold.

property threshold

Threshold for rejecting weights.

property boxes_per_dim

Tuple, number of sorting boxes per dimension.

property verbose

Show some more particles info.

property equation_params

Parameters appearing in model equation due to Struphy normalization.

property initial_condition: KineticBackground

Kinetic initial condition

property f_init

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

property u_init

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

property f0: Maxwellian

property control_variate

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

property domain_array

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

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

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

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

property mpi_dims_mask

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

property nprocs

Number of MPI processes in each dimension.

property n_mks_load

Array of number of markers on each process at loading stage

property markers

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

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

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

property holes

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

property ghost_particles

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

property markers_wo_holes

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

property markers_wo_holes_and_ghost

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

property domain

From struphy.geometry.domains.

property equil

From struphy.fields_background.equils.

property projected_equil

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

property lost_markers

Array containing the last infos of removed markers

property n_lost_markers

Number of removed particles.

property index

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

property valid_mks

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

update_valid_mks()

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

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

property weights0

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

property marker_ids

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

property is_volume_form

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

Type:

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

property spatial

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

property f_coords_index

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

property f_jacobian_coords_index

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

property f_coords

Coordinates of the distribution function.

property args_markers

Collection of mandatory arguments for pusher kernels.

property f_jacobian_coords

Coordinates of the velocity jacobian determinant of the distribution fuction.

property sorting_boxes

property tesselation

Tesselation of the current process domain.

classmethod ker_dct()

Available smoothing kernels, numbers must be multiplies of 100.

draw_markers(sort: bool = True, verbose: bool = True)

Drawing markers according to the volume density \(s^\textrm{vol}_{\textnormal{in}}\). In Struphy, the initial marker distribution \(s^\textrm{vol}_{\textnormal{in}}\) is always of the form

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

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

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

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

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

The following choices are available in Struphy:

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

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

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

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

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

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

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

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

Then,

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

So then,

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

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

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)

• verbose (bool) – Show info on screen.

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

Sorts markers according to MPI domain decomposition.

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

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

Parameters:

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

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

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

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

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

Computes the initial weights

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

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

Parameters:

• bckgr_params (dict) – Kinetic background parameters.

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

update_weights()

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

reset_marker_ids()

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

binning(components: tuple[bool], bin_edges: tuple[ndarray], divide_by_jac: bool = True)

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

Parameters:

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

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

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

Returns:

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

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

show_distribution_function(components, bin_edges)

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

Parameters:

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

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

apply_kinetic_bc(newton=False)

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

Parameters:

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

particle_refilling()

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

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

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

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

gyro_transfer(outside_inds)

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

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

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

Then changes the direction of the perpendicular velocity,

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

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

Parameters:

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

Returns:

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

Return type:

xp.array (bool)

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

Bases: object

Boxes used for the sorting of the particles.

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

Parameters:

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

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

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

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

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

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

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

• comm (Intracomm) – MPI communicator or None.

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

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

property nx

property ny

property nz

property comm

property box_index

property boxes

property neighbours

property communicate

property is_domain_boundary

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

property bc_sph

List of boundary conditions for sph evaluation in each direction.

property bc_sph_index_shifts

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

put_particles_in_boxes()

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

check_and_assign_particles_to_boxes()

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

do_sort(use_numpy_argsort=False)

Assign the particles to boxes and then sort them.

remove_ghost_particles()

prepare_ghost_particles()

Markers for boundary conditions and MPI communication.

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

determine_markers_in_box(list_boxes)

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

get_destinations_box()

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

self_communication_boxes()

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

communicate_boxes(verbose=False)

sendrecv_all_to_all_boxes()

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

sendrecv_markers_boxes()

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

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

Density function as 0-form.

Parameters:

• eta1 (array_like) – Logical evaluation points (flat or meshgrid evaluation).

• eta2 (array_like) – Logical evaluation points (flat or meshgrid evaluation).

• eta3 (array_like) – Logical evaluation points (flat or meshgrid evaluation).

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

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

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

• kernel_type (str) – Name of the smoothing kernel to be used.

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

• fast (bool) – True: box-based evaluation, False: naive evaluation.

Returns:

• out (array-like) – Same size as eta1.

• ——-

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

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

\[b(\boldsymbol \eta_i) = \frac 1N \sum_k \beta_k W_h(\boldsymbol \eta_i - \boldsymbol \eta_k)\,.\]

The coefficients \(\beta_k\) must be stored at self.markers[k, index]. The possible choices for \(W_h\) are listed in sph_smoothing_kernels and in ker_dct().

Parameters:

• eta1 (array_like) – Logical evaluation points.

• eta2 (array_like) – Logical evaluation points.

• eta3 (array_like) – Logical evaluation points.

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

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

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

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

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

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

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

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

update_holes()

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

update_ghost_particles()

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

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

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

Parameters:

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

Returns:

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

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

sendrecv_get_destinations(send_inds)

Determine to which process particles have to be sent.

Parameters:

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

Returns:

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

Return type:

array[int]

sendrecv_all_to_all(send_info)

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

Parameters:

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

Returns:

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

Return type:

array[int]

sendrecv_markers(recv_info, hole_inds_after_send)

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

Parameters:

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

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

class struphy.pic.base.Tesselation(tiles_pb: int | float, *, comm: Intracomm | None = None, domain_array: ndarray | None = None, sorting_boxes: SortingBoxes | None = None)

Bases: object

Make a tesselation of the simulation domain into tiles of equal size.

Parameters:

• tiles_pb (int) – Number of equally sized tiles per box defined in sorting boxes (there is 1 box if sorting_boxes=None). This is equal to particels per box (ppb) when used for SPH markers.

• comm (Intracomm) – MPI communicator.

• domain_array (xp.ndarray) – A 2d array[float] of shape (comm.Get_size(), 9) holding info on the domain decomposition.

• sorting_boxes (Particles.SortingBoxes) – Box info for SPH evaluations.

get_tiles()

Compute tesselation of a single sorting box.

draw_markers()

Draw markers on the tile midpoints.

cell_averages(fun, n_quad=None)

Compute cell averages of fun over all tiles on current process.

Parameters:

fun (callable) – Some callable function.

property tiles_pb

Number of equally sized tiles per sorting box.

property n_tiles

Total number of tiles on current process.

property nt_per_dim

3-list of number of equally sized tiles per sorting box per direction.

property starts

3-list of domain starts (left boundaries) on current process.

property ends

3-list of domain ends (right boundaries) on current process.

property tile_breaks

3-list of tile break points within the single sorting box [0.0, sorting_box_width], in each direction.

property tile_midpoints

3-list of tile midpoints within the single sorting box [0.0, sorting_box_width], in each direction.

property tile_volume

Volume of a single tile.

property tile_quad_pts

3-list of quadrature points (GL) within a single tile, in each direction.

property tile_quad_wts

3-list of quadrature weights (GL) within a single tile, in each direction.

property rank

Current process rank.

property boxes_per_dim

Sorting boxes per direction.

property box_widths

3-list of sorting box widths in each direction.

property dims_mask

Boolean array of size 3; entry is True if direction participates in tesselation.

class struphy.pic.particles.Particles6D(**kwargs)

Bases: Particles

A class for initializing particles in models that use the full 6D phase space.

The numpy marker array is as follows:

index	0 | 1 | 2 |	3 | 4 | 5 |	6	7	8	>=9
value	position (eta)	velocities	weight	s0	w0	buffer

classmethod default_background()

The default background (of type Maxwellian).

property vdim

Dimension of the velocity space.

property n_cols_diagnostics

Number of the diagnostics columns.

property n_cols_aux

Number of the auxiliary columns.

property coords

Coordinates of the Particles6D, \((v_1, v_2, v_3)\).

svol(eta1, eta2, eta3, *v)

Sampling density function as volume form.

Parameters:

• eta1 (array_like) – Logical evaluation points.

• eta2 (array_like) – Logical evaluation points.

• eta3 (array_like) – Logical evaluation points.

• *v (array_like) – Velocity evaluation points.

Returns:

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

• ——-

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

Sampling density function as 0 form.

Parameters:

• eta1 (array_like) – Logical evaluation points.

• eta2 (array_like) – Logical evaluation points.

• eta3 (array_like) – Logical evaluation points.

• *v (array_like) – 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()

Calculate each markers’ guiding center constants of motions and assign them into diagnostics columns of marker array:

diagnostics index	0 | 1 | 2 |	3 |	4 |	5 |	6 |
value	guiding_center	energy	magn. moment	can. momentum	para. velocity

Only equilibrium magnetic field is considered.

class struphy.pic.particles.DeltaFParticles6D(**kwargs)

Bases: Particles6D

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

classmethod default_background()

The default background (of type Maxwellian).

set_n_to_zero(background: Maxwellian | SumKineticBackground)

class struphy.pic.particles.Particles5D(projected_equil: ProjectedFluidEquilibriumWithB, **kwargs)

Bases: Particles

A class for initializing particles in guiding-center, drift-kinetic or gyro-kinetic models that use the 5D phase space.

The numpy marker array is as follows:

index	0 | 1 | 2 |	3	4	5	6	7	>=8
value	position (eta)	v_parallel	v_perp	weight	s0	w0	buffer

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.

classmethod default_background()

The default background (of type Maxwellian).

property vdim

Dimension of the velocity space.

property n_cols_diagnostics

Number of the diagnostics columns.

property n_cols_aux

Number of the auxiliary columns.

property magn_bckgr

Fluid equilibrium with B.

property absB0_h

Discrete 0-form coefficients of |B_0|.

property unit_b1_h

Discrete 1-form coefficients of B/|B|.

property epsilon

One of equation params, epsilon

property coords

Coordinates of the Particles5D, \((v_\parallel, \mu)\).

property derham

Discrete Deram complex.

svol(eta1, eta2, eta3, *v)

Sampling density function as volume-form.

Parameters:

• eta1 (array_like) – Logical evaluation points.

• eta2 (array_like) – Logical evaluation points.

• eta3 (array_like) – Logical evaluation points.

• *v (array_like) – Velocity evaluation points.

Returns:

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

• ——-

s3(eta1, eta2, eta3, *v)

Sampling density function as 3-form.

Parameters:

• eta1 (array_like) – Logical evaluation points.

• eta2 (array_like) – Logical evaluation points.

• eta3 (array_like) – Logical evaluation points.

• *v (array_like) – Velocity evaluation points.

Returns:

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

• ——-

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

Sampling density function as 0-form.

Parameters:

• eta1 (array_like) – Logical evaluation points.

• eta2 (array_like) – Logical evaluation points.

• eta3 (array_like) – Logical evaluation points.

• v_parallel (array_like) – Velocity evaluation points.

• v_perp (array_like) – 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, verbose: bool = True)

Drawing markers according to the volume density \(s^\textrm{vol}_{\textnormal{in}}\). In Struphy, the initial marker distribution \(s^\textrm{vol}_{\textnormal{in}}\) is always of the form

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

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

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

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

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

The following choices are available in Struphy:

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

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

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

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

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

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

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

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

Then,

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

So then,

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

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

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)

• verbose (bool) – Show info on screen.

save_constants_of_motion()

Calculate each markers’ energy and canonical toroidal momentum and assign them into diagnostics columns of marker array:

diagnostics index	0 |	1 |	2 |
value	energy	magn. moment	can. momentum

Only equilibrium magnetic field is considered.

save_magnetic_energy(PBb)

Calculate magnetic field energy at each particles’ position and assign it into markers[:,self.first_diagnostics_idx].

Parameters:

b2 (BlockVector) – Finite element coefficients of the time-dependent magnetic field.

save_magnetic_background_energy()

Evaluate \(mu_p |B_0(\boldsymbol \eta_p)|\) for each marker. The result is stored at markers[:, self.first_diagnostics_idx,].

save_magnetic_moment()

Calculate magnetic moment of each particles and assign it into markers[:,self.first_diagnostics_idx,+1].

class struphy.pic.particles.Particles3D(**kwargs)

Bases: Particles

A class for initializing particles in 3D configuration space.

The numpy marker array is as follows:

index	0 | 1 | 2 |	3	4	5	>=6
value	position (eta)	weight	s0	w0	buffer

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.

classmethod default_background()

The default background (of type Maxwellian).

property vdim

Dimension of the velocity space.

property n_cols_diagnostics

Number of the diagnostics columns.

property n_cols_aux

Number of the auxiliary columns.

property coords

Coordinates of the Particles3D.

svol(eta1, eta2, eta3)

Sampling density function as volume form.

Parameters:

• eta1 (array_like) – Logical evaluation points.

• eta2 (array_like) – Logical evaluation points.

• eta3 (array_like) – Logical evaluation points.

• *v (array_like) – Velocity evaluation points.

Returns:

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

• ——-

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

Sampling density function as 0 form.

Parameters:

• eta1 (array_like) – Logical evaluation points.

• eta2 (array_like) – Logical evaluation points.

• eta3 (array_like) – Logical evaluation points.

• *v (array_like) – 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.

• ——-

class struphy.pic.particles.ParticlesSPH(**kwargs)

Bases: Particles

A class for initializing particles in SPH models.

The numpy marker array is as follows:

index	0 | 1 | 2 |	3 | 4 | 5 |	6	7	8	>=9
value	position (eta)	velocities	weight	s0	w0	buffer

Parameters:

• name (str) – Name of the particle species.

• **params (dict) – Parameters for markers, see Particles.

classmethod default_background()

The default background (of type Maxwellian).

property vdim

Dimension of the velocity space.

property n_cols_diagnostics

Number of the diagnostics columns.

property n_cols_aux

Number of the auxiliary columns.

property coords

Coordinates of the Particles6D, \((v_1, v_2, v_3)\).

svol(eta1, eta2, eta3, *v)

Sampling density function as volume form.

Parameters:

• eta1 (array_like) – Logical evaluation points.

• eta2 (array_like) – Logical evaluation points.

• eta3 (array_like) – Logical evaluation points.

• *v (array_like) – Velocity evaluation points.

Returns:

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

• ——-

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

Sampling density function as 0 form.

Parameters:

• eta1 (array_like) – Logical evaluation points.

• eta2 (array_like) – Logical evaluation points.

• eta3 (array_like) – Logical evaluation points.

• *v (array_like) – 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.

• ——-

Pusher modules

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

Bases: object

Class for solving particle ODEs

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

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

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

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

Notes

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

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

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

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

Parameters:

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

• kernel (pyccelized function) – The pusher kernel.

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

• args_domain (DomainArguments) – Mapping infos.

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

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

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

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

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

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

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

• verbose (bool) – Whether to print some info or not.

property particles

Particle object.

property kernel

The pyccelized pusher kernel.

property init_kernels

A dict of kernels for initial spline evaluation before iteration.

property eval_kernels

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

property args_kernel

Optional arguments for kernel.

property args_domain

Mandatory Domain arguments.

property n_stages

Number of stages of the pusher.

property maxiter

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

property tol

Iteration terminates when residual<tol.

property mpi_sort

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

property verbose

Print more info.

Particle-to-grid coupling

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

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

Bases: object

Struphy accumulation (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} \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} \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 evaluated at the particle position \(\boldsymbol \eta_p\), and \(A^{\mu,\nu}_p\) and \(B^\mu_p\) are particle-dependent “filling functions”, to be defined in the module accum_kernels.

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 6D accumulation kernels and 5D accumulation kernels for details.

property particles

Particle object.

property kernel: Pyccelkernel

The accumulation kernel.

property derham

Discrete Derham complex on the logical unit cube.

property args_domain

Mapping info for evaluating metric coefficients.

property space_id

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

property form

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

property symmetry

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

property operators

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

property vectors

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

property accfilter

Callable filters

init_control_variate(mass_ops)

Set up the use of noise reduction by control variate.

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

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

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

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

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

Bases: object

Same as 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.

property particles

Particle object.

property kernel: Pyccelkernel

The accumulation kernel.

property derham

Discrete Derham complex on the logical unit cube.

property args_domain

Mapping arguments.

property space_id

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

property form

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

property vectors

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

property accfilter

Callable filters

init_control_variate(mass_ops)

Set up the use of noise reduction by control variate.

show_accumulated_spline_field(mass_ops, eta_direction=0)

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

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

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