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 <memory>
       8             : #include <pup.h>
       9             : #include <string>
      10             : #include <vector>
      11             : 
      12             : #include "DataStructures/Tensor/TypeAliases.hpp"
      13             : #include "Domain/Tags.hpp"
      14             : #include "Domain/Tags/FaceNormal.hpp"
      15             : #include "Elliptic/BoundaryConditions/BoundaryCondition.hpp"
      16             : #include "Elliptic/BoundaryConditions/BoundaryConditionType.hpp"
      17             : #include "Elliptic/Systems/Xcts/FluxesAndSources.hpp"
      18             : #include "Options/String.hpp"
      19             : #include "Utilities/Gsl.hpp"
      20             : #include "Utilities/Serialization/CharmPupable.hpp"
      21             : #include "Utilities/TMPL.hpp"
      22             : 
      23             : /// \cond
      24             : class DataVector;
      25             : /// \endcond
      26             : 
      27             : namespace Xcts::BoundaryConditions {
      28             : 
      29             : /*!
      30             :  * \brief Impose Robin boundary conditions at the outer boundary
      31             :  *
      32             :  * Impose $\partial_r(r\psi)=0$, $\partial_r(\alpha\psi)=0$, and
      33             :  * $\partial_r(r\beta_\mathrm{excess}^j)=0$ on the boundary. These Robin
      34             :  * boundary conditions incur an error of order $1/R^2$, where $R$ is the outer
      35             :  * radius of the domain. This allows to place the outer boundary much closer
      36             :  * to the center than for Dirichlet boundary conditions
      37             :  * (Xcts::BoundaryConditions::Flatness), which incur an error of order $1/R$.
      38             :  *
      39             :  * The Robin boundary conditions are imposed as Neumann-type as follows:
      40             :  *
      41             :  * \begin{align*}
      42             :  * (n^i \partial_i \psi)_N &= -\frac{\psi}{r}, \\
      43             :  * (n^i \partial_i (\alpha\psi))_N &= -\frac{\alpha\psi}{r}, \\
      44             :  * (n^i (L\beta)^{ij})_N = n^i (L\beta)^{ij} - \left[
      45             :  *   \frac{\beta_\mathrm{excess}^j}{r}
      46             :  *   + \frac{1}{3} n^j \frac{n^i \beta_\mathrm{excess}^i}{r}
      47             :  *   + n^i \partial_i \beta_\mathrm{excess}^j
      48             :  *   + \frac{1}{3} n^j n^i n_k \partial_i \beta_\mathrm{excess}^k
      49             :  * \right]
      50             :  * \end{align*}
      51             :  *
      52             :  * Here, the condition on the longitudinal shift is derived by imposing the
      53             :  * Robin boundary condition
      54             :  * $n^i\partial_i\beta_\mathrm{excess}^j=-\beta_\mathrm{excess}^j/r$ only on the
      55             :  * normal component of the shift gradient. To do this we can use the projection
      56             :  * operator $P_{ij}=\delta_{ij}-n_i n_j$ to set the shift gradient to
      57             :  * $\partial_i\beta_\mathrm{excess}^j=P_{ik}\partial_k\beta_\mathrm{excess}^j
      58             :  * -n_i \beta_\mathrm{excess}^j/r$ and then apply the longitudinal operator.
      59             :  *
      60             :  * \tparam EnabledEquations The subset of XCTS equations that are being solved
      61             :  */
      62             : template <Xcts::Equations EnabledEquations>
      63           1 : class Robin : public elliptic::BoundaryConditions::BoundaryCondition<3> {
      64             :  private:
      65           0 :   using Base = elliptic::BoundaryConditions::BoundaryCondition<3>;
      66             : 
      67             :  public:
      68           0 :   using options = tmpl::list<>;
      69           0 :   static constexpr Options::String help =
      70             :       "Impose Robin boundary conditions at the outer boundary. "
      71             :       "They incur an error of order 1/R^2, where R is the outer radius.";
      72             : 
      73           0 :   Robin() = default;
      74           0 :   Robin(const Robin&) = default;
      75           0 :   Robin& operator=(const Robin&) = default;
      76           0 :   Robin(Robin&&) = default;
      77           0 :   Robin& operator=(Robin&&) = default;
      78           0 :   ~Robin() override = default;
      79             : 
      80             :   /// \cond
      81             :   explicit Robin(CkMigrateMessage* m) : Base(m) {}
      82             :   using PUP::able::register_constructor;
      83             :   WRAPPED_PUPable_decl_template(Robin);
      84             :   /// \endcond
      85             : 
      86           0 :   std::unique_ptr<domain::BoundaryConditions::BoundaryCondition> get_clone()
      87             :       const override;
      88             : 
      89           0 :   std::vector<elliptic::BoundaryConditionType> boundary_condition_types()
      90             :       const override;
      91             : 
      92           0 :   using argument_tags =
      93             :       tmpl::list<domain::Tags::Coordinates<3, Frame::Inertial>,
      94             :                  domain::Tags::FaceNormal<3, Frame::Inertial>>;
      95           0 :   using volume_tags = tmpl::list<>;
      96             : 
      97           0 :   void apply(gsl::not_null<Scalar<DataVector>*> conformal_factor_minus_one,
      98             :              gsl::not_null<Scalar<DataVector>*> n_dot_conformal_factor_gradient,
      99             :              const tnsr::i<DataVector, 3>& deriv_conformal_factor,
     100             :              const tnsr::I<DataVector, 3>& x,
     101             :              const tnsr::i<DataVector, 3>& face_normal) const;
     102             : 
     103           0 :   void apply(
     104             :       gsl::not_null<Scalar<DataVector>*> conformal_factor_minus_one,
     105             :       gsl::not_null<Scalar<DataVector>*> lapse_times_conformal_factor_minus_one,
     106             :       gsl::not_null<Scalar<DataVector>*> n_dot_conformal_factor_gradient,
     107             :       gsl::not_null<Scalar<DataVector>*>
     108             :           n_dot_lapse_times_conformal_factor_gradient,
     109             :       const tnsr::i<DataVector, 3>& deriv_conformal_factor,
     110             :       const tnsr::i<DataVector, 3>& deriv_lapse_times_conformal_factor,
     111             :       const tnsr::I<DataVector, 3>& x,
     112             :       const tnsr::i<DataVector, 3>& face_normal) const;
     113             : 
     114           0 :   void apply(
     115             :       gsl::not_null<Scalar<DataVector>*> conformal_factor_minus_one,
     116             :       gsl::not_null<Scalar<DataVector>*> lapse_times_conformal_factor_minus_one,
     117             :       gsl::not_null<tnsr::I<DataVector, 3>*> shift_excess,
     118             :       gsl::not_null<Scalar<DataVector>*> n_dot_conformal_factor_gradient,
     119             :       gsl::not_null<Scalar<DataVector>*>
     120             :           n_dot_lapse_times_conformal_factor_gradient,
     121             :       gsl::not_null<tnsr::I<DataVector, 3>*> n_dot_longitudinal_shift_excess,
     122             :       const tnsr::i<DataVector, 3>& deriv_conformal_factor,
     123             :       const tnsr::i<DataVector, 3>& deriv_lapse_times_conformal_factor,
     124             :       const tnsr::iJ<DataVector, 3>& deriv_shift_excess,
     125             :       const tnsr::I<DataVector, 3>& x,
     126             :       const tnsr::i<DataVector, 3>& face_normal) const;
     127             : 
     128           0 :   using argument_tags_linearized =
     129             :       tmpl::list<domain::Tags::Coordinates<3, Frame::Inertial>,
     130             :                  domain::Tags::FaceNormal<3, Frame::Inertial>>;
     131           0 :   using volume_tags_linearized = tmpl::list<>;
     132             : 
     133           0 :   void apply_linearized(
     134             :       gsl::not_null<Scalar<DataVector>*> conformal_factor_correction,
     135             :       gsl::not_null<Scalar<DataVector>*>
     136             :           n_dot_conformal_factor_gradient_correction,
     137             :       const tnsr::i<DataVector, 3>& deriv_conformal_factor_correction,
     138             :       const tnsr::I<DataVector, 3>& x,
     139             :       const tnsr::i<DataVector, 3>& face_normal) const;
     140             : 
     141           0 :   void apply_linearized(
     142             :       gsl::not_null<Scalar<DataVector>*> conformal_factor_correction,
     143             :       gsl::not_null<Scalar<DataVector>*>
     144             :           lapse_times_conformal_factor_correction,
     145             :       gsl::not_null<Scalar<DataVector>*>
     146             :           n_dot_conformal_factor_gradient_correction,
     147             :       gsl::not_null<Scalar<DataVector>*>
     148             :           n_dot_lapse_times_conformal_factor_gradient_correction,
     149             :       const tnsr::i<DataVector, 3>& deriv_conformal_factor_correction,
     150             :       const tnsr::i<DataVector, 3>&
     151             :           deriv_lapse_times_conformal_factor_correction,
     152             :       const tnsr::I<DataVector, 3>& x,
     153             :       const tnsr::i<DataVector, 3>& face_normal) const;
     154             : 
     155           0 :   void apply_linearized(
     156             :       gsl::not_null<Scalar<DataVector>*> conformal_factor_correction,
     157             :       gsl::not_null<Scalar<DataVector>*>
     158             :           lapse_times_conformal_factor_correction,
     159             :       gsl::not_null<tnsr::I<DataVector, 3>*> shift_excess_correction,
     160             :       gsl::not_null<Scalar<DataVector>*>
     161             :           n_dot_conformal_factor_gradient_correction,
     162             :       gsl::not_null<Scalar<DataVector>*>
     163             :           n_dot_lapse_times_conformal_factor_gradient_correction,
     164             :       gsl::not_null<tnsr::I<DataVector, 3>*>
     165             :           n_dot_longitudinal_shift_excess_correction,
     166             :       const tnsr::i<DataVector, 3>& deriv_conformal_factor_correction,
     167             :       const tnsr::i<DataVector, 3>&
     168             :           deriv_lapse_times_conformal_factor_correction,
     169             :       const tnsr::iJ<DataVector, 3>& deriv_shift_excess_correction,
     170             :       const tnsr::I<DataVector, 3>& x,
     171             :       const tnsr::i<DataVector, 3>& face_normal) const;
     172             : };
     173             : 
     174             : template <Xcts::Equations LocalEnabledEquations>
     175           0 : bool operator==(const Robin<LocalEnabledEquations>& lhs,
     176             :                 const Robin<LocalEnabledEquations>& rhs);
     177             : 
     178             : template <Xcts::Equations EnabledEquations>
     179           0 : bool operator!=(const Robin<EnabledEquations>& lhs,
     180             :                 const Robin<EnabledEquations>& rhs);
     181             : 
     182             : }  // namespace Xcts::BoundaryConditions