Line data    Source code

       1            0 : // Distributed under the MIT License.
       2              : // See LICENSE.txt for details.
       3              : 
       4              : #pragma once
       5              : 
       6              : #include "DataStructures/DataBox/Prefixes.hpp"
       7              : #include "DataStructures/TaggedTuple.hpp"
       8              : #include "DataStructures/Tensor/TypeAliases.hpp"
       9              : #include "Evolution/Systems/Burgers/Tags.hpp"
      10              : #include "Options/String.hpp"
      11              : #include "PointwiseFunctions/AnalyticData/AnalyticData.hpp"
      12              : #include "PointwiseFunctions/InitialDataUtilities/InitialData.hpp"
      13              : #include "Utilities/Serialization/CharmPupable.hpp"
      14              : #include "Utilities/TMPL.hpp"
      15              : 
      16              : /// \cond
      17              : class DataVector;
      18              : 
      19              : namespace PUP {
      20              : class er;
      21              : }  // namespace PUP
      22              : /// \endcond
      23              : 
      24              : namespace Burgers::AnalyticData {
      25              : /*!
      26              :  * \brief Analytic data (with an "exact" solution known) that is periodic over
      27              :  * the interval \f$[0,2\pi]\f$.
      28              :  *
      29              :  * The initial data is given by:
      30              :  *
      31              :  * \f{align}{
      32              :  *   u(x, 0) = \sin(x)
      33              :  * \f}
      34              :  *
      35              :  * At future times the analytic solution can be found by solving the
      36              :  * transcendental equation \cite Harten19973
      37              :  *
      38              :  * \f{align}{
      39              :  *   \label{eq:transcendental burgers periodic}
      40              :  *   \mathcal{F}=\sin\left(x-\mathcal{F}t\right)
      41              :  * \f}
      42              :  *
      43              :  * on the interval \f$x\in(0,\pi)\f$. The solution from \f$x\in(\pi,2\pi)\f$ is
      44              :  * given by \f$\mathcal{F}(x, t)=-\mathcal{F}(2\pi-x,t)\f$. The transcendental
      45              :  * equation \f$(\ref{eq:transcendental burgers periodic})\f$ can be solved with
      46              :  * a Newton-Raphson iterative scheme. Since this can be quite sensitive to the
      47              :  * initial guess we implement this solution as analytic data. The python code
      48              :  * below can be used to compute the analytic solution if desired.
      49              :  *
      50              :  * At time \f$1\f$ the solution develops a discontinuity at \f$x=\pi\f$ followed
      51              :  * by the amplitude of the solution decaying over time.
      52              :  *
      53              :  * \note We have rescaled \f$x\f$ and \f$t\f$ by \f$\pi\f$ compared to
      54              :  * \cite Harten19973.
      55              :  *
      56              :  * \code{py}
      57              :    import numpy as np
      58              :    from scipy.optimize import newton
      59              : 
      60              :    # x_grid is a np.array of positions at which to evaluate the solution
      61              :    def burgers_periodic(x_grid, time):
      62              :        def opt_fun(F, x, t):
      63              :            return np.sin((x - F * t)) - F
      64              : 
      65              :        results = []
      66              :        for i in range(len(x_grid)):
      67              :            x = x_grid[i]
      68              :            greater_than_pi = False
      69              :            if x > np.pi:
      70              :                x = x - np.pi
      71              :                x = -x
      72              :                x = x + np.pi
      73              :                greater_than_pi = True
      74              : 
      75              :            guess = 0.0
      76              :            if len(results) > 0:
      77              :                if results[-1] < 0.0:
      78              :                    guess = -results[-1]
      79              :                else:
      80              :                    guess = results[-1]
      81              :            res = newton(lambda F: opt_fun(F, x, time), x0=guess)
      82              : 
      83              :            if greater_than_pi:
      84              :                results.append(-res)
      85              :            else:
      86              :                results.append(res)
      87              : 
      88              :        return np.asarray(results)
      89              :  * \endcode
      90              :  */
      91            1 : class Sinusoid : public evolution::initial_data::InitialData,
      92              :                  public MarkAsAnalyticData {
      93              :  public:
      94            0 :   using options = tmpl::list<>;
      95            0 :   static constexpr Options::String help{
      96              :       "A solution that is periodic over the interval [0,2pi]. The solution "
      97              :       "starts as a sinusoid: u(x,0) = sin(x) and develops a "
      98              :       "discontinuity at x=pi and t=1."};
      99              : 
     100            0 :   Sinusoid() = default;
     101            0 :   Sinusoid(const Sinusoid&) = default;
     102            0 :   Sinusoid& operator=(const Sinusoid&) = default;
     103            0 :   Sinusoid(Sinusoid&&) = default;
     104            0 :   Sinusoid& operator=(Sinusoid&&) = default;
     105            0 :   ~Sinusoid() override = default;
     106              : 
     107            0 :   auto get_clone() const
     108              :       -> std::unique_ptr<evolution::initial_data::InitialData> override;
     109              : 
     110              :   /// \cond
     111              :   explicit Sinusoid(CkMigrateMessage* msg);
     112              :   using PUP::able::register_constructor;
     113              :   WRAPPED_PUPable_decl_template(Sinusoid);
     114              :   /// \endcond
     115              : 
     116              :   template <typename T>
     117            0 :   Scalar<T> u(const tnsr::I<T, 1>& x) const;
     118              : 
     119            0 :   tuples::TaggedTuple<Tags::U> variables(const tnsr::I<DataVector, 1>& x,
     120              :                                          tmpl::list<Tags::U> /*meta*/) const;
     121              : 
     122              :   // NOLINTNEXTLINE(google-runtime-references)
     123            0 :   void pup(PUP::er& p) override;
     124              : };
     125              : 
     126            0 : bool operator==(const Sinusoid& /*lhs*/, const Sinusoid& /*rhs*/);
     127              : 
     128            0 : bool operator!=(const Sinusoid& lhs, const Sinusoid& rhs);
     129              : }  // namespace Burgers::AnalyticData