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 the outflow lasted (a valve surge, a blowdown): 5.0e-4 with WENO and
  4.5e-4 with van Leer in ErgunValveSurgeTest (Kozeny-Carman), against
  1.1e-5 with upwind. The schemes clamped every face mole fraction to [0, 1], the bed-end faces
  included; those take the upwind cell's value and are not rescaled to its
  sum, so an end cell holding pure nitrogen at y = 1 + 1e-9 sent out a face
  that summed to one instead of the cell's sum, and the cell's departure from
  one was no longer carried out of the bed. Only the reconstructed values are
  clamped now; the upwind faces are left as they are, as under Upwind. The
  same clamp made the logged balances disagree with the solver (the stream
  logs integrate the end cell's value, the solver used the clamped face):
  the DAC charge with WENO had its total mole balance off by 3.45e-2 mol
  (now 2.7e-5, the upwind value) and the outlet cell's sum by 4.0e-5 (now
  9e-15); the temperature purge with van Leer logged 2.4e-5 mol of a
  nitrogen that is in neither the feed nor the bed (now 1.3e-8). The drift
  also reached 1.6e-6 (WENO) in the steep-front breakthrough of
  AdvectionSchemeConservation at reltol 1e-8; every scheme now keeps about
  1e-12 there. Upwind is unchanged (its faces were never clamped). Tests:
  AdvectionSchemeConservation.EveryFaceSumsToItsUpwindCellSum (the
  discrete identity, face by face); the valve-surge test run with every
  scheme (ErgunValveSurgeTest.*WENO, *VanLeer), which bounds the outlet
  cell's sum at the end of each feed step by reltol; UserDataTest.DAC and
  TemperaturePurge run with every scheme (*WENO, *VanLeer), with their
  logged balances bounded at the upwind level; and the steep-front sum
  check, now at reltol 1e-8 with a 1e-9 bound. Found by the WENO-everywhere
  experiment of the verification programme (WP20-E1).

Fixed — the ColouredCheck Jacobian diagnostic no longer changes the run

• In JacobianMode::ColouredCheck the band evaluation built for the
  comparison learned every coupling it found into the pattern without
  recolouring, and its count was discarded: a coupling the coloured pattern
  missed stayed missing from the coloured Jacobian for the rest of the run,
  no audit ever counted it (NumMisses() stayed 0), and the check run took a
  different path from the Coloured run it was meant to observe. The
  comparison now only observes (no learning, and no residual-call,
  out-of-band or pattern counts, and a failed comparison evaluation does not
  fail the step), so a ColouredCheck run is the Coloured run exactly and
  reports a missed coupling until an audit adds it. Test:
  ColouredJacobian.ColouredCheckObservesTheColouredRun (it starts from a
  pattern learned from one evaluation, the new testing switch
  IDASettings::jacobian_learn_once_for_testing_).
• With the WENO and van Leer schemes the pattern learned in the first cycle
  missed couplings of the columns in the three cells at each end (e.g. the
  first cell's mole fractions on the inlet ghost cell's other mole
  fractions); the audits added them later. Those columns now also take the
  interior pattern, so the coloured Jacobian equals the band DQ Jacobian entry
  for entry for all three schemes (ColouredJacobian.MatchesBandDQ_*_WENO /
  _VanLeer). Results of WENO / van Leer runs in the default Coloured mode
  move at integrator-tolerance level (they are now bit-identical to BandDQ
  runs); upwind runs are unchanged, bit for bit.

Added — constant gas heat capacity and viscosity; heat of adsorption as −ΔU + RT

• fluid.gas_heat_capacity_mode: coolprop (default, unchanged) or
  constant, with a molar value per component in
  fluid.gas_heat_capacity[<component>] [J/(mol·K)] (a new component starts
  at 29.1, 7R/2). The constants replace the CoolProp values that
  Fluid::Initialise caches; the molar mixing rule is unchanged, so equal
  values give that heat capacity at every composition. The Gas phase
  adsorbed-phase heat capacity follows them.
• fluid.viscosity_mode: polynomial (default, unchanged) or constant,
  with one mixture viscosity fluid.viscosity [Pa·s] in every cell.
• fluid.heat_of_adsorption_convention = "isosteric_rt": the stored heat
  of adsorption is −ΔU, the internal energy of adsorption, and the heat
  released per mole adsorbed is −ΔU + R·T at the local adsorbed-phase
  temperature (LTE bed or LTNE solid temperature). This is the source term
  of Yancy-Caballero et al. (2020) / PEESEgroup PSA, −(ΔU − RT); a constant
  isosteric value fixes it at one temperature.
• With the three switches a deck reproduces models that use fixed gas
  properties (C_p,g = 30.7 J/(mol·K), μ = 1.72e-5 Pa·s in Yancy-Caballero
  2020 and Haghpanah 2013). Defaults are unchanged, so every example deck
  and validation result is unchanged. Fluid is class version 4; older
  decks load with the defaults. Non-positive constants fail the deck
  checks. Registry, Python (sim.fluid.gas_heat_capacity_mode,
  gas_heat_capacity[...], viscosity_mode, viscosity) and the script
  export carry the new fields. Docs: equation-of-state.md (constant gas
  properties), energy-balance.md (heat-of-adsorption conventions). Tests:
  FluidPropertyModes* and the IsostericRT cases of
  ComponentEnergySourceLog.

Added — valves can read the bed at the end cell

• Solver setting settings.valve_bed_pressure (IDASettings::valve_bed_pressure_,
  "Valve_Bed_Pressure"; IDASettings cereal version 7): face (the default,
  unchanged behaviour) evaluates the bed side of every valve stream's linear
  law at the bed-end face, end_cell at the last cell, as toPSAil does. A
  connected valve then reads its partner's end-cell pressure, which every
  stream and pressure boundary now logs (StreamData::GetEndCellPressure).
  Decks from before load with face.

Added — sensitivity studies saved with the deck

• A deck now saves its sensitivity study's set-up (UserData::sensitivity_study,
  Sensitivity_Study in the file; UserData cereal version 4): inputs and their
  distributions, outputs, run policy, seed, and the screening and propagation
  sizes. Results are not saved. Decks from before load with no study.
• adsim/solver/sensitivity_study.h holds Distribution, Input, RunPolicy
  (moved from sensitivity.h, which includes it; the API is unchanged) and the
  new Study. sensitivity.h adds StudyOutputs and CheckStudy(deck, study).
• Python: Study.from_deck(sim, **overrides), Study.save_to(sim) and
  Study.setup(); the study's sizes are constructor arguments, the defaults of
  screen(), propagate() and Screening.important().
• CLI: skarstrom sensitivity deck.json without --input runs the deck's
  saved study; options given override it.

