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