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 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 10⁷ of the humidity, and a Jacobian formed one step earlier can damp the first correction by a factor 10⁶ 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, H2OPAIR in the case settings. At 0, the default, liquid water is not a separate species: the state H2O holds all water and the part above the saturation concentration H2OEQ is liquid (the @H2O = min(H2O, H2OEQ) line), which asserts saturation directly. At 1 you get FACSIMILE's own treatment instead: a separate liquid species H2OLIQ and the fast pair it used, condensation first order in the vapour and evaporation zeroth order at 1E6·H2OEQ, switched off once the liquid is gone. The two hold H2O at H2OEQ by different means and agree when the answer has converged — on the 500-year 13g case both give 65.94 g of water at rtol 10⁻⁸ — 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 puts H2OPAIR back 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 NO₂ (44) and HNO₂ (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 O₂, H₂, NH₃, HNO₃, HNO₂ and H₂O₂ 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.

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 two families do different amounts of it. NDF and BDF do all three: it damps the derivative so that a species at zero cannot be pushed below it, folds the size of any violation into the error test, and projects the accepted step back. The Julia ports project only. 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 Julia ports

Seven more solvers, ported into JavaScript from DifferentialEquations.jl: FBDF, QNDF, QBDF, Rodas5P, RadauIIA5, KenCarp4 and TRBDF2. 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 pair does.

Which does what. 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. QBDF stands in the same relation to BDF: κ set to zero throughout, which is what turns a numerical differentiation formula back into a plain backward differentiation formula, so running the pair either way 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. Every one of them that finishes agrees with NDF on the water left after 500 years — 65.935 g against 65.933 g — so the physics is not in question. What differs is the cost, and by a lot:

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 — across all thirty-nine presets the total is three. None of the ported methods has the escape, because none of them has it in Julia either. FBDF is the one that does not need it. Rodas5P and RadauIIA5 need more than that and do not get through the case at 10^-30 at all; raise the absolute tolerance to 10^-20 and Rodas5P solves it in a few hundred steps. 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:

From 10^-7 downwards every method agrees on 65.94 g and stays there. Above it they scatter between 36 and 66 g, and not even monotonically — NDF is closer at 10^-5 than at 10^-6. 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. So the 65.933 g against 65.935 g above says the ported solvers and this one take the same path at the page's defaults, not that either figure is right to four significant digits. If the number matters, tighten the relative tolerance to 10^-7 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 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.

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 interp state. ramp(x) is max(x, 0) (FACSIMILE's RAMP) and step(x) is 1 for positive x.

# the anoxic corrosion of the model
4 H2O = 4 H2, rf = c2/3*(1 - f1)*f2      # 3Fe + 4H2O = Fe3O4 + 4H2
# 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)

Units

Time in seconds, concentrations in mol/cm³, temperature in K, energies in cal/mol; bimolecular rate constants in cm³ mol⁻¹ s⁻¹ (the literature values in cm³ molecule⁻¹ s⁻¹ are multiplied by NA). Dose rates are given in Gy/h and converted to the 100 eV cm⁻³ s⁻¹ unit the G-values need.

Solver settings

• Tolerances. The error of each species is compared with rtol·max(|y|, atol/rtol). The Python port used atol = 10⁻¹⁰⁰ and below, i.e. purely relative control down to nothing; 10⁻³⁰ mol/cm³ leaves species below that unresolved, which is where the physics stops mattering.
• Per-species absolute tolerance. One SPECIES VALUE a 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 is max(|y|, atol/rtol), so at rtol 10⁻⁵ and atol 10⁻³⁰ anything above 10⁻²⁵ is judged relatively and its atol never enters: setting water's to 10⁻⁴⁰ or to 10⁻¹² gives the same 5888 steps, to the step. And relaxing the fast trace species — the obvious thing to do — backfires here. Taking the five ions E, HP, O2M, H3OP and NOP from 10⁻³⁰ to 10⁻²⁰ 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⁻²² 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, leaves NDF slightly worse at 6451 against 5888, and doubles QNDF's work. Across the built-in scenarios NDF takes 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.
• Method. NDF solves a 500-year case here in a second or two, and is the one to use. 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. The Julia port of FBDF is about as quick and is an independent check that costs nothing. 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 and keeps the sparse path only when it is cheaper than a dense LU.
• 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. 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.

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:

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, from the Ecolego web port), 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 follows the numerical differentiation formulas of Shampine and Reichelt, as ported for ecolego-js; thermodynamic data from the CHEMKIN data base (SAND87-8215B).