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" // IWYU pragma: keep 7 : #include "DataStructures/Tensor/TypeAliases.hpp" 8 : #include "Evolution/Systems/Burgers/Tags.hpp" // IWYU pragma: keep 9 : #include "Options/String.hpp" 10 : #include "PointwiseFunctions/AnalyticData/AnalyticData.hpp" 11 : #include "PointwiseFunctions/InitialDataUtilities/InitialData.hpp" 12 : #include "Utilities/Serialization/CharmPupable.hpp" 13 : #include "Utilities/TMPL.hpp" 14 : #include "Utilities/TaggedTuple.hpp" 15 : 16 : /// \cond 17 : class DataVector; 18 : // IWYU pragma: no_forward_declare Tensor 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