Subcritical radiative shock

The experiment

Physical setup

This test evolves gas and radiation together. The intended subcritical regime has a radiative precursor with T_gas,pre < T_downstream; radiation and gas temperatures need not agree in the precursor or at the shock spike.

Ensman's original calculation used a moving spherical piston at 6 × 10⁵ cm/s in a shell spanning 8 × 10¹¹ to 8.7 × 10¹¹ cm. Initially ρ = 7.78 × 10⁻¹⁰ g/cm³, the absorption coefficient is 3.115 × 10⁻¹⁰ cm⁻¹, and the 300 equal-mass zones begin at 10 K at the outer edge, rising by 0.25 K per zone inward. The gas uses γ = 5/3 and fully ionized hydrogen (mean molecular weight μ = 0.5).

Local Cartesian adaptation. The box spans −3.5 × 10¹⁰ to +3.5 × 10¹⁰ cm. Gas moves toward a reflecting left wall at v_x = −6 × 10⁵ cm/s, with fixed inflow at the right face and zero-gradient transverse boundaries. Temperature varies linearly from 85 K at the wall to 10 K at the inflow. Absorption follows χ = 3.115 × 10⁻¹⁰(ρ / 7.78 × 10⁻¹⁰ g/cm³) cm⁻¹. The default duration is 38300 s. A stationary planar wall replaces the spherical moving piston; the M1 moment closure replaces the original transfer treatment.

These cases enable gas dynamics and radiation coupling on a Cartesian grid; gravity is disabled. Uniform runs are the default; adaptive runs are identified in their recorded setup. Each completed run records its actual configuration, time, units, and executable hash below. The configured finest equivalent cells per side are N = INX × 2^level.

The comparison

Comparison and interpretation

No exact analytic profile is supplied for the coupled problem. Ensman (1994), Figure 9 at t = 3.83 × 10⁴ s, is an original numerical comparison target. Its curves are not embedded here, digitized, or used to manufacture error norms.

Approximate analytic shock scales

Let U be the imposed wall-frame inflow speed and R = k_B/(μm_H) = 1.65022 × 10⁸ erg/(g K). A strong, gas-pressure-dominated adiabatic shock with negligible upstream pressure has

ρ₂/ρ₁ = 4, u₁ = 4U/3, s = U/3
T₂ ≈ 2(γ − 1)u₁²/[R(γ + 1)²] = U²/(3R)

Here u₁ is the upstream speed relative to the shock and s is the outward shock speed in the wall frame. For U = 6 km/s, this gives s ≈ 2 km/s and T₂ ≈ 727 K. These are approximate conservation-law scales, not an analytic transient profile or an error threshold. Upstream heating, finite Mach number, and radiative energy losses modify the jump.

What this test measures

Inspect the radiation precursor, the difference between gas and radiation temperatures, the gas temperature overshoot, and the shock speed. A subcritical shock has gas immediately upstream cooler than the downstream equilibrium gas. The physical regime must be assessed after sufficient evolution; the case name alone does not demonstrate that the regime has been reached.

Measured results

All resolutions

Level	Grid	Batch status
2	32³	Ready

These are numerical solutions without a supplied analytic reference. No error norms or convergence orders are inferred from the profiles. Short runs used for implementation checks do not establish the published late-time shock or blast structure.

32³ cells · actual final time 3.830000e+04 s

Central cell layer at z=dx/2; ray at y=z=dx/2.

Numerical Cartesian adaptation. No independent reference profile or quantitative agreement with the original spherical benchmark is established.

M1 RK1 (Euler). Physical-light-speed radiation equations (alpha = 1); c = 29979245800 cm/s.

Recorded setup and run provenance

Origin	Octo-TIGER application
Domain	Cube of side 7.000000e+10 cm; centered at the origin
Resolution	32³ cells; finest Δx = 2.187500e+09 cm
Final comparison time	3.830000e+04 s
Radiation model	M1
Radiation integration	M1 RK1 (first order in time)
Flux semantics	Evolved M1 flux moment
Radiation backend	cpu
Physical light speed	2.997925e+10 cm/s
Reduced-speed factor α	1.000000e+00
Radiation transport speed	2.997925e+10 cm/s
Periodic boundaries	off
Hydrodynamics	on
Gravity	off
Implicit radiation source	on
Uniform grid	on
CFL	0.4
Timestep cap	63.833333333333336 s
Output interval	319.16666666666669 s