Added — cross-code regression against the PEESE model and toPSAil

• core/test/crosscode/: twelve single-column breakthrough and purge cases
  (CO2/N2 on 13X with concentration- and pressure-form Langmuir and Henry
  isotherms, 1-50 % feeds, purge, adiabatic, a 25 % Ergun pressure drop,
  CO2/CH4 on LiX, a steel wall with ambient loss) run in the PEESEgroup/PSA
  column model (MIT) and in toPSAil (GPL-3.0; numerical outputs only), frozen
  as data with their decks and replayed by CrossCodePeese.* and
  CrossCodeToPSAil.* (about 15 s). Bounds come from the measured, attributed
  differences: isothermal cases agree within 0.6 s on t50 at the same grid;
  the adiabatic offsets (up to 14.5 s) are the CoolProp gas heat capacity and
  the constant heat of adsorption. The folder README has the attribution.

Fixed — a connection replays a bottom pressure boundary's outflow with the right sign

• PressureBoundaryCondition::Log did not record the bed end it logged at,
  so a pressure boundary's StreamData always said TOP. For a pressure
  boundary at the BOTTOM of its bed, StreamData::GetOutflow then returned
  gas leaving the bed as a negative outflow, and a no-valve stream connected
  to it (e.g. a rinse or heavy-reflux inlet fed from a blowdown or evacuation)
  drew that gas out of its own bed instead of taking it in: the connection
  destroyed twice the transferred moles, and with a step-change blowdown IDA
  failed at the start of the receiving step. The log now records the
  boundary's location, as a stream's log does. A TOP pressure boundary was not
  affected, and no bundled deck, cycle template or paper case connects a
  stream to a BOTTOM pressure boundary, so no published result moves. Tests:
  ConnectedStreamTest.PressureBoundaryDonorLogRecordsItsBedEnd,
  ConnectionSign.*.

Added — verification: extreme operating corners

• core/test/validation/test_extreme_conditions.cpp: very high and very low
  pressure ratios, near-zero flow, very short and very long steps and high
  temperatures each converge or fail cleanly, with finite and physical states.

Added — verification: thermodynamic bound on the separation energy

• core/test/validation/test_thermodynamic_bound.cpp: at cyclic steady state
  the logged pump energy of a vacuum-swing and a pressure-swing cycle exceeds
  the reversible pump work, the bed's exergy change and the minimum work of
  separation of the products, in that order.

Added — verification: parameter limits and degenerate mixtures

• core/test/validation/test_parameter_limits.cpp: the rate at which the
  model reaches its limits (LDF k to infinity and to zero, heat of adsorption
  to zero, wall heat capacity and bed-wall coefficient to infinity, Langmuir b
  to zero, pressure ratio to one, trace feed), and two degenerate mixtures
  (identical components keep their feed ratio; zero affinity is inert).

Added — verification: metamorphic relations

• core/test/validation/test_metamorphic_relations.cpp: equivalent decks must
  give the same solution to within the integrator's measured error, for the
  component order, a packing split into identical layers, a step cut in two,
  inputs written in other units, and the step a cyclic steady state is started
  from (with a connection whose lag changes under the rotation).
  BreakthroughCase gains n_layers.

Added — verification: evacuation, pressure equalisation and buffer tank

• core/test/validation/test_unit_operations.cpp checks constant-speed
  evacuation (P = P0 exp(-S t / V a R T)), two-bed pressure equalisation
  through connected valve steps at the cyclic steady state (mole balance and
  the exp(-2 cv t / 1e5 V a) decay of the pressure difference; with a
  Langmuir isotherm, conservation of the two beds' inventory), and a
  buffer-tank fill (it delivers the set percentage of the vented moles), each
  without and with Henry adsorption.

Added — verification: thermal fronts

• core/test/validation/test_thermal_fronts.cpp checks the speed of an
  inert-bed thermal front against the energy-balance jump condition
  (including the gas heat-capacity term), and the steady exponential axial
  temperature profile of a flowing bed losing heat through the wall, with
  first-order convergence to it.

Added — verification: blowdown composition path

• core/test/validation/test_blowdown_composition_path.cpp checks that a
  binary on linear isotherms follows the equilibrium-theory (Knaebel-Hill)
  composition path y(P) during a valve blowdown, with the departure from
  the path equal to the LDF lag (it halves when k doubles).

Added — verification: binary competitive Langmuir fronts

• core/test/validation/test_competitive_langmuir_fronts.cpp compares a
  binary extended-Langmuir breakthrough (constant separation factor) with
  equilibrium theory: the roll-up plateau of the weaker component from the
  coherence condition, and the times of both fronts from the jump conditions.

Added — verification: Langmuir desorption simple wave

• core/test/validation/test_desorption_simple_wave.cpp checks that a
  Langmuir desorption front spreads as the equilibrium-theory simple wave:
  the outlet times of three levels inside the wave converge at first order to
  the method-of-characteristics times, and their Richardson extrapolation
  agrees to 0.1 %.

Added — verification: LDF-limited front against the Thomas solution

• core/test/validation/test_ldf_thomas.cpp compares the whole outlet curve of
  a Henry-isotherm, LDF-limited, dispersion-free front with the
  Anzelius/Thomas closed form, and checks that the grid error is exactly the
  upwind scheme's: the variance minus L dx R^2/u^2 is the Thomas variance,
  the curve error falls at first order and, after the closed-form O(dx)
  correction, at second order.

Changed — CI skips the builds for docs-only pull requests

• A new changes job lists a pull request's files. The C++ (Linux, macOS,
  Windows), Python wheel, Emscripten and subproject jobs run only if one of
  them is code: anything except Markdown, docs/, mkdocs.yml,
  requirements-docs.txt, .gitattributes, CITATION.cff and issue
  templates. docs/python.md (its example is run by python/tests),
  python/README.md, THIRD_PARTY_NOTICES.md, LICENSE and NOTICE count as
  code. The docs and workflow-lint jobs always run, and pushes to main,
  nightly and manual runs build everything.
• .gitattributes: CHANGELOG.md uses git's union merge driver, so two
  pull requests adding entries at the top of Unreleased no longer conflict in
  local merges.

Fixed — connected streams conserve moles when the donor's flow reverses

• A connected no-valve stream (e.g. a velocity-0 equalisation consumer)
  replayed the magnitude of the connected stream's logged flow. When the
  donor's flow reversed (a valve outlet drawing gas back into its bed), the
  backflow reached the receiving step as more inflow and the connection
  created gas. The replay now keeps the sign: what left the donor's bed enters
  the receiver and what flowed back into it is drawn out (StreamData::GetOutflow;
  the log records the bed end it was taken at). Test:
  ConnectedStreamTest.ReplayKeepsTheSignOfTheDonorFlow.
