module radiative cooling – add optically thin radiative cooling More...

Data Types
interface	eos_scalar2_func
	Scalar EoS inverse, e.g. fleint_from_T(log_nH, log_T) More...

interface	get_2var_subr

interface	get_subr1

type	rc_fluid

Functions/Subroutines
subroutine	radiative_cooling_init_params (phys_gamma, he_abund)
	Radiative cooling initialization.

subroutine	voigt_escape_init_table ()
	Build the Voigt escape probability lookup table. Called once (guarded by voigt_table_ready flag). Uses 64-point Gauss-Legendre quadrature on [0, x_max] to evaluate E(tau_0) = (2 phi(0)/tau_0) * integral_0^{x_max} [1-exp(-tau_0 g(x))] dx where g(x) = phi(x)/phi(0) for the Voigt profile H(a,x) truncated at x_max.

double precision function	voigt_profile_ratio (a, x)
	Voigt profile ratio phi(x)/phi(0) using Humlicek (1982) Region I/II approx. For the small-a regime (a < 0.01), this simplifies to: H(a,x)/H(a,0) ≈ exp(-x^2) + a*sqrt(pi)/x^2 for \|x\| > ~2 We use the exact Gaussian core + Lorentzian wing decomposition.

double precision function	voigt_escape_lookup (tau)
	Look up the Voigt escape probability for a given tau. Uses linear interpolation in log10(tau) space.

subroutine	voigt_gauss_legendre (a, b, n, x, w)
	Gauss-Legendre quadrature nodes and weights on [a,b]. Uses the Golub-Welsch algorithm for n points.

subroutine	radiative_cooling_init (fl, read_params)

subroutine	create_y_ppl (fl)

subroutine	getvar_cooling (ixil, ixol, w, x, coolrate, fl)

subroutine	getvar_cooling_exact (qdt, ixil, ixol, wct, w, x, coolrate, fl)

subroutine	radiative_cooling_add_source (qdt, ixil, ixol, wct, wctprim, w, x, qsourcesplit, active, fl)

subroutine	floortemperature (qdt, ixil, ixol, wct, w, x, fl)

subroutine	radiative_cooling_taper (ixd, x_ndim, rho_val, te_val, fl, factor)

subroutine	get_cool_equi (qdt, ixil, ixol, wct, w, x, fl, res)

subroutine	cool_exact (qdt, ixil, ixol, wct, wctprim, w, x, fl)

subroutine	calc_l_extended (tpoint, lpoint, fl)

double precision function	lowfip_fraction (tpoint, fl)

subroutine	findl (tpoint, lpoint, fl)

subroutine	findy (tpoint, ypoint, fl)

subroutine	findt (tpoint, ypoint, fl)

subroutine	finddldt (tpoint, dlpoint, fl)

subroutine	build_y_mod_table (fl)
	===================================================================

double precision function	simpson_composite (f, n, h)
	Composite Simpson's rule on (N+1) equally spaced samples (N even). N must be a positive even integer; h is the step size.

double precision function	boole_composite (f, n, h)
	Composite Boole's rule on (N+1) equally spaced samples (N a multiple of 4). Each 4-step block contributes (2h/45)*(7 f0 + 32 f1 + 12 f2 + 32 f3 + 7 f4).

double precision function	invert_row_bisect (y_row, t_grid, ncool, y_target)
	Bisection on a single Y_mod row to find log10(T) such that Y_mod_row(i) = y_target. The row is monotonically increasing in decreasing i (since Y(ncool)=0 grows toward Y(1)=Y_max). Returns log10 T (code units).

double precision function	findy_mod (te_loc, nh_loc, fl)
	Forward Y_mod lookup: bilinear interpolation in (log10 nH, log10 T) on the precomputed Y_mod table. Both axes are uniform in log space.

double precision function	findt_mod (y_target, nh_loc, fl)
	Inverse Y_mod lookup: given Y_target and nH, return T such that Y_mod(log10 nH, log10 T) = Y_target. Bisection on the cooling-table i index using values interpolated linearly between the two adjacent nH rows (same convention as the classical findT). Saturates at the table extremes when Y_target falls outside [Y(tcoolmax), Y(tcoolmin)].

