Simulation contracts¶

These interfaces define the reusable simulation lifecycle. Read the architecture and application pipeline for ownership and control flow before using the exact declarations.

`pfc::Model`¶

class Model

The Model class represents the physics model for simulations in OpenPFC.

The Model class is responsible for introducing the physics to the simulation model. Users can override the initialize and step functions to define their own functionality and implement specific physics for their simulations.

The initialize function is called at the beginning of the simulation and is used to allocate arrays, pre-calculate operators used in time integration, and perform other necessary initialization tasks.

The step function is called sequentially during the time integration process. OpenPFC currently uses a semi-implicit time integration scheme, where the linear part is solved implicitly and the non-linear part is solved explicitly.

Public Functions

~Model() = default: Destroy the Model object.

Model(const Model&) = delete: Disable copy constructor.

Model &operator=(const Model&) = delete: Disable copy assignment operator.

inline bool is_rank0() const noexcept

Check if current MPI rank is 0.

Useful for conditional output (only rank 0 prints).

if (model.is_rank0()) {
    std::cout << "Status message\n";
}

Note

This function is thread-safe and has zero overhead (defined in header)

Returns:: true if rank is 0, false otherwise

inline MPI_Comm mpi_comm() const noexcept

MPI communicator used for rank-0 checks (is_rank0())

Matches the communicator passed at construction; align with Simulator and FFT when using a subcommunicator or non-default world communicator.

inline const World &get_world() const noexcept

Get the world object associated with the model.

Returns:: Reference to the World object

inline const ModelFieldRegistry &fields() const noexcept

Read-only access to the named field registry.

Use for helpers, tests, or custom get_allocated_memory_bytes() that should account for registered_field_storage_bytes(), and for introspection via list_field_names().

inline virtual void prepare_for_field_modifiers()

Hook run by the Simulator just before field modifiers (initial / boundary conditions) or result writers touch the registered fields.

Default: no-op. A model whose authoritative field lives on a device (a GPU model with a registered host mirror) overrides this to copy device -> host so modifiers and writers observe current data. Paired with finalize_after_field_modifiers().

Note

Fixes audit 4.1: previously the App/Simulator path applied JSON initial conditions to a GPU model’s host mirror and never propagated them to the device, so App-driven GPU runs integrated from an unseeded device buffer. The Simulator now brackets modifier application and result writing with these hooks.

inline virtual void finalize_after_field_modifiers(): Hook run by the Simulator just after field modifiers have mutated the registered (host) fields, to propagate changes back to the model’s authoritative storage (host -> device). Default: no-op.

inline bool has_real_field(std::string_view field_name) const noexcept

Check if the model has a real-valued field with the given name.

Parameters:: field_name – Name of the field to check
Returns:: True if the field exists, False otherwise

inline const RealField &get_real_field(std::string_view name) const

Retrieve a registered real-valued field by name (const version)

Returns a const reference to a previously registered field for reading.

See also

get_real_field() for non-const version

Parameters:: name – Name of the field to retrieve
Throws:: std::out_of_range – if field name not registered
Returns:: Const reference to the RealField object

inline const ComplexField &get_complex_field(std::string_view name) const

Retrieve a registered complex-valued field by name (const version)

Returns a const reference to a previously registered complex field for reading.

See also

get_complex_field() for non-const version

Parameters:: name – Name of the field to retrieve
Throws:: std::out_of_range – if field name not registered
Returns:: Const reference to the ComplexField object

inline void add_field(const std::string &name, RealField &field)

Add a field to the model.

Parameters:

name – Name of the field
field – Reference to the RealField object representing the field

inline void add_field(const std::string &name, ComplexField &field)

Add a field to the model.

Parameters:

name – Name of the field
field – Reference to the ComplexField object representing the field

inline virtual Field &get_field()

Get a reference to the default primary unknown field.

Deprecated:: Prefer get_real_field("default") or another registered name. The legacy free function pfc::get_field(Model&) was removed from OpenPFC; use get_real_field(model, "default") (or the overload on Simulator) at new call sites.

Note

Out-of-tree models that still override this hook for a single primary field should migrate to explicit register_real_field("default", …) (or a domain-specific name) and call sites that use get_real_field. Behavior is rank-local like other field accessors; MPI consistency of the underlying DiscreteField data is unchanged.

Returns:: Reference to the RealField called “default”

`pfc::Simulator`¶

class Simulator

Base class for time-stepping simulations with forward Euler integration.