{
  "background": 0.0,
  "benchmark_time": 38300.0,
  "c": 29979245800.0,
  "case": "subcritical_shock",
  "cells": 32,
  "command": [
    "/ddnB/work/dmarce1/octotiger-m1-fld-20261014/campaign/control/launch-solver",
    "/ddnB/work/dmarce1/octotiger-m1-fld-20261014/build/octotiger",
    "--config_file=/ddnB/work/dmarce1/octotiger-m1-fld-20261014/campaign/runs/l2-subcritical_shock-m1/subcritical_shock/l2/run.ini",
    "--hpx:threads=64"
  ],
  "comparison_signature": {
    "atomic_mass": "1",
    "atomic_number": "1",
    "cfl": "0.4",
    "code_to_cm": "1",
    "code_to_g": "1",
    "code_to_s": "1",
    "correct_am_hydro": "off",
    "disable_analytic": "on",
    "disable_diagnostics": "on",
    "eos": "IDEAL",
    "gravity": "off",
    "hard_dt": "63.833333333333336",
    "hydro": "on",
    "hydro_device_kernel_type": "OFF",
    "hydro_host_kernel_type": "LEGACY",
    "inflow_bc": "off",
    "monopole_device_kernel_type": "OFF",
    "multipole_device_kernel_type": "OFF",
    "n_species": "1",
    "omega": "0",
    "periodic": "off",
    "problem": "RADIATION_SUBCRITICAL_SHOCK",
    "rad_implicit": "on",
    "radiation": "on",
    "radiation_algorithm_version": 1,
    "radiation_integrator": "euler",
    "radiation_model": "m1",
    "radiation_reduced_c_factor": "1",
    "reflect_bc": "off",
    "stop_time": "38300",
    "unigrid": "on",
    "xscale": "35000000000"
  },
  "config": {
    "atomic_mass": "1",
    "atomic_number": "1",
    "cfl": "0.4",
    "code_to_cm": "1",
    "code_to_g": "1",
    "code_to_s": "1",
    "correct_am_hydro": "off",
    "datadir": "/ddnB/work/dmarce1/octotiger-m1-fld-20261014/campaign/runs/l2-subcritical_shock-m1/subcritical_shock/l2/",
    "disable_analytic": "on",
    "disable_diagnostics": "on",
    "disable_output": "off",
    "eos": "IDEAL",
    "gravity": "off",
    "hard_dt": "63.833333333333336",
    "hydro": "on",
    "hydro_device_kernel_type": "OFF",
    "hydro_host_kernel_type": "LEGACY",
    "inflow_bc": "off",
    "max_level": "2",
    "min_level": "2",
    "monopole_device_kernel_type": "OFF",
    "multipole_device_kernel_type": "OFF",
    "n_species": "1",
    "odt": "319.16666666666669",
    "omega": "0",
    "periodic": "off",
    "problem": "RADIATION_SUBCRITICAL_SHOCK",
    "rad_implicit": "on",
    "radiation": "on",
    "radiation_backend": "cpu",
    "radiation_integrator": "euler",
    "radiation_model": "m1",
    "radiation_reduced_c_factor": "1",
    "reflect_bc": "off",
    "stop_time": "38300",
    "unigrid": "on",
    "xscale": "35000000000"
  },
  "conservation": [
    {
      "boundary": 1.7720122575424458e+25,
      "boundary_fraction": 0.5771086978222407,
      "field": "er",
      "final": 1.5399642032259397e+29,
      "final_fraction": 5015.353207822257,
      "history_end_time": 38300.0,
      "history_scope": "full_run",
      "history_start_time": 0.0,
      "initial": 3.070500001523553e+25,
      "initial_fraction": 1.0,
      "integral_units": "erg",
      "max_abs_residual": 9.58774139420672e+15,
      "max_normalized_error": 3.1225342418008057e-10,
      "normalization_basis": "abs(initial_energy)",
      "normalization_scale": 3.070500001523553e+25,
      "normalization_time": 0.0,
      "normalized_error": 2.308956512744449e-10,
      "normalized_residual": -2.308956512744449e-10,
      "raw_change": 1.5396571532257874e+29,
      "residual": -7.089650975899648e+15,
      "source": 1.5398343544516126e+29,
      "source_fraction": 5014.930316520311,
      "zero_initial_flux_relative_tolerance": 1.4210854715202004e-14
    },
    {
      "boundary": -1.3150101704590293e+44,
      "boundary_fraction": -5353404045784.618,
      "field": "fx",
      "final": 3.421872834434828e+38,
      "final_fraction": 13930438.172679234,
      "history_end_time": 38300.0,
      "history_scope": "full_run",
      "history_start_time": 0.0,
      "initial": -2.4564000012188425e+31,
      "initial_fraction": -1.0,
      "integral_units": "erg cm/s",
      "max_abs_residual": 1.9807040628566084e+28,
      "max_normalized_error": 0.0008063442687973464,
      "normalization_basis": "abs(initial_component)",
      "normalization_scale": 2.4564000012188425e+31,
      "normalization_time": 0.0,
      "normalized_error": 0.0,
      "normalized_residual": 0.0,
      "raw_change": 3.421873080074828e+38,
      "residual": 0.0,
      "source": -1.3150067485859492e+44,
      "source_fraction": -5353390115345.446,
      "zero_initial_flux_relative_tolerance": 1.4210854715202004e-14
    },
    {
      "boundary": 2.0937414811633094e+21,
      "boundary_fraction": 2.2745382914722038e-15,
      "field": "fy",
      "final": 0.0,
      "final_fraction": 0.0,
      "history_end_time": 38300.0,
      "history_scope": "full_run",
      "history_start_time": 0.0,
      "initial": 0.0,
      "initial_fraction": 0.0,
      "integral_units": "erg cm/s",
      "max_abs_residual": 1.2526726334734016e+22,
      "max_normalized_error": 1.3608422516095426e-14,
      "normalization_basis": "c * abs(initial_energy)",
      "normalization_scale": 9.205127427457496e+35,
      "normalization_time": 0.0,
      "normalized_error": 2.2745382914722038e-15,
      "normalized_residual": 2.2745382914722038e-15,
      "raw_change": 0.0,
      "residual": 2.0937414811633094e+21,
      "source": 0.0,
      "source_fraction": 0.0,
      "zero_initial_flux_relative_tolerance": 1.4210854715202004e-14
    },
    {
      "boundary": 1.2128770257259158e+22,
      "boundary_fraction": 1.3176102506827747e-14,
      "field": "fz",
      "final": 0.0,
      "final_fraction": 0.0,
      "history_end_time": 38300.0,
      "history_scope": "full_run",
      "history_start_time": 0.0,
      "initial": 0.0,
      "initial_fraction": 0.0,
      "integral_units": "erg cm/s",
      "max_abs_residual": 2.2828116007193212e+22,
      "max_normalized_error": 2.479934817534455e-14,
      "normalization_basis": "c * abs(initial_energy)",
      "normalization_scale": 9.205127427457496e+35,
      "normalization_time": 0.0,
      "normalized_error": 1.3176102506827747e-14,
      "normalized_residual": 1.3176102506827747e-14,
      "raw_change": 0.0,
      "residual": 1.2128770257259158e+22,
      "source": 0.0,
      "source_fraction": 0.0,
      "zero_initial_flux_relative_tolerance": 1.4210854715202004e-14
    }
  ],
  "dx": 2187500000.0,
  "executable_sha256": "ca38f12e338eed0b7074dbb15c70c089a8d31d1c12aa87f5b044c42824721f0b",
  "inx": 8,
  "length": 70000000000.0,
  "level": 2,
  "movie_capture": {
    "hard_dt": 63.833333333333336,
    "odt": 319.1666666666667,
    "requested_snapshots": 121,
    "steps_per_output_check": 5
  },
  "origin": "Octo-TIGER application",
  "physical_c": 29979245800.0,
  "radiation_algorithm_version": 1,
  "radiation_backend": "cpu",
  "radiation_equations": "physical-c",
  "radiation_flux_semantics": "evolved_m1_flux",
  "radiation_integrator": "euler",
  "radiation_model": "m1",
  "radiation_reduced_c_factor": 1.0,
  "reference_kind": "numerical-only",
  "status": "complete",
  "time": 38300.0,
  "transport_c": 29979245800.0,
  "units": {
    "chi": "1/cm",
    "energy_density": "erg/cm^3",
    "flux": "erg/(cm^2 s)",
    "length": "cm",
    "luminosity": "erg/s",
    "mass": "g",
    "source": "erg/(cm^3 s)",
    "system": "CGS",
    "time": "s"
  },
  "validation_scope": "Numerical Cartesian adaptation. No independent reference profile or quantitative agreement with the original spherical benchmark is established."
}

