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