Experimental

Pyrolysis.Experimental

Namespace for experimental or unstable code paths. The names below remain functional and are accessible as Pyrolysis.Experimental.X, but they are deliberately removed from the top-level public API. These interfaces may change between releases; prefer the stable top-level exports when possible.

Contents:

• P1 quasi-steady radiation (residual wrapper + G kernels). Lives in RadiationP1.jl — built on top of the unified (def, ws) pair and the unified StateLayout; the quasi-steady G field lives in the RadiationCache, not in the ODE state vector.

• h-AMR (h-refinement, h-coarsening, AMR controller, adjacent-cell merging, surface-depletion utilities). Will move to a dedicated AdaptivityH extension once the active-cell desync invariants are reworked.

AMRController{EI,QM}

Controls adaptive mesh refinement operations.

Type Parameters

• EI: Error indicator type
• QM: Quality metric type

Fields

• error_indicator: Error indicator for refinement decisions
• quality_metric: Mesh quality metric to enforce
• θ_refine: Refinement threshold (refine if η > θ_refine * max(η))
• θ_coarsen: Coarsening threshold (coarsen if η < θ_coarsen * max(η))
• h_min: Minimum allowed cell size
• h_max: Maximum allowed cell size
• max_level: Maximum refinement levels
• check_every: Timesteps between AMR checks
• n_refinements: Counter for total refinements performed
• n_coarsenings: Counter for total coarsenings performed
• last_check_time: Last simulation time when AMR was checked

AMRController(; kwargs...)

Construct an AMR controller with keyword arguments.

Keyword Arguments

• error_indicator: Error indicator (default: GradientIndicator(:T))
• quality_metric: Quality metric (default: AspectRatioMetric())
• θ_refine: Refinement threshold (default: 0.7)
• θ_coarsen: Coarsening threshold (default: 0.1)
• h_min: Minimum cell size (default: 1e-6)
• h_max: Maximum cell size (default: 1e-2)
• max_level: Maximum refinement level (default: 4)
• check_every: Check frequency in timesteps (default: 10)

DirichletRadiationBC{FG} <: AbstractRadiationBC

Prescribed incident radiation intensity at a boundary.

Sets the value of G directly at the boundary: Gboundary = Gprescribed

This is rarely used in practice but can be useful for:

• Coupling to external radiation field calculations
• Validation against analytical solutions
• Specialized boundary treatments

Fields

• G::FG: Prescribed G value [W/m²]. Can be constant or time-dependent G(t).

Physical Meaning

For any surface at temperature Tb in volumetric equilibrium: G = 4σTb⁴ (blackbody, standard P1 approximation)

Note: In P1, the volumetric equilibrium is ALWAYS blackbody (ε=1 locally). Surface emissivity ε < 1 only appears in boundary conditions (Marshak BC) where it represents partial reflection at the boundary, not in G itself.

Example

# Isothermal black surface at 500 K
bc = DirichletRadiationBC(G = 4 * STEFAN_BOLTZMANN * 500.0^4)

# Time-varying prescribed G
bc = DirichletRadiationBC(G = t -> 1000.0 * (1 + 0.5*sin(2π*t/60)))

Note

Using DirichletRadiationBC requires careful consideration of the physics. In most pyrolysis applications, MarshakRadiationBC is preferred as it properly couples external sources with surface emission.

See Also

• MarshakRadiationBC for typical exposed surfaces
• DirichletThermalBC for prescribed temperature

MarshakRadiationBC{FI, FT, Fε} <: AbstractRadiationBC

Marshak boundary condition for P1 radiation at a diffuse surface.

Models radiation exchange between a surface and its surroundings, including external radiation sources (e.g., cone heater, fire). This is the primary BC for exposed surfaces in pyrolysis simulations.

Fields

• I_external::FI: External incident radiation flux [W/m²]. Can be constant or time-dependent function I_external(t). This represents radiation from external sources like cone heaters, radiant panels, or flames.
• T_ambient::FT: Ambient temperature [K] for radiation exchange with surroundings. Can be constant or time-dependent T_ambient(t).
• emissivity::Fε: Surface emissivity ε ∈ [0,1]. Can be constant or temperature-dependent emissivity(T). For pyrolysis materials:
  • Virgin polymers: 0.85-0.95
  • Char/carbon: 0.90-0.98

Physical Model