• Vysyaraju 2024 PSA_Exp1 was the one case affected: PER took in 29 % more
  gas than PED released, so the conservation-based He recovery read 17 % (product
  He / fed He was 69 %). With the fix the two recovery measures agree (63.1 %),
  He purity is 99.36 %, and the tolerances are re-measured (see the paper's
  README); the deck's step cap is 0.5 s. The case still fails on macOS and
  Windows (it did before the fix too): its single-bed equalisation lags PER's
  pressure by one cycle, and which state the cycle settles in depends on
  round-off. Linux passes.
• CI: one run per pull request (pushes to feature branches no longer start a
  second run), superseded runs are cancelled, the paper replications run on
  pull requests (full suite on Linux and macOS, fast subset on Windows), Windows
  builds with Ninja + ccache (windows-ninja-release preset), the Emscripten
  builds use ccache, and ctest's per-test timings are cached so the longest
  tests start first.

Fixed — face values of logged cell variables on non-uniform meshes

• EquationLog::GetFaceVariable / GetComponentFaceVariable interpolate
  with the right weights. The cell-centre-to-face interpolation weighted
  each neighbour by its own width, (y_L dx_L + y_R dx_R) / (dx_L + dx_R),
  which biases the face value towards the larger cell. It is now
  (y_L dx_R + y_R dx_L) / (dx_L + dx_R), exact for a linear profile, so a
  face between cells of different widths (between layers with different cell
  sizes, or between a void cell and the bed) sits where it should. The two
  boundary faces now take the plain mean of the ghost and end cells. That is
  the face value the boundary conditions impose (P_face = (P_0 + P_1) / 2),
  and it matters because the ghost cell copies the width of the first or last
  bed cell, not of a void end cell. A uniform mesh is unchanged: every
  example deck, the paper validation and docs/validation.md are unaffected.
  Tests: EquationLogFaces.* (hand-built log, non-uniform mesh, linear
  profile exact at every interior face) and
  VoidRegionSolverTest.FaceReconstructedBoundaryMolesMatchTheStreamLogs
  (with voids, the face-rebuilt boundary moles now equal the stream logs to
  round-off; they were 1e-5 to 2e-5 off).
  What moves: the feed/product KPIs (purity, recovery, productivity,
  energy) do not use these face values and do not change. The total
  mass-balance CSS check (css.method 2), which reads the boundary-face
  pressure and temperature, changes for beds with void regions. Its error had
  a spurious floor. On Perez 2019 case A with both voids enabled (1e-4 m³,
  void cell 0.124 m against 0.01 m bed cells), the floor was about 3.46e-3
  (tolerance 5e-3) and is now 2.4e-5, the same order as without voids
  (1.9e-5); CSS is still reached at cycle 75. With 3e-4 m³ voids the old
  floor of 1.26e-2 never passed the tolerance (no CSS in 200 cycles); CSS is
  now reached at cycle 198. Face profiles plotted by front ends also move at
  layer and void interfaces.

Fixed — every cycle step starts the integrator afresh

• A step can no longer step over a short pulse at its start. Cycle steps
  ran back to back on one IDA history, so the first internal step of a step
  inherited the size of the previous step's last, tstop-clipped, step. The
  boundary condition's algebraic jump is out of the error test
  (suppress_alg_) and an implicit step samples the forcing only at its end,
  so a step whose forcing is a pulse shorter than that inherited step was
  integrated as if the pulse were not there. On Vysyaraju2024_PSA_exp1 the
  PER step replays the PED blowdown (2.8 mol/s decaying in about 0.1 s): when
  BLO happened to end on a 0.2 s remainder instead of a 0.015 s one, PER's
  first step jumped the whole pulse, the bed missed its equalisation gas and
  the run restarted its approach to cyclic steady state (PRES uptake
  0.1497 → 0.284 → −0.002 → 0.071 mol …). Which cycle tripped hung on
  round-off: natively at cycle 3 (and at cycle 48 with a reltol 1e-4 away),
  in the WebAssembly build at cycle 51, so the web solver disagreed with the
  native one at the last cycle (relative error 0.72). The engine before
  #36/#37 tripped too, but early (cycle 6 natively), and had re-converged by
  cycle 60. Now
  SolverData::RestartForNewStep, called by the blocking and the chunked
  solver at the start of every cycle step, restarts IDA from a consistent
  initial condition at the new boundary conditions with fresh history and
  hin = init_step (falling back to the unchanged end-of-step point if the
  consistent IC fails). Native and WebAssembly runs of the deck now agree to
  3e-5 per cycle, and native runs no longer depend on the reltol's last
  digits. The run's IDA statistics stay cumulative across re-initialisations
  (they used to restart from zero at each retry restart). Tests:
  StepBoundaryRestart.*. Costs about 6 % run time on the deck. Other
  decks' endpoint totals move by at most 3.4e-4 relative (Nguyen 2023 VII);
  Vysyaraju 2024's PRES uptake at cycle 60 by 1.2e-3.

Fixed — gas-mixture heat capacity mixed on the molar basis

• The gas heat capacity of a mixture is mass-consistent. UpdateHeatCapacity
  mixed the per-kilogram pure-component heat capacities that Fluid caches
  from CoolProp by mole fraction, C_p,g = Σ y_i C_p,i, which is right only
  when every component has the same molar mass. It now uses
  C_p,g = Σ y_i M_i C_p,i / Σ y_i M_i [J/(kg·K)] (the molar mixture heat
  capacity over the mixture molar mass), so ρ_g C_p,g = C Σ y_i M_i C_p,i.
  For 50/50 H₂/CO₂ the gas heat capacity was about 174 J/(mol·K) and is now
  about 33; for CO₂/N₂ 15/85 it was 1.3 % high. Every user of the gas heat
  capacity takes this one value: the LTE and LTNE energy balances, the
  convective enthalpy terms, the thermal dispersion and heat-transfer
  correlations, and the logged heat capacity. The rule lives in
  Fluid::MixtureHeatCapacity / Fluid::MixtureHeatCapacityFromSums, with
  Fluid::MolarHeatCapacity (C_p,i M_i) and
  Fluid::AdsorbedPhaseHeatCapacity (the Gas phase adsorbed-phase option,
  which was already per component and is unchanged) replacing the two copies
  in the equations. A single-component gas, or a mixture with all but one
  mole fraction zero, is unchanged. The stream "enthalpy flows" and the pump
  and compressor energy do not use the gas heat capacity (they are molar flow
  × temperature and polytropic relations) and are unchanged. Tests:
  GasHeatCapacityTest.* and GasHeatCapacitySolverTest.*.
  Results move for decks with light components: Casas 2012 (H₂/CO₂)
  breakthrough times fall by 0.2–2.2 % (Fig. 4 50/50: H₂ 87.98 → 87.40 s,
  CO₂ 445.1 → 436.0 s; Fig. 6 75/25: H₂ 133.2 → 131.5 s, CO₂ 359.3 →
  351.4 s), the Streb 2021 He-carrier breakthroughs move (thermocouple RMSE
  down 2–18 %, the A-series outlet and front-time errors up), Nguyen 2022
  front times fall by up to 1.5 %, Vysyaraju 2024 He recovery rises from
  17.54 to 17.91 %; CO₂/N₂ and O₂/N₂ cycles move in the fourth digit.

