Numerics

Numerics

This page describes how the solver turns the governing equations into numbers: the spatial discretisation, the time integrator and its settings, how the Jacobian is built, and how cyclic steady state is detected and accelerated. Everything here is taken from the code: core/src/equations/adsorbent_equations.cpp (residuals), core/src/equations/equation.cpp (unknowns and layout), core/src/solver/coupled_solver.cpp and core/include/adsim/solver/ida_settings.h (integration), core/src/solver/jacobian.cpp (Jacobian) and core/src/solver/cyclic_steady_state.cpp (CSS).

---

Spatial discretisation: cell-centred finite volumes

The column is divided into cells along its axis. Each layer contributes its Number of cells of equal width \(\Delta z = L/N\); optional bottom and top void regions add one lumped cell each; and one ghost cell at each end carries the boundary condition. Layers of different cell widths can sit next to each other.

The cell-centred unknowns, in every cell, are:

Unknown	Residual it carries
total concentration \(C\)	ideal gas law, \(C - P/(RT) = 0\) (Equation of State)
pressure \(P\)	overall mass balance (Mass Balance)
gas temperature \(T\)	bed energy balance (Energy Balance)
solid temperature \(T_s\)	solid energy balance in LTNE; \(T_s = T\) in LTE
wall temperature \(T_w\)	wall energy balance
mole fraction \(y_i\) of every component	component mass balance
adsorbed concentration \(q_{i,k}\) of every component on every adsorbent \(k\)	the kinetic rate model, only in the cells of layers that contain adsorbent \(k\)

In an isothermal run the three temperatures are removed from the unknowns, and void cells carry no solid temperature or loadings.

Where results are reported. A ghost-cell value is not a value anywhere in the bed (an inlet temperature ghost holds \(2T_{in} - T_1\), so that the face between it and the first cell sits at \(T_{in}\)), so the results never report one. Axial profiles of cell quantities are given at the cell centres, and of face quantities (velocity, flow rates) at the faces. The bed ends report the physical boundary value: the composition, temperature and pressure an inlet imposes, or the bed-end face value at an outlet, which is what the boundary-stream logs record (Results, adsim/solver/bed_series.h).

Staggered face velocities. The superficial velocity \(u\) is not stored at cell centres but on the faces between neighbouring cells, including the faces between the ghost cells and the bed. Each face carries the momentum balance as its residual, so the flux through every face is computed from the pressures on its two sides without interpolation. The finite-volume balances then use these face velocities directly: for example the convective term of the component balance in cell \(n\) is the difference of the face fluxes

\[ F_{n+1/2} - F_{n-1/2},\qquad F_{n+1/2} = \frac{T_n}{P_n}\thinspace \frac{y_{i,n+1/2}\thinspace P_{n+1/2}\thinspace u_{n+1/2}}{T_{n+1/2}} \]

with \(P\) and \(T\) at a face taken as the mean of the two cells. This staggering avoids the checkerboard pressure modes a collocated grid can develop, and makes the scheme conservative: what leaves one cell through a face enters its neighbour.

Dispersion and conduction fluxes are central differences. At a face between cells of different width or different dispersion coefficient, the face coefficient is the cell-width-weighted harmonic mean of the two cell values, so a sharp change of packing is handled consistently.

Boundary faces (Danckwerts). No dispersive or conductive flux crosses the two boundary faces of the bed: the flux through a boundary face is the convective one alone, evaluated with the boundary value. At an inlet the component flux is therefore exactly \(u\thinspace C\thinspace y_{i,in}\) and the enthalpy flux \(\rho\thinspace c_p\thinspace u\thinspace T_{in}\), the amounts the stream logs, the performance metrics and the cyclic-steady-state mass balance record; this is the finite-volume form of the Danckwerts inlet condition \(u\thinspace C\thinspace y_{in} = u\thinspace C\thinspace y(0^+) - \varepsilon_b D_{ax} C\thinspace \partial y/\partial z\rvert_{0^+}\), and the first cell's composition can differ from the feed's. At an outlet or a closed end the ghost cell mirrors the first cell (zero gradient), so there was no flux to drop. The species ghost cell holds the boundary composition itself; the temperature ghost cell is set so that the mean of the ghost and the first cell, the face temperature, is the boundary temperature, and the energy balance takes that face temperature as the upstream value of an inlet face.

Momentum: Darcy and Ergun as face residuals

For a face between cells \(n\) and \(n+1\) the two half-cells act in series, each with its own viscous resistance \(\mu\thinspace \Delta z/2K\) and, for Ergun, inertial resistance \(\beta\rho\thinspace \Delta z/2\). Darcy faces (Kozeny–Carman or Blake–Kozeny, no inertial term) use the explicit residual \(u - \Delta P / A = 0\). Ergun faces are solved implicitly: instead of imposing the root of the quadratic \(\Delta P = A u + B\lvert u\rvert u\), which behaves like a square root near zero flow, each face carries the residual

\[ r_u = u - \frac{\Delta P - B\thinspace \lvert u\rvert\thinspace u}{A} = 0 \]

which is smooth through flow reversal, linear in the pressures and has a slope of at least 1 in \(u\). Its solution is exactly the Ergun root. The boundary face of a pressure or closed boundary condition, which fixes that face's pressure drop, instead imposes the explicit Ergun root of the boundary pressure drop: when a step switch makes that face's velocity jump, the implicit residual's slope in \(u\) changes by an order of magnitude and IDA's corrector (whose Jacobian predates the jump) fails to converge, while the explicit root is linear in \(u\) and depends only on the first-cell pressure, which is continuous across the switch. Flow boundaries, which fix the face velocity, keep the implicit residual. With this split Ergun decks cross step switches as cleanly as Darcy ones. The Momentum Balance page gives the details, including void regions and flow boundaries.

