Skip to content

Programming Guide

PSuM Programming Guide

If you are not yet familiar with basic C++ templates, lambdas, and SYCL syntax, you can start with C++ in PSuM. This document only covers "how to write programs using PSuM"; for internal implementation details, see Framework Design.

This document is intended for users who can already compile PSuM examples and are ready to write their own simulations. If your environment is not yet set up, consider Docker Deployment first (if you have control over your environment, such as a personal computer, this is much simpler than manual setup). If Docker is not suitable, read Dependency Installation and Build System Instructions to configure manually.

Environment Setup

Before writing code, you typically need to load the project-generated environment variables:

source env_load.sh

env_load.sh is usually generated by env_scan.sh and sets up paths for AdaptiveCpp, Eigen, CUDA, UMFPACK, and other dependencies. To disable CUDA or adjust the architecture, override USE_CUDA, CUDA_ARCH, and other options in config.mk.local rather than modifying config.mk.

To verify the environment is working, you can run:

cd test
bash runCheck.sh

This test typically takes a few minutes; please be patient.

Alternatively, enter an example directory and run make. The code in this guide assumes the above environment is already loaded.

Code Entry Point and Namespaces

Most user programs only need to include the main entry point:

#include <psum/psum.hpp>

using namespace psum::prelude;

psum::prelude aggregates the namespaces of commonly used modules, including tags, serialization, particle containers, fields, boundaries, field solvers, collisions, random numbers, and various SYCL utilities. Using prelude keeps your code concise.

If you prefer explicit control over dependencies, you can include individual module entry points:

#include <psum/field.hpp>
#include <psum/particle_container.hpp>
#include <psum/field_solver.hpp>

and use the corresponding module namespaces:

using namespace psum::field;
using namespace psum::particle_container;
using namespace psum::field_solver;

You should prefer the public entry points under include/psum/. Headers under src/ reflect the project's internal organization. While many implementations currently exist in header-only form, they are not necessarily stable user-facing entry points.

Hello PSuM!

Here is a minimal program: create a compute device, write data into a device-side array, and copy it back to the host.

#include <psum/psum.hpp>
#include <iostream>

using namespace psum::prelude;

int main() {
    sycl::queue q{sycl::default_selector_v};
    std::cout << q.get_device().get_info<sycl::info::device::name>() << std::endl;

    device_vector<double> values(q, 16);
    values.for_each([&](sycl::handler& h) {
        return [=](double& x) {
            x = 3.14;
        };
    });

    auto host_values = values.to_host();
    std::cout << host_values[0] << std::endl;
}

for_each uses the "double lambda" pattern common in PSuM: the outer lambda prepares resources on the host side, and the inner lambda executes in parallel on the device. The inner lambda should typically use [=] value capture. For a more detailed explanation, see C++ in PSuM.

Defining Particle Types

PSuM uses a tag system to describe particle attributes. Use tagged_struct:

using Particle = tagged_struct<
    tag_bind<property::position, Eigen::RowVector3d>,
    tag_bind<property::velocity, Eigen::RowVector3d>,
    tag_bind<property::weight, double>,
    tag_bind<property::random_seed, uint32_t>
>;

Access attributes using get<Tag>(particle):

Particle p;
get<property::position>(p) = Eigen::RowVector3d(0.0, 0.0, 0.0);
get<property::velocity>(p) = Eigen::RowVector3d::Zero();
get<property::weight>(p) = 1.0;
get<property::random_seed>(p) = 1;

Common built-in tags include:

  • property::position
  • property::velocity
  • property::acceleration
  • property::mass
  • property::charge
  • property::weight
  • property::random_seed
  • property::species_id
  • property::internal_energy

To define a custom tag, simply inherit from abstract_tag and provide tag_name:

struct init_position : tag::foundation::abstract_tag {
    inline const static std::string tag_name = "init_position";
};

If a tag needs to restrict the bound type, you can implement check<T>(). The compile-time lookup and type verification details of the tag system are internal mechanisms; see Framework Design.

A common pitfall is name conflicts between local variables and tag names. For example, a local variable mass will shadow property::mass. In practice, it is recommended to keep the property:: prefix.

Using device_vector

device_vector<T> is a device-side dynamic array, similar to std::vector<T>, but with data stored in SYCL device memory.

Common construction methods:

std::vector<Particle> host_particles(100);
device_vector<Particle> from_host(q, host_particles);

device_vector<Particle> buffer(q, 1000);  // reserved capacity: 1000

Data transfer between host and device:

auto host_copy = from_host.to_host();
from_host.copy(host_particles);

Device-side iteration:

from_host.for_each([&](sycl::handler& h) {
    return [=](Particle& p) {
        get<property::velocity>(p) = Eigen::RowVector3d::Zero();
    };
});

To fill new elements on the device side, you must first obtain an accessor in the outer lambda:

device_vector<Particle> selected(q, 1000);

from_host.for_each([&](sycl::handler& h) {
    auto out = selected.get_access(h);
    return [=](Particle& p) {
        if (get<property::position>(p).x() > 0.0) {
            out.push_back(p);
        }
    };
});

push_back is thread-safe, but capacity does not grow automatically. If the number of writes exceeds the reserved capacity, the container enters an overflow state, and subsequent host-side operations will surface this error as early as possible.

Using particle_group

particle_group is the primary particle container in PSuM. It is built on top of device_vector and additionally supports lazy deletion and space reuse.

Use particle_group if:

  • Data requires frequent deletion
  • Device-side filling is not needed

When defining a particle group, you must specify the particle type and a validator:

using ParticleGroup = particle_group<Particle, pos_x_nan_is_invalid>;

