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