Advection schemes

The mole fraction at a face, \(y_{i,n+1/2}\), has to be reconstructed from the cell values. Three schemes are implemented; they change only this face value, never the conserved cell values, so total species mass is conserved whichever is used:

• Upwind (AdvectionScheme::Upwind): the face takes the value of the upstream cell. First order, monotone and robust; it smears sharp fronts over a few cells, so it needs many cells (below). It was the default up to engine 0.4, and a deck saved before the setting existed (Solver_Settings version 1) still runs with it.
• WENO (AdvectionScheme::WENO, the default): the third-order WENO scheme of Jiang and Shu (two stencils, \(k = 2\)) used by Haghpanah et al. (2013) and Yancy-Caballero et al. (2020). For flow in the \(+z\) direction it blends the central interpolation \(\tfrac12(y_n+y_{n+1})\) and the upwind extrapolation \(\tfrac32 y_n - \tfrac12 y_{n-1}\) with linear weights 2/3 and 1/3, \(y_{n+1/2} = \dfrac{\alpha_0\thinspace\tfrac12(y_n+y_{n+1}) + \alpha_1\thinspace(\tfrac32 y_n - \tfrac12 y_{n-1})}{\alpha_0+\alpha_1}\), \(\alpha_0 = \dfrac{2/3}{(\beta_0+\varepsilon)^2}\), \(\alpha_1 = \dfrac{1/3}{(\beta_1+\varepsilon)^2}\), with smoothness indicators \(\beta_0 = (y_{n+1}-y_n)^2\), \(\beta_1 = (y_n-y_{n-1})^2\) and \(\varepsilon = 10^{-12}\) (the reference codes write \((\Delta y + 10^{-10})^4\), which is the same away from flat profiles). With the linear weights this is the third-order upwind-biased reconstruction \(-\tfrac16 y_{n-1} + \tfrac56 y_n + \tfrac13 y_{n+1}\); near a front the weight shifts to the smoother stencil. On a non-uniform grid (layers of different cell size) both candidates are the linear interpolation or extrapolation in \(z\), and \(\beta\) is the slope times the cell width.
• van Leer (AdvectionScheme::VanLeer): a limited second-order reconstruction \(y_{n+1/2} = y_n + \tfrac12\thinspace \phi(r)\thinspace (y_{n+1}-y_n)\) with \(r = (y_n - y_{n-1})/(y_{n+1}-y_n)\) and the van Leer limiter \(\phi(r) = (r + \lvert r\rvert)/(1 + \lvert r\rvert)\). The ratio is evaluated as \(r = ab/(b^2 + 10^{-24})\) with \(a = y_n - y_{n-1}\), \(b = y_{n+1} - y_n\), which is the same for \(y\) and \(1-y\) and goes smoothly to zero on a flat profile.

Boundaries. The ghost cell holds a value at the bed-end face, not at a cell centre: the inlet composition (Dirichlet), or a copy of the end cell at an outlet. For the higher-order schemes the two bed-end faces therefore take the upwind value, exactly as Upwind does, so the convective flux through a bed end is the stream's own composition times its flow and equals the stream integral the KPIs and the cycle mass balance are built from. In WENO the first interior face next to an inlet extrapolates from that face value, half a cell away: \(2y_1 - y_0\) with \(\beta_1 = \left(2(y_1-y_0)\right)^2\), as in the reference implementations; van Leer doubles the upwind slope in \(r\) there in the same way. A void region is a well-mixed tank, so both schemes treat it in the same way (its faces are upwind and its value sits on the face of the next cell). The reconstructed face value is bounded to \([\min(0, y_n, y_{n+1}),\ \max(1, y_n, y_{n+1})]\), the physical range widened to the two cells that share the face, so an overshoot cannot inject a negative or super-unity composition: WENO's upwind extrapolation overshoots where a cell is a local extremum (the cells 0.45, 0, 1 give a face value of \(-0.17\)), while the van Leer face lies between its two cells and is never changed. For cells inside [0, 1] this is the plain [0, 1] clamp. It is widened because a component absent from the feed and the bed holds only round-off around zero (and a pure one round-off around one): clamped to [0, 1], the face of a cell marginally below zero was zero, with a zero derivative, so negative round-off could not be carried out of the bed while positive round-off was, the absent component drifted to a few percent of its absolute tolerance, and the face switching between clamped and free during the Newton iteration made it fail to converge. WENO then took 25 to 80 times longer than Upwind on such decks (van Leer 4 to 15 times), nearly all of it in steps retried after Newton convergence failures; with the widened bound it takes the same steps as Upwind. A face value between the two cells it separates is not an overshoot, so allowing it keeps round-off behaving as under Upwind. The faces that take the upwind value (the bed ends and the faces of void cells) are not bounded at all, exactly as under Upwind: a bed-end face is not rescaled (below), so clamping it would make the outlet face of a cell holding, say, \(y_{N_2} = 1 + 10^{-9}\) sum to one instead of to the cell's sum, and that cell's departure of the sum would no longer be carried out of the bed. Every face is single-valued, so each species is conserved whatever the scheme.

