Line data Source code
1 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 : #include "Evolution/Systems/RadiationTransport/M1Grey/Tags.hpp" // IWYU pragma: keep
11 : #include "PointwiseFunctions/GeneralRelativity/Tags.hpp"
12 : #include "PointwiseFunctions/Hydro/Tags.hpp"
13 : #include "Utilities/Gsl.hpp" // for not_...
14 : #include "Utilities/TMPL.hpp" // for EXPAND_PACK_LEFT_TO...
15 :
16 : // IWYU pragma: no_forward_declare Tensor
17 :
18 : /// \cond
19 : class DataVector;
20 : /// \endcond
21 :
22 : namespace RadiationTransport {
23 : namespace M1Grey {
24 :
25 : // Implementation of the M1 closure for an
26 : // individual species
27 : namespace detail {
28 : void compute_closure_impl(
29 : gsl::not_null<Scalar<DataVector>*> closure_factor,
30 : gsl::not_null<tnsr::II<DataVector, 3, Frame::Inertial>*> pressure_tensor,
31 : gsl::not_null<Scalar<DataVector>*> comoving_energy_density,
32 : gsl::not_null<Scalar<DataVector>*> comoving_momentum_density_normal,
33 : gsl::not_null<tnsr::i<DataVector, 3, Frame::Inertial>*>
34 : comoving_momentum_density_spatial,
35 : const Scalar<DataVector>& energy_density,
36 : const tnsr::i<DataVector, 3, Frame::Inertial>& momentum_density,
37 : const tnsr::I<DataVector, 3, Frame::Inertial>& fluid_velocity,
38 : const Scalar<DataVector>& fluid_lorentz_factor,
39 : const tnsr::ii<DataVector, 3, Frame::Inertial>& spatial_metric,
40 : const tnsr::II<DataVector, 3, Frame::Inertial>& inv_spatial_metric);
41 : } // namespace detail
42 :
43 : template <typename NeutrinoSpeciesList>
44 0 : struct ComputeM1Closure;
45 : /*!
46 : * Compute the 2nd moment of the neutrino distribution function
47 : * in the inertial frame (pressure tensor) from the 0th (energy density)
48 : * and 1st (momentum density) moments. The M1 closure sets
49 : * \f{align}{
50 : * P_{ij} = d_{\rm thin} P_{ij,\rm thin} + d_{\rm thick} P_{ij,\rm thick}
51 : * \f}
52 : * with \f$P_{ij}\f$ the pressure tensor in the inertial frame, and
53 : * \f$P_{ij,thin/thick}\f$ its value assuming an optically thick / optically
54 : * thin medium.
55 : *
56 : * Following the algorithm described in Appendix A of \cite Foucart2015vpa,
57 : * we choose the Minerbo closure (\cite Minerbo1978) which sets
58 : * \f{align}{
59 : * d_{\rm thin} = & 1.5\chi-0.5\\
60 : * d_{\rm thick} = & 1.5-1.5\chi \\
61 : * \chi = & \frac{1}{3} + \xi^2 \frac{6-2\xi+6\xi^2}{15} \\
62 : * \xi = & H^aH_a/J^2
63 : * \f}
64 : * with \f$J\f$ the energy density in the
65 : * frame comoving with the fluid, and \f$H^a\f$ the momentum density
66 : * in that frame.
67 : * The optically thick closure is
68 : * \f{align}{
69 : * K^{ab,\rm thick} = \frac{J_{\rm thick}}{3} (g^{ab}+u^a u^b)
70 : * \f}
71 : * with \f$ K^{ab,\rm thick} \f$ the pressure tensor in the frame comoving with
72 : * the fluid, and \f$ u^a \f$ the fluid 4-velocity, while the optically thin
73 : * closure is
74 : * \f{align}{
75 : * P^{ij,\rm thin} = \frac{S^i S^j}{S^k S_k} E
76 : * \f}
77 : * with \f$ E \f$, \f$ S^{i}
78 : * \f$ the 0th and 1st moments in the inertial frame (see paper for the
79 : * computation of \f$P_{ij,thick}\f$ from the other moments).
80 : *
81 : * The main step of this function is to solve the equation
82 : * \f{align}{
83 : * \frac{\xi^2 J^2 - H^aH_a}{E^2} = 0
84 : * \f}
85 : * for \f$\xi\f$, with \f$J\f$ and \f$H^a\f$ being functions
86 : * of \f$\xi\f$ and of the input variables (the 0th and 1st moments
87 : * in the inertial frame). This is done by separating \f$H^a\f$ as
88 : * \f{align}{
89 : * H_a = - H_t t_a - H_v v_a - H_f F_a
90 : * \f}
91 : * (with \f$v^a\f$ the fluid 3-velocity and \f$t^a\f$ the unit normal
92 : * to the slice) and then writing each of the
93 : * variables \f$J\f$, \f$H_n\f$, \f$H_v\f$, \f$H_f\f$ as
94 : * \f{align}{
95 : * X = X_0 + d_{\rm thin} X_{\rm thin} + d_{\rm thick} X_{\rm thick}
96 : * \f}
97 : * so that evaluating
98 : * \f{align}{
99 : * \frac{\xi^2 J^2 - H^aH_a}{E^2} = 0
100 : * \f}
101 : * for a given \f$\xi\f$ only requires recomputing \f$d_{\rm thin,thick}\f$
102 : * and their derivatives with respect to \f$\xi\f$.
103 : * We perform the root-finding using a Newton-Raphson algorithm, with the
104 : * accuracy set by the variable root_find_number_of_digits (6 significant digits
105 : * at the moment).
106 : *
107 : * The function returns the closure factors \f$\xi\f$ (to be used as initial
108 : * guess for this function at the next step), the pressure tensor \f$P_{ij}\f$,
109 : * and the neutrino moments in the frame comoving with the fluid.
110 : * The momentum density in the frame comoving with the fluid
111 : * is decomposed into its normal component \f$ H^a t_a\f$, and its spatial
112 : * components \f$ \gamma_{ia} H^a\f$.
113 : */
114 : template <typename... NeutrinoSpecies>
115 1 : struct ComputeM1Closure<tmpl::list<NeutrinoSpecies...>> {
116 0 : using return_tags =
117 : tmpl::list<Tags::ClosureFactor<NeutrinoSpecies>...,
118 : Tags::TildeP<Frame::Inertial, NeutrinoSpecies>...,
119 : Tags::TildeJ<NeutrinoSpecies>...,
120 : Tags::TildeHNormal<NeutrinoSpecies>...,
121 : Tags::TildeHSpatial<Frame::Inertial, NeutrinoSpecies>...>;
122 :
123 0 : using argument_tags =
124 : tmpl::list<Tags::TildeE<Frame::Inertial, NeutrinoSpecies>...,
125 : Tags::TildeS<Frame::Inertial, NeutrinoSpecies>...,
126 : hydro::Tags::SpatialVelocity<DataVector, 3>,
127 : hydro::Tags::LorentzFactor<DataVector>,
128 : gr::Tags::SpatialMetric<DataVector, 3>,
129 : gr::Tags::InverseSpatialMetric<DataVector, 3>>;
130 :
131 0 : static void apply(
132 : const gsl::not_null<typename Tags::ClosureFactor<
133 : NeutrinoSpecies>::type*>... closure_factor,
134 : const gsl::not_null<typename Tags::TildeP<
135 : Frame::Inertial, NeutrinoSpecies>::type*>... tilde_p,
136 : const gsl::not_null<
137 : typename Tags::TildeJ<NeutrinoSpecies>::type*>... tilde_j,
138 : const gsl::not_null<
139 : typename Tags::TildeHNormal<NeutrinoSpecies>::type*>... tilde_hn,
140 : const gsl::not_null<typename Tags::TildeHSpatial<
141 : Frame::Inertial, NeutrinoSpecies>::type*>... tilde_hi,
142 : const typename Tags::TildeE<Frame::Inertial,
143 : NeutrinoSpecies>::type&... tilde_e,
144 : const typename Tags::TildeS<Frame::Inertial,
145 : NeutrinoSpecies>::type&... tilde_s,
146 : const tnsr::I<DataVector, 3>& spatial_velocity,
147 : const Scalar<DataVector>& lorentz_factor,
148 : const tnsr::ii<DataVector, 3>& spatial_metric,
149 : const tnsr::II<DataVector, 3>& inv_spatial_metric) {
150 : EXPAND_PACK_LEFT_TO_RIGHT(detail::compute_closure_impl(
151 : closure_factor, tilde_p, tilde_j, tilde_hn, tilde_hi, tilde_e, tilde_s,
152 : spatial_velocity, lorentz_factor, spatial_metric, inv_spatial_metric));
153 : }
154 : };
155 :
156 : } // namespace M1Grey
157 : } // namespace RadiationTransport
|