\Large{\bfseries Integration Schemes for Simple Cooling}}

Note: This page is based on/inspired by An Exact Integration Scheme for Radiative Cooling in Hydrodynamical Simulations. Townsend, R.H.D. 2009, ApJS, 181, 391, to be attached to this page.

Background

% \section{Background} % In many of the astrophysical systems studied with $AstroBEAR$, an important piece of physics is the loss of thermal energy via optically-thin emission, a process known as $radiative\ cooling$. \\ \\ The ``optically thin'' qualifier is used to represent systems where the radiating energy can be assumed to simply leave the system. In contrast, ``optically thick'' systems, in which the photons interact with the system before exiting, require treatment known as $radiation\ transport$, a complicated bit of physics and code development which is not covered in this page. \\ \\ Optically thin radiative cooling is employed in many computational codes. Typically, codes treat the effects in an $operator\ split$ approach: one may write the evolution of the pressure $P$ as \begin{equation}\tag{1} \frac{dP}{dt}\bigg|_{total}=\frac{dP}{dt}\bigg|_{ad} + \frac{dP}{dt}\bigg|_{cool} \end{equation} where \begin{align} \frac{dP}{dt}\bigg|_{ad} &= \frac{\gamma P}{\rho} \frac{d\rho}{dt} \tag{1a}\ ,\qquad \text{and}\\ \frac{dP}{dt}\bigg|_{cool} &= -(\gamma-1)n_e n_H \Lambda(T)\ . \tag{1b} \end{align} In the above, $t$ is time, $\gamma$ is ratio of specific heats, $n_e$ and $n_H$ number densities of electrons and hydrogen atoms, respectively, $\Lambda(T)$ the cooling function, and $\rho$ mass density ($\rho\equiv \mu n$, where $\mu$ is mean molecular weight and $n$ number density). (Note that for an ideal gas, $P=n k_B T=\rho\mu^{\text{-}1} k_B T$ where $k_B$ is Boltzmann's constant.) \\ \\ The exact treatment or approximation of the cooling function $\Lambda(T)$ determines what methods may be used to evaluate (1b) above.

Simplest form: power law

Main wikipage: TestSuite/RadiativeShocks.

% \subsection{Simplest form: power law} % The simplest form of the cooling function is that of a power law, \begin{align} \Lambda(T) = \alpha T^\beta\ . \tag{2} \end{align} For certain regimes this may be exact, for example temperatures above $10^8$ K, where the only contributor is assumed to be Bremsstrahlung radiation. \\ \\ However, for temperatures less than this, one would not expect the cooling function to be a simple function of temperature. For example, the cooling curve of Dalgarno \& McCray (DM; 1972), based on assumptions including solar abundances of elements, displays a fair amount of structure in the range $10^4\le T\le10^8$ K: %\begin{figure}[!htbp] %\plotone{dmcooling.png} %\caption{\label{}} %\end{figure}

his then implies that a more complicated functional form of $\Lambda(T)$ will be required in these temperature ranges. \\ \\ Note that temperatures below $10^4$ K are not handled well with the DM cooling curve. This results from the fact that molecular physics become increasingly important as the temperature drops enough for molecules to make up a large fraction of the gas. Molecular physics is not very well implemented in the DM curve and requires additional handling (see below).

Piecewise-power-law form

% \subsection{Piecewise-power-law form} % A common method to approximate curves like the above is to assume the cooling function may be represented by a series of power laws $\Lambda_k(T)$, e.g. \begin{align} \Lambda(T) = \Lambda_k\left(\frac{T}{T_k}\right)^{\beta_k} \qquad T_k\le T \le T_{k+1}\ . \tag{3} \end{align} This admits a well-defined cooling rate at any point in the range of cooling functions $k$, allowing direct calculation of $\Lambda_k(T)$ at any point in this range.

Piecewise-linear-reconstruction form

% \subsection{Piecewise-linear-reconstruction form} % As a final alternative, one may read in a table of discrete points of an arbitrary cooling function, \begin{align} \Lambda_i(T_i),\quad \{i=1\ldots N\}\ . \tag{4} \end{align} Then, linear interpolation may be employed to get values of the cooling function at temperatures located between $T_i$ and $T_{i+1}$ in the table.

Three Schemes for Solving the Cooling Equation

% \section{Three Schemes for Solving the Cooling Equation} % A semidiscrete, piecewise-constant finite-volume discretization of Eq.~(1b) for cell $i$ and timestep $\Delta t$ from time $n$ to time $n+1$ yields \begin{align} \frac{T_i^{n+1}-T_i^n}{\Delta t} = -\frac{(\gamma-1)\rho_i\mu}{k_B\mu_e\mu_H} \Lambda(T_i)\ . \tag{5} \end{align} where $\mu_e$ and $\mu_H$ are electron and hydrogen mean molecular weights. This equation is ``semidiscrete'' because the cooling function $\Lambda(T_i)$ is not yet discretized in time. How one chooses to do this dictates whether the discretization is $explicit$ or $implicit$. After discussing these two methods, we will revisit Eq.~(1b) and discuss the $exact\ integration$ (EI) method.

Explicit scheme

% \subsection{Explicit scheme} % One obtains the explicit scheme for solving Eq.~(5) by substituting the present temperature $T_i^n$ into the cooling function. Rearranging---and defining a parameter $t_{cool}$---we reach an explicit expression for the updated temperature $T_i^{n+1}$, to first order Eq.~(5) becomes \begin{align} T_i^{n+1} = T_i^n\left[1-\frac{\Delta t}{t_{cool}^n}\right]\ , \tag{6} \end{align} where \begin{align} t_{cool}^n = \frac{k_B \mu_e\mu_H T_i^n}{(\gamma-1)\rho_i\mu \Lambda(T_i^n)}\ . \tag{7} \end{align} One solves Eq.~(7) then by using one of the prescriptions given in Eq's (2)--(4) to get a value of $\Lambda(T_i^n)$ and directly calculate the new temperature $T_i^{n+1}$. \\ \\ Note that the above, being first-order, is diffuse. Higher-order discretizations in time and/or space (where the concept of order is based upon Taylor expansion) produce increasingly accurate results. However, each increase in order requires an additional calculation of $\Lambda(T)$, and hence results in multiples of the computational time of the first-order method. \\ \\ Note also that if $\Delta t\ge t_{cool}^n$, the above may perform poorly. This can be ameliorated by subcycling: reducing the timestep $\Delta t$ in order to successively solve Eq.~(7) $M$ times over timesteps $\Delta t_m\equiv \Delta t/M$. This strategy is successful because it essentially decreases the cooling strength, allowing the method to converge more quickly than the cooling timescale.

The explicit scheme in AstroBEAR

Click for a discussion of the way the explicit scheme is implemented in AstroBEAR

Aside: At present (11/2010), the source term algorithm does not distinguish between stiff and non-stiff source terms. A ``stiff'' source term is one in which the timescale associated with the source is much smaller than the hydrodynamic timescale, i.e. \begin{align} \Delta t_{src}\ll \Delta t_{hydro} \equiv \frac{c_s}{\Delta x}C \ , \tag{8} \end{align} where $c_s$ is the sound speed $(c_s\equiv \sqrt{\gamma P/\rho})$, $\Delta x$ the cell size, and $C$ the CFL number.] The algorithm in AstroBEAR treats the cooling term as a stiff source term. The basic algorithm is as follows below. (Note that this method is carried out individually on every cell in the computational domain, admitting a possible avenue of improvement involving some sort of binning and global search-and-replacement.) \begin{enumerate} \item Protect incoming values of q. This may not be necessary. \item Using initial q, get $\Delta q/\Delta t$, internal energy and sound speed to establish both the initial timestep $\Delta t^0$ and a quantity used to measure the error in the solution ``qScale'', i.e. \begin{align*} qScale=\begin{pmatrix} \rho \\ c_s \\ c_s \\ e \end{pmatrix}\ ,\quad \frac{\Delta q}{\Delta t} = \begin{pmatrix} 0 \\ 0 \\ 0 \\ -S_c\Lambda(T) \rho^2 \\ \end{pmatrix} \end{align*} where $S_c$ is a scale factor to convert to computational units, and $\Lambda(T)$ is calculated via linear interpolation of a table as in Eq.~(4) above. The initial timestep for cell $i$ at timestep $n$ is then \begin{align*} \Delta t_i^0 &= \text{MAX}\left( \text{MINVAL}(\left|\frac{qScale}{\Delta q/\Delta t}\right|), 10^{-8} \right) \\ &= \text{MAX}\left( \frac{e_i}{S_c \Lambda(T_i^n) \rho_i^2}, 10^{-8} \right) \end{align*} \item The code then enters a loop which updates $qScale$ above based on the previous iteration, makes a call to STIFF, and updates the timestep for the next iteration based on the result of STIFF. \end{enumerate}

