Extensions

This document describes the built-in extensions available in π-PIC, their physical models, parameters, and usage examples.

Overview

Extensions add physics modules to simulations via custom handlers. They are registered using sim.add_handler() and operate on particles or fields during the simulation advance.

Extensions

QED Extensions

qed_gonoskov2015

Description: QED photon emission and pair production using the locally constant field approximation (LCFA) with rejection sampling.

Reference: A. Gonoskov et al., PRE (2015), arXiv:1412.6426

Physics:

• Photon emission (synchrotron radiation) from electrons/positrons

• Pair production from photons

• Based on synchrotron functions with efficient approximations

Parameters:

handler = qed.handler(electron_type, positron_type, photon_type,
                     probability_threshold=1e-3,
                     probability_subcycle=0.1)

• electron_type (int): Type index for electrons (use sim.get_type_index('electron')).

• positron_type (int): Type index for positrons (use sim.get_type_index('positron')).

• photon_type (int): Type index for photons (use sim.get_type_index('photon')).

• probability_threshold (float): Below this threshold, event weights are increased and events are rarified. Default: 1e-3.

• probability_subcycle (float): Maximum estimated probability per QED event within a single substep. Default: 0.1.

Usage:

import pipic.extensions.qed_gonoskov2015 as qed

# Must have electron, positron, and photon types defined
sim.add_particles(name='electron', ...)
sim.add_particles(name='positron', ...)
sim.add_particles(name='photon', charge=0, mass=0, ...)

# Get type indices
electron_type = sim.get_type_index('electron')
positron_type = sim.get_type_index('positron')
photon_type = sim.get_type_index('photon')

handler_ptr = qed.handler(electron_type, positron_type, photon_type)
sim.add_handler(name='qed_gonoskov2015',
               subject='electron, positron, photon',
               handler=handler_ptr)

qed_volokitin2023

Description: Optimized QED event generator using minimal rate computations per event (faster equivalent to qed_gonoskov2015).

Reference: V. Volokitin et al., JCS 74, 102170 (2023) - DOI: 10.1016/j.jocs.2023.102170

Physics:

• Photon emission (Compton scattering) from electrons/positrons

• Pair production (Breit-Wheeler process) from photons

• Based on Fast_QED module from pyHiChi

• Optimized algorithm minimizing rate computations

Developer: Joel Magnusson (joel.magnusson@physics.gu.se), based on pyHiChi implementation

Parameters:

handler = qed.handler(electron_type, positron_type, photon_type)

• electron_type (int): Type index for electrons (use sim.get_type_index('electron')).

• positron_type (int): Type index for positrons (use sim.get_type_index('positron')).

• photon_type (int): Type index for photons (use sim.get_type_index('photon')).

Usage:

import pipic.extensions.qed_volokitin2023 as qed

# Must have electron, positron, and photon types defined
sim.add_particles(name='electron', ...)
sim.add_particles(name='positron', ...)
sim.add_particles(name='photon', charge=0, mass=0, ...)

# Get type indices
electron_type = sim.get_type_index('electron')
positron_type = sim.get_type_index('positron')
photon_type = sim.get_type_index('photon')

handler_ptr = qed.handler(electron_type, positron_type, photon_type)
sim.add_handler(name='qed_volokitin2023',
               subject='electron, positron, photon',
               handler=handler_ptr)

Notes:

• Equivalent physics to qed_gonoskov2015 but with optimized performance.

• Recommended over qed_gonoskov2015 for production simulations.

• See examples/qed_volokitin2023_test.py for usage example.

Particle Management

downsampler_gonoskov2022

Description: Agnostic conservative downsampling for dynamically reducing particle count while preserving distributions and conserved quantities.

Reference: A. Gonoskov, CPC 271, 108200 (2022) - DOI: 10.1016/j.cpc.2021.108200, arXiv:1607.03755

Physics:

• Reduces particle ensemble without introducing flattening or systematic distribution changes

• Preserves conserved quantities: total weight, energy, momentum, CIC contributions

• Operates on subsets of particles contributing to nearby grid nodes (2/4/8 nodes depending on dimensionality)

• Can be configured for multiple particle types via add_assignment()

Developer: Arkady Gonoskov (arkady.gonoskov@physics.gu.se)

Parameters:

handler = downsampler.handler(ensemble_data, type_index,
                             preserve_energy=True,
                             preserve_momentum=True,
                             preserve_cic_weight=True,
                             cap=15,
                             target_ratio=1.0)