Face temperature. The energy balance's convective face temperature is reconstructed with the same scheme, the same stencils and the same upwind faces (bed ends, void cells) as the mole fractions, so the temperature has the scheme's order of accuracy too. Two details differ. The temperature ghost cell does not hold the face value, so wherever a mole-fraction reconstruction reads the ghost (the upstream value of a bed-end face, and the far cell of the first interior face) the temperature reads the face temperature, the mean of the ghost and the first cell, which is the boundary temperature at an inlet. The face temperature is not bounded. And the reconstruction acts on \(T/T_c\), with \(T_c\) the upwind cell's temperature, scaled back afterwards: WENO's \(\varepsilon\) and van Leer's \(10^{-24}\) are absolute and sized for mole fractions, and on kelvin they would let noise-level differences of a nearly flat temperature profile switch the weights. The energy balance is written in advective form, \(\rho c_p\thinspace \bar u\thinspace (T_{n+1/2} - T_{n-1/2})\) with the cell's own \(\rho c_p\) and mean face velocity, so the face value enters it as a difference, as before. Before this was introduced (see the changelog) the face temperature was always the upwind value, whatever the scheme, which kept the temperature first order (see the order of accuracy below).

Sum of the face mole fractions. Every species and the total have their own balance, so nothing in the equations pulls the cell sum \(\sum_i y_i\) back to one; the species fluxes must add up to the flux the total balance carries. Upwind faces sum to the upwind cell's sum, so a departure of the cell sum (at the level of the solver tolerance) is simply carried along and out of the bed. The nonlinear reconstructions, built component by component, do not keep that sum, and worse, linearised about a steep front they weight the downwind cell, so a departure of the cell sum grows instead of being carried: before this was handled, van Leer drove the cell sum to 90 in a CO2/N2 breakthrough on 13X. The interior faces of WENO and van Leer are therefore rescaled so that their mole fractions sum to the upwind cell's sum (blended through zero flow like the upwind value), and the bed-end faces take the upwind cell's values unclamped, so every face sums to its upwind cell's sum. Summing the species balances of a cell and subtracting its sum \(\Sigma_n\) times the total balance gives \(\varepsilon_t \Delta z\, \mathrm{d}\Sigma_n/\mathrm{d}t = -(\Sigma_e - \Sigma_n) F_e + (\Sigma_w - \Sigma_n) F_w + \dots\), with \(\Sigma_e\), \(\Sigma_w\) the sums of the face mole fractions and \(F\) the face fluxes: with upwind sums the departure of \(\Sigma\) is carried from cell to cell and out of the bed. The sum then behaves exactly as under upwinding while each component keeps its high-order profile; each face stays single-valued, so each species is still conserved. This couples every mole fraction across the stencil, which widens the Jacobian band (below).

Smoothed upwinding at zero flow. A hard switch of the upwind cell when \(u\) changes sign makes the residual discontinuous, and Newton iterations chatter when a face sits at zero flow (a closed bed end, a bed that has just reached the boundary pressure). The face value therefore blends the two candidates with the weight \(w = \tfrac12\left(1 + \tanh(u/u_s)\right)\), \(u_s = 10^{-6}\) m/s. For \(\lvert u\rvert > 20\thinspace u_s\) this is exactly the hard upwind choice. Flow boundaries blend their inlet composition and temperature in the same way as the boundary flow passes through zero.

Choosing the scheme and the number of cells. A grid study of every example deck (each run on 10 to 160 cells per layer, Richardson extrapolation from the three finest grids) found that recovery within 0.1 percentage point of its grid-converged value needs 200-600 or more upwind cells on the cyclic decks (Perez2019 A-E 200 to over 1000, Nguyen2023 170-340, Ward2024 360, Ferre2026 610), but 25-110 WENO cells (Perez2019 23-65, Nguyen2023 40-76, Ward2024 94, Ferre2026 110). A WENO run costs about 1.5 times an upwind run on the same grid (more steps and a wider Jacobian band), so at equal accuracy WENO is 3-10 times cheaper; that is why it is the default. Fronts of a strongly adsorbed component converge on 20-40 cells under every scheme, while the front of a weakly adsorbed one (H₂ or CH₄ in a 1.2 m column) needs 100-190 WENO cells. The example decks use grid-converged counts with WENO: 80 cells for the 0.3 m VSA columns (Perez2019, Nguyen2023, Vysyaraju2024), 100 for Ward2024 and Arora2021 PVSA, 120 for Ferre2026, 160 for Arora2021 PSA (its 7 % net O₂ recovery is very grid-sensitive) and for the Casas2012 and Streb2021 breakthroughs; the short breakthrough columns (Nguyen2022, Ward2022) and the illustrative TSA are converged on their 30-60 cells. Upwind remains the scheme of choice where a first-order reference is wanted (a tanks-in-series model, an analytic solution with the upwind numerical dispersion), and van Leer is slightly more accurate per cell than WENO at the same cost. Three caveats: WENO is not strictly monotone, so cell mole fractions and loadings can undershoot zero (or overshoot one) by up to a few times the absolute tolerance of 10⁻⁸ (−4×10⁻⁸ in a Perez2019_A run) where upwind stays at round-off; its smoothness weights use an absolute \(\varepsilon\), so they depend on the scale of the profile and two identical components fed in different amounts drift from their feed ratio by about 10⁻⁵ (upwind keeps it to round-off); and sharper fronts make the approach to cyclic steady state slower on decks whose CO₂ inventory settles slowly (light-product pressurisation cycles), which the example decks offset by starting from a feed-loaded bed. Whatever the scheme, check a result by re-running it on twice the cells.