Numerical gas temperature, radiation temperature, and density. Central cell layer at z=dx/2; ray at y=z=dx/2. · PDF

Numerical profiles; no reference curves are supplied. Central cell layer at z=dx/2; ray at y=z=dx/2. · PDF

Saved-state movie

Subcritical Shock · L2 · M1 RK1

121 stored states; 121 shown. Physical time: 0–38300 s. Playback: 60 s.

Zero-order hold; no solution interpolation. Playback holds each saved numerical state over its time interval and includes endpoint holds.

All 121 snapshot times

Snapshot times: 0, 319.166667, 638.333333, 957.5, 1276.66667, 1915, 2234.16667, 2553.33333, 2872.5, 3191.66667, 3510.83333, 3830, 4149.16667, 4468.33333, 4787.5, 5106.66667, 5425.83333, 5745, 6064.16667, 6383.33333, 6702.5, 7021.66667, 7340.83333, 7660, 7979.16667, 8298.33333, 8617.5, 8936.66667, 9255.83333, 9575, 9894.16667, 10213.3333, 10532.5, 10851.6667, 11170.8333, 11490, 11809.1667, 12128.3333, 12447.5, 12766.6667, 13085.8333, 13405, 13724.1667, 14043.3333, 14362.5, 14681.6667, 15000.8333, 15320, 15639.1667, 15958.3333, 16277.5, 16596.6667, 16915.8333, 17235, 17554.1667, 17873.3333, 18192.5, 18511.6667, 18830.8333, 19150, 19469.1667, 19788.3333, 20107.5, 20426.6667, 20745.8333, 21065, 21384.1667, 21703.3333, 22022.5, 22341.6667, 22660.8333, 22980, 23299.1667, 23618.3333, 23937.5, 24256.6667, 24575.8333, 24895, 25214.1667, 25533.3333, 25852.5, 26171.6667, 26490.8333, 26810, 27129.1667, 27448.3333, 27767.5, 28086.6667, 28405.8333, 28725, 29044.1667, 29363.3333, 29682.5, 30001.6667, 30320.8333, 30640, 30959.1667, 31278.3333, 31597.5, 31916.6667, 32235.8333, 32555, 32874.1667, 33193.3333, 33512.5, 33831.6667, 34150.8333, 34470, 34789.1667, 35108.3333, 35427.5, 35746.6667, 36065.8333, 36385, 36704.1667, 37023.3333, 37342.5, 37661.6667, 37980.8333, 38300, 38300 s.