o, what is STIFF? STIFF attempts to perform a 4th-order temporal integration of q over time step $\Delta t$. It performs 3 protections along the way; if any fails the timestep is multiplied by a factor of 0.10 and it tries again. If the protections do not fail, an error $\varepsilon$ is estimated and a relative error, given by \begin{align*} \varepsilon_{rel} = \frac{\varepsilon}{\text{qScale}\cdot\text{prec}} \end{align*} where ``prec'' is a precision tolerance set in physics.data, normally at $10^{-3}$. If $\varepsilon_{rel}$, the timestep $\Delta t$ is reduced and the algorithm cycles. If $\varepsilon_{rel}$, the $\Delta t$ is returned as well as a prediction for the next $\Delta t$. The STIFF algorithm then exits, so that a new qScale may be calculated before being called again.

Implicit scheme

onversely, instead of subcycling one may consider the implicit form of Eq.~(5), in which we substitute the updated temperature $T_i^{n+1}$ into the cooling function, which yields \begin{align} t_i^{n+1} &= T_i^n \left[ 1-\frac{\Delta t}{t_{cool}^{n+1}}\right] \notag\\ &= T_i^n \left[1-\frac{\Lambda(T_i^{n+1})}{\Lambda(T_i^n)}\frac{\Delta t}{t_{cool}^n} \right] \ ,\tag{9} \end{align} where we define the updated cooling time $t_{cool}^{n+1}$ analogously to Eq.~(7). This has the advantage of not depending upon the present cooling function and therefore not susceptible to the same deficiencies as the explicit Eq.~(6). \\ \\ As you can see, both the left and right sides of the above depend upon the updated temperature $T_i^{n+1}$, the left explicitly and the right $implicitly$. The equation therefore cannot be solved directly, instead requiring some manner of iterative solve to calculate the updated temperature. We briefly discuss one method below.

An implicit method example: Secant/bisection

% \subsubsection{An implicit method: Secant/bisection} % asdf