Temperature, pressure and dose rate
This model reports none of these. The first three charts are the canister model’s; use Your own selection below.
Water, relative humidity, oxygen and corrosion
This model reports none of these. The first three charts are the canister model’s; use Your own selection below.
Amounts of the main species
This model reports none of these. The first three charts are the canister model’s; use Your own selection below.
Your own selection
Tick series above to draw them here.
What this page does
It integrates the gas-phase radiolysis and steel-corrosion model of an intact spent-fuel canister that was written in FACSIMILE for SKB (report TR-22-15) and later ported to Python. The reaction system, the constants and the case settings are a text file you can edit in the Model tab; the settings column shows the values from its <SETTINGS> section, and the Scenario list fills them with the cases of the study. Choosing one also says where that case is defined: cases 1–22 are Table 3-1 of SKB TR-22-15, with the section that presents its results; the rest are variants that come with the FACSIMILE model files rather than the report, and say so instead of borrowing a citation. The reference travels with the run into the HDF5 file and the spreadsheet.
Everything runs in your browser. The text is compiled into JavaScript for the right-hand side and for its Jacobian, and the system of stiff ordinary differential equations is solved with a variable-order multistep method of orders 1–5 — the numerical differentiation formulas, or the plain backward differentiation formulas.
The button at the right end of the tabs gives the page the whole window, without the site’s header and footer; the kvot mark at the left of the tabs then goes to the home page, and the switch beside the button changes between light and dark. The page remembers the choice.
The Jacobian
A stiff solver needs the matrix df/dy repeatedly. Instead of estimating it with finite differences (one evaluation of the whole model per species), it is derived from the equations when the text is compiled: every expression is emitted together with its derivative with respect to each species it depends on, so a rate constant that contains the third-body concentration M contributes to the columns of the six species in M, and a switch such as ramp contributes its true slope rather than the artefact a difference across its kink gives. Which entries can be non-zero is known before any number is computed, so the matrix is stored and factorised as a sparse matrix (compressed column storage, Gilbert–Peierls LU with a reverse Cuthill–McKee ordering). Whether the sparse factorisation is used is decided from its measured fill: for this model of 63 species the factor fills in to more than a third of the dense matrix, and a dense LU is then cheaper. The Jacobian tab shows the pattern, and Check Jacobian compares the analytic entries with finite differences.
Two details of the Newton iteration matter for this model. The corrections are computed in variables scaled by max(|y|, atol/rtol), because the concentrations span forty orders of magnitude and an unscaled factorisation gets the trace species wrong. And at least two Newton iterations are taken per step: the relative-humidity switch on the corrosion rate turns on over one part in 107 of the humidity, and a Jacobian formed one step earlier can damp the first correction by a factor 106 while the residual is far from zero.
Differences from the FACSIMILE original
- Rate constants are evaluated at the current temperature; FACSIMILE computed them once at the initial temperature.
- Constants carry more digits (Avogadro's number, the gas constants, 365.25 days per year).
- Liquid water is a setting,
H2OPAIRin the case settings. At 0, the default, liquid water is not a separate species: the stateH2Oholds all water and the part above the saturation concentrationH2OEQis liquid (the@H2O = min(H2O, H2OEQ)line), which asserts saturation directly. At 1 you get FACSIMILE's own treatment instead: a separate liquid speciesH2OLIQand the fast pair it used, condensation first order in the vapour and evaporation zeroth order at1E6·H2OEQ, switched off once the liquid is gone. The two holdH2OatH2OEQby different means and agree when the answer has converged — on the 500-year 13g case both give 65.94 g of water atrtol10−8 — but not before it has: at the tolerance this page starts with they read 65.93 g and 63.00 g, for the reason given under the solver comparison. Choosing a scenario putsH2OPAIRback to 0, because that is the model the published cases were run with. - A negative concentration is read as zero in the rate laws (the same as the Python port), and the solver projects a species that a step has left negative back onto zero.
- The gas density uses FACSIMILE's molar masses for NO2 (44) and HNO2 (45) so that the ppm outputs match the reports.
Against the Python port the engine agrees to 0.1 % or better for the amounts of O2, H2, NH3, HNO3, HNO2 and H2O2 while oxygen is consumed and the gas chemistry develops. In the phase where corrosion is throttled by the humidity switch the two integrators follow the same water curve to within a few per cent, and both agree with the FACSIMILE run of case 13g on the water left after 500 years. A 500-year case takes one or two seconds here; the Python port needed minutes because it differenced its Jacobian.
Which settings apply
The solver settings on the left are not all read by every method, and the ones that are not are hidden with a line underneath saying which. Maximum and minimum order mean nothing to a one-step method such as Rodas5P or KenCarp4 — only the multistep formulas have an order to cap. A Rosenbrock method has no nonlinear iteration to give a Newton tolerance to, and its Jacobian is part of the method rather than an optimisation, so it is re-formed every step and there is no age to set. RadauIIA5 measures its error against Hairer's own transformed tolerances, in a norm of its own. An explicit method — Dormand–Prince, Tsit5, Vern7 — forms no matrix and reads no Jacobian, and FBDF by GMRES forms no matrix either. Norm control is the NDF's alone, and Auto reads what the NDF reads, since the NDF is what it hands the stiff part of a run to.
Three settings deserve a word because they are not all-or-nothing.
- Reading negatives as zero is not a solver setting at all. It changes the generated rate laws, so it is built into the derivative and the Jacobian that every solver is handed, and it applies to all of them identically.
- Keeping every species non-negative is three things at once, and the methods do different amounts of it. NDF, with its BDF formulas or without, and Rosenbrock 2-3 do all three: they damp the derivative so that a species at zero cannot be pushed below it, fold the size of any violation into the error test, and project the accepted step back. Dormand–Prince does the last two, since an explicit method's stages would straddle the kink a damped derivative makes. The Julia ports project only, and Auto does that while it is explicit and all three once it is on the NDF. The panel says which you are getting.
- The Newton tolerance κ is how closely each stage is solved, as a share of one unit of the integration tolerance. It is not a cosmetic knob: at 1/100, which is what SciML uses, the left-over iteration error in an ESDIRK's stages is enough to corrupt its embedded error estimate, and the step-size controller then acts on a corrupted number. The default here is 1/1000, which on the standard stiff test problems is both more accurate and cheaper.
What each method reads is declared in the code that hands the settings over, which is the only place the answer stays honest.
The other methods
The method menu is Kompartment's, under each method's own name. Three of them come with the solver core this page shares with Kompartment and the reactive transport page: Auto, Rosenbrock 2-3 and Dormand–Prince 4-5. The rest are ported into JavaScript from DifferentialEquations.jl: its default algorithm, DefaultODEAlgorithm, and FBDF, FBDF by GMRES, QNDF, Rodas5P, Rosenbrock23, RadauIIA5, KenCarp4, TRBDF2, Tsit5 and Vern7. They are a package of their own — the sources are under resources/js/ode/julia/ with their own tests and README — and they are plain JavaScript: nothing is downloaded and they work offline exactly as the built-in NDF does.
Auto starts on an explicit method — Tsit5, or Vern7 below a relative tolerance of 10−5 — and watches every step for stiffness, as DifferentialEquations.jl's default does. Once it finds it, the NDF takes the rest of the run, and every segment after the next event as well. A model that never turns stiff stays explicit throughout. The explicit start works at a tenth of the tolerances and lands on every corner of the tables the model reads at the clock — here the temperature profile — because its error test is on the ends of its steps and a long step across a change of slope can be right at both ends and wrong between; for the same reason it reads the model's output times off its own interpolant rather than off a cubic between its steps. A model with algebraic variables is the NDF's from the start. DefaultODEAlgorithm is the same start as DifferentialEquations.jl has it, switching to a stiff method of its own — Rosenbrock23 up to 50 species, FBDF up to 500, FBDF by GMRES above — and back when the stiffness passes; it is there to compare with Julia, and Kompartment measured Auto the faster of the two at every size of model.
Which does what. Rosenbrock 2-3 is Shampine and Reichelt's Rosenbrock pair, one Jacobian and one factorisation a step and no Newton iteration, and Rosenbrock23 is the same method as DifferentialEquations.jl runs it. Dormand–Prince 4-5, Tsit5 and Vern7 are explicit Runge–Kutta pairs of orders 5, 5 and 7: no Jacobian at all, cheap on a smooth model that is not stiff, and on a stiff one held at the edge of their stability. FBDF by GMRES solves FBDF's Newton iterations without forming a matrix, which is what DifferentialEquations.jl's default does above 500 states. QNDF is the same method as this page's own NDF — the same numerical differentiation formulas, the same κ coefficients and the same backward-difference machinery — written by other people, so it is the most direct check there is on the built-in solver. Both have a BDF formulas switch under Advanced settings, which sets κ to zero throughout, and QNDF then runs as QBDF. That is what turns a numerical differentiation formula back into a plain backward differentiation formula, so running either with the switch on and off shows exactly what those terms are worth on a given problem. FBDF is a close relative: a multistep formula of variable order that reuses one matrix factorisation across many steps, which is what makes it cheap on a system this size. RadauIIA5 is fully implicit and suffers the least from stiffness, so it is the one to believe when two others disagree. KenCarp4 and TRBDF2 are diagonally implicit and cheap at moderate and loose tolerances. Rosenbrock methods like Rodas5P take one linear solve per stage and no nonlinear iteration at all, so there is nothing that can fail to converge.
On this model, in particular. The 500-year case 13g at the page's defaults, timed in Node on the development machine:
| solver | steps | rejected | seconds | water left (g) |
|---|---|---|---|---|
| NDF | 5 665 | 1 600 | 0.4 | 42.18 |
| BDF | 5 519 | 1 483 | 0.4 | 65.20 |
| Auto | 4 990 | 1 295 | 0.4 | 59.45 |
| QNDF | 14 358 | 22 002 | 10 | 65.94 |
| KenCarp4 | 12 701 | 30 846 | 15 | 65.94 |
| TRBDF2 | 31 335 | 65 838 | 31 | 65.94 |
| DefaultODEAlgorithm | 104 134 | 130 442 | 70 | 65.93 |
| FBDF | 109 545 | 139 832 | 73 | 65.94 |
Auto took four Tsit5 steps before it found the model stiff, and DefaultODEAlgorithm fourteen before it went to FBDF. The rest do not finish the case at all. Rosenbrock 2-3 uses up its two million steps 140 years in, Rodas5P and Rosenbrock23 theirs about 13 years in, after ten minutes, and RadauIIA5 had not finished after fifteen. FBDF by GMRES, Dormand–Prince, Tsit5 and Vern7 do not get past the first microsecond: the explicit methods, and GMRES with nothing but a diagonal scaling to help it, are held there by the stiffness of the fast ions.
Why the rejections. The absolute tolerance here is 10-30 and the trace species sit forty orders of magnitude below the main ones, so the error test is asking for the impossible on components that hold nothing. This page's NDF gets through by accepting, up to five times in a row, a step at its smallest allowed size that has failed the error test anyway. That is a departure from the published method, not a feature of it: reaching the same point, a variable-order multistep code raises a tolerance-not-met error and returns the partial solution. Without the departure, three of the presets here stop early; with it, they finish, and the footer says in bold how many such steps were taken. The ported methods are handed the same escape, which is not theirs in Julia. Loosening the absolute tolerance, or choosing the root-mean-square norm instead of the maximum, makes all of them far cheaper and is worth trying before concluding a method is unsuitable.
Read the table for cost, not for accuracy. The water left after 500 years is not converged at the relative tolerance this page starts with. Run the 13g case at a sequence of tolerances and it goes:
| rtol | 10-4 | 10-5 | 10-6 | 10-7 | 10-8 | 10-9 |
|---|---|---|---|---|---|---|
| NDF | 46.53 | 42.18 | 40.26 | 49.71 | 65.94 | 65.94 |
| BDF | 19.76 | 65.20 | 65.93 | 58.74 | 65.94 | 65.94 |
From 10-8 downwards both agree on 65.94 g and stay there, and the ported methods that finish the case land on it already at the defaults. Above it the NDF and the BDF scatter between 20 and 66 g, and not even monotonically. Water at 500 years sits on the far side of the corrosion switch and the relative humidity that drives it; small differences early decide which side of that threshold the run comes down on, and an agreement at a loose tolerance is partly luck. If the number matters, tighten the relative tolerance to 10-8 and check it stops moving.
How much of a run is kept
The charts, the table and the exported files hold the solver's own steps, up to twenty thousand of them. Past that the run is thinned as it goes — every second point, then every fourth — so that a solver taking hundreds of thousands of steps cannot fill the machine's memory. The two ends of the run are always kept, and the footer says when it happened: one in 32 beside the step count.
The model tab
Colour the syntax paints comments, section headings, numbers, function names and the words that mean something where they stand — kf and rate in a reaction, down and stop in an event, log and lin in a time list. It is a second copy of the text drawn underneath the box you are typing in, which is the only way a plain text area can be coloured while keeping its own undo, selection and caret; turn it off for a very large model, or if the colours are in the way. The choice is remembered between visits.
The panel on the left
Drag the edge between the panel and this area to give either more of the page; double-click it to put it back. Each heading folds away what is under it. The width and what is folded are remembered between visits. Run, Stop and Check Jacobian sit at the foot of the panel and stay there whatever is scrolled past them, with the status line beside them.
The model file
A section starts with a line <NAME>. # starts a comment; FACSIMILE's * comment lines and trailing ; are accepted. Names are case-sensitive; function names are not.
| Section | Lines |
|---|---|
<SETTINGS> | NAME = value # label (unit). A value may be an expression of earlier settings and constants, or the name of a table (a profile). These lines are what the settings column on the left shows: typing in a field rewrites its line here, and editing the line here shows up in the field. There is one copy of a setting and it is the line, so the file you save is the case you ran. |
<CONSTANTS> | NAME = expression, evaluated once. |
<SPECIES> | Names, in the order you want them listed. Species that appear only in reactions are appended. |
<TABLE NAME> | Two numbers per line, x y, with increasing x; used as interp(NAME, x), linear and clamped at the ends. |
<THERMO> | NAME a1 … a7, CHEMKIN coefficients; defines DGNAME, the Gibbs energy in cal/mol at the temperature T. |
<INITIAL> | Evaluated at t = 0 in order. A species name sets its initial concentration (mol/cm³); another name becomes a constant of the run. Equations that depend on time only (T, H2OEQ, …) may be used. |
<EQUATIONS> | NAME = expression, evaluated at every derivative evaluation in order; may use t (seconds), the species, the constants and earlier equations. @X = expression replaces the species X in everything below it (the effective concentration), while the state keeps the total. |
<REACTIONS> | A + 2 B = C, kf = k (mass action) with optionally kb = k or keq = K (then kb = kf/keq); rf = r and rb = r give absolute rates instead. FACSIMILE's %kf%kb : A + B = C ; and = r : A = ; are accepted too. Reactants or products may be empty. Adding rate = NAME, or writing the name before the first % in FACSIMILE's own form, reports the net rate of that reaction as a quantity of its own. |
<EVENTS> | expression, NAME = value, …: when the expression crosses zero from below, the named settings or constants are changed, the run constants are recomputed and the integration restarts there. down or both beside the assignments chooses which way the crossing is to be taken, and stop ends the run there instead of carrying on — FACSIMILE's % on a WHEN clause. An event must do one or the other. Writing the trigger as expression = v1 v2 v3 fires it wherever the expression passes each of those values, which is FACSIMILE's valuelist, and the value forms of <TIMES> work here too. once makes it fire one time only, which is the whole difference between FACSIMILE's WHEN and its WHENEVER. |
<ALGEBRAIC> | NAME : expression: the solver varies NAME to hold the expression at zero. NAME is a state variable like a species, but has no differential equation — its value is whatever satisfies the constraint — so nothing may declare it in <SPECIES> or change it in a reaction, though a reaction may use it as a catalyst or a third body. Give it a starting value in <INITIAL>; the solver moves it onto its constraint before the run and reports how far it had to. |
<TIMES [unit]> | Times at which the run is also to be reported, in s (the default), min, h, d or y. A line is a list of values, first + step * count, or from .. to log|lin count. The Table tab then offers at the model's output times, and the CSV and the Excel file follow it; the values there are the cubic through the two accepted steps on either side and their derivatives, not the solver's own interpolant, so that every method on the menu gives the same grid. |
<OUTPUTS> | NAME = expression, for the table and the charts; may use earlier outputs. state(X) is the raw state of a species (all water, for H2O), and deriv(X) is its time derivative — what the FACSIMILE model carried in five dummy variables to report its apparent G-value. deriv() belongs here and nowhere else: a derivative inside a rate law would not be a differential equation. |
Expressions use + - * / ** (also ^ and FACSIMILE's @ for powers), parentheses, numbers such as 1.5E-12 or 2.0D+14, and the functions exp log ln log10 sqrt abs min max pow ramp step stepf sign sin cos tan atan tanh amod interp state deriv. ramp(x) is max(x, 0) (FACSIMILE's RAMP) and step(x) is 1 for positive x. stepf, sign, amod and the trigonometric functions are FACSIMILE's, spelled as it spells them (artan is accepted for atan); the one difference to watch is at zero, where FACSIMILE's stepf is 1 and the step this page already had is 0. There is no random-number function: a derivative that changed from one evaluation to the next would not define a differential equation for any solver to integrate.
# the anoxic corrosion of the model, with its net rate named 4 H2O = 4 H2, rf = c2/3*(1 - f1)*f2, rate = RANOX # 3Fe + 4H2O = Fe3O4 + 4H2 (mol cm-3 s-1) # a reversible reaction with an equilibrium constant from Gibbs energies 2 OH = H2O + O, kf = 1.5E+09*T**1.14*exp(-50/T), keq = exp(-(-2*DGOH+DGH2O+DGO)/RT) # FACSIMILE's own form, where the name goes before the first % RAB%K1%K2 : A + B = C
A named rate is reported beside the outputs, on the fourth chart, in the table, in the CSV, in the Excel file and in /Rates of the HDF5 file. It is the net rate of that one reaction — forward less backward, or the absolute rate where one was given — so a species' derivative is the sum of the rates that reach it, each times its stoichiometric coefficient. That is what the built-in model's two corrosion fluxes now are:
3 O2 = , rf = c1/4*f1*f2, rate = ROXID # three O2 consumed per turn dDUMO2 = 3*ROXID # in <OUTPUTS> GHNO3 = (deriv(HNO3) + deriv(HNO2) + deriv(NO2) + deriv(NO3) + deriv(NO))/DOSER
Algebraic variables
A quantity given by an equation rather than by a rate of change goes in <ALGEBRAIC>, and the model becomes a differential-algebraic system: M y′ = f, with a zero on the diagonal of the mass matrix wherever a variable has a constraint instead of a derivative. Nothing about the method changes. The backward differentiation formulas were extended to these systems by Gear, they are what DASSL and SUNDIALS IDA are built on, and the formulas this page uses solve them in their original home too; substituting the difference quotient for y′ leaves a system whose iteration matrix is M − hγJ, which is non-singular whenever the constraints determine the algebraic variables. A model with no <ALGEBRAIC> section hands over a mass matrix of ones and is solved by exactly the same arithmetic as before.
# Robertson's problem, the usual index-1 test, in this form <SPECIES> Y1 Y2 <ALGEBRAIC> Y3 : Y1 + Y2 + Y3 - 1 <REACTIONS> Y1 = Y2, kf = 0.04 Y2 + Y3 = Y1 + Y3, kf = 1.0E4 # Y3 is a catalyst here, which is allowed Y2 = , rf = 3.0E7*Y2*Y2
Three things are worth knowing.
- The start is solved for you. A differential-algebraic system has to begin from a state that satisfies its own constraints, and the value you write in
<INITIAL>generally will not. FACSIMILE asked the author for one close enough for its iteration to work; what that gives when the guess is off is not a warning but a stall. Take the circle of the FACSIMILE User Guide, y₁′ = y₂ with y₁² + y₂² = 1, started at its own documented y₂ = −0.9 where the constraint wants −0.866: the first step moves the algebraic variable by 0.034, the error test measures that jump and rejects the step, and halving the step does not shrink the jump, because the constraint has to hold at the new point whatever the step is. The step size collapses and the run never leaves zero. So the constraints are solved here first, by Newton on the algebraic variables with the differential ones held where you put them, and the footer says how far they moved. Your starting value is a first guess, and only has to be close enough for that Newton to find the right root. - Only the NDF, with its BDF formulas or without — and Auto, which hands such a model to the NDF from the start. The other methods take no mass matrix and would read a constraint residual as a rate of change, so they are refused rather than allowed to produce a wrong answer quietly.
- Index 1. Each constraint must determine its own variable. If it does not, the iteration matrix is singular and the page says which variable and why, rather than reporting an arithmetic failure.
A worked example is in the repository at resources/data/facsimile-langmuir.fac: a gas consumed at a fixed rate with a surface whose coverage is in equilibrium with it at every instant. Open it with Open…, or drop it on the page. It uses all of this at once — an algebraic variable, a named rate, an output grid, an event value list and an event that stops the run — and everything in it has an exact answer to check against.
There is an example of the trade-off in the built-in model already. Setting H2OPAIR to 1 holds water at saturation with FACSIMILE's evaporation and condensation pair at a rate constant of a million, which is the usual way of avoiding a constraint; the default substitutes the saturation expression at compile time instead. An algebraic variable is the third way, and imposes it exactly.
Events at a list of values
An event trigger written as expression = v1 v2 v3 fires wherever the expression passes each of those values, and takes the same forms as <TIMES>, so TIMY = 0 + 50 * 10 is every fifty years for five hundred. Adding once makes an event fire one time only. Those two together are FACSIMILE's WHEN and WHENEVER: the same command, differing only in whether the clause survives its first firing. Restarting the integration at an event, which FACSIMILE asked for with RESTART, is what this page always does.
One limitation is worth stating, because it is shared with FACSIMILE and with most solvers and is invisible when it bites. A crossing is found from the sign of the expression at the two ends of a step, and located by searching the solver's own interpolant between them. When the trigger moves faster than the solution does the step is long while the expression is not, and two things go wrong: an excursion that goes out and comes back inside one step has the same sign at both ends and is not seen at all, and a crossing that is seen can be located on an interpolant too coarse to place it, so the reported time is wrong. On sin(t) against 0.9, left to itself, the second crossing comes out at 3.01 where it belongs at 2.02. The remedy for both is the maximum step size under the advanced solver settings: capping it makes the solver look more often and interpolate over less.
Units
Time in seconds, concentrations in mol/cm³, temperature in K, energies in cal/mol; bimolecular rate constants in cm³ mol−1 s−1 (the literature values in cm³ molecule−1 s−1 are multiplied by NA). Dose rates are given in Gy/h and converted to the 100 eV cm−3 s−1 unit the G-values need.
Solver settings
- Tolerances. The error of each species is compared with
rtol·max(|y|, atol/rtol). The Python port usedatol= 10−100 and below, i.e. purely relative control down to nothing; 10−30 mol/cm³ leaves species below that unresolved, which is where the physics stops mattering. - Per-species absolute tolerance. One
SPECIES VALUEa line, for the species that should not be judged against the single number above; everything unlisted keeps it. Every solver on the menu reads it — the error test has always weighed each species separately, and this is the only part of that weight which was one number for all sixty-three. Two things are worth knowing before reaching for it. It does nothing for a species that is never small. The weight ismax(|y|, atol/rtol), so atrtol10−5 andatol10−30 anything above 10−25 is judged relatively and itsatolnever enters: setting water's to 10−40 or to 10−12 gives the same 5888 steps, to the step. And relaxing the fast trace species — the obvious thing to do — backfires here. Taking the five ionsE,HP,O2M,H3OPandNOPfrom 10−30 to 10−20 costs 71 335 steps against 5888, with the Jacobian re-formed 24 234 times against 698. Raising a tolerance past a species' own magnitude makes it invisible to the error test, the step controller then asks for steps the Newton iteration cannot solve, and the run spends itself failing to converge and cutting back. On this model the error test on the fast ions is what was keeping the step inside the Newton's reach. The cost climbs gently to about 10−22 and then runs away. - Let the absolute tolerance follow the solution upwards. After every accepted step, each species' absolute tolerance is raised to
rtol × |y|if that is larger, and never lowered again. Each component is then judged against the largest it has ever been, rather than against a floor chosen before the run started. The case for it on a model like this one is direct: a radical that rose to 10-5 and has since decayed to 10-40 is otherwise still being held to an absolute tolerance of 10-30, thirty-five orders below anything it ever was, and the step size pays for it. The case against is that the tolerance only ever loosens, so a species that genuinely needs absolute accuracy after having once been large will not get it. It is off by default and it is worth trying, because what it does depends on the method: on the 500-year case it takes FBDF from 4943 steps to 2897 and gives the same answer, leavesNDFslightly worse at 6451 against 5888, and doubles QNDF's work. Across the built-in scenariosNDFtakes fewer steps in total with it on, but one of them stops early. - Accept failing steps at the floor. What to do when a step at the smallest size the clock can represent fails the error test anyway. There is nothing left to try: a shorter step does not exist. Set to 0 the run stops there and hands back what it had, which is what the published method does — it raises a tolerance-not-met error and returns — and what every solver in DifferentialEquations.jl does. Above 0, that many such steps in a row are accepted instead, and the footer reports how many were taken. The case for accepting is that the error test is not always asking a sensible question: after a restart from an interpolated state, a species with a lifetime of femtoseconds is off its steady state by more than the tolerance, and no step size mends that — the implicit step itself does. The case against is that you are accepting a step known to be inaccurate. The default here is 5, because at 0 three of the built-in scenarios stop early; across all thirty-nine, the escape is used three times in total.
- Error norm. How the error of one step is turned from a number per species into a single number to compare with the tolerance. max takes the largest, after each species is scaled by its own tolerance, so the worst-resolved one decides the step; that is what the numerical differentiation formulas were published with and what this page defaults to. rms takes the root mean square instead, which is what CVODE and SUNDIALS do. The difference is not small here: with one species out by exactly one unit of tolerance and the other sixty-two perfect, max reads 1 and rejects the step while rms reads 0.126 and accepts it — a factor of √63, about 8. So rms gives longer steps and a looser answer, and on a model whose trace species sit forty orders of magnitude below the main ones, it is the trace species that the choice is about.
- First step. The step to try first, at the start of the run and again after every event, where the run is restarted; 0, the default, lets the solver choose it from the derivative there. A step that is too long is rejected and cut like any other, so a guess costs steps rather than accuracy. Auto hands it to its explicit start, and the NDF chooses its own once the run has turned stiff.
- Norm control. The NDF's error test weighs each species against its own size; with norm control it weighs the whole vector at once, the 2-norm of the error over the 2-norm of the solution, as MATLAB's NormControl does. The error norm then does not arise, and of the absolute tolerances only the first is read. On a model whose species span forty orders of magnitude that leaves every trace species uncontrolled, which is why it is off; it is there to match a result worked out that way.
- Method. NDF solves a 500-year case here in under a second, and is the one to use; Auto costs the same, since this model is stiff from the first second and it hands the run to the NDF almost at once. BDF is the same integrator with the κ terms switched off: a little more stable and a little less accurate per step, which on this model costs about 6 % more steps for the same answer, and is the honest control on what those terms buy. It gets through thirty-eight of the thirty-nine presets; the one it misses stops at the moment the water supply is cut off, where restarting from the interpolated state leaves a trace species further from its steady state than a plain BDF will tolerate. Letting the absolute tolerance follow the solution gets it through that one, in fewer steps than NDF takes. QNDF is the cheapest independent check among the ported methods, at about twenty-five times the NDF's time. Any run reports its progress in simulated years as it goes, and Stop ends it at once.
- Iteration matrix. auto measures the fill of one sparse factorisation, then tries a sparse LU that keeps its pivots and uses it wherever it costs no more than the sparse or dense LU would. On this model it needs about 6,000 multiply-adds per factorisation, against 85,000 for the dense LU. It chooses its pivots from the values once, among entries at least a tenth of the largest in their column, in an order that keeps the factor small. Then it repeats the same eliminations on every later matrix, checking each pivot: it must be at least a hundredth of the largest entry below it. A pivot that fails is chosen again from that point on, and a matrix it cannot factor goes to the dense or sparse LU with partial pivoting. On this model, that makes each step about 30 % cheaper, and it is as accurate: measured against a run at a tolerance of 1e-9, the two give the same errors. The individual steps are not the same, because on this model a change at the level of round-off, from either factorisation, moves the step sequence by several hundred steps either way.
- Jacobian. finite differences differences the model through the analytic pattern, one evaluation per group of columns that share no row (45 groups for 63 species), for comparison.
- Non-negative. Three things at once: the derivative of a species at or below zero may not take it lower, a step that leaves a species negative by more than the tolerance is rejected, and what is left below zero is projected back. On by default: with it every scenario of the study integrates in a second or two, whereas a trace species that has drifted a little below zero is exactly where the Newton iteration of a stiff solver comes to grief (the derivative of the clamped rate law is zero on that side). The Python port ran without it and paid in tolerance and time.
Files
Save as .fac writes the model text; Open… reads one back, and so does dropping the file anywhere on the page. A dropped file is taken as a model only when it ends in .fac, .txt or .in, and only one at a time, so that something dropped by mistake does not replace the text in the editor. Download Excel writes SETTINGS, DATA (the outputs) and STATES (the species) sheets in the layout of the Python port; CSV writes the table's columns at every solver step. Where the model has a <TIMES> section and the Table tab is set to those times, all three write that grid instead, and the buttons say so.
Sending a run to the HDF5 Browser
View in HDF5 Browser on the Charts tab writes the run as an HDF5 file and opens it in the HDF5 Browser in a second tab. The bytes are handed straight across between the two tabs, so nothing is written to disk and nothing is uploaded. Download HDF5 on the Table tab writes the same file, for h5py or anything else that reads HDF5.
/Settings and /Constants hold the case and what it works out to as a dataset each, carrying its own unit and description — the same two attributes every other dataset in the file has, read the same way from the comment beside the value in the model text. They were attributes on the group to begin with, which read as an empty group in a browser that lists children; then both, which made the group itself a wall of sixty name–value pairs. One dataset each is what a reader can actually open.
The file carries the run and the case that produced it:
/time | the integration clock, in hours, which everything is plotted against |
/Results | one series per <OUTPUTS> line, with the unit read off its comment |
/Equations | one per <EQUATIONS> line: the rate constants, the Gibbs energies, the corrosion switches |
/Species | one per species, in mol/cm³ |
/Settings, /Constants | the case, and what it works out to, as attributes |
/Model | the model text and the reactions, so the file says what produced it |
/Events | when an event fired and what it changed, where the case has one |
Selecting one of the three groups draws all of its series together, which is what the IndexLists attribute on it is for. The writer is a small one (resources/js/kvot-hdf5-write.js), so the page does not have to fetch 4.7 MB of WebAssembly in order to write a few hundred kilobytes; what it writes is read back by h5py in the tests.
Sources: the FACSIMILE model files skbcanister11b/13g/16x.fac and the Python port skbcanister.py of the SKB TR-22-15 study; the built-in solver implements the numerical differentiation formulas of Shampine and Reichelt; thermodynamic data from the CHEMKIN data base (SAND87-8215B).