Line data Source code
1 0 : // Distributed under the MIT License.
2 : // See LICENSE.txt for details.
3 :
4 : #pragma once
5 :
6 : #include <cstdint>
7 :
8 : #include "DataStructures/Tensor/TypeAliases.hpp"
9 : #include "Domain/Tags.hpp"
10 : #include "Evolution/DgSubcell/Tags/Coordinates.hpp"
11 : #include "Evolution/Systems/GrMhd/ValenciaDivClean/TagsDeclarations.hpp"
12 : #include "NumericalAlgorithms/LinearOperators/PartialDerivatives.hpp"
13 : #include "PointwiseFunctions/GeneralRelativity/TagsDeclarations.hpp"
14 : #include "PointwiseFunctions/Hydro/Tags.hpp"
15 : #include "PointwiseFunctions/Hydro/TagsDeclarations.hpp"
16 : #include "Utilities/TMPL.hpp"
17 :
18 : /// \cond
19 : class DataVector;
20 : template <size_t Dim>
21 : class Mesh;
22 :
23 : namespace Spectral {
24 : enum class Quadrature : uint8_t;
25 : } // namespace Spectral
26 : namespace Tags {
27 : template <typename>
28 : struct Source;
29 : } // namespace Tags
30 : /// \endcond
31 :
32 : namespace grmhd {
33 : namespace ValenciaDivClean {
34 : namespace detail {
35 : void cartoon_sources_impl(
36 : gsl::not_null<tnsr::i<DataVector, 3, Frame::Inertial>*> source_tilde_s,
37 : gsl::not_null<tnsr::I<DataVector, 3, Frame::Inertial>*> source_tilde_b,
38 :
39 : const Scalar<DataVector>& pressure_star,
40 : const tnsr::i<DataVector, 3, Frame::Inertial>& magnetic_field_one_form,
41 : const Scalar<DataVector>& magnetic_field_dot_spatial_velocity,
42 :
43 : const tnsr::i<DataVector, 3, Frame::Inertial>& tilde_s,
44 : const tnsr::I<DataVector, 3, Frame::Inertial>& tilde_b,
45 : const Scalar<DataVector>& tilde_phi,
46 : const tnsr::I<DataVector, 3, Frame::Inertial>& spatial_velocity,
47 : const Scalar<DataVector>& lorentz_factor, const Scalar<DataVector>& lapse,
48 : const tnsr::I<DataVector, 3, Frame::Inertial>& shift,
49 : const tnsr::ii<DataVector, 3, Frame::Inertial>& spatial_metric,
50 : const tnsr::II<DataVector, 3, Frame::Inertial>& inv_spatial_metric,
51 : const Scalar<DataVector>& sqrt_det_spatial_metric,
52 : const tnsr::I<DataVector, 3, Frame::Inertial>& inertial_coords,
53 : Spectral::Quadrature cartoon_quadrature);
54 :
55 : void sources_impl(
56 : gsl::not_null<Scalar<DataVector>*> source_tilde_tau,
57 : gsl::not_null<tnsr::i<DataVector, 3, Frame::Inertial>*> source_tilde_s,
58 : gsl::not_null<tnsr::I<DataVector, 3, Frame::Inertial>*> source_tilde_b,
59 : gsl::not_null<Scalar<DataVector>*> source_tilde_phi,
60 :
61 : gsl::not_null<tnsr::I<DataVector, 3, Frame::Inertial>*> tilde_s_up,
62 : gsl::not_null<tnsr::II<DataVector, 3, Frame::Inertial>*> densitized_stress,
63 : gsl::not_null<Scalar<DataVector>*> h_rho_w_squared_plus_b_squared,
64 :
65 : const Scalar<DataVector>& magnetic_field_dot_spatial_velocity,
66 : const Scalar<DataVector>& magnetic_field_squared,
67 : const Scalar<DataVector>& one_over_w_squared,
68 : const Scalar<DataVector>& pressure_star,
69 : const tnsr::I<DataVector, 3, Frame::Inertial>&
70 : trace_spatial_christoffel_second,
71 :
72 : const Scalar<DataVector>& tilde_d, const Scalar<DataVector>& tilde_ye,
73 : const Scalar<DataVector>& tilde_tau,
74 : const tnsr::i<DataVector, 3, Frame::Inertial>& tilde_s,
75 : const tnsr::I<DataVector, 3, Frame::Inertial>& tilde_b,
76 : const Scalar<DataVector>& tilde_phi, const Scalar<DataVector>& lapse,
77 : const Scalar<DataVector>& sqrt_det_spatial_metric,
78 : const tnsr::II<DataVector, 3, Frame::Inertial>& inv_spatial_metric,
79 : const tnsr::i<DataVector, 3, Frame::Inertial>& d_lapse,
80 : const tnsr::iJ<DataVector, 3, Frame::Inertial>& d_shift,
81 : const tnsr::ijj<DataVector, 3, Frame::Inertial>& d_spatial_metric,
82 : const tnsr::I<DataVector, 3, Frame::Inertial>& spatial_velocity,
83 : const Scalar<DataVector>& lorentz_factor,
84 : const tnsr::I<DataVector, 3, Frame::Inertial>& magnetic_field,
85 :
86 : const Scalar<DataVector>& rest_mass_density,
87 : const Scalar<DataVector>& electron_fraction,
88 : const Scalar<DataVector>& pressure,
89 : const Scalar<DataVector>& specific_internal_energy,
90 : const tnsr::ii<DataVector, 3, Frame::Inertial>& extrinsic_curvature,
91 : double constraint_damping_parameter);
92 : } // namespace detail
93 :
94 : /*!
95 : * \brief Compute the source terms for the flux-conservative Valencia
96 : * formulation of GRMHD with divergence cleaning, coupled with electron
97 : * fraction.
98 : *
99 : * A flux-conservative system has the generic form:
100 : * \f[
101 : * \partial_t U_i + \partial_m F^m(U_i) = S(U_i)
102 : * \f]
103 : *
104 : * where \f$F^a()\f$ denotes the flux of a conserved variable \f$U_i\f$ and
105 : * \f$S()\f$ denotes the source term for the conserved variable.
106 : *
107 : * For the Valencia formulation:
108 : * \f{align*}
109 : * S({\tilde D}) = & 0\\
110 : * S({\tilde S}_i) = & \frac{1}{2} \alpha {\tilde S}^{mn} \partial_i \gamma_{mn}
111 : * + {\tilde S}_m \partial_i \beta^m - ({\tilde D} + {\tilde \tau}) \partial_i
112 : * \alpha \\
113 : * S({\tilde \tau}) = & \alpha {\tilde S}^{mn} K_{mn} - {\tilde S}^m \partial_m
114 : * \alpha \\
115 : * S({\tilde B}^i) = & -{\tilde B}^m \partial_m \beta^i + {\tilde \Phi}
116 : * \gamma^{im} \partial_m \alpha + \alpha {\tilde \Phi} \left( \frac{1}{2}
117 : * \gamma^{il} \gamma^{jk} - \gamma^{ij} \gamma^{lk} \right) \partial_l
118 : * \gamma_{jk} \\ S({\tilde \Phi}) = & {\tilde B}^k \partial_k \alpha - \alpha K
119 : * {\tilde \Phi}
120 : * - \alpha \kappa {\tilde \Phi}
121 : * \f}
122 : *
123 : * where
124 : * \f{align*}
125 : * {\tilde S}^i = & {\tilde S}_m \gamma^{im} \\
126 : * {\tilde S}^{ij} = & \sqrt{\gamma} \left[ \left(h \rho W^2 + B^n B_n \right)
127 : * v^i v^j + \left(p + p_m \right) \gamma^{ij} - B^n v_n \left( B^i v^j + B^j
128 : * v^i \right) - \frac{B^i B^j}{W^2} \right]
129 : * \f}
130 : *
131 : * where \f${\tilde D}\f$, \f${\tilde S}_i\f$, \f${\tilde \tau}\f$, \f${\tilde
132 : * B}^i\f$, and \f${\tilde \Phi}\f$ are a generalized mass-energy density,
133 : * momentum density, specific internal energy density, magnetic field, and
134 : * divergence cleaning field. Furthermore, \f$\gamma\f$ is the determinant of
135 : * the spatial metric \f$\gamma_{ij}\f$, \f$\rho\f$ is the rest mass density,
136 : * \f$W\f$ is the Lorentz factor, \f$h\f$ is the specific enthalpy, \f$v^i\f$ is
137 : * the spatial velocity, \f$B^k\f$ is the magnetic field, \f$p\f$ is the
138 : * pressure, \f$p_m = \frac{1}{2} \left[ \left( B^n v_n \right)^2 + B^n B_n /
139 : * W^2 \right]\f$ is the magnetic pressure, \f$\alpha\f$ is the lapse,
140 : * \f$\beta^i\f$ is the shift, \f$K\f$ is the trace of the extrinsic curvature
141 : * \f$K_{ij}\f$, and \f$\kappa\f$ is a damping parameter that damps violations
142 : * of the divergence-free (no-monopole) condition \f$\Phi = \partial_i {\tilde
143 : * B}^i = 0\f$ .
144 : *
145 : * When using the cartoon method, the flux of vector-quantity conserved
146 : * variables result in terms that are not flux-like but can be interpreted as
147 : * sources. For spherical symmetry these extra source terms are
148 : * \f{align*}
149 : * S_{\mathrm{spherical}, \tilde{S}_x} &=
150 : * \frac{2}{x} \alpha \sqrt{\gamma} (p + p_m) \\
151 : * S_{\mathrm{spherical}, \tilde{B}^x} &=
152 : * \frac{1}{x}\alpha (\gamma^{yy} + \gamma^{zz}) \tilde{\Phi}
153 : * \f}
154 : *
155 : * and for axial symmetry they are
156 : * \f{align*}
157 : * S_{\mathrm{axial}, \tilde{S}_x} &=
158 : * \frac{1}{x} \left( \alpha \sqrt{\gamma} (p + p_m)
159 : * + \tilde{S}_z v^z_\mathrm{tr}
160 : * - \frac{\alpha}{W} b_z \tilde{B}^z \right) \\
161 : * S_{\mathrm{axial}, \tilde{S}_z} &=
162 : * \frac{1}{x} \left( -\tilde{S}_x v^z_\mathrm{tr}
163 : * + \frac{\alpha}{W} b_x \tilde{B}^z \right) \\
164 : * S_{\mathrm{axial}, \tilde{B}^x} &=
165 : * \frac{1}{x} \left( \alpha \gamma^{zz}\tilde{\Phi}
166 : * - \alpha v^z \tilde{B}^z
167 : * + \tilde{B}^z v^z_\mathrm{tr} \right) \\
168 : * S_{\mathrm{axial}, \tilde{B}^z} &=
169 : * \frac{1}{x} \left( -\alpha \gamma^{zx} \tilde{\Phi}
170 : * + \alpha v^x \tilde{B}^z
171 : * - \tilde{B}^x v^z_\mathrm{tr} \right),
172 : * \f}
173 : *
174 : * where \f$v^i_\mathrm{tr} = \alpha v^i - \beta^i\f$ is the transport
175 : * velocity. These terms are associated with the $x$ component of
176 : * $\tilde{S}_j$ and $\tilde{B}^j$ because that is the divergence direction we
177 : * have cartoonified.
178 : *
179 : * \note For the electron fraction side, the source term is currently set to
180 : * \f$S(\tilde{Y}_e) = 0\f$ where the conserved variable \f$\tilde{Y}_e\f$ is a
181 : * generalized electron fraction. Implementing the source term using neutrino
182 : * scheme is in progress (Last update : Oct 2022).
183 : *
184 : */
185 1 : struct ComputeSources {
186 0 : using return_tags =
187 : tmpl::list<::Tags::Source<grmhd::ValenciaDivClean::Tags::TildeTau>,
188 : ::Tags::Source<grmhd::ValenciaDivClean::Tags::TildeS<>>,
189 : ::Tags::Source<grmhd::ValenciaDivClean::Tags::TildeB<>>,
190 : ::Tags::Source<grmhd::ValenciaDivClean::Tags::TildePhi>>;
191 :
192 0 : using argument_tags =
193 : tmpl::list<grmhd::ValenciaDivClean::Tags::TildeD,
194 : grmhd::ValenciaDivClean::Tags::TildeYe,
195 : grmhd::ValenciaDivClean::Tags::TildeTau,
196 : grmhd::ValenciaDivClean::Tags::TildeS<>,
197 : grmhd::ValenciaDivClean::Tags::TildeB<>,
198 : grmhd::ValenciaDivClean::Tags::TildePhi,
199 : hydro::Tags::SpatialVelocity<DataVector, 3>,
200 : hydro::Tags::MagneticField<DataVector, 3>,
201 : hydro::Tags::RestMassDensity<DataVector>,
202 : hydro::Tags::ElectronFraction<DataVector>,
203 : hydro::Tags::SpecificInternalEnergy<DataVector>,
204 : hydro::Tags::LorentzFactor<DataVector>,
205 : hydro::Tags::Pressure<DataVector>, gr::Tags::Lapse<DataVector>,
206 : gr::Tags::Shift<DataVector, 3>,
207 : ::Tags::deriv<gr::Tags::Lapse<DataVector>, tmpl::size_t<3>,
208 : Frame::Inertial>,
209 : ::Tags::deriv<gr::Tags::Shift<DataVector, 3>, tmpl::size_t<3>,
210 : Frame::Inertial>,
211 : gr::Tags::SpatialMetric<DataVector, 3>,
212 : ::Tags::deriv<gr::Tags::SpatialMetric<DataVector, 3>,
213 : tmpl::size_t<3>, Frame::Inertial>,
214 : gr::Tags::InverseSpatialMetric<DataVector, 3>,
215 : gr::Tags::SqrtDetSpatialMetric<DataVector>,
216 : gr::Tags::ExtrinsicCurvature<DataVector, 3>,
217 : grmhd::ValenciaDivClean::Tags::ConstraintDampingParameter,
218 : evolution::dg::subcell::Tags::Coordinates<3, Frame::Inertial>,
219 : domain::Tags::Mesh<3>>;
220 :
221 0 : static void apply(
222 : gsl::not_null<Scalar<DataVector>*> source_tilde_tau,
223 : gsl::not_null<tnsr::i<DataVector, 3, Frame::Inertial>*> source_tilde_s,
224 : gsl::not_null<tnsr::I<DataVector, 3, Frame::Inertial>*> source_tilde_b,
225 : gsl::not_null<Scalar<DataVector>*> source_tilde_phi,
226 : const Scalar<DataVector>& tilde_d, const Scalar<DataVector>& tilde_ye,
227 : const Scalar<DataVector>& tilde_tau,
228 : const tnsr::i<DataVector, 3, Frame::Inertial>& tilde_s,
229 : const tnsr::I<DataVector, 3, Frame::Inertial>& tilde_b,
230 : const Scalar<DataVector>& tilde_phi,
231 : const tnsr::I<DataVector, 3, Frame::Inertial>& spatial_velocity,
232 : const tnsr::I<DataVector, 3, Frame::Inertial>& magnetic_field,
233 : const Scalar<DataVector>& rest_mass_density,
234 : const Scalar<DataVector>& electron_fraction,
235 : const Scalar<DataVector>& specific_internal_energy,
236 : const Scalar<DataVector>& lorentz_factor,
237 : const Scalar<DataVector>& pressure, const Scalar<DataVector>& lapse,
238 : const tnsr::I<DataVector, 3, Frame::Inertial>& shift,
239 : const tnsr::i<DataVector, 3, Frame::Inertial>& d_lapse,
240 : const tnsr::iJ<DataVector, 3, Frame::Inertial>& d_shift,
241 : const tnsr::ii<DataVector, 3, Frame::Inertial>& spatial_metric,
242 : const tnsr::ijj<DataVector, 3, Frame::Inertial>& d_spatial_metric,
243 : const tnsr::II<DataVector, 3, Frame::Inertial>& inv_spatial_metric,
244 : const Scalar<DataVector>& sqrt_det_spatial_metric,
245 : const tnsr::ii<DataVector, 3, Frame::Inertial>& extrinsic_curvature,
246 : double constraint_damping_parameter,
247 : const tnsr::I<DataVector, 3, Frame::Inertial>& inertial_coords,
248 : const Mesh<3>& dg_mesh);
249 : };
250 : } // namespace ValenciaDivClean
251 : } // namespace grmhd
|