The scheme is a solver setting stored in the deck ("Advection_Scheme" in Solver_Settings: 0 = Upwind, 1 = WENO, 2 = van Leer) and in the C++ API (IDASettings::advection_scheme_). It is also the parameter settings.advection_scheme (an integer) in the Python parameter registry.

Order of accuracy. The method of manufactured solutions measures the order of the whole discretisation with every term active at once: a smooth, time-dependent solution for \(P\), \(y\), \(u\), \(q\) (and \(T\)) is substituted into the continuous equations, the remainder is added to the residuals as a source through a verification-only hook (UserData::manufactured_sources_, adsim/equations/manufactured_sources.h: not serialised, not reachable from a deck, Python or the web, and with no source set the results are bit-identical), and the discrete solution is compared with the chosen one on 20 to 160 cells (verification). Isothermally, with dispersion, LDF kinetics on a Langmuir isotherm and Darcy flow, every field converges at order 1 with Upwind and order 2 with WENO and van Leer: the central dispersion, pressure and Darcy terms and the cell-centre sampling are second order, so WENO's third-order interior reconstruction shows only on coarse grids (interior error order 2.5 to 2.8 there). The bed-end faces take the upwind value in every scheme, which is first order wherever the profile has a slope at the bed end (as at the outlet of a breakthrough, where it caps the outlet-curve order at 1); with the zero end gradients the Danckwerts conditions impose on a dispersed profile it is second order. With the energy balance on, the temperature converges at the same orders (1 with Upwind; 2.0 to 2.25 with WENO and 1.95 to 2.03 with van Leer on 20 to 160 cells, still 2.0 from 320 to 640 cells), and the other fields keep their isothermal orders. When the face temperature was always the upwind value it converged at first order for every scheme, and on fine grids that error pulled the other fields of WENO and van Leer towards first order too.

Time integration: SUNDIALS IDA

After discretisation the model is a differential-algebraic system \(F(t, \mathbf{y}, \dot{\mathbf{y}}) = 0\), integrated with IDA from SUNDIALS, a variable-order, variable-step backward differentiation formula (BDF) solver. At each step IDA predicts the new state from a polynomial through the previous points, then corrects it with a Newton iteration on the implicit BDF equations.

Settings (Solver → Solver Settings; IDASettings):

Setting	Default	Meaning
Max BDF order	5	highest order IDA may use (1–5)
Relative tolerance	10⁻⁶	one scalar for every unknown
Initial step	10⁻¹⁰ s	first step attempted, and after every restart
Min step	10⁻¹⁵ s	IDA fails rather than go below it
Max step	0 = no limit	cap on the internal step size
Max nonlinear iterations / step	15	Newton iterations per step attempt
Nonlinear convergence coef	0.1	safety factor on the Newton convergence test
Max convergence failures	200	Newton failures allowed in one step
Max error-test failures	40	local-error failures allowed in one step
Max number of steps	100 000	IDA's limit between two output times

Absolute tolerances are a vector, one value per kind of unknown, fixed in the code: \(10^{-8}\) for the total concentration [mol m⁻³], the mole fractions and the loadings [mol kg⁻¹]; \(10^{-6}\) for the pressure [Pa] and the face velocities [m s⁻¹]; \(10^{-4}\) K for the three temperatures. IDA accepts a step when the weighted norm of its local error estimate, with weights \(1/(\text{rtol}\thinspace \lvert y\rvert + \text{atol})\), is below 1.

Differential and algebraic unknowns. IDA is told which unknowns are algebraic (pressure, the face velocities and every ghost-cell value, which only appear through boundary conditions) and which are differential. This is what IDACalcIC needs.

Consistent initial conditions. The initial conditions you enter set the differential unknowns but rarely satisfy the algebraic equations (the boundary conditions, the momentum balance, the ideal gas law). Before the first step, IDACalcIC in IDA_YA_YDP_INIT mode solves for the algebraic unknowns and the time derivatives of the differential ones so that \(F = 0\) holds, and the corrected state is what gets logged at \(t = 0\). If it fails, the run stops with Initial condition calculation failed, which almost always means the initial conditions are far from what the first step imposes (see Running the Solver).

Steps of a cycle. Each cycle step sets IDA's stop time to the end of the step and integrates to it one internal step at a time (IDA_ONE_STEP), so a step boundary is never stepped over. The new step's boundary conditions take effect when the step starts. Results are written at up to 200 points per step (max_log_points_per_step_), evenly spaced in time and filled in from IDA's interpolating polynomial, plus the first and last point of every step; the boundary flows that feed the KPIs are logged at every accepted step. Only the last Max cycles to keep cycles (default 10) are retained for plotting.

Recovery from failures. If IDA fails on a recoverable error (corrector or error-test failures, a residual or linear-solver error, too much work), the step is retried up to four times. Each retry restarts IDA from the last accepted point with a fresh history, zeroed algebraic derivatives and a new IDACalcIC, which clears the large derivative estimates an abrupt boundary switch can leave in the history. The log reports each retry; after the fourth the run fails with the IDA error and the unknown with the largest local error.

The linear solver and the Jacobian

