Line data Source code
1 0 : // Distributed under the MIT License.
2 : // See LICENSE.txt for details.
3 :
4 : #pragma once
5 :
6 : #include <cstddef>
7 : #include <iosfwd>
8 : #include <memory>
9 : #include <optional>
10 : #include <string>
11 : #include <vector>
12 :
13 : #include "DataStructures/DataVector.hpp"
14 : #include "DataStructures/Tensor/TypeAliases.hpp"
15 : #include "Options/Auto.hpp"
16 : #include "Options/String.hpp"
17 : #include "PointwiseFunctions/InitialDataUtilities/InitialData.hpp"
18 : #include "Utilities/Gsl.hpp"
19 : #include "Utilities/Serialization/CharmPupable.hpp"
20 :
21 : /// \cond
22 : namespace Options {
23 : class Option;
24 : template <typename T>
25 : struct create_from_yaml;
26 : } // namespace Options
27 : namespace PUP {
28 : class er;
29 : } // namespace PUP
30 : /// \endcond
31 :
32 : namespace evolution::initial_data {
33 :
34 : /// How the noise amplitude passed to `add_noise_to_tensor` is interpreted.
35 1 : enum class NoiseAmplitudeType : uint8_t {
36 : /// The noise is drawn from \f$[-A,\,A]\f$.
37 : Absolute,
38 : /// The noise is drawn from \f$[-A M,\,A M]\f$, where \f$M\f$ is the
39 : /// maximum absolute value over all stored components and grid points of the
40 : /// tensor being perturbed.
41 : Relative
42 : };
43 :
44 0 : std::ostream& operator<<(std::ostream& os, NoiseAmplitudeType value);
45 :
46 : /*!
47 : * \brief Add uniform random noise to all stored components of \p tensor.
48 : *
49 : * Each independent stored component receives noise drawn independently from
50 : * \f$\text{Uniform}[-A,\,A]\f$. For `NoiseAmplitudeType::Absolute`,
51 : * \f$A\f$ is \p amplitude. For `NoiseAmplitudeType::Relative`, \f$A\f$ is
52 : * \p amplitude times the maximum absolute value over all stored components and
53 : * grid points of \p tensor, so a tensor that is identically zero receives no
54 : * noise. The RNG seed
55 : * for component \f$i\f$ is `hash(element_seed, component_offset + i)`, so
56 : * calling this function with different \p component_offset values keeps
57 : * separate tensor fields independent even when \p element_seed is the same.
58 : *
59 : * \param tensor Tensor to perturb in place.
60 : * \param amplitude Sets \f$\pm A\f$ of the noise interval (absolute or
61 : * relative, depending on \p amplitude_type).
62 : * \param amplitude_type Whether \p amplitude is absolute or relative to the
63 : * maximum absolute value of \p tensor.
64 : * \param element_seed Per-element seed, typically constructed by hashing a
65 : * base seed with the inertial coordinates of the element's first grid point.
66 : * \param component_offset Offset added to the component index before hashing,
67 : * used to keep different tensor fields independent.
68 : */
69 : template <typename TensorType>
70 1 : void add_noise_to_tensor(gsl::not_null<TensorType*> tensor, double amplitude,
71 : NoiseAmplitudeType amplitude_type, size_t element_seed,
72 : size_t component_offset);
73 :
74 : /*!
75 : * \brief Combine \p base_seed with the first grid point of \p inertial_coords
76 : * to produce a per-element RNG seed.
77 : *
78 : * Because most domains have elements have distinct first grid points, this
79 : * gives independent noise on almost every element even when `base_seed` is the
80 : * same.
81 : */
82 : template <size_t Dim>
83 1 : size_t make_element_seed(
84 : size_t base_seed,
85 : const tnsr::I<DataVector, Dim, Frame::Inertial>& inertial_coords);
86 :
87 : /*!
88 : * \brief Wraps any analytic `InitialData` and adds uniform random noise to
89 : * selected evolution variables after the inner solution has initialized them.
90 : *
91 : * \details The inner solution first sets all evolution variables in the normal
92 : * way. Then each variable whose `db::tag_name` appears in `Variables` (or
93 : * every variable if `All` is listed) has independent noise from
94 : * \f$\text{Uniform}[-A,\,A]\f$ added to every grid point and every independent
95 : * tensor component. If `AmplitudeType` is `Absolute`, \f$A\f$ is `Amplitude`.
96 : * If `AmplitudeType` is `Relative`, \f$A\f$ is `Amplitude` times the maximum
97 : * absolute value of the variable on the element, taken over all of its
98 : * independent components and grid points. Relative noise is therefore suited
99 : * to systems whose variables differ greatly in magnitude. Because the maximum
100 : * is computed element by element, \f$A\f$ can differ slightly between
101 : * elements, and a variable that is identically zero on an element (e.g.
102 : * `Pi` for Minkowski) receives no noise there; use `Absolute` to perturb such
103 : * variables.
104 : *
105 : * Valid variable names are the `db::tag_name`s of the system's evolution
106 : * variables. Examples:
107 : * - **ScalarWave**: `Psi`, `Pi`, `Phi`
108 : * - **GeneralizedHarmonic**: `SpacetimeMetric`, `Pi`, `Phi`
109 : *
110 : * Different elements receive independent noise because the RNG seed is mixed
111 : * with the inertial coordinates of each element's first grid point. Set
112 : * `Seed` to a fixed integer for reproducible results.
113 : *
114 : * \note Only analytic initial data (analytic solutions or analytic data) are
115 : * supported as the inner `Solution`.
116 : *
117 : * \note Nesting `WithNoise` inside `WithNoise` is not supported.
118 : *
119 : * \note Systems with primitive and conservative variables (e.g.
120 : * ValenciaDivClean, NewtonianEuler) are not yet supported and trigger a
121 : * `static_assert` in `evolution::Initialization::Actions::SetVariables`.
122 : *
123 : * \note **GeneralizedHarmonic systems**: the initialization phase
124 : * `InitializeInitialDataDependentQuantities` runs
125 : * `gh::gauges::SetPiAndPhiFromConstraints` *after* `WithNoise` applies noise.
126 : * That action unconditionally overwrites `Phi` with the numerical spatial
127 : * derivative of `SpacetimeMetric` and recomputes `Pi` from the gauge source
128 : * function and the current geometry. As a result, noise added to `Pi` or `Phi`
129 : * is effectively erased. Only noise added to `SpacetimeMetric` persists and
130 : * propagates into the other variables: it appears in `Phi` amplified by
131 : * roughly one over the grid spacing (due to differentiation), and in `Pi`
132 : * through the geometric quantities (lapse, shift) extracted from the metric.
133 : *
134 : * Example input-file snippet for ScalarWave:
135 : * \code{.yaml}
136 : * InitialData:
137 : * WithNoise:
138 : * Amplitude: 1.0e-6
139 : * AmplitudeType: Absolute
140 : * Seed: 42
141 : * Variables: [Psi, Pi]
142 : * Solution:
143 : * PlaneWave:
144 : * ...
145 : * \endcode
146 : */
147 1 : class WithNoise : public evolution::initial_data::InitialData {
148 : public:
149 : /// The analytic initial data to wrap.
150 1 : struct Solution {
151 0 : using type = std::unique_ptr<evolution::initial_data::InitialData>;
152 0 : static constexpr Options::String help =
153 : "The analytic initial data to wrap. Must be an analytic solution or "
154 : "analytic data (not numeric initial data).";
155 : };
156 :
157 : /// Amplitude \f$A\f$ of the noise interval \f$[-A,\,A]\f$.
158 1 : struct Amplitude {
159 0 : using type = double;
160 0 : static constexpr Options::String help =
161 : "Amplitude A of the noise interval [-A, A] applied "
162 : "independently to each grid point and tensor component.";
163 : };
164 :
165 : /// Whether `Amplitude` is absolute or relative to the maximum absolute
166 : /// value of each variable on each element.
167 1 : struct AmplitudeType {
168 0 : using type = NoiseAmplitudeType;
169 0 : static constexpr Options::String help =
170 : "'Absolute': noise is drawn from [-A, A] with A = Amplitude. "
171 : "'Relative': A = Amplitude * M, where M is the maximum absolute value "
172 : "of the variable over all of its components and grid points on the "
173 : "element.";
174 : };
175 :
176 : /// Seed for the random number generator.
177 1 : struct Seed {
178 0 : using type = Options::Auto<size_t, Options::AutoLabel::None>;
179 0 : static constexpr Options::String help =
180 : "Seed for the noise RNG. Set to 'None' for a random seed chosen at "
181 : "runtime (runs will not be reproducible).";
182 : };
183 :
184 : /// Names of evolution variables to perturb, or `[All]` for all.
185 1 : struct Variables {
186 0 : using type = std::vector<std::string>;
187 0 : static constexpr Options::String help =
188 : "List of variable names to add noise to, or [All] to noise every "
189 : "evolution variable. Valid names depend on the system; see class "
190 : "documentation.";
191 : };
192 :
193 0 : using options =
194 : tmpl::list<Solution, Amplitude, AmplitudeType, Seed, Variables>;
195 0 : static constexpr Options::String help =
196 : "Wraps any analytic initial data and adds uniform random noise to "
197 : "selected evolution variables after the inner solution sets them.";
198 :
199 0 : WithNoise() = default;
200 0 : WithNoise(const WithNoise& rhs);
201 0 : WithNoise& operator=(const WithNoise& rhs);
202 0 : WithNoise(WithNoise&&) = default;
203 0 : WithNoise& operator=(WithNoise&&) = default;
204 0 : ~WithNoise() override = default;
205 :
206 0 : WithNoise(std::unique_ptr<evolution::initial_data::InitialData> solution,
207 : double amplitude, NoiseAmplitudeType amplitude_type,
208 : std::optional<size_t> seed, std::vector<std::string> variables);
209 :
210 : /// \cond
211 : explicit WithNoise(CkMigrateMessage* msg);
212 : using PUP::able::register_constructor;
213 : WRAPPED_PUPable_decl_template(WithNoise); // NOLINT
214 : /// \endcond
215 :
216 0 : std::unique_ptr<evolution::initial_data::InitialData> get_clone()
217 : const override;
218 :
219 : // NOLINTNEXTLINE(google-runtime-references)
220 0 : void pup(PUP::er& p) override;
221 :
222 0 : const evolution::initial_data::InitialData& solution() const {
223 : return *solution_;
224 : }
225 : /// Unwrap to the inner analytic solution, recursively resolving nested
226 : /// `WithNoise` wrappers.
227 1 : const evolution::initial_data::InitialData& unwrap() const override {
228 : return solution_->unwrap();
229 : }
230 0 : double amplitude() const { return amplitude_; }
231 0 : NoiseAmplitudeType amplitude_type() const { return amplitude_type_; }
232 0 : size_t seed() const { return seed_; }
233 0 : const std::vector<std::string>& variables() const { return variables_; }
234 :
235 0 : friend bool operator==(const WithNoise& lhs, const WithNoise& rhs);
236 :
237 0 : friend bool operator!=(const WithNoise& lhs, const WithNoise& rhs);
238 :
239 : private:
240 0 : std::unique_ptr<evolution::initial_data::InitialData> solution_;
241 0 : double amplitude_{0.};
242 0 : NoiseAmplitudeType amplitude_type_{NoiseAmplitudeType::Absolute};
243 0 : size_t seed_{0};
244 0 : std::vector<std::string> variables_;
245 : };
246 :
247 : } // namespace evolution::initial_data
248 :
249 : template <>
250 0 : struct Options::create_from_yaml<evolution::initial_data::NoiseAmplitudeType> {
251 : template <typename Metavariables>
252 0 : static evolution::initial_data::NoiseAmplitudeType create(
253 : const Options::Option& options) {
254 : return create<void>(options);
255 : }
256 : };
257 :
258 : template <>
259 0 : evolution::initial_data::NoiseAmplitudeType
260 : Options::create_from_yaml<evolution::initial_data::NoiseAmplitudeType>::create<
261 : void>(const Options::Option& options);
|