Arrows show the independently evolved M1 radiation flux. Arrows have equal length and show direction only.

Movies show saved numerical states; no applicable analytic reference is recorded.

Central z=0+ radiation energy slice with faint cell mesh and equal-length flux direction arrows. Arrows show the independently evolved M1 radiation flux. Gray density contours and a matching density/velocity panel. Synchronized gas and radiation temperature profiles along the central x line appear below the maps. Fixed color limits; actual stored Silo times held without solution interpolation. Movies show saved numerical states; no applicable analytic reference is recorded. The maps show the full simulation domain without cropping. The temperature profile also spans the full domain.

Download MP4 · Rendering and timing metadata · Full run report

All campaign movies

Data and reproducibility

Effective configuration Run metadata Raw conservation segment Solver log Numerical slice CSV Data unitsFinal numerical Silo (local file; excluded from staged site)

Reading the numerical profiles

Gas and radiation temperatures are distinct diagnostics; the radiation temperature is defined by E = a_rT_r⁴. The central ray lies at y = z = Δx/2 and is offset from the exact origin by half a cell in each transverse direction. The two-dimensional plane is z = Δx/2.

Check positivity, the recorded-frame radiation flux envelope, conservation with the appropriate boundary budget, thermal equilibration, shock speed and precursor structure, and sensitivity to resolution and timestep. A finite, positive short run is an implementation check, not evidence that a subcritical or supercritical state has developed. Comparisons to published figures also require matching geometry, opacities, boundaries, and physical time.

Sources & provenance

Where this problem comes from

Ensman (1994), “Test Problems for Radiation and Radiation-Hydrodynamics Codes,” ApJ 424, 275, §§5.1–5.2, pp. 283–285, describes the original radiative shocks. Original paper PDF. The original spherical geometry, piston boundary, and time must be matched before claiming quantitative agreement with its figures.

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.

Batch settings and source hashes Batch status