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