The Marshak condition closes the P1 half-range fluxes at a gray diffuse surface: the inward half-range flux equals the transmitted external radiosity plus the reflected outgoing half-range flux,

G/4 + q/2 = ε(I_external + σT_amb⁴) + (1-ε)(G/4 - q/2)

Solving for the net flux into the domain:

q_rad = c · (4σT_amb⁴ + 4·I_external - G_surf),    c = ε/(2(2-ε))

This accounts for:

• External radiation input (weight 4c = 2ε/(2-ε) on I_external — the same half-range weight the ambient term carries)
• Ambient radiation exchange (c·4σT_amb⁴)
• Partial reflection (1-ε) of the interior radiation back into the domain

The surface temperature does not appear: the medium's own emission is volumetric (4σαT⁴) in P1, and the surface returns radiation only by reflection. See radiation_bc_flux and _marshak_exchange_coefficient for the implementation and derivation.

Examples

# Cone calorimeter with 25 kW/m² heater (constant)
bc = MarshakRadiationBC(
	I_external = 25e3,      # W/m²
	T_ambient = 300.0,      # K
	emissivity = 0.9
)

# Time-varying external flux (ramp-up)
bc = MarshakRadiationBC(
	I_external = t -> 50e3 * (1 - exp(-t/10)),
	T_ambient = 300.0,
	emissivity = 0.9
)

# Temperature-dependent emissivity
bc = MarshakRadiationBC(
	I_external = 35e3,
	T_ambient = 300.0,
	emissivity = T -> 0.85 + 0.0001 * (T - 300)  # Increases with T
)

See Also

• ReflectiveRadiationBC for back surfaces
• RadiativeBC (thermal BC) for surface-only radiation without P1

P1Residual{D, W, RC, RBC}

Callable struct holding everything needed to evaluate the P1 residual: the immutable def, the mutable ws, the radiation cache (which owns the quasi-steady G field), and the radiation BCs.

The state vector is the unified layout from ws.layout — build u0 with create_state_vector(ws.layout, ws.mesh). Read the G field via p.rad_cache.G or get_radiation_intensity.

The constructor requires def.options.radiation_model == P1_QUASI_STEADY (use to_problem_def(problem; radiation_model = P1_QUASI_STEADY)); this is what keeps the incident radiation single-counted between the Marshak BC and the thermal boundary conditions.

Construct explicitly in Pyrolysis.Experimental; the core solve path keeps P1 out of the unified residual until it gets a dedicated radiation trait.

RadiationBCSet{B1, B2}

Type-stable container for radiation boundary conditions (top, bottom).

This mirrors BoundaryConditionSet for thermal/mass BCs and avoids Dict lookups inside hot radiation kernels.

RadiationBCSet(; top=ReflectiveRadiationBC(), bottom=ReflectiveRadiationBC())

Keyword constructor with defaults.

ReflectiveRadiationBC <: AbstractRadiationBC

Reflective (adiabatic) boundary condition for radiation.

Imposes zero normal gradient on the incident radiation field: ∂G/∂n = 0

This is appropriate for:

• Insulated back faces (no radiation escape)
• Symmetry planes
• Reflective substrates
• Impermeable boundaries where radiation cannot escape

Physical Interpretation

A reflective BC means all radiation incident on the boundary is reflected back into the domain. No net radiative flux crosses the boundary.

Example

# Back face of sample (insulated)
radiation_bcs = RadiationBCSet(
	top = MarshakRadiationBC(I_external = 25e3, T_ambient = 300.0),
	bottom = ReflectiveRadiationBC()  # Adiabatic
)

See Also

• MarshakRadiationBC for exposed surfaces
• AdiabaticBC (thermal BC) for thermal insulation

_marshak_driving_field(bc::MarshakRadiationBC, t) -> G_ext + 4·I_ext

The driving incident field of the Marshak relation q = c·(G_ext + 4·I_ext − G): the ambient equilibrium field Gext = 4σTamb⁴ plus the external incident flux at its half-range weight 4 (so that c·4·Iext = 2ε/(2−ε)·Iext).

_marshak_effective_coefficient(bc, T_cell, D_cell, δ) -> c_eff

Effective Marshak exchange coefficient seen by the cell-centered value GN of the boundary cell: the surface exchange conductance c = ε/(2(2−ε)) in series with the half-cell diffusive conductance β = Drad/δ (δ = distance from the cell center to the boundary face):

