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