Variables
double precision, public	lb_cool_accum = 0.0d0
	Per-rank cooling-only compute accumulator for lb_diagnose. Sums the wall time spent inside radiative_cooling_add_source across all blocks and substages within one full advance() call. Reset by mod_advance::advance at the start of each step. Note: this is a subset of lb_compute_accum (which already encloses the whole advect1 block loop); it isolates the cooling Newton/exact-solver kernel cost from the surrounding finite-volume work.

Detailed Description

module radiative cooling – add optically thin radiative cooling

only uses the (Townsend) exact integration method, can be used in HD, ffhd, MHD, twofl

Assumptions: full ionized plasma dominated by H and He, ionization equilibrium Formula: Q=-n_H*n_e*f(T), positive f(T) function is pre-computed and tabulated or a piecewise power law Uses the various cooling tables stored in mod_radloss_tables.t

Function/Subroutine Documentation

◆ boole_composite()

double precision function mod_radiative_cooling::boole_composite	(	double precision, dimension(0:n), intent(in)	f,
		integer, intent(in)	n,
		double precision, intent(in)	h
	)

Composite Boole's rule on (N+1) equally spaced samples (N a multiple of 4). Each 4-step block contributes (2h/45)*(7 f0 + 32 f1 + 12 f2 + 32 f3 + 7 f4).

Definition at line 1970 of file mod_radiative_cooling.t.

◆ build_y_mod_table()

subroutine mod_radiative_cooling::build_y_mod_table ( type(rc_fluid), intent(inout) fl )

===================================================================

Variable-c_V Townsend extension (Y_mod)

The original Townsend (2009) exact integration scheme assumes a constant heat capacity c_V = ρR/(γ-1). For LTE plasmas with H/He ionisation, c_V is a strong function of temperature in the recombination zone (Ibañez 1985, 1992 thermal-instability buffering).