• ensemble_data (int): Ensemble data pointer from sim.ensemble_data().

• type_index (int): Type index of particles to downsample.

• preserve_energy (bool): Preserve energy of particle subsets. Default: True.

• preserve_momentum (bool): Preserve momentum of particle subsets. Default: True.

• preserve_cic_weight (bool): Preserve CIC contributions to grid nodes. Default: True.

• cap (int): Threshold for number of particles per subset (effectively per cell). Default: 15.

• target_ratio (float): Target particle count after downsampling as fraction of cap (≤1.0). Default: 1.0.

Additional Configuration:

downsampler.add_assignment(type_index, preserve_energy=True,
                          preserve_momentum=True, preserve_cic_weight=True,
                          cap=15, target_ratio=1.0)

• Use to configure downsampling for additional particle types after initialization.

• Parameters same as handler() except ensemble_data not needed.

Usage:

import pipic.extensions.downsampler_gonoskov2022 as downsampler

# Get ensemble data pointer
ensemble_ptr = sim.ensemble_data()

# Get type index for electrons
electron_type = sim.get_type_index('electron')

# Initialize handler
handler_ptr = downsampler.handler(ensemble_ptr, electron_type,
                                 preserve_energy=True,
                                 preserve_momentum=True,
                                 preserve_cic_weight=True,
                                 cap=15,
                                 target_ratio=1.0)

# Add additional particle types if needed
positron_type = sim.get_type_index('positron')
downsampler.add_assignment(positron_type, cap=15, target_ratio=1.0)

# Register handler (must use 'cells' as subject)
sim.add_handler(name='downsampler',
               subject='cells',
               handler=handler_ptr)

Notes:

• Subject must be set to ‘cells’, not specific particle types.

• Downsampling applied when particle count in a cell exceeds cap.

• Reduces particles to target_ratio * cap when triggered.

• See examples/downsampler_gonoskov2022_test.py for usage example.

Radiation Reaction

landau_lifshitz

Description: Account for radiation reaction using leading terms of the Landau-Lifshitz model.

Physics:

• Classical radiation reaction force on charged particles

• Based on Landau-Lifshitz equation (leading order)

• Accounts for energy loss due to radiation in strong electromagnetic fields

Developer: Joel Magnusson (joel.magnusson@physics.gu.se)

Parameters:

handler = landau_lifshitz.handler()

• No parameters required for initialization.

Usage:

import pipic.extensions.landau_lifshitz as ll

# Initialize handler
handler_ptr = ll.handler()

# Register handler for electrons and positrons
sim.add_handler(name='landau_lifshitz',
               subject='electron, positron',
               handler=handler_ptr)

Notes:

• Applies classical radiation reaction force to particles.

• Most relevant for high-intensity laser-plasma interactions where radiation damping is significant.

• Can be applied to any charged particle types.

Field Generation

focused_pulse

Description: Inject a focused Gaussian laser pulse via a field handler (paraxial approximation).

Configuration functions:

focused_pulse.set_box(nx, ny, nz, xmin, ymin, zmin, xmax, ymax, zmax)
focused_pulse.set_center(x, y, z)
focused_pulse.set_path(x, y, z)          # Propagation direction
focused_pulse.set_e_axis(x, y, z)        # Polarization direction
focused_pulse.set_theta_max(theta_max)   # Angular aperture
focused_pulse.set_l_size(l_size)         # Characteristic length
focused_pulse.set_shape(shape)           # Temporal shape function
focused_pulse.add_pulse()                # Add configured pulse
focused_pulse.clear_pulse()              # Clear all pulses

Usage:

import pipic.extensions.focused_pulse as fp

fp.set_box(sim.nx, sim.ny, sim.nz,
          sim.xmin, sim.ymin, sim.zmin,
          sim.xmax, sim.ymax, sim.zmax)
fp.set_center(0, 0, 0)
fp.set_path(1, 0, 0)
fp.set_e_axis(0, 1, 0)
fp.set_theta_max(30.0)
fp.set_l_size(3e-4)
fp.set_shape(lambda t: ...)
fp.add_pulse()

field_handler_ptr = fp.field_loop_cb()
sim.add_handler(name='focused_pulse',
               subject='cells',
               field_handler=field_handler_ptr)

Notes:

• Only use as a field handler (no particle handler).

• Paraxial approximation; keep beam waist and angle consistent with grid.

Boundary Handling

