Line data    Source code

       1           0 : // Distributed under the MIT License.
       2             : // See LICENSE.txt for details.
       3             : 
       4             : #pragma once
       5             : 
       6             : #include <memory>
       7             : #include <optional>
       8             : #include <pup.h>
       9             : #include <string>
      10             : #include <type_traits>
      11             : 
      12             : #include "DataStructures/DataBox/Prefixes.hpp"
      13             : #include "DataStructures/DataVector.hpp"
      14             : #include "DataStructures/Tensor/Tensor.hpp"
      15             : #include "DataStructures/Variables.hpp"
      16             : #include "Evolution/BoundaryConditions/Type.hpp"
      17             : #include "Evolution/Systems/ScalarWave/BoundaryConditions/BoundaryCondition.hpp"
      18             : #include "Evolution/Systems/ScalarWave/Tags.hpp"
      19             : #include "NumericalAlgorithms/LinearOperators/PartialDerivatives.hpp"
      20             : #include "Options/Options.hpp"
      21             : #include "Options/String.hpp"
      22             : #include "PointwiseFunctions/AnalyticData/Tags.hpp"
      23             : #include "PointwiseFunctions/AnalyticSolutions/AnalyticSolution.hpp"
      24             : #include "Utilities/Gsl.hpp"
      25             : #include "Utilities/Serialization/CharmPupable.hpp"
      26             : #include "Utilities/TMPL.hpp"
      27             : 
      28             : /// \cond
      29             : namespace domain::Tags {
      30             : template <size_t Dim, typename Frame>
      31             : struct Coordinates;
      32             : }  // namespace domain::Tags
      33             : /// \endcond
      34             : 
      35             : namespace ScalarWave::BoundaryConditions {
      36             : namespace detail {
      37             : /// The type of spherical radiation boundary condition to impose
      38             : enum class ConstraintPreservingSphericalRadiationType {
      39             :   /// Impose \f$(\partial_t + \partial_r)\Psi=0\f$
      40             :   Sommerfeld,
      41             :   /// Impose \f$(\partial_t + \partial_r + r^{-1})\Psi=0\f$
      42             :   FirstOrderBaylissTurkel,
      43             :   /// Imposes a second-order Bayliss-Turkel boundary condition
      44             :   SecondOrderBaylissTurkel
      45             : };
      46             : 
      47             : ConstraintPreservingSphericalRadiationType
      48             : convert_constraint_preserving_spherical_radiation_type_from_yaml(
      49             :     const Options::Option& options);
      50             : }  // namespace detail
      51             : 
      52             : /*!
      53             :  * \brief Constraint-preserving spherical radiation boundary condition that
      54             :  * seeks to avoid ingoing constraint violations and radiation.
      55             :  *
      56             :  * The constraint-preserving part of the boundary condition imposes the
      57             :  * following condition on the time derivatives of the characteristic fields:
      58             :  *
      59             :  * \f{align*}{
      60             :  *   d_tw^\Psi&\to d_tw^{\Psi}+\lambda_{\Psi}n^i\mathcal{C}_i, \\
      61             :  *   d_tw^0_i&\to d_tw^{0}_i+\lambda_{0}n^jP^k_i\mathcal{C}_{ik},
      62             :  * \f}
      63             :  *
      64             :  * where
      65             :  *
      66             :  * \f{align*}{
      67             :  * P^k{}_i=\delta^k_i-n^kn_i
      68             :  * \f}
      69             :  *
      70             :  * projects a tensor onto the spatial surface to which \f$n_i\f$ is normal, and
      71             :  * \f$d_t w\f$ is the evolved to characteristic field transformation applied to
      72             :  * the time derivatives of the evolved fields. That is,
      73             :  *
      74             :  * \f{align*}{
      75             :  * d_t w^\Psi&=\partial_t \Psi, \\
      76             :  * d_t w_i^0&=(\delta^k_i-n^k n_i)\partial_t \Phi_k, \\
      77             :  * d_t w^{\pm}&=\partial_t\Pi\pm n^k\partial_t\Phi_k - \gamma_2\partial_t\Psi.
      78             :  * \f}
      79             :  *
      80             :  * The constraints are defined as:
      81             :  *
      82             :  * \f{align*}{
      83             :  *  \mathcal{C}_i&=\partial_i\Psi - \Phi_i=0, \\
      84             :  *  \mathcal{C}_{ij}&=\partial_{[i}\Phi_{j]}=0
      85             :  * \f}
      86             :  *
      87             :  * Radiation boundary conditions impose a condition on \f$\Pi\f$ or its time
      88             :  * derivative. We denote the boundary condition value of the time derivative of
      89             :  * \f$\Pi\f$ by \f$\partial_t\Pi^{\mathrm{BC}}\f$. With this, we can impose
      90             :  * boundary conditions on the time derivatives of the evolved variables as
      91             :  * follows:
      92             :  *
      93             :  * \f{align*}{
      94             :  * \partial_{t} \Psi&\to\partial_{t}\Psi +
      95             :  *                    \lambda_\Psi n^i \mathcal{C}_i, \\
      96             :  * \partial_{t}\Pi&\to\partial_{t}\Pi-\left(\partial_t\Pi
      97             :  *                  - \partial_t\Pi^{\mathrm{BC}}\right)
      98             :  *                  +\gamma_2\lambda_\Psi n^i \mathcal{C}_i
      99             :  *                  =\partial_t\Pi^{\mathrm{BC}}
     100             :  *                  +\gamma_2\lambda_\Psi n^i \mathcal{C}_i, \\
     101             :  * \partial_{t}\Phi_i&\to\partial_{t}\Phi_i+ 2 \lambda_0n^j \mathcal{C}_{ji}.
     102             :  * \f}
     103             :  *
     104             :  * Below we assume the normal vector \f$n^i\f$ is the radial unit normal vector.
     105             :  * That is, we assume the outer boundary is spherical. A Sommerfeld
     106             :  * \cite Sommerfeld1949 radiation condition is given by
     107             :  *
     108             :  * \f{align*}{
     109             :  *  \partial_t\Psi=n^i\Phi_i
     110             :  * \f}
     111             :  *
     112             :  * Or, assuming that \f$\partial_tn^i=0\f$ (or is very small),
     113             :  *
     114             :  * \f{align*}{
     115             :  *  \partial_t\Pi^{\mathrm{BC}}=n^i\partial_t\Phi_i
     116             :  * \f}
     117             :  *
     118             :  * The Bayliss-Turkel \cite BaylissTurkel boundary conditions are given by:
     119             :  *
     120             :  * \f{align*}{
     121             :  *  \prod_{l=1}^m\left(\partial_t + \partial_r + \frac{2l-1}{r}\right)\Psi=0
     122             :  * \f}
     123             :  *
     124             :  * The first-order form is
     125             :  *
     126             :  * \f{align*}{
     127             :  *  \partial_t\Pi^{\mathrm{BC}}=n^i\partial_t\Phi_i + \frac{1}{r}\partial_t\Psi,
     128             :  * \f}
     129             :  *
     130             :  * assuming \f$\partial_t n^i=0\f$ and \f$\partial_t r=0\f$.
     131             :  *
     132             :  * The second-order boundary condition is given by,
     133             :  *
     134             :  * \f{align*}{
     135             :  *  \partial_t\Pi^{\mathrm{BC}}
     136             :  *   &=\left(\partial_t\partial_r + \partial_r\partial_t +
     137             :  *     \partial_r^2+\frac{4}{r}\partial_t
     138             :  *     +\frac{4}{r}\partial_r + \frac{2}{r^2}\right)\Psi \notag \\
     139             :  *     &=n^i(\partial_t\Phi_i-\partial_i\Pi) + n^i n^j\partial_i\Phi_j -
     140             :  *     \frac{4}{r}\Pi+\frac{4}{r}n^i\Phi_i + \frac{2}{r^2}\Psi \notag \\
     141             :  *     &=n^i(\gamma_2(\partial_i\Psi - \Phi_i)-2\partial_i\Pi)
     142             :  *     + n^i n^j\partial_i\Phi_j -
     143             :  *     \frac{4}{r}\Pi+\frac{4}{r}n^i\Phi_i + \frac{2}{r^2}\Psi,
     144             :  * \f}
     145             :  *
     146             :  * where \f$\partial_t\f$ is the derivative with respect to the inertial frame
     147             :  * time, and we are assuming \f$\partial_t n^i=0\f$ and \f$\partial_t r=0\f$.
     148             :  *
     149             :  * The moving mesh can be accounted for by using
     150             :  *
     151             :  * \f{align*}{
     152             :  * \partial_t r = \frac{1}{r}x^i\delta_{ij}\partial_t x^j
     153             :  * \f}
     154             :  *
     155             :  * \note It is not clear if \f$\partial_t\Phi_i\f$ should be replaced by
     156             :  * \f$-\partial_i\Pi\f$, which is the evolution equation but without the
     157             :  * constraint.
     158             :  *
     159             :  * \note On a moving mesh the characteristic speeds change according to
     160             :  * \f$\lambda\to\lambda-v^i_gn_i\f$ where \f$v^i_g\f$ is the mesh velocity.
     161             :  *
     162             :  * \note For the scalar wave system \f$\lambda_0 = \lambda_\psi\f$
     163             :  *
     164             :  * \warning The boundary conditions are implemented assuming the outer boundary
     165             :  * is spherical. It might be possible to generalize the condition to
     166             :  * non-spherical boundaries by using \f$x^i/r\f$ instead of \f$n^i\f$, but this
     167             :  * hasn't been tested.
     168             :  *
     169             :  * \warning The received time derivatives are in the logical frame, not the
     170             :  * inertial frame and would need to first be corrected by subtracting the mesh
     171             :  * velocity. Similarly, the time derivative correction is in the logical frame.
     172             :  * We instead substitute the evolution equations directly in since the scalar
     173             :  * wave system is quite simple.
     174             :  */
     175             : template <size_t Dim>
     176           1 : class ConstraintPreservingSphericalRadiation final
     177             :     : public BoundaryCondition<Dim> {
     178             :  public:
     179           0 :   struct TypeOptionTag {
     180           0 :     using type = detail::ConstraintPreservingSphericalRadiationType;
     181           0 :     static std::string name() { return "Type"; }
     182           0 :     static constexpr Options::String help{
     183             :         "Whether to impose Sommerfeld, first-order Bayliss-Turkel, or "
     184             :         "second-order Bayliss-Turkel spherical radiation boundary conditions."};
     185             :   };
     186             : 
     187           0 :   using options = tmpl::list<TypeOptionTag>;
     188           0 :   static constexpr Options::String help{
     189             :       "Constraint-preserving spherical radiation boundary conditions setting "
     190             :       "the time derivatives of Psi, Phi, and Pi to avoid incoming constraint "
     191             :       "violations, and imposing radiation boundary conditions."};
     192             : 
     193           0 :   ConstraintPreservingSphericalRadiation(
     194             :       detail::ConstraintPreservingSphericalRadiationType type);
     195             : 
     196           0 :   ConstraintPreservingSphericalRadiation() = default;
     197             :   /// \cond
     198             :   ConstraintPreservingSphericalRadiation(
     199             :       ConstraintPreservingSphericalRadiation&&) = default;
     200             :   ConstraintPreservingSphericalRadiation& operator=(
     201             :       ConstraintPreservingSphericalRadiation&&) = default;
     202             :   ConstraintPreservingSphericalRadiation(
     203             :       const ConstraintPreservingSphericalRadiation&) = default;
     204             :   ConstraintPreservingSphericalRadiation& operator=(
     205             :       const ConstraintPreservingSphericalRadiation&) = default;
     206             :   /// \endcond
     207           0 :   ~ConstraintPreservingSphericalRadiation() override = default;
     208             : 
     209           0 :   explicit ConstraintPreservingSphericalRadiation(CkMigrateMessage* msg);
     210             : 
     211           0 :   WRAPPED_PUPable_decl_base_template(
     212             :       domain::BoundaryConditions::BoundaryCondition,
     213             :       ConstraintPreservingSphericalRadiation);
     214             : 
     215           0 :   auto get_clone() const -> std::unique_ptr<
     216             :       domain::BoundaryConditions::BoundaryCondition> override;
     217             : 
     218           0 :   static constexpr evolution::BoundaryConditions::Type bc_type =
     219             :       evolution::BoundaryConditions::Type::TimeDerivative;
     220             : 
     221           0 :   void pup(PUP::er& p) override;
     222             : 
     223           0 :   using dg_interior_evolved_variables_tags =
     224             :       tmpl::list<ScalarWave::Tags::Psi, ScalarWave::Tags::Pi,
     225             :                  ScalarWave::Tags::Phi<Dim>>;
     226           0 :   using dg_interior_temporary_tags =
     227             :       tmpl::list<domain::Tags::Coordinates<Dim, Frame::Inertial>,
     228             :                  Tags::ConstraintGamma2>;
     229             :   // using dg_interior_dt_vars_tags = tmpl::list<>;
     230           0 :   using dg_interior_deriv_vars_tags = tmpl::list<
     231             :       ::Tags::deriv<ScalarWave::Tags::Psi, tmpl::size_t<Dim>, Frame::Inertial>,
     232             :       ::Tags::deriv<ScalarWave::Tags::Pi, tmpl::size_t<Dim>, Frame::Inertial>,
     233             :       ::Tags::deriv<ScalarWave::Tags::Phi<Dim>, tmpl::size_t<Dim>,
     234             :                     Frame::Inertial>>;
     235           0 :   using dg_gridless_tags = tmpl::list<>;
     236             : 
     237           0 :   std::optional<std::string> dg_time_derivative(
     238             :       gsl::not_null<Scalar<DataVector>*> dt_psi_correction,
     239             :       gsl::not_null<Scalar<DataVector>*> dt_pi_correction,
     240             :       gsl::not_null<tnsr::i<DataVector, Dim, Frame::Inertial>*>
     241             :           dt_phi_correction,
     242             : 
     243             :       const std::optional<tnsr::I<DataVector, Dim, Frame::Inertial>>&
     244             :           face_mesh_velocity,
     245             :       const tnsr::i<DataVector, Dim, Frame::Inertial>& normal_covector,
     246             : 
     247             :       const Scalar<DataVector>& psi, const Scalar<DataVector>& pi,
     248             :       const tnsr::i<DataVector, Dim, Frame::Inertial>& phi,
     249             : 
     250             :       const tnsr::I<DataVector, Dim, Frame::Inertial>& coords,
     251             :       const Scalar<DataVector>& gamma2,
     252             : 
     253             :       const tnsr::i<DataVector, Dim, Frame::Inertial>& d_psi,
     254             :       const tnsr::i<DataVector, Dim, Frame::Inertial>& d_pi,
     255             :       const tnsr::ij<DataVector, Dim, Frame::Inertial>& d_phi) const;
     256             : 
     257             :  private:
     258           0 :   detail::ConstraintPreservingSphericalRadiationType type_{
     259             :       detail::ConstraintPreservingSphericalRadiationType::Sommerfeld};
     260             : };
     261             : }  // namespace ScalarWave::BoundaryConditions
     262             : 
     263             : template <>
     264             : struct Options::create_from_yaml<
     265             :     ScalarWave::BoundaryConditions::detail::
     266             :         ConstraintPreservingSphericalRadiationType> {
     267             :   template <typename Metavariables>
     268             :   static typename ScalarWave::BoundaryConditions::detail::
     269             :       ConstraintPreservingSphericalRadiationType
     270             :       create(const Options::Option& options) {
     271             :     return ScalarWave::BoundaryConditions::detail::
     272             :         convert_constraint_preserving_spherical_radiation_type_from_yaml(
     273             :             options);
     274             :   }
     275             : };