Sinusoid.hpp
1 // 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
8 #include "Evolution/Systems/Burgers/Tags.hpp" // IWYU pragma: keep
9 #include "Options/Options.hpp"
10 #include "Utilities/TMPL.hpp"
11 #include "Utilities/TaggedTuple.hpp"
12 
13 /// \cond
14 class DataVector;
15 // IWYU pragma: no_forward_declare Tensor
16 namespace PUP {
17 class er;
18 } // namespace PUP
19 /// \endcond
20 
21 namespace Burgers {
22 namespace AnalyticData {
23 /*!
24  * \brief Analytic data (with an "exact" solution known) that is periodic over
25  * the interval \f$[0,2\pi]\f$.
26  *
27  * The initial data is given by:
28  *
29  * \f{align}{
30  * u(x, 0) = \sin(x)
31  * \f}
32  *
33  * At future times the analytic solution can be found by solving the
34  * transcendental equation \cite Harten19973
35  *
36  * \f{align}{
37  * \label{eq:transcendental burgers periodic}
38  * \mathcal{F}=\sin\left(x-\mathcal{F}t\right)
39  * \f}
40  *
41  * on the interval \f$x\in(0,\pi)\f$. The solution from \f$x\in(\pi,2\pi)\f$ is
42  * given by \f$\mathcal{F}(x, t)=-\mathcal{F}(2\pi-x,t)\f$. The transcendental
43  * equation \f$(\ref{eq:transcendental burgers periodic})\f$ can be solved with
44  * a Newton-Raphson iterative scheme. Since this can be quite sensitive to the
45  * initial guess we implement this solution as analytic data. The python code
46  * below can be used to compute the analytic solution if desired.
47  *
48  * At time \f$1\f$ the solution develops a discontinuity at \f$x=\pi\f$ followed
49  * by the amplitude of the solution decaying over time.
50  *
51  * \note We have rescaled \f$x\f$ and \f$t\f$ by \f$\pi\f$ compared to
52  * \cite Harten19973.
53  *
54  * \code{py}
55  import numpy as np
56  from scipy.optimize import newton
57 
58  # x_grid is a np.array of positions at which to evaluate the solution
59  def burgers_periodic(x_grid, time):
60  def opt_fun(F, x, t):
61  return np.sin((x - F * t)) - F
62 
63  results = []
64  for i in range(len(x_grid)):
65  x = x_grid[i]
66  greater_than_pi = False
67  if x > np.pi:
68  x = x - np.pi
69  x = -x
70  x = x + np.pi
71  greater_than_pi = True
72 
73  guess = 0.0
74  if len(results) > 0:
75  if results[-1] < 0.0:
76  guess = -results[-1]
77  else:
78  guess = results[-1]
79  res = newton(lambda F: opt_fun(F, x, time), x0=guess)
80 
81  if greater_than_pi:
82  results.append(-res)
83  else:
84  results.append(res)
85 
86  return np.asarray(results)
87  * \endcode
88  */
89 class Sinusoid {
90  public:
91  using options = tmpl::list<>;
92  static constexpr OptionString help{
93  "A solution that is periodic over the interval [0,2pi]. The solution "
94  "starts as a sinusoid: u(x,0) = sin(x) and develops a "
95  "discontinuity at x=pi and t=1."};
96 
97  Sinusoid() = default;
98  Sinusoid(const Sinusoid&) noexcept = default;
99  Sinusoid& operator=(const Sinusoid&) noexcept = default;
100  Sinusoid(Sinusoid&&) noexcept = default;
101  Sinusoid& operator=(Sinusoid&&) noexcept = default;
102  ~Sinusoid() noexcept = default;
103 
104  template <typename T>
105  Scalar<T> u(const tnsr::I<T, 1>& x) const noexcept;
106 
107  tuples::TaggedTuple<Tags::U> variables(const tnsr::I<DataVector, 1>& x,
108  tmpl::list<Tags::U> /*meta*/) const
109  noexcept;
110 
111  // clang-tidy: no pass by reference
112  void pup(PUP::er& p) noexcept; // NOLINT
113 };
114 
115 bool operator==(const Sinusoid& /*lhs*/, const Sinusoid& /*rhs*/) noexcept;
116 
117 bool operator!=(const Sinusoid& lhs, const Sinusoid& rhs) noexcept;
118 } // namespace AnalyticData
119 } // namespace Burgers
Options.hpp
Burgers
Items related to evolving the Burgers equation .
Definition: Characteristics.hpp:25
Burgers::AnalyticData::Sinusoid
Analytic data (with an "exact" solution known) that is periodic over the interval .
Definition: Sinusoid.hpp:89
tuples::TaggedTuple
An associative container that is indexed by structs.
Definition: TaggedTuple.hpp:272
DataVector
Stores a collection of function values.
Definition: DataVector.hpp:42
Scalar
Tensor< T, Symmetry<>, index_list<> > Scalar
Definition: TypeAliases.hpp:21
TypeAliases.hpp
Prefixes.hpp
OptionString
const char *const OptionString
The string used in option structs.
Definition: Options.hpp:30
TMPL.hpp