The Simulator class provides a framework for time-dependent field evolution using an embedded forward Euler time-integration method. The class owns specific lifecycle stages (pre-step preparation, RHS evaluation, post-step updates, output generation, checkpointing) and coordinates boundary condition application and MPI halo exchange with these stages.

Time-integration assumptions

Embedded integrator: forward Euler (first-order explicit)
One RHS evaluation per time step
No substepping; each call to advance() evolves fields by exactly dt
Time step size dt is fixed and enforced by an external scheduler; the simulator does not adapt dt
Field update follows: phi_next = phi + dt * RHS(phi, t)

Lifecycle stage ownership

The Simulator base class owns the following stages (implemented via virtual methods that subclasses may override):

Pre-step preparation: anything required before RHS evaluation (e.g., loading boundary conditions)
RHS evaluation: computes the right-hand side of the evolution equation
Post-step updates: applies time-stepped field updates after RHS is computed
Output generation: writes field data, diagnostics, or visualization files
Checkpointing: saves simulation state for restart

Boundary/halo synchronization

Boundary conditions and MPI halo exchanges must be synchronized with the integration stages. The expected ordering is:

Apply boundary conditions to fill halo regions
Compute RHS using the synchronized field state
Post-step updates complete before the next synchronization

Contract for substituting alternative integrators

To swap in a different time-integration scheme (e.g., Runge-Kutta), subclasses must:

Override advance() to implement the multi-stage algorithm
Provide intermediate storage for stage values if needed
Call boundary synchronization at appropriate stage boundaries
Maintain the same pre- and post-step hooks for compatibility with output and checkpoint scheduling
Preserve the contract that dt is externally fixed (no adaptive dt)
Document the new scheme’s stage count, RHS evaluations per step, and any additional synchronization points

Note

This base class does not own the time-stepping loop or schedule; external code determines when advance() is called and how to drive the simulation to completion or termination.

Public Functions

inline Simulator(Model &model, Time &time, MPI_Comm mpi_comm = MPI_COMM_WORLD)

Construct a new Simulator object.

Parameters:

model – The model to simulate.
time – The time object to use for simulation.
mpi_comm – Communicator passed to field modifiers (MPI-IO collectives, etc.)

inline Model &get_model()

Get the model object.

Returns:: Model&

inline const World &get_world() const noexcept

Get the decomposition object.

Get the world object

Returns:: const Decomposition&
Returns:: const World&

inline fft::IFFT &get_fft()

Get the FFT object.

Returns:: fft::IFFT&

inline Time &get_time()

Get the time object.

Returns:: Time&

inline bool is_rank0() const noexcept: Rank 0 in mpi_comm() (same communicator passed to field modifiers)

inline const ResultsWriterMap &results_writers() const noexcept

Inspect registered results writers (read-only; for tests and tooling)

Writers are still added only via add_results_writer.

See also

pfc::results_writers

inline const std::vector<std::unique_ptr<FieldModifier>> &get_initial_conditions() const

Gets the initial conditions of the simulation.

This function returns a const reference to the vector of unique pointers to FieldModifier objects that represent the initial conditions of the simulation.

Returns:: A const reference to the vector of unique pointers to FieldModifier objects.

inline const std::vector<std::unique_ptr<FieldModifier>> &get_boundary_conditions() const

Gets the boundary conditions of the simulation.

This function returns a const reference to the vector of unique pointers to FieldModifier objects that represent the boundary conditions of the simulation.

Returns:: A const reference to the vector of unique pointers to FieldModifier objects.

inline void set_result_counter(int result_counter)

See also

pfc::set_result_counter

inline int get_result_counter() const

See also

pfc::get_result_counter

inline void write_results()

Increment the result counter and invoke every registered writer.

Implemented by write_scheduled_simulator_results. For tests that only need the write loop, call pfc::write_results_for_registered_fields with a stub model and writer map instead of wiring a full Simulator.

See also

write_scheduled_simulator_results

See also

pfc::write_results

See also

results_writers()

See also

pfc::write_results_for_registered_fields

inline void begin_integrator_step()

Prologue of one integrator step (same ordering as step())

When Time::get_increment() == 0, applies initial conditions, boundary conditions, and optionally writes results if Time::do_save(). Then always calls Time::next() and applies boundary conditions at the new time.

Call this once per iteration before the physics update, then invoke Model::step (or a model-specific step(simulator, model) overload), then call end_integrator_step().

See also

end_integrator_step()

See also

step_with_physics()

See also

step()

inline void end_integrator_step()

Epilogue of one integrator step: write results if at a save point.

See also

begin_integrator_step()

template<class PhysicsFn> inline void step_with_physics(PhysicsFn &&physics_fn)

