Line data    Source code

       1           1 : // Distributed under the MIT License.
       2             : // See LICENSE.txt for details.
       3             : 
       4             : /// \file
       5             : /// Defines MathFunctions::Gaussian.
       6             : 
       7             : #pragma once
       8             : 
       9             : #include <array>
      10             : #include <memory>
      11             : #include <pup.h>
      12             : 
      13             : #include "Options/String.hpp"
      14             : #include "PointwiseFunctions/MathFunctions/MathFunction.hpp"  // IWYU pragma: keep
      15             : #include "Utilities/Serialization/CharmPupable.hpp"
      16             : #include "Utilities/TMPL.hpp"
      17             : 
      18             : /// \cond
      19             : class DataVector;
      20             : /// \endcond
      21             : 
      22             : namespace MathFunctions {
      23             : template <size_t VolumeDim, typename Fr>
      24             : class Gaussian;
      25             : 
      26             : /*!
      27             :  *  \ingroup MathFunctionsGroup
      28             :  *  \brief 1D Gaussian \f$f = A \exp\left(-\frac{(x-x_0)^2}{w^2}\right)\f$
      29             :  *
      30             :  *  \details Input file options are: Amplitude, Width, and Center. The function
      31             :  *  takes input of type `double` or `DataVector` and returns
      32             :  *  the same type as the input type.
      33             :  */
      34             : template <typename Fr>
      35           1 : class Gaussian<1, Fr> : public MathFunction<1, Fr> {
      36             :  public:
      37           0 :   struct Amplitude {
      38           0 :     using type = double;
      39           0 :     static constexpr Options::String help = {"The amplitude."};
      40             :   };
      41             : 
      42           0 :   struct Width {
      43           0 :     using type = double;
      44           0 :     static constexpr Options::String help = {"The width."};
      45           0 :     static type lower_bound() { return 0.; }
      46             :   };
      47             : 
      48           0 :   struct Center {
      49           0 :     using type = double;
      50           0 :     static constexpr Options::String help = {"The center."};
      51             :   };
      52           0 :   using options = tmpl::list<Amplitude, Width, Center>;
      53             : 
      54           0 :   static constexpr Options::String help = {
      55             :       "Computes a Gaussian about an arbitrary coordinate center with given "
      56             :       "width and amplitude"};
      57             : 
      58           0 :   WRAPPED_PUPable_decl_base_template(SINGLE_ARG(MathFunction<1, Fr>),
      59             :                                      Gaussian);  // NOLINT
      60             : 
      61           0 :   explicit Gaussian(CkMigrateMessage* /*unused*/) {}
      62             : 
      63           0 :   Gaussian(double amplitude, double width, double center);
      64           0 :   Gaussian(double amplitude, double width, const std::array<double, 1>& center);
      65             : 
      66           0 :   Gaussian() = default;
      67           0 :   std::unique_ptr<MathFunction<1, Fr>> get_clone() const override;
      68             : 
      69           1 :   double operator()(const double& x) const override;
      70           0 :   DataVector operator()(const DataVector& x) const override;
      71             :   using MathFunction<1, Fr>::operator();
      72             : 
      73           1 :   double first_deriv(const double& x) const override;
      74           0 :   DataVector first_deriv(const DataVector& x) const override;
      75             :   using MathFunction<1, Fr>::first_deriv;
      76             : 
      77           1 :   double second_deriv(const double& x) const override;
      78           0 :   DataVector second_deriv(const DataVector& x) const override;
      79             :   using MathFunction<1, Fr>::second_deriv;
      80             : 
      81           1 :   double third_deriv(const double& x) const override;
      82           0 :   DataVector third_deriv(const DataVector& x) const override;
      83             :   using MathFunction<1, Fr>::third_deriv;
      84             : 
      85           0 :   bool operator==(const MathFunction<1, Fr>& other) const override;
      86           0 :   bool operator!=(const MathFunction<1, Fr>& other) const override;
      87             : 
      88           0 :   bool operator==(const Gaussian<1, Fr>& other) const;
      89           0 :   bool operator!=(const Gaussian<1, Fr>& other) const;
      90             :   // NOLINTNEXTLINE(google-runtime-references)
      91           0 :   void pup(PUP::er& p) override;
      92             : 
      93             :  private:
      94           0 :   double amplitude_{};
      95           0 :   double inverse_width_{};
      96           0 :   double center_{};
      97             : 
      98             :   template <typename T>
      99           0 :   T apply_call_operator(const T& x) const;
     100             :   template <typename T>
     101           0 :   T apply_first_deriv(const T& x) const;
     102             :   template <typename T>
     103           0 :   T apply_second_deriv(const T& x) const;
     104             :   template <typename T>
     105           0 :   T apply_third_deriv(const T& x) const;
     106             : };
     107             : 
     108             : /*!
     109             :  * \ingroup MathFunctionsGroup
     110             :  * \brief Gaussian \f$f = A \exp\left(-\frac{(x-x_0)^2}{w^2}\right)\f$
     111             :  *
     112             :  * \details Input file options are: Amplitude, Width, and Center. The function
     113             :  * takes input coordinates of type `tnsr::I<T, VolumeDim, Fr>`, where `T` is
     114             :  * e.g. `double` or `DataVector`, `Fr` is a frame (e.g. `Frame::Inertial`), and
     115             :  * `VolumeDim` is the dimension of the spatial volume, i.e. 2 or 3. (The case of
     116             :  * VolumeDim == 1 is handled specially by Gaussian<1, T, Frame::Inertial>.)
     117             :  */
     118             : template <size_t VolumeDim, typename Fr>
     119           1 : class Gaussian : public MathFunction<VolumeDim, Fr> {
     120             :  public:
     121           0 :   struct Amplitude {
     122           0 :     using type = double;
     123           0 :     static constexpr Options::String help = {"The amplitude."};
     124             :   };
     125             : 
     126           0 :   struct Width {
     127           0 :     using type = double;
     128           0 :     static constexpr Options::String help = {"The width."};
     129           0 :     static type lower_bound() { return 0.; }
     130             :   };
     131             : 
     132           0 :   struct Center {
     133           0 :     using type = std::array<double, VolumeDim>;
     134           0 :     static constexpr Options::String help = {"The center."};
     135             :   };
     136           0 :   using options = tmpl::list<Amplitude, Width, Center>;
     137             : 
     138           0 :   static constexpr Options::String help = {
     139             :       "Computes a Gaussian about an arbitrary coordinate center with given "
     140             :       "width and amplitude"};
     141             : 
     142           0 :   WRAPPED_PUPable_decl_base_template(SINGLE_ARG(MathFunction<VolumeDim, Fr>),
     143             :                                      Gaussian);  // NOLINT
     144             : 
     145           0 :   explicit Gaussian(CkMigrateMessage* /*unused*/) {}
     146             : 
     147           0 :   Gaussian(double amplitude, double width,
     148             :            const std::array<double, VolumeDim>& center);
     149           0 :   Gaussian() = default;
     150           0 :   std::unique_ptr<MathFunction<VolumeDim, Fr>> get_clone() const override;
     151             : 
     152           1 :   Scalar<double> operator()(
     153             :       const tnsr::I<double, VolumeDim, Fr>& x) const override;
     154           1 :   Scalar<DataVector> operator()(
     155             :       const tnsr::I<DataVector, VolumeDim, Fr>& x) const override;
     156             : 
     157           1 :   tnsr::i<double, VolumeDim, Fr> first_deriv(
     158             :       const tnsr::I<double, VolumeDim, Fr>& x) const override;
     159           1 :   tnsr::i<DataVector, VolumeDim, Fr> first_deriv(
     160             :       const tnsr::I<DataVector, VolumeDim, Fr>& x) const override;
     161             : 
     162           1 :   tnsr::ii<double, VolumeDim, Fr> second_deriv(
     163             :       const tnsr::I<double, VolumeDim, Fr>& x) const override;
     164           1 :   tnsr::ii<DataVector, VolumeDim, Fr> second_deriv(
     165             :       const tnsr::I<DataVector, VolumeDim, Fr>& x) const override;
     166             : 
     167           1 :   tnsr::iii<double, VolumeDim, Fr> third_deriv(
     168             :       const tnsr::I<double, VolumeDim, Fr>& x) const override;
     169           1 :   tnsr::iii<DataVector, VolumeDim, Fr> third_deriv(
     170             :       const tnsr::I<DataVector, VolumeDim, Fr>& x) const override;
     171             : 
     172           0 :   bool operator==(const MathFunction<VolumeDim, Fr>& other) const override;
     173           0 :   bool operator!=(const MathFunction<VolumeDim, Fr>& other) const override;
     174             : 
     175           0 :   bool operator==(const Gaussian<VolumeDim, Fr>& other) const;
     176           0 :   bool operator!=(const Gaussian<VolumeDim, Fr>& other) const;
     177             :   // NOLINTNEXTLINE(google-runtime-references)
     178           0 :   void pup(PUP::er& p) override;
     179             : 
     180             :  private:
     181           0 :   double amplitude_{};
     182           0 :   double inverse_width_{};
     183           0 :   std::array<double, VolumeDim> center_{};
     184             : 
     185             :   template <typename T>
     186           0 :   tnsr::I<T, VolumeDim, Fr> centered_coordinates(
     187             :       const tnsr::I<T, VolumeDim, Fr>& x) const;
     188             :   template <typename T>
     189           0 :   Scalar<T> apply_call_operator(
     190             :       const tnsr::I<T, VolumeDim, Fr>& centered_coords) const;
     191             :   template <typename T>
     192           0 :   tnsr::i<T, VolumeDim, Fr> apply_first_deriv(
     193             :       const tnsr::I<T, VolumeDim, Fr>& centered_coords,
     194             :       const Scalar<T>& gaussian) const;
     195             :   template <typename T>
     196           0 :   tnsr::ii<T, VolumeDim, Fr> apply_second_deriv(
     197             :       const tnsr::I<T, VolumeDim, Fr>& centered_coords,
     198             :       const Scalar<T>& gaussian,
     199             :       const tnsr::i<T, VolumeDim, Fr>& d_gaussian) const;
     200             :   template <typename T>
     201           0 :   tnsr::iii<T, VolumeDim, Fr> apply_third_deriv(
     202             :       const tnsr::I<T, VolumeDim, Fr>& centered_coords,
     203             :       const Scalar<T>& gaussian, const tnsr::i<T, VolumeDim, Fr>& d_gaussian,
     204             :       const tnsr::ii<T, VolumeDim, Fr>& d2_gaussian) const;
     205             : };
     206             : 
     207             : template <size_t VolumeDim, typename Fr>
     208           0 : bool operator!=(const Gaussian<VolumeDim, Fr>& lhs,
     209             :                 const Gaussian<VolumeDim, Fr>& rhs) {
     210             :   return not(lhs == rhs);
     211             : }
     212             : }  // namespace MathFunctions
     213             : 
     214             : /// \cond
     215             : template <size_t VolumeDim, typename Fr>
     216             : PUP::able::PUP_ID MathFunctions::Gaussian<VolumeDim, Fr>::my_PUP_ID =
     217             :     0;  // NOLINT
     218             : 
     219             : template <typename Fr>
     220             : PUP::able::PUP_ID MathFunctions::Gaussian<1, Fr>::my_PUP_ID = 0;  // NOLINT
     221             : /// \endcond