Line data Source code
1 0 : // Distributed under the MIT License. 2 : // See LICENSE.txt for details. 3 : 4 : #pragma once 5 : 6 : #include <cstddef> 7 : 8 : #include "DataStructures/VariablesTag.hpp" 9 : #include "Evolution/Systems/GrMhd/ValenciaDivClean/BoundaryConditions/BoundaryCondition.hpp" 10 : #include "Evolution/Systems/GrMhd/ValenciaDivClean/BoundaryCorrections/BoundaryCorrection.hpp" 11 : #include "Evolution/Systems/GrMhd/ValenciaDivClean/Characteristics.hpp" 12 : #include "Evolution/Systems/GrMhd/ValenciaDivClean/ConservativeFromPrimitive.hpp" 13 : #include "Evolution/Systems/GrMhd/ValenciaDivClean/NewmanHamlin.hpp" 14 : #include "Evolution/Systems/GrMhd/ValenciaDivClean/PrimitiveFromConservative.hpp" 15 : #include "Evolution/Systems/GrMhd/ValenciaDivClean/Tags.hpp" 16 : #include "Evolution/Systems/GrMhd/ValenciaDivClean/TimeDerivativeTerms.hpp" 17 : #include "PointwiseFunctions/GeneralRelativity/Tags.hpp" 18 : #include "PointwiseFunctions/Hydro/Tags.hpp" 19 : #include "Utilities/TMPL.hpp" 20 : 21 : /// \ingroup EvolutionSystemsGroup 22 : /// \brief Items related to general relativistic magnetohydrodynamics (GRMHD) 23 : namespace grmhd { 24 : /// Items related to the Valencia formulation of ideal GRMHD with divergence 25 : /// cleaning coupled with electron fraction. 26 : /// 27 : /// References: 28 : /// - Numerical 3+1 General Relativistic Magnetohydrodynamics: A Local 29 : /// Characteristic Approach \cite Anton2006 30 : /// - GRHydro: a new open-source general-relativistic magnetohydrodynamics code 31 : /// for the Einstein toolkit \cite Moesta2014 32 : /// - Black hole-neutron star mergers with a hot nuclear equation of state: 33 : /// outflow and neutrino-cooled disk for a low-mass, high-spin case 34 : /// \cite Deaton2013 35 : /// - A Simflowny-based finite-difference code for high-performance computing 36 : /// in Relativity \cite Palenzuela2018 37 : /// 38 : namespace ValenciaDivClean { 39 : 40 : /*! 41 : * \brief Ideal general relativistic magnetohydrodynamics (GRMHD) system with 42 : * divergence cleaning coupled with electron fraction. 43 : * 44 : * We adopt the standard 3+1 form of metric 45 : * \f{align*} 46 : * ds^2 = -\alpha^2 dt^2 + \gamma_{ij}(dx^i + \beta^i dt)(dx^j + \beta^j dt) 47 : * \f} 48 : * where \f$\alpha\f$ is the lapse, \f$\beta^i\f$ is the shift vector, and 49 : * \f$\gamma_{ij}\f$ is the spatial metric. 50 : * 51 : * Primitive variables of the system are 52 : * - rest mass density \f$\rho\f$ 53 : * - electron fraction \f$Y_e\f$ 54 : * - the spatial velocity \f{align*} v^i = 55 : * \frac{1}{\alpha}\left(\frac{u^i}{u^0} + \beta^i\right) \f} measured by the 56 : * Eulerian observer 57 : * - specific internal energy density \f$\epsilon\f$ 58 : * - magnetic field \f$B^i = -^*F^{ia}n_a = -\alpha (^*F^{0i})\f$ 59 : * - the divergence cleaning field \f$\Phi\f$ 60 : * 61 : * with corresponding derived physical quantities which frequently appear in 62 : * equations: 63 : * - The transport velocity \f$v^i_{tr} = \alpha v^i - \beta^i\f$ 64 : * - The Lorentz factor \f{align*} W = -u^an_a = \alpha u^0 = 65 : * \frac{1}{\sqrt{1-\gamma_{ij}v^iv^j}} = \sqrt{1+\gamma^{ij}u_iu_j} = 66 : * \sqrt{1+\gamma^{ij}W^2v_iv_j} \f} 67 : * - The specific enthalpy \f$h = 1 + \epsilon + p/\rho\f$ where \f$p\f$ is the 68 : * pressure specified by a particular equation of state (EoS) 69 : * - The comoving magnetic field \f$b^b = -^*F^{ba} u_a\f$ in the ideal MHD 70 : * limit 71 : * \f{align*} 72 : * b^0 & = W B^i v_i / \alpha \\ 73 : * b^i & = (B^i + \alpha b^0 u^i) / W \\ 74 : * b^2 & = b^ab_a = B^2/W^2 + (B^iv_i)^2 75 : * \f} 76 : * - Augmented enthalpy density \f$(\rho h)^* = \rho h + b^2\f$ and augmented 77 : * pressure \f$p^* = p + b^2/2\f$ which include contributions from magnetic 78 : * field 79 : * 80 : * \note We are using the electromagnetic variables with the scaling convention 81 : * that the factor \f$4\pi\f$ does not appear in Maxwell equations and the 82 : * stress-energy tensor of the EM fields (geometrized Heaviside-Lorentz units). 83 : * To recover the physical value of magnetic field in the usual CGS Gaussian 84 : * unit, the conversion factor is 85 : * \f{align*} 86 : * \sqrt{4\pi}\frac{c^4}{G^{3/2}M_\odot} \approx 8.35 \times 10^{19} 87 : * \,\text{Gauss} 88 : * \f} 89 : * For example, magnetic field $10^{-5}$ with the code unit corresponds to the 90 : * $8.35\times 10^{14}\,\text{G}$ in the CGS Gaussian unit. See also 91 : * documentation of hydro::units::cgs::gauss_unit for details. 92 : * 93 : * The GRMHD equations can be written in a flux-balanced form 94 : * \f[ 95 : * \partial_t U+ \partial_i F^i(U) = S(U). 96 : * \f] 97 : * 98 : * Evolved (conserved) variables \f$U\f$ are 99 : * \f{align*}{ 100 : * U = \sqrt{\gamma}\left[\,\begin{matrix} 101 : * D \\ 102 : * D Y_e \\ 103 : * S_j \\ 104 : * \tau \\ 105 : * B^j \\ 106 : * \Phi \\ 107 : * \end{matrix}\,\right] \equiv \left[\,\,\begin{matrix} 108 : * \tilde{D} \\ 109 : * \tilde{Y_e} \\ 110 : * \tilde{S}_j \\ 111 : * \tilde{\tau} \\ 112 : * \tilde{B}^j \\ 113 : * \tilde{\Phi} \\ 114 : * \end{matrix}\,\,\right] = \sqrt{\gamma} \left[\,\,\begin{matrix} 115 : * \rho W \\ 116 : * \rho W Y_e \\ 117 : * (\rho h)^* W^2 v_j - \alpha b^0 b_j \\ 118 : * (\rho h)^* W^2 - p^* - (\alpha b^0)^2 - \rho W \\ 119 : * B^j \\ 120 : * \Phi \\ 121 : * \end{matrix}\,\,\right] 122 : * \f} 123 : * 124 : * where \f${\tilde D}\f$, \f${\tilde Y}_e\f$,\f${\tilde S}_i\f$, \f${\tilde 125 : * \tau}\f$, \f${\tilde B}^i\f$, and \f${\tilde \Phi}\f$ are a generalized 126 : * mass-energy density, electron fraction, momentum density, specific internal 127 : * energy density, magnetic field, and divergence cleaning field. Also, 128 : * \f$\gamma\f$ is the determinant of the spatial metric \f$\gamma_{ij}\f$. 129 : * 130 : * Corresponding fluxes \f$F^i(U)\f$ are 131 : * 132 : * \f{align*} 133 : * F^i({\tilde D}) &= {\tilde D} v^i_{tr} \\ 134 : * F^i({\tilde Y}_e) &= {\tilde Y}_e v^i_{tr} \\ 135 : * F^i({\tilde S}_j) &= {\tilde S}_j v^i_{tr} + \alpha \sqrt{\gamma} p^* 136 : * \delta^i_j - \alpha b_j \tilde{B}^i / W \\ 137 : * F^i({\tilde \tau}) &= {\tilde \tau} v^i_{tr} + \alpha \sqrt{\gamma} p^* v^i 138 : * - \alpha^2 b^0 \tilde{B}^i / W \\ 139 : * F^i({\tilde B}^j) &= {\tilde B}^j v^i_{tr} - v^j_{tr} {\tilde B}^i 140 : * + \alpha \gamma^{ij} {\tilde \Phi} \\ 141 : * F^i({\tilde \Phi}) &= \alpha {\tilde B^i} - {\tilde \Phi} \beta^i 142 : * \f} 143 : * 144 : * and source terms \f$S(U)\f$ are 145 : * 146 : * \f{align*} 147 : * S({\tilde D}) &= 0 \\ 148 : * S({\tilde Y}_e) &= 0 \\ 149 : * S({\tilde S}_j) &= \frac{\alpha}{2} {\tilde S}^{mn} \partial_j \gamma_{mn} 150 : * + {\tilde S}_k \partial_j \beta^k - ({\tilde D} 151 : * + {\tilde \tau}) \partial_j \alpha \\ 152 : * S({\tilde \tau}) &= \alpha {\tilde S}^{mn} K_{mn} 153 : * - {\tilde S}^m \partial_m \alpha \\ 154 : * S({\tilde B}^j) &= \Phi \partial_k (\alpha \sqrt{\gamma}\gamma^{jk}) \\ 155 : * S({\tilde \Phi}) &= {\tilde B}^k \partial_k \alpha - \alpha K 156 : * {\tilde \Phi} - \alpha \kappa {\tilde \Phi} 157 : * \f} 158 : * 159 : * with 160 : * \f{align*} 161 : * {\tilde S}^{ij} = & \sqrt{\gamma} \left[ (\rho h)^* W^2 v^i v^j + p^* 162 : * \gamma^{ij} - \gamma^{ik}\gamma^{jl}b_kb_l \right] 163 : * \f} 164 : * 165 : * where \f$K\f$ is the trace of the extrinsic curvature \f$K_{ij}\f$, and 166 : * \f$\kappa\f$ is a damping parameter that damps violations of the 167 : * divergence-free (no-monopole) condition \f$\Phi = \partial_i {\tilde B}^i = 168 : * 0\f$. 169 : * 170 : * \note The fluxes and sources of ${\tilde B}$ follow eq. 32 of 171 : * \cite Palenzuela2018 rather than eq. 87 of \cite Moesta2014 172 : * 173 : * \note On the electron fraction side, the source term is currently set to 174 : * \f$S(\tilde{Y}_e) = 0\f$. Implementing the source term using neutrino scheme 175 : * is in progress (Last update : Oct 2022). 176 : * 177 : */ 178 1 : struct System { 179 0 : static constexpr bool is_in_flux_conservative_form = true; 180 0 : static constexpr bool has_primitive_and_conservative_vars = true; 181 0 : static constexpr size_t volume_dim = 3; 182 : 183 0 : using boundary_conditions_base = BoundaryConditions::BoundaryCondition; 184 0 : using boundary_correction_base = BoundaryCorrections::BoundaryCorrection; 185 : 186 0 : using variables_tag = ::Tags::Variables< 187 : tmpl::list<Tags::TildeD, Tags::TildeYe, Tags::TildeTau, Tags::TildeS<>, 188 : Tags::TildeB<>, Tags::TildePhi>>; 189 0 : using flux_variables = 190 : tmpl::list<Tags::TildeD, Tags::TildeYe, Tags::TildeTau, Tags::TildeS<>, 191 : Tags::TildeB<>, Tags::TildePhi>; 192 0 : using non_conservative_variables = tmpl::list<>; 193 0 : using gradient_variables = tmpl::list<>; 194 0 : using primitive_variables_tag = 195 : ::Tags::Variables<hydro::grmhd_tags<DataVector>>; 196 0 : using spacetime_variables_tag = 197 : ::Tags::Variables<gr::tags_for_hydro<volume_dim, DataVector>>; 198 0 : using flux_spacetime_variables_tag = ::Tags::Variables< 199 : tmpl::list<gr::Tags::Lapse<DataVector>, gr::Tags::Shift<DataVector, 3>, 200 : gr::Tags::SpatialMetric<DataVector, 3>, 201 : gr::Tags::SqrtDetSpatialMetric<DataVector>, 202 : gr::Tags::InverseSpatialMetric<DataVector, 3>>>; 203 : 204 0 : using compute_volume_time_derivative_terms = TimeDerivativeTerms; 205 : 206 0 : using conservative_from_primitive = ConservativeFromPrimitive; 207 : template <typename OrderedListOfPrimitiveRecoverySchemes> 208 0 : using primitive_from_conservative = 209 : PrimitiveFromConservative<OrderedListOfPrimitiveRecoverySchemes>; 210 : 211 0 : using compute_largest_characteristic_speed = 212 : Tags::ComputeLargestCharacteristicSpeed; 213 : 214 0 : using inverse_spatial_metric_tag = 215 : gr::Tags::InverseSpatialMetric<DataVector, volume_dim>; 216 : }; 217 : } // namespace ValenciaDivClean 218 : } // namespace grmhd