x_reflector_c

Description: Reflect particles in a finite region along x.

Parameters:

handler = x_reflector.handler(location, thickness)

• location (float): Center position of reflective region (cm).

• thickness (float): Thickness of reflective region (cm).

Usage:

import pipic.extensions.x_reflector_c as reflector

handler_ptr = reflector.handler(location=1e-3, thickness=1e-5)
sim.add_handler(name='x_reflector',
               subject='electron, positron',
               handler=handler_ptr)

Behavior:

• Particles within the region with \(p_x > 0\) have momentum reversed: \(p_x \to -p_x\).

• Region spans location - thickness/2 to location + thickness/2.

x_converter_c

Description: Convert particles to another type when traversing an x-plane region.

Parameters:

handler = x_converter.handler(location, thickness, typeTo)

• location (float): Center position of conversion region (cm).

• thickness (float): Thickness of conversion region (cm).

• typeTo (int): Target particle type index.

Usage:

import pipic.extensions.x_converter_c as converter

photon_idx = sim.get_type_index('photon')
handler_ptr = converter.handler(location=0, thickness=1e-4, typeTo=photon_idx)
sim.add_handler(name='x_converter',
               subject='electron',
               handler=handler_ptr)

Notes:

• Converted particles keep position, momentum, and weight.

absorbing_boundaries

Description: Absorb particles and damp fields near boundaries.

Parameters (particle handler):

handler = absorber.handler(ensemble_data, simulation_box, characteristic_wavelength,
                          density_profile=-1, boundary_size=-1.0, axis='x',
                          fall=-1.0, temperature=0.0, particles_per_cell=1.0,
                          remove_particles_every=10, moving_window_velocity=0.0,
                          moving_window_direction='x')

• ensemble_data (int): Ensemble data pointer from sim.ensemble_data().

• simulation_box (int): Pointer from sim.simulation_box().

• characteristic_wavelength (float): Characteristic wavelength (cm).

• density_profile (int): Density profile callback address. Default: -1.

• boundary_size (float): Layer thickness (cm). Default: -1.0 (auto).

• axis (char): Boundary axis. Default: ‘x’.

• fall (float): Field damping coefficient. Default: -1.0.

• temperature (float): Temperature for re-injected particles. Default: 0.0.

• particles_per_cell (float): Density for re-injection. Default: 1.0.

• remove_particles_every (int): Removal frequency. Default: 10.

• moving_window_velocity (float): Moving window speed (cm/s). Default: 0.0.

• moving_window_direction (char): Moving window axis. Default: ‘x’.

Parameters (field handler):

field_handler = absorber.field_handler(simulation_box, characteristic_wavelength,
                                      boundary_size=-1.0, axis='x', fall=-1.0)

• simulation_box (int): Simulation box pointer.

• characteristic_wavelength (float): Characteristic wavelength (cm).

• boundary_size (float): Absorption layer thickness (cm). Default: -1.0.

• axis (char): Boundary axis. Default: ‘x’.

• fall (float): Field damping coefficient. Default: -1.0.

Usage:

import pipic.extensions.absorbing_boundaries as absorber

ensemble_ptr = sim.ensemble_data()
simbox_ptr = sim.simulation_box()

handler_ptr = absorber.handler(ensemble_ptr, simbox_ptr,
                               characteristic_wavelength=0.8e-4,
                               boundary_size=20*dx,
                               axis='x')
field_handler_ptr = absorber.field_handler(simbox_ptr,
                                           characteristic_wavelength=0.8e-4,
                                           boundary_size=20*dx,
                                           axis='x')

sim.add_handler(name='absorber',
               subject='all_types',
               handler=handler_ptr,
               field_handler=field_handler_ptr)

Behavior:

• Particle weights reduced in absorption layer; optional re-injection.

• Fields damped near boundaries to suppress reflections.

Moving Window

moving_window

Description: Move the simulation window while maintaining plasma injection.

Parameters (particle handler):

handler = mw.handler(simulation_box, particles_per_cell,
                    temperature, density_profile,
                    thickness=-1, velocity=lightVelocity, axis='x', angle=0)

• simulation_box (int): Pointer from sim.simulation_box().

• particles_per_cell (float): Particle density for re-injection.

• temperature (float): Temperature for re-injected particles (erg).

• density_profile (int): Density profile callback address.

• thickness (float): Boundary thickness (cm). Default: -1.

• velocity (float): Window velocity (cm/s). Default: lightVelocity.

