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 <cstdint>
9 : #include <optional>
10 : #include <unordered_map>
11 : #include <utility>
12 :
13 : #include "DataStructures/ApplyMatrices.hpp"
14 : #include "DataStructures/DataBox/DataBox.hpp"
15 : #include "DataStructures/DataBox/PrefixHelpers.hpp"
16 : #include "DataStructures/DataBox/Prefixes.hpp"
17 : #include "Domain/Amr/Flag.hpp"
18 : #include "Domain/Amr/Helpers.hpp"
19 : #include "Domain/Amr/Info.hpp"
20 : #include "Domain/Amr/Tags/Flags.hpp"
21 : #include "Domain/Structure/ElementId.hpp"
22 : #include "Domain/Tags.hpp"
23 : #include "Evolution/Initialization/Tags.hpp"
24 : #include "NumericalAlgorithms/Spectral/Mesh.hpp"
25 : #include "NumericalAlgorithms/Spectral/Projection.hpp"
26 : #include "Parallel/Tags/ArrayIndex.hpp"
27 : #include "ParallelAlgorithms/Amr/Protocols/Projector.hpp"
28 : #include "Time/AdaptiveSteppingDiagnostics.hpp"
29 : #include "Time/ChooseLtsStepSize.hpp"
30 : #include "Time/Slab.hpp"
31 : #include "Time/StepChoosers/StepChooser.hpp"
32 : #include "Time/Tags/AdaptiveSteppingDiagnostics.hpp"
33 : #include "Time/Tags/HistoryEvolvedVariables.hpp"
34 : #include "Time/Tags/StepChoosers.hpp"
35 : #include "Time/Tags/Time.hpp"
36 : #include "Time/Tags/TimeStep.hpp"
37 : #include "Time/Tags/TimeStepId.hpp"
38 : #include "Time/Tags/TimeStepper.hpp"
39 : #include "Time/Time.hpp"
40 : #include "Time/TimeStepId.hpp"
41 : #include "Time/TimeSteppers/LtsTimeStepper.hpp"
42 : #include "Time/TimeSteppers/TimeStepper.hpp"
43 : #include "Utilities/Algorithm.hpp"
44 : #include "Utilities/ErrorHandling/Error.hpp"
45 : #include "Utilities/Gsl.hpp"
46 : #include "Utilities/TMPL.hpp"
47 : #include "Utilities/TaggedTuple.hpp"
48 :
49 : namespace Initialization {
50 :
51 : namespace detail {
52 : inline Time initial_time(const bool time_runs_forward,
53 : const double initial_time_value,
54 : const double initial_slab_size) {
55 : const Slab initial_slab =
56 : time_runs_forward
57 : ? Slab::with_duration_from_start(initial_time_value,
58 : initial_slab_size)
59 : : Slab::with_duration_to_end(initial_time_value, initial_slab_size);
60 : return time_runs_forward ? initial_slab.start() : initial_slab.end();
61 : }
62 :
63 : template <typename TimeStepper>
64 : void set_next_time_step_id(const gsl::not_null<TimeStepId*> next_time_step_id,
65 : const Time& initial_time,
66 : const bool time_runs_forward,
67 : const TimeStepper& time_stepper) {
68 : *next_time_step_id = TimeStepId(
69 : time_runs_forward,
70 : -static_cast<int64_t>(time_stepper.number_of_past_steps()), initial_time);
71 : }
72 : } // namespace detail
73 :
74 : /// \ingroup InitializationGroup
75 : /// \brief Initialize items related to time stepping
76 : ///
77 : /// \details See the type aliases defined below for what items are added to the
78 : /// GlobalCache and DataBox and how they are initialized
79 : ///
80 : /// Since the evolution has not started yet, initialize the state
81 : /// _before_ the initial time. So `Tags::TimeStepId` is undefined at this point,
82 : /// and `Tags::Next<Tags::TimeStepId>` is the initial time.
83 : template <typename Metavariables, typename TimeStepperBase>
84 1 : struct TimeStepping {
85 : /// Tags for constant items added to the GlobalCache. These items are
86 : /// initialized from input file options.
87 1 : using const_global_cache_tags = tmpl::conditional_t<
88 : TimeStepperBase::local_time_stepping,
89 : tmpl::list<::Tags::ConcreteTimeStepper<TimeStepperBase>,
90 : ::Tags::StepChoosers>,
91 : tmpl::list<::Tags::ConcreteTimeStepper<TimeStepperBase>>>;
92 :
93 : /// Tags for mutable items added to the GlobalCache. These items are
94 : /// initialized from input file options.
95 1 : using mutable_global_cache_tags = tmpl::list<>;
96 :
97 : /// Tags for items fetched by the DataBox and passed to the apply function
98 1 : using argument_tags =
99 : tmpl::list<::Tags::Time, Tags::InitialTimeDelta,
100 : Tags::InitialSlabSize<TimeStepperBase::local_time_stepping>,
101 : ::Tags::ConcreteTimeStepper<TimeStepperBase>>;
102 :
103 : /// Tags for simple DataBox items that are initialized from input file options
104 1 : using simple_tags_from_options =
105 : tmpl::list<::Tags::Time, Tags::InitialTimeDelta,
106 : Tags::InitialSlabSize<TimeStepperBase::local_time_stepping>>;
107 :
108 : /// Tags for simple DataBox items that are default initialized.
109 1 : using default_initialized_simple_tags =
110 : tmpl::push_back<StepChoosers::step_chooser_simple_tags<
111 : Metavariables, TimeStepperBase::local_time_stepping>,
112 : ::Tags::TimeStepId, ::Tags::AdaptiveSteppingDiagnostics>;
113 :
114 : /// Tags for items in the DataBox that are mutated by the apply function
115 1 : using return_tags =
116 : tmpl::list<::Tags::Next<::Tags::TimeStepId>, ::Tags::TimeStep,
117 : ::Tags::Next<::Tags::TimeStep>>;
118 :
119 : /// Tags for mutable DataBox items that are either default initialized or
120 : /// initialized by the apply function
121 1 : using simple_tags =
122 : tmpl::append<default_initialized_simple_tags, return_tags>;
123 :
124 : /// Tags for immutable DataBox items (compute items or reference items) added
125 : /// to the DataBox.
126 1 : using compute_tags = time_stepper_ref_tags<TimeStepperBase>;
127 :
128 : /// Given the items fetched from a DataBox by the argument_tags when using
129 : /// LTS, mutate the items in the DataBox corresponding to return_tags
130 1 : static void apply(const gsl::not_null<TimeStepId*> next_time_step_id,
131 : const gsl::not_null<TimeDelta*> time_step,
132 : const gsl::not_null<TimeDelta*> next_time_step,
133 : const double initial_time_value,
134 : const double initial_dt_value,
135 : const double initial_slab_size,
136 : const LtsTimeStepper& time_stepper) {
137 : const bool time_runs_forward = initial_dt_value > 0.0;
138 : const Time initial_time = detail::initial_time(
139 : time_runs_forward, initial_time_value, initial_slab_size);
140 : detail::set_next_time_step_id(next_time_step_id, initial_time,
141 : time_runs_forward, time_stepper);
142 : *time_step = choose_lts_step_size(initial_time, initial_dt_value);
143 : *next_time_step = *time_step;
144 : }
145 :
146 : /// Given the items fetched from a DataBox by the argument_tags, when not
147 : /// using LTS, mutate the items in the DataBox corresponding to return_tags
148 1 : static void apply(const gsl::not_null<TimeStepId*> next_time_step_id,
149 : const gsl::not_null<TimeDelta*> time_step,
150 : const gsl::not_null<TimeDelta*> next_time_step,
151 : const double initial_time_value,
152 : const double initial_dt_value,
153 : const double initial_slab_size,
154 : const TimeStepper& time_stepper) {
155 : const bool time_runs_forward = initial_dt_value > 0.0;
156 : const Time initial_time = detail::initial_time(
157 : time_runs_forward, initial_time_value, initial_slab_size);
158 : detail::set_next_time_step_id(next_time_step_id, initial_time,
159 : time_runs_forward, time_stepper);
160 : *time_step = (time_runs_forward ? 1 : -1) * initial_time.slab().duration();
161 : *next_time_step = *time_step;
162 : }
163 : };
164 :
165 : /// \brief Initialize/update items related to time stepping after an AMR change
166 : template <size_t Dim>
167 1 : struct ProjectTimeStepping : tt::ConformsTo<amr::protocols::Projector> {
168 0 : using return_tags =
169 : tmpl::list<::Tags::TimeStepId, ::Tags::Next<::Tags::TimeStepId>,
170 : ::Tags::TimeStep, ::Tags::Next<::Tags::TimeStep>, ::Tags::Time,
171 : ::Tags::AdaptiveSteppingDiagnostics>;
172 0 : using argument_tags = tmpl::list<Parallel::Tags::ArrayIndex>;
173 :
174 0 : static void apply(
175 : const gsl::not_null<TimeStepId*> /*time_step_id*/,
176 : const gsl::not_null<TimeStepId*> /*next_time_step_id*/,
177 : const gsl::not_null<TimeDelta*> /*time_step*/,
178 : const gsl::not_null<TimeDelta*> /*next_time_step*/,
179 : const gsl::not_null<double*> /*time*/,
180 : const gsl::not_null<AdaptiveSteppingDiagnostics*>
181 : /*adaptive_stepping_diagnostics*/,
182 : const ElementId<Dim>& /*element_id*/,
183 : const std::pair<Mesh<Dim>, Element<Dim>>& /*old_mesh_and_element*/) {
184 : // Do not change anything for p-refinement
185 : }
186 :
187 : template <typename... Tags>
188 0 : static void apply(const gsl::not_null<TimeStepId*> time_step_id,
189 : const gsl::not_null<TimeStepId*> next_time_step_id,
190 : const gsl::not_null<TimeDelta*> time_step,
191 : const gsl::not_null<TimeDelta*> next_time_step,
192 : const gsl::not_null<double*> time,
193 : const gsl::not_null<AdaptiveSteppingDiagnostics*>
194 : adaptive_stepping_diagnostics,
195 : const ElementId<Dim>& element_id,
196 : const tuples::TaggedTuple<Tags...>& parent_items) {
197 : *time_step_id = get<::Tags::TimeStepId>(parent_items);
198 : *next_time_step_id = get<::Tags::Next<::Tags::TimeStepId>>(parent_items);
199 : *time_step = get<::Tags::TimeStep>(parent_items);
200 : *next_time_step = get<::Tags::Next<::Tags::TimeStep>>(parent_items);
201 : *time = get<::Tags::Time>(parent_items);
202 :
203 : // Since AdaptiveSteppingDiagnostics are reduced over all elements, we
204 : // set the slab quantities to the same value over all children, and the
205 : // step quantities to belong to the first child
206 : const auto& parent_diagnostics =
207 : get<::Tags::AdaptiveSteppingDiagnostics>(parent_items);
208 : const auto& parent_amr_flags =
209 : get<amr::Tags::Info<Dim>>(parent_items).flags;
210 : const auto& parent_id =
211 : get<Parallel::Tags::ArrayIndexImpl<ElementId<Dim>>>(parent_items);
212 : auto children_ids = amr::ids_of_children(parent_id, parent_amr_flags);
213 : if (element_id == children_ids.front()) {
214 : *adaptive_stepping_diagnostics = parent_diagnostics;
215 : } else {
216 : adaptive_stepping_diagnostics->number_of_slabs =
217 : parent_diagnostics.number_of_slabs;
218 : adaptive_stepping_diagnostics->number_of_slab_size_changes =
219 : parent_diagnostics.number_of_slab_size_changes;
220 : }
221 : }
222 :
223 : template <typename... Tags>
224 0 : static void apply(
225 : const gsl::not_null<TimeStepId*> time_step_id,
226 : const gsl::not_null<TimeStepId*> next_time_step_id,
227 : const gsl::not_null<TimeDelta*> time_step,
228 : const gsl::not_null<TimeDelta*> next_time_step,
229 : const gsl::not_null<double*> time,
230 : const gsl::not_null<AdaptiveSteppingDiagnostics*>
231 : adaptive_stepping_diagnostics,
232 : const ElementId<Dim>& /*element_id*/,
233 : const std::unordered_map<ElementId<Dim>, tuples::TaggedTuple<Tags...>>&
234 : children_items) {
235 : const auto slowest_child =
236 : alg::min_element(children_items, [](const auto& a, const auto& b) {
237 : const auto& time_step_a = get<::Tags::TimeStep>(a.second);
238 : const auto& time_step_b = get<::Tags::TimeStep>(b.second);
239 : ASSERT(time_step_a.is_positive() == time_step_b.is_positive(),
240 : "Elements are not taking time steps in the same direction!");
241 : return time_step_a.is_positive() ? (time_step_a < time_step_b)
242 : : (time_step_a > time_step_b);
243 : });
244 : const auto& slowest_child_items = (*slowest_child).second;
245 : *time_step_id = get<::Tags::TimeStepId>(slowest_child_items);
246 : *next_time_step_id =
247 : get<::Tags::Next<::Tags::TimeStepId>>(slowest_child_items);
248 : *time_step = get<::Tags::TimeStep>(slowest_child_items);
249 : *next_time_step = get<::Tags::Next<::Tags::TimeStep>>(slowest_child_items);
250 : *time = get<::Tags::Time>(slowest_child_items);
251 : const auto& slowest_child_diagnostics =
252 : get<::Tags::AdaptiveSteppingDiagnostics>(slowest_child_items);
253 :
254 : adaptive_stepping_diagnostics->number_of_slabs =
255 : slowest_child_diagnostics.number_of_slabs;
256 : adaptive_stepping_diagnostics->number_of_slab_size_changes =
257 : slowest_child_diagnostics.number_of_slab_size_changes;
258 : for (const auto& [_, child_items] : children_items) {
259 : *adaptive_stepping_diagnostics +=
260 : get<::Tags::AdaptiveSteppingDiagnostics>(child_items);
261 : }
262 : }
263 : };
264 :
265 : /// \ingroup InitializationGroup
266 : /// \brief Initialize time-stepper items
267 : ///
268 : /// DataBox changes:
269 : /// - Adds:
270 : /// * `db::add_tag_prefix<Tags::dt, variables_tag>`
271 : /// * `Tags::HistoryEvolvedVariables<variables_tag, dt_variables_tag>`
272 : /// - Removes: nothing
273 : /// - Modifies: nothing
274 : ///
275 : /// \note HistoryEvolvedVariables is allocated, but needs to be initialized
276 : template <typename Metavariables>
277 1 : struct TimeStepperHistory {
278 0 : static constexpr size_t dim = Metavariables::volume_dim;
279 0 : using variables_tag = typename Metavariables::system::variables_tag;
280 0 : using dt_variables_tag = db::add_tag_prefix<::Tags::dt, variables_tag>;
281 :
282 0 : using const_global_cache_tags = tmpl::list<>;
283 0 : using mutable_global_cache_tags = tmpl::list<>;
284 0 : using simple_tags_from_options = tmpl::list<>;
285 0 : using simple_tags =
286 : tmpl::list<dt_variables_tag,
287 : ::Tags::HistoryEvolvedVariables<variables_tag>>;
288 0 : using compute_tags = tmpl::list<>;
289 :
290 0 : using argument_tags =
291 : tmpl::list<::Tags::TimeStepper<TimeStepper>, domain::Tags::Mesh<dim>>;
292 0 : using return_tags = simple_tags;
293 :
294 0 : static void apply(
295 : const gsl::not_null<typename dt_variables_tag::type*> dt_vars,
296 : const gsl::not_null<TimeSteppers::History<typename variables_tag::type>*>
297 : history,
298 : const TimeStepper& time_stepper, const Mesh<dim>& mesh) {
299 : // Will be overwritten before use
300 : dt_vars->initialize(mesh.number_of_grid_points());
301 :
302 : // All steppers we have that need to start at low order require
303 : // one additional point per order, so this is the order that
304 : // requires no initial past steps.
305 : const size_t starting_order =
306 : time_stepper.order() - time_stepper.number_of_past_steps();
307 : history->integration_order(starting_order);
308 : }
309 : };
310 :
311 : /// \brief Initialize/update items related to time stepper history after an AMR
312 : /// change
313 : template <typename Metavariables>
314 1 : struct ProjectTimeStepperHistory : tt::ConformsTo<amr::protocols::Projector> {
315 0 : static constexpr size_t dim = Metavariables::volume_dim;
316 0 : using variables_tag = typename Metavariables::system::variables_tag;
317 0 : using dt_variables_tag = db::add_tag_prefix<::Tags::dt, variables_tag>;
318 0 : using history_tag = ::Tags::HistoryEvolvedVariables<variables_tag>;
319 :
320 0 : using return_tags = tmpl::list<dt_variables_tag, history_tag>;
321 0 : using argument_tags =
322 : tmpl::list<domain::Tags::Mesh<dim>, Parallel::Tags::ArrayIndex>;
323 :
324 0 : static void apply(
325 : const gsl::not_null<typename dt_variables_tag::type*> dt_vars,
326 : const gsl::not_null<typename history_tag::type*> history,
327 : const Mesh<dim>& new_mesh, const ElementId<dim>& /*element_id*/,
328 : const std::pair<Mesh<dim>, Element<dim>>& old_mesh_and_element) {
329 : const auto& old_mesh = old_mesh_and_element.first;
330 : if (old_mesh == new_mesh) {
331 : return; // mesh was not refined, so no projection needed
332 : }
333 : const auto projection_matrices =
334 : Spectral::p_projection_matrices(old_mesh, new_mesh);
335 : const auto& old_extents = old_mesh.extents();
336 : history->map_entries(
337 : [&projection_matrices, &old_extents](const auto entry) {
338 : *entry = apply_matrices(projection_matrices, *entry, old_extents);
339 : });
340 : dt_vars->initialize(new_mesh.number_of_grid_points());
341 : }
342 :
343 : template <typename... Tags>
344 0 : static void apply(
345 : const gsl::not_null<typename dt_variables_tag::type*> /*dt_vars*/,
346 : const gsl::not_null<typename history_tag::type*> /*history*/,
347 : const Mesh<dim>& /*new_mesh*/, const ElementId<dim>& /*element_id*/,
348 : const tuples::TaggedTuple<Tags...>& /*parent_items*/) {
349 : ERROR("h-refinement not implemented yet");
350 : }
351 :
352 : template <typename... Tags>
353 0 : static void apply(
354 : const gsl::not_null<typename dt_variables_tag::type*> /*dt_vars*/,
355 : const gsl::not_null<typename history_tag::type*> /*history*/,
356 : const Mesh<dim>& /*new_mesh*/, const ElementId<dim>& /*element_id*/,
357 : const std::unordered_map<ElementId<dim>, tuples::TaggedTuple<Tags...>>&
358 : /*children_items*/) {
359 : ERROR("h-refinement not implemented yet");
360 : }
361 : };
362 : } // namespace Initialization
|