Fixed — bed-variable results at the right positions and boundary values

• Simulation::Series / sim.series and Simulation::Profile /
  sim.profile read the log at the right places. They mapped the raw cell
  rows, ghost cells included, onto the face positions, so every profile and
  every z query sat half a cell downstream, bottom returned the inlet ghost
  cell and top the last cell. In an adiabatic CO2/N2 breakthrough (feed
  313.15 K, 30 cells, t = 500 s) series("temperature", location="bottom")
  read 281.3 K, the ghost's 2 T_feed − T_1 (as low as 269 K during the run),
  and top pressure read 100004 Pa against the 100000 Pa set pressure. Now:
  • cell variables are at the cell centres and face variables (velocity,
    volumetric and molar flow rate, approach to fluidisation) at the faces;
    profiles leave out the ghost cells (AxialProfile::faces / prof.faces
    says which grid z is). Simulation::Positions / sim.positions are
    the cell centres, as documented; the new FacePositions /
    sim.face_positions are the faces.
  • top and bottom are the physical boundary values. Pressure,
    temperature and mole fraction are the logged values of the step's
    endpoint at that end, identical to stream_series (several streams are
    mixed by flow): at an inlet the value the boundary imposes (feed
    composition, feed temperature after any compressor, set pressure), at an
    outlet the bed-end face value. A closed end, and every other variable,
    gives the face value (the mean of the ghost and the end cell; a mole
    fraction's ghost, which is the face value itself; the nearest packed cell
    for adsorbent-only variables). A z between a bed end and the first cell
    centre interpolates between the two. The time base is unchanged.
  • The same logic is the new public header adsim/solver/bed_series.h
    (ExtractBedTimeSeries, CellCentrePositions, IsFaceRow);
    ExtractBulkTimeSeries / ExtractComponentTimeSeries
    (parameterisation.h) use it and take an optional BedBoundarySources
    (boundary log and steps) for the endpoint values. Front ends that fit to
    or plot these series move with it.