Each Newton iteration solves a linear system with the iteration matrix \(J = \partial F/\partial \mathbf{y} + c_j\thinspace \partial F/\partial \dot{\mathbf{y}}\). The unknowns are ordered cell by cell (all unknowns of cell 0, then cell 1, and so on), so every equation couples only to its own cell and its neighbours and \(J\) is banded. Cross-variable couplings reach one cell either side, so a half-bandwidth of twice the number of unknowns per cell covers the upwind scheme. The WENO and van Leer stencils reach two cells, and their interior face mole fractions are rescaled over the components, so a mole fraction there reads every mole fraction two cells away; for these schemes the half-bandwidth is twice the number of unknowns per cell plus the number of components minus one. SUNDIALS' banded direct LU solver (SUNLinSol_Band) factorises it.

The Jacobian itself is approximated by difference quotients: perturb a group of unknowns, re-evaluate the residuals, and read off the columns. How the groups are chosen is the Jacobian mode (JacobianMode, a solver implementation detail that is not saved in decks; the bench_solver tool selects it with --jacobian ida|band|coloured):

Mode	What it does
IdaInternal	IDA's own banded difference-quotient Jacobian: columns a full bandwidth apart are perturbed together, so each Jacobian costs one residual evaluation per column of the band (upper + lower half-bandwidth + 1).
BandDQ	Skarstrom's own replica of IDA's banded algorithm, with exactly the same increments and column groups, so every entry is bit-identical to IdaInternal. Used as the reference for the coloured mode.
Coloured (default)	Difference quotients over a learned sparsity pattern. Columns whose rows do not overlap are perturbed together, which needs far fewer residual evaluations per Jacobian than the band stride.
ColouredCheck	Coloured, plus a full band evaluation at every Jacobian for comparison; a diagnostic that reports any mismatch. The comparison only observes (it learns nothing and changes no counter), so the run is exactly the Coloured run and a missed coupling shows up as mismatches until an audit adds it.

How the Coloured mode stays safe:

• The pattern is learned from band evaluations: at three randomised "generic" states near the current one at the start of each step of the first cycle (so couplings that vanish only at special states such as \(u = 0\) or \(y = 0\) still show up), and at the first ten real evaluations of every step in the first cycle, since each step has its own boundary conditions.
• Interior columns share a translation-invariant pattern keyed by variable and cell offset; the three cells at each end keep their own exact patterns, because boundary conditions overwrite those residuals, and also take the interior pattern. A boundary column's own pattern is learned from one column only, and with WENO or van Leer some of its couplings did not show up in the first cycle's learning (for example the first cell's mole fractions on the inlet ghost cell's other mole fractions, through the face rescaling); the interior pattern has seen the same coupling all along the bed. With it, the coloured Jacobian equals the band DQ Jacobian entry for entry for all three advection schemes on the decks test_coloured_jacobian.cpp runs.
• Every 50th Jacobian is a full band evaluation that audits the pattern; any coupling found outside it is added and the columns are recoloured, so a missed dependency corrects itself.
• The Jacobian requested after a Newton convergence failure is always a full band evaluation, so in rare regimes (a boundary whose flow reverses, for example) the solver is never less robust than IDA's own Jacobian.

The benchmark output of bench_solver reports the number of colours, audits and pattern misses for each deck.

Cyclic steady state: detection and acceleration

After the last step of every cycle the solver applies the deck's cyclic steady state rule (CyclicSteadyState::CheckCSS); if it holds, the run stops and the status becomes Cyclic steady state reached at cycle n. Each method has a single-cycle test (below), and the rule holds when that test has passed on Consecutive cycles cycles in a row (default 3) and at least Minimum cycles cycles have run. The check reads only the equation log, so the logs keep at least the cycles it needs.

Method	Single-cycle test
None	never; the run does every cycle
Molar Profile	\(\lVert \mathbf{C}^{\text{curr}} - \mathbf{C}^{\text{prev}}\rVert_2 \le \varepsilon_{abs} + \varepsilon_{rel}\lVert \mathbf{C}^{\text{curr}}\rVert_2\), where \(\mathbf{C}\) is the total-concentration profile over all cells at the end of the first step of the current and the previous cycle
Mass Balance (Total)	\(\lvert n_{in} - n_{out}\rvert / n_{in} \le \varepsilon_{mb}\) over the whole cycle, summing every face of every step by the sign of its net flow, and \(\Delta_{\text{state}} \le \varepsilon_{\text{state}}\)
Mass Balance (Per-Component)	the balance for each component that is fed, taking the worst, and \(\Delta_{\text{state}} \le \varepsilon_{\text{state}}\)

