Line data Source code
1 0 : // Distributed under the MIT License.
2 : // See LICENSE.txt for details.
3 :
4 : #pragma once
5 :
6 : #include <array>
7 : #include <cstddef>
8 : #include <limits>
9 :
10 : #include "DataStructures/DataVector.hpp"
11 : #include "DataStructures/TaggedTuple.hpp"
12 : #include "DataStructures/Tensor/Tensor.hpp"
13 : #include "DataStructures/Tensor/TypeAliases.hpp"
14 : #include "Options/Context.hpp"
15 : #include "Options/String.hpp"
16 : #include "PointwiseFunctions/AnalyticData/AnalyticData.hpp"
17 : #include "PointwiseFunctions/Hydro/EquationsOfState/EquationOfState.hpp"
18 : #include "PointwiseFunctions/Hydro/EquationsOfState/IdealFluid.hpp"
19 : #include "PointwiseFunctions/Hydro/Tags.hpp"
20 : #include "PointwiseFunctions/InitialDataUtilities/InitialData.hpp"
21 : #include "Utilities/MakeArray.hpp"
22 : #include "Utilities/TMPL.hpp"
23 :
24 : /// \cond
25 : namespace PUP {
26 : class er;
27 : } // namespace PUP
28 : /// \endcond
29 :
30 : namespace NewtonianEuler::AnalyticData {
31 : /*!
32 : * \brief A cylindrical or spherical Sod explosion \cite Toro2009 \cite Sod19781
33 : *
34 : * Common initial conditions are:
35 : *
36 : * \f{align*}{
37 : * (\rho, v^i, p) =
38 : * &\left\{
39 : * \begin{array}{ll}
40 : * (1 ,0, 1) & \mathrm{if} \; r \le 0.5 \\
41 : * (0.125 ,0, 0.1) & \mathrm{if} \; r > 0.5
42 : * \end{array}\right.
43 : * \f}
44 : *
45 : * where \f$r\f$ is the cylindrical (2d) or spherical (3d) radius. This test
46 : * problem uses an adiabatic index of 1.4. A reference solution can be computed
47 : * in 1d by solving the Newtonian Euler equations in cylindrical or spherical
48 : * symmetry. Note that the inner and outer density and pressure, as well as
49 : * where the initial discontinuity is can be chosen arbitrarily.
50 : */
51 : template <size_t Dim>
52 1 : class SodExplosion : public evolution::initial_data::InitialData,
53 : public MarkAsAnalyticData {
54 : public:
55 : static_assert(Dim > 1, "Sod explosion is a 2d and 3d problem.");
56 0 : using equation_of_state_type = EquationsOfState::IdealFluid<false>;
57 :
58 : /// Initial radius of the discontinuity
59 1 : struct InitialRadius {
60 0 : using type = double;
61 0 : static constexpr Options::String help = {
62 : "The initial radius of the discontinuity."};
63 0 : static type lower_bound() { return 0.0; }
64 : };
65 :
66 0 : struct InnerMassDensity {
67 0 : using type = double;
68 0 : static constexpr Options::String help = {"The inner mass density."};
69 0 : static type lower_bound() { return 0.0; }
70 : };
71 :
72 0 : struct InnerPressure {
73 0 : using type = double;
74 0 : static constexpr Options::String help = {"The inner pressure."};
75 0 : static type lower_bound() { return 0.0; }
76 : };
77 :
78 0 : struct OuterMassDensity {
79 0 : using type = double;
80 0 : static constexpr Options::String help = {"The outer mass density."};
81 0 : static type lower_bound() { return 0.0; }
82 : };
83 :
84 0 : struct OuterPressure {
85 0 : using type = double;
86 0 : static constexpr Options::String help = {"The outer pressure."};
87 0 : static type lower_bound() { return 0.0; }
88 : };
89 :
90 0 : using options = tmpl::list<InitialRadius, InnerMassDensity, InnerPressure,
91 : OuterMassDensity, OuterPressure>;
92 :
93 0 : static constexpr Options::String help = {
94 : "Cylindrical or spherical Sod explosion."};
95 :
96 0 : SodExplosion() = default;
97 0 : SodExplosion(const SodExplosion& /*rhs*/) = default;
98 0 : SodExplosion& operator=(const SodExplosion& /*rhs*/) = default;
99 0 : SodExplosion(SodExplosion&& /*rhs*/) = default;
100 0 : SodExplosion& operator=(SodExplosion&& /*rhs*/) = default;
101 0 : ~SodExplosion() override = default;
102 :
103 0 : auto get_clone() const
104 : -> std::unique_ptr<evolution::initial_data::InitialData> override;
105 :
106 : /// \cond
107 : explicit SodExplosion(CkMigrateMessage* msg);
108 : using PUP::able::register_constructor;
109 : WRAPPED_PUPable_decl_template(SodExplosion);
110 : /// \endcond
111 :
112 0 : SodExplosion(double initial_radius, double inner_mass_density,
113 : double inner_pressure, double outer_mass_density,
114 : double outer_pressure, const Options::Context& context = {});
115 :
116 : /// Retrieve a collection of hydrodynamic variables at position x
117 : template <typename... Tags>
118 1 : tuples::TaggedTuple<Tags...> variables(
119 : const tnsr::I<DataVector, Dim, Frame::Inertial>& x,
120 : tmpl::list<Tags...> /*meta*/) const {
121 : return {tuples::get<Tags>(variables(x, tmpl::list<Tags>{}))...};
122 : }
123 :
124 0 : const equation_of_state_type& equation_of_state() const {
125 : return equation_of_state_;
126 : }
127 :
128 : // NOLINTNEXTLINE(google-runtime-references)
129 0 : void pup(PUP::er& /*p*/) override;
130 :
131 : private:
132 0 : tuples::TaggedTuple<hydro::Tags::RestMassDensity<DataVector>> variables(
133 : const tnsr::I<DataVector, Dim, Frame::Inertial>& x,
134 : tmpl::list<hydro::Tags::RestMassDensity<DataVector>> /*meta*/) const;
135 :
136 : tuples::TaggedTuple<
137 : hydro::Tags::SpatialVelocity<DataVector, Dim, Frame::Inertial>>
138 0 : variables(const tnsr::I<DataVector, Dim, Frame::Inertial>& x,
139 : tmpl::list<hydro::Tags::SpatialVelocity<
140 : DataVector, Dim, Frame::Inertial>> /*meta*/) const;
141 :
142 0 : tuples::TaggedTuple<hydro::Tags::Pressure<DataVector>> variables(
143 : const tnsr::I<DataVector, Dim, Frame::Inertial>& x,
144 : tmpl::list<hydro::Tags::Pressure<DataVector>> /*meta*/) const;
145 :
146 : tuples::TaggedTuple<hydro::Tags::SpecificInternalEnergy<DataVector>>
147 0 : variables(
148 : const tnsr::I<DataVector, Dim, Frame::Inertial>& x,
149 : tmpl::list<hydro::Tags::SpecificInternalEnergy<DataVector>> /*meta*/)
150 : const;
151 :
152 : template <size_t SpatialDim>
153 : friend bool
154 0 : operator==( // NOLINT (clang-tidy: readability-redundant-declaration)
155 : const SodExplosion<SpatialDim>& lhs, const SodExplosion<SpatialDim>& rhs);
156 :
157 0 : double initial_radius_ = std::numeric_limits<double>::signaling_NaN();
158 0 : double inner_mass_density_ = std::numeric_limits<double>::signaling_NaN();
159 0 : double inner_pressure_ = std::numeric_limits<double>::signaling_NaN();
160 0 : double outer_mass_density_ = std::numeric_limits<double>::signaling_NaN();
161 0 : double outer_pressure_ = std::numeric_limits<double>::signaling_NaN();
162 0 : equation_of_state_type equation_of_state_{};
163 : };
164 :
165 : template <size_t Dim>
166 0 : bool operator!=(const SodExplosion<Dim>& lhs, const SodExplosion<Dim>& rhs);
167 : } // namespace NewtonianEuler::AnalyticData
|