ParticleGroup particles(q);
particles.insert(host_particles);

pos_x_nan_is_invalid determines whether a particle is valid based on whether position.x() is NaN. When iterating over a particle_group, invalid particles are automatically skipped:

particles.for_each([&](sycl::handler& h) {
    return [=](Particle& p) {
        if (get<property::position>(p).x() > 1.0) {
            ParticleGroup::validator::make_invalid(p);
        }
    };
});

Deletion only marks particles as invalid and does not immediately move memory. Afterwards, you can use:

particles.compress();
particles.shrink(particles.size() * 1.2);

size() returns the number of valid particles, and capacity() returns the allocated capacity. The dual-array structure and compression algorithm of particle_group are described in Framework Design.

Using Grids and Fields

The field system is responsible for storing physical quantities such as charge density, electric potential, and electric field on a regular grid.

Creating a grid:

grid2D grid({0.0, 0.0}, {1.0, 1.0}, {64, 64});

There are two common types of fields:

  • device_field: Device-side field, suitable for reading and writing inside kernel functions.
  • host_field: Host-side field, more convenient for direct manipulation on the host.

For particle-field interaction scenarios, you should generally use device_field.

Example:

node_field2D<double> rho(q, grid);
node_field2D<double> phi(q, grid);

When using fields, pay attention to three things: dimension, position type (node or cell-centered), and numeric type.

Interpolation and Deposition

In PIC, there are two fundamental operations between particles and fields:

  • interp: Interpolate from the grid field to particle positions.
  • add_back: Deposit particle weights back onto the grid.

Core functions provided by PSuM:

auto value = interp(pos, field_acc);
auto grad = interp_diff(pos, field_acc);
add_back(pos, weight, field_acc);

These functions are typically used inside particle iteration:

particles.for_each([&](sycl::handler& h) {
    auto rho_acc = rho.get_access(h);
    return [=](Particle& p) {
        auto pos = get<property::position>(p);
        auto w = get<property::weight>(p);
        add_back(pos, w, rho_acc);  // use rho_acc, not rho!
    };
});

Different field types support different interpolation/deposition methods.

If you only need nearest-neighbor versions, you can use interp_nearest and add_back_nearest. For details on interpolation kernels, tensor product expansion, and atomic deposition, see Framework Design.

Solving the Poisson Equation

A typical electrostatic PIC program can be summarized as:

for each time step:
    rho.setZero()
    particles.for_each(deposit charge to rho)
    solver.solve(phi, rho)
    particles.for_each(push particles by -interp_diff(phi); apply boundary conditions)
    execute collisions

Particle boundaries can sometimes be quite simple and written directly in the push function; particle collisions can sometimes be ignored, but the field solver is essential. Poisson solvers are organized by dimension:

  • Poisson_solver_1d
  • Poisson_solver_2d
  • Poisson_solver_3d

The typical usage flow is:

1. Create a grid
2. Create fixed (Dirichlet) or mixed (Robin) boundary lists
3. solver.init(grid, boundaries, ...)
4. Call solver.solve(phi, source) each time step

Fixed boundaries (Dirichlet) directly specify potential values; mixed boundaries (Robin) handle conditions involving both function values and normal derivatives. The solver can use the default native backend, or you can select an alternative backend during initialization.

Matrix assembly, ghost node elimination, and the backend plugin mechanism are internal implementation details; see Framework Design.

Handling Particle Boundaries

The particle boundary system handles the intersection of particle trajectories with geometric surfaces. The usage approach is:

1. Load or generate triangle mesh geometry
2. Assign material types to triangles
3. Create a boundary_router
4. After particle push, call router_acc.deal(...) with old and new positions

Typical call pattern inside a kernel:

particles.for_each([&](sycl::handler& h) {
    auto router_acc = router.get_access(h);
    return [=](Particle& p) {
        auto old_pos = get<property::position>(p);
        // update position here
        auto new_pos = get<property::position>(p);
        router_acc.deal(
            old_pos.x(), old_pos.y(), old_pos.z(),
            new_pos.x(), new_pos.y(), new_pos.z(),
            p
        );
    };
});

Built-in material behaviors include absorption, specular reflection, diffuse reflection, and Maxwell reflection. Many reflection-type behaviors also require the particle to have a random_seed attribute.

Handling Collisions

PSuM's collision module mainly includes:

  • MCC: Monte Carlo Collision between particles and background species.
  • DPMCC: Deferred Pairing Collision that requires explicit pairing.

Users need to understand three concepts:

  • Collision model: Describes the incident species and collision channels.
  • Species context: Binds particle groups, background density fields, and new particle buffers.
  • Execution function: Invokes the collision model within a time step and cleans up buffers afterward.

Simplified flow:

prepare species contexts
execute_mcc_model<Model>(ctx, dt)
execute_all_clean_buffers<SpeciesTuple>(ctx)

Custom collision cross-sections, channel selection, and pairing algorithms are typically too lengthy for a main tutorial. Complete usage examples should be provided through examples.md or specific application documentation.

Saving, Restoring, and Configuration

Basic save/load uses mas_file:

mas_file fp("state.mas", mas_file::replaceMode);
save(fp, "particles", particles);
load(fp, "particles", particles);

When managing multiple objects, you can use object_manager. It is well-suited for checkpoint/restart: objects are registered with the object_manager, and on the next run they are restored from an existing MAS file.

Configuration files can be loaded with json_loader:

json_loader loader;
loader.load_json("config.json");
double dt = loader["solver"].obj<double>("dt");

Serialization of device_vector and particle_group involves device information; you cannot assume unconditional exchange between different hardware.