Line data Source code
1 0 : // Distributed under the MIT License.
2 : // See LICENSE.txt for details.
3 :
4 : #pragma once
5 :
6 : #include <algorithm>
7 : #include <cmath>
8 : #include <cstddef>
9 : #include <cstdint>
10 : #include <optional>
11 : #include <type_traits>
12 : #include <unordered_map>
13 : #include <utility>
14 :
15 : #include "DataStructures/DataBox/PrefixHelpers.hpp"
16 : #include "DataStructures/DataBox/Prefixes.hpp"
17 : #include "DataStructures/TaggedTuple.hpp"
18 : #include "DataStructures/TaggedVariant.hpp"
19 : #include "Domain/Amr/Helpers.hpp"
20 : #include "Domain/Structure/ChildSize.hpp"
21 : #include "Domain/Structure/DirectionalIdMap.hpp"
22 : #include "Domain/Structure/Element.hpp"
23 : #include "Domain/Structure/ElementId.hpp"
24 : #include "Evolution/DiscontinuousGalerkin/MortarInfo.hpp"
25 : #include "Evolution/DiscontinuousGalerkin/TimeSteppingPolicy.hpp"
26 : #include "Evolution/Initialization/Tags.hpp"
27 : #include "NumericalAlgorithms/Spectral/Mesh.hpp"
28 : #include "NumericalAlgorithms/Spectral/Projection.hpp"
29 : #include "ParallelAlgorithms/Amr/Protocols/Projector.hpp"
30 : #include "Time/AdaptiveSteppingDiagnostics.hpp"
31 : #include "Time/ChangeSlabSize/Tags.hpp"
32 : #include "Time/ChooseLtsStepSize.hpp"
33 : #include "Time/History.hpp"
34 : #include "Time/LtsMode.hpp"
35 : #include "Time/Slab.hpp"
36 : #include "Time/StepChoosers/StepChooser.hpp"
37 : #include "Time/Tags/AdaptiveSteppingDiagnostics.hpp"
38 : #include "Time/Tags/HistoryEvolvedVariables.hpp"
39 : #include "Time/Tags/LtsMode.hpp"
40 : #include "Time/Tags/StepNumberWithinSlab.hpp"
41 : #include "Time/Tags/Time.hpp"
42 : #include "Time/Tags/TimeStep.hpp"
43 : #include "Time/Tags/TimeStepId.hpp"
44 : #include "Time/Tags/TimeStepper.hpp"
45 : #include "Time/Time.hpp"
46 : #include "Time/TimeStepId.hpp"
47 : #include "Time/TimeSteppers/LtsTimeStepper.hpp"
48 : #include "Time/TimeSteppers/TimeStepper.hpp"
49 : #include "Utilities/Algorithm.hpp"
50 : #include "Utilities/ErrorHandling/Assert.hpp"
51 : #include "Utilities/ErrorHandling/Error.hpp"
52 : #include "Utilities/Gsl.hpp"
53 : #include "Utilities/Literals.hpp"
54 : #include "Utilities/MakeArray.hpp"
55 : #include "Utilities/ProtocolHelpers.hpp"
56 : #include "Utilities/SetNumberOfGridPoints.hpp"
57 : #include "Utilities/TMPL.hpp"
58 :
59 : /// \cond
60 : namespace Parallel::Tags {
61 : template <typename Index>
62 : struct ArrayIndex;
63 : } // namespace Parallel::Tags
64 : namespace Tags {
65 : template <typename Tag>
66 : struct StepperErrors;
67 : } // namespace Tags
68 : namespace amr::Tags {
69 : template <size_t VolumeDim>
70 : struct Info;
71 : } // namespace amr::Tags
72 : namespace domain::Tags {
73 : template <size_t VolumeDim>
74 : struct Element;
75 : template <size_t VolumeDim>
76 : struct Mesh;
77 : } // namespace domain::Tags
78 : namespace evolution::dg::Tags {
79 : template <size_t Dim>
80 : struct MortarInfo;
81 : } // namespace evolution::dg::Tags
82 : /// \endcond
83 :
84 : namespace Initialization {
85 : /// \ingroup InitializationGroup
86 : /// \brief Initialize items related to time stepping
87 : ///
88 : /// \details See the type aliases defined below for what items are added to the
89 : /// GlobalCache and DataBox and how they are initialized
90 : ///
91 : /// Since the evolution has not started yet, initialize the state
92 : /// _before_ the initial time. So `Tags::TimeStepId` is undefined at this point,
93 : /// and `Tags::Next<Tags::TimeStepId>` is the initial time.
94 : template <typename Metavariables, typename TimeStepperBase,
95 : bool WithControlSystems, bool AllowLocalTimeStepping>
96 1 : struct TimeStepping {
97 : /// Tags for constant items added to the GlobalCache. These items are
98 : /// initialized from input file options.
99 1 : using const_global_cache_tags = tmpl::list<
100 : ::Tags::ConcreteTimeStepper<TimeStepperBase, WithControlSystems>,
101 : tmpl::conditional_t<AllowLocalTimeStepping, ::Tags::LtsMode,
102 : ::Tags::LtsModeForced<LtsMode::Off>>>;
103 :
104 : /// Tags for mutable items added to the GlobalCache. These items are
105 : /// initialized from input file options.
106 1 : using mutable_global_cache_tags = tmpl::list<>;
107 :
108 : /// Tags for items fetched by the DataBox and passed to the apply function
109 1 : using argument_tags =
110 : tmpl::list<::Tags::Time, Tags::InitialTimeDelta, Tags::InitialSlabSize,
111 : ::Tags::TimeStepper<TimeStepperBase>, ::Tags::LtsMode>;
112 :
113 : /// Tags for simple DataBox items that are initialized from input file options
114 1 : using simple_tags_from_options =
115 : tmpl::list<::Tags::Time, Tags::InitialTimeDelta, Tags::InitialSlabSize>;
116 :
117 : /// Tags for simple DataBox items that are default initialized.
118 1 : using default_initialized_simple_tags =
119 : tmpl::push_back<StepChoosers::step_chooser_simple_tags<Metavariables>,
120 : ::Tags::TimeStepId, ::Tags::StepNumberWithinSlab,
121 : ::Tags::AdaptiveSteppingDiagnostics>;
122 :
123 : /// Tags for items in the DataBox that are mutated by the apply function
124 1 : using return_tags =
125 : tmpl::list<::Tags::Next<::Tags::TimeStepId>, ::Tags::TimeStep,
126 : ::Tags::ChangeSlabSize::SlabSizeGoal>;
127 :
128 : /// Tags for mutable DataBox items that are either default initialized or
129 : /// initialized by the apply function
130 1 : using simple_tags =
131 : tmpl::append<default_initialized_simple_tags, return_tags>;
132 :
133 : /// Tags for immutable DataBox items (compute items or reference items) added
134 : /// to the DataBox.
135 1 : using compute_tags =
136 : time_stepper_ref_tags<TimeStepperBase, WithControlSystems>;
137 :
138 0 : static void apply(const gsl::not_null<TimeStepId*> next_time_step_id,
139 : const gsl::not_null<TimeDelta*> time_step,
140 : const gsl::not_null<double*> slab_size_goal,
141 : const double initial_time_value,
142 : const double initial_dt_value,
143 : const double initial_slab_size,
144 : const TimeStepper& time_stepper, const LtsMode lts_mode) {
145 : // Ignore the sign of initial_slab_size for user convenience. GTS
146 : // evolutions can set the initial slab size and time step the same
147 : // even if they are negative.
148 :
149 : const bool time_runs_forward = initial_dt_value > 0.0;
150 : const Slab initial_slab =
151 : time_runs_forward
152 : ? Slab::with_duration_from_start(initial_time_value,
153 : std::abs(initial_slab_size))
154 : : Slab::with_duration_to_end(initial_time_value,
155 : std::abs(initial_slab_size));
156 : const auto initial_time =
157 : time_runs_forward ? initial_slab.start() : initial_slab.end();
158 :
159 : *next_time_step_id =
160 : TimeStepId(time_runs_forward,
161 : -static_cast<int64_t>(time_stepper.number_of_past_steps()),
162 : initial_time);
163 : *slab_size_goal =
164 : (time_runs_forward ? 1.0 : -1.0) * std::abs(initial_slab_size);
165 :
166 : if (lts_mode == LtsMode::Off) {
167 : if (std::abs(initial_slab_size) != std::abs(initial_dt_value)) {
168 : ERROR_NO_TRACE(
169 : "InitialSlabSize and InitialTimeStep must be equal with "
170 : "LocalTimeStepping off.");
171 : }
172 : *time_step = (time_runs_forward ? 1 : -1) * initial_slab.duration();
173 : } else {
174 : *time_step = choose_lts_step_size(initial_time, initial_dt_value);
175 : }
176 : }
177 : };
178 :
179 : /// \brief Initialize/update items related to time stepping after an AMR change
180 : template <size_t Dim>
181 1 : struct ProjectTimeStepping : tt::ConformsTo<amr::protocols::Projector> {
182 0 : using return_tags =
183 : tmpl::list<::Tags::TimeStepId, ::Tags::Next<::Tags::TimeStepId>,
184 : ::Tags::TimeStep, ::Tags::Time, ::Tags::StepNumberWithinSlab,
185 : ::Tags::AdaptiveSteppingDiagnostics,
186 : ::Tags::ChangeSlabSize::SlabSizeGoal>;
187 0 : using argument_tags = tmpl::list<Parallel::Tags::ArrayIndex<ElementId<Dim>>>;
188 :
189 0 : static void apply(
190 : const gsl::not_null<TimeStepId*> /*time_step_id*/,
191 : const gsl::not_null<TimeStepId*> /*next_time_step_id*/,
192 : const gsl::not_null<TimeDelta*> /*time_step*/,
193 : const gsl::not_null<double*> /*time*/,
194 : const gsl::not_null<uint64_t*> /*step_number_within_slab*/,
195 : const gsl::not_null<AdaptiveSteppingDiagnostics*>
196 : /*adaptive_stepping_diagnostics*/,
197 : const gsl::not_null<double*> /*slab_size_goal*/,
198 : const ElementId<Dim>& /*element_id*/,
199 : const std::pair<Mesh<Dim>, Element<Dim>>& /*old_mesh_and_element*/) {
200 : // Do not change anything for p-refinement
201 : }
202 :
203 : template <typename... Tags>
204 0 : static void apply(const gsl::not_null<TimeStepId*> time_step_id,
205 : const gsl::not_null<TimeStepId*> next_time_step_id,
206 : const gsl::not_null<TimeDelta*> time_step,
207 : const gsl::not_null<double*> time,
208 : const gsl::not_null<uint64_t*> step_number_within_slab,
209 : const gsl::not_null<AdaptiveSteppingDiagnostics*>
210 : adaptive_stepping_diagnostics,
211 : const gsl::not_null<double*> slab_size_goal,
212 : const ElementId<Dim>& element_id,
213 : const tuples::TaggedTuple<Tags...>& parent_items) {
214 : *time_step_id = get<::Tags::TimeStepId>(parent_items);
215 : *next_time_step_id = get<::Tags::Next<::Tags::TimeStepId>>(parent_items);
216 : *time_step = get<::Tags::TimeStep>(parent_items);
217 : *time = get<::Tags::Time>(parent_items);
218 : *slab_size_goal = get<::Tags::ChangeSlabSize::SlabSizeGoal>(parent_items);
219 : *step_number_within_slab = get<::Tags::StepNumberWithinSlab>(parent_items);
220 :
221 : // Since AdaptiveSteppingDiagnostics are reduced over all elements, we
222 : // set the slab quantities to the same value over all children, and the
223 : // step quantities to belong to the first child
224 : const auto& parent_diagnostics =
225 : get<::Tags::AdaptiveSteppingDiagnostics>(parent_items);
226 : const auto& parent_amr_flags =
227 : get<amr::Tags::Info<Dim>>(parent_items).flags;
228 : const auto& parent_id =
229 : get<Parallel::Tags::ArrayIndex<ElementId<Dim>>>(parent_items);
230 : auto children_ids = amr::ids_of_children(parent_id, parent_amr_flags);
231 : if (element_id == children_ids.front()) {
232 : *adaptive_stepping_diagnostics = parent_diagnostics;
233 : } else {
234 : adaptive_stepping_diagnostics->number_of_slabs =
235 : parent_diagnostics.number_of_slabs;
236 : adaptive_stepping_diagnostics->number_of_slab_size_changes =
237 : parent_diagnostics.number_of_slab_size_changes;
238 : }
239 : }
240 :
241 : template <typename... Tags>
242 0 : static void apply(
243 : const gsl::not_null<TimeStepId*> time_step_id,
244 : const gsl::not_null<TimeStepId*> next_time_step_id,
245 : const gsl::not_null<TimeDelta*> time_step,
246 : const gsl::not_null<double*> time,
247 : const gsl::not_null<uint64_t*> step_number_within_slab,
248 : const gsl::not_null<AdaptiveSteppingDiagnostics*>
249 : adaptive_stepping_diagnostics,
250 : const gsl::not_null<double*> slab_size_goal,
251 : const ElementId<Dim>& /*element_id*/,
252 : const std::unordered_map<ElementId<Dim>, tuples::TaggedTuple<Tags...>>&
253 : children_items) {
254 : const auto slowest_child =
255 : alg::min_element(children_items, [](const auto& a, const auto& b) {
256 : const auto& time_step_a = get<::Tags::TimeStep>(a.second);
257 : const auto& time_step_b = get<::Tags::TimeStep>(b.second);
258 : ASSERT(time_step_a.is_positive() == time_step_b.is_positive(),
259 : "Elements are not taking time steps in the same direction!");
260 : return time_step_a.is_positive() ? (time_step_a < time_step_b)
261 : : (time_step_a > time_step_b);
262 : });
263 : const auto& slowest_child_items = (*slowest_child).second;
264 : *time_step_id = get<::Tags::TimeStepId>(slowest_child_items);
265 : *next_time_step_id =
266 : get<::Tags::Next<::Tags::TimeStepId>>(slowest_child_items);
267 : *time_step = get<::Tags::TimeStep>(slowest_child_items);
268 : *time = get<::Tags::Time>(slowest_child_items);
269 : *slab_size_goal =
270 : get<::Tags::ChangeSlabSize::SlabSizeGoal>(slowest_child_items);
271 : *step_number_within_slab =
272 : get<::Tags::StepNumberWithinSlab>(slowest_child_items);
273 : const auto& slowest_child_diagnostics =
274 : get<::Tags::AdaptiveSteppingDiagnostics>(slowest_child_items);
275 :
276 : adaptive_stepping_diagnostics->number_of_slabs =
277 : slowest_child_diagnostics.number_of_slabs;
278 : adaptive_stepping_diagnostics->number_of_slab_size_changes =
279 : slowest_child_diagnostics.number_of_slab_size_changes;
280 : for (const auto& [_, child_items] : children_items) {
281 : *adaptive_stepping_diagnostics +=
282 : get<::Tags::AdaptiveSteppingDiagnostics>(child_items);
283 : }
284 : }
285 : };
286 :
287 : /// \ingroup InitializationGroup
288 : /// \brief Initialize time-stepper items
289 : ///
290 : /// DataBox changes:
291 : /// - Adds:
292 : /// * `db::add_tag_prefix<Tags::dt, variables_tag>`
293 : /// * `Tags::HistoryEvolvedVariables<variables_tag>`
294 : /// - Removes: nothing
295 : /// - Modifies: nothing
296 : ///
297 : /// \note HistoryEvolvedVariables is allocated, but needs to be initialized
298 : /// @{
299 : template <typename System,
300 : template <typename> typename CacheTagPrefix = std::type_identity_t,
301 : typename = tmpl::conditional_t<
302 : tt::is_a_v<tmpl::list, typename System::variables_tag>,
303 : typename System::variables_tag,
304 : tmpl::list<typename System::variables_tag>>>
305 1 : struct TimeStepperHistory;
306 :
307 : template <typename System, template <typename> typename CacheTagPrefix,
308 : typename... VariablesTags>
309 0 : struct TimeStepperHistory<System, CacheTagPrefix,
310 : tmpl::list<VariablesTags...>> {
311 0 : using const_global_cache_tags = tmpl::list<>;
312 0 : using mutable_global_cache_tags = tmpl::list<>;
313 0 : using simple_tags_from_options = tmpl::list<>;
314 0 : using simple_tags =
315 : tmpl::list<db::add_tag_prefix<::Tags::dt, VariablesTags>...,
316 : ::Tags::HistoryEvolvedVariables<VariablesTags>...>;
317 0 : using compute_tags = tmpl::list<>;
318 :
319 0 : using argument_tags =
320 : tmpl::list<CacheTagPrefix<::Tags::TimeStepper<TimeStepper>>,
321 : VariablesTags...>;
322 0 : using return_tags = simple_tags;
323 :
324 0 : static void apply(
325 : const gsl::not_null<typename db::add_tag_prefix<
326 : ::Tags::dt, VariablesTags>::type*>... dt_vars,
327 : const gsl::not_null<
328 : TimeSteppers::History<typename VariablesTags::type>*>... history,
329 : const TimeStepper& time_stepper,
330 : const typename VariablesTags::type&... vars) {
331 : // Will be overwritten before use
332 : expand_pack((set_number_of_grid_points(dt_vars, vars), 0)...);
333 :
334 : // All steppers we have that need to start at low order require
335 : // one additional point per order, so this is the order that
336 : // requires no initial past steps.
337 : const size_t starting_order =
338 : visit(
339 : []<typename Tag>(
340 : const std::pair<tmpl::type_<Tag>, typename Tag::type&&> order) {
341 : if constexpr (std::is_same_v<Tag,
342 : TimeSteppers::Tags::FixedOrder>) {
343 : return order.second;
344 : } else {
345 : return order.second.minimum;
346 : }
347 : },
348 : time_stepper.order()) -
349 : time_stepper.number_of_past_steps();
350 : expand_pack((history->integration_order(starting_order), 0)...);
351 : }
352 : };
353 : /// @}
354 :
355 : /// \brief Initialize/update items related to time stepper history after an AMR
356 : /// change
357 : ///
358 : /// \note `Tags::TimeStep` and `Tags::Next<Tags::TimeStepId>` are not
359 : /// initially set by this projector. They are only updated if the
360 : /// time stepper must be restarted because of LTS h-refinement.
361 : template <typename Metavariables>
362 1 : struct ProjectTimeStepperHistory : tt::ConformsTo<amr::protocols::Projector> {
363 0 : static constexpr size_t dim = Metavariables::volume_dim;
364 0 : using variables_tag = typename Metavariables::system::variables_tag;
365 0 : using dt_variables_tag = db::add_tag_prefix<::Tags::dt, variables_tag>;
366 0 : using history_tag = ::Tags::HistoryEvolvedVariables<variables_tag>;
367 :
368 0 : using return_tags =
369 : tmpl::list<dt_variables_tag, history_tag,
370 : ::Tags::Next<::Tags::TimeStepId>, ::Tags::TimeStep>;
371 0 : using argument_tags = tmpl::list<
372 : domain::Tags::Mesh<dim>, Parallel::Tags::ArrayIndex<ElementId<dim>>,
373 : ::Tags::TimeStepper<TimeStepper>, evolution::dg::Tags::MortarInfo<dim>>;
374 :
375 0 : static void apply(
376 : const gsl::not_null<typename dt_variables_tag::type*> dt_vars,
377 : const gsl::not_null<typename history_tag::type*> history,
378 : const gsl::not_null<TimeStepId*> /*next_time_step_id*/,
379 : const gsl::not_null<TimeDelta*> /*time_step*/, const Mesh<dim>& new_mesh,
380 : const ElementId<dim>& /*element_id*/, const TimeStepper& /*time_stepper*/,
381 : const DirectionalIdMap<dim, evolution::dg::MortarInfo<dim>>&
382 : /*mortar_info*/,
383 : const std::pair<Mesh<dim>, Element<dim>>& old_mesh_and_element) {
384 : const auto& old_mesh = old_mesh_and_element.first;
385 : if (old_mesh == new_mesh) {
386 : return; // mesh was not refined, so no projection needed
387 : }
388 : history->map_entries([&old_mesh, &new_mesh](const auto entry) {
389 : *entry = Spectral::project(*entry, old_mesh, new_mesh,
390 : make_array<dim>(Spectral::SegmentSize::Full),
391 : make_array<dim>(Spectral::SegmentSize::Full));
392 : });
393 : dt_vars->initialize(new_mesh.number_of_grid_points());
394 : }
395 :
396 : template <typename... Tags>
397 0 : static void apply(
398 : const gsl::not_null<typename dt_variables_tag::type*> dt_vars,
399 : const gsl::not_null<typename history_tag::type*> history,
400 : const gsl::not_null<TimeStepId*> next_time_step_id,
401 : const gsl::not_null<TimeDelta*> time_step, const Mesh<dim>& new_mesh,
402 : const ElementId<dim>& element_id, const TimeStepper& time_stepper,
403 : const DirectionalIdMap<dim, evolution::dg::MortarInfo<dim>>& mortar_info,
404 : const tuples::TaggedTuple<Tags...>& parent_items) {
405 : dt_vars->initialize(new_mesh.number_of_grid_points());
406 : ASSERT(element_id.refinement_levels() == make_array<dim>(0_st) or
407 : not mortar_info.empty(),
408 : "Element has no neighbors but is not a full block");
409 : const auto time_stepping_policy = time_stepping_policy_for_h_refinement(
410 : mortar_info, &get<evolution::dg::Tags::MortarInfo<dim>>(parent_items),
411 : {});
412 : switch (time_stepping_policy) {
413 : case evolution::dg::TimeSteppingPolicy::Conservative: {
414 : if (time_stepper.number_of_past_steps() != 0) {
415 : ERROR_NO_TRACE(
416 : "Cannot perform h-refinement with LTS steppers requiring "
417 : "initialization.");
418 : }
419 : const auto integrator_order = time_stepper.order();
420 : if (variants::holds_alternative<TimeSteppers::Tags::FixedOrder>(
421 : integrator_order)) {
422 : *history = typename history_tag::type{
423 : get<TimeSteppers::Tags::FixedOrder>(integrator_order)};
424 : return;
425 : }
426 : const auto start_order =
427 : get<TimeSteppers::Tags::VariableOrder>(integrator_order).minimum;
428 : *history = typename history_tag::type{start_order};
429 :
430 : const auto reduced_step = restart_time_step(parent_items, start_order);
431 : if (abs(reduced_step) < abs(*time_step)) {
432 : *time_step = reduced_step;
433 : *next_time_step_id = time_stepper.next_time_id(
434 : get<::Tags::TimeStepId>(parent_items), *time_step);
435 : }
436 : break;
437 : }
438 : case evolution::dg::TimeSteppingPolicy::EqualRate: {
439 : const auto& parent_id =
440 : get<domain::Tags::Element<dim>>(parent_items).id();
441 : const auto& parent_mesh = get<domain::Tags::Mesh<dim>>(parent_items);
442 : const auto child_sizes = domain::child_size(element_id.segment_ids(),
443 : parent_id.segment_ids());
444 : transform(history, get<history_tag>(parent_items),
445 : [&](const auto& source_entry) {
446 : return Spectral::project(
447 : source_entry, parent_mesh, new_mesh,
448 : make_array<dim>(Spectral::SegmentSize::Full),
449 : child_sizes);
450 : });
451 : break;
452 : }
453 : default:
454 : ERROR("Unhandled time-stepping policy: " << time_stepping_policy);
455 : }
456 : }
457 :
458 : template <typename... Tags>
459 0 : static void apply(
460 : const gsl::not_null<typename dt_variables_tag::type*> dt_vars,
461 : const gsl::not_null<typename history_tag::type*> history,
462 : const gsl::not_null<TimeStepId*> next_time_step_id,
463 : const gsl::not_null<TimeDelta*> time_step, const Mesh<dim>& new_mesh,
464 : const ElementId<dim>& element_id, const TimeStepper& time_stepper,
465 : const DirectionalIdMap<dim, evolution::dg::MortarInfo<dim>>& mortar_info,
466 : const std::unordered_map<ElementId<dim>, tuples::TaggedTuple<Tags...>>&
467 : children_items) {
468 : dt_vars->initialize(new_mesh.number_of_grid_points());
469 : ASSERT(element_id.refinement_levels() == make_array<dim>(0_st) or
470 : not mortar_info.empty(),
471 : "Element has no neighbors but is not a full block");
472 : const auto time_stepping_policy =
473 : time_stepping_policy_for_h_refinement<Tags...>(mortar_info, {},
474 : {&children_items});
475 : switch (time_stepping_policy) {
476 : case evolution::dg::TimeSteppingPolicy::Conservative: {
477 : if (time_stepper.number_of_past_steps() != 0) {
478 : ERROR_NO_TRACE(
479 : "Cannot perform h-refinement with LTS steppers requiring "
480 : "initialization.");
481 : }
482 : const auto integrator_order = time_stepper.order();
483 : if (variants::holds_alternative<TimeSteppers::Tags::FixedOrder>(
484 : integrator_order)) {
485 : *history = typename history_tag::type{
486 : get<TimeSteppers::Tags::FixedOrder>(integrator_order)};
487 : return;
488 : }
489 : const auto start_order =
490 : get<TimeSteppers::Tags::VariableOrder>(integrator_order).minimum;
491 : *history = typename history_tag::type{start_order};
492 :
493 : for (const auto& [child, child_items] : children_items) {
494 : const auto reduced_step = restart_time_step(child_items, start_order);
495 : if (abs(reduced_step) < abs(*time_step)) {
496 : *time_step = reduced_step;
497 : *next_time_step_id = time_stepper.next_time_id(
498 : get<::Tags::TimeStepId>(children_items.begin()->second),
499 : *time_step);
500 : }
501 : }
502 : break;
503 : }
504 : case evolution::dg::TimeSteppingPolicy::EqualRate: {
505 : bool first_child = true;
506 : for (const auto& [child_id, child_items] : children_items) {
507 : const auto& child_mesh = get<domain::Tags::Mesh<dim>>(child_items);
508 : const auto child_sizes = domain::child_size(child_id.segment_ids(),
509 : element_id.segment_ids());
510 : if (first_child) {
511 : transform(history, get<history_tag>(child_items),
512 : [&](const auto& source_entry) {
513 : return Spectral::project(
514 : source_entry, child_mesh, new_mesh, child_sizes,
515 : make_array<dim>(Spectral::SegmentSize::Full));
516 : });
517 : first_child = false;
518 : } else {
519 : transform_mutate(
520 : history, get<history_tag>(child_items),
521 : [&](const auto dest_entry, const auto& source_entry) {
522 : *dest_entry += Spectral::project(
523 : source_entry, child_mesh, new_mesh, child_sizes,
524 : make_array<dim>(Spectral::SegmentSize::Full));
525 : });
526 : }
527 : }
528 : break;
529 : }
530 : default:
531 : ERROR("Unhandled time-stepping policy: " << time_stepping_policy);
532 : }
533 : }
534 :
535 : private:
536 : template <typename... Tags>
537 : static evolution::dg::TimeSteppingPolicy
538 0 : time_stepping_policy_for_h_refinement(
539 : const DirectionalIdMap<dim, evolution::dg::MortarInfo<dim>>&
540 : element_infos,
541 : const std::optional<gsl::not_null<
542 : const DirectionalIdMap<dim, evolution::dg::MortarInfo<dim>>*>>&
543 : parent_infos,
544 : const std::optional<gsl::not_null<const std::unordered_map<
545 : ElementId<dim>, tuples::TaggedTuple<Tags...>>*>>& children_items) {
546 : std::optional<evolution::dg::TimeSteppingPolicy> policy{};
547 : const auto process_infos =
548 : [&policy](const DirectionalIdMap<dim, evolution::dg::MortarInfo<dim>>&
549 : infos) {
550 : for (const auto& info : infos) {
551 : const auto mortar_policy = info.second.time_stepping_policy();
552 : ASSERT(not policy.has_value() or *policy == mortar_policy,
553 : "Inconsistent policies: "
554 : << *policy << " and " << mortar_policy
555 : << "\nWe currently require all mortars to have the "
556 : "same policy when doing h-refinement. This might "
557 : "be relaxed when we add an LTS policy other than "
558 : "Conservative.");
559 : policy.emplace(mortar_policy);
560 : }
561 : };
562 :
563 : process_infos(element_infos);
564 : if (parent_infos.has_value()) {
565 : process_infos(**parent_infos);
566 : }
567 : if (children_items.has_value()) {
568 : if constexpr (sizeof...(Tags) > 0) {
569 : for (const auto& child : **children_items) {
570 : process_infos(
571 : get<evolution::dg::Tags::MortarInfo<dim>>(child.second));
572 : }
573 : } else {
574 : ERROR("Children but no child data");
575 : }
576 : }
577 :
578 : ASSERT(policy.has_value(),
579 : "Found no mortars, either before or after h-refinement. A mortar "
580 : "should have been created or destroyed, so this is not possible.");
581 : return *policy;
582 : }
583 :
584 : template <typename... Tags>
585 0 : static TimeDelta restart_time_step(
586 : const tuples::TaggedTuple<Tags...>& old_items, const size_t start_order) {
587 : const auto& old_history = get<history_tag>(old_items);
588 : if (old_history.integration_order() == start_order) {
589 : return get<::Tags::TimeStep>(old_items);
590 : }
591 :
592 : const auto& time_step_id = get<::Tags::TimeStepId>(old_items);
593 : const auto& errors = get<::Tags::StepperErrors<variables_tag>>(old_items);
594 : if (not errors[1].has_value()) {
595 : ERROR_NO_TRACE(
596 : "Evolutions performing h-refinement with variable-order local "
597 : "time-stepping must use ErrorControl for step size choosing.");
598 : }
599 : ASSERT(errors[1]->errors[start_order - 1].has_value(),
600 : "Start-order estimate not available.");
601 :
602 : // At this low an order, we are almost certainly step-size-limited
603 : // by accuracy rather than stability, so ignore things like the
604 : // change to the grid spacing.
605 : return choose_lts_step_size(
606 : time_step_id.step_time(),
607 : errors[1]->step_size.value() *
608 : pow(1.0 / std::max(*errors[1]->errors[start_order - 1], 1e-14),
609 : 1.0 / static_cast<double>(start_order)));
610 : }
611 : };
612 : } // namespace Initialization
|