One full step with a custom physics body (same ordering as step())

Equivalent to begin_integrator_step(); physics_fn(); end_integrator_step(); with physics_fn typically calling pfc::step(model, time.get_current()) or a model-specific overload such as step(simulator, concrete_model).

Template Parameters:: PhysicsFn – nullary callable

`pfc::Time`¶

class Time

Public Functions

inline Time(const std::array<double, 3> &time)

Construct a new Time object with the specified time interval and default save interval.

The save interval is set to the same value as the time step.

Parameters:: time – An array containing the start time, end time, and time step in that order

inline Time(const std::array<double, 3> &time, double saveat)

Construct a new Time object with the specified time interval and default save interval (euler integrator).

The save interval is set to the same value as the time step.

Parameters:: time – An array containing the start time, end time, and time step in that order

inline double get_t0() const

Get the start time.

Returns:: The start time

inline double get_t1() const

Get the end time.

Returns:: The end time

inline double get_dt() const

Get the time step.

Returns:: The time step

inline void set_dt(double dt)

Set the time step.

Sets a new time step size for adaptive time-stepping. The time step must be positive to ensure meaningful time integration.

This method enables runtime adjustment of the recommended / policy dt. It does not change accepted simulation time; adaptive drivers must advance time via commit_attempt (or fixed-step next).

See also

get_dt() - query the current time step

See also

set_increment() - fixed-dt reconstruction helper (not adaptive commit)

Note

Changing dt does not rewrite past accepted time. Use begin_attempt / commit_attempt for clipped attempts.

Parameters:: dt – [in] The new time step size (must be > 0)
Throws:: std::invalid_argument – If dt <= 0
Post:: get_dt() returns the new dt value
Post:: get_accepted_time() is unchanged

inline int get_increment() const

Get the current time increment.

Returns:: Current increment

inline int get_step_count() const

Get the number of steps completed so far.

Self-documenting alternative to get_increment() for integrator authors querying progress (checkpoint/restart, adaptive step rejection logic, output scheduling): “increment” reads as a small delta, not a count, whereas “step count” states the intent directly. Returns the same underlying value as get_increment().

Returns:: Current step count

inline void increment_step_success()

Record a successful adaptive step attempt.

Increments the accepted-attempt counter. Unlike next(), this does not advance simulation time; it only updates attempt statistics for adaptive time-stepping algorithms.

Note

When used with TimeStateGuard, call this before commit() (or outside an uncommitted guard); otherwise the guard restores the counter on destruction.

Post:: get_accepted_steps() returns previous value + 1

inline void increment_step_rejection()

Record a rejected adaptive step attempt.

Increments the rejected-attempt counter. Unlike next(), this does not advance simulation time; it only updates attempt statistics.

Note

TimeStateGuard also restores accepted/rejected counters unless commit() was called. Call this after an uncommitted guard’s scope ends (or outside the guard) so the rejection count is kept.

Post:: get_rejected_steps() returns previous value + 1

inline int get_accepted_steps() const

Get the number of accepted adaptive step attempts.

Returns:: Accepted step count (starts at 0)

inline int get_rejected_steps() const

Get the number of rejected adaptive step attempts.

Returns:: Rejected step count (starts at 0)

inline int get_stage() const

Get the current stage index within this time step.

Stage tracking is independent of get_increment(): stages range from 0 to get_stage_count() - 1 and describe progress within a single multi-stage step (e.g. RK2/RK4), while the increment counts completed time steps.

Returns:: Current stage index (0-based)

inline void set_stage(int stage)

Set the current stage index within this time step.

Multi-stage steppers should call this before computing each stage.

Parameters:: stage – [in] New stage index (must satisfy 0 <= stage < get_stage_count())
Throws:: std::invalid_argument – If stage is out of range

inline int get_stage_count() const

Get the total number of stages in the current time step.

Returns:: Stage count (always >= 1)

inline void set_stage_count(int stage_count)

Set the total number of stages for this time step.

Multi-stage steppers should call this once during initialization to configure how many stages each step has (e.g. 2 for RK2, 4 for RK4).

Parameters:: stage_count – [in] New stage count (must be >= 1)
Throws:: std::invalid_argument – If stage_count < 1

inline double get_accepted_time() const noexcept

Get the accepted simulation time (read-only)

Returns the stored accepted clock. Unchanged for the duration of an active attempt (begin_attempt … commit_attempt / reject_attempt).

Returns:: Accepted simulation time

inline double get_current() const

Get the current (accepted) simulation time.

