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