SpECTRE Documentation Coverage Report
Current view: top level - Evolution/Systems/GeneralizedHarmonic/GaugeSourceFunctions - DampedHarmonic.hpp Hit Total Coverage
Commit: 6bd9cb26e64b74498abb59750d6afe6a33bff976 Lines: 6 28 21.4 %
Date: 2024-05-21 21:00:13
Legend: Lines: hit not hit

          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             : #include <limits>
       8             : 
       9             : #include "DataStructures/DataBox/Tag.hpp"
      10             : #include "DataStructures/Tensor/TypeAliases.hpp"
      11             : #include "Domain/Tags.hpp"
      12             : #include "Evolution/Systems/GeneralizedHarmonic/GaugeSourceFunctions/Gauges.hpp"
      13             : #include "Evolution/Systems/GeneralizedHarmonic/Tags.hpp"
      14             : #include "PointwiseFunctions/GeneralRelativity/TagsDeclarations.hpp"
      15             : #include "Utilities/ContainerHelpers.hpp"
      16             : #include "Utilities/TMPL.hpp"
      17             : 
      18             : /// \cond
      19             : class DataVector;
      20             : namespace gsl {
      21             : template <class T>
      22             : class not_null;
      23             : }  // namespace gsl
      24             : /// \endcond
      25             : 
      26             : namespace gh {
      27             : namespace gauges {
      28             : /*!
      29             :  * \brief Damped harmonic gauge source function and its spacetime derivative.
      30             :  *
      31             :  * \details The gauge condition has been taken from \cite Szilagyi2009qz and
      32             :  * \cite Deppe2018uye. We provide both a "rollon" version
      33             :  * (`damped_harmonic_rollon`), and a "non-rollon" version (`damped_harmonic`).
      34             :  * In the non-rollon version the rollon function \f$R(t)=1\f$.
      35             :  *
      36             :  * \warning Only the non-rollon version can be used with a moving mesh.
      37             :  *
      38             :  * The covariant form of the source function \f$H_a\f$ is written as:
      39             :  *
      40             :  * \f{align*}
      41             :  * H_a :=  [1 - R(t)] H_a^\mathrm{init} +
      42             :  *  [\mu_{L1} \mathrm{log}(\sqrt{g}/N) + \mu_{L2} \mathrm{log}(1/N)] t_a
      43             :  *   - \mu_S g_{ai} N^i / N
      44             :  * \f}
      45             :  *
      46             :  * where \f$N, N^k\f$ are the lapse and shift respectively, \f$t_a\f$ is the
      47             :  * unit normal one-form to the spatial slice, and \f$g_{ab}\f$ is
      48             :  * the spatial metric (obtained by projecting the spacetime metric onto the
      49             :  * 3-slice, i.e. \f$g_{ab} = \psi_{ab} + t_a t_b\f$). The prefactors are:
      50             :  *
      51             :  * \f{align*}
      52             :  *  \mu_{L1} &= A_{L1} R(t) W(x^i) \mathrm{log}(\sqrt{g}/N)^{e_{L1}}, \\
      53             :  *  \mu_{L2} &= A_{L2} R(t) W(x^i) \mathrm{log}(1/N)^{e_{L2}}, \\
      54             :  *  \mu_{S} &= A_{S} R(t) W(x^i) \mathrm{log}(\sqrt{g}/N)^{e_{S}},
      55             :  * \f}
      56             :  *
      57             :  * temporal roll-on function \f$ R(t)\f$ is:
      58             :  *
      59             :  * \f{align*}
      60             :  * \begin{array}{ll}
      61             :  *     R(t) & = 0, & t< t_{0} \\
      62             :  *          & = 1 - \exp[-((t - t_{0})/ \sigma_t)^4], & t\geq t_{0} \\
      63             :  * \end{array}
      64             :  * \f}
      65             :  *
      66             :  * and the spatial weight function is:
      67             :  *
      68             :  * \f{align*}
      69             :  * W(x^i) = \exp[-(r/\sigma_r)^2].
      70             :  * \f}
      71             :  *
      72             :  * This weight function can be written with multiple constant factors in the
      73             :  * exponent in literature \cite Deppe2018uye, but we absorb them all into
      74             :  * \f$ \sigma_r\f$ here. The coordinate \f$ r\f$ is the Euclidean radius
      75             :  * in Inertial coordinates.
      76             :  *
      77             :  * Note that for the last three terms in \f$H_a\f$ (with \f$ X = \{L1, L2, S\}
      78             :  * \f$):
      79             :  *   - Amplitude factors \f$ A_{X} \f$ are taken as input here as `amp_coef_X`
      80             :  *   - Exponents \f$ e_X\f$ are taken as input here as `exp_X`.
      81             :  *   - Spatial weight function \f$W\f$ is specified completely by
      82             :  *     \f$\sigma_r\f$, which is taken as input here as `sigma_r`.
      83             :  *
      84             :  * Also computes spacetime derivatives, i.e. \f$\partial_a H_b\f$, of the damped
      85             :  * harmonic source function H. Using notation from damped_harmonic_h(), we
      86             :  * rewrite the same as:
      87             :  *
      88             :  * \f{align*}
      89             :  * \partial_a H_b =& \partial_a T_1 + \partial_a T_2 + \partial_a T_3, \\
      90             :  * H_a =& T_1 + T_2 + T_3,
      91             :  * \f}
      92             :  *
      93             :  * where:
      94             :  *
      95             :  * \f{align*}
      96             :  * T_1 =& [1 - R(t)] H_a^\mathrm{init}, \\
      97             :  * T_2 =& [\mu_{L1} \mathrm{log}(\sqrt{g}/N) + \mu_{L2} \mathrm{log}(1/N)] t_a,
      98             :  * \\
      99             :  * T_3 =& - \mu_S g_{ai} N^i / N.
     100             :  * \f}
     101             :  *
     102             :  * Derivation:
     103             :  *
     104             :  * \f$\blacksquare\f$ For \f$ T_1 \f$, the derivatives are:
     105             :  * \f{align*}
     106             :  * \partial_a T_1 = (1 - R(t))
     107             :  * \partial_a H_b^\mathrm{init}
     108             :  *                - H_b^\mathrm{init} \partial_a R.
     109             :  * \f}
     110             :  *
     111             :  * \f$\blacksquare\f$ Write \f$ T_2 \equiv (\mu_1 + \mu_2) t_b \f$. Then:
     112             :  *
     113             :  * \f{align*}
     114             :  * \partial_a T_2 =& (\partial_a \mu_1 + \partial_a \mu_2) t_b \\
     115             :  *               +& (\mu_1 + \mu_2) \partial_a t_b,
     116             :  * \f}
     117             :  *
     118             :  * where
     119             :  *
     120             :  * \f{align*}
     121             :  * \partial_a t_b =& \left(-\partial_a N, 0, 0, 0\right) \\
     122             :  *
     123             :  * \partial_a \mu_1
     124             :  *  =& \partial_a [A_{L1} R(t) W(x^i) \mathrm{log}(\sqrt{g}/N)^{e_{L1} +
     125             :  * 1}], \\
     126             :  *  =& A_{L1} R(t) W(x^i) \partial_a [\mathrm{log}(\sqrt{g}/N)^{e_{L1} +
     127             :  * 1}] \\
     128             :  *   +& A_{L1} \mathrm{log}(\sqrt{g}/N)^{e_{L1} + 1} \partial_a [R(t)
     129             :  * W(x^i)],\\
     130             :  *
     131             :  * \partial_a \mu_2
     132             :  *  =& \partial_a [A_{L2} R(t) W(x^i) \mathrm{log}(1/N)^{e_{L2} + 1}], \\
     133             :  *  =& A_{L2} R(t) W(x^i) \partial_a [\mathrm{log}(1/N)^{e_{L2} + 1}] \\
     134             :  *     +& A_{L2} \mathrm{log}(1/N)^{e_{L2} + 1} \partial_a [R(t) W(x^i)],
     135             :  * \f}
     136             :  *
     137             :  * where \f$\partial_a [R W] = \left(\partial_0 R(t), \partial_i
     138             :  * W(x^j)\right)\f$.
     139             :  *
     140             :  * \f$\blacksquare\f$ Finally, the derivatives of \f$ T_3 \f$ are:
     141             :  *
     142             :  * \f[
     143             :  * \partial_a T_3 = -\partial_a(\mu_S/N) g_{bi} N^i
     144             :  *                  -(\mu_S/N) \partial_a(g_{bi}) N^i
     145             :  *                  -(\mu_S/N) g_{bi}\partial_a N^i,
     146             :  * \f]
     147             :  *
     148             :  * where
     149             :  *
     150             :  * \f{align*}
     151             :  * \partial_a(\mu_S / N) =& (1/N)\partial_a \mu_S
     152             :  *                       - \frac{\mu_S}{N^2}\partial_a N, \,\,\mathrm{and}\\
     153             :  * \partial_a \mu_S =& \partial_a [A_S R(t) W(x^i)
     154             :  * \mathrm{log}(\sqrt{g}/N)^{e_S}], \\
     155             :  *                  =& A_S R(t) W(x^i) \partial_a
     156             :  * [\mathrm{log}(\sqrt{g}/N)^{e_S}] \\
     157             :  *                  +& A_S \mathrm{log}(\sqrt{g} / N)^{e_S} \partial_a [R(t)
     158             :  * W(x^i)].
     159             :  * \f}
     160             :  */
     161             : template <size_t SpatialDim, typename Frame>
     162           1 : void damped_harmonic_rollon(
     163             :     gsl::not_null<tnsr::a<DataVector, SpatialDim, Frame>*> gauge_h,
     164             :     gsl::not_null<tnsr::ab<DataVector, SpatialDim, Frame>*> d4_gauge_h,
     165             :     const tnsr::a<DataVector, SpatialDim, Frame>& gauge_h_init,
     166             :     const tnsr::ab<DataVector, SpatialDim, Frame>& dgauge_h_init,
     167             :     const Scalar<DataVector>& lapse,
     168             :     const tnsr::I<DataVector, SpatialDim, Frame>& shift,
     169             :     const Scalar<DataVector>& sqrt_det_spatial_metric,
     170             :     const tnsr::II<DataVector, SpatialDim, Frame>& inverse_spatial_metric,
     171             :     const tnsr::abb<DataVector, SpatialDim, Frame>& d4_spacetime_metric,
     172             :     const Scalar<DataVector>& half_pi_two_normals,
     173             :     const tnsr::i<DataVector, SpatialDim, Frame>& half_phi_two_normals,
     174             :     const tnsr::aa<DataVector, SpatialDim, Frame>& spacetime_metric,
     175             :     const tnsr::iaa<DataVector, SpatialDim, Frame>& phi, double time,
     176             :     const tnsr::I<DataVector, SpatialDim, Frame>& coords, double amp_coef_L1,
     177             :     double amp_coef_L2, double amp_coef_S, int exp_L1, int exp_L2, int exp_S,
     178             :     double rollon_start_time, double rollon_width, double sigma_r);
     179             : 
     180             : /*!
     181             :  * \copydoc damped_harmonic_rollon()
     182             :  */
     183             : template <size_t SpatialDim, typename Frame>
     184           1 : void damped_harmonic(
     185             :     gsl::not_null<tnsr::a<DataVector, SpatialDim, Frame>*> gauge_h,
     186             :     gsl::not_null<tnsr::ab<DataVector, SpatialDim, Frame>*> d4_gauge_h,
     187             :     const Scalar<DataVector>& lapse,
     188             :     const tnsr::I<DataVector, SpatialDim, Frame>& shift,
     189             :     const Scalar<DataVector>& sqrt_det_spatial_metric,
     190             :     const tnsr::II<DataVector, SpatialDim, Frame>& inverse_spatial_metric,
     191             :     const tnsr::abb<DataVector, SpatialDim, Frame>& d4_spacetime_metric,
     192             :     const Scalar<DataVector>& half_pi_two_normals,
     193             :     const tnsr::i<DataVector, SpatialDim, Frame>& half_phi_two_normals,
     194             :     const tnsr::aa<DataVector, SpatialDim, Frame>& spacetime_metric,
     195             :     const tnsr::iaa<DataVector, SpatialDim, Frame>& phi,
     196             :     const tnsr::I<DataVector, SpatialDim, Frame>& coords, double amp_coef_L1,
     197             :     double amp_coef_L2, double amp_coef_S, int exp_L1, int exp_L2, int exp_S,
     198             :     double sigma_r);
     199             : 
     200             : /*!
     201             :  * \brief Impose damped harmonic gauge.
     202             :  *
     203             :  * \see `damped_harmonic()`
     204             :  */
     205           1 : class DampedHarmonic final : public GaugeCondition {
     206             :  public:
     207             :   /// The width of the Gaussian for the spatial decay of the damped harmonic
     208             :   /// gauge.
     209           1 :   struct SpatialDecayWidth {
     210           0 :     using type = double;
     211           0 :     static constexpr Options::String help{
     212             :         "Spatial width (sigma_r) of weight function (W(x^i)) used in the "
     213             :         "damped harmonic gauge."};
     214             :   };
     215             :   /// The amplitudes for the L1, L2, and S terms, respectively, for the damped
     216             :   /// harmonic gauge.
     217           1 :   struct Amplitudes {
     218           0 :     using type = std::array<double, 3>;
     219           0 :     static constexpr Options::String help{
     220             :         "Amplitudes [A_{L1}, A_{L2}, A_{S}] for the damped harmonic gauge."};
     221             :   };
     222             :   /// The exponents for the L1, L2, and S terms, respectively, for the damped
     223             :   /// harmonic gauge.
     224           1 :   struct Exponents {
     225           0 :     using type = std::array<int, 3>;
     226           0 :     static constexpr Options::String help{
     227             :         "Exponents [e_{L1}, e_{L2}, e_{S}] for the damped harmonic gauge."};
     228             :   };
     229             : 
     230           0 :   static constexpr Options::String help{
     231             :       "Apply damped harmonic/damped wave gauge."};
     232             : 
     233           0 :   using options = tmpl::list<SpatialDecayWidth, Amplitudes, Exponents>;
     234             : 
     235           0 :   DampedHarmonic(double width, const std::array<double, 3>& amps,
     236             :                  const std::array<int, 3>& exps);
     237             : 
     238           0 :   DampedHarmonic() = default;
     239           0 :   DampedHarmonic(const DampedHarmonic&) = default;
     240           0 :   DampedHarmonic& operator=(const DampedHarmonic&) = default;
     241           0 :   DampedHarmonic(DampedHarmonic&&) = default;
     242           0 :   DampedHarmonic& operator=(DampedHarmonic&&) = default;
     243           0 :   ~DampedHarmonic() override = default;
     244             : 
     245             :   /// \cond
     246             :   explicit DampedHarmonic(CkMigrateMessage* msg);
     247             :   using PUP::able::register_constructor;
     248             :   WRAPPED_PUPable_decl_template(DampedHarmonic);  // NOLINT
     249             :   /// \endcond
     250             : 
     251           0 :   std::unique_ptr<GaugeCondition> get_clone() const override;
     252             : 
     253             :   template <size_t SpatialDim>
     254           0 :   void gauge_and_spacetime_derivative(
     255             :       gsl::not_null<tnsr::a<DataVector, SpatialDim, Frame::Inertial>*> gauge_h,
     256             :       gsl::not_null<tnsr::ab<DataVector, SpatialDim, Frame::Inertial>*>
     257             :           d4_gauge_h,
     258             :       const Scalar<DataVector>& lapse,
     259             :       const tnsr::I<DataVector, SpatialDim, Frame::Inertial>& shift,
     260             :       const Scalar<DataVector>& sqrt_det_spatial_metric,
     261             :       const tnsr::II<DataVector, SpatialDim, Frame::Inertial>&
     262             :           inverse_spatial_metric,
     263             :       const tnsr::abb<DataVector, SpatialDim, Frame::Inertial>&
     264             :           d4_spacetime_metric,
     265             :       const Scalar<DataVector>& half_pi_two_normals,
     266             :       const tnsr::i<DataVector, SpatialDim, Frame::Inertial>&
     267             :           half_phi_two_normals,
     268             :       const tnsr::aa<DataVector, SpatialDim, Frame::Inertial>& spacetime_metric,
     269             :       const tnsr::iaa<DataVector, SpatialDim, Frame::Inertial>& phi,
     270             :       double time,
     271             :       const tnsr::I<DataVector, SpatialDim, Frame::Inertial>& inertial_coords)
     272             :       const;
     273             : 
     274             :   // NOLINTNEXTLINE(google-runtime-references)
     275           0 :   void pup(PUP::er& p) override;
     276             : 
     277             :  private:
     278           0 :   double spatial_decay_width_{std::numeric_limits<double>::signaling_NaN()};
     279           0 :   std::array<double, 3> amplitudes_{
     280             :       {std::numeric_limits<double>::signaling_NaN(),
     281             :        std::numeric_limits<double>::signaling_NaN(),
     282             :        std::numeric_limits<double>::signaling_NaN()}};
     283           0 :   std::array<int, 3> exponents_{{std::numeric_limits<int>::max(),
     284             :                                  std::numeric_limits<int>::max(),
     285             :                                  std::numeric_limits<int>::max()}};
     286             : };
     287             : 
     288             : namespace DampedHarmonicGauge_detail {
     289             : // Used in the test
     290             : double roll_on_function(double time, double t_start, double sigma_t);
     291             : 
     292             : double time_deriv_of_roll_on_function(double time, double t_start,
     293             :                                       double sigma_t);
     294             : }  // namespace DampedHarmonicGauge_detail
     295             : }  // namespace gauges
     296             : }  // namespace gh

Generated by: LCOV version 1.14