• axis (char): Movement axis. Default: ‘x’.

• angle (float): Window angle (radians). Default: 0.

Parameters (field handler):

field_handler = mw.field_handler(simulation_box, thickness=-1,
                                 velocity=lightVelocity, axis='x', angle=0)

• simulation_box (int): Simulation box pointer.

• thickness (float): Boundary thickness (cm). Default: -1.

• velocity (float): Window velocity (cm/s). Default: lightVelocity.

• axis (char): Movement axis. Default: ‘x’.

• angle (float): Window angle (radians). Default: 0.

Usage:

import pipic.extensions.moving_window as mw
from pipic import consts, types

@cfunc(types.add_particles_callback)
def density_profile(r, data_double, data_int):
    return 1e18

simbox_ptr = sim.simulation_box()

handler_ptr = mw.handler(simbox_ptr,
                        particles_per_cell=10,
                        temperature=1e-6*consts.electron_mass*consts.light_velocity**2,
                        density_profile=density_profile.address,
                        velocity=consts.light_velocity,
                        axis='x')
field_handler_ptr = mw.field_handler(simbox_ptr,
                                    velocity=consts.light_velocity,
                                    axis='x')

sim.add_handler(name='moving_window',
               subject='all_types',
               handler=handler_ptr,
               field_handler=field_handler_ptr)

Behavior:

• Shifts simulation domain along axis at velocity.

• Removes particles exiting the rear; optionally re-injects at front via density_profile.

Extension Development

π-PIC offers the possibility to develop custom extensions in Python, C/C++, Fortran, and other languages that can produce a callable function for Python. Extensions can modify, add, and remove particles based on local field state, and can also modify field state. This enables implementing physics like ionization, radiation reaction, QED processes, etc.

Overview

Extensions connect to π-PIC via two types of an interfaces: one for particles and one for fields. The particle interface is based on the cellInterface structure, which provides direct data access without performance loss. The field access is realized through the fieldLoop() method defined for every field_solver. The interfaces are defined in:

• Python: /pipic/interfaces/cellinterface.py

• C/C++: src/interfaces.h

Extensions must provide Handler and/or fieldHandler functions that process particle subsets and/or electromagnetic fields. The Handler is called during advance() for each subset of particles in each cell for each specified type, while the fieldHandler is called for each gridpoint.

Key capabilities:

• Add new particles within the same cell (buffered until next loop), enabling implementation of particle creation mechanisms

• Process particles by type or process all cells via subject='cells' for applications like ionization

• Pass custom data via data_double and data_int parameters for storing extension-specific state across calls

Example registration:

import my_extension
from pipic import addressof

# Initialize handlers
handler_ptr = my_extension.handler(param1, param2)
field_handler_ptr = my_extension.field_handler(param3)

# Optional: allocate custom data arrays
data_double = np.zeros((10,), dtype=np.double)
data_int = np.zeros((5,), dtype=np.intc)

# Register both particle and field handlers
sim.add_handler(name=my_extension.name,
               subject='electron, positron',     # Apply particle handler to these types
               handler=handler_ptr,               # Processes particles
               field_handler=field_handler_ptr,   # Processes fields (optional)
               data_double=addressof(data_double),  # Custom data (optional)
               data_int=addressof(data_int))        # Custom data (optional)

Python Extensions

Python extensions use Numba C callbacks for high performance while maintaining Python flexibility.

Extension Structure:

1. Set extension name (must match filename):

   name = "my_extension"
   

2. Allocate shared data (if needed, must be thread-safe):

   from numba import int32, float64
   from ctypes import c_double, c_int
   
   DataDouble = (c_double * 10)()  # Array of 10 doubles
   DataInt = (c_int * 5)()         # Array of 5 integers
   

3. Implement Handler using CellInterface:

   from numba import cfunc
   from pipic.types import handler_callback
   from pipic.interfaces import CellInterface
   
   @cfunc(handler_callback)
   def Handler(CI_I, CI_D, CI_F, CI_P, CI_NP, data_double, data_int):
       C = CellInterface(CI_I, CI_D, CI_F, CI_P, CI_NP)
   
       # Process each particle in the subset
       for ip in range(C.particleSubsetSize):
           P = C.Particle(ip)
   
           # Access electromagnetic field at particle position
           E, B = C.interpolateField(P.r)
   
           # Modify particle momentum/position
           P.p.x += ...
   
           # Add new particle if needed
           if C.particleBufferSize < C.particleBufferCapacity:
               newP = C.newParticle(C.particleBufferSize)
               newP.r = P.r
               newP.p = ...
               newP.w = P.w
               newP.id = target_type_index
               C.particleBufferSize += 1
   