Returns the stored accepted time (m_accepted_time), clamped to t1 if needed. Prefer get_accepted_time for adaptive drivers; this alias preserves the historical name used by Simulator and writers.

See also

get_accepted_time() - explicit accepted-time accessor

See also

get_increment() - get the number of steps taken

See also

next() - advance to next time step (fixed dt)

See also

done() - check if t_current >= t1

Note

Accepted time advances on next or commit_attempt only. set_dt does not rewrite this value.

Returns:: Current simulation time, clamped to [t0, t1]

inline double get_saveat() const

Get the time interval for saving data.

Returns:: The save interval

inline IntegratorMethod method() const noexcept

Get the integrator method.

Returns:: The integrator method

inline void set_increment(int increment)

Set the current time increment (fixed-dt reconstruction helper).

Sets the committed step count and reconstructs accepted time as min(t0 + increment * dt, t1). This preserves restart / rewind helpers that assume a constant dt. It is not the adaptive commit path — use commit_attempt after a clipped attempt instead.

The increment must be non-negative so that get_current() stays in [t0, t1] and done() / do_save() remain meaningful.

Parameters:: increment – The current time increment (number of completed steps from t0, i.e. same convention as after next()).
Throws:: std::invalid_argument – if increment < 0
Post:: get_increment() == increment
Post:: get_accepted_time() == min(t0 + increment * dt, t1)

inline void set_saveat(double saveat)

Set the time interval for saving data.

Parameters:: saveat – The save interval

inline void next()

Advance to the next time step (fixed dt)

Advances accepted time by dt (clamped to t1) and increments the step counter by 1. For adaptive clipped intervals, use begin_attempt / commit_attempt instead.

See also

get_increment() - query current step number

See also

get_current() - accepted simulation time

See also

done() - check completion status

See also

commit_attempt() - advance by a clipped attempted interval

Post:: get_increment() returns previous value + 1
Post:: get_accepted_time() returns previous value + dt (clamped to t1)

inline double clip_attempt_dt(double candidate_dt) const

Clip a candidate step interval against terminal and output bounds.

Returns an attempted dt such that accepted_time + attempted_dt does not exceed t1. When saveat > 0, also does not pass the next output-alignment time without landing on it.

Alignment uses the same tolerance family as do_save (1e-9): next_save = ceil((accepted + 1e-9) / saveat) * saveat. When both the terminal bound and saveat constrain the step, the most restrictive (minimum) wins. If next_save >= t1, only the terminal bound matters. When saveat <= 0, output-alignment clipping is skipped.

Parameters:: candidate_dt – [in] Proposed step size (must be > 0)
Throws:: std::invalid_argument – If candidate_dt <= 0
Returns:: Clipped attempted interval

inline bool attempt_active() const noexcept: Whether an attempt transaction is currently open.

inline double get_attempted_dt() const

Clipped interval for the active attempt.

Throws:: std::logic_error – If no attempt is active

inline void begin_attempt(double candidate_dt)

Begin an attempt: clip candidate_dt and leave accepted time unchanged.

Parameters:

candidate_dt – [in] Proposed step size (must be > 0)

Throws:

std::logic_error – If an attempt is already active
std::invalid_argument – If candidate_dt <= 0 (via clip_attempt_dt)

Post:

attempt_active() is true

Post:

get_accepted_time() is unchanged

Post:

get_attempted_dt() == clip_attempt_dt(candidate_dt)

inline void commit_attempt()

Commit the active attempt: advance accepted time by attempted_dt.

Does not auto-call increment_step_success (counters stay caller-owned).

Throws:: std::logic_error – If no attempt is active
Post:: get_accepted_time() advanced by the attempted interval (clamped to t1)
Post:: get_increment() increased by 1
Post:: attempt_active() is false

inline void reject_attempt()

Reject the active attempt without changing accepted time.

Does not auto-call increment_step_rejection (counters stay caller-owned).

Throws:: std::logic_error – If no attempt is active
Post:: get_accepted_time() and get_increment() unchanged
Post:: attempt_active() is false

inline operator double() const

Conversion operator to retrieve the current time as a double value.

Returns:: The current time as a double value

Friends

inline friend std::ostream &operator<<(std::ostream &os, const Time &t)

Overloaded stream insertion operator to print the Time object.

Parameters:

os – The output stream
t – The Time object to be printed

Returns:

The output stream

`pfc::FieldModifier`¶

class FieldModifier

Subclassed by pfc::Constant, pfc::FileReader, pfc::FixedBC, pfc::MovingBC, pfc::RandomSeeds, pfc::SeedGrid, pfc::SingleSeed