c_eff = c·β/(c + β) = 1 / (1/c + δ/D_rad)

so the applied flux q = c_eff·(G_ext + 4·I_ext − G_N) is algebraically identical to the exact face relation q = c·(G_ext + 4·I_ext − G_face) closed with the discrete diffusive flux q = Drad·(Gface − GN)/δ. This removes the O(δ) error of evaluating the Marshak relation at the cell center. Used by both the flux application (`applymarshakbc!`) and the quasi-steady Newton Jacobian — the two must stay consistent.

_marshak_exchange_coefficient(bc::MarshakRadiationBC, T_surf) -> c

The Marshak G-exchange coefficient c = ε/(2(2−ε)) (Modest, Radiative Heat Transfer, 3rd ed., Eq. 16.48), with ε evaluated at T_surf.

adapt_mesh!(mesh::Mesh1D{NC}, controller::AMRController;
			material=nothing) where NC -> NamedTuple

Perform full adaptive mesh refinement cycle.

material is forwarded to the coarsening pass for the energy-conserving temperature merge; without it sibling merges fall back to volume-weighted-T averaging (not energy-conserving for mixed compositions).

Algorithm

1. Compute error indicators for all active cells
2. Mark cells for refinement (η > θrefine * max(η), h > hmin, level < max_level)
3. Enforce 2:1 balance on marked cells
4. Perform h-refinement on marked cells
5. Recompute indicators on refined mesh
6. Mark sibling pairs for coarsening (both η < θcoarsen * max(η), combined h < hmax)
7. Perform h-coarsening on marked pairs
8. Check mesh quality; apply smoothing if needed

Returns

NamedTuple with:

• n_refined: Number of cells refined
• n_coarsened: Number of cell pairs coarsened
• quality_ok: Whether mesh passes quality check

amr_statistics(controller::AMRController) -> NamedTuple

Get AMR statistics from controller.

apply_radiation_divergence!(dT::Vector{Float64}, mesh, rad_cache, prop_cache)

Add radiation flux divergence contribution to temperature RHS.

dT[i] += (1/ρcₚ[i]) × (1/V[i]) × (q_rad[R] - q_rad[L]) × A

This implements the ∇·(D_rad∇G) term in the energy equation.

Arguments

• dT: Temperature RHS vector [K/s] (modified in-place)
• mesh: Mesh with geometry
• rad_cache: Radiation cache with computed q_rad fluxes
• prop_cache: Property cache with ρcₚ values

Notes

The radiation flux qrad must be computed first using `computeradiation_fluxes!`. Sign convention: positive flux = into domain at boundaries, in +z direction for interior.

Example

# First compute radiation fluxes
update_radiation_properties!(rad_cache, mesh, prop_cache)
compute_radiation_fluxes!(rad_cache, mesh, prop_cache, radiation_bcs, t)

# Then add to temperature RHS
apply_radiation_divergence!(dT, mesh, rad_cache, prop_cache)

can_coarsen(mesh::Mesh1D, cell_id::Int) -> Bool

Check if a cell can be coarsened (merged with sibling).

A cell can be coarsened if:

1. It has a parent (is not a root cell)
2. It is a leaf (has no children)
3. Its sibling is also a leaf

check_p1_validity(mesh, prop_cache::PropertyCache; 
				  warn::Bool=true) -> NamedTuple

Check if P1 approximation is valid for the current material state.

Returns a NamedTuple with:

• τ_total: Total optical thickness
• τ_min: Minimum cell optical thickness
• τ_max: Maximum cell optical thickness
• regime: String describing optical regime
• valid: Bool indicating if P1 is appropriate

Arguments

• mesh: Mesh with geometry
• prop_cache: Property cache with α values
• warn: If true, print warning for marginal cases

Example

validity = check_p1_validity(mesh, prop_cache)
println("Optical regime: $(validity.regime)")
println("P1 valid: $(validity.valid)")

coarsen_cells!(mesh::Mesh1D{NC}, cell_a_id::Int, cell_b_id::Int; material=nothing) -> Int

Merge two sibling cells back into their parent.

Arguments

• mesh: The mesh to modify
• cell_a_id: First child cell index
• cell_b_id: Second child cell index (sibling)
• material=nothing: Optional material for energy-conserving temperature merge. If provided, uses ρcₚ-weighted energy conservation. Otherwise, uses simple volume-weighted averaging.

