OpenPFC / Build a simulation
Let the error choose the next step
Compare memoryless and PI control, keep rejected candidates separate, and land on a known change in decay rate.
Download runnable lesson (.zip)
Includes the model, inputs, required headers and a short run guide.
A fixed step must be small enough for the most demanding part of a run. An adaptive method asks how much error the current step introduces. This refactors the earlier adaptive example around the deliberately simple model \(y'=-k(t)y\), with \(y(0)=1\), \(k=50\) before \(t_s=0.15\) and \(k=1\) afterwards. Its exact solution is
\[y(t)=\exp[-50\min(t,t_s)-\max(0,t-t_s)]. \tag{1}\]
Two approximations from shared stages
The native Bogacki–Shampine factory supplies a third-order result and an embedded second-order result. Their difference is \(O(h^3)\). The controller scales it component by component:
\[r=\max_i\frac{|e_i|}{\mathrm{atol}+\mathrm{rtol}\max(|y_i|,|y_i^*|)}. \tag{2}\]
All MPI ranks agree on the maximum. A step is admitted when the metric is finite and at most one. The memoryless next-step factor is proportional to \(r^{-1/3}\). The native PI controller also uses the previous accepted metric; its exact gains and safety convention are documented in the public controller header. A local error tolerance is not a guarantee on the final global error. The control-theory viewpoint is developed by Söderlind (2002); the concrete gain and safety convention used here follows the OpenPFC contract above.
A candidate is not an accepted state
while (!time.done()) {
if (++attempts > 10000) throw std::runtime_error("Attempt limit reached");
const double t0 = time.get_accepted_time();
double proposed = time.get_dt();
if (t0 < transition) proposed = std::min(proposed, transition - t0);
rate = t0 < transition ? 50 : 1; // One smooth phase per attempted interval.
time.begin_attempt(proposed);
const double dt = time.get_attempted_dt();
const auto candidate = stepper.attempt(t0, dt, y);
const auto decision = controller.decide_from_embedded_error(
dt, stepper.error(), y, stepper.u_high(), comm);
if (decision.accepted) pfc::sim::steppers::commit_step_attempt(y, candidate);
controller.apply(time, decision);The stepper keeps its candidate separate from y. Only an admitted attempt copies that candidate into the accepted vector. controller.apply then commits or rejects the clock and stores the next proposed step size. Rejected attempts do not move accepted time or change y.
The rate switch is known in advance, so the driver clips a step to land on it. The rate stays fixed throughout that attempted interval, including its endpoint stage. The next interval uses the other rate. This avoids asking a smooth Runge–Kutta error model to infer an unresolved discontinuity.
Make a controlled comparison
openpfc run model.cpp input.json --source-root .. --launcher=noneBoth controllers run at the same tolerance, used for both atol and rtol. The other inputs set the first proposed dt, transition and final end. Each rank writes results/attempts_rank_<rank>.csv, including rejections, accepted time, the accepted value, the exact solution and RHS counts. The scalar state is replicated across ranks here; the reduction illustrates collective decisions without claiming distributed ODE storage.
Tighten the tolerance from \(10^{-4}\) to \(10^{-7}\). Compare final error and work; do not infer that one controller is always better from one trajectory. Short intervals can also result from clipping to the rate switch or final time. A small last step is not, by itself, evidence of a new physical transient. The checks require exact rejection rollback, event alignment, consistent rank decisions, a more accurate tighter run, and four RHS calls per attempt. The current stepper does not reuse an FSAL stage.
The complete C++ model and starting input are available with this lesson. It also uses the shared lesson.hpp; keep its relative location when saving the files.