← Octo-TIGER / All test problemsSteady diffusion

Equilibrium sphere

A smooth spherical emitter balances outward radiation transport. An error-function solution describes its steady diffusion equilibrium.

The experiment

Physical setup

A Gaussian light source is centered in a stationary medium of prescribed opacity χ. It emits total luminosity ℒ over all space, with a smooth width w:

q(r) = ℒ exp(−r²/w²) / (π3/2 w³),   r = |x|

The domain starts from the steady diffusion solution on top of background Ebg. Boundaries are nonperiodic: ghost cells are fixed to cell averages of that same analytic equilibrium. The incoming background is part of this boundary state. The gas remains fixed; the source adds radiation energy and the medium damps flux.

The standard CGS templates use w = 6 × 10⁹ cm, χ = 1.6666666667 × 10−10 cm−1, Ebg = 1 erg/cm³, and ℒ = 2.698132122 × 10²⁹ erg/s. This is a smooth Gaussian emitter. The saved settings below give the actual values used.

The standard runner uses a uniform cubic grid with hydrodynamics and gravity disabled. Its CGS box spans −3 × 10¹⁰ to +3 × 10¹⁰ cm on each axis, with c = 2.99792458 × 10¹⁰ cm/s and default final time 4 s. Opt-in adaptive runs are identified below; N = INX × 2level then describes their configured finest equivalent resolution. Each completed run below lists its recorded parameters.

The comparison

Comparison and interpretation

In steady isotropic diffusion, D = c/(3χ), ∇ · F = q, and F = −D∇E. Integrating the enclosed Gaussian luminosity gives

E(r) = Ebg + [3χℒ/(4πc)] erf(r/w)/r
F(x) = [ℒ/(4πr³)] [erf(s) − 2s e−s²/√π] x,   s = r/w

The apparent singularities are removable: E(0) = Ebg + 3χℒ/(2π3/2cw) and F(0) = 0. The implementation uses series expansions near the origin.

Initial, boundary, and comparison states are cell averages evaluated with four-point Gauss–Legendre quadrature on each axis (64 samples per cell). The smooth-field averages therefore have quadrature error. The Gaussian source has an exact cell integral, evaluated as a product of three error-function differences:

q̄ = (ℒ/Δx³) ∏d=x,y,z ½ {erf[(xd + Δx/2)/w] − erf[(xd − Δx/2)/w]}

The prescribed source adds Δt q̄ to energy and damps flux at rate cχ. Its numerical integration depends on the recorded method, as described below. The steady reference itself is independent of time.

What this test measures

This is a diffusion-limit equilibrium test. Nonlinear M1 has an anisotropic pressure at finite flux, so the diffusion solution need not remain exactly stationary in that model. Errors can reflect closure differences, transport/source splitting, discretization, and quadrature. The plots measure equilibrium drift and Cartesian directional effects as resolution changes.

Energy injected inside the finite cube should balance its outward transport and the change in stored energy. ℒ is the all-space luminosity; use the measured in-domain source budget rather than assuming the finite cube receives exactly ℒ. Symmetry can make net flux integrals very small, which is why their conservation normalization includes an energy-based scale.

Measured results

All resolutions

No completed run is available for this problem in this batch. Its description is available above; results will appear here as runs finish.

Momentum plots are omitted for this fixed-material problem. The prescribed medium acts as an external momentum reservoir, so radiation momentum exchange does not imply that a compensating gas momentum was evolved.

Reading the results

Errors, budgets, and movies

Every norm uses the full three-dimensional domain, with cell-volume weighting. Here e = numerical − reference and V is the domain volume:

L1 = Σ |e| ΔV / V,   L2 = √(Σ e² ΔV / V),   L∞ = max |e|
p = ln(errorcoarse / errorfine) / ln(Δxcoarse / Δxfine)

Orders require two completed runs with matching physical time, configuration, units, executable, and origin. An exact zero error has no logarithmic order. The central slice and line profiles are visual diagnostics: the report slice is the cell layer at z = Δx/2, and its line profile is at y = z = Δx/2.