Returns

• Parent cell index

Algorithm

1. Verify cells are siblings and both are leaves
2. Compute conservative average of child states
3. Update parent state
4. Restore parent face connectivity
5. Mark children as inactive
6. Mark parent as active

Energy Conservation Note

When material is provided, temperature is computed via energy conservation: Etotal = Va × ρcₚa × Ta + Vb × ρcₚb × Tb Tparent = Etotal / (Vparent × ρcₚ_parent)

This is critical when merging cells with different compositions (e.g., solid vs gas), as simple volume-weighted averaging can introduce significant temperature errors.

Note

This does NOT remove the child cells or intermediate node from the mesh arrays. They are simply marked inactive. This avoids index invalidation.

compute_p1_residual!(du, u, def, t, ws, p1) -> nothing

Unified residual augmented with P1 quasi-steady volumetric radiation.

u is an ordinary unified-layout state vector (T, ξ, z when ALE, χ when WithChi) — G is not a state variable. p1 provides rad_cache and radiation_bcs (a P1Residual or any object with those fields).

Steps:

1. Evaluate the unified residual!. def must carry radiation_model = P1_QUASI_STEADY (enforced), which zeroes the boundary RadiativeFluxBC/RadiativeBC fluxes and the core volumetric radiation source — P1 owns all radiation energy.
2. Mirror u into the legacy Mesh1D (the P1 kernels are Float64-only and read cell temperatures and live geometry from the mesh).
3. Update radiation properties (α, D_rad) and solve the quasi-steady P1 system for G. This happens on every evaluation — G is an algebraic field, never frozen.
4. Add S_rad/ρcₚ = α(G − 4σT⁴)/ρcₚ to the energy rows of du.

compute_radiation_diffusion_coefficient(α::Float64, σ_s::Float64=0.0) -> Float64

Compute the P1 radiation diffusion coefficient.

D_rad = 1 / (3 · max(α + σ_s, 1.0e-10))

For pyrolysis materials, scattering is typically negligible (σs ≈ 0): Drad = 1 / (3α)

Arguments

• α: Absorption coefficient [m⁻¹].
• σ_s: Scattering coefficient [m⁻¹] (default: 0). Typically zero for pyrolysis.

Returns

• D_rad: Radiation diffusion coefficient [m]

Physical Interpretation

Drad represents the "mean free path" for radiation diffusion. Large Drad (small α) means radiation travels far before being absorbed. Small D_rad (large α) means radiation is absorbed locally.

Transparent limit

Below an extinction of P1_MIN_EXTINCTION the coefficient is clamped to the large finite value 1/(3·P1_MIN_EXTINCTION) rather than diverging (or being zeroed): a transparent cell then behaves as a near-perfect conductor of G, which keeps the quasi-steady Newton system regular. P1 itself is not a valid model in this regime (τ < 0.1) — see check_p1_validity.

Examples

# Typical virgin polymer (α ≈ 1000-2000 m⁻¹)
D_rad = compute_radiation_diffusion_coefficient(1500.0)  # ≈ 2.2×10⁻⁴ m

# Char layer (α ≈ 5000-10000 m⁻¹)
D_rad = compute_radiation_diffusion_coefficient(8000.0)  # ≈ 4.2×10⁻⁵ m

# With scattering (rare for pyrolysis)
D_rad = compute_radiation_diffusion_coefficient(5000.0, 1000.0)

compute_radiation_flux_divergence!(rad_cache::RadiationCache, mesh, prop_cache::PropertyCache)

Compute the net volumetric radiation source from the P1 model.

S_rad = α(G - 4σT⁴)

where:

• α is volumetric absorption coefficient [1/m]
• G is incident radiation (the quasi-steady P1 field) [W/m²]
• σ is Stefan-Boltzmann constant
• T is temperature [K]

Sign Convention

• Positive S_rad: net absorption (G > 4σT⁴) → material heats up
• Negative S_rad: net emission (G < 4σT⁴) → material cools down

Physics Derivation

Integrating the RTE over all solid angles gives ∇·qrad = κ(4σT⁴ - G). The energy equation has -∇·qrad = κ(G - 4σT⁴) = Srad as a source term. The result is stored in `radcache.S_rad` (not the flux divergence itself).

