Line data    Source code

       1           0 : // Distributed under the MIT License.
       2             : // See LICENSE.txt for details.
       3             : 
       4             : #pragma once
       5             : 
       6             : #include "DataStructures/Tensor/TypeAliases.hpp"  // IWYU pragma: keep
       7             : #include "Evolution/Systems/GrMhd/ValenciaDivClean/TagsDeclarations.hpp"  // IWYU pragma: keep
       8             : #include "NumericalAlgorithms/LinearOperators/PartialDerivatives.hpp"  //  IWYU pragma: keep
       9             : #include "PointwiseFunctions/GeneralRelativity/TagsDeclarations.hpp"  // IWYU pragma: keep
      10             : #include "PointwiseFunctions/Hydro/Tags.hpp"
      11             : #include "PointwiseFunctions/Hydro/TagsDeclarations.hpp"  // IWYU pragma: keep
      12             : #include "Utilities/TMPL.hpp"
      13             : 
      14             : // IWYU pragma: no_include "Evolution/Systems/GrMhd/ValenciaDivClean/Tags.hpp"
      15             : // IWYU pragma: no_include "PointwiseFunctions/GeneralRelativity/Tags.hpp"
      16             : // IWYU pragma: no_include "PointwiseFunctions/Hydro/Tags.hpp"
      17             : // IWYU pragma: no_include "Utilities/Gsl.hpp"
      18             : 
      19             : // IWYU pragma: no_forward_declare Tensor
      20             : // IWYU pragma: no_forward_declare Tags::deriv
      21             : 
      22             : /// \cond
      23             : class DataVector;
      24             : 
      25             : namespace Tags {
      26             : template <typename>
      27             : struct Source;
      28             : }  // namespace Tags
      29             : /// \endcond
      30             : 
      31             : namespace grmhd {
      32             : namespace ValenciaDivClean {
      33             : namespace detail {
      34             : void sources_impl(
      35             :     gsl::not_null<Scalar<DataVector>*> source_tilde_tau,
      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             :     gsl::not_null<Scalar<DataVector>*> source_tilde_phi,
      39             : 
      40             :     gsl::not_null<tnsr::I<DataVector, 3, Frame::Inertial>*> tilde_s_up,
      41             :     gsl::not_null<tnsr::II<DataVector, 3, Frame::Inertial>*> densitized_stress,
      42             :     gsl::not_null<Scalar<DataVector>*> h_rho_w_squared_plus_b_squared,
      43             : 
      44             :     const Scalar<DataVector>& magnetic_field_dot_spatial_velocity,
      45             :     const Scalar<DataVector>& magnetic_field_squared,
      46             :     const Scalar<DataVector>& one_over_w_squared,
      47             :     const Scalar<DataVector>& pressure_star,
      48             :     const tnsr::I<DataVector, 3, Frame::Inertial>&
      49             :         trace_spatial_christoffel_second,
      50             : 
      51             :     const Scalar<DataVector>& tilde_d, const Scalar<DataVector>& tilde_ye,
      52             :     const Scalar<DataVector>& tilde_tau,
      53             :     const tnsr::i<DataVector, 3, Frame::Inertial>& tilde_s,
      54             :     const tnsr::I<DataVector, 3, Frame::Inertial>& tilde_b,
      55             :     const Scalar<DataVector>& tilde_phi, const Scalar<DataVector>& lapse,
      56             :     const Scalar<DataVector>& sqrt_det_spatial_metric,
      57             :     const tnsr::II<DataVector, 3, Frame::Inertial>& inv_spatial_metric,
      58             :     const tnsr::i<DataVector, 3, Frame::Inertial>& d_lapse,
      59             :     const tnsr::iJ<DataVector, 3, Frame::Inertial>& d_shift,
      60             :     const tnsr::ijj<DataVector, 3, Frame::Inertial>& d_spatial_metric,
      61             :     const tnsr::I<DataVector, 3, Frame::Inertial>& spatial_velocity,
      62             :     const Scalar<DataVector>& lorentz_factor,
      63             :     const tnsr::I<DataVector, 3, Frame::Inertial>& magnetic_field,
      64             : 
      65             :     const Scalar<DataVector>& rest_mass_density,
      66             :     const Scalar<DataVector>& electron_fraction,
      67             :     const Scalar<DataVector>& pressure,
      68             :     const Scalar<DataVector>& specific_internal_energy,
      69             :     const tnsr::ii<DataVector, 3, Frame::Inertial>& extrinsic_curvature,
      70             :     double constraint_damping_parameter);
      71             : }  // namespace detail
      72             : 
      73             : /*!
      74             :  * \brief Compute the source terms for the flux-conservative Valencia
      75             :  * formulation of GRMHD with divergence cleaning, coupled with electron
      76             :  * fraction.
      77             :  *
      78             :  * A flux-conservative system has the generic form:
      79             :  * \f[
      80             :  * \partial_t U_i + \partial_m F^m(U_i) = S(U_i)
      81             :  * \f]
      82             :  *
      83             :  * where \f$F^a()\f$ denotes the flux of a conserved variable \f$U_i\f$ and
      84             :  * \f$S()\f$ denotes the source term for the conserved variable.
      85             :  *
      86             :  * For the Valencia formulation:
      87             :  * \f{align*}
      88             :  * S({\tilde D}) = & 0\\
      89             :  * S({\tilde S}_i) = & \frac{1}{2} \alpha {\tilde S}^{mn} \partial_i \gamma_{mn}
      90             :  * + {\tilde S}_m \partial_i \beta^m - ({\tilde D} + {\tilde \tau}) \partial_i
      91             :  * \alpha \\
      92             :  * S({\tilde \tau}) = & \alpha {\tilde S}^{mn} K_{mn} - {\tilde S}^m \partial_m
      93             :  * \alpha \\
      94             :  * S({\tilde B}^i) = & -{\tilde B}^m \partial_m \beta^i + {\tilde \Phi}
      95             :  * \gamma^{im} \partial_m \alpha + \alpha {\tilde \Phi} \left( \frac{1}{2}
      96             :  * \gamma^{il} \gamma^{jk} - \gamma^{ij} \gamma^{lk} \right) \partial_l
      97             :  * \gamma_{jk} \\ S({\tilde \Phi}) = & {\tilde B}^k \partial_k \alpha - \alpha K
      98             :  * {\tilde \Phi}
      99             :  * - \alpha \kappa {\tilde \Phi}
     100             :  * \f}
     101             :  *
     102             :  * where
     103             :  * \f{align*}
     104             :  * {\tilde S}^i = & {\tilde S}_m \gamma^{im} \\
     105             :  * {\tilde S}^{ij} = & \sqrt{\gamma} \left[ \left(h \rho W^2 + B^n B_n \right)
     106             :  * v^i v^j + \left(p + p_m \right) \gamma^{ij} - B^n v_n \left( B^i v^j + B^j
     107             :  * v^i \right) - \frac{B^i B^j}{W^2} \right]
     108             :  * \f}
     109             :  *
     110             :  * where \f${\tilde D}\f$, \f${\tilde S}_i\f$, \f${\tilde \tau}\f$, \f${\tilde
     111             :  * B}^i\f$, and \f${\tilde \Phi}\f$ are a generalized mass-energy density,
     112             :  * momentum density, specific internal energy density, magnetic field, and
     113             :  * divergence cleaning field.  Furthermore, \f$\gamma\f$ is the determinant of
     114             :  * the spatial metric \f$\gamma_{ij}\f$, \f$\rho\f$ is the rest mass density,
     115             :  * \f$W\f$ is the Lorentz factor, \f$h\f$ is the specific enthalpy, \f$v^i\f$ is
     116             :  * the spatial velocity, \f$B^k\f$ is the magnetic field, \f$p\f$ is the
     117             :  * pressure, \f$p_m = \frac{1}{2} \left[ \left( B^n v_n \right)^2 + B^n B_n /
     118             :  * W^2 \right]\f$ is the magnetic pressure, \f$\alpha\f$ is the lapse,
     119             :  * \f$\beta^i\f$ is the shift, \f$K\f$ is the trace of the extrinsic curvature
     120             :  * \f$K_{ij}\f$, and \f$\kappa\f$ is a damping parameter that damps violations
     121             :  * of the divergence-free (no-monopole) condition \f$\Phi = \partial_i {\tilde
     122             :  * B}^i = 0\f$ .
     123             :  *
     124             :  * \note For the electron fraction side, the source term is currently set to
     125             :  * \f$S(\tilde{Y}_e) = 0\f$ where the conserved variable \f$\tilde{Y}_e\f$ is a
     126             :  * generalized electron fraction. Implementing the source term using neutrino
     127             :  * scheme is in progress (Last update : Oct 2022).
     128             :  *
     129             :  */
     130           1 : struct ComputeSources {
     131           0 :   using return_tags =
     132             :       tmpl::list<::Tags::Source<grmhd::ValenciaDivClean::Tags::TildeTau>,
     133             :                  ::Tags::Source<grmhd::ValenciaDivClean::Tags::TildeS<>>,
     134             :                  ::Tags::Source<grmhd::ValenciaDivClean::Tags::TildeB<>>,
     135             :                  ::Tags::Source<grmhd::ValenciaDivClean::Tags::TildePhi>>;
     136             : 
     137           0 :   using argument_tags =
     138             :       tmpl::list<grmhd::ValenciaDivClean::Tags::TildeD,
     139             :                  grmhd::ValenciaDivClean::Tags::TildeYe,
     140             :                  grmhd::ValenciaDivClean::Tags::TildeTau,
     141             :                  grmhd::ValenciaDivClean::Tags::TildeS<>,
     142             :                  grmhd::ValenciaDivClean::Tags::TildeB<>,
     143             :                  grmhd::ValenciaDivClean::Tags::TildePhi,
     144             :                  hydro::Tags::SpatialVelocity<DataVector, 3>,
     145             :                  hydro::Tags::MagneticField<DataVector, 3>,
     146             :                  hydro::Tags::RestMassDensity<DataVector>,
     147             :                  hydro::Tags::ElectronFraction<DataVector>,
     148             :                  hydro::Tags::SpecificInternalEnergy<DataVector>,
     149             :                  hydro::Tags::LorentzFactor<DataVector>,
     150             :                  hydro::Tags::Pressure<DataVector>, gr::Tags::Lapse<DataVector>,
     151             :                  ::Tags::deriv<gr::Tags::Lapse<DataVector>, tmpl::size_t<3>,
     152             :                                Frame::Inertial>,
     153             :                  ::Tags::deriv<gr::Tags::Shift<DataVector, 3>, tmpl::size_t<3>,
     154             :                                Frame::Inertial>,
     155             :                  gr::Tags::SpatialMetric<DataVector, 3>,
     156             :                  ::Tags::deriv<gr::Tags::SpatialMetric<DataVector, 3>,
     157             :                                tmpl::size_t<3>, Frame::Inertial>,
     158             :                  gr::Tags::InverseSpatialMetric<DataVector, 3>,
     159             :                  gr::Tags::SqrtDetSpatialMetric<DataVector>,
     160             :                  gr::Tags::ExtrinsicCurvature<DataVector, 3>,
     161             :                  grmhd::ValenciaDivClean::Tags::ConstraintDampingParameter>;
     162             : 
     163           0 :   static void apply(
     164             :       gsl::not_null<Scalar<DataVector>*> source_tilde_tau,
     165             :       gsl::not_null<tnsr::i<DataVector, 3, Frame::Inertial>*> source_tilde_s,
     166             :       gsl::not_null<tnsr::I<DataVector, 3, Frame::Inertial>*> source_tilde_b,
     167             :       gsl::not_null<Scalar<DataVector>*> source_tilde_phi,
     168             :       const Scalar<DataVector>& tilde_d, const Scalar<DataVector>& tilde_ye,
     169             :       const Scalar<DataVector>& tilde_tau,
     170             :       const tnsr::i<DataVector, 3, Frame::Inertial>& tilde_s,
     171             :       const tnsr::I<DataVector, 3, Frame::Inertial>& tilde_b,
     172             :       const Scalar<DataVector>& tilde_phi,
     173             :       const tnsr::I<DataVector, 3, Frame::Inertial>& spatial_velocity,
     174             :       const tnsr::I<DataVector, 3, Frame::Inertial>& magnetic_field,
     175             :       const Scalar<DataVector>& rest_mass_density,
     176             :       const Scalar<DataVector>& electron_fraction,
     177             :       const Scalar<DataVector>& specific_internal_energy,
     178             :       const Scalar<DataVector>& lorentz_factor,
     179             :       const Scalar<DataVector>& pressure, const Scalar<DataVector>& lapse,
     180             :       const tnsr::i<DataVector, 3, Frame::Inertial>& d_lapse,
     181             :       const tnsr::iJ<DataVector, 3, Frame::Inertial>& d_shift,
     182             :       const tnsr::ii<DataVector, 3, Frame::Inertial>& spatial_metric,
     183             :       const tnsr::ijj<DataVector, 3, Frame::Inertial>& d_spatial_metric,
     184             :       const tnsr::II<DataVector, 3, Frame::Inertial>& inv_spatial_metric,
     185             :       const Scalar<DataVector>& sqrt_det_spatial_metric,
     186             :       const tnsr::ii<DataVector, 3, Frame::Inertial>& extrinsic_curvature,
     187             :       double constraint_damping_parameter);
     188             : };
     189             : }  // namespace ValenciaDivClean
     190             : }  // namespace grmhd