4. Create initialization function returning Handler address:

   def handler(param1, param2, ...):
       # Configure extension with parameters
       DataDouble[0] = param1
       DataInt[0] = param2
       return Handler.address
   

5. Provide data accessors (if using shared data):

   from ctypes import addressof
   
   def data_double():
       return addressof(DataDouble)
   
   def data_int():
       return addressof(DataInt)
   

Usage example:

import my_extension

handler_ptr = my_extension.handler(param1=1.0, param2=42)
sim.add_handler(name=my_extension.name,
               subject='electron',
               handler=handler_ptr,
               data_double=my_extension.data_double(),
               data_int=my_extension.data_int())

Reference: See x_reflector_py.py and its usage in x_reflector_py_test.py.

C/C++ Extensions

C/C++ extensions provide maximum flexibility and performance. Unlike Python extensions using Numba, C/C++ extensions are compiled independently and can achieve full CPU performance. The following outlines both local development and the recommended contribution workflow.

Extension Structure:

1. Include required headers:

   #include "interfaces.h"
   #include <pybind11/pybind11.h>
   #include "pybind11/stl.h"
   #include <pybind11/operators.h>
   

2. Set extension name and global static parameters (used for module export):

   const string name = "my_extension";
   static double globalParameter;
   

3. Implement Handler function processing particles:

   // Particle handler: processes particles of specified types
   void Handler(int *I, double *D, double *F, double *P,
               double *NP, double *dataDouble, int *dataInt) {
       cellInterface CI(I, D, F, P, NP);
   
       // Generate new particles when CI.particleTypeIndex == -1 (cell-level handler)
       if (CI.particleTypeIndex == -1) {
           // Called once per cell, independent of particle types
           // Use this for particle generation/injection
   
           if(CI.particleBufferSize < CI.particleBufferCapacity) {
               particle *newP = CI.newParticle(CI.particleBufferSize);
               newP->r = ...;  // Set position
               newP->p = ...;  // Set momentum
               newP->w = ...;  // Set weight
               newP->id = target_type_index;
               CI.particleBufferSize++;
           }
       }
       // Modify/remove particles when CI.particleTypeIndex == particleTypeIndex
       else if (CI.particleTypeIndex == particleTypeIndex) {
           // Called for each particle subset of the specified type
           // Use this for particle modification or removal
   
           for(int ip = 0; ip < CI.particleSubsetSize; ip++) {
               particle *P = CI.Particle(ip);
   
               // Interpolate EM field at particle position
               double3 E, B;
               CI.interpolateField(P->r, E, B);
   
               // Modify particle momentum. Note: do not modify positions since it can cause buffer overflow.
               P->p.x += ...;
   
               // Remove particle by setting weight to zero
               // P->w = 0.0;
           }
       }
   }
   

4. Implement initialization function returning Handler address:

   int64_t handler(double param1, int param2) {
       globalParameter = param1;
       return (int64_t)Handler;
   }
   

5. Implement fieldHandler function processing electromagnetic fields:

   // Field handler: processes electromagnetic fields
   void fieldHandler(int *ind, double *r, double *E, double *B,
                    double *dataDouble, int *dataInt) {
       // Modify field at position r and grid index ind
       E[0] *= damping_factor;
       // ...
   }
   

6. Implement initialization function returning fieldHandler address:

   int64_t field_handler(double param1) {
       globalParameter = param1;
       return (int64_t)fieldHandler;
   }
   

7. Export via pybind11:

   namespace py = pybind11;
   PYBIND11_MODULE(_my_extension, object) {
       object.attr("name") = name;
       object.def("handler", &handler,
                 py::arg("param1"), py::arg("param2"));
       object.def("field_handler", &field_handler,
                 py::arg("param1"));
   }
   

Thread Management in Handlers

Handlers are executed in parallel using OpenMP, with multiple threads processing different cells simultaneously. Any data structures that are modified during Handler execution must be thread-safe to avoid race conditions and data corruption.

Example: threadHandler struct

Random number generation (RNG) requires maintaining state (the current seed value) that changes after each call, otherwise multiple OpenMP threads access the same RNG state simultaneously. The solution is to provide each OpenMP thread with its own independent RNG instance via thread-local storage and deterministic seeding based on CI.rngSeed.