Note: For gray medium in LTE (Kirchhoff's law), the absorption and emission coefficients are equal, so no separate emissivity factor appears here. Emissivity only enters at boundaries via Marshak BCs.

compute_radiation_fluxes!(rad_cache::RadiationCache,
						  mesh::Mesh1D{NC},
						  prop_cache::PropertyCache,
						  radiation_bcs::RadiationBCSet,
						  t) where NC

Compute radiation diffusion fluxes at all faces using P1 approximation.

For interior faces: qrad[f] = -Drad_face × (G[R] - G[L]) / Δx

For boundary faces, the Marshak or other BC is applied.

Arguments

• rad_cache: Radiation cache (updated in-place: q_rad)
• mesh: Mesh with geometry
• prop_cache: Property cache with α, ε values
• radiation_bcs: RadiationBCSet with radiation boundary conditions
• t: Current simulation time [s] (supports dual numbers for AD)

Notes

• Type-generic to support automatic differentiation

compute_radiation_source_terms!(Q_emission::Vector{Float64},
								Q_absorption::Vector{Float64},
								mesh::Mesh1D{NC},
								rad_cache::RadiationCache,
								prop_cache::PropertyCache) where NC

Compute volumetric emission and absorption source terms for P1 equation.

Emission Source (blackbody emission into radiation field):

Q_emission[i] = 4σ × ε[i] × α[i] × T[i]⁴    [W/m³]

Absorption (removal from radiation field):

Q_absorption[i] = α[i] × G[i]                [W/m³]

The net radiation source for the energy equation is: Qnet = Qabsorption - Q_emission

Note: These are the sources for the G equation. For the T equation, the radiation contribution comes through ∇·(D_rad ∇G).

Arguments

• Q_emission: Output vector for emission sources [W/m³]
• Q_absorption: Output vector for absorption sources [W/m³]
• mesh: Mesh with current temperature state
• rad_cache: Radiation cache with G values
• prop_cache: Property cache with α, ε values

create_default_controller(; h_target::Float64=1e-4) -> AMRController

Create an AMR controller with default settings tuned for pyrolysis.

The default indicator combines temperature gradient and curvature (curvature targets the sharp fronts that reaction zones produce). A ReactionRateIndicator cannot be used here: adapt_mesh! recomputes indicators after topology changes, when any precomputed rates are stale, so that indicator is standalone-only.

Arguments

• h_target: Target cell size in reaction zones

default_radiation_bcs(I_external::Float64=0.0; T_ambient::Float64=300.0)

Create default radiation boundary conditions for a 1D domain.

Returns a RadiationBCSet with:

• Top (boundary 1): MarshakRadiationBC with specified external flux
• Bottom (boundary 2): ReflectiveRadiationBC (adiabatic)

Arguments

• I_external: External radiation flux at top surface [W/m²]
• T_ambient: Ambient temperature for radiation exchange [K]
• emissivity: Surface emissivity (default 0.9)

Example

# Default for cone calorimeter at 50 kW/m²
radiation_bcs = default_radiation_bcs(50e3, T_ambient=300.0)

enforce_2to1_balance!(mesh::Mesh1D, marked::Vector{Int}, max_level::Int) -> Vector{Int}

Ensure 2:1 refinement level balance by marking additional cells.

Adjacent cells should differ by at most 1 refinement level.

Returns

• Expanded vector of cells to refine (including balance cells)

estimate_gradient(mesh::Mesh1D, cell_id::Int, φ::Float64, var) -> Float64

Estimate gradient of a field at a cell for interpolation.

get_radiation_intensity(p::P1Residual, i::Int) -> Float64

Quasi-steady incident radiation G in cell i [W/m²], from the most recent residual evaluation. G is an algebraic field stored in the radiation cache — it is not part of the ODE state vector.

get_sibling(mesh::Mesh1D, cell_id::Int) -> Int

Get the sibling cell ID (the other child of the same parent). Returns 0 if no sibling exists.

hr_adapt!(mesh::Mesh1D{NC};
		  h_indicator::AbstractErrorIndicator=GradientIndicator(:T),
		  r_monitor::MonitorFunction=GradientMonitor(),
		  θ_refine::Float64=0.7, θ_coarsen::Float64=0.1,
		  h_min::Float64=1e-6, h_max::Float64=1e-2,
		  max_level::Int=4, r_iterations::Int=5,
		  material=nothing) where NC

Perform combined h-refinement and r-refinement.

material is forwarded to the coarsening pass for the energy-conserving temperature merge; without it sibling merges fall back to volume-weighted-T averaging (not energy-conserving for mixed compositions).

Strategy (from Design Doc Section 4.6)

1. Perform h-adaptivity first (change topology)
2. Perform r-adaptivity second (smooth and optimize)

initialize_g_from_temperature!(rad_cache, mesh, prop_cache)

Initialize G field from local equilibrium with temperature.

For local thermodynamic equilibrium in P1: G = 4σT⁴ (blackbody)

The standard P1 approximation uses blackbody emission (no emissivity factor) in the volumetric source term. At equilibrium: αG = 4σαT⁴ → G = 4σT⁴

Arguments

• rad_cache: Radiation cache (G is updated in-place)
• mesh: Mesh with temperature state
• prop_cache: Property cache (not used, kept for API compatibility)

mark_cells_for_coarsening(mesh::Mesh1D, η::Vector{Float64},
						  θ_coarsen::Float64, h_max::Float64)
						  -> Vector{Tuple{Int,Int}}

Mark sibling cell pairs for coarsening based on error indicators.

Arguments

• mesh: The mesh
• η: Error indicator values
• θ_coarsen: Coarsening threshold (coarsen if η < θ_coarsen * max(η))
• h_max: Maximum allowed cell size

Returns

• Vector of (cella, cellb) tuples representing sibling pairs to coarsen

mark_cells_for_refinement(mesh::Mesh1D, η::Vector{Float64}, 
						  θ_refine::Float64, h_min::Float64, max_level::Int)
						  -> Vector{Int}

Mark cells for refinement based on error indicators.

Arguments

• mesh: The mesh
• η: Error indicator values
• θ_refine: Refinement threshold (refine if η > θ_refine * max(η))
• h_min: Minimum allowed cell size
• max_level: Maximum refinement level

Returns

• Vector of cell indices to refine

maybe_adapt!(mesh::Mesh1D{NC}, controller::AMRController, 
			 t::Float64, step::Int) where NC -> Union{Nothing,NamedTuple}

Conditionally perform adaptation based on check frequency.

Returns

• nothing if adaptation was skipped
• Result of adapt_mesh! if adaptation was performed

optical_thickness(mesh::Mesh1D{NC}, prop_cache::PropertyCache) where NC -> Float64

Compute the optical thickness τ = α × L for the domain.

This is a key parameter for assessing P1 validity:

• τ > 10: Optically thick (P1 excellent)
• 1 < τ < 10: Intermediate (P1 good)
• 0.1 < τ < 1: Optically thin (P1 fair)
• τ < 0.1: Transparent (P1 poor)

Arguments

• mesh: Mesh with geometry
• prop_cache: Property cache with α values

Returns

• τ: Optical thickness (dimensionless)

Example

τ = optical_thickness(mesh, prop_cache)
if τ < 0.1
	@warn "Material is nearly transparent (τ = $τ). P1 may be inaccurate."
end

optical_thickness_cell(mesh::Mesh1D{NC}, prop_cache::PropertyCache, i::Int) where NC -> Float64

Compute the optical thickness of a single cell: τi = αi × Δx_i

This is useful for identifying which cells are optically thick vs thin.

perform_h_coarsening!(mesh::Mesh1D{NC}, pairs::Vector{Tuple{Int,Int}};
					  material=nothing) where NC

Perform h-coarsening on all marked sibling pairs.

Arguments

• mesh: The mesh to modify
• pairs: Vector of (cella, cellb) tuples to coarsen
• material=nothing: Material for the energy-conserving temperature merge. Without it every pair falls back to volume-weighted-T averaging, which does not conserve energy when the siblings have different compositions — pass the material whenever one is available.

Notes

• After coarsening, call update_geometry!(mesh) to refresh caches

perform_h_refinement!(mesh::Mesh1D{NC}, marked::Vector{Int}) where NC

Perform h-refinement on all marked cells.

Arguments

• mesh: The mesh to modify
• marked: Vector of cell indices to refine

Notes

• Cells are refined in order of decreasing index to avoid invalidating indices
• After refinement, call update_geometry!(mesh) to refresh caches

radiation_bc_flux(bc::DirichletRadiationBC, args...)

For Dirichlet BCs, the flux is determined by the interior gradient. This returns zero as a placeholder - actual implementation needs the interior G value.

radiation_bc_flux(bc::MarshakRadiationBC, G_surf, T_surf, α_surf, ε_surf, t)

Compute radiation flux at boundary with Marshak condition, evaluated at the surface (face) value of G.

The P1 Marshak condition for a gray diffuse surface with emissivity ε:

q = c · (G_ext + 4·I_ext - G_surf),    c = ε / (2(2-ε))

where:

• Iext is external incident radiation [W/m²], carrying its half-range weight 4c = 2ε/(2-ε) — the same weight the ambient term carries (c·Gext = c·4σTamb⁴ ≡ 4c·σTamb⁴). It is 2.0 at ε = 1 and 1.64 at ε = 0.9; a weight of 1 under-injects the external flux by that factor.
• Gext = 4σTamb⁴ is the equilibrium incident radiation from ambient
• G_surf is the P1 incident radiation at the surface
• ε is surface emissivity (1-ε is reflected back into domain)

Derived from the P1 half-range fluxes qin = G/4 + q/2, qout = G/4 − q/2 with the surface closure qin = ε(Iext + σTamb⁴) + (1−ε)·qout (see the kernel comments above _marshak_exchange_coefficient).

This formulation ensures:

1. At equilibrium (G = 4σT⁴, Tamb = Tsurf, I_ext = 0): q = 0 ✓
2. External and ambient radiation enter with consistent half-range weights
3. Surface partially reflects (1-ε) of interior radiation back

The production solver applies this relation at the boundary-cell center through the half-cell-corrected coefficient (_marshak_effective_coefficient); this face-value form is the exact pointwise relation.

Arguments

• bc: Marshak boundary condition
• G_surf: Incident radiation at the surface (face) [W/m²]
• T_surf: Surface temperature [K] (used only to evaluate ε(T))
• α_surf: Surface absorption coefficient [m⁻¹] (not used)
• ε_surf: Surface emissivity [-] (not used — the BC's own emissivity governs)
• t: Current time [s] (supports dual numbers for AD)

Returns

• Radiation flux [W/m²], positive = into domain

Notes

• Type-generic to support automatic differentiation with Dual numbers
• Uses BC emissivity, not εsurf from propcache, for consistency

radiation_bc_flux(::ReflectiveRadiationBC, args...)

Zero flux for reflective boundary (∂G/∂n = 0).

radiation_energy_balance(rad_cache::RadiationCache, mesh) -> NamedTuple

Compute radiation energy balance for the domain.

Returns:

• total_emission: ∫ Q_emission dV [W/m²] (integrated over depth)
• total_absorption: ∫ Q_absorption dV [W/m²]
• net_source: totalabsorption - totalemission [W/m²]
• flux_top: Radiation flux at top boundary [W/m²]
• flux_bottom: Radiation flux at bottom boundary [W/m²]
• balance_error: (fluxtop − fluxbottom) + net_source — zero at a converged quasi-steady state

The quasi-steady G equation integrates to ∫(Qabs − Qem)dV + ∮q·n̂ dA = 0 with q_rad storing the physical +z flux at every face (boundary faces included), so the discrete residual is (flux_top − flux_bottom) + net_source. The two terms cancel — they do not equal each other.

refine_cell!(mesh::Mesh1D{NC}, cell_id::Int) -> Tuple{Int,Int}

Split a single cell into two children.

Arguments

• mesh: The mesh to modify
• cell_id: Index of cell to refine

Returns

• Tuple of (childaid, childbid)

Algorithm

1. Create new node at cell center
2. Create two child cells (left and right)
3. Create new interior face between children
4. Update face connectivity
5. Interpolate state conservatively to children
6. Mark parent as inactive

reset_statistics!(controller::AMRController)

Reset AMR statistics counters.

resize_p1_residual!(p::P1Residual, n_nodes_new, n_cells_new, n_faces_new)

Resize the P1 residual after AMR / depletion: walks the workspace (which rebuilds ws.layout) and the radiation cache.

should_adapt(controller::AMRController, t::Float64, step::Int) -> Bool

Check if adaptation should be performed at current time/step.

Arguments

• controller: AMR controller
• t: Current simulation time
• step: Current timestep number

Returns

• true if adaptation should be performed

solve_quasi_steady_radiation!(rad_cache, mesh, prop_cache, radiation_bcs, t;
							  tol=1e-10, maxiter=20, damping=1.0)

Solve the quasi-steady P1 radiation equation for G.

Since radiation equilibrates much faster than conduction/reaction time scales (τrad ≈ 10⁻¹² s vs τcond ≈ 10² s), we can assume dG/dt = 0:

∇·(D_rad∇G) = αG - 4σαT⁴

This is a nonlinear algebraic equation (due to T⁴) solved using Newton-Raphson.

Arguments

• rad_cache: Radiation cache (G is updated in-place)
• mesh: Mesh with temperature state
• prop_cache: Property cache (must be updated before calling)
• radiation_bcs: RadiationBCSet with top/bottom radiation BCs
• t: Current time [s] (supports dual numbers for AD)
• tol: Convergence tolerance (default 1e-10)
• maxiter: Maximum iterations (default 20)
• damping: Damping factor for updates (default 1.0)

Returns

• converged::Bool: True if converged within tolerance
• iterations::Int: Number of iterations taken

Algorithm

Newton-Raphson iteration:

1. Assemble residual R[i] = ∇·(D_rad∇G) - αG + 4σαT⁴
2. Assemble Jacobian J[i,j] = ∂R[i]/∂G[j]
3. Solve J × ΔG = -R
4. Update G = G + damping × ΔG
5. Check convergence: ||ΔG|| < tol

For 1D, the Jacobian is tridiagonal and can be solved efficiently.

Notes

• Type-generic to support automatic differentiation

Example

# First update radiation properties
update_radiation_properties!(rad_cache, mesh, prop_cache)

# Then solve for quasi-steady G
converged, iters = solve_quasi_steady_radiation!(
	rad_cache, mesh, prop_cache, radiation_bcs, t)

if !converged
	@warn "Quasi-steady radiation did not converge"
end

thomas_solve!(dl, d, du, x, d_work)

Solve tridiagonal system Ax = b in-place using Thomas algorithm (TDMA).

The system has the form: [d₁ du₁ 0 ... 0 ] [x₁] [b₁] [dl₁ d₂ du₂ ... 0 ] [x₂] [b₂] [0 dl₂ d₃ ... 0 ] [x₃] = [b₃] [... ... ] [..] [..] [0 ... dl_{n-1} dₙ ] [xₙ] [bₙ]

Arguments

• dl::Vector{Float64}: Lower diagonal (length n-1), unchanged on exit
• d::Vector{Float64}: Main diagonal (length n), unchanged on exit
• du::Vector{Float64}: Upper diagonal (length n-1), unchanged on exit
• x::Vector{Float64}: On entry: RHS vector b. On exit: solution x
• d_work::Vector{Float64}: Work array (length n) for modified diagonal

Notes

• This is allocation-free when d_work is pre-allocated
• Uses partial pivoting implicitly (Thomas algorithm is stable for diagonally dominant systems like the P1 radiation equation)
• The original diagonals dl, d, du are preserved

update_radiation_properties!(rad_cache::RadiationCache, mesh,
							 prop_cache::PropertyCache)

Update radiation-specific properties in the cache. This is the single method of this function — the transparent limit is handled exclusively by compute_radiation_diffusion_coefficient (large finite D, never 0).

Computes:

1. Cell-centered D_rad[i] = 1/(3α[i])
2. Face-centered Dradface[f] = harmonic mean weighted by distance

The harmonic mean is used for D_rad at faces because radiation transport follows resistance-in-series (like conduction). This ensures correct flux at interfaces between materials with different optical properties.

Arguments

• rad_cache: Radiation cache to update
• mesh: Mesh with geometry information
• prop_cache: Property cache with α values (must be updated first)

verify_coarsening_conservation(V_a, state_a, V_b, state_b, V_p, state_p, NC) -> NamedTuple

Verify that coarsening conserved mass and energy.

Arguments

• Pre-coarsening: Va, statea, Vb, stateb (children)
• Post-coarsening: Vp, statep (parent)
• NC: Number of components

Returns

• NamedTuple with fields: masserror, energyerror

verify_refinement_conservation(mesh::Mesh1D{NC}, cell_id::Int) -> NamedTuple

Verify that refinement conserved mass and energy.

Returns

• NamedTuple with fields: masserror, energyerror