\(\Delta_{\text{state}}\) (CyclicSteadyState::EndOfCycleStateChange) compares the last logged row of the cycle's last step with the previous cycle's: the largest of \(\operatorname{RMS}(x^{c} - x^{c-1}) / s\) over the adsorbed loading of each component (\(s\) the RMS total loading, at least the solver's absolute tolerance on loadings, 10⁻⁸ mol/kg, so the round-off loadings of an inert packing read as zero change), the mole fractions (\(s = 1\)) and the fluid, solid and wall temperatures (\(s\) their RMS). The boundary flux in the balance is reconstructed from the logged face values and integrated over the logged points, so the balance carries a small systematic error of its own; the state test does not, and it does not cross zero the way the balance does on its way to CSS. The Cyclic Steady State page gives the reasoning and the measured stopping cycles.

Aitken (Graves–Morris) acceleration. The approach to CSS is often a slow geometric decay: each cycle closes a fixed fraction of the remaining gap. With CSS acceleration method set to Aitken (Graves-Morris), the solver extrapolates the end-of-cycle state towards the fixed point (AccelerateCyclicSteadyState::Accelerate):

1. From the third cycle on, it collects three successive natural end-of-cycle states \(x_0, x_1, x_2\).
2. It computes one scalar coefficient shared by all fields,

\[ \mu = \frac{\sum_f w_f\thinspace r_{2,f}\cdot d_f}{\sum_f w_f\thinspace d_f\cdot d_f},\qquad r_1 = x_1 - x_0,\thickspace r_2 = x_2 - x_1,\thickspace d = r_2 - r_1,\thickspace w_f = \frac{1}{\overline{x_f^2}} \]

where the weights stop fields of large magnitude (temperatures) from dominating fields of small magnitude (mole fractions), and moves every field to \(x^{*} = x_2 - \mu\thinspace r_2\). Sharing one coefficient keeps the moved fields consistent with one another. 3. Only the slow differential fields are moved: the adsorbed loadings, the solid and wall temperatures and the mole fractions. The total concentration and gas temperature are left alone (they fix the pressure through the gas law and relax within the next step), and the algebraic unknowns are recomputed. 4. It checks the three samples before jumping. For a single geometric mode of decay rate \(\rho\), \(r_2 = \rho\thinspace r_1\) and \(\mu = \rho/(\rho - 1)\), which is exact. The solver - skips the jump if the steps are not shrinking (\(\lVert r_2\rVert_w \ge \lVert r_1\rVert_w\)) or there is no curvature to fit; - skips it if \(\mu\) has the wrong sign for the single-mode rate \(\hat\rho = (r_1\cdot r_2)_w/(r_1\cdot r_1)_w\): a monotone decay (\(\hat\rho > 0\)) needs \(\mu < 0\). A positive \(\mu\) with steps that point the same way (typically a front drifting by the same amount every cycle) would move the state back along its own trajectory. A positive \(\mu\) belongs to an oscillating mode (\(\rho < 0\), \(0 < \mu < 1/2\)); one of 1 or more has no contracting mode behind it and is skipped too, so the limit below applies to negative \(\mu\) only; - limits the jump to twenty natural steps (\(\lvert\mu\rvert \le 20\), the coefficient of a decay rate of 0.95). A slower mode asks for more (\(\mu \approx -28\) at \(\rho = 0.965\), \(-140\) at 0.993); the shortened jump undershoots the fixed point, multiplying the remaining error along that mode by \(1 - 20(1-\rho)/\rho\) (0.27 at \(\rho = 0.965\), 0.86 at 0.993). The limit also bounds what a jump does to faster modes still in the samples. It was chosen on measured decks: with limits up to 20 the KPIs at the CSS stop stayed within the error the stopping rule leaves on a plain run, while a limit of 50 stopped Perez 2019 F (10 cells) 1.6 times further from its CSS than the plain run.

Each skipped window is logged as CSS acceleration skipped (reason), and at the high log level each window's fit as CSS acceleration fit: mu = …, step ratio …. After the jump, negative values are clipped to zero and the mole fractions renormalised to sum to one. 5. The extrapolated state is first integrated for a short probe (the smaller of 2 s and a tenth of the first step) in a throwaway copy of the integrator. Only if that succeeds is it committed to the live run, with a fresh IDACalcIC; otherwise it is discarded and the log says CSS acceleration skipped (re-seed not integrable). 6. The three samples are then discarded, so the next extrapolation uses three fresh cycles. Jumps are therefore at least three cycles apart and the bed relaxes between them.

The gain depends on the deck. For the Perez 2019 VPSA tutorial deck (Perez2019_D1: WENO, 80 cells, feed-loaded start), acceleration brings CSS from cycle 94 to cycle 48, with purity within 1.3×10⁻⁵ and recovery within 5.7×10⁻⁵ of the plain run's. Nguyen 2023 experiment X on 20 cells (CSSAcceleration.SlowGeometricModeReachesCssSooner), whose approach to CSS is one slow mode of rate 0.984, reaches CSS at cycle 84 instead of 349; its purity and recovery differ from the plain run's by about 10⁻⁶, while the plain stop itself is 5×10⁻⁵ and 1.1×10⁻⁴ from the converged values. On the test cycle of CSSAcceleration.ReducesCyclesToCSS* acceleration reaches CSS in 24 to 26 cycles instead of 39 with each advection scheme. The shipped decks with Aitken and Anderson are compared on the Cyclic Steady State page.

Anderson acceleration. With the method set to Anderson, the solver combines several earlier cycle-to-cycle steps instead of fitting one geometric mode (AccelerateCyclicSteadyState::AccelerateAnderson). It moves the same slow fields, in the same weighted norm, through the same projection and the same throwaway-integrator probe as Aitken.

1. It keeps up to five pairs \((x_j, G(x_j))\) of a cycle's start and end state (\(G\) is the cycle-to-cycle map). A pair is kept only when \(x_j\) is known exactly: when the previous cycle ended without a jump, so \(x_j\) is that cycle's natural end state. The cycle after a jump gives no pair, because the probe may have rejected the jump and the solver then continued from the natural state. The older pairs remain valid samples of \(G\) and are kept, so the method is not restarted. Pairs start from the third cycle.
2. With residuals \(f_j = G(x_j) - x_j\), the newest pair \(n\) and the older pairs \(j\), it solves the small least-squares problem $$ \gamma = \arg\min_\gamma \Bigl\lVert f_n - \sum_j \gamma_j (f_n - f_j) \Bigr\rVert_w $$ by column-pivoted QR, dropping the oldest pairs until the columns are independent, and jumps to \(x^{*} = G(x_n) - \sum_j \gamma_j \bigl(G(x_n) - G(x_j)\bigr)\). For a cycle map that is linear with \(k\) error modes, \(k\) independent steps put \(x^{*}\) on the fixed point.
3. It skips the jump when the fitted residual is more than half of \(\lVert f_n\rVert_w\) (little gain, and a late jump delays the stopping rule, which needs consecutive quiet cycles), and limits it to twenty natural steps \(\lVert f_n\rVert_w\), as for Aitken. Jumps are at most every other cycle.

The depth (four earlier steps), the jump spacing and the skip threshold (one half) were chosen on Perez 2019 A and E and Nguyen 2023 X. With a threshold of 0.9, a depth of three or five, or a jump every third cycle, stopped Perez 2019 E at cycle 55 to 64 instead of 48, and Perez 2019 A between cycles 34 and 37 instead of 35. With depth four, thresholds of 0.7 and 0.9 stopped Perez 2019 A and E as late or later than one half, and 0.9 stopped Nguyen 2023 X at cycle 234 instead of 190. The stops and their KPI errors are tabulated on the Cyclic Steady State page.

Earlier versions skipped every window with \(\lvert\mu\rvert > 5\) (a decay rate above 0.83) instead of limiting the jump, and accepted a \(\mu\) of either sign. Once a run settled on a slow mode no window passed again (with WENO or van Leer, acceleration stopped after two jumps on the test cycle), and in a drifting-front phase the accepted positive \(\mu\) moved the state backwards: Perez 2019 D1 with WENO reached CSS at cycle 250 accelerated against 237 plain.

Convergence of the results

A KPI reported at cyclic steady state carries three numerical errors: the grid, the integrator tolerance and the distance from CSS. The verification tests measure each one on a synthetic cycle; this section summarises them and gives the grid study of three example decks.

Grid. Each deck was run to CSS on 10, 20, 40, 80 and 160 cells (python scripts/verification/grid_convergence.py, which writes the full table, with the productivity too; CSS state tolerance tightened to 10⁻⁶, everything else as in the deck). From the three finest grids, the observed order is \(p = \log_2\big((f_{40}-f_{80})/(f_{80}-f_{160})\big)\), Roache's grid convergence index of the 160-cell value is \(\text{GCI} = 1.25\thinspace\lvert(f_{160}-f_{80})/f_{160}\rvert/(2^p-1)\), and the Richardson extrapolation \(f_{160} + (f_{160}-f_{80})/(2^p-1)\) estimates the grid-converged value \(f^{*}\). The last column is the grid at which the error \(\lvert f_N - f^{*}\rvert \approx \lvert f_{160} - f^{*}\rvert\thinspace(160/N)^p\) falls to 0.1 percentage points.

The same analysis is in the engine as adsim/solver/grid_convergence.h, for values on three or more grids at any constant refinement ratio \(r\) (the formulas above with \(\log_r\) and \(r^p\) in place of \(\log_2\) and \(2^p\)). AnalyseGridConvergence(cells, values) returns the observed order (of the three finest grids, and of every consecutive triple), the GCI, the extrapolated value and the cells for a target error (default 10⁻³, 0.1 points), with a status: monotone; oscillatory (the differences change sign) or divergent (a difference does not shrink), both outside the asymptotic range, so no order is reported; grid independent (the two finest values equal); or too few levels. The building blocks ObservedOrder, RichardsonExtrapolation, GridConvergenceIndex and CellsForTargetError are public too. The cycle-KPI grid test (core/test/validation/test_grid_convergence.cpp) uses it; the script above keeps its own Python copy of the formulas until the analysis is exposed to Python.

Deck	Scheme	KPI	10	20	40	80	160	\(p\)	GCI	\(f^{*}\)	Cells for 0.1 pt
Perez2019_A	Upwind	CO₂ recovery	0.5772	0.5874	0.5921	0.5945	0.5957	0.96	0.26 %	0.5969	200
	WENO		0.5872	0.5934	0.5956	0.5965	0.5968	1.49	0.036 %	0.5970	50
	van Leer		0.5897	0.5951	0.5963	0.5967	0.5969	1.40	0.023 %	0.5970	35
	Upwind	CO₂ purity	0.9942	0.9947	0.9949	0.9950	0.9950	1.14	0.004 %	0.9951	10
Ward2024	Upwind	CO₂ recovery	0.8936	0.9054	0.9120	0.9149	0.9165	0.86	0.27 %	0.9185	360
	WENO		0.9040	0.9129	0.9166	0.9180	0.9187	1.14	0.074 %	0.9192	95
	van Leer		0.9097	0.9157	0.9176	0.9185	0.9189	0.96	0.062 %	0.9193	70
	Upwind	CO₂ purity	0.9643	0.9628	0.9618	0.9613	0.9610	0.97	0.033 %	0.9608	40
Illustrative_TSA_13X	Upwind	CO₂ recovery	0.6602	0.6634	0.6649	0.6655	0.6658	1.12	0.05 %	0.6661	50
	WENO		0.6648	0.6657	0.6660	0.6662	0.6662	(3.2)	< 0.001 %	0.6662	20
	van Leer		0.6651	0.6658	0.6661	0.6662	0.6662	–	< 0.001 %	0.6662	15
	Upwind	CO₂ purity	0.8333	0.8343	0.8348	0.8350	0.8351	1.04	0.016 %	0.8352	20

What the study shows:

• Every deck converges monotonically, and upwind at first order (\(p\) = 0.86 to 1.14). WENO and van Leer have errors two to nine times smaller at the same grid, but on Ward2024 their order also falls to about 1 on the finest grids: the two bed-end faces take the upwind value in every scheme (see Advection schemes), so a KPI read at a bed end is formally first order whatever the interior reconstruction; how soon that shows depends on the gradient at the product end. On the synthetic cycle of the tests, whose product end is almost flat, WENO and van Leer converge at \(p \approx 2.5\) on 20 to 80 cells.
• Recovery is the KPI that needs the grid; purity much less so. At the 30 upwind cells these decks used up to engine 0.4, the recovery of Perez2019_A is about 0.6 points and that of Ward2024 about 0.8 points below its grid-converged value (from the fitted error), while the purities are within 0.02 and 0.13 points. This matches the open-source benchmark, in which the reference code's 10-cell recovery rose from 91.6 % to 96.3 % at 100 cells.
• The grid each deck needs. For recovery within 0.1 points: Perez2019_A about 200 cells with upwind or 35 to 50 with WENO / van Leer; Ward2024 about 360 with upwind or 70 to 95 with the higher-order schemes; Illustrative_TSA_13X 50 with upwind, 15 to 20 with WENO / van Leer. For a design study comparing small differences between runs, a higher-order scheme on 40 to 80 cells is far cheaper than upwind on several hundred. Every example deck has since been studied the same way; they now run WENO on the counts listed under Advection schemes (Perez2019_A 80, Ward2024 100, Illustrative_TSA_13X 30). For a deck of your own, run the script before relying on a KPI to a fraction of a point.
• On the TSA deck the higher-order schemes are converged to about 10⁻⁵ by 80 cells, which is the level of the integrator and CSS tolerances of the study: their 80-to-160 differences are noise, so \(p\) there is not meaningful (van Leer's is reported as oscillatory).

Integrator tolerance. IDA's error-per-step control gives a global error proportional to \(\text{reltol}^{k/(k+1)}\) for BDF order \(k\), a slope between 0.5 and 1 on a log-log plot: measured 0.68 on a breakthrough and 0.77 for the bed state of a cycle. Once \(\text{reltol}\thinspace\lvert y\rvert\) falls below the absolute tolerance, the error stops falling and stays at about the absolute tolerance (10⁻⁸ for mole fractions and loadings). The KPIs converge more slowly (slope 0.42 for the recovery of the test cycle), because they integrate the logged boundary flows by the trapezoid rule and the flow transient after a step opens a pressure boundary lasts only the pressure diffusion time of one cell (about \(\mu\thinspace\varepsilon_b\thinspace\Delta z^2/(K P) \approx 2\times10^{-5}\) s), which only IDA's own steps resolve. At the default reltol of 10⁻⁶ the recovery error of the test cycle (with upwind) is 6×10⁻⁵, a quarter of the change from 20 to 40 cells; tighten reltol only when comparing runs that differ by less than that. With the default WENO scheme the grid error is much smaller (a 20 to 40 cell change of 1.8×10⁻⁵ on the same cycle), so there the reltol error, about 5×10⁻⁵, is the larger of the two; both are far below the 0.1-point level to which the example decks are grid-converged.

Maximum step. max_step caps IDA's step; it never relaxes the error test, so results do not depend on it beyond the integrator tolerance: in the tests a cap above IDA's own steps leaves a run bit-for-bit the same, and a 0.1 s cap (twice the steps) moves profiles and KPIs by less than the error of a 100 times looser reltol. Decks set it for robustness, not accuracy: Vysyaraju2024_PSA_exp1 uses 0.5 s because with a 1 s cap IDA stalled at its minimum step in one equalisation step; with 0.5 s (or reltol 10⁻⁷ at 1 s) the run completes and the two agree to 0.01 points in purity and recovery.

Cyclic steady state. Near CSS the end-of-cycle map contracts by a fixed factor \(\rho\) per cycle, so the state change decays geometrically and every KPI converges at the same rate. The error a KPI still carries when the stopping rule fires at a state change \(\Delta \le \varepsilon\) is then \(g\thinspace\Delta\thinspace\rho/(1-\rho)\), with \(g\) the KPI change per unit state change: proportional to the state tolerance and to \(1/(1-\rho)\). On the test cycle with its energy balance on (\(\rho = 0.966\)) the default rule stops at cycle 108 with a recovery error of 5.95×10⁻⁵ against 5.94×10⁻⁵ predicted; isothermal (\(\rho = 0.106\)) each decade of state tolerance costs one cycle and buys a decade of KPI accuracy. Aitken acceleration stops 28 cycles earlier on the same cycle, on KPIs 1.4 to 2.8×10⁻⁵ from the CSS value, inside the error the rule allows a plain run.

The chunked solver for cooperative front ends

The C++ API (and the Python package) can run a whole cycle loop in one call (CoupledSolver::IntegrateCycle). A single-threaded host such as a browser worker cannot do that: a long blocking call could never receive a Cancel message. For such hosts the engine provides CoupledSolverChuncked, which performs the same steps, retries, CSS checks and acceleration with the same integrator code, but in slices: each call to Integrate(max_time_ms) advances the current step until the wall-clock budget is used or the step ends, then returns. A WebAssembly front end can call it with a short budget (e.g. 50 ms) in a loop, yielding between slices so that progress, log messages and cancel requests get through. Results are the same as a native run up to floating-point rounding differences between compilers.