Public Functions

inline void set_field_names(std::vector<std::string> names): Declare every field this modifier may write (for Simulator checks).

inline virtual void set_mpi_comm(MPI_Comm): Optional MPI communicator for modifiers that use collectives (MPI-IO, reductions). Default is a no-op; FileReader and MovingBC override.

inline virtual void apply(const SimulationContext &simulation_context, Model &model, double time)

Apply the field modification with explicit simulation context.

The simulator invokes this overload so modifiers can use simulation_context (e.g. mpi_comm()) without relying solely on set_mpi_comm(). The default implementation ignores the context and calls apply(Model&, double).

Modifiers that need MPI collectives should override this method and use simulation_context.mpi_comm(). They may still override apply(Model&, double) to wrap a default SimulationContext for direct/test calls.

Note

Contract (substitutability): Production runs use this overload. Prefer implementing one core body (e.g. a private apply_impl(...)) and having both apply(SimulationContext,...) and apply(Model&,double) forward to it so direct unit tests and the simulator stay consistent. If you override only apply(Model&,double), the context overload’s default still delegates there—ensure any MPI-aware logic is reachable from that path or override the context overload as well.

virtual void apply(Model &model, double time) = 0

Apply the field modification to the model (pure virtual)

This is the main interface method that derived classes must implement to define their modification logic. The method receives full mutable access to the Model and current simulation time, allowing arbitrary modifications.

Implementation Responsibilities:

Retrieve field(s) via get_real_field(model, name) or get_complex_field(model, name)
Access geometry via pfc::get_world(model) and pfc::get_fft(model)
Modify field values according to modifier’s purpose
Handle MPI parallelism (operate on local subdomain)

Typical Implementation Pattern:

void apply(pfc::Model& model, double time) override {
  // 1. Get field to modify
  auto& field = get_real_field(model, get_field_name());

  // 2. Get geometry information
  const auto& world = pfc::get_world(model);
  const auto& fft = pfc::get_fft(model);
  auto inbox = pfc::fft::get_inbox(fft);

  // 3. Loop over local subdomain
  int idx = 0;
  for (int k = inbox.low[2]; k <= inbox.high[2]; k++) {
    for (int j = inbox.low[1]; j <= inbox.high[1]; j++) {
      for (int i = inbox.low[0]; i <= inbox.high[0]; i++) {
        // Compute modification based on position and/or time
        auto pos = pfc::world::to_coords(world, Int3{i, j, k});
        field[idx++] = compute_value(pos, time);
      }
    }
  }
}

See also

Model::get_real_field() for field access

See also

get_world(const Model&) for domain geometry

See also

get_fft(Model&) for subdomain bounds

Note

For initial conditions, time is typically 0.0

Note

For boundary conditions, time reflects current simulation time

Note

Method is called on every MPI rank; each rank operates on its subdomain

Note

The simulator’s entry point is apply(SimulationContext,...); see its documentation for how to keep this overload and that one equivalent.

Warning

Ensure modifications maintain physical correctness and don’t violate model invariants (e.g., mass conservation if required)

Parameters:

model – Mutable reference to the Model containing fields to modify
time – Current simulation time (useful for time-dependent BCs)

Pre:

Model must have the field specified by get_field_name() registered

Post:

Field values are modified according to modifier’s logic

virtual ~FieldModifier() = default

Destructor for the FieldModifier class.

The destructor is declared as default, allowing proper destruction of derived classes.

`pfc::ResultsWriter`¶

class ResultsWriter

Subclassed by pfc::BinaryWriter, pfc::VTKWriter

Public Functions

virtual ~ResultsWriter() noexcept(false) = default

Virtual destructor.

Marked noexcept(false) to allow subclasses to throw on MPI cleanup failures (fail-closed policy).

virtual MPI_Status write(int increment, const ComplexField &data) = 0

Write a complex-valued field to file at specified time step.

Writes the local portion of a ComplexField (complex doubles) to the output file. Useful for storing Fourier coefficients or k-space data.

See also

FFT::forward() - produces ComplexField from RealField

Note

Complex field size is typically ~50% of real field size (r2c symmetry).

Parameters:

increment – [in] Time step or frame number
data – [in] Local complex field data

Returns:

MPI_Status Information about the write operation

Simulation contracts¶

pfc::Model¶

pfc::Simulator¶

pfc::Time¶

pfc::FieldModifier¶

pfc::ResultsWriter¶

`pfc::Model`¶

`pfc::Simulator`¶

`pfc::Time`¶

`pfc::FieldModifier`¶

`pfc::ResultsWriter`¶