adam_prism_fnl_pic_conserving_object

ADAM, PRISM charge-conserving PIC current (londrillo, totorica, esirkepov_Q<N>), FNL backend.

Device port of the CPU charge-conserving current chain of prism_cpu_object (integrate_rk_ssp_pic_charge_conserving and the helpers it calls). Host-side twins of the CPU matrix builders/factorization live here (private) so the CPU backend is left untouched; the arithmetic of every kernel follows the CPU routine named in its header line.

Multi-block and multi-rank (single refinement level): every solver works on a virtual block, the whole domain seen as one NX x NY x NZ block (NX = nbx*ni, ...). The conserving operators are global along grid lines (the inverse of the londrillo line matrix is dense, the matched face flux is a prefix sum from the domain boundary, the P/Q filters reach up to 7/5 cells), so a single-block solve of the whole domain is the reference result: on the virtual block the phase-1 single-block kernels run unchanged and the multi-block result equals the single-block one by construction. The real blocks are touched only at the two interfaces: (i) inputs, Esirkepov directly in global cell indices, the matched rho through the regular block deposit + reduce_ghost_local_gpu + gather of the block interiors; (ii) outputs, J (and the filtered rho) scattered from the virtual block to every real cell, ghosts included, so interface ghosts get the neighbour values and physical ghosts the single-block odd/zero values. With one block the virtual block is the block itself (same origin, spacing and cell centers): phase-1 arithmetic. Several ranks: the virtual block is replicated on every rank. Each rank adds the sources of its own particles on the virtual cells (block interiors and in-domain ghosts, so that no ghost fold is needed), the complete sources are summed over the ranks (reduce_virtual, MPI_Allreduce staged through the host), every rank runs the same virtual solve and scatters the result on its own blocks. The per-particle cleanup sums the particle currents of each rank, then reduces the virtual current.

Device layouts (FNL convention, block index first):

• real q-like fields: [nb, 1-ngc:ni+ngc, 1-ngc:nj+ngc, 1-ngc:nk+ngc, nv];
• virtual fields: the same with nb=1 and NX, NY, NZ (H, scratches, line buffers, the virtual current jv_gpu);
• directional charge changes H: [1, NX, NY, NZ, 3] (interior only, the CPU solvers never read H ghosts);
• particles: [particle, variable] (coalesced across particles).

Like the CPU reference, a particle whose support touches the domain boundary stops the run (zero-flux closure, no particle boundary conditions).

Source: src/app/prism/fnl/adam_prism_fnl_pic_conserving_object.F90

Dependencies

Contents

• prism_fnl_pic_conserving_object
• destroy
• initialize
• prepare
• allocate_gather
• initialize_coherent_gather
• deposit_matched_charge_box
• coherent_field_weighting
• copy_field
• copy_particles
• build_stage_target
• build_centered_trajectory
• esirkepov_displacement
• deposit_matched_current_box
• store_box_current
• combine_stage_source
• solve_dispatch
• solve_lines
• solve_modified
• solve_modified_with_cleanup
• deposit_modified_charge
• impose_current_ghosts
• reduce_virtual4
• reduce_virtual5
• set_particle_box
• reduce_virtual_box
• grow_virtual_pack_buffer
• grow_virtual_host_buffer
• build_virtual_geometry
• scatter_current
• impose_virtual_ghosts
• impose_virtual_ghosts_periodic
• check_source_support
• zero_h
• zero_current
• filter_axis
• box_filter_axis
• transposed_filter_axis
• line_cell
• shape_weights_dev
• check_particle_error
• modified_kernels
• odd_mirror_index
• build_centered_reconstruction_matrix
• factorize_matrix_pivot
• esirkepov_shape_weight_dev
• particle_shape_weight_dev
• effective_gaussian_weight_dev
• gaussian_weight_dev
• support_radius

Variables

Name	Type	Attributes	Description
virtual_host_buffer	real(kind=R8P)	allocatable, target, save
virtual_pack_gpu	real(kind=R8P)	pointer, save
CONSERVING_SOLVER_LONDRILLO	integer(kind=I4P)	parameter	Esirkepov face fluxes + inverse reconstruction.
CONSERVING_SOLVER_MODIFIED	integer(kind=I4P)	parameter	Matched current, totorica or esirkepov_Q (P/Q kernels).
SHAPE_NGP	integer(kind=I4P)	parameter	NGP particle shape.
SHAPE_CIC	integer(kind=I4P)	parameter	CIC particle shape.
SHAPE_TSC	integer(kind=I4P)	parameter	TSC particle shape.
SHAPE_CUBIC	integer(kind=I4P)	parameter	Cubic particle shape.
SHAPE_QUARTIC	integer(kind=I4P)	parameter	Quartic particle shape.
SHAPE_QUINTIC	integer(kind=I4P)	parameter	Quintic particle shape.
SHAPE_SEXTIC	integer(kind=I4P)	parameter	Sextic particle shape.
SHAPE_GAUSSIAN	integer(kind=I4P)	parameter	Gaussian particle shape.
MAX_SPAN	integer(kind=I4P)	parameter	Max 1D cells touched by one particle trajectory (gang-private weights).
LQ_SPAN	integer(kind=I4P)	parameter	Box cells widened by the Q taps (box current, nq <= 10).
LP_SPAN	integer(kind=I4P)	parameter	Box cells widened by the P taps (box current, radius <= 7).
ERR_OUTSIDE	integer(kind=I4P)	parameter	Particle outside the domain.
ERR_CHARGE	integer(kind=I4P)	parameter	Particle charge changed during the trajectory.
ERR_BOUNDARY	integer(kind=I4P)	parameter	Particle shape touches the domain boundary.
ERR_FILTER	integer(kind=I4P)	parameter	esirkepov_Q P filter support reaches the domain boundary.
ERR_SPAN	integer(kind=I4P)	parameter	Trajectory support wider than MAX_SPAN (FNL-only limit).
ERR_GAUSS	integer(kind=I4P)	parameter	Empty Gaussian deposition support.
ERR_PU	integer(kind=I4P)	parameter	B-spline weights do not sum to one on the deposition stencil.
TAU_TAIL	real(kind=R8P)	parameter	Relative tail threshold of the particle current cleanup (CPU twin).
TAU_PU	real(kind=R8P)	parameter	Max |raw B-spline weight sum - 1| before the closure (CPU twin).