• stream_series(..., "pressure"): the bed-end face pressure of an
  endpoint (a pressure boundary's set value) is now readable
  (BoundaryVariableType::PRESSURE, appended to kBoundaryVariableTypeList),
  and BoundaryLog::GetStreamData returns an endpoint's logged data.
• Fitted or validated numbers that compared bed-end or z series (the
  example scripts' outlet curves and thermocouples) move slightly: outlet
  compositions and temperatures were already the end cell, so mainly
  thermocouple positions (half a cell) and inlet-end series change.

Fixed — component mass balance with axial dispersion (Danckwerts inlets)

• Inlet faces are Danckwerts. With non-zero axial dispersion the component
  balance did not close at cyclic steady state: the species residual let a
  dispersive flux eb·D·C·(y₁ − y_in)/Δz cross the inlet face on top of the
  convective u·C·y_in, because the inlet ghost cell holds the feed composition.
  Totals still balanced (the dispersive fluxes sum to zero over components) but
  each component did not, by an amount proportional to D and nearly independent
  of the grid, and the stream logs, KPIs and the CSS mass-balance check never
  saw the extra species (reported: CO₂ out/in 1.0035 at D = 3.1e-4 m²/s and
  1.032 at 3e-3 m²/s). No dispersive flux now crosses the two boundary faces of
  the bed (UpdateMoleFraction), so an inlet face carries exactly u·C·y_in and
  the first cell may differ from the feed; at outlets and closed ends the ghost
  mirrors the first cell, so nothing changes there. The energy balance had the
  same defect: conduction across the inlet face, and an upwind value at the
  inlet face equal to the mirrored ghost temperature 2·T_in − T₁ rather than
  T_in. Both are fixed (UpdateTemperature, UpdateFluidTemperature). The WENO
  and van Leer reconstructions take the ghost value at the boundary faces
  (they averaged ghost and first cell, which at an inlet is not the feed).
  GetTotalComponentMolesIn (CSS mass-balance check, mass-balance panels)
  reconstructs the boundary flux without the dispersive term, matching the
  residual.
• Flow boundaries on Darcy layers (Kozeny–Carman, Blake–Kozeny) deliver exactly the stream flow.
  FlowBoundaryCondition::Apply computed the boundary-face flow of a Darcy
  layer with the ghost cell's mobility alone, K/(μ_ghost·Δz), while the face
  velocity that carries the gas into the bed obeys the face's Darcy law with
  the two half-cells in series. Whenever the ghost and first-cell viscosities
  differ (feed composition vs the bed gas, or the mirrored ghost temperature)
  the bed took in a different flow from the one the stream logs record (≈1 %
  in the new test; several per cent in the Streb 2021 decks, whose beds start
  in helium). The condition now matches the stream flow against the face's own
  velocity unknown, as Ergun faces (and the momentum-balance docs) already did.
• New regression tests DispersionMassBalance.*
  (core/test/validation/test_dispersion_mass_balance.cpp): a dispersive
  (D = 3e-3 m²/s) single-bed cycle with a flow-condition inlet at the bottom and
  a pressure-condition inlet at the top. Per step and component, the moles the
  stream logs record entering the bed equal the bed's inventory change (before
  the fix: errors of 0.4–50 %), and at cyclic steady state out/in of every
  component is 1 within 5e-4, the log quadrature error (before: O₂ 1.053, N₂ 0.983).
• Verification tests adjusted: OpenSystemEnergyBalanceConvergesUnderGridRefinement
  asserted first-order convergence of an energy closure error that came from the
  inlet ghost upwinding; the error is now 1.5e-3 → 6.6e-5 over the four grids
  (at least first order is asserted, and < 2e-4 on the finest grid instead of
  < 5e-3). EndOfCycleProfileConverges allows an absolute 2e-4 mol/m³ wobble at
  the ~1e-3 mol/m³ floor the gaps level off at.
• Validation numbers that moved (example decks re-saved; tolerances
  re-measured where they no longer held, each documented in the paper's
  README):
  • Ward 2022 (species fix; the bed had taken in ~2 % more CO₂ than fed): CO₂
    front time vs experiment +1.9/+1.7/+0.1 % → +4.2/+3.3/+1.5 % at 50/100/200
    cm³/min, outlet-curve NRMSE 0.18/0.13/0.029 → 0.33/0.24/0.12; thermocouple
    RMSE 1.19/2.06/4.01 → 1.18/1.18/2.76 K.
  • Streb 2021 (flow-face fix): CH₄ fronts ~3 % earlier; CH₄ front-time error
    A2 7.6 → 10.5 %, A3 0.3 → 3.1 %, A5 12.1 → 14.2 %, A6 8.7 → 11.7 %,
    A7 7.3 → 10.2 %, B1 0.7 → 2.1 %; A3 y_CH₄ NRMSE 0.21 → 0.50; A6 y_CO₂ NRMSE
    0.048 → 0.10; A2 thermocouple RMSE 14.1 → 18.6 K.
  • Nguyen 2023: CO₂ recovery of IX and X no longer exceeds 100 %
    (100.35/100.78 → 98.27/97.64 %); purity X 80.84 → 79.67 %, IX 84.69 →
    84.28 %; VIII productivity 0.1566 → 0.1545 (−1.4 %).
  • Within their bounds: Perez 2019 A–F recovery down by up to 1.1 pp
    (purity ≤ 0.2 pp); Arora 2021 O₂ recovery at cycle 50 4.24 → 4.53 % (PSA)
    and 19.88 → 20.74 % (PVSA); Ferré 2026 CH₄ recovery up by up to 0.75 pp;
    Casas 2012 H₂ breakthrough times ≤ 1.2 % earlier; Nguyen 2022 outlet-curve
    NRMSE improves (Fig. 4 50 %: 0.061 → 0.042); Vysyaraju 2024 He recovery
    17.32 → 17.54 %; Ward & Pini 2024 ≤ 0.2 %. docs/validation.md and
    validation_summary.json are regenerated.

Changed

• Cyclic steady state: the mass-balance methods no longer stop by accident.
  On the way to CSS the per-cycle balance error swings through zero (a bed
  alternately gaining and losing a component), so the old single-cycle test
  with its 0.5 % band passed early: on the Yancy-Caballero (2020) 13X cycle at
  10 cells it stopped at the fifth cycle with purity and recovery 4–13 × 10⁻³
  from their converged values, and on Perez2019_A…_F 26 to 166 cycles
  early (KPIs up to 1.8 × 10⁻³ off). Mass Balance (Total) and Mass Balance
  (Per-Component) now also require the end-of-cycle bed state to have
  stopped changing: the largest relative RMS change since the previous cycle
  of the component loadings (scaled by the total loading), the mole fractions
  and the fluid, solid and wall temperatures
  (CyclicSteadyState::EndOfCycleStateChange) must be within the new
  css.state_tolerance (default 2 × 10⁻⁵; 0 switches it off). Every
  method's test must now hold on css.consecutive_cycles cycles in a row
  (default 3) before CSS is declared. The balance tolerance is unchanged: a
  tighter one cannot be relied on, because the balance carries a small
  systematic error of its own (up to 1 × 10⁻⁴ at CSS on Perez2019_E).
  With the defaults the Yancy cycle stops at cycles 59 and 100 (C1, C2,
  0-based) instead of 4, within 1.4 × 10⁻⁵ of the 150-cycle KPIs, and the
  example decks within 4.7 × 10⁻⁵ of their converged KPIs (Perez2019_F, which
  closes only ~1.5 % of its gap per cycle, 1.3 × 10⁻⁴); the stopping cycles
  and errors of 14 decks are tabulated in docs/cyclic-steady-state.md.
  CyclicSteadyState is class version 2: decks saved before it load with the
  defaults, so existing decks that use a mass-balance method run longer and
  report converged (different) purity and recovery, and Molar Profile decks
  stop at least two cycles later. css.state_tolerance = 0 with
  css.consecutive_cycles = 1 restores the old test. The solver's logs keep at
  least the cycles the check reads (one more than consecutive_cycles), even
  when settings.max_cycles_to_keep is smaller. CSSResult, CSSRecord and
  Python's sim.css_history() gain state_change, state_tolerance,
  passing_cycles and required_cycles; the solver log at level medium
  prints them after every cycle. Every paper case that runs to CSS moved
  (docs/validation.md): Perez 2019 A–F by up to 0.17 pp (A's recovery
  against the experiment crossed its bound, 0.294 → 0.417 pp; the bound is now
  0.5, see the Perez README), Nguyen 2023 I–X by up to 0.22 pp, Ferre 2026 by
  up to 0.32 pp, Ward & Pini 2024 by 0.04 pp and Vysyaraju 2024's He recovery
  by 0.64 pp (17.32 → 16.69 %); all within their bounds. Every bundled deck is
  re-saved with the two new fields, and the CSS decks with their new results.
• Adsorbent library: one recommended entry per material (catalogue schema 3).
  library.json groups its entries into 22 materials (id, name, family,
  recommended entry and why, aliases). A material's other entries are
  variants, each naming the product, form or data range that sets it apart;
  labels are the material name plus the variant ("Zeolite 13X (Zeochem Z10-02,
  CO2/N2 to 1 bar)") instead of "Material (Author year)". Each isotherm records
  its data_quality (measured with tabulated data, measured fit only, refit of
  other authors' data, copied from an earlier paper, not yet checked). The
  validator, the native library tests and the generated docs page (now grouped
  by material) enforce it. family moved from the entry to its material.
• AddLibraryAdsorbent, Simulation.add_library_adsorbent and
  skarstrom.adsorbent_path accept a material id (e.g. "zeolite-13x") and add
  or return that material's recommended entry; AdsorbentLibrary gains
  materials, FindMaterial, Resolve and Variants. Entry keys are unchanged.

Fixed

• WENO and van Leer advection (advection_scheme = weno | van_leer) are fixed; Upwind (the
  default, and every bundled deck) is unchanged. Compared line by line with the WENO of
  Haghpanah et al. (2013) / Yancy-Caballero et al. (2020) (PEESEgroup/PSA CycleSteps/WENO.m):
  (1) WENO's upwind-side candidate was the face value one face upstream, (y[n-1] + y[n])/2,
  instead of the extrapolation 3/2 y[n] - 1/2 y[n-1], which made the scheme first order;
  (2) both schemes gave the bed-end faces the average of the ghost and the first cell, but the
  ghost holds the inlet composition itself, so the bed received (y_feed + y_1)/2 while the
  stream logs, KPIs and cycle mass balance counted y_feed (CO2 balance 0.94-0.98 on the
  Yancy-Caballero 13X cycle at 10 cells with dispersion off, where upwind closes to 1.0000);
  (3) built component by component, the limited faces did not sum to the cell sum, and
  linearised about a steep front they weight the downwind cell, so a departure of the cell sum
  of mole fractions from one grew without bound (van Leer reached a sum of 90 in a CO2/N2
  breakthrough). Now the bed-end faces take the upwind value; the first interior face next to
  an inlet extrapolates from the inlet value half a cell away (2 y_1 - y_0, as the reference;
  van Leer doubles that slope likewise); void cells are treated as well mixed; WENO uses the
  cell widths on non-uniform grids; and the interior faces are rescaled so their mole fractions
  sum to the upwind cell's sum, which couples the mole fractions two cells apart, so the
  Jacobian half-bandwidth for these schemes grows from 2 neq to 2 neq + components - 1.
  On the Yancy-Caballero 13X cycle at 10 cells, WENO now gives a CO2 recovery of 90.0 % for
  the minimum-energy case (PSASimulator.jl WENO reference 89.7 %; was 77.5 %), and its CO2
  balance with dispersion off closes to 1e-5 at every grid. Tests:
  AdvectionSchemeConservation.EachSpeciesIsConservedThroughTheBedEnds and
  AdvectionSchemeConservation.CellMoleFractionsKeepSummingToOneAtASteepFront.
• CreateBoundaryValue and CreateInletTemperatureProfile know the double-exponential
  profile. The web document layer builds profiles through the first factory, so a deck
  with a double-exponential pressure boundary (the Perez 2019 decks) lost it to a null
  profile, i.e. a constant 1 bar, in the browser. A unit test now checks every listed
  kind has a factory entry.

Added

• A Blake–Kozeny momentum model (MomentumModel::BlakeKozeny = 2, builder kind and
  registry choice blake_kozeny, Python layer.set_momentum("blake_kozeny") /
  add_layer(..., momentum="blake_kozeny"), display name "Darcy (Blake-Kozeny)"):
  Darcy's law with K = ε³ d_p² / (150 (1 − ε)²), i.e. Ergun's viscous term without the
  inertial one. It is the pressure-drop law of Haghpanah et al. (2013) and Estupiñan
  Perez et al. (2019, SI Eq. S5). It takes the Darcy face path everywhere, like
  Kozeny–Carman. The stored values 0 (Kozeny–Carman) and 1 (Ergun) are unchanged, so
  existing decks load as before; a deck that uses the new value needs this engine.
• An optional constant molecular diffusivity for the macropore LDF coefficients
  (macropore, macropore_capacity, macropore_knudsen): registry leaf
  adsorbent[...].kinetics[<component>].molecular_diffusion_coefficient [m2/s], Python
  set_rate_model(..., molecular_diffusion_coefficient=1.6e-5) / kinetics[c].molecular_diffusion_coefficient.
  Above 0 it replaces the binary-diffusion D_m computed from T and P (which scales as
  1/P, so at 3-10 kPa evacuation pressures the coefficient was 10-33 times the value
  of the papers that use a constant D_m, e.g. Haghpanah et al. 2013 and Estupiñan Perez
  et al. 2019, Table 1). The default 0 keeps the previous behaviour bit for bit. Class
  versions bumped (macropore and macropore + Knudsen to 3, capacity factor to 2); older
  decks load with 0. Haynes is unchanged: its film term uses D_m with the Sherwood
  correlation, so a constant there is not a drop-in.
• Ideal Adsorbed Solution Theory (IAST) as a per-adsorbent mixing rule
  (AdsorptionSource::mixing_rule_, MixingRule::None|IAST; registry
  adsorbent[<id>].mixing_rule with choices none, iast; builder
  SetMixingRule and kind list MixingRuleKinds(); Python ads.mixing_rule,
  ads.set_mixing_rule(kind)). Under iast the mixture loadings come from the
  components' pure isotherms (any kinds, mixed freely), solved per cell by a
  bracketed Newton iteration on the reduced spreading pressure, warm-started
  from the cell's last solve. Closed-form spreading pressures for most kinds,
  Gauss–Legendre quadrature for Toth, Toth 2 and Redlich–Peterson. Absent
  components (y ≤ 0) take the infinite-dilution limit, so the loading is
  smooth through zero mole fraction. The evaluation functions are public in
  adsim/equations/iast.h and in Python as skarstrom.evaluate_iast. The
  solver log warns when IAST is combined with an extended (competitive) kind.
  The capacity-factor macropore coefficient (macropore_capacity), which
  evaluates q* itself, takes the IAST mixture loading under the IAST rule.
  AdsorptionSource is at cereal version 1; older decks load with none and
  run exactly as before. Docs: the isotherms page.
• IsCompetitiveIsothermKind(type) and PairedIsothermKind(type, out) in
  adsim/equations/isotherm_kinds.h: which kinds couple components through a
  shared denominator, and each one's pure counterpart with the same parameters
  (and back). The IAST warning uses them and names the pure kind to switch to.
• A double exponential profile (DoubleExponentialBoundaryValue,
  BoundaryValueType::DOUBLE_EXPONENTIAL, builder kind double_exponential, Python
  sk.double_exponential(end, rate_1, rate_2, weight, start=)): a two-term decay
  towards the end value, v(t) = v1 + (v0 − v1)[w e^(−r1 t) + (1 − w) e^(−r2 t)], the
  form Estupiñan Perez et al. (2019) fitted to their rig's evacuation pressure. The
  rates are entered in 1/s. Registry leaves <head>_double_exponential_start/_end/ _rate_1/_rate_2/_weight and pressure_double_exponential_end/_rate_1/_rate_2/_weight.
• A capacity-factor macropore LDF coefficient (MassTransferCoefficientMacroporeCapacity,
  kind macropore_capacity): k = (15 ε_p D_m / (τ r_p²)) · c / (ρ_p q*), the
  Haghpanah et al. (2013) form the Rajendran-group papers use. The ratio is evaluated
  from the isotherm at max(c, c_floor), with the isotherm re-evaluated at the floored
  concentration, so it stays finite where c and q* both vanish; capped at 1. The plain
  macropore kind is unchanged (it never applied the factor; its ratio = 1.0 stub is
  gone).

Changed

• Native WebAssembly exceptions by default. A new CMake option, ADSIM_EMSCRIPTEN_EXCEPTIONS, chooses the C++
  exception scheme of an Emscripten build: wasm (the new default, -fwasm-exceptions) or js (-fexceptions,
  the previous behaviour). The whole engine tree is compiled with the chosen flag, and Adsim::adsim now passes it to
  the consumer's compile and link. On the Skarstrom web solver, native exceptions ran 10 % faster on six example decks
  (4–17 %) with a 12 % smaller .wasm; they need Chrome 95, Firefox 100, Safari 15.2 or Node 17. Consumers: drop
  your own exception flags. -sNO_DISABLE_EXCEPTION_CATCHING=1 is rejected by Emscripten together with
  -fwasm-exceptions; keep the old scheme with -DADSIM_EMSCRIPTEN_EXCEPTIONS=js. The subproject consumer test now
  checks that an exception thrown inside the engine reaches the consumer's catch, and CI runs it with each scheme.
• CI builds the Emscripten consumer with Emscripten 6.0.13 (was 4.0.10).
• Perez 2019 example decks (Perez2019_A/B/C1/D1/E/F) now use the paper's model
  forms, after an audit of the deck against the paper and its SI
  (core/test/papers/perez_2019/README.md, provenance and attribution tables): the
  two-term exponential evacuation pressure of Sect. 4.3 (0.019 and 0.101 1/s; the split
  of the terms is not printed and an equal split is used, chosen on A), wall ends at
  ambient (SI Table S2) and the capacity-factor LDF of SI Eq. S4. The four-step
  recoveries, 8-10 pp above the experiments before, are now within 0.7-1.1 pp for A and
  B; the attribution shows the evacuation pressure history was the cause. The decks keep
  Kozeny-Carman: the paper calls the pressure drop negligible, Ergun is worth 0.3 pp, and
  it fails the initial-condition calculation of the fast-feed connected configuration in
  ConfigD1_ConnectedBC (an open engine item). examples/perez_2019.py builds the same
  decks.

Changed

• Ward 2022 example decks (Ward2022_Q50ccm/Q100ccm/Q200ccm) now use the
  paper's model forms, after an audit of the deck against the paper text
  (core/test/papers/ward_2022/README.md, provenance and attribution tables):
  helium as the carrier (the CO2/He experiment), the paper's dispersion
  coefficient (Eq. 4 with D_m = 0.605e-5 m2/s, times the factor 4 the paper
  fitted), -dU_1 + RT as the heat of adsorption (32.4 kJ/mol), the wall ends
  held at ambient (Table 4) and the Ergun pressure drop. The CO2 front now
  arrives 0.1-1.9 % late instead of 4-6 % early, the curve NRMSE is 3-18 %
  instead of 34-38 %, and the temperature rise is within 0.9-2.9 K of the
  measurement (it was 8-11 K high). The Ward test's temperature-rise
  indicator follows the paper's Table 9 (against the feed temperature).
  examples/ward_2022.py builds the same decks.

Added

• A log-linear profile (LogLinearBoundaryValue, BoundaryValueType::LOG_LINEAR,
  builder kind log_linear, Python sk.log_linear(end, start=, time=)) for pressures,
  flows and inlet temperatures: v(t) = v0 (v1 / v0)^(t / t_r), a geometric ramp whose
  relative rate of rise d(ln v)/dt is constant. A pressurisation from vacuum with a
  linear ramp puts most of the relative pressure rise, and with it the inlet pressure
  gradient and approach to fluidisation, in the first seconds; the log-linear ramp
  spreads it over the step. It holds its end value after the ramp, and ramps linearly
  when the start or end value is zero or negative. Registry leaves
  <head>_log_linear_start/_end/_time and pressure_log_linear_end.

Changed

• Solver residual, byte-identical: the competitive concentration-based
  Langmuir isotherms (ExtendedLangmuirIsotherm_Concentration,
  ExtendedDualSiteLangmuirIsotherm_Concentration, the kinds most example
  decks use) gain the single-pass UpdateAll the pressure-based ones had, so
  their shared denominators and exponentials are computed once per cell
  instead of once per component (Illustrative TSA: 17% fewer instructions for
  the run, the isotherm 2.6x cheaper). Every competitive isotherm's
  per-component Update is now that pass for one component, which also
  serves the capacity-factor MTC that re-runs the isotherm. The binary
  diffusion model memoises the Neufield collision integral and T^1.5 per
  cell on the exact temperature bits (TemperatureCache), since most residual
  evaluations repeat the previous temperature in most cells (15% fewer
  instructions in that model).
• Solver residual, byte-identical: the Langmuir and dual-site Langmuir isotherm
  functions (and their concentration-based variants) compute each exp(q/(R T))
  factor once instead of twice (GCC's default -fmath-errno does not merge the
  duplicate calls, so Linux builds paid 4 exp per cell per component), and a
  constant mass-transfer coefficient is evaluated once per layer instead of once
  per cell (ExpressionParam::IsRunConstant()). -17% instructions per residual
  call on FiveStepProductPurge; CSV output identical on six example decks.
• Solver: two settings are now on by default (IDASettings, transient):
  • suppress_alg_ with pressure_differential_: the pressure is tagged a
    differential unknown (its residual carries dP/dt) and IDA's local error
    test leaves out the algebraic entries, i.e. the ghost cells (mirror
    relations with no time derivative) and the face velocity. Their jumps at
    a boundary-condition switch used to fail the error test and force tiny
    restart steps; the cyclic example decks now take 40-60% fewer steps
    (1.3-2.1x wall) at the same accuracy against a reltol = 1e-8 reference.
  • alias_lte_solid_temperature_ and derived_concentration_: with local
    thermal equilibrium the solid temperature is a copy of the fluid
    temperature instead of an unknown, and the gas concentration is derived
    from P and T (C = P / RT; no residual reads dC/dt). Two unknowns per cell
    fewer, a narrower band, one colour fewer (1.1-1.25x).
    The equations solved are unchanged; the integration path is not, so
    results move in their last digits (the example decks are re-saved). The
    per-step boundary/outlet logs have fewer points than before, so a value read
    from them at an exact time is interpolated over a longer step. Set the
    switches to false for the previous behaviour.

Changed

• Solver: the banded Jacobian matrix is sized from the coupling reach a
  difference-quotient probe measures at the initial state (12 / 16 instead of
  2 neq = 20 / 20 diagonals at ten equations per cell), about 45% less LU
  work; results are byte-identical. If a Jacobian audit later finds a
  coupling outside the band, the solver widens it to the full reach and
  restarts from the current point (logged at medium level).
  IDASettings::band_from_probe_ switches it off.

Added

• Transient (not serialised) IDASettings switches for solver experiments,
  all off by default: alias_lte_solid_temperature_ and
  derived_concentration_ (drop Ts and C from the DAE unknowns),
  suppress_alg_ (IDASetSuppressAlg), delta_cj_lsetup_
  (IDASetDeltaCjLSetup) and pressure_differential_ (tag P as a
  differential unknown). bench_solver gets --unknowns, --band,
  --suppress-alg, --dcj, --p-diff and --reltol; its report adds the
  band widths, out-of-band count and whether the band was widened.
• Engine internals de-duplicated (no behaviour or deck-format change): the
  twenty concrete isotherm classes share one ParamIsotherm<Params, kType> base and
  only declare their Update (archives are byte-identical: the parameter map keeps
  its name and each class its cereal_class_version); the four flow boundary value
  kinds keep their profile_ and serialisation in IFlowBoundaryValue; the
  per-stream GetComponentMolarFlow / GetEnthalpyFlow, which were copied into
  every stream kind, live once in IBoundaryStream with a BedInletTemperature
  hook for the pump-driven kinds' compression heating; EquationLog / BoundaryLog
  getters and face interpolation, the KPI bed-accumulation integrals and the
  compressor / vacuum-pump cost structs each have a single implementation. The
  unused CalculateAccumulation_Volume is gone. Test fixtures that were copied
  across files (the validation PSA cycle, the feed/hold integration deck, the
  Ogata-Banks solution) are shared headers now.

What's Changed

• Ergun: boundary faces use the interior neighbour's gas state; warn once per run about constant-pressure step changes by @samaffleck in #9
• Solver log: drop the end-of-run IDA statistics line by @samaffleck in #10
• Start 0.2.0 development after the 0.1.0 release by @samaffleck in #11
• Install consumer: request the tree's own version, not a fixed 0.1 by @samaffleck in #12
• Adsorbent library: drop the legacy Zeolite 13X entry and bulk density by @samaffleck in #13
• Sensitivity analysis and KPI error bars by @samaffleck in #14
• Give the polynomial pump efficiency its P^2 term by @samaffleck in #15
• Release 0.2.0 by @samaffleck in #16
• Upgrade cibuildwheel and pip to address certificate authority issue by @samaffleck in #17
• Allow expressions in inlet and initial compositions by @samaffleck in #18
• Wall pressure design, operating-limit warnings and warning switches by @samaffleck in #19
• Release 0.3.0 by @samaffleck in #20
• Size the Jacobian band from a coupling probe; solver experiment switches by @samaffleck in #21
• Solver: suppress algebraic error-test entries and reduce the DAE unknowns by default by @samaffleck in #22
• Add a log-linear boundary value profile by @samaffleck in #23
• Residual: one exp per temperature factor in the Langmuir isotherms; constant MTC evaluated once per layer by @samaffleck in #24
• Validation: Ward 2022 on the paper's model forms, provenance audits for Ward 2022 and Streb 2021 by @samaffleck in #25
• Perez 2019 on the paper's model forms: two-term evacuation profile and capacity-factor LDF by @samaffleck in #26
• Residual: single-pass UpdateAll for the concentration-based competitive Langmuir isotherms; memoised Neufield table by @samaffleck in #27
• Profile factories: add the double-exponential kind by @samaffleck in #28
• Remove duplicated code in the engine and the test suite by @samaffleck in #30
• Solve the Ergun flow-boundary face explicitly at step switches by @samaffleck in #32
• Add Ideal Adsorbed Solution Theory as a per-adsorbent mixing rule by @samaffleck in #33
• Adsorbent library: one recommended entry per material (catalogue schema 3) by @samaffleck in #35
• Blake-Kozeny momentum model and constant D_m for the macropore LDF coefficients by @samaffleck in #34
• Close the component mass balance with axial dispersion (Danckwerts inlets) by @samaffleck in #37
• Fix the WENO and van Leer schemes; CSS requires a settled bed state by @samaffleck in #36
• Emscripten: native WebAssembly exceptions by default; CI on Emscripten 6.0.13 by @samaffleck in #38
• Molar-basis gas heat capacity; series and profiles at the right positions by @samaffleck in #39
• Start every cycle step with a fresh IDA history (WASM/native parity) by @samaffleck in #40
• Fix face interpolation weights in EquationLog::GetFaceVariable by @samaffleck in #42
• Mole-conserving connected streams; faster CI (one run per PR, cached Windows build, paper tests before merge) by @samaffleck in #41
• Union merge for CHANGELOG.md; CI skips builds for docs-only PRs by @samaffleck in #43
• WP15-E7: list every verification-suite file in the README Files table by @samaffleck in #44
• WP15-E1: LDF-limited front against the Anzelius/Thomas solution by @samaffleck in #46
• WP15-E2: Langmuir desorption simple wave against the method of characteristics by @samaffleck in #47
• WP15-E3: binary competitive Langmuir fronts against equilibrium theory by @samaffleck in #48
• WP15-E4 (part): blowdown composition path of a binary on linear isotherms by @samaffleck in #49
• WP15-E5: thermal front speed and steady wall-loss temperature profile by @samaffleck in #50
• WP15-E6: evacuation, two-bed pressure equalisation and buffer-tank fill by @samaffleck in #51
• WP17-E3: metamorphic relations (component order, layers, split steps, units, cycle phase) by @samaffleck in #52
• WP17-E1, WP17-E2: parameter limits (rates of approach) and degenerate mixtures by @samaffleck in #53
• WP17-E4: thermodynamic bound on the separation energy at CSS by @samaffleck in #54
• WP17-E5: extreme operating corners converge or fail cleanly by @samaffleck in #55
• Fix: connected streams replay a bottom pressure boundary's outflow with the wrong sign by @samaffleck in #56
• WP18-E2: cross-code regression against the PEESE model and toPSAil by @samaffleck in #57
• Save a sensitivity study's set-up with the deck by @samaffleck in #59
• Constant gas heat capacity and viscosity; heat of adsorption as −ΔU + RT by @samaffleck in #45
• Fix: the ColouredCheck Jacobian diagnostic learns the pattern without recolouring by @samaffleck in #60
• WP19-E3: valves can read the bed side of their law at the end cell by @samaffleck in #61
• WP19-E4: extended Langmuir-Freundlich isotherm (toPSAil form) by @samaffleck in #62
• WP16-E1..E4, E7: convergence verification of cycle results (grid, integrator tolerance, CSS) by @samaffleck in #58
• Fix: WENO / van Leer outlet-cell mole-fraction sum drifts during fast outflow by @samaffleck in #63
• WP16-E6: method-of-manufactured-solutions hook and order test by @samaffleck in #65
• Own ramp time for pressure ramps; user exponential factor by @samaffleck in #68
• Fix: breakthrough decks stall on fine grids (Ward2022_Q50ccm at 320 cells) by @samaffleck in #66
• Keep the KPIs of every cycle, not only the logged ones by @samaffleck in #69
• Fix: the CSS state test never passes for an inert bed by @samaffleck in #71
• Release automatically when a version change reaches main by @samaffleck in #70
• Release 0.4.0 by @github-actions[bot] in #72

New Contributors

• @github-actions[bot] made their first contribution in #72

Full Changelog: v0.1.0...v0.4.0