Skarstrom 0.6.0

Added — regression test for CSS acceleration on a thermal cycle

• A slow-tier test checks that Aitken acceleration shortens the approach to
  cyclic steady state on the thermal four-step LPP cycle (87 cycles without
  it, 39 with it) and that both runs stop at the same KPIs. The jump limit
  that replaced the |mu| <= 5 cap already gives the reduction; earlier,
  both runs took the same number of cycles.

Fixed — library re-saves under ctest --parallel

• The tests that re-save data/adsorbents/*.json now write through the
  suite's atomic-write helper (test::SaveAdsorbentToJsonAtomically): the
  rename over the target is retried while a reader holds the file open
  (Windows), and a failed write throws instead of leaving the old file in
  place unnoticed. The library files are unchanged.

Fixed — flaky tests reading example decks under ctest --parallel

• The tests that re-save examples/decks/*.json now write each deck to a
  temporary file in the same directory and rename it over the deck, so a test
  loading the deck at the same time always reads a complete file instead of
  truncated JSON (rapidjson internal assertion failure: IsObject()). On
  Windows, where the replace fails while a reader holds the file open, the
  rename is retried for up to about 5 s. The deck contents are unchanged.

Fixed — class versions visible wherever a deck is saved or loaded

• Every serialised engine class now declares its CEREAL_CLASS_VERSION in the
  header that holds its inline serialize template, not in its .cpp. A
  translation unit that saved or loaded a deck without seeing the declaration
  used version 0, so which version a save wrote (and which branch a load took)
  depended on which instantiation the linker kept. Saved decks do not change.
  A new unit test (CerealClassVersions.*) pins the version each class writes.

Fixed — class versions of the fluid, connections and CSS acceleration

• adsim::Fluid, adsim::Connection and adsim::AccelerateCyclicSteadyState
  now declare their CEREAL_CLASS_VERSION in their headers too, so every
  serialised engine class does. Saved decks do not change; the
  CerealClassVersions.* test now pins these three classes.

Changed — changelog entries are fragments in changelog.d/

• A pull request adds its changelog entry as a new file
  changelog.d/<topic>.md instead of editing CHANGELOG.md, so pull
  requests in flight together no longer conflict on the top of Unreleased
  (GitHub ignores the union merge driver). scripts/release.py reads the
  fragments as part of Unreleased; the release pull request moves them into
  CHANGELOG.md and deletes them. New commands: release.py fragments
  (check) and release.py collect. CONTRIBUTING.md, RELEASING.md and the pull
  request template describe it.

Added — cycle state: restart, warm start and one-cycle runs

• New public header adsim/solver/cycle_state.h: CycleState, the full state
  at a cycle boundary. It holds the DAE unknowns and their derivatives in the
  solver's order, a field map, and the differential entries as a flat vector
  (Differential / SetDifferential). It also holds every endpoint's logged
  history, which connected inlets replay, and the buffer-tank summaries. It
  records the deck structure it belongs to: components, layers with cell
  counts, steps, endpoints, step times and the number of cycles behind it.
• Simulation::FinalState() is the state at the end of the last completed
  run. RunCycle(state) runs exactly one cycle from a state, so FinalState()
  then gives the cycle map applied to it; RunCycle() with no state steps on
  from the last final state. SetInitialState(state) sets a warm start that
  later runs start from instead of the initial conditions. CheckState(state)
  gives the reason a state cannot start the deck: other components, layers or
  cell counts, another step structure, a changed state layout, or changed step
  times when connections or buffer tanks read the histories. A warm-started
  run counts the state's cycles towards the CSS minimum number of cycles.
• Python: skarstrom.CycleState; Simulation.final_state(), run_cycle(),
  set_initial_state(), warm_start(), clear_initial_state(),
  initial_state and check_state().
• The state is kept for the session only and is never written to the deck
  file. Runs without a state are unchanged.

Added — moving a cycle state to another grid

• RemapCycleState (adsim/solver/cycle_state.h) and
  Simulation::RemapState transfer a CycleState onto a deck with other
  cells per layer. Inside each layer, every per-cell field is remapped by
  overlap-weighted cell averages. The gas is remapped as component
  concentrations, so every component's gas and adsorbed inventory in every
  layer is kept to rounding. Connection histories and buffer-tank summaries
  are carried over unchanged. Refining by a whole-number ratio and coarsening
  back gives the original state.
• A warm start or RunCycle from a state on another grid now remaps it,
  where it was refused before. Different components, layers, void regions or
  step structure are still refused with the reason.

Added — electrical energy KPI

• The specific energy KPI is now also reported under its explicit name,
  electrical_energy (compressor and vacuum-pump work over the cycle per kg
  of net product, J/kg): a KPIRecord field and Python attribute, a key of
  KPIRecord.as_dict() / kpi_table() / the CLI's --json KPIs, a
  sensitivity output (Component.electrical_energy), an appended
  electrical energy [J/kg] column of the CLI's kpis.csv, and the
  ELECTRICAL_ENERGY enumerator of KPIType and KPIVariableType (the
  same enumerator as ENERGY_CONSUMPTION). The KPI's display label is now
  "Electrical Energy", and skarstrom run prints it as
  electrical energy [J/kg].

Changed — energy_consumption is documented as the electrical energy

• energy_consumption keeps its name, value and meaning in KPI records,
  sensitivity outputs (still the default), optimiser objectives, the cost
  estimate and the existing kpis.csv columns; it is documented as the
  specific electrical energy. Step energy sources (heaters, coolers) are
  still not included. Results do not change.

Fixed — extended Langmuir–Freundlich trace loadings

• extended_langmuir_freundlich no longer multiplies each competitive term by
  the smooth step y² / (y² + 1e-10) (with a softplus on y). That step biased
  trace loadings by about 1e-10 / y²: 1e-4 at y = 1e-3, 1 % at 1e-4, about
  50 % at 1e-5. The power (P y / P_ref)^n is now continued below y = 1e-5 by
  the same C1 quadratic the competitive Sips kinds use, so the isotherm is
  exact for y ≥ 1e-5 in both the solver and the single-component evaluation,
  with a finite Henry slope below and exactly zero loading at y = 0. Results
  of runs with trace components on this kind move; no bundled deck uses it.

Fixed — energy-balance docs equation

• The LTE energy-balance equation in the docs now matches the engine: the adsorption heat term no longer carries a stray $RT/P$ factor, its sign follows the stored heat of adsorption (positive for exothermic adsorption, so adsorption heats the bed), and the equation now shows the adsorbed-phase heat capacity in the accumulation term and the step energy source.

Added — GAB and BET water isotherms

• Two isotherm kinds for water and other condensable vapours, in the relative
  pressure x = p_i / p_sat(T): gab (Guggenheim-Anderson-de Boer with
  C = C0 exp(dH_C/RT), K = K0 exp(dH_K/RT), q_m = q_m0 exp(dH_m/RT); BET is
  K0 = 1, dH_K = 0) and gab_young (the layer-energy form of Young et al.
  2021 and Low et al. 2023, with the papers' own parameters). The saturation
  pressure comes from CoolProp's saturation ancillary (water: within 0.02 % of
  IAPWS-95). The loading is exact for 0 <= x <= min(1, 0.95/K) and continued
  smoothly beyond (the Henry line for x < 0, the tangent above), so it stays
  finite at and above saturation. Not competitive: water adsorbs on its own.
• adsim::IsostericHeat(isotherm, component, T, P, y): the isosteric heat
  implied by any isotherm's temperature dependence, to help choose the
  constant heat of adsorption the energy balance uses.
• The GAB kinds work under the IAST mixing rule (closed-form spreading
  pressure).

Fixed — IAST with an isotherm kind it does not know

• An isotherm kind without an IAST form now raises an error instead of being
  treated as inert (which silently gave that component zero loading).

Added — grid convergence analysis in the engine

• New public header adsim/solver/grid_convergence.h: observed order of
  accuracy, Richardson extrapolation, Roache's grid convergence index and the
  cell count for a target error, for a result (a cycle KPI, say) on three or
  more grids at any constant refinement ratio. It reports whether the grids
  are in the asymptotic range: monotone, oscillatory (no order reported),
  divergent, grid independent or too few levels. The cycle-KPI grid
  convergence test now uses it (same numbers).

Added — dense finite-difference check of the solver Jacobian

• New integration test JacobianDenseCheck: at states in the middle of each
  step of real decks, the Jacobian IDA receives (coloured, and the band
  evaluations that learn and audit its pattern) is compared entry by entry
  with a dense forward difference over the whole matrix, with no colouring
  and no band assumption. It covers upwind, WENO and van Leer, thermal and
  isothermal runs, Kozeny-Carman and Ergun momentum, a wall energy balance,
  a breakthrough, a PSA and a TSA cycle, and each deck's grid as well as 10
  and 20 cells. Every entry agrees to rounding (relative error about 2e-16);
  a guard case shows that the comparison detects an incomplete pattern.
• BandedJacobian::SetEvaluationObserverForTesting: a read-only test hook
  that sees every assemb...

Skarstrom 0.4.0

Added — automatic releases

• Merging a version change into main now releases it: release.yml runs on
  a push to main that changes the version in CMakeLists.txt, creates the
  tag v<version> itself, uses the version's CHANGELOG section as the
  GitHub Release notes and publishes to PyPI (TestPyPI for pre-releases). Before,
  only a hand-pushed tag released, and 0.2.0 and 0.3.0 were never tagged or
  published. A pushed tag still works. A final version is published only if
  CHANGELOG.md has its section with nothing left under Unreleased, and
  CITATION.cff carries it.
• New workflow prepare-release.yml keeps one "Release X.Y.Z" pull request
  open (branch bot/release) whenever Unreleased has entries, rebuilt from
  main on every changelog change. The next version comes from the Unreleased
  subsections: Added, Changed, Removed or Deprecated give a minor release,
  Fixed or Security alone a patch release; major releases and pre-releases are
  chosen by hand (Run workflow with bump / suffix).
• scripts/release.py (current, next, prepare, check, notes) does the
  version arithmetic and the three-file release change, with tests in
  scripts/test_release.py (run by CI's docs job). The version on main is
  now the last released version; it is no longer bumped for development.
  RELEASING.md and CONTRIBUTING.md describe the new flow.

Added — KPIs of every cycle (KPIHistory)

• The solvers record the KPIs of every tracked component when each cycle
  completes (adsim::KPIHistory, adsim/solver/kpi_history.h;
  CoupledSolver::GetKPIHistory, CoupledSolverChuncked::GetKPIHistory), so a
  run's KPIs cover all its cycles, not only the last
  settings.max_cycles_to_keep the logs hold. Simulation::KPIs() (Python
  sim.kpis(), the CLI's kpis and kpis.csv) now lists every cycle:
  cycles still in the log are computed from it as before, evicted ones come
  from the history. New: Simulation::KPICycles(),
  Simulation::KPIHistoryRecords(), and adsim::ComputeCycleKPIs, the
  per-cycle computation both use. KPIRecord moved from simulation.h to
  kpi_history.h (which simulation.h includes).

Added — own ramp time for pressure ramps; exponential factor

• Linear, exponential and log-linear pressure profiles (pressure boundaries
  and valve streams) take an optional ramp time of their own:
  use_step_time_ (deck field Use_Step_Time, default on) keeps the ramp
  over the whole step as before; off, the profile's ramp_time_ /
  step_time_ is used and must not exceed the step time
  (CheckPressureProfile, run by Cycle::PerformChecks). Registry:
  endpoint[i].pressure_<kind>_use_step_time, and pressure_<kind>_time when
  it is off. Python: sk.linear(end, time=10) for a pressure.
• The exponential profile's factor is a parameter: factor_ (deck field
  Factor, registry <head>_exponential_factor, Python
  sk.exponential(..., factor=)), v(t) = v1 + (v0 - v1) F^(-t / t_r).
  The default 100 is the old fixed value (99 % at the ramp time); F must be
  greater than 1 (CheckProfileFactor).
• Class versions: LinearBoundaryValue 4, ExponentialBoundaryValue 4,
  LogLinearBoundaryValue 2. Older decks load with the step time and F = 100,
  so their results are unchanged.

Added — method-of-manufactured-solutions hook and order test

• Verification hook UserData::manufactured_sources_
  (adsim/equations/manufactured_sources.h): source terms, functions of z
  and t, for the overall balance, each component balance, each loading's
  rate law, the momentum faces (a pressure-gradient source on Darcy faces
  and interior Ergun faces) and the LTE energy balance, added to the
  residuals in the scaling of the neighbouring physical terms. Testing only:
  not serialised and not reachable from decks, the parameter registry,
  Python or the web. With no source set the results are bit-identical
  (bench_solver logs of two example decks are byte-identical to the
  previous residual code) and the residual cost is unchanged.
  UpdateResiduals records the integrator time in UserData::residual_time_
  for it.
• core/test/validation/test_manufactured_solutions.cpp: a smooth,
  time-dependent manufactured solution with dispersion, LDF kinetics on a
  Langmuir isotherm, Kozeny-Carman flow and (thermal case) the energy
  balance, sources derived by hand and checked by finite differences.
  Observed orders on 20-160 cells: upwind 1.0 in every field; WENO and
  van Leer 2 (WENO interior 2.5-2.8 before the second-order central terms
  take over); with the energy balance the temperature is first order for
  every scheme (the convective face temperature is always upwind). About 4 s.

Fixed — the CSS state test never passed for an inert bed

• For a bed of inert packing (an adsorbent with no isotherm) the adsorbed
  loadings are round-off, 1e-26 to 1e-22 mol/kg of either sign. The state
  test of the mass-balance CSS methods
  (CyclicSteadyState::EndOfCycleStateChange) scales the loading change by
  the RMS total loading and skipped the field only when that scale was below
  1e-300, so it divided noise by noise and reported a change of 0.77 to 51
  on a bed whose state repeats: a mass-balance method with the state test on
  (the default) never stopped.
  The loading scale is now floored at the solver's absolute tolerance on the
  loading unknowns (1e-8 mol/kg); a loading below it is not resolved. Every
  adsorbing deck has an RMS loading far above the floor, so its state change
  is unchanged to the last bit. Test:
  CyclicSteadyState.InertBedReachesCSSUnderTheStateTest (the shared PSA
  cycle with inert packing stops at the minimum-cycle floor, cycle 4). Found
  by the dynamic-tank CSS test of the verification programme (WP19-E2a).

Added — convergence verification of cycle results (grid, integrator tolerance, CSS)

• CycleKpiGridConvergence (core/test/validation/test_grid_convergence.cpp):
  purity, recovery and specific energy of the shared PSA cycle at cyclic
  steady state on 10, 20, 40 and 80 cells, for each advection scheme.
  Convergence is monotone; the observed order (Richardson, three finest grids)
  is 1.55 for upwind (inside the [1, 2] its error expansion allows) and 2.4
  to 2.5 for WENO and van Leer; Roache's GCI of the 80-cell upwind KPIs is
  0.06 to 0.12 % and contains the higher-order solutions.
• IntegratorTolerance (test_integrator_tolerance.cpp): the error falls with
  reltol at slope 0.68 (breakthrough) and 0.77 (cycle state), inside the
  [0.5, 1] of BDF error-per-step control, and floors at the absolute
  tolerance; max_step changes results by less than the tolerance error. The
  KPIs converge with reltol at slope 0.42, because the boundary flows are
  integrated over the logged points and the flow transient after a pressure
  switch (~2e-5 s) is resolved only by the integrator's own steps.
• CSSConvergence (test_css_convergence.cpp): the approach to CSS is
  geometric (rate 0.106 isothermal, 0.966 with the energy balance), the KPIs
  follow it, the error the stopping rule leaves is g Delta rho / (1 - rho)
  (measured 6.11e-5 in recovery at the default state tolerance, 6.21e-5
  predicted), and Aitken acceleration stops 25 cycles earlier on the same KPIs.
• scripts/verification/grid_convergence.py writes the grid study of the
  example decks (observed order, GCI, extrapolated KPIs); the table and the grid
  each studied deck needs are on the Numerics page ("Convergence of the
  results").

Added — extended Langmuir-Freundlich isotherm (toPSAil form)

• New competitive isotherm kind extended_langmuir_freundlich
  (IsothermType::EXTENDED_LANGMUIR_FREUNDLICH,
  ExtendedLangmuirFreundlichIsotherm): q_i = (k1 + k2 T) z_i / (1 + sum z_j),
  z_i = k3 exp(k4 / T) (P y_i / P_ref)^{n_i}, 1/n_i = k5 + k6 / T (Yang 1987,
  toPSAil's calcIsothermExtLangFreu.m; P_ref = 1e5 Pa keeps toPSAil's
  bar-based k3). No existing kind maps it unless k2 = k6 = 0. Registry,
  Python isotherm_kinds(), IAST (an exact Sips term) and the docs know it;
  the fitter does not offer it yet. Near y = 0 each term is regularised with a
  step at 1e-5 (loading changed by 1e-10 / y^2: 1e-4 at y = 1e-3).

Fixed — breakthrough decks no longer stall on fine grids

• Ward2022_Q50ccm ran on 10-160 cells but failed on 320 with every
  advection scheme ("too much work": 50 000 steps, 24 000 Newton convergence
  failures, 72 000 Jacobians) as the CO2 front left the column. The
  difference-quotient Jacobian perturbed the pressure by IDA's increment
  floor reltol |P| + abstol (0.1 Pa at 1 bar). A pressure outlet imposes the
  outlet-face velocity as the explicit Ergun root of the end-cell pressure
  drop, which is linear only over A² / (4B), a pressure scale proportional to
  the cell width (0.12 Pa on 60 cells, 0.02 Pa on 320): the 0.1 Pa secant
  was half the true slope there, and the modified Newton iteration diverged
  unless the step was tiny. The pressure columns of the BandDQ, Coloured
  and ColouredCheck Jacobians now use sqrt(eps) max(|P|, |h P'|) floored at
  the pressure's absolute tolerance (1.5e-3 Pa at 1 bar); every other column
  keeps IDA's increment. The deck now runs on 320 cells in 2 169 steps with
  24 convergence failures (3 s), and on 640 cells in 2 047 steps. Results of
  every deck move at integrator-tolerance level, since the Jacobian's
  pressure columns change in the last digits: across the paper validation
  metrics, front times by at most 2e-5 relative and curve RMSEs by at most
  1e-5 absolute, and every paper test still passes.
  JacobianMode::IdaInternal (IDA's own Jacobian) is unchanged. Test:
  FineGridBreakthrough.Ward2022Q50ccmConvergesTo320Cells.

Fixed — WENO / van Leer outlet-cell mole-fraction sum drifts during fast outflow

• With AdvectionScheme::WENO or VanLeer the cell at a bed end that gas
  leaves through could see its mole fractions drift from summing to one
  while...

Skarstrom 0.1.0

What's Changed

• Fix first CI run: map-order-dependent isotherm sums (macOS), workflow lint by @samaffleck in #1
• Build CoolProp for the native architecture (fix Linux aarch64 wheels) by @samaffleck in #2
• Python object API, builder, LLNL/pint units, KPI recycle role, LPP runs to CSS by @samaffleck in #3
• Validation: refresh the Perez 2019 and Nguyen 2023 findings for the CSS runs by @samaffleck in #4
• Release 0.1.0: wheels for CPython 3.10–3.15 by @samaffleck in #6
• Seed the adsorbed-phase heat capacity in AddComponentToDeck by @samaffleck in #7
• Test that the edit cascades leave nothing for Initialise to seed by @samaffleck in #8

New Contributors

• @samaffleck made their first contribution in #1

Full Changelog: https://github.com/samaffleck/skarstrom/commits/v0.1.0