Numerical uncertainty

A number computed by a simulation (a growth rate, a mean, a value at a point) depends on the resolution and on the time step. cn.uncertainty reports it with an error estimate from studies that you declare: three resolutions and three time steps, nothing chosen for you.

A value with its error

import numpy as np

import conatus as cn

def u_at_origin(N, dt):
    """u(0, T) of the heat equation u_t = nu u_xx, u(x, 0) = exp(cos x)."""
    m = cn.Model(cn.Box(x=cn.periodic(2 * np.pi)), name="heat", nu=0.5)
    m.scalar("u")
    m.equations("dt(u) - nu*dx(dx(u)) = 0")
    sim = m.simulate(N=N, scheme="RK222", dt=dt)
    sim.initial(u=lambda x: np.exp(np.cos(x)))
    sim.run(until=0.4)
    return float(sim.grid("u")[0])                   # x = 0 is the first grid point

rep = cn.uncertainty(u_at_origin, N=[8, 12, 16], dt=[0.04, 0.02, 0.01], scheme_order=2,
                     name="u(0, T)")
lo, hi = rep.interval
print(f"u(0, T) = {rep.value:.8f} +- {rep.error:.1e}   trusted: {rep.trusted}")
print(f"observed order in time: {rep.temporal.order:.2f} (scheme order 2)")

u(0, T) = 2.32103549 +- 3.8e-06   trusted: True
observed order in time: 2.01 (scheme order 2)
exact value 2.32103854: inside the interval: True

The last line is the check made after the example: the exact value of this problem (a series of Bessel functions) lies inside the reported interval.

What is computed

quantity(N=..., dt=...) is run five times: the three time steps at the finest resolution, and the three resolutions at the finest time step (the finest run is shared).

• Time: the observed order of convergence \(p\) from the three values; the Richardson value is shown beside the result, and the error estimate is \(1.25\,|q_R - q_3|\), with the safety factor of Roache’s grid-convergence index (Roache 1994). With scheme_order, an observed order more than 0.3 away from it is flagged: the time steps are not yet in the asymptotic range.

• Space: spectral convergence is geometric, \(q(N) = q_\infty + C e^{-aN}\) at the resolutions actually given; the ratio \(\rho\) of successive differences fixes \(a\), and the estimate is \(1.25\,|d_2|/(e^{a h_2}-1)\) with \(h_2\) the last increment of \(N\) (equal increments: \(1.25\,|d_2|\,\rho/(1-\rho)\)). Differences that do not decrease geometrically, and tuple resolutions refined in different proportions, are flagged.

• Result: the value at the finest resolution and time step (the Richardson value is never substituted for it), the sum of both estimates, and the interval. rep.trusted is false as soon as a study is flagged (also without scheme_order: the asymptotic range is then unchecked), and rep.explain() prints every number and assumption.

The time steps must have a constant ratio and the resolutions must increase, the finest last; anything else is refused.

The mean of a statistically steady run

cn.time_average(series, dt=..., discard=k) gives the mean of a correlated time series after its first k samples (the transient is yours to judge: discard has no default), with its standard error from the integrated autocorrelation time. A batch-means estimate is printed beside it as a check, and a series shorter than about 50 correlation times is flagged.

import numpy as np
import conatus as cn

rng = np.random.default_rng(0)                  # stands for the Nusselt number of a run,
nusselt = np.empty(20000)                       # one value every 0.01: a correlated series
nusselt[0] = 6.0                                # with a transient
for i in range(1, nusselt.size):
    nusselt[i] = 4.9 + 0.98 * (nusselt[i - 1] - 4.9) + 0.02 * rng.standard_normal()
rep = cn.time_average(nusselt, dt=0.01, discard=2000)
print(rep.explain())        # mean, standard error, correlation time, batch-means check