Derived Types

prism_fnl_pic_conserving_object

Charge-conserving PIC current, device state and kernels.

Components

Name	Type	Attributes	Description
enabled	logical		Conserving current active.
prepared	logical		Geometry-dependent caches built.
solver	integer(kind=I4P)		Selected solver (CONSERVING_SOLVER_*).
shape	integer(kind=I4P)		Particle shape (SHAPE_*).
filter_deposition	logical		Binomial filter on deposition.
tail_cleanup	logical		Matched-current per-particle tail cleanup.
is_periodic	logical		Periodic particle closure along each axis (host pic%bc).
box_current	logical		Matched current and charge built per particle on its box (periodic faces).
sigma	real(kind=R8P)		Gaussian width.
cutoff_sigma	real(kind=R8P)		Gaussian cutoff in sigma units.
gaussian_support_cells	integer(kind=I4P)		Gaussian support radius in cells.
fdv_order	integer(kind=I4P)		FD order.
hs	integer(kind=I4P)		FD half stencil.
ni	integer(kind=I4P)		Virtual block sizes (whole domain), ghost cells.
nj	integer(kind=I4P)		Virtual block sizes (whole domain), ghost cells.
nk	integer(kind=I4P)		Virtual block sizes (whole domain), ghost cells.
ngc	integer(kind=I4P)		Virtual block sizes (whole domain), ghost cells.
nb	integer(kind=I4P)		Virtual blocks (always 1).
blocks_number	integer(kind=I4P)		Virtual blocks (always 1).
bni	integer(kind=I4P)		Real block sizes.
bnj	integer(kind=I4P)		Real block sizes.
bnk	integer(kind=I4P)		Real block sizes.
bnb	integer(kind=I4P)		Real allocated blocks.
bblocks	integer(kind=I4P)		Real actual blocks (cache key).
vemin	real(kind=R8P)		Virtual block origin (domain minimum).
vdx	real(kind=R8P)		Cell size (single refinement level).
nv	integer(kind=I4P)		Field variables number.
particle_number	integer(kind=I4P)		Particles number owned by this rank.
particle_number_global	integer(kind=I4P)		Particles number, all ranks.
particle_capacity	integer(kind=I4P)		Particle arrays capacity of this rank (issue #112).
nrk	integer(kind=I4P)		RK stages number.
esirkepov_q	integer(kind=I4P)		Q filter coefficients of esirkepov_Q.
solver_name	character(len=:)	allocatable	Input name of the solver (messages).
radius	integer(kind=I4P)		Modified P/Q kernels geometry.
qfirst	integer(kind=I4P)		Modified P/Q kernels geometry.
nq	integer(kind=I4P)		Modified P/Q kernels geometry.
shift	real(kind=R8P)		Base-shape position shift: dx/2 for the odd Q (dual nodes), 0 otherwise.
p	real(kind=R8P)		Modified P kernel (cell centers / transverse lines).
qw	real(kind=R8P)		Modified Q kernel (normal faces to centers).
residual_max	real(kind=R8P)		Max current-solver residual of the current step.
coherent	logical		Coherent matched gather (field_weighting_model totorica or esirkepov_Q).
gradius	integer(kind=I4P)		Coherent gather: taps of its P kernel.
gshift	real(kind=R8P)		Coherent gather: dx/2 for the odd esirkepov_Q (dual nodes), 0 otherwise.
vbox_lo	integer(kind=I4P)		Particle box of the boxed reduce_virtual, lower virtual cells (issue #111).
vbox_hi	integer(kind=I4P)		Particle box of the boxed reduce_virtual, upper virtual cells.
vbox_empty	logical		No rank holds particles: a boxed reduce_virtual has nothing to sum.
off_gpu	integer(kind=I4P)	pointer	Real block origin in the virtual block, cells [bnb,3].
vx_gpu	real(kind=R8P)	pointer	Virtual x cell centers [NX].
vy_gpu	real(kind=R8P)	pointer	Virtual y cell centers [NY].
vz_gpu	real(kind=R8P)	pointer	Virtual z cell centers [NZ].
jv_gpu	real(kind=R8P)	pointer	Virtual current [1,ghosts...,3].
p_gpu	real(kind=R8P)	pointer	P kernel [-7:7].
qw_gpu	real(kind=R8P)	pointer	Q kernel [0:9].
fv1_gpu	real(kind=R8P)	pointer	FV1_CC(1:hs,hs), face reconstruction coefficients.
h_gpu	real(kind=R8P)	pointer	Stage directional charge changes [nb,ni,nj,nk,3,nrk].
q_next_gpu	real(kind=R8P)	pointer	Particle stage target [np,8].
q_ref_gpu	real(kind=R8P)	pointer	Particle reference (step start / t0 minus) [np,8].
q_hist_gpu	real(kind=R8P)	pointer	Particle stage targets history [np,8,nrk].
q_mixed_gpu	real(kind=R8P)	pointer	Shifted positions, odd esirkepov_Q rho [np,8].
work_a_gpu	real(kind=R8P)	pointer	Interior scratch [nb,ni,nj,nk].
work_b_gpu	real(kind=R8P)	pointer	Interior scratch [nb,ni,nj,nk].
line_gpu	real(kind=R8P)	pointer	Line-major scratch for line solves/scans.
gather_gpu	real(kind=R8P)	pointer	Real-block field copy, external-field gather [bnb,...,nv].
gp_gpu	real(kind=R8P)	pointer	Coherent gather P kernel [-7:7].
gq_gpu	real(kind=R8P)	pointer	Coherent gather: E, B filtered by P^T [bnb,...,6].
q_dual_gpu	real(kind=R8P)	pointer	Coherent gather: stage particles at x-dx/2 [np,8].
jp_gpu	real(kind=R8P)	pointer	Single-particle virtual current, cleanup [1,ghosts...,3].
src_gpu	real(kind=R8P)	pointer	Single-particle source, cleanup [nb,ni,nj,nk,3].
src_stage_gpu	real(kind=R8P)	pointer	Single-particle stage sources [nb,ni,nj,nk,3,nrk].
mat_x_gpu	real(kind=R8P)	pointer	Factored x-line matrices [ni,ni,nb].
mat_y_gpu	real(kind=R8P)	pointer	Factored y-line matrices [nj,nj,nb].
mat_z_gpu	real(kind=R8P)	pointer	Factored z-line matrices [nk,nk,nb].
piv_x_gpu	integer(kind=I4P)	pointer	x-line pivots [ni,nb].
piv_y_gpu	integer(kind=I4P)	pointer	y-line pivots [nj,nb].
piv_z_gpu	integer(kind=I4P)	pointer	z-line pivots [nk,nb].

Type-Bound Procedures

Name	Attributes	Description
destroy	pass(self)	Free device data.
initialize	pass(self)	Initialize from host configuration.
prepare	pass(self)	Build geometry-dependent caches.
allocate_gather	pass(self)	Allocate the external-field gather copy.
initialize_coherent_gather	pass(self)	Kernels and buffers of the coherent matched gather.
coherent_field_weighting	pass(self)	Gather E, B with P^T and the base shape.
copy_field	pass(self)	Copy a q-like device field.
copy_particles	pass(self)	Copy a device particle array.
build_stage_target	pass(self)	q_next = q_pic + dt sum_r c_r k_r.
build_centered_trajectory	pass(self)	t=0 virtual trajectory q_pic -/+ dt/2 v.
esirkepov_displacement	pass(self)	Esirkepov directional charge change.
combine_stage_source	pass(self)	H_s = (H_s - sum_r c_r H_r)/c_s.
solve_dispatch	pass(self)	Select the configured current solver.
solve_lines	pass(self)	Londrillo line solves.
solve_modified	pass(self)	Matched current (totorica, esirkepov_Q).
solve_modified_with_cleanup	pass(self)	Per-particle modified current with tail cleanup.
deposit_modified_charge	pass(self)	Unfiltered deposit + P_c filter of rho.
deposit_matched_current_box	pass(self)	Matched current per particle on its box, folded (periodic).
store_box_current	pass(self)	Real q(J) = virtual stage current of the box solver.
impose_current_ghosts	pass(self)	Odd/zero current ghosts on the physical faces of real blocks.
reduce_virtual4	pass(self)	Sum a rank-4 virtual-block array over the ranks.
reduce_virtual5	pass(self)	Sum a rank-5 virtual-block array over the ranks.
reduce_virtual		Sum a virtual-block array over the ranks.
set_particle_box	pass(self)	Global box of the virtual cells the particles reach.
build_virtual_geometry	pass(self)	Virtual block of the whole domain.
scatter_current	pass(self)	Real q(J) = virtual J, ghosts included.
impose_virtual_ghosts	pass(self)	Odd/zero current ghosts of the virtual block.
impose_virtual_ghosts_periodic	pass(self)	Virtual current ghosts with periodic axes.
deposit_matched_charge_box	pass(self)	Matched rho per particle, folded (periodic).
check_source_support	pass(self)	Particle source inside its trajectory support box.
transposed_filter_axis	pass(self)	P^T convolution of the coherent gather along one axis.

Subroutines

destroy

Free device data owned by this object.

subroutine destroy(self)

Arguments

Name	Type	Intent	Attributes	Description
self	class(prism_fnl_pic_conserving_object)	inout		Conserving current object.

Call graph

initialize

Initialize from the host configuration and allocate the step-invariant device buffers. Twin of the CPU validation in prism_cpu_object%initialize plus the allocations of integrate_rk_ssp_pic_charge_conserving. The virtual block (whole domain) is built here: the blocks must already be in place (uniform refinement included).

subroutine initialize(self, pic, grid, field, nb, nv, nrk, fdv_order, hs)

Arguments

Name	Type	Intent	Attributes	Description
self	class(prism_fnl_pic_conserving_object)	inout		Conserving current object.
pic	type(prism_pic_object)	in		Host PIC object.
grid	type(grid_object)	in		Grid.
field	type(field_object)	in		Host field (block geometry).
nb	integer(kind=I4P)	in		Allocated real blocks.
nv	integer(kind=I4P)	in		Field variables number.
nrk	integer(kind=I4P)	in		RK stages number.
fdv_order	integer(kind=I4P)	in		FD order.
hs	integer(kind=I4P)	in		FD half stencil.

Call graph

prepare

Build the geometry-dependent caches: virtual block geometry and the factored line matrices. Twin of CPU ensure_londrillo_current_solver_cache (host factorization, one upload). The matrices span whole virtual lines (the domain extent), one per direction (single refinement level).

subroutine prepare(self, field, grid, force)

Arguments

Name	Type	Intent	Attributes	Description
self	class(prism_fnl_pic_conserving_object)	inout		Conserving current object.
field	type(field_object)	in		Host field.
grid	type(grid_object)	in		Grid.
force	logical	in	optional	Rebuild even if already prepared.

Call graph

allocate_gather

Allocate the field copy used to gather with external fields (CPU twin: q_gather).

subroutine allocate_gather(self)

Arguments

Name	Type	Intent	Attributes	Description
self	class(prism_fnl_pic_conserving_object)	inout		Conserving current object.

Call graph

initialize_coherent_gather

Kernels and buffers of the coherent gather of totorica and esirkepov_Q<N> (pic%set_coherent_gather already ran). CPU twin: the gather_* state of prism_pic_object.

subroutine initialize_coherent_gather(self, pic)

Arguments

Name	Type	Intent	Attributes	Description
self	class(prism_fnl_pic_conserving_object)	inout		Conserving current object.
pic	type(prism_pic_object)	in		Host PIC object.

Call graph

deposit_matched_charge_box

Matched charge rho = P S rho^P built particle by particle and folded into the virtual block along the periodic axes, summed over the ranks, then scattered to every real cell (ghosts beyond a periodic face take the periodic image, physical ghosts 0). CPU twin: deposit_matched_charge_box. Used with periodic faces only.

P is linear and shift invariant, so it commutes with the sum over the particles: the base shape S rho^P is deposited particle by particle, then P is applied once along every axis on the virtual block (box_filter_axis, circular along the periodic axes). Parallel mapping: one gang per particle, the three 1D weights computed sequentially, the deposit one vector lane per cell, atomic.

subroutine deposit_matched_charge_box(self, pic_fnl, q, nv)

Arguments

Name	Type	Intent	Attributes	Description
self	class(prism_fnl_pic_conserving_object)	inout		Conserving current object.
pic_fnl	type(prism_fnl_pic_object)	inout		Device PIC helper.
q	real(kind=R8P)	inout		Real field.
nv	integer(kind=I4P)	in		Charge density slot.

Call graph

coherent_field_weighting

Gather E and B with the transpose of the matched charge operator: P^T on the fields, then the unfiltered base shape (on the dual nodes x_c+dx/2 for the odd esirkepov_Q<N>), the adjoint of deposit_modified_charge. CPU twin: coherent_field_weighting of prism_pic_object.

Multi-block, multi-rank: as the charge deposit, P^T runs on the virtual block. Each field component is copied from the real interiors (every virtual cell has one owner block on one rank, the others add zero), summed over the ranks, filtered, and scattered back to every real cell, ghosts included (interface ghosts get the neighbour values, physical ghosts 0, the zero extension of the single-block filter). Along a periodic axis P^T is circular and the ghosts beyond the face take the periodic image.

subroutine coherent_field_weighting(self, pic_fnl, field_fnl, field, grid, pic_fields_gpu, q, q_pic_gpu, nv)

Arguments

Name	Type	Intent	Attributes	Description
self	class(prism_fnl_pic_conserving_object)	inout		Conserving current object.
pic_fnl	type(prism_fnl_pic_object)	inout		Device PIC helper.
field_fnl	type(field_fnl_object)	in		Device field helper.
field	type(field_object)	in		Host field.
grid	type(grid_object)	in		Grid.
pic_fields_gpu	real(kind=R8P)	inout		Fields at the particles [np,6].
q	real(kind=R8P)	in		Real field.
q_pic_gpu	real(kind=R8P)	in		Gather particles [np,8].
nv	integer(kind=I4P)	in		Field variables number.

Call graph

copy_field

dst = src, full real-block q-like device field (ghosts included).

subroutine copy_field(self, src, dst)

Arguments

Name	Type	Intent	Attributes	Description
self	class(prism_fnl_pic_conserving_object)	in		Conserving current object.
src	real(kind=R8P)	in		Source.
dst	real(kind=R8P)	inout		Destination.

Call graph

copy_particles

dst = src, device particle arrays [np,8].

subroutine copy_particles(self, src, dst)

Arguments

Name	Type	Intent	Attributes	Description
self	class(prism_fnl_pic_conserving_object)	in		Conserving current object.
src	real(kind=R8P)	in		Source.
dst	real(kind=R8P)	inout		Destination.

Call graph

build_stage_target

q_next(1:6) = q_pic(1:6) + dt sum_{r<=s} coeff(r) k_r(1:6), q_next(7:8) = q_pic(7:8). CPU twin: the q_pic_next loop of integrate_rk_ssp_pic_charge_conserving (same r order).

subroutine build_stage_target(self, q_pic, q_pic_rk, s, coeff, dt)

Arguments

Name	Type	Intent	Attributes	Description
self	class(prism_fnl_pic_conserving_object)	inout		Conserving current object.
q_pic	real(kind=R8P)	in		Particles at step start [np,8].
q_pic_rk	real(kind=R8P)	in		Stage right-hand sides [np,8,nrk].
s	integer(kind=I4P)	in		Current stage.
coeff	real(kind=R8P)	in		Stage coefficients (1:s used).
dt	real(kind=R8P)	in		Time step.

Call graph

build_centered_trajectory

q_ref = q_pic - dt/2 v, q_next = q_pic + dt/2 v (positions only). CPU twin: initialize_pic_conserving_current_time_zero.

subroutine build_centered_trajectory(self, q_pic, dt)

Arguments

Name	Type	Intent	Attributes	Description
self	class(prism_fnl_pic_conserving_object)	inout		Conserving current object.
q_pic	real(kind=R8P)	in		Particles [np,8].
dt	real(kind=R8P)	in		Virtual time step.

Call graph

esirkepov_displacement

Symmetric Esirkepov density decomposition. CPU twin: compute_pic_esirkepov_displacement.

Parallel mapping: one gang per particle, one vector lane per cell of the particle's trajectory box. The six 1D weight vectors (3 directions x 2 endpoints) are computed once per gang into gang-private arrays (CPU recomputes them inside the triple loop for the non-Gaussian/non-quartic shapes: same values, fewer evaluations). Different particles may overlap, hence the atomic updates.

Cells are virtual-block (global) indices: a trajectory box may straddle any number of real blocks, only the domain boundary (and MAX_SPAN) limits it.

subroutine esirkepov_displacement(self, field_fnl, q_ref, q_to, h, p_first, p_last)

Arguments

Name	Type	Intent	Attributes	Description
self	class(prism_fnl_pic_conserving_object)	in		Conserving current object.
field_fnl	type(field_fnl_object)	in		Device field helper (unused: virtual geometry).
q_ref	real(kind=R8P)	in		Trajectory start [np,8].
q_to	real(kind=R8P)	in		Trajectory end [np,8].
h	real(kind=R8P)	inout		Directional charge change [1,NX,NY,NZ,3].
p_first	integer(kind=I4P)	in		First particle.
p_last	integer(kind=I4P)	in		Last particle.

Call graph

deposit_matched_current_box

Matched current (totorica, esirkepov_Q<N>) of the trajectories q_ref -> q_to, J_d = Q_d P_e1 P_e2 F_d^E, folded into the virtual block along the periodic axes. CPU twin: deposit_matched_current (the separable factors are explained there). Used with periodic faces only: h receives the stage current in place of the stage source, combined over the stages by combine_stage_source as the sources are (the solve is linear) and written to the real blocks by store_box_current.

Q and P are finite shift-invariant stencils, so they commute with the sum over the particles: each particle deposits only its Esirkepov current F_d^E (face prefix sum of the weight change along d, times the transverse weights), then Q_d P_e1 P_e2 are applied once on the grid (box_filter_axis, circular along the periodic axes). The prefix sum stays per particle, where it is compact (the weight change of one particle sums to zero) and keeps the net current through a periodic face; a prefix sum on the grid would make J nonlocal and lose that k = 0 mode. The support of each particle's current is unchanged, (rr+nq) x (rr+2 radius)^2 cells per direction.

Parallel mapping: one gang per particle. The prefix sums are computed sequentially into a gang-private array, one column per direction and written once (the zero-flux closure past the last weight change is applied when they are read: nvfortran 25.11 -fast can lose a second store into a private array). The deposit is one vector lane per face and transverse cell of the Esirkepov box, atomic (particles overlap). Cells are unwrapped, center emin+(i-1/2)dx, and folded modulo n along the periodic axes; face f (between cells f and f+1) is stored at cell f.

subroutine deposit_matched_current_box(self, q_ref, q_to, h, p_first, p_last, dt)

Arguments

Name	Type	Intent	Attributes	Description
self	class(prism_fnl_pic_conserving_object)	inout		Conserving current object.
q_ref	real(kind=R8P)	in		Trajectory start [np,8].
q_to	real(kind=R8P)	in		Trajectory end [np,8].
h	real(kind=R8P)	inout		Current [1,NX,NY,NZ,3].
p_first	integer(kind=I4P)	in		First particle.
p_last	integer(kind=I4P)	in		Last particle.
dt	real(kind=R8P)	in		Time step.

Call graph

store_box_current

Real q(J) = virtual stage current of the box solver (deposit_matched_current_box + combine_stage_source), ghosts included: the virtual ghosts are the periodic images along the periodic axes and the odd mirror along the others (impose_virtual_ghosts_periodic), then scatter_current. CPU twin: store_matched_current.

subroutine store_box_current(self, q, var_jx, hq)

Arguments

Name	Type	Intent	Attributes	Description
self	class(prism_fnl_pic_conserving_object)	inout		Conserving current object.
q	real(kind=R8P)	inout		Real field.
var_jx	integer(kind=I4P)	in		Jx slot in q.
hq	real(kind=R8P)	in		Virtual current [1,NX,NY,NZ,3].

Call graph

combine_stage_source

H_s = (H_s - sum_{r<s} coeff(r) H_r) / coeff(s). CPU twin: the h_q SSP combination of the conserving integrator.

subroutine combine_stage_source(self, h, s, coeff)

Arguments

Name	Type	Intent	Attributes	Description
self	class(prism_fnl_pic_conserving_object)	in		Conserving current object.
h	real(kind=R8P)	inout		Stage sources [nb,ni,nj,nk,3,nrk].
s	integer(kind=I4P)	in		Current stage.
coeff	real(kind=R8P)	in		Stage coefficients (1:s used).

Call graph

solve_dispatch

Select the configured current reconstruction. CPU twin: solve_pic_charge_conserving_current_dispatch. The solve runs on the virtual block (jv_gpu), then J is scattered to the real blocks, ghosts included.

subroutine solve_dispatch(self, field_fnl, q, var_jx, hq, dt)

Arguments

Name	Type	Intent	Attributes	Description
self	class(prism_fnl_pic_conserving_object)	inout		Conserving current object.
field_fnl	type(field_fnl_object)	in		Device field helper.
q	real(kind=R8P)	inout		Real target field.
var_jx	integer(kind=I4P)	in		Jx slot in q (Jy, Jz follow).
hq	real(kind=R8P)	in		Directional source [1,NX,NY,NZ,3].
dt	real(kind=R8P)	in		Time step.

Call graph

solve_lines

Londrillo line solves along x, y, z: Esirkepov face prefix sum + inverse reconstruction. CPU twins: solve_pic_charge_conserving_current_londrillo, build_esirkepov_face_line, update_esirkepov_face_residual, solve_factored_matrix_pivot.

Parallel mapping: one thread per grid line (b, t1, t2), sequential LU substitution along the line against the cached factored matrix of its block (all lanes of a block read the same matrix entry: broadcast). The line is staged in a line-major scratch line(line + nl*(m-1)), so at each substitution step adjacent lanes touch adjacent words (coalesced), instead of the strided q(b,i,j,k,v) layout (stride nb*(ni+2ngc) for x-lines). The CPU flux array is not stored: the residual pass rebuilds the prefix sum in the same order.

subroutine solve_lines(self, field_fnl, q, var_jx, hq, dt)

Arguments

Name	Type	Intent	Attributes	Description
self	class(prism_fnl_pic_conserving_object)	inout		Conserving current object.
field_fnl	type(field_fnl_object)	in		Device field helper.
q	real(kind=R8P)	inout		Target field.
var_jx	integer(kind=I4P)	in		Jx slot in q.
hq	real(kind=R8P)	in		Directional source [nb,ni,nj,nk,3].
dt	real(kind=R8P)	in		Time step.

Call graph

solve_modified

J_d = Q_d P_transverse F_d^E. CPU twins: solve_pic_charge_conserving_current_modified, filter_modified_axis, modified_current_line, build_esirkepov_face_line.

Parallel mapping: (1) the two transverse P filters are fully parallel separable convolutions (one thread per cell, zero extension outside the block interior, taps summed in the CPU order); (2) the face prefix sum F^E is a scan along each line, one thread per line, written line-major into the scratch; (3) the Q convolution F^E -> J is again one thread per cell, reading the line-major scratch with line-fastest (coalesced) indexing. Only steps (2) is sequential and it is O(n) per line.

subroutine solve_modified(self, field_fnl, q, var_jx, hq, dt)

Arguments

Name	Type	Intent	Attributes	Description
self	class(prism_fnl_pic_conserving_object)	inout		Conserving current object.
field_fnl	type(field_fnl_object)	in		Device field helper.
q	real(kind=R8P)	inout		Target field.
var_jx	integer(kind=I4P)	in		Jx slot in q.
hq	real(kind=R8P)	in		Directional source [nb,ni,nj,nk,3].
dt	real(kind=R8P)	in		Time step.

Call graph

solve_modified_with_cleanup

Per-particle modified current with tail cleanup, then accumulation. CPU twins: solve_modified_current_with_cleanup, cleanup_modified_particle_current, report_particle_tail_error.

The particle loop stays on the host, as on the CPU (each particle needs a whole-grid modified solve, so the work per particle is already a set of grid-parallel kernels). The support box of each particle is obtained with min/max reductions (no device-to-host array copy). The tail check reports the largest offending ratio instead of the first offending cell. Everything runs on the virtual block (support boxes in global cells, so a box may straddle real blocks); the accumulated current is scattered to the real blocks at the end.

subroutine solve_modified_with_cleanup(self, pic_fnl, field_fnl, q, var_jx, q_ref, q_hist, active_stage, alph, beta, dt, mixed_source)

Arguments

Name	Type	Intent	Attributes	Description
self	class(prism_fnl_pic_conserving_object)	inout		Conserving current object.
pic_fnl	type(prism_fnl_pic_object)	in		Device PIC helper.
field_fnl	type(field_fnl_object)	in		Device field helper.
q	real(kind=R8P)	inout		Real target.
var_jx	integer(kind=I4P)	in		Jx slot in q.
q_ref	real(kind=R8P)	in		Trajectory start [np,8].
q_hist	real(kind=R8P)	in		Stage targets [np,8,>=stages].
active_stage	integer(kind=I4P)	in		Active stage (0 = t0 virtual step).
alph	real(kind=R8P)	in		SSP alpha coefficients (host).
beta	real(kind=R8P)	in		SSP beta coefficients (host).
dt	real(kind=R8P)	in		Time step.
mixed_source	real(kind=R8P)	in	optional	Mixed stage source.

Call graph

deposit_modified_charge

Deposit the unfiltered particle shape, then apply P_c in all directions. CPU twin: deposit_pic_modified_charge. The neighbour list must be current (caller runs the grid index first, as on the CPU).

Multi-block: the regular deposit lands in the real blocks (ghost spill folded back with reduce_ghost_local_gpu), the complete interiors are gathered into the virtual block, filtered there (P reaches 7 cells, beyond ngc), and the result is scattered back to every real cell: interior and interface ghosts get the filtered values, physical ghosts 0 (as the single-block CPU, which zeroes the rho ghosts).

subroutine deposit_modified_charge(self, pic_fnl, field_fnl, field, grid, q, nv)

Arguments

Name	Type	Intent	Attributes	Description
self	class(prism_fnl_pic_conserving_object)	inout		Conserving current object.
pic_fnl	type(prism_fnl_pic_object)	inout		Device PIC helper.
field_fnl	type(field_fnl_object)	in		Device field helper.
field	type(field_object)	in		Host field.
grid	type(grid_object)	in		Grid.
q	real(kind=R8P)	inout		Real field.
nv	integer(kind=I4P)	in		Charge density slot.

Call graph

impose_current_ghosts

Odd mirror (zero normal face flux) current ghosts on the PHYSICAL faces of the real blocks; interface ghosts are left to the ghost exchange. CPU twin: impose_odd_current_ghosts (single block: every face is physical). Every ghost write reads a cell that is interior in the component's own direction, so all three components are imposed in one race-free kernel.

subroutine impose_current_ghosts(self, q, var_jx)

Arguments

Name	Type	Intent	Attributes	Description
self	class(prism_fnl_pic_conserving_object)	in		Conserving current object.
q	real(kind=R8P)	inout		Real field.
var_jx	integer(kind=I4P)	in		Jx slot in q.

Call graph

reduce_virtual4

Sum a rank-4 virtual-block device array over the ranks, in place: every rank adds the sources of its own particles, the sum is the source of all the particles. Staged through the host: an MPI_Allreduce on device buffers is not portable (the UCC/UCX collectives of HPC-X read device memory as host memory on WSL, SIGSEGV). No-op on one rank. Collective.

boxed (issue #111): sum only the particle box of the last set_particle_box, packed on the device. The caller guarantees that outside the box either every rank holds zero (sources: the sum is zero there too) or the values are never read (gathered fields). A box larger than half the array takes the whole-array path.

subroutine reduce_virtual4(self, a, boxed)

Arguments

Name	Type	Intent	Attributes	Description
self	class(prism_fnl_pic_conserving_object)	in		Conserving current object.
a	real(kind=R8P)	inout		Device array [1,NX,NY,NZ].
boxed	logical	in	optional	Sum only the particle box.

Call graph

reduce_virtual5

Rank-5 twin of reduce_virtual4 (a boxed sum runs component by component).

subroutine reduce_virtual5(self, a, boxed)

Arguments

Name	Type	Intent	Attributes	Description
self	class(prism_fnl_pic_conserving_object)	in		Conserving current object.
a	real(kind=R8P)	inout		Device array [1,NX,NY,NZ,nc].
boxed	logical	in	optional	Sum only the particle box.

Call graph

set_particle_box

Global box of the virtual cells that the sources and the gathers of the particles of every rank can reach (issue #111), for the boxed reduce_virtual: the cells of qa (and qb, the other end of the trajectories) widened by the base-shape support, the widest P/Q tap and a margin of three cells (the dual-node shift, the face storage), clamped to the domain; along a periodic axis a box that would cross the face takes the whole axis. Empty when no rank holds particles. Collective.

subroutine set_particle_box(self, qa, qb, np)

Arguments

Name	Type	Intent	Attributes	Description
self	class(prism_fnl_pic_conserving_object)	inout		Conserving current object.
qa	real(kind=R8P)	in		Particles [np,8].
qb	real(kind=R8P)	in	optional	Other trajectory ends [np,8].
np	integer(kind=I4P)	in		Particles of this rank.

Call graph

reduce_virtual_box

Sum the box lo:hi of a rank-4 virtual-block device array over the ranks, in place: packed on the device, staged through the host as reduce_virtual4, unpacked. Collective.

subroutine reduce_virtual_box(lo, hi, a)

Arguments

Name	Type	Intent	Attributes	Description
lo	integer(kind=I4P)	in		Box lower virtual cells.
hi	integer(kind=I4P)	in		Box upper virtual cells.
a	real(kind=R8P)	inout		Device array [1,NX,NY,NZ].

Call graph

grow_virtual_pack_buffer

Grow the device packing buffer of reduce_virtual_box to n values at least (it never shrinks).

subroutine grow_virtual_pack_buffer(n)

Arguments

Name	Type	Intent	Attributes	Description
n	integer(kind=I4P)	in		Values needed.

Call graph

grow_virtual_host_buffer

Grow the host staging buffer of reduce_virtual to n values at least (it never shrinks).

subroutine grow_virtual_host_buffer(n)

Arguments

Name	Type	Intent	Attributes	Description
n	integer(kind=I8P)	in		Values needed.

Call graph

build_virtual_geometry

Build the virtual block of the whole domain: sizes, origin, spacing, cell centers and the origin of every local real block. Requirements (error otherwise): one refinement level (all blocks of all ranks with the same spacing), real blocks of all ranks tiling the domain exactly once. With one block the virtual block is the block itself. The virtual cell centers are the real ones of the blocks that hold them (gathered over the ranks): the Esirkepov weights use the field geometry of the run, identical on every rank.

subroutine build_virtual_geometry(self, field, grid)

Arguments

Name	Type	Intent	Attributes	Description
self	class(prism_fnl_pic_conserving_object)	inout		Conserving current object.
field	type(field_object)	in		Host field.
grid	type(grid_object)	in		Grid.

Call graph

scatter_current

Real q(J) = virtual J at the same global cell, for every real cell including ghosts: interface ghosts receive the neighbour values, physical ghosts the virtual (single-block) odd/zero ghosts. Other variables are untouched.

subroutine scatter_current(self, q, var_jx)

Arguments

Name	Type	Intent	Attributes	Description
self	class(prism_fnl_pic_conserving_object)	in		Conserving current object.
q	real(kind=R8P)	inout		Real field.
var_jx	integer(kind=I4P)	in		Jx slot in q.

Call graph

impose_virtual_ghosts

Odd mirror current ghosts of the virtual block (all its faces are physical). CPU twin: impose_odd_current_ghosts.

subroutine impose_virtual_ghosts(self, q, var_jx)

Arguments

Name	Type	Intent	Attributes	Description
self	class(prism_fnl_pic_conserving_object)	in		Conserving current object.
q	real(kind=R8P)	inout		Virtual field.
var_jx	integer(kind=I4P)	in		Jx slot in q.

Call graph

impose_virtual_ghosts_periodic

Current ghosts of the virtual block with periodic axes: along a periodic axis every component takes the periodic image, along the others the normal component is the odd mirror (impose_virtual_ghosts). One kernel per axis, in order, so that the edge and corner ghosts take the images of already filled ghosts. CPU twin: impose_odd_current_ghosts.

subroutine impose_virtual_ghosts_periodic(self, q, var_jx)

Arguments

Name	Type	Intent	Attributes	Description
self	class(prism_fnl_pic_conserving_object)	in		Conserving current object.
q	real(kind=R8P)	inout		Virtual field.
var_jx	integer(kind=I4P)	in		Jx slot in q.

Call graph

check_source_support

Stop if the directional source of one particle has nonzero cells outside the support box built from its trajectory end points. Box and source must describe the same trajectory; a mismatch would otherwise surface downstream as a misleading "current outside theoretical support" tail error. CPU twin: the same check in solve_modified_current_with_cleanup.

subroutine check_source_support(self, h, slo, shi, particle)

Arguments

Name	Type	Intent	Attributes	Description
self	class(prism_fnl_pic_conserving_object)	in		Conserving current object.
h	real(kind=R8P)	in		Particle source [1,NX,NY,NZ,3].
slo	integer(kind=I4P)	in		Support box, lower cell.
shi	integer(kind=I4P)	in		Support box, upper cell.
particle	integer(kind=I4P)	in		Particle index (report only).

Call graph

zero_h

h = 0 over the interior directional source.

subroutine zero_h(self, h)

Arguments

Name	Type	Intent	Attributes	Description
self	class(prism_fnl_pic_conserving_object)	in		Conserving current object.
h	real(kind=R8P)	inout		Directional source [nb,ni,nj,nk,3].

Call graph

zero_current

q(J) = 0 over the full extent (CPU: q(var_Jx:var_Jz,:,:,:😅 = 0).

subroutine zero_current(self, q, var_jx)

Arguments

Name	Type	Intent	Attributes	Description
self	class(prism_fnl_pic_conserving_object)	in		Conserving current object.
q	real(kind=R8P)	inout		Field.
var_jx	integer(kind=I4P)	in		Jx slot in q.

Call graph

filter_axis

Separable P convolution along one axis, zero extension outside the block interior. CPU twin: filter_modified_axis (taps summed from -radius to +radius, as on the CPU).

subroutine filter_axis(self, input, output, axis)

Arguments

Name	Type	Intent	Attributes	Description
self	class(prism_fnl_pic_conserving_object)	in		Conserving current object.
input	real(kind=R8P)	in		Input [nb,ni,nj,nk].
output	real(kind=R8P)	inout		Output [nb,ni,nj,nk].
axis	integer(kind=I4P)	in		Filter axis.

Call graph

box_filter_axis

Convolution along one axis of the box current, output_c = sum_{t=0}^{ntaps-1} w(woff+t) input_{c+first+t}: circular along a periodic axis (the fold of deposit_matched_current_box), zero extension otherwise (no box reaches the non-periodic faces, ERR_FILTER). Q: w = qw, woff = 1, first = qfirst, ntaps = nq (faces to cells, face f stored at cell f); P: w = p, woff = 8-radius, first = -radius, ntaps = 2 radius+1.

subroutine box_filter_axis(self, input, output, axis, w, woff, first, ntaps)

Arguments

Name	Type	Intent	Attributes	Description
self	class(prism_fnl_pic_conserving_object)	in		Conserving current object.
input	real(kind=R8P)	in		Input [nb,ni,nj,nk].
output	real(kind=R8P)	inout		Output [nb,ni,nj,nk].
axis	integer(kind=I4P)	in		Filter axis.
w	real(kind=R8P)	in	pointer	Taps (device).
woff	integer(kind=I4P)	in		Index in w of the first tap.
first	integer(kind=I4P)	in		Offset of the first tap.
ntaps	integer(kind=I4P)	in		Taps number.

Call graph

transposed_filter_axis

Transpose of filter_axis with the coherent gather kernel: output_m = sum_s p(s) input_{m-s}, zero extension (circular along a periodic axis). CPU twin: transposed_filter_axis of prism_pic_object.

subroutine transposed_filter_axis(self, input, output, axis)

Arguments

Name	Type	Intent	Attributes	Description
self	class(prism_fnl_pic_conserving_object)	in		Conserving current object.
input	real(kind=R8P)	in		Input [nb,ni,nj,nk].
output	real(kind=R8P)	inout		Output [nb,ni,nj,nk].
axis	integer(kind=I4P)	in		Filter axis.

Call graph

line_cell

Cell (i,j,k) of the m-th element of the line (t1,t2) along dir.

Attributes: pure

subroutine line_cell(dir, m, t1, t2, i, j, k)

Arguments

Name	Type	Intent	Attributes	Description
dir	integer(kind=I4P)	in
m	integer(kind=I4P)	in
t1	integer(kind=I4P)	in
t2	integer(kind=I4P)	in
i	integer(kind=I4P)	out
j	integer(kind=I4P)	out
k	integer(kind=I4P)	out

Call graph

shape_weights_dev

1D Esirkepov weights of one endpoint on the trajectory box. CPU twins: compute_esirkepov_gaussian_shape_1d, closed_shape_1d (unfiltered B-splines, closed weights of bspline_closed_weights), esirkepov_shape_weight_1d.

Attributes: pure

subroutine shape_weights_dev(shp, filter, sigma, cutoff_limit, radius, x, ic, span, dx, xc, w, perr)

Arguments

Name	Type	Intent	Attributes	Description
shp	integer(kind=I4P)	in		Shape.
filter	logical	in		Binomial filter on the Esirkepov shape.
sigma	real(kind=R8P)	in		Gaussian width.
cutoff_limit	real(kind=R8P)	in		Gaussian cutoff.
radius	integer(kind=I4P)	in		Support radius.
x	real(kind=R8P)	in		Particle coordinate.
ic	integer(kind=I4P)	in		Particle cell, local box index.
span	integer(kind=I4P)	in		Box cells.
dx	real(kind=R8P)	in		Cell size.
xc	real(kind=R8P)	in		Cell centers of the box.
w	real(kind=R8P)	out		Weights on the box.
perr	integer(kind=I4P)	inout		Error code.

Call graph

check_particle_error

Translate a device particle error code into the CPU error message.

subroutine check_particle_error(ierr)

Arguments

Name	Type	Intent	Attributes	Description
ierr	integer(kind=I4P)	in		Error code.

Call graph

modified_kernels

P = A Q (base-shape nodes to centers) and Q (face prefix sums to centers) of esirkepov_Q<N>. CPU twin: modified_kernels (offsets and the half-cell shift of the odd filters documented there).

subroutine modified_kernels(esirkepov_q, p, qw, radius, qfirst, nq)

Arguments

Name	Type	Intent	Attributes	Description
esirkepov_q	integer(kind=I4P)	in
p	real(kind=R8P)	out
qw	real(kind=R8P)	out
radius	integer(kind=I4P)	out
qfirst	integer(kind=I4P)	out
nq	integer(kind=I4P)	out

Call graph

odd_mirror_index

Odd mirror of a line index. CPU twin: odd_mirror_index.

subroutine odd_mirror_index(idx, n, col, sgn)

Arguments

Name	Type	Intent	Attributes	Description
idx	integer(kind=I4P)	in
n	integer(kind=I4P)	in
col	integer(kind=I4P)	out
sgn	real(kind=R8P)	out

Call graph

build_centered_reconstruction_matrix

Face-reconstruction matrix R such that FD1_CC J = (R_{i+1/2}-R_{i-1/2})/dx. CPU twin of the same name.

subroutine build_centered_reconstruction_matrix(hs, n, a)

Arguments

Name	Type	Intent	Attributes	Description
hs	integer(kind=I4P)	in
n	integer(kind=I4P)	in
a	real(kind=R8P)	inout

Call graph

factorize_matrix_pivot

LU with partial pivoting, in place. CPU twin of the same name.

subroutine factorize_matrix_pivot(a, pivot, n)

Arguments

Name	Type	Intent	Attributes	Description
a	real(kind=R8P)	inout
pivot	integer(kind=I4P)	out
n	integer(kind=I4P)	in

Call graph

Functions

esirkepov_shape_weight_dev

Optionally binomial-filtered 1D particle shape. CPU twin: esirkepov_shape_weight_1d.

Attributes: pure

Returns: real(kind=R8P)

function esirkepov_shape_weight_dev(shp, filter, x, xc, dx) result(w)

Arguments

Name	Type	Intent	Attributes	Description
shp	integer(kind=I4P)	in
filter	logical	in
x	real(kind=R8P)	in
xc	real(kind=R8P)	in
dx	real(kind=R8P)	in

Call graph

particle_shape_weight_dev

1D particle shape (polynomial forms of the CPU). CPU twin: particle_shape_weight_1d (Gaussian not reached here).

Attributes: pure

Returns: real(kind=R8P)

function particle_shape_weight_dev(shp, x, xc, dx) result(w)

Arguments

Name	Type	Intent	Attributes	Description
shp	integer(kind=I4P)	in
x	real(kind=R8P)	in
xc	real(kind=R8P)	in
dx	real(kind=R8P)	in

Call graph

effective_gaussian_weight_dev

Optionally binomial-filtered 1D Gaussian weight. Same expressions as prism_pic_object%effective_gaussian_weight.

Attributes: pure

Returns: real(kind=R8P)

function effective_gaussian_weight_dev(r, shift, cutoff_limit, filter) result(weight)

Arguments

Name	Type	Intent	Attributes	Description
r	real(kind=R8P)	in
shift	real(kind=R8P)	in
cutoff_limit	real(kind=R8P)	in
filter	logical	in

Call graph

gaussian_weight_dev

Unnormalized 1D Gaussian with compact cutoff.

Attributes: pure

Returns: real(kind=R8P)

function gaussian_weight_dev(r, cutoff_limit) result(weight)

Arguments

Name	Type	Intent	Attributes	Description
r	real(kind=R8P)	in
cutoff_limit	real(kind=R8P)	in

Call graph

support_radius

Integer support radius of the particle shape. CPU twin: particle_weighting_radius.

Returns: integer(kind=I4P)

function support_radius(self) result(rx)

Arguments

Name	Type	Intent	Attributes	Description
self	class(prism_fnl_pic_conserving_object)	in		Conserving current object.

Call graph