The derivation defines a modified TEF Ỹ(T; ρ) ≡ ∫_T^{T_ref} c_V(T'; ρ) / (n_H n_e(T'; ρ) Λ(T')) dT' with c_V the LTE volumetric heat capacity including the ionisation-energy reservoir. The integral is recast via the change of variables u = e_int/n_H so that the discrete construction never finite-differences c_V — both integrand and energy update are driven by the same eint_from_T table, guaranteeing bit-consistent energy conservation.

Must be called after eos_finalise() so that eoseint_from_T, eosT (forward), and eosneOnH are all in code units. The hook is bind_eos_to_source() in mod_hd_eos.t / mod_mhd_eos.t.

Algorithm: change of variables from T to u = e_int/n_H. For each log10 nH grid point j (taken from the eoseint_from_T table's nH axis):

Cache u_i = e_int/n_H at every cooling-curve T_i.
Walk i = ncool-1 down to 1 and accumulate Y_mod(j, i) = Y_mod(j, i+1) + ∫_{u_i}^{u_{i+1}} du / (n_e Λ) using a composite Simpson (3-point, O(h^4)) or Boole (5-point, O(h^6)) rule with N_sub sub-intervals.
The integrand evaluates n_e via y_from_nH_eint and Λ via findL (with T from T_from_nH_eint).

Per-row inverse table T_mod_inv(j, k) is also built for the 'table' inverse method (alternative to bisection).

Definition at line 1795 of file mod_radiative_cooling.t.

Here is the call graph for this function:

◆ calc_l_extended()

subroutine mod_radiative_cooling::calc_l_extended	(	double precision, intent(in)	tpoint,
		double precision, intent(out)	lpoint,
		type(rc_fluid), intent(in)	fl
	)

Definition at line 1567 of file mod_radiative_cooling.t.

◆ cool_exact()

subroutine mod_radiative_cooling::cool_exact	(	double precision, intent(in)	qdt,
		integer, intent(in)	ixi,
		integer, intent(in)	l,
		integer, intent(in)	ixo,
			l,
		double precision, dimension(ixi^s,1:nw), intent(in)	wct,
		double precision, dimension(ixi^s,1:nw), intent(in)	wctprim,
		double precision, dimension(ixi^s,1:nw), intent(inout)	w,
		double precision, dimension(ixi^s,1:ndim), intent(in)	x,
		type(rc_fluid), intent(in)	fl
	)

Always compute the CLASSICAL Townsend advance first. This gives Tlocal2 using the same path as ionE=false. If the resulting dT is small (saturated y region), the classical kinetic form is exact and we're done. If dT is large (recombination zone), we redo the advance with the Y_mod path to capture variable-c_V dynamics from ionisation buffering.

Upgrade to Y_mod only in the recombination zone (large dT). Guard: if findY_mod overflows (cooling curve near-zero at this T), keep the classical result – variable c_V is irrelevant where Lambda ~ 0.

Recombination zone: use table-based de for variable c_V.

Saturated y or non-ionE: classical Townsend kinetic form.

Definition at line 1381 of file mod_radiative_cooling.t.

Here is the call graph for this function:

◆ create_y_ppl()

subroutine mod_radiative_cooling::create_y_ppl ( type(rc_fluid) fl )

Definition at line 964 of file mod_radiative_cooling.t.

◆ finddldt()

subroutine mod_radiative_cooling::finddldt	(	double precision, intent(in)	tpoint,
		double precision, intent(out)	dlpoint,
		type(rc_fluid), intent(in)	fl
	)

Definition at line 1719 of file mod_radiative_cooling.t.

◆ findl()

subroutine mod_radiative_cooling::findl	(	double precision, intent(in)	tpoint,
		double precision, intent(out)	lpoint,
		type(rc_fluid), intent(in)	fl
	)

Definition at line 1609 of file mod_radiative_cooling.t.

◆ findt()

subroutine mod_radiative_cooling::findt	(	double precision, intent(out)	tpoint,
		double precision, intent(in)	ypoint,
		type(rc_fluid), intent(in)	fl
	)

Definition at line 1672 of file mod_radiative_cooling.t.

◆ findt_mod()

double precision function mod_radiative_cooling::findt_mod	(	double precision, intent(in)	y_target,
		double precision, intent(in)	nh_loc,
		type(rc_fluid), intent(in)	fl
	)

Inverse Y_mod lookup: given Y_target and nH, return T such that Y_mod(log10 nH, log10 T) = Y_target. Bisection on the cooling-table i index using values interpolated linearly between the two adjacent nH rows (same convention as the classical findT). Saturates at the table extremes when Y_target falls outside [Y(tcoolmax), Y(tcoolmin)].

Definition at line 2067 of file mod_radiative_cooling.t.

◆ findy()

subroutine mod_radiative_cooling::findy	(	double precision, intent(in)	tpoint,
		double precision, intent(out)	ypoint,
		type(rc_fluid), intent(in)	fl
	)

Definition at line 1634 of file mod_radiative_cooling.t.

◆ findy_mod()

double precision function mod_radiative_cooling::findy_mod	(	double precision, intent(in)	te_loc,
		double precision, intent(in)	nh_loc,
		type(rc_fluid), intent(in)	fl
	)

Forward Y_mod lookup: bilinear interpolation in (log10 nH, log10 T) on the precomputed Y_mod table. Both axes are uniform in log space.

Definition at line 2032 of file mod_radiative_cooling.t.

◆ floortemperature()

subroutine mod_radiative_cooling::floortemperature	(	double precision, intent(in)	qdt,
		integer, intent(in)	ixi,
		integer, intent(in)	l,
		integer, intent(in)	ixo,
			l,
		double precision, dimension(ixi^s,1:nw), intent(in)	wct,
		double precision, dimension(ixi^s,1:nw), intent(inout)	w,
		double precision, dimension(ixi^s,1:ndim), intent(in)	x,
		type(rc_fluid), intent(in)	fl
	)

Parameters

[in]	ixi	This will need revisiting in LTE
[in]	l	This will need revisiting in LTE
[in]	ixo	This will need revisiting in LTE
[in]	l	This will need revisiting in LTE

Definition at line 1170 of file mod_radiative_cooling.t.

◆ get_cool_equi()

subroutine mod_radiative_cooling::get_cool_equi	(	double precision, intent(in)	qdt,
		integer, intent(in)	ixi,
		integer, intent(in)	l,
		integer, intent(in)	ixo,
			l,
		double precision, dimension(ixi^s,1:nw), intent(in)	wct,
		double precision, dimension(ixi^s,1:nw), intent(inout)	w,
		double precision, dimension(ixi^s,1:ndim), intent(in)	x,
		type(rc_fluid), intent(in)	fl,
		double precision, dimension(ixi^s), intent(out)	res
	)

Always classical Townsend first. Upgrade to Y_mod only in the recombination zone (large ΔT) — saturated y uses the identical formula as ionE=false, removing per-substep asymmetry.

Definition at line 1262 of file mod_radiative_cooling.t.

Here is the call graph for this function:

◆ getvar_cooling()

subroutine mod_radiative_cooling::getvar_cooling	(	integer, intent(in)	ixi,
		integer, intent(in)	l,
		integer, intent(in)	ixo,
			l,
		double precision, dimension(ixi^s,1:nw)	w,
		double precision, dimension(ixi^s,1:ndim), intent(in)	x,
		double precision, dimension(ixi^s), intent(out)	coolrate,
		type(rc_fluid), intent(in)	fl
	)

Definition at line 988 of file mod_radiative_cooling.t.

Here is the call graph for this function:

◆ getvar_cooling_exact()

subroutine mod_radiative_cooling::getvar_cooling_exact	(	double precision, intent(in)	qdt,
		integer, intent(in)	ixi,
		integer, intent(in)	l,
		integer, intent(in)	ixo,
			l,
		double precision, dimension(ixi^s, 1:nw), intent(in)	wct,
		double precision, dimension(ixi^s, 1:nw)	w,
		double precision, dimension(ixi^s, 1:ndim), intent(in)	x,
		double precision, dimension(ixi^s), intent(out)	coolrate,
		type(rc_fluid), intent(in)	fl
	)

Always classical Townsend first. Upgrade to Y_mod only where ionisation buffering matters (large ΔT, recombination zone).

Recombination zone: table-based de for variable c_V.

Saturated y or non-ionE: identical Townsend kinetic form.

Definition at line 1034 of file mod_radiative_cooling.t.

Here is the call graph for this function:

◆ invert_row_bisect()

double precision function mod_radiative_cooling::invert_row_bisect	(	double precision, dimension(ncool), intent(in)	y_row,
		double precision, dimension(ncool), intent(in)	t_grid,
		integer, intent(in)	ncool,
		double precision, intent(in)	y_target
	)

Bisection on a single Y_mod row to find log10(T) such that Y_mod_row(i) = y_target. The row is monotonically increasing in decreasing i (since Y(ncool)=0 grows toward Y(1)=Y_max). Returns log10 T (code units).

Definition at line 1990 of file mod_radiative_cooling.t.

◆ lowfip_fraction()

double precision function mod_radiative_cooling::lowfip_fraction	(	double precision, intent(in)	tpoint,
		type(rc_fluid), intent(in)	fl
	)

Definition at line 1582 of file mod_radiative_cooling.t.

Here is the call graph for this function:

◆ radiative_cooling_add_source()

subroutine mod_radiative_cooling::radiative_cooling_add_source	(	double precision, intent(in)	qdt,
		integer, intent(in)	ixi,
		integer, intent(in)	l,
		integer, intent(in)	ixo,
			l,
		double precision, dimension(ixi^s,1:nw), intent(in)	wct,
		double precision, dimension(ixi^s,1:nw), intent(in)	wctprim,
		double precision, dimension(ixi^s,1:nw), intent(inout)	w,
		double precision, dimension(ixi^s,1:ndim), intent(in)	x,
		logical, intent(in)	qsourcesplit,
		logical, intent(inout)	active,
		type(rc_fluid), intent(in)	fl
	)

Definition at line 1142 of file mod_radiative_cooling.t.

Here is the call graph for this function:

◆ radiative_cooling_init()

subroutine mod_radiative_cooling::radiative_cooling_init	(	type(rc_fluid), intent(inout)	fl,
		external subroutine(type(rc_fluid), intent(inout) fl)	read_params
	)

Definition at line 446 of file mod_radiative_cooling.t.

Here is the call graph for this function:

◆ radiative_cooling_init_params()

subroutine mod_radiative_cooling::radiative_cooling_init_params	(	double precision, intent(in)	phys_gamma,
		double precision, intent(in)	he_abund
	)

Radiative cooling initialization.

Definition at line 262 of file mod_radiative_cooling.t.

◆ radiative_cooling_taper()

subroutine mod_radiative_cooling::radiative_cooling_taper	(	integer, intent(in)	ix,
		integer, intent(in)	d,
		double precision, intent(in)	x_ndim,
		double precision, intent(in)	rho_val,
		double precision, intent(in)	te_val,
		type(rc_fluid), intent(in)	fl,
		double precision, intent(out)	factor
	)

Parameters

[in]	ix	Compute multiplicative taper factor for radiative cooling. Returns 1.0 when no tapering applies; < 1.0 near boundaries, at high density, or in optically thick regions (escape probability).
[in]	d	Compute multiplicative taper factor for radiative cooling. Returns 1.0 when no tapering applies; < 1.0 near boundaries, at high density, or in optically thick regions (escape probability).

Definition at line 1191 of file mod_radiative_cooling.t.

Here is the call graph for this function:

◆ simpson_composite()

double precision function mod_radiative_cooling::simpson_composite	(	double precision, dimension(0:n), intent(in)	f,
		integer, intent(in)	n,
		double precision, intent(in)	h
	)

Composite Simpson's rule on (N+1) equally spaced samples (N even). N must be a positive even integer; h is the step size.

Definition at line 1953 of file mod_radiative_cooling.t.

◆ voigt_escape_init_table()

subroutine mod_radiative_cooling::voigt_escape_init_table

Build the Voigt escape probability lookup table. Called once (guarded by voigt_table_ready flag). Uses 64-point Gauss-Legendre quadrature on [0, x_max] to evaluate E(tau_0) = (2 phi(0)/tau_0) * integral_0^{x_max} [1-exp(-tau_0 g(x))] dx where g(x) = phi(x)/phi(0) for the Voigt profile H(a,x) truncated at x_max.

Definition at line 275 of file mod_radiative_cooling.t.

Here is the call graph for this function:

◆ voigt_escape_lookup()

double precision function mod_radiative_cooling::voigt_escape_lookup ( double precision, intent(in) tau )

Look up the Voigt escape probability for a given tau. Uses linear interpolation in log10(tau) space.

Definition at line 365 of file mod_radiative_cooling.t.

Here is the call graph for this function:

◆ voigt_gauss_legendre()

subroutine mod_radiative_cooling::voigt_gauss_legendre	(	double precision, intent(in)	a,
		double precision, intent(in)	b,
		integer, intent(in)	n,
		double precision, dimension(n), intent(out)	x,
		double precision, dimension(n), intent(out)	w
	)

Gauss-Legendre quadrature nodes and weights on [a,b]. Uses the Golub-Welsch algorithm for n points.

Definition at line 403 of file mod_radiative_cooling.t.

◆ voigt_profile_ratio()

double precision function mod_radiative_cooling::voigt_profile_ratio	(	double precision, intent(in)	a,
		double precision, intent(in)	x
	)

Voigt profile ratio phi(x)/phi(0) using Humlicek (1982) Region I/II approx. For the small-a regime (a < 0.01), this simplifies to: H(a,x)/H(a,0) ≈ exp(-x^2) + a*sqrt(pi)/x^2 for |x| > ~2 We use the exact Gaussian core + Lorentzian wing decomposition.

Definition at line 339 of file mod_radiative_cooling.t.

Here is the call graph for this function:

Variable Documentation

◆ lb_cool_accum

double precision, public mod_radiative_cooling::lb_cool_accum = 0.0d0

Per-rank cooling-only compute accumulator for lb_diagnose. Sums the wall time spent inside radiative_cooling_add_source across all blocks and substages within one full advance() call. Reset by mod_advance::advance at the start of each step. Note: this is a subset of lb_compute_accum (which already encloses the whole advect1 block loop); it isolates the cooling Newton/exact-solver kernel cost from the surrounding finite-volume work.

Definition at line 24 of file mod_radiative_cooling.t.

Data Types

Functions/Subroutines

Variables

Detailed Description

Function/Subroutine Documentation

◆ boole_composite()

◆ build_y_mod_table()

Variable-c_V Townsend extension (Y_mod)

◆ calc_l_extended()

◆ cool_exact()

◆ create_y_ppl()

◆ finddldt()

◆ findl()

◆ findt()

◆ findt_mod()

◆ findy()

◆ findy_mod()

◆ floortemperature()

◆ get_cool_equi()

◆ getvar_cooling()

◆ getvar_cooling_exact()

◆ invert_row_bisect()

◆ lowfip_fraction()

◆ radiative_cooling_add_source()

◆ radiative_cooling_init()

◆ radiative_cooling_init_params()

◆ radiative_cooling_taper()

◆ simpson_composite()

◆ voigt_escape_init_table()

◆ voigt_escape_lookup()

◆ voigt_gauss_legendre()

◆ voigt_profile_ratio()

Variable Documentation

◆ lb_cool_accum