Physics of Radiation Transfer

Spectral Intensity

Spectral Intensity

Typically when we discuss the radiation field we use the spectral intensity \left ( \nu, \mathbf{x}, \Omega \right ) which is a function of frequency, position, and direction. This is very similar to the phase space density used in deriving the fluid equations \left ( \mathbf{x}, \mathbf{v} \right ) except that

• light always travels at \(c\), so the velocity dependence is just a direction dependence.
• Furthermore, photons can have different frequencies, so there is an extra dimension to the phase space.
• Instead of storing the phase space density of photons, the spectral intensity is the phase space density of energy flux…

Going between photon number and energy just involves a factor of \( h \nu \) and going from energy density to energy flux density just involves a factor of \(c\) so we have

I \left ( \nu, \mathbf{x}, \Omega, \right ) = h \nu c f \left ( \nu, \mathbf{x}, \Omega, \right )

This can also be seen by considering the differential energy

E = I \left ( \nu, \mathbf{x}, \Omega, \right ) d\nu d\Omega dA dt = h \nu f \left ( \nu, \mathbf{x}, \Omega, \right ) d\nu d\Omega d

where the number of photons traveling normal to the surface dA that cross the surface dA in time dt is just the number of photons in the volume dV = dA c dt (assuming the photons are headed normal to dA)…

so we also have

E = h \nu f \left ( \nu, \mathbf{x}, \Omega, \right ) d\nu d\Omega dA c d

which gives

\left ( \nu, \mathbf{x}, \Omega, \right ) = h \nu c f \left ( \nu, \mathbf{x}, \Omega, \right )

Deriving the Transport Equation

Deriving the Transport Equation

If we consider the Boltzmann transport equation for photons of a specific frequency \(f_\nu\) we have

frac{\partial}{\partial t} f_\nu + \mathbf{v} \cdot \nabla f + \mathbf{F} \cdot \nabla_p f = \left ( \frac{\partial f}{\partial t} \right )_{coll}

Now photons don't experience body forces, always travel at the speed of light, and in general the "collision term" consists of photon emission and absorptions… so we have

frac{\partial}{\partial t} f_\nu + c \mathbf{n} \cdot \nabla f_\nu = A_{\nu} - \chi_\nu f_\nu c

where \(A_{\nu} \) is the emission rate of photons of frequency \(\nu\) and the mean free path length is given by \(\chi_\nu = \sigma_\nu n\) where \(\sigma_\nu\) is the particle scattering cross section and \(n\) is the number density of particles.

Now if we multiply through by \( h\nu\) we have

frac{\partial}{c \partial t} I_\nu + \mathbf{n} \cdot \nabla I_\nu = \eta_\nu - \sigma_\nu I_\nu

where \(\eta_\nu = h\nu A_{\nu}\) is the radiative power

If we solve the transport equation along a characteristic

left [ \mathbf{x} \left ( s \right ) , t \left ( s \right ) \right ] = \left [ \mathbf{x0} + \mathbf{n} s, \frac{s}{c} \right ]

we have

frac{d I_\nu}{d s} = \frac{\partial I_\nu}{\partial x^i} \frac{\partial{x^i}}{\partial s} + \frac{\partial I_\nu}{\partial t}\frac{\partial t}{\partial s} = \mathbf{n} \cdot \nabla I_\nu + \frac{1}{c}\frac{\partial I_\nu}{\partial t} = \eta_\nu(s) - \chi_\nu(s) I_\nu (s)

where \(f(s) = f(\mathbf{x}(s), t(s)) = f \left ( \mathbf{x_0}+\mathbf{n} s, \frac{s}{c} \right ) \) and then we can divide through by \(\chi_\nu(s)\) we get

frac{dI_\nu}{\chi_\nu(s) ds} = \frac{\eta_\nu(s)}{\chi_\nu(s)} - I_\nu(s) = S_\nu(s) - I_\nu(s

Now if we define \(d\tau_\nu = \chi_\nu(s) ds \) which gives

tau_\nu(s) = \int\limits_0^s \chi_\nu(s') ds

and

(\tau_\nu) = \int\limits_0^\tau_\nu \frac{1}{\chi_\nu} d\tau'_\nu

we can write the transport equation in the simplest form

frac{dI_\nu}{d\tau_\nu} = S_\nu(\tau_\nu) - I_\nu(\tau_\nu

although the RHS is now more difficult to evaluate as

f \left ( \tau_\nu \right ) = f \left ( s \left ( \tau_\nu \right ) \right ) = f \left ( \mathbf{x} \left (s \left ( \tau_\nu \right ) \right ), t \left ( s \left ( \tau_\nu \right ) \right ) \right )

Also if we include scattering then the source function can depend on the mean radiative flux \(\frac{cE}{4 \pi}\) and the transport equation becomes an integro-differential equation that must be solved iteratively…

There are also a few important dimensionless numbers to consider:

tau = l \kappa=\frac{l}{\lambda_p	beta = \frac{u}{c
tau << 1	streaming limit
tau >> 1 \mbox{, } \beta \tau <<	static diffusion limit
tau >> 1 \mbox{, } \beta \tau >>	dynamic diffusion limit

Equations of Radiation Hydrodynamics

Equations of Radiation Hydrodynamics

Some of what follows is taken from ​Krumholz et al. 2007

frac{\partial \rho}{\partial t} + \nabla \cdot \left ( \rho \mathbf{v} \right ) =

frac{\partial }{\partial t} \left ( \rho \mathbf{v} \right ) + \nabla \cdot \left ( \rho \mathbf{vv} \right ) = -\nabla P+\mathbf{G

frac{\partial e}{\partial t} + \nabla \cdot \left [ \left ( e + P \right ) \mathbf{v} \right ] = c G^

frac{\partial E}{\partial t} + \nabla \cdot \mathbf{F} = -c G_0

frac{1}{c^2} \frac{\partial \mathbf{F}}{\partial t} + \nabla \cdot \mathbf{P} = -\mathbf{G

where the moments of the specific intensity are defined as

Radiation Energy moments	Corresponding fluid moments
E =\int_0^\infty d \nu \int d \Omega I(\mathbf{n}, \nu	rho=\int d \mathbf{v} f(\mathbf{v}
mathbf{F}=\int_0^\infty d \nu \int d \Omega \mathbf{n} I(\mathbf{n}, \nu	rho \mathbf{v} =\int d \mathbf{v} \mathbf{v} f(\mathbf{v}
\mathbf{P}=\int_0^\infty d \nu \int d \Omega \mathbf{nn} I(\mathbf{n}, \nu	mathbf{P} =\int d \mathbf{v} \mathbf{vv} f(\mathbf{v}

and the radiation 4-force density is given by

G^0 = \int_0^\infty d \nu \int d \Omega \left [ \kappa (\mathbf{n}, v)I(\mathbf{n}, \nu)-\eta(\mathbf{n},v) \right

\mathbf{G} = \int_0^\infty d \nu \int d \Omega \left [ \kappa (\mathbf{n}, v)I(\mathbf{n}, \nu)-\eta(\mathbf{n},v) \right ] \mathbf{n

If we had a closure relation for the radiation pressure then we could solve this system. For gas particles, collisions tend to produce a Boltzmann Distribution which is isotropic and gives a pressure tensor that is a multiple of the identity tensor. Photons do not "collide" with each other and they all have the same velocity 'c' but in various directions. If the field were isotropic than \(P^{ij}=\delta{ij} 1/3 E\) but in general \(P^{ij}=f{ij} E\) where 'f' is the Eddington Tensor.

Simplifiying assumptions

Simplifying assumptions

• If the flux spectrum of the radiation is direction-independent then we can write the radiation four-force density in terms of the moments of the radiation field

^0=\gamma [ \gamma^2 \kappa_{0E} + (1-\gamma^2)\kappa_{0F}]E-\gamma \kappa_{0P}a_RT_0^4-\gamma(\mathbf{v} \cdot \mathbf{F}/c^2)[\kappa_{0F}-2\gamma^2(\kappa_{0F}-\kappa_{0E})] - \gamma^3(\kappa_{0F}-\kappa+{0E})(\mathbf{vv}) : \mathbf{P}/c^

mathbf{G} = \gamma \kappa_{0F}(\mathbf{F}/c)-\gamma \kappa_{0P}a_RT_0^4 (\mathbf{v}/c) - \left [ \gamma^3( \kappa_{0F} - \kappa_{0E})( \mathbf{v}/c)E + \gamma \kappa_{0F}(\mathbf{v}/c) \cdot \mathbf{P} \right ] + \gamma^3 ( \kappa_{0F} - \kappa_{0E}) \left [ 2 \mathbf{v} \cdot \mathbf{F}/c^3 - ( \mathbf{vv}) : \mathbf{P}/c^3 \right ] \mathbf{v}

where

kappa_{0P	comoving-frame Planck function weighted opacity
kappa_{0E	comoving-frame radiation energy weighted opacity
kappa_{0F	comoving-frame radiation flux weighted opacity

• If the radiation has a blackbody spectrum then \(\kappa_{0E}=\kappa_{0P}\)
• If the radiation is optically thick, then

mathbf{F_0}(\nu_0) \propto -\nabla E_0(\nu_0)/\kappa_0(\nu_0) \propto -[\partial B(\nu_0, T_0) / \partial T_0](\nabla T_0)/\kappa_0(\nu_0

which implies that kappa_{0F}^{-1}=\kappa_{0R}^{-1}=\frac{\int_0^\infty d \nu_0 \kappa_0(\nu_0)^{-1}[\partial B(\nu_0,T_0)/\partial T_0]}{\int_0^\infty d \nu_0 [\partial B(\nu_0, T_0)/\partial T_0]

• In the optically thin regime, \(|\mathbf{F}_0(\nu_0)| \rightarrow cE_0(\nu_0)\) so we would have \(\kappa_{0F}=\kappa_{0E}\) however assuming a blackbody temperature in the optically thin limit may be any more accurate than assuming that \(\kappa_{0F}=\kappa_{0R}\)

Flux limited diffusion

Flux limited diffusion

The flux limited diffusion approximation drops the radiation momentum equation in favor of \(\mathbf{F}_0=-\frac{c\lambda}{\kappa_{0R}}\nabla E_0\) where \(\lambda\) is the flux-limiter and is given by

lambda = \frac{1}{R} \left ( \coth R - \frac{1}{R} \right )

and

=\frac{|\nabla E_0|}{\kappa_{0R}E_0

which corresponds to a pressure tensor

mathbf{P}_0=\frac{E_0}{2}[(1-R_2)\mathbf{I}+(3R_2-1)\mathbf{n}_0\mathbf{n}_0

where

_2=\lambda+\lambda^2R^

Here are \(\lambda®\) and \(R_2®\) plotted

If we Lorentz boost the comoving terms into the lab frame and keep terms necessary to maintain accuracy we get:

^0=\kappa_{0P} \left ( E-\frac{4 \pi B}{c} \right ) + \left ( \frac{\lambda}{c} \right ) \left ( 2 \frac{\kappa_{0P}}{\kappa_{0R}} - 1 \right ) \mathbf{v} \cdot \nabla E - \frac{\kappa_{0P}}{c^2} E \left [ \frac{3-R_2}{2}v^2 + \frac{3R_2-1}{2}(\mathbf{v} \cdot \mathbf{n})^2 \right ] + \frac{1}{2} \left ( \frac{v}{c} \right ) ^2 \kappa_{0P} \left ( E - \frac{4 \pi B}{c} \right ) mathbf{G} = -\lambda \nabla E + \kappa_{0P} \frac{\mathbf{v}}{c} \left ( E - \frac{4 \pi B}{c} \right ) - \frac{1}{2} \left ( \frac{v}{c} \right ) ^2 \lambda \nabla E + 2 \lambda \left ( \frac{\kappa_{0P}}{\kappa_{0R}} - 1 \right ) \frac{(\mathbf{v} \cdot \nabla E )\mathbf{v}}{c^2

Numerics of Flux Limited Diffusion

If we plug the expressions for the radiation 4-momentum back into the gas equations and keep terms necessary to maintain accuracy we get:

frac{\partial }{\partial t} \left ( \rho \mathbf{v} \right ) + \nabla \cdot \left ( \rho \mathbf{vv} \right ) = -\nabla P\color{green}{-\lambda \nabla E frac{\partial e}{\partial t} + \nabla \cdot \left [ \left ( e + P \right ) \mathbf{v} \right ] = \color{red}{-\kappa_{0P}(4 \pi B-cE)} \color{green}{+\lambda \left ( 2 \frac{\kappa_{0P}}{\kappa_{0R}}-1 \right ) \mathbf{v} \cdot \nabla E} \color{blue}{-\frac{3-R_2}{2}\kappa_{0P}\frac{v^2}{c}E frac{\partial E}{\partial t} \color{red}{ - \nabla \cdot \frac{c\lambda}{\kappa_{0R}} \nabla E} = \color{red}{\kappa_{0P} (4 \pi B-cE)} \color{green}{-\lambda \left(2\frac{\kappa_{0P}}{\kappa_{0R}}-1 \right )\mathbf{v}\cdot \nabla E} \color{green}{-\nabla \cdot \left ( \frac{3-R_2}{2}\mathbf{v}E \right )} \color{blue}{+\frac{3-R_2}{2}\kappa_{0P}\frac{v^2}{c}E}

For static diffusion, the terms in blue with \(\frac{v^{2}{c}\) can be dropped and the system can be split into the usual hydro update (black), radiative source terms (green), and a coupled implicit solve (red) for the radiation energy density and thermal energy density (ie temperature).}

Operator Splitting 1

Operator splitting

Krumholz et al. perform Implicit Radiative, Explicit Hydro, Explicit Radiative

In AstroBEAR this would look like:

• Initialization
  • Prolongate, \(\rho\), \(\rho\mathbf{v}\), \(e\), \(E\), and \(\dot{E}^I\)
• Step 1
  • Overlap \(\rho\), \(\rho \mathbf{v} \) , \(e\), \(E\) and do physical BC's
  • Do IR which updates \(e_0\) and \(E_0\) using \(\rho_1\), \(e_1\), \(E_1\), and \(\dot{E}^I_1\)
  • Update \(E_{2\mbox{mbc}}\) using \(\dot{E}^{I_{2\mbox{mbc}}\)}
  • Update \(e_{2\mbox{mbc}}\) using \(E_{2\mbox{mbc}}\), \(\dot{E}^{I_{2\mbox{mbc}}\), and \(e_{2\mbox{mbc}}\)}
  • Update \(\dot{E}^{I_0\) using pre IR and post IR \(E_0\)}
  • Ghost \(e_{2\mbox{mbc}}\), \(E_{\mbox{mbc}+1}\), \(\dot{E}^{I_{\mbox{mbc}+1}\)}
  • Do first EH_mbc
  • Do ER_mbc —- Terms with \(\nabla E\) can be done without ghosting since EH did not change \(E\). The \(\nabla \cdot \mathbf{v}E\) term needs time centered face centered velocities which can be stored during the hydro update.
  • Store \(\dot{E}\) in child arrays to be prolongated
• Step 2
  • Overlap \(\rho\), \(\rho \mathbf{v} \) , \(e\), \(E\) and do physical BC's
  • Do IR which updates \(e_0\) and \(E_0\) using \(\rho_1\), \(e_1\), \(E_1\), and \(\dot{E}^I_1\)
  • Update \(\dot{E}^{I_0\) using pre IR and post IR \(E_0\)}
  • Update \(E_{1}\) using \(\dot{E}^I_{1}\)
  • Ghost \(e_{mbc}\), \(E_{1}\), \(\dot{E}^I_{1}\)
  • Do second EH₀
  • Do ER_mbc —- Terms with \(\nabla E\) can be done without ghosting since EH did not
  • Do ER₀ —- Terms with grad E can be done without ghosting since EH did not change \(E\). The \(\nabla \cdot \mathbf{v}E\) term needs time centered face centered velocities which can be stored during the hydro update.

Explicit Update 1

Explicit Update 1

The extra terms in the explicit update due to radiation energy are as follows:

frac{\partial }{\partial t} \left ( \rho \mathbf{v} \right ) =-\lambda \nabla frac{\partial e}{\partial t} = \lambda \left ( 2 \frac{\kappa_{0P}}{\kappa_{0R}}-1 \right ) \mathbf{v} \cdot \nabla E frac{\partial E}{\partial t} = -\lambda \left(2\frac{\kappa_{0P}}{\kappa_{0R}}-1 \right )\mathbf{v}\cdot \nabla E -\nabla \cdot \left ( \frac{3-R_2}{2}\mathbf{v}E \right )

These can be discretized as follows:

^{n+1}_i=p^n_i - \frac{1}{2}\frac{\Delta t}{\Delta x} \lambda_{i} \left ( E^n_{i+1}-E^n_{i-1} \right )

^{n+1}_i=e^n_i + \frac{1}{2}\frac{\Delta t}{\Delta x} \lambda_i \left ( 2 \frac{\kappa^n_{i,0P}}{\kappa^n_{i,0R}}-1 \right ) \left ( v^n_i \left ( E^n_{i+1}-E^n_{i-1} \right ) \right )

^{n+1}_i=E^n_i - \frac{\Delta t}{\Delta x} \left ( \frac{\lambda_i}{2} \left ( 2 \frac{\kappa^n_{i,0P}}{\kappa^n_{i,0R}}-1 \right ) \left ( v^n_i \left ( E^n_{i+1}-E^n_{i-1} \right ) \right ) + \left ( F_{i+1/2}-F_{i-1/2} \right ) \right

where _{i+1/2}=\frac{3-R_{2,i+1/2}}{8} \left(v_{i}^n+v_{i+1}^n \right ) \left ( E^n_i+E^n_{i+1} \right )

where

_{2,i+1/2} = \lambda_{i+1/2}+\lambda_{i+1/2}^2 R_{i+1/2}^ and _{i+1/2} = \frac{\left | E^n_{i+1}-E^n_{i} \right | }{2 \kappa_{0R,i+1/2} \left ( E^n_i+E^n_{i+1} \right )

and lambda_{i+1/2} = \frac{1}{R_{i+1/2}} \left ( \coth R_{i+1/2} - \frac{1}{R_{i+1/2}} \right )

and kappa_{0R,i+1/2} = \frac{\kappa^n_{0R,i}+\kappa^n_{0R,i+1}}{2} \mbox{ and } \lambda_{i+1/2} = \frac{1}{R_{i+1/2}} \left ( \coth R_{i+1/2} - \frac{1}{R_{i+1/2}} \right )

and

lambda_{i} = \frac{1}{R_{i}} \left ( \coth R_{i} - \frac{1}{R_{i}} \right )

and _{i} = \frac{\left | E^n_{i+1}-E^n_{i-1} \right | }{2 \kappa_{0R,i} E^n_{i}

Implicit Update 1

Implicit Update 1

For now we will assume that \(\kappa_{0P}\) and \(\kappa_{0R}\) are constant over the implicit update and we will treat the energy as the total internal energy ignoring kinetic and magnetic contributions. In this case we can solve the radiation energy equations:

frac{\partial E}{\partial t} = \nabla \cdot \frac{c\lambda}{\kappa_{0R}} \nabla E + \kappa_{0P} (4 \pi B(T)-cE) = \nabla \cdot \mathbf{F} + \kappa_{0P} (4 \pi B(T)-cE frac{\partial e}{\partial t} = - \kappa_{0P} (4 \pi B(T)-cE

where mathbf{F} = \frac{c\lambda}{\kappa_{0R}} \nabla

Expanding about e₀

Of course even if the opacity is independent of energy and radiation energy, the above combined system of equations is still non-linear due to the dependence on Temperature of the Planck Function \(B(T)\)

However we can expand the Plank Function about \(e_0\)

(T) = B \left ( T_0+dT \right ) = B \left ( T_0 \right ) + \left . \frac{\partial B}{\partial T} \right | _{T_0} \frac{\partial T}{\partial e} de = \frac{c}{4 \pi} a_R \left ( T_0^4 + 4T_0^3\Gamma de \right ) = B_0 \left ( 1 + 4\Gamma \frac{e-e_0}{T_0} \right )

where

Gamma = \frac{\partial T}{\partial e} = \frac{(\gamma-1)}{n k_B

Then the system of equations becomes

frac{\partial E}{\partial t} = \nabla \cdot \frac{c\lambda}{\kappa_{0R}} \nabla E + \kappa_{0P} \left [4 \pi B_0 \left ( 1 + 4\Gamma \frac{e-e_0}{T_0} \right )-cE \right ] frac{\partial e}{\partial t} = - \kappa_{0P} \left [ 4 \pi B_0 \left ( 1 + 4\Gamma \frac{e-e_0}{T_0} \right )-cE \right ]

which will be accurate as long as \(\left | 4\Gamma \frac{e-e_0}{T_0} \right | < \xi << 1\) or \(\Delta e = \left | e-e_0 \right | < \xi \frac{T_0}{4 \Gamma}\)

We can calculate the time scale for this to be true using the evolution equation for the energy density

Delta e = \Delta t \kappa_{0P} \left | 4 \pi B_0 -cE \right | < \xi \frac{T_0}{4 \Gamma

which gives Delta t < \xi \frac{T_0}{4 \Gamma \kappa_{0P} \left | 4 \pi B_0 - cE \right |

Implicit Discretization 1

We can now discretize the equations

^{n+1}_i-E^{n}_i = \left [ \alpha_{i+1/2} \left ( E^{*}_{i+1}-E^{*}_{i} \right ) - \alpha_{i-1/2} \left ( E^{*}_{i}-E^{*}_{i-1} \right ) \right ] - \epsilon E^{*}_i + \phi e^{*}_i + \theta_i - \phi_i e^n_i) ^{n+1}_i-e^{n}_i = \epsilon_i E^{*}_i - \phi_i e^{*}_i - \left(\theta_i-\phi_i e^n_i \right )

where the diffusion coefficient is given by

alpha_{i+1/2}=\frac{\Delta t}{\Delta x^2} \frac{c \lambda_{i+1/2}}{\kappa_{0R,i+1/2}} \mbox{ where } \kappa_{0R,i+1/2} = \frac{\kappa^n_{0R,i}+\kappa^n_{0R,i+1}}{2} \mbox{ and } \lambda_{i+1/2} = \frac{1}{R_{i+1/2}} \left ( \coth R_{i+1/2} - \frac{1}{R_{i+1/2}} \right ) and where _{i+1/2} = \frac{\left | E^n_{i+1}-E^n_{i} \right | }{2 \kappa_{0R,i+1/2} \left ( E^n_i+E^n_{i+1} \right )

and

epsilon_i=c\Delta t \kappa_{0P,i

represents the number of absorption/emissions during the time step

and

phi_i = \epsilon_i \frac{4 \pi}{c} B \left ( T^n_i \right ) \left ( \frac{4\Gamma}{T^n_i} \right )

theta_i = \epsilon_i \frac{4 \pi}{c} B \left ( T^n_i \right )

and we can think of the radiative flux as

frac{\Delta t}{\Delta x}\mathbf{F}_{i+1/2} = \alpha_{i+1/2} \left ( E^{*}_{i+1} - E^{*}_i \right )

Time Discretization

Now all the terms on the right hand side that are linear in E or e have been written as E^* or e^* because there are different ways to approximate E^* (e^*). For Backward Euler we have ^*_i = E^{n+1}_ and for Crank Nicholson we have ^*_i = \frac{1}{2} \left ( E^{n+1}_i + E^n_i \right ) or we can parameterize the solution ^*_i = \psi E^{n+1}_i + \bar{\psi}E^n_ where \(\bar{\psi} = 1-\psi\)

Backward Euler has \(\psi=1\) and Crank Nicholson has \(\psi=½\)

Forward Euler has \(\psi=0\)

In any event in 1D we have the following matrix coefficients

left [ 1 + \psi \left( \alpha_{i+1/2} + \alpha_{i-1/2} + \epsilon_i \right ) \right ] E^{n+1}_i - \left ( \psi \alpha_{i+1/2} \right ) E^{n+1}_{i+1} - \left ( \psi \alpha_{i-1/2} \right ) E^{n+1}_{i-1} - \left ( \psi \phi_i \right ) e^{n+1}_i=\left [ 1 - \bar{\psi} \left( \alpha_{i+1/2} + \alpha_{i-1/2} + \epsilon_i \right ) \right ] E^n_i + \left ( \bar{\psi} \alpha_{i+1/2} \right ) E^{n}_{i+1} + \left ( \bar{\psi} \alpha_{i-1/2} \right ) E^{n}_{i-1} - \psi\phi_i e^n_i + \theta_

left ( 1 +\psi \phi_i \right ) e^{n+1}_i - \left ( \psi \epsilon_i \right )E^{n+1}_i =\left ( 1 +\psi\phi_i \right ) e^n_i + \left ( \bar{\psi} \epsilon_i \right ) E^n_i-\theta_i

Now since the second equation has no spatial dependence, we can solve it for

color{purple}{e^{n+1}_i = e^n_i + \frac{1}{ 1 +\psi \phi_i}\left \{ \left ( \psi \epsilon_i \right )E^{n+1}_i + \left ( \bar{\psi} \epsilon_i \right ) E^n_i-\theta_i \right \}}

and plug the result into the first equation to get a matrix equation involving only one variable.

color{purple}{\left [ 1 + \psi \left( \alpha_{i+1/2} + \alpha_{i-1/2} + \frac{\epsilon_i}{ 1 +\psi \phi_i}\right ) \right ] E^{n+1}_i - \left ( \psi \alpha_{i+1/2} \right ) E^{n+1}_{i+1} - \left ( \psi \alpha_{i-1/2} \right ) E^{n+1}_{i-1} =\left [ 1 - \bar{\psi} \left( \alpha_{i+1/2} + \alpha_{i-1/2} +\frac{\epsilon_i }{ 1 +\psi \phi_i} \right ) \right ] E^n_i + \left ( \bar{\psi} \alpha_{i+1/2} \right ) E^{n}_{i+1} + \left ( \bar{\psi} \alpha_{i-1/2} \right ) E^{n}_{i-1} + \frac{\theta_i}{ 1 +\psi \phi_i}

This equation decouples and can be solved on it's own, and then the solution plugged back into the second equation to solve for the new energy.

2D etc…

For 2D or 3D we have more connections to add to the matrix elements but it is very straight forward… There will be additional alpha terms for each dimension, but everything else stays the same.

Initial solution vector

For the initial solution vector, we can just use Edot from the parent update (or last time step if we are on the coarse grid) to guess E, and then we can solve for the new e given our guess at the new E using the same time stepping (Backward Euler, Crank Nicholson, etc…).

^{n+1}_i = E^{n}_i+\dot{E}^{I}_i \Delta ^{n+1}_i = e^n_i + \frac{1}{1+\psi \phi_i} \left \{ \left ( \psi \epsilon_i \right )E^{n+1}_i + \left ( \bar{\psi} \epsilon_i \right ) E^n_i-\theta_i \right \}

Coarse-Fine Boundaries

Since we are doing our implicit solves first, we can use time interpolated solutions for the implicit solve for non-refined ghost zones. To do this we just need Edot. The opacities etc… in the ghost zones can be obtained from the hydro terms.

And the radiative implicit heating in coarse ghost cells can be done along with the initial solution vector so they are available for the hydro update.

Physical Boundary Conditions

Open (Free streaming) boundaries

We would like the radiation to leave at the free streaming limit. So frac{c \lambda}{\kappa_{0R}} \nabla E = \mathbf{F} = cE\mathbf{n} = \frac{c \lambda_g}{\kappa_g}\frac{ \left ( E-E_g \right ) }{\Delta x

Clearly if we set \(E_g = 0\) and \(\lambda_g =\kappa_g \Delta x\) we should get the correct flux.

This corresponds to an \(\alpha = c \frac{\Delta t}{\Delta x}\)

So we would just modify \(\alpha\) and zero out the matrix coefficient to the ghost zone

Constant Slope Boundary

Here we want the flux to be constant so energy does not pile up near the boundary. If we cancel all derivative terms on both sides of the cell, this will effectively match the incoming flux with the outgoing flux. This can also be done by setting \(\alpha_g = \alpha_i = 0\)

Periodic Boundary

This is the same as internal zones - it just maps the neighbor cell to be across the domain. Hypre has built in functionality for this under for the Struct Interface

User defined radiation field/Coarse Fine boundary

This will be the boundary at internal coarse-fine boundaries, but could also be used at the physical boundary if the radiation energy were specified.

Reflecting/ZeroSlope Boundary

Reflecting boundary should be fairly straightforward. This an be achieved by setting \(\alpha_g = 0\) which zeros out any flux - and has the same effect as setting \(E^{*}_g=E{*}_i\) or \(E^{n+1}_g=E{n+1}_i \mbox{ and } E^{n}_g=E{n}_i\)

Constant radiative flux

To have a constant radiative flux we must zero out terms involving the gradient and just add \(F_0 \frac{\Delta t}{\Delta x}\) in the source vector

Summary

Numerical value	Boundary	^{n+1}_{i+1	^{n+1}_	^{n+1}_{i-1
0	RAD_FREE_STREAMING		psi c \frac{\Delta t}{\Delta x	\psi\alpha_{i-1/2	\bar{\psi}c \frac{\Delta t}{\Delta x} E^n_i + \bar{\psi}\alpha_{i-1/2} \left ( E^n_{i-1}-E^n_i \right
1	RAD_EXTRAPOLATE_FLUX
2	INTERNAL/PERIODIC	\psi\alpha_{i+1/2	psi\alpha_{i+1/2	\psi\alpha_{i-1/2	bar{\psi}\alpha_{i+1/2} \left ( E^n_{i+1}-E^n_i \right ) + \bar{\psi}\alpha_{i-1/2} \left ( E^n_{i-1}-E^n_i \right
3	RAD_REFLECTING			\psi\alpha_{i-1/2	\bar{\psi}\alpha_{i-1/2} \left ( E^n_{i-1}-E^n_i \right
4	RAD_USERDEFINED_BOUNDARY/AMR_BOUNDARY		psi\alpha_{i+1/2	\psi\alpha_{i-1/2	bar{\psi}\alpha_{i+1/2} \left ( E^n_{i+1}-E^n_i \right ) - \psi \alpha_{i+1/2} E^{n+1}_{i+1} + \bar{\psi}\alpha_{i-1/2} \left ( E^n_{i-1}-E^n_i \right
5	RAD_USERDEFINED_FLUX			\psi\alpha_{i-1/2	_0 \frac{\Delta t}{\Delta x} + \bar{\psi}\alpha_{i-1/2} \left ( E^n_{i-1}-E^n_i \right

Alternative Splitting Method

While the previous method for splitting the equation technically works, it is likely to lead to large differences in the radiation field and the energy fields at coarse fine boundaries since there is no radiative flux to coarsen from the explicit updates to keep the various AMR levels consistent. There is one term that could potentially be coarsened and applied like any flux, but it is not the dominant term. As a result it might be better to treat the entire radiative energy equation implicitly (though keeping the velocity constant)

frac{\partial }{\partial t} \left ( \rho \mathbf{v} \right ) + \nabla \cdot \left ( \rho \mathbf{vv} \right ) = -\nabla P\color{green}{-\lambda \nabla E frac{\partial e}{\partial t} + \nabla \cdot \left [ \left ( e + P \right ) \mathbf{v} \right ] = \color{red}{-\kappa_{0P}(4 \pi B-cE)} \color{magenta}{+\lambda \left ( 2 \frac{\kappa_{0P}}{\kappa_{0R}}-1 \right ) \mathbf{v} \cdot \nabla E} \color{magenta}{-\frac{3-R_2}{2}\kappa_{0P}\frac{v^2}{c}E frac{\partial E}{\partial t} \color{red}{ - \nabla \cdot \frac{c\lambda}{\kappa_{0R}} \nabla E} = \color{red}{\kappa_{0P} (4 \pi B-cE)} \color{magenta}{-\lambda \left(2\frac{\kappa_{0P}}{\kappa_{0R}}-1 \right )\mathbf{v}\cdot \nabla E} \color{magenta}{-\nabla \cdot \left ( \frac{3-R_2}{2}\mathbf{v}E \right )} \color{magenta}{+\frac{3-R_2}{2}\kappa_{0P}\frac{v^2}{c}E}

It is possible that with included the terms in magenta in a semi-implicit method, the dynamic diffusion regime may be stable… In any event, it costs very little to add all of the terms in magenta to the implicit solve (using the old velocity). Then the momentum update can be done explicitly - though using a time centered radiation field.

Operator Splitting 2

Operator splitting 2

In AstroBEAR this would look like:

• Initialization
  • Prolongate, \(\rho\), \(\rho\mathbf{v}\), \(e\), \(E\), and \(\dot{E}\)
• Step 1
  • Overlap \(\rho\), \(\rho \mathbf{v} \) , \(e\), \(E\) and do physical BC's
  • Do IR which updates \(e_0\) and \(E_0\) using \(\rho_1\), \(e_1\), \(E_1\), and \(\dot{E}_1\)
  • Update \(E_{2\mbox{mbc}}\) using \(\dot{E}_{2\mbox{mbc}}\)
  • Update \(e_{2\mbox{mbc}}\) using \(E_{2\mbox{mbc}}\), \(\dot{E}_{2\mbox{mbc}}\), and \(e_{2\mbox{mbc}}\)
  • Update \(\dot{E}_0\) using pre IR and post IR \(E_0\)
  • Ghost \(e_{2\mbox{mbc}}\), \(E_{\mbox{mbc}+1}\), \(\dot{E}_{\mbox{mbc}+1}\)
  • Do first EH_mbc
  • Do ER_mbc —- Terms with \(\nabla E\) can be done without ghosting since EH did not change \(E\).
  • Store \(\dot{E}\) in child arrays to be prolongated
• Step 2
  • Overlap \(\rho\), \(\rho \mathbf{v} \) , \(e\), \(E\) and do physical BC's
  • Do IR which updates \(e_0\) and \(E_0\) using \(\rho_1\), \(e_1\), \(E_1\), and \(\dot{E}_1\)
  • Update \(\dot{E}_0\) using pre IR and post IR \(E_0\)
  • Update \(E_{1}\) using \(\dot{E}_{1}\)
  • Ghost \(e_{mbc}\), \(E_{1}\), \(\dot{E}_{1}\)
  • Do second EH₀
  • Do ER₀ —- Terms with \(\nabla E\) can be done without ghosting since EH did not change \(E\).

Explicit Update 2

Explicit Update 2

The extra terms in the explicit update due to radiation energy are as follows:

frac{\partial }{\partial t} \left ( \rho \mathbf{v} \right ) =-\lambda \nabla

These can be discretized as follows:

^{n+1}_i=p^n_i - \frac{1}{4}\frac{\Delta t}{\Delta x} \lambda_{i} \left ( \left ( E^{n}_{i+1}+E^{n+1}_{i+1} \right ) - \left ( E^n_{i-1}+E^{n+1}_{i-1} \right ) \right )

lambda_{i} = \frac{1}{R_{i}} \left ( \coth R_{i} - \frac{1}{R_{i}} \right )

and _{i} = \frac{\left | E^n_{i+1}-E^n_{i-1} \right | }{2 \kappa_{0R,i} E^n_{i}

Implicit Update 2

Implicit Update 2

For now we will assume that \(\kappa_{0P}\) and \(\kappa_{0R}\) are constant over the implicit update and we will treat the energy as the total internal energy ignoring kinetic and magnetic contributions. In this case we can solve the radiation energy equations:

frac{\partial e}{\partial t} = -\kappa_{0P}(4 \pi B-cE) + \lambda \left ( 2 \frac{\kappa_{0P}}{\kappa_{0R}}-1 \right ) \mathbf{v} \cdot \nabla E -\frac{3-R_2}{2}\kappa_{0P}\frac{v^2}{c}E frac{\partial E}{\partial t} = \nabla \cdot \frac{c\lambda}{\kappa_{0R}} \nabla E + \kappa_{0P} (4 \pi B(T)-cE) -\lambda \left(2\frac{\kappa_{0P}}{\kappa_{0R}}-1 \right )\mathbf{v}\cdot \nabla E + \frac{3-R_2}{2}\kappa_{0P}\frac{v^2}{c}E -\nabla \cdot \left ( \frac{3-R_2}{2}\mathbf{v}E \right )

which we can also write as

frac{\partial e}{\partial t} = f \left ( e \right ) + g \left ( E,\nabla E \right ) frac{\partial E}{\partial t} = \nabla \cdot \frac{c\lambda}{\kappa_{0R}} \nabla E -\nabla \cdot \left ( \frac{3-R_2}{2}\mathbf{v}E \right ) - f \left ( e \right ) - g \left (e,E, \nabla E \right ) where

f \left ( e \right ) = -4 \pi \kappa_{0P} B(T(e))

and

g \left (e,E, \nabla E \right ) = \kappa_{0P}(cE) + \lambda \left ( 2 \frac{\kappa_{0P}}{\kappa_{0R}}-1 \right ) \mathbf{v} \cdot \nabla E -\frac{3-R_2}{2}\kappa_{0P}\frac{v^2}{c}E

Now we can linearize f about e₀ f \left ( e \right ) = f \left ( e_0 \right ) + \frac{\partial f}{\partial e} \left ( e - e_0 \right )

so that the first equation can be written as

frac{\partial e}{\partial t} = f \left ( e_0 \right ) + g \left ( E,\nabla E \right ) + \frac{\partial f}{\partial e} \left ( e - e_0 \right )

and then discretized as

^{n+1}_i-e^{n}_i = \Delta t \left ( f \left ( e^n_i \right ) + g \left (E^*,\nabla E^* \right ) \right ) + \psi \phi e^n_i - \psi \phi e^{n+1}_i

where

phi = -\Delta t \frac{\partial f}{\partial e} = 4 \pi \kappa_{0P} \Delta t \frac{\partial B(T(e))}{\partial e

which can be solved for ^{n+1}_i = e^{n}_i + \frac{\Delta t}{1 + \psi \phi} \left ( f \left ( e^n_i \right ) + g \left (E^*,\nabla E^* \right ) \right )

or

^{n+1}_i = e^{n}_i + \frac{1}{1 + \psi \phi_i} \left ( -\theta_i + \epsilon_i E^*_i \omega_i v^n_{x,i} \left ( E^*_{i+1}-E^*_{i-1} \right ) - \xi_i E^* \right )

Then if we take the semi-discretized equation for E

frac{\partial E}{\partial t} = \nabla \cdot \frac{c \lambda}{\kappa_{0R}} \nabla E - \nabla \cdot \left ( \frac{3-R_2}{2} \mathbf{v} E \right ) - f \left ( e^n_i \right ) - g \left ( E, \nabla E \right ) - \frac{1}{\Delta t} \left ( \psi \phi_i e^n_i - \psi \phi_i e^{n+1}_i \right )

and then plugin the solution for eⁿ⁺¹_i we get

frac{\partial E}{\partial t} = \nabla \cdot \frac{c\lambda}{\kappa_{0R}} \nabla E -\nabla \cdot \left ( \frac{3-R_2}{2}\mathbf{v}E \right ) - f \left (e^n_i \right ) - g \left (E, \nabla E \right ) - \frac{1}{\Delta t} \left ( \psi \phi_i e^n_i - \psi \phi_i e^n_i - \frac{\psi \phi_i }{1+\psi \phi_i} \left ( f \left ( e^n_i \right ) + g \left (E,\nabla E \right ) \right ) \right )

which simplifies to

frac{\partial E}{\partial t} = \nabla \cdot \frac{c\lambda}{\kappa_{0R}} \nabla E -\nabla \cdot \left ( \frac{3-R_2}{2}\mathbf{v}E \right ) - \frac{1}{1+\psi \phi_i} \left ( f \left ( e^n_i \right ) + g \left ( E,\nabla E \right ) \right )

Now we have 1 equation with 1 variable that we can solve implicitly using hypre, and then we can use \(E^{{n+1}\) and \(E}n\) to construct \(E^{*\) which we can plug into the equation for \(e}{n+1}\)

Expanding f about e₀

Expanding

(T(e)) = B \left ( T_0+dT(e) \right ) = B \left ( T_0 \right ) + \left . \frac{\partial B}{\partial T} \right | _{T_0} \frac{\partial T}{\partial e} de = \frac{c}{4 \pi} a_R \left ( T_0^4 + 4T_0^3\Gamma de \right ) = B_0 \left ( 1 + 4\Gamma \frac{e-e_0}{T_0} \right )

where

Gamma = \frac{\partial T}{\partial e} = \frac{(\gamma-1)}{n k_B

and we can identify phi = -\Delta t \frac{\partial f}{\partial e} = 4 \pi \kappa_{0P} \Delta t \frac{\partial B}{\partial e} = 16 \pi \kappa_{0P} \Delta t B_0 \frac{\Gamma}{T_0

Then the equation for e becomes

frac{\partial e}{\partial t} = - \kappa_{0P} \left [ 4 \pi B_0 \left ( 1 + 4\Gamma \frac{e-e_0}{T_0} \right )-cE \right ] + \lambda \left ( 2 \frac{\kappa_{0P}}{\kappa_{0R}}-1 \right ) \mathbf{v} \cdot \nabla E -\frac{3-R_2}{2}\kappa_{0P}\frac{v^2}{c}

which will be accurate as long as \(4\Gamma \frac{\left | e-e_0 \right | }{T_0} < \xi << 1\) or \(\left | \Delta e \right | = \left | e-e_0 \right | < \xi \frac{T_0}{4 \Gamma}\)

We can calculate the time scale for this to be true using the evolution equation for the energy density

left | e-e_0 \right | = \left | \Delta e \right | = \Delta t \left | \frac{\partial e}{\partial t} \right | < \xi \frac{T_0}{4 \Gamma

which gives Delta t < \xi \frac{T_0}{4 \Gamma \left | \frac{\partial e}{\partial t} \right | } = \xi \frac{T_0}{4 \Gamma \left | - \kappa_{0P} \left [ 4 \pi B_0 \left ( 1 + 4\Gamma \frac{e-e_0}{T_0} \right )-cE \right ] + \lambda \left ( 2 \frac{\kappa_{0P}}{\kappa_{0R}}-1 \right ) \mathbf{v} \cdot \nabla E -\frac{3-R_2}{2}\kappa_{0P}\frac{v^2}{c}E \right | }

or in discretized form

elta t < \xi \frac{T^n_i \Delta t}{4 \Gamma \left | \left ( -\theta_i + \epsilon_i E^*_i + \omega_i v^n_{x,i} \left ( E^*_{i+1}-E^*_{i-1} \right ) - \xi_i E^* \right ) \right |}

Implicit Discretization 2

Now we can discretize (E*, \nabla E*) = \kappa_{0P}cE^* + \lambda \left ( 2 \frac{\kappa_{0P}}{\kappa_{0R}}-1 \right ) \mathbf{v} \cdot \nabla E^* -\frac{3-R_2}{2}\kappa_{0P}\frac{v^2}{c}E^

as

= \epsilon \left ( \psi E^{n+1}_i + \bar{\psi} E^n_i \right ) + \omega v^n_x \left ( \psi E^{n+1}_{i+1} - \psi E^{n+1}_{i-1} + \bar{\psi} E^n_{i+1}- \bar{\psi} E^n_{i-1} \right ) - \xi \left ( \psi E^{n+1}_i + \bar{\psi} E^{n}_i \right )

which along with the other terms gives

begin{eqnarray} E^{n+1}_i-E^{n}_i & = & \left [ \alpha^n_{i+1/2} \left ( \psi E^{n+1}_{i+1} + \bar{\psi} E^{n}_{i+1}- \psi E^{n+1}_{i} - \bar{\psi} E^n_{i} \right ) - \alpha^n_{i-1/2} \left ( \psi E^{n+1}_{i} + \bar{\psi} E^{n}_i - \psi E^{n+1}_{i-1} - \bar{\psi}E^{n}_{i-1} \right ) \right ] \\ & - & \left [ \zeta^n_{i+1/2} v^n_{x,i+1/2} \left ( \psi E^{n+1}_{i+1} + \bar{\psi} E^{n}_{i+1} + \psi E^{n+1}_{i} + \bar{\psi} E^n_{i} \right ) - \zeta^n_{i-1/2} v^n_{x,i-1/2}\left ( \psi E^{n+1}_{i} + \bar{\psi} E^n_i + \psi E^{n+1}_{i-1} + \bar{\psi}E^{n}_{i-1} \right ) \right ] \\ & - & \frac{1}{1+\psi \phi} \left [ - \theta + \epsilon \left ( \psi E^{n+1}_i + \bar{\psi} E^n_i \right ) + \omega v^n_x \left ( \psi E^{n+1}_{i+1} + \bar{\psi} E^n_{i+1} - \psi E^{n+1}_{i-1} - \bar{\psi} E^n_{i-1} \right ) - \xi \left ( \psi E^{n+1}_i + \bar{\psi} E^{n}_i \right ) \right ] \\ \end{eqnarray}

where the diffusion coefficient is given by

alpha_{i+1/2}=\frac{\Delta t}{\Delta x^2} \frac{c \lambda_{i+1/2}}{\kappa_{0R,i+1/2}}

and where

zeta_{i+1/2}= \frac{\Delta t}{\Delta x}\frac{3-R_{2,i+1/2}}{4}

and

epsilon_i=c\Delta t \kappa_{0P,i

and

phi_i- = \Delta t 4 \pi \kappa_{0P,i} B \left ( T^n_i \right ) \left ( \frac{4\Gamma}{T^n_i} \right )

and

theta_i = - \Delta t f(e^n_i) = \Delta t 4 \pi \kappa_{0P,i} B \left ( T^n_i \right )

and

omega_{i} = \frac{\lambda_i \Delta t}{\Delta x} \left ( \frac{\kappa_{0P,i}}{\kappa_{0R,i}}-\frac{1}{2} \right

and

xi_i = \Delta t \frac{3-R_{2,i}}{2}\kappa_{0P,i}\frac{v_i^2}{c} and

where \kappa_{0R,i+1/2} = \frac{\kappa_{0R,i}+\kappa_{0R,i+1}}{2}

and

_{i+1/2} = \frac{\left | E^n_{i+1}-E^n_{i} \right | }{2 \kappa_{0R,i+1/2} \left ( E^n_i+E^n_{i+1} \right )

and

lambda_{i+1/2} = \frac{1}{R_{i+1/2}} \left ( \coth R_{i+1/2} - \frac{1}{R_{i+1/2}} \right )

and

_{2,i+1/2} = \lambda_{i+1/2}+\lambda_{i+1/2}^2 R_{i+1/2}^

and

_{i} = \frac{\left | E^n_{i+1}-E^n_{i-1} \right | }{2 \kappa_{0R,i} E^n_{i}

and

lambda_{i} = \frac{1}{R_{i}} \left ( \coth R_{i} - \frac{1}{R_{i}} \right )

and

_{2,i} = \lambda_{i}+\lambda_{i}^2 R_{i}^

Which we can arrange into the following form

begin{eqnarray} & \left ( 1 + \psi \left ( \alpha^n_{i+1/2} + \alpha^n_{i-1/2} + \zeta^n_{i+1/2} v^n_{x,i+1/2} - \zeta^n_{i-1/2} v^n_{x,i-1/2} + \frac{\epsilon - \xi}{1+\psi \phi} \right ) \right ) E^{n+1}_i \\ - & \left ( \psi \left ( \alpha^n_{i+1/2} - \zeta^n_{i+1/2} v^n_{x,i+1/2} - \frac{\omega v^n_x}{1+\psi \phi} \right ) \right ) E^{n+1}_{i+1} \\ - & \left ( \psi \left ( \alpha^n_{i-1/2} + \zeta^n_{i-1/2} v^n_{x,i-1/2} + \frac{\omega v^n_x}{1+\psi \phi} \right ) \right ) E^{n+1}_{i-1} \\ = & \left ( 1 - \bar{\psi} \left ( \alpha^n_{i+1/2} + \alpha^n_{i-1/2} + \zeta^n_{i+1/2} v^n_{x,i+1/2} - \zeta^n_{i-1/2} v^n_{x,i-1/2} + \frac{\epsilon - \xi}{1+\psi \phi} \right ) \right ) E^{n}_i \\ + & \left ( \bar{\psi} \left ( \alpha^n_{i+1/2} - \zeta^n_{i+1/2} v^n_{x,i+1/2} - \frac{\omega v^n_x}{1+\psi \phi} \right ) \right ) E^{n}_{i+1} \\ + & \left ( \bar{\psi} \left ( \alpha^n_{i-1/2} + \zeta^n_{i-1/2} v^n_{x,i-1/2} + \frac{\omega v^n_x}{1+\psi \phi} \right ) \right ) E^{n}_{i-1} \\ + & \frac{\theta}{1+\psi \phi} \\ \end{eqnarray}

2D etc…

For 2D or 3D we have more connections to add to the matrix elements but it is very straight forward… There will be additional terms for each dimension - and the velocity components (v_x) will change, but everything else stays the same.

Initial solution vector

For the initial solution vector, we can just use Edot from the parent update (or last time step if we are on the coarse grid) to guess E ^{n+1}_i = E^{n}_i+\dot{E}^{n}_i \Delta

Coarse-Fine Boundaries

Since we are doing our implicit solves first, we can use time interpolated solutions for the implicit solve for non-refined ghost zones. To do this we just need Edot. The opacities etc… in the ghost zones can be obtained from the hydro terms.

And the radiative implicit heating in coarse ghost cells can be done along with the initial solution vector so they are available for the hydro update.

Physical Boundary Conditions

For each boundary type we need to specify \(E^{n_g, E}{n+1}_g, \kappa_g, \mbox { and } v_g\) in terms of other quantities (ie including Eⁿ⁺¹_i). The v_g and kappa_g terms should come from the hydro boundary conditions so we just need an equation for Eⁿ_g

Open (Free streaming) boundaries

We would like the radiation to leave at the free streaming limit. But we could have to modify the opacity in the ghost zone to keep the radiation energy positive

frac{c \lambda}{\kappa_{0R}} \nabla E = \mathbf{F} = cE\mathbf{n} = \frac{c \lambda_g}{\kappa_g}\frac{ \left ( E-E_g \right ) }{\Delta x

Clearly if we set \(E_g = 0\) and \(\kappa_g << \frac{1}{\Delta x} \) we should get \(\lambda = \kappa_g \Delta x\) and a free streaming flux of \(c E \mathbf{n}\)

which also corresponds to an \(\alpha = c \frac{\Delta t}{\Delta x}\)

Constant Slope Boundary

_g = 2 E_i - E_{i+1

Periodic Boundary

This is the same as internal zones - it just maps the neighbor cell to be across the domain. Hypre has built in functionality for this under for the Struct Interface

User defined radiation field/Coarse Fine boundary

This will be the boundary at internal coarse-fine boundaries, but could also be used at the physical boundary if the radiation energy were specified.

Reflecting/ZeroSlope Boundary

_g = E_

Constant radiative flux

_g = E_i - F_0 \frac{\kappa_g \Delta x}{c \lambda_g