Conservation is assessed using the measured domain integrals and the signed cumulative boundary and source budgets:

A labeled checkpoint continuation uses its checkpoint as t_ref, with counters and fixed normalization starting there. Its reported residuals cover that resumed segment only; original histories remain separate.

R(t) = Q(t) − Q(0) + B(t) − S(t)

B is outward transport and S is the source change. Conservation plots and this page's budget tables show fractions of the initial totals: Q(t)/D, [Q(0) − B(t) + S(t)]/D, and the signed residual R(t)/D. The denominator D = |Q(0)| is fixed for the entire history. It never grows with the current total or accumulated transport/source budgets. A negative initial flux therefore starts at −1, and the plotted fraction retains its sign.

A component whose initial net flux is zero or negligible cannot be divided by its own initial total. It instead uses the fixed scale D = c |E₀|, where E₀ is the initial integrated radiation energy. This fallback is labeled on its panels and in the table. “Negligible” means |Q(0)|/(c |E₀|) ≤ 64ε ≈ 1.421085 × 10−14, with ε the double-precision machine epsilon. This prevents roundoff left by cancellation in a symmetric field from becoming the denominator. For the initially cold Marshak case, the prescribed bath sets D = Ebath V for energy and D = c Ebath V for flux, with Ebath = 4/c. These scales remain fixed even though radiation enters the initially empty domain.

Raw histories and dimensional summary columns remain in physical units: energy integrals in erg, physical-flux integrals in erg cm/s for CGS runs. Dividing a flux integral by c² gives radiation momentum. Fraction-history CSV files are dimensionless except for their time column; summary files record each denominator and its basis. Maximum fractional residuals use the same fixed denominator as the final residual.

For runs labeled RSLA, these budgets check the modified radiation equations. They do not establish conservation of physical gas-plus-radiation energy. With reduced-speed factor α, the closed-system coupled invariants are Egas + Er/α and gas momentum + F/(αc²), with physical c in the closure. Compare multiple α values before interpreting an RSLA result as an approximation to physical-c evolution.

Movies show the numerical Silo states with a fixed color range for each movie. Recorded states are held for their actual simulation-time gaps; intermediate solutions are not interpolated. The 3D option shows the domain's exterior surface. Movies and plots retain their recorded physical units.

Sources & provenance

Where this problem comes from

This smooth Gaussian-emitter test is defined by equilibriumSphere, sphereAverage, and bulbSourceAverage in octotiger/test_problems/radiation/profiles.hpp. Its reference follows by solving the steady diffusion equation for that source. The available implementation does not identify it as a verbatim published sphere benchmark, so the provenance here is the suite's analytic construction.

Radiation method context

M1 RK1 evolves radiation energy and flux moments with the retained first-order radiation update. FLD uses backward-Euler radiation diffusion/material exchange with the Levermore–Pomraning limiter. Moving cases add first-order conservative gas transport and radiation force/work; their reconstructed lab-frame flux includes (E I+P) v as well as −D∇E. The scalar implicit transport solve uses composite-grid multigrid and an iterative linear solver. Both time integrators are first order.

The suite includes Gaussian scattering, Marshak heating, a Gaussian emitter, explicit thin positive-opacity wave/front comparisons, and moving-gas shock/blast problems on CPU at physical light speed. A configured case does not establish that a solve converged; only validated completed runs appear as results. Flux diagnostics in FLD have no separate M1 momentum conservation test. A difference between paired M1 and FLD solutions is a model comparison, not an analytic error. ADER2 has been removed from the current solver; historical result files retain their original provenance.

Skinner & Ostriker (2013) and Hanawa & Audit (2014) provide context for the M1 equations. The implementation and saved executable/source hashes identify the actual numerical methods used in this batch.

Reproducing this batch

Descriptions document the reviewed source implementation. Older or edited solver versions can use different profile constants or averaging rules; these are not all encoded in run.ini. Use the run metadata, executable hash, and saved batch source hashes to identify the implementation behind a result. A source path or descriptive formula alone is not proof of the executable's contents.