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 <pup.h>
       8             : 
       9             : #include "DataStructures/CachedTempBuffer.hpp"
      10             : #include "DataStructures/ComplexDataVector.hpp"
      11             : #include "DataStructures/DataBox/Prefixes.hpp"
      12             : #include "DataStructures/DataVector.hpp"
      13             : #include "DataStructures/Tensor/Tensor.hpp"
      14             : #include "Elliptic/Systems/Poisson/Tags.hpp"
      15             : #include "NumericalAlgorithms/LinearOperators/PartialDerivatives.hpp"
      16             : #include "Options/String.hpp"
      17             : #include "PointwiseFunctions/InitialDataUtilities/AnalyticSolution.hpp"
      18             : #include "Utilities/ContainerHelpers.hpp"
      19             : #include "Utilities/Gsl.hpp"
      20             : #include "Utilities/TMPL.hpp"
      21             : #include "Utilities/TaggedTuple.hpp"
      22             : 
      23             : namespace Poisson::Solutions {
      24             : 
      25             : namespace detail {
      26             : template <typename DataType, size_t Dim>
      27             : struct LorentzianVariables {
      28             :   using Cache = CachedTempBuffer<
      29             :       Tags::Field<DataType>,
      30             :       ::Tags::deriv<Tags::Field<DataType>, tmpl::size_t<Dim>, Frame::Inertial>,
      31             :       ::Tags::Flux<Tags::Field<DataType>, tmpl::size_t<Dim>, Frame::Inertial>,
      32             :       ::Tags::FixedSource<Tags::Field<DataType>>>;
      33             : 
      34             :   // NOLINTNEXTLINE(cppcoreguidelines-avoid-const-or-ref-data-members)
      35             :   const tnsr::I<DataVector, Dim>& x;
      36             :   double constant;
      37             :   double complex_phase;
      38             : 
      39             :   void operator()(gsl::not_null<Scalar<DataType>*> field,
      40             :                   gsl::not_null<Cache*> cache,
      41             :                   Tags::Field<DataType> /*meta*/) const;
      42             :   void operator()(gsl::not_null<tnsr::i<DataType, Dim>*> field_gradient,
      43             :                   gsl::not_null<Cache*> cache,
      44             :                   ::Tags::deriv<Tags::Field<DataType>, tmpl::size_t<Dim>,
      45             :                                 Frame::Inertial> /*meta*/) const;
      46             :   void operator()(gsl::not_null<tnsr::I<DataType, Dim>*> flux_for_field,
      47             :                   gsl::not_null<Cache*> cache,
      48             :                   ::Tags::Flux<Tags::Field<DataType>, tmpl::size_t<Dim>,
      49             :                                Frame::Inertial> /*meta*/) const;
      50             :   void operator()(gsl::not_null<Scalar<DataType>*> fixed_source_for_field,
      51             :                   gsl::not_null<Cache*> cache,
      52             :                   ::Tags::FixedSource<Tags::Field<DataType>> /*meta*/) const;
      53             : };
      54             : }  // namespace detail
      55             : 
      56             : /*!
      57             :  * \brief A Lorentzian solution to the Poisson equation
      58             :  *
      59             :  * \details This implements the Lorentzian solution
      60             :  * \f$u(\boldsymbol{x})=\left(1+r^2\right)^{-\frac{1}{2}}\f$ to the
      61             :  * three-dimensional Poisson equation
      62             :  * \f$-\Delta u(\boldsymbol{x})=f(\boldsymbol{x})\f$, where
      63             :  * \f$r^2=x^2+y^2+z^2\f$. The corresponding source is
      64             :  * \f$f(\boldsymbol{x})=3\left(1+r^2\right)^{-\frac{5}{2}}\f$.
      65             :  *
      66             :  * If `DataType` is `ComplexDataVector`, the solution is multiplied by
      67             :  * `exp(i * complex_phase)` to rotate it in the complex plane. This allows to
      68             :  * use this solution for the complex Poisson equation.
      69             :  *
      70             :  * \note Corresponding 1D and 2D solutions are not implemented yet.
      71             :  */
      72             : template <size_t Dim, typename DataType = DataVector>
      73           1 : class Lorentzian : public elliptic::analytic_data::AnalyticSolution {
      74             :   static_assert(
      75             :       Dim == 3,
      76             :       "This solution is currently implemented in 3 spatial dimensions only");
      77             : 
      78             :  public:
      79           0 :   struct PlusConstant {
      80           0 :     using type = double;
      81           0 :     static constexpr Options::String help{"Constant added to the solution."};
      82             :   };
      83             : 
      84           0 :   struct ComplexPhase {
      85           0 :     using type = double;
      86           0 :     static constexpr Options::String help{
      87             :         "Phase 'phi' of a complex exponential 'exp(i phi)' that rotates the "
      88             :         "solution in the complex plane."};
      89             :   };
      90             : 
      91           0 :   using options = tmpl::flatten<tmpl::list<
      92             :       PlusConstant,
      93             :       tmpl::conditional_t<std::is_same_v<DataType, ComplexDataVector>,
      94             :                           ComplexPhase, tmpl::list<>>>>;
      95           0 :   static constexpr Options::String help{
      96             :       "A Lorentzian solution to the Poisson equation."};
      97             : 
      98           0 :   Lorentzian() = default;
      99           0 :   Lorentzian(const Lorentzian&) = default;
     100           0 :   Lorentzian& operator=(const Lorentzian&) = default;
     101           0 :   Lorentzian(Lorentzian&&) = default;
     102           0 :   Lorentzian& operator=(Lorentzian&&) = default;
     103           0 :   ~Lorentzian() override = default;
     104             : 
     105           0 :   explicit Lorentzian(const double constant, const double complex_phase = 0.)
     106             :       : constant_(constant), complex_phase_(complex_phase) {
     107             :     ASSERT((std::is_same_v<DataType, ComplexDataVector> or complex_phase == 0.),
     108             :            "The complex phase is only supported for ComplexDataVector.");
     109             :   }
     110             : 
     111           0 :   double constant() const { return constant_; }
     112           0 :   double complex_phase() const { return complex_phase_; }
     113             : 
     114           0 :   std::unique_ptr<elliptic::analytic_data::AnalyticSolution> get_clone()
     115             :       const override {
     116             :     return std::make_unique<Lorentzian>(*this);
     117             :   }
     118             : 
     119             :   /// \cond
     120             :   explicit Lorentzian(CkMigrateMessage* m)
     121             :       : elliptic::analytic_data::AnalyticSolution(m) {}
     122             :   using PUP::able::register_constructor;
     123             :   WRAPPED_PUPable_decl_template(Lorentzian);  // NOLINT
     124             :   /// \endcond
     125             : 
     126             :   template <typename... RequestedTags>
     127           0 :   tuples::TaggedTuple<RequestedTags...> variables(
     128             :       const tnsr::I<DataVector, Dim>& x,
     129             :       tmpl::list<RequestedTags...> /*meta*/) const {
     130             :     using VarsComputer = detail::LorentzianVariables<DataType, Dim>;
     131             :     typename VarsComputer::Cache cache{get_size(*x.begin())};
     132             :     const VarsComputer computer{x, constant_, complex_phase_};
     133             :     return {cache.get_var(computer, RequestedTags{})...};
     134             :   }
     135             : 
     136           0 :   void pup(PUP::er& p) override {
     137             :     elliptic::analytic_data::AnalyticSolution::pup(p);
     138             :     p | constant_;
     139             :     p | complex_phase_;
     140             :   }
     141             : 
     142             :  private:
     143           0 :   double constant_ = std::numeric_limits<double>::signaling_NaN();
     144           0 :   double complex_phase_ = std::numeric_limits<double>::signaling_NaN();
     145             : };
     146             : 
     147             : /// \cond
     148             : template <size_t Dim, typename DataType>
     149             : PUP::able::PUP_ID Lorentzian<Dim, DataType>::my_PUP_ID = 0;  // NOLINT
     150             : /// \endcond
     151             : 
     152             : template <size_t Dim, typename DataType>
     153           0 : bool operator==(const Lorentzian<Dim, DataType>& lhs,
     154             :                 const Lorentzian<Dim, DataType>& rhs) {
     155             :   return lhs.constant() == rhs.constant() and
     156             :          lhs.complex_phase() == rhs.complex_phase();
     157             : }
     158             : 
     159             : template <size_t Dim, typename DataType>
     160           0 : bool operator!=(const Lorentzian<Dim, DataType>& lhs,
     161             :                 const Lorentzian<Dim, DataType>& rhs) {
     162             :   return not(lhs == rhs);
     163             : }
     164             : 
     165             : }  // namespace Poisson::Solutions