M1Closure.hpp
Go to the documentation of this file.
1 // Distributed under the MIT License.
2 // See LICENSE.txt for details.
3 
4 ///\file
5 /// Defines functions to compute the M1 closure
6 
7 #pragma once
8 
9 #include "DataStructures/Tensor/TypeAliases.hpp" // IWYU pragma: keep
10 
11 // IWYU pragma: no_forward_declare Tensor
12 
13 /// \cond
14 class DataVector;
15 namespace gsl {
16 template <class>
17 class not_null;
18 } // namespace gsl
19 /// \endcond
20 
21 namespace RadiationTransport {
22 namespace M1Grey {
23 
24 /*!
25  * Compute the 2nd moment of the neutrino distribution function
26  * in the inertial frame (pressure tensor) from the 0th (energy density)
27  * and 1st (momentum density) moments. The M1 closure sets
28  * \f{align}{
29  * P_{ij} = d_{\rm thin} P_{ij,\rm thin} + d_{\rm thick} P_{ij,\rm thick}
30  * \f}
31  * with \f$P_{ij}\f$ the pressure tensor in the inertial frame, and
32  * \f$P_{ij,thin/thick}\f$ its value assuming an optically thick / optically
33  * thin medium.
34  *
35  * Following the algorithm described in Appendix A of \cite Foucart2015vpa,
36  * we choose the Minerbo closure (\cite Minerbo1978) which sets
37  * \f{align}{
38  * d_{\rm thin} = & 1.5\chi-0.5\\
39  * d_{\rm thick} = & 1.5-1.5\chi \\
40  * \chi = & \frac{1}{3} + \xi^2 \frac{6-2\xi+6\xi^2}{15} \\
41  * \xi = & H^aH_a/J^2
42  * \f}
43  * with \f$J\f$ the energy density in the
44  * frame comoving with the fluid, and \f$H^a\f$ the momentum density
45  * in that frame.
46  * The optically thick closure is
47  * \f{align}{
48  * K^{ab,\rm thick} = \frac{J_{\rm thick}}{3} (g^{ab}+u^a u^b)
49  * \f}
50  * with \f$ K^{ab,\rm thick} \f$ the pressure tensor in the frame comoving with
51  * the fluid, and \f$ u^a \f$ the fluid 4-velocity, while the optically thin
52  * closure is
53  * \f{align}{
54  * P^{ij,\rm thin} = \frac{S^i S^j}{S^k S_k} E
55  * \f}
56  * with \f$ E \f$, \f$ S^{i}
57  * \f$ the 0th and 1st moments in the inertial frame (see paper for the
58  * computation of \f$P_{ij,thick}\f$ from the other moments).
59  *
60  * The main step of this function is to solve the equation
61  * \f{align}{
62  * \frac{\xi^2 J^2 - H^aH_a}{E^2} = 0
63  * \f}
64  * for \f$\xi\f$, with \f$J\f$ and \f$H^a\f$ being functions
65  * of \f$\xi\f$ and of the input variables (the 0th and 1st moments
66  * in the inertial frame). This is done by separating \f$H^a\f$ as
67  * \f{align}{
68  * H_a = - H_t t_a - H_v v_a - H_f F_a
69  * \f}
70  * (with \f$v^a\f$ the fluid 3-velocity and \f$t^a\f$ the unit normal
71  * to the slice) and then writing each of the
72  * variables \f$J\f$, \f$H_n\f$, \f$H_v\f$, \f$H_f\f$ as
73  * \f{align}{
74  * X = X_0 + d_{\rm thin} X_{\rm thin} + d_{\rm thick} X_{\rm thick}
75  * \f}
76  * so that evaluating
77  * \f{align}{
78  * \frac{\xi^2 J^2 - H^aH_a}{E^2} = 0
79  * \f}
80  * for a given \f$\xi\f$ only requires recomputing \f$d_{\rm thin,thick}\f$
81  * and their derivatives with respect to \f$\xi\f$.
82  * We perform the root-finding using a Newton-Raphson algorithm, with the
83  * accuracy set by the variable root_find_number_of_digits (6 significant digits
84  * at the moment).
85  *
86  * The function returns the closure factors \f$\xi\f$ (to be used as initial
87  * guess for this function at the next step), the pressure tensor \f$P_{ij}\f$,
88  * and the neutrino moments in the frame comoving with the fluid.
89  * The momentum density in the frame comoving with the fluid
90  * is decomposed into its normal component \f$ H^a t_a\f$, and its spatial
91  * components \f$ \gamma_{ia} H^a\f$.
92  */
93 void M1Closure(
94  gsl::not_null<Scalar<DataVector>*> closure_factor,
95  gsl::not_null<tnsr::II<DataVector, 3, Frame::Inertial>*> pressure_tensor,
96  gsl::not_null<Scalar<DataVector>*> comoving_energy_density,
97  gsl::not_null<Scalar<DataVector>*> comoving_momentum_density_normal,
98  gsl::not_null<tnsr::i<DataVector, 3, Frame::Inertial>*>
99  comoving_momentum_density_spatial,
100  const Scalar<DataVector>& energy_density,
101  const tnsr::i<DataVector, 3, Frame::Inertial>& momentum_density,
102  const tnsr::I<DataVector, 3, Frame::Inertial>& fluid_velocity,
103  const Scalar<DataVector>& fluid_lorentz_factor,
104  const tnsr::ii<DataVector, 3, Frame::Inertial>& spatial_metric,
105  const tnsr::II<DataVector, 3, Frame::Inertial>&
106  inv_spatial_metric) noexcept;
107 
108 } // namespace M1Grey
109 } // namespace RadiationTransport
Implementations from the Guideline Support Library.
Definition: ConservativeFromPrimitive.hpp:10
Namespace for all radiation transport algorithms.
Definition: M1Closure.cpp:45
tnsr::ii< DataType, SpatialDim, Frame > spatial_metric(const tnsr::aa< DataType, SpatialDim, Frame > &spacetime_metric) noexcept
Compute spatial metric from spacetime metric.
Definition: ComputeSpacetimeQuantities.cpp:43
Defines a list of useful type aliases for tensors.
Stores a collection of function values.
Definition: DataVector.hpp:42
Tensor< T, Symmetry<>, index_list<> > Scalar
Scalar type.
Definition: TypeAliases.hpp:21
Require a pointer to not be a nullptr
Definition: ConservativeFromPrimitive.hpp:12
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