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