struct threadHandler {
    mt19937 rng;
    std::uniform_real_distribution<double> U1;
    std::normal_distribution<double> N1;
    double working_array[100];  // Temporary computation buffers

    threadHandler(): U1(0, 1.0), N1(0, 1.0) {}
    double random() { return U1(rng); }
};

static vector<threadHandler> Thread;

void Handler(...) {
    cellInterface CI(...);
    Thread.resize(omp_get_max_threads());
    threadHandler &th = Thread[CI.threadNum];  // Get thread-local instance
    th.rng.seed(CI.rngSeed);  // Deterministic per-cell seeding

    // Use th.random(), th.working_array[i], etc.
}

Example from extension.cpp:

See extension.cpp for a minimal example using threadHandler.

Build extension (C/C++)

Two build strategies are available for extensions: local build which creates a local executable, or integrated build which makes the extension callable within the π-PIC framework.

Local build

1. Prepare build environment (from a working directory):

   mkdir my_extension && cd my_extension
   git clone https://github.com/hi-chi/pipic.git
   cd pipic
   
   # Optional: create and switch to development branch
   git checkout -b my_extension
   
   # Optional: reinstall pipic
   pip uninstall pipic
   python3 -m pip install .
   

2. Set up extension folder:

   cd src/extensions
   mkdir my_extension && cd my_extension
   

3. Copy necessary files:

   git clone https://github.com/pybind/pybind11
   cp ../../primitives.h .
   cp ../../interfaces.h .
   cp ../../ensemble.h .    # For particle management
   cp ../../services.h .    # For logging
   cp ../../CMakeLists.txt .
   

4. Configure CMakeLists.txt (replace pipic references):

   Change from:

   set(pipic
       pipic.cpp)
   pybind11_add_module(_pipic${pipic})
   

   To:

   set(my_extension
       my_extension.cpp)
   pybind11_add_module(_my_extension ${my_extension})
   

   Optional: Remove FFTW if not needed:

   # Remove these lines from CMakeLists.txt:
   # find_package(PkgConfig REQUIRED)
   # pkg_search_module(FFTW REQUIRED fftw3 IMPORTED_TARGET)
   # Remove #include <fftw3.h> from primitives.h
   

5. Add your extension files: my_extension.cpp. See extension.cpp for a complete example.

6. Compile locally:

   cmake .
   make
   

7. Test in Python:

   import my_extension
   from pipic import addressof
   
   # Initialize extension
   handler_ptr = my_extension.handler(param1=1.0, param2=42)
   field_handler_ptr = my_extension.field_handler(param1=0.1)
   
   # Register handlers
   sim.add_handler(name=my_extension.name,
                  subject='electron',
                  handler=handler_ptr,
                  field_handler=field_handler_ptr)
   

Practical Example:

The tutorial plasma_oscillation_extension_development.ipynb demonstrates developing an absorbing boundary extension.

Integrated build

After successful development and testing:

1. Update pybind11 module name (add underscore prefix to avoid conflicts):

   PYBIND11_MODULE(_my_extension, object) {
   

2. Organize files in src/extensions/my_extension/:

   src/extensions/my_extension/
   ├── my_extension.cpp
   └── my_extension.h       (if needed)
   

3. Create example my_extension_test.py in examples/ folder:

   from pipic.extensions import my_extension
   import pipic
   from pipic import consts
   
   sim = pipic.init(...)
   handler_ptr = my_extension.handler(...)
   sim.add_handler(name=my_extension.name,
                  subject='electron',
                  handler=handler_ptr)
   

4. Register in pipic/extensions/__init__.py

5. Test clean build:

   pip uninstall pipic
   python3 -m pip install .
   python3 examples/my_extension_test.py
   

6. Commit and push:

   git add src/extensions/my_extension/
   git add examples/my_extension_test.py
   git add pipic/extensions/__init__.py
   git commit -m "Add my_extension: ..."
   git push origin my_extension
   # Create pull request on GitHub
   

7. Document in EXTENSIONS.md with physics description, parameters, and usage.

Performance Tips

• Use C/C++ for performance-critical extensions; numba callbacks are slower but more convenient for prototyping.

• Minimize allocations inside hot loops (particle/cell processing).

• Use thread-local RNG via CI.rngSeed for deterministic random number generation.

• Disable threading during debugging: sim.advance(..., use_omp=False).