Namespaces | Functions
RadiationTransport::M1Grey Namespace Reference

Namespace for the grey-M1 radiation transport scheme. More...

Namespaces

 Tags
 Tags for the evolution of neutrinos using a grey M1 scheme.
 

Functions

void M1Closure (const gsl::not_null< Scalar< DataVector > *> closure_factor, const gsl::not_null< tnsr::II< DataVector, 3, Frame::Inertial > *> pressure_tensor, const gsl::not_null< Scalar< DataVector > *> comoving_energy_density, const gsl::not_null< Scalar< DataVector > *> comoving_momentum_density_normal, const gsl::not_null< tnsr::i< DataVector, 3, Frame::Inertial > *> comoving_momentum_density_spatial, const Scalar< DataVector > &energy_density, const tnsr::i< DataVector, 3, Frame::Inertial > &momentum_density, const tnsr::I< DataVector, 3, Frame::Inertial > &fluid_velocity, const Scalar< DataVector > &fluid_lorentz_factor, const tnsr::ii< DataVector, 3, Frame::Inertial > &spatial_metric, const tnsr::II< DataVector, 3, Frame::Inertial > &inv_spatial_metric) noexcept
 
void M1Closure (gsl::not_null< Scalar< DataVector > *> closure_factor, gsl::not_null< tnsr::II< DataVector, 3, Frame::Inertial > *> pressure_tensor, gsl::not_null< Scalar< DataVector > *> comoving_energy_density, gsl::not_null< Scalar< DataVector > *> comoving_momentum_density_normal, gsl::not_null< tnsr::i< DataVector, 3, Frame::Inertial > *> comoving_momentum_density_spatial, const Scalar< DataVector > &energy_density, const tnsr::i< DataVector, 3, Frame::Inertial > &momentum_density, const tnsr::I< DataVector, 3, Frame::Inertial > &fluid_velocity, const Scalar< DataVector > &fluid_lorentz_factor, const tnsr::ii< DataVector, 3, Frame::Inertial > &spatial_metric, const tnsr::II< DataVector, 3, Frame::Inertial > &inv_spatial_metric) noexcept
 

Detailed Description

Namespace for the grey-M1 radiation transport scheme.

Function Documentation

◆ M1Closure()

void RadiationTransport::M1Grey::M1Closure ( gsl::not_null< Scalar< DataVector > *>  closure_factor,
gsl::not_null< tnsr::II< DataVector, 3, Frame::Inertial > *>  pressure_tensor,
gsl::not_null< Scalar< DataVector > *>  comoving_energy_density,
gsl::not_null< Scalar< DataVector > *>  comoving_momentum_density_normal,
gsl::not_null< tnsr::i< DataVector, 3, Frame::Inertial > *>  comoving_momentum_density_spatial,
const Scalar< DataVector > &  energy_density,
const tnsr::i< DataVector, 3, Frame::Inertial > &  momentum_density,
const tnsr::I< DataVector, 3, Frame::Inertial > &  fluid_velocity,
const Scalar< DataVector > &  fluid_lorentz_factor,
const tnsr::ii< DataVector, 3, Frame::Inertial > &  spatial_metric,
const tnsr::II< DataVector, 3, Frame::Inertial > &  inv_spatial_metric 
)
noexcept

Compute the 2nd moment of the neutrino distribution function in the inertial frame (pressure tensor) from the 0th (energy density) and 1st (momentum density) moments. The M1 closure sets

\begin{align} P_{ij} = d_{\rm thin} P_{ij,\rm thin} + d_{\rm thick} P_{ij,\rm thick} \end{align}

with \(P_{ij}\) the pressure tensor in the inertial frame, and \(P_{ij,thin/thick}\) its value assuming an optically thick / optically thin medium.

Following the algorithm described in Appendix A of [7], we choose the Minerbo closure ([21]) which sets

\begin{align} d_{\rm thin} = & 1.5\chi-0.5\\ d_{\rm thick} = & 1.5-1.5\chi \\ \chi = & \frac{1}{3} + \xi^2 \frac{6-2\xi+6\xi^2}{15} \\ \xi = & H^aH_a/J^2 \end{align}

with \(J\) the energy density in the frame comoving with the fluid, and \(H^a\) the momentum density in that frame. The optically thick closure is

\begin{align} K^{ab,\rm thick} = \frac{J_{\rm thick}}{3} (g^{ab}+u^a u^b) \end{align}

with \( K^{ab,\rm thick} \) the pressure tensor in the frame comoving with the fluid, and \( u^a \) the fluid 4-velocity, while the optically thin closure is

\begin{align} P^{ij,\rm thin} = \frac{S^i S^j}{S^k S_k} E \end{align}

with \( E \), \( S^{i} \) the 0th and 1st moments in the inertial frame (see paper for the computation of \(P_{ij,thick}\) from the other moments).

The main step of this function is to solve the equation

\begin{align} \frac{\xi^2 J^2 - H^aH_a}{E^2} = 0 \end{align}

for \(\xi\), with \(J\) and \(H^a\) being functions of \(\xi\) and of the input variables (the 0th and 1st moments in the inertial frame). This is done by separating \(H^a\) as

\begin{align} H_a = - H_t t_a - H_v v_a - H_f F_a \end{align}

(with \(v^a\) the fluid 3-velocity and \(t^a\) the unit normal to the slice) and then writing each of the variables \(J\), \(H_n\), \(H_v\), \(H_f\) as

\begin{align} X = X_0 + d_{\rm thin} X_{\rm thin} + d_{\rm thick} X_{\rm thick} \end{align}

so that evaluating

\begin{align} \frac{\xi^2 J^2 - H^aH_a}{E^2} = 0 \end{align}

for a given \(\xi\) only requires recomputing \(d_{\rm thin,thick}\) and their derivatives with respect to \(\xi\). We perform the root-finding using a Newton-Raphson algorithm, with the accuracy set by the variable root_find_number_of_digits (6 significant digits at the moment).

The function returns the closure factors \(\xi\) (to be used as initial guess for this function at the next step), the pressure tensor \(P_{ij}\), and the neutrino moments in the frame comoving with the fluid. The momentum density in the frame comoving with the fluid is decomposed into its normal component \( H^a t_a\), and its spatial components \( \gamma_{ia} H^a\).