Line data Source code
1 0 : // Distributed under the MIT License.
2 : // See LICENSE.txt for details.
3 :
4 : #pragma once
5 :
6 : #include <array>
7 : #include <cstddef>
8 : #include <memory>
9 : #include <optional>
10 : #include <tuple>
11 : #include <utility>
12 : #include <vector>
13 :
14 : #include "ControlSystem/Tags/FunctionsOfTimeInitialize.hpp"
15 : #include "DataStructures/DataBox/DataBox.hpp"
16 : #include "DataStructures/TaggedTuple.hpp"
17 : #include "DataStructures/Tensor/IndexType.hpp"
18 : #include "Domain/CoordinateMaps/CoordinateMap.hpp"
19 : #include "Domain/CoordinateMaps/Identity.hpp"
20 : #include "Domain/CoordinateMaps/Tags.hpp"
21 : #include "Domain/CreateInitialElement.hpp"
22 : #include "Domain/Creators/Tags/Domain.hpp"
23 : #include "Domain/Creators/Tags/FunctionsOfTime.hpp"
24 : #include "Domain/Creators/Tags/InitialExtents.hpp"
25 : #include "Domain/Creators/Tags/InitialRefinementLevels.hpp"
26 : #include "Domain/Domain.hpp"
27 : #include "Domain/ElementMap.hpp"
28 : #include "Domain/MinimumGridSpacing.hpp"
29 : #include "Domain/Structure/CreateInitialMesh.hpp"
30 : #include "Domain/Structure/Direction.hpp"
31 : #include "Domain/Structure/Element.hpp"
32 : #include "Domain/Structure/ElementId.hpp"
33 : #include "Domain/Structure/IsValidDgMesh.hpp"
34 : #include "Domain/Tags.hpp"
35 : #include "Domain/Tags/NeighborMesh.hpp"
36 : #include "Domain/TagsTimeDependent.hpp"
37 : #include "Evolution/DiscontinuousGalerkin/Initialization/QuadratureTag.hpp"
38 : #include "Evolution/TagsDomain.hpp"
39 : #include "NumericalAlgorithms/Spectral/Basis.hpp"
40 : #include "NumericalAlgorithms/Spectral/LogicalCoordinates.hpp"
41 : #include "NumericalAlgorithms/Spectral/Mesh.hpp"
42 : #include "NumericalAlgorithms/Spectral/Quadrature.hpp"
43 : #include "Parallel/AlgorithmExecution.hpp"
44 : #include "Parallel/GlobalCache.hpp"
45 : #include "Parallel/Tags/ArrayIndex.hpp"
46 : #include "ParallelAlgorithms/Amr/Protocols/Projector.hpp"
47 : #include "ParallelAlgorithms/Initialization/MutateAssign.hpp"
48 : #include "Utilities/ErrorHandling/Error.hpp"
49 : #include "Utilities/Gsl.hpp"
50 : #include "Utilities/TMPL.hpp"
51 :
52 : /// \cond
53 : namespace Frame {
54 : struct Inertial;
55 : } // namespace Frame
56 : /// \endcond
57 :
58 : namespace evolution::dg::Initialization {
59 :
60 : /// \ingroup InitializationGroup
61 : /// \brief Initialize items related to the basic structure of the element
62 : ///
63 : /// \details See the type aliases defined below for what items are added to the
64 : /// GlobalCache and DataBox and how they are initialized
65 :
66 : template <typename Metavariables, bool UseControlSystems = false>
67 1 : struct Domain {
68 0 : static constexpr size_t dim = Metavariables::volume_dim;
69 :
70 : /// Tags for constant items added to the GlobalCache. These items are
71 : /// initialized from input file options.
72 1 : using const_global_cache_tags = tmpl::list<::domain::Tags::Domain<dim>>;
73 :
74 : /// Tags for mutable items added to the GlobalCache. These items are
75 : /// initialized from input file options.
76 1 : using mutable_global_cache_tags = tmpl::list<tmpl::conditional_t<
77 : UseControlSystems, ::control_system::Tags::FunctionsOfTimeInitialize,
78 : ::domain::Tags::FunctionsOfTimeInitialize>>;
79 :
80 : /// Tags for simple DataBox items that are initialized from input file options
81 1 : using simple_tags_from_options =
82 : tmpl::list<::domain::Tags::InitialExtents<dim>,
83 : ::domain::Tags::InitialRefinementLevels<dim>,
84 : evolution::dg::Tags::Quadrature>;
85 :
86 : /// Tags for simple DataBox items that are default initialized.
87 1 : using default_initialized_simple_tags = tmpl::list<>;
88 :
89 : /// Tags for items fetched by the DataBox and passed to the apply function
90 1 : using argument_tags =
91 : tmpl::append<const_global_cache_tags, simple_tags_from_options,
92 : tmpl::list<::Parallel::Tags::ArrayIndex<ElementId<dim>>>>;
93 :
94 : /// Tags for items in the DataBox that are mutated by the apply function
95 1 : using return_tags =
96 : tmpl::list<::domain::Tags::Mesh<dim>, ::domain::Tags::Element<dim>,
97 : ::domain::Tags::ElementMap<dim, Frame::Grid>,
98 : ::domain::CoordinateMaps::Tags::CoordinateMap<dim, Frame::Grid,
99 : Frame::Inertial>,
100 : ::domain::Tags::NeighborMesh<dim>>;
101 :
102 : /// Tags for mutable DataBox items that are either default initialized or
103 : /// initialized by the apply function
104 1 : using simple_tags =
105 : tmpl::append<default_initialized_simple_tags, return_tags>;
106 :
107 : /// Tags for immutable DataBox items (compute items or reference items) added
108 : /// to the DataBox.
109 1 : using compute_tags = tmpl::list<
110 : ::domain::Tags::LogicalCoordinates<dim>,
111 : // Compute tags for Frame::Grid quantities
112 : ::domain::Tags::MappedCoordinates<
113 : ::domain::Tags::ElementMap<dim, Frame::Grid>,
114 : ::domain::Tags::Coordinates<dim, Frame::ElementLogical>>,
115 : ::domain::Tags::InverseJacobianCompute<
116 : ::domain::Tags::ElementMap<dim, Frame::Grid>,
117 : ::domain::Tags::Coordinates<dim, Frame::ElementLogical>>,
118 : // Compute tag to retrieve functions of time from global cache.
119 : Parallel::Tags::FromGlobalCache<
120 : tmpl::conditional_t<UseControlSystems,
121 : ::control_system::Tags::FunctionsOfTimeInitialize,
122 : ::domain::Tags::FunctionsOfTimeInitialize>,
123 : Metavariables>,
124 : // Compute tags for Frame::Inertial quantities
125 : ::domain::Tags::CoordinatesMeshVelocityAndJacobiansCompute<
126 : ::domain::CoordinateMaps::Tags::CoordinateMap<dim, Frame::Grid,
127 : Frame::Inertial>>,
128 : ::domain::Tags::GridToInertialInverseJacobian<dim>,
129 : ::domain::Tags::GridToInertialJacobian<dim>,
130 :
131 : ::domain::Tags::InertialFromGridCoordinatesCompute<dim>,
132 : ::domain::Tags::ElementToInertialInverseJacobian<dim>,
133 : ::domain::Tags::DetInvJacobianCompute<dim, Frame::ElementLogical,
134 : Frame::Inertial>,
135 : ::domain::Tags::InertialMeshVelocityCompute<dim>,
136 : evolution::domain::Tags::DivMeshVelocityCompute<dim>,
137 : // Compute tags for other mesh quantities
138 : ::domain::Tags::MinimumGridSpacingCompute<dim, Frame::Inertial>>;
139 :
140 : /// Given the items fetched from a DataBox by the argument_tags, mutate
141 : /// the items in the DataBox corresponding to return_tags
142 1 : static void apply(
143 : const gsl::not_null<Mesh<dim>*> mesh,
144 : const gsl::not_null<Element<dim>*> element,
145 : const gsl::not_null<ElementMap<dim, Frame::Grid>*> element_map,
146 : const gsl::not_null<std::unique_ptr<
147 : ::domain::CoordinateMapBase<Frame::Grid, Frame::Inertial, dim>>*>
148 : grid_to_inertial_map,
149 : const gsl::not_null<DirectionalIdMap<dim, Mesh<dim>>*> neighbor_mesh,
150 : const ::Domain<dim>& domain,
151 : const std::vector<std::array<size_t, dim>>& initial_extents,
152 : const std::vector<std::array<size_t, dim>>& initial_refinement,
153 : const Spectral::Quadrature& i1_quadrature,
154 : const ElementId<dim>& element_id) {
155 : *element = ::domain::create_initial_element(element_id, domain.blocks(),
156 : initial_refinement);
157 : const Spectral::Basis i1_basis{Spectral::Basis::Legendre};
158 : *mesh = ::domain::create_initial_mesh(initial_extents, *element, i1_basis,
159 : i1_quadrature);
160 : if (not ::domain::is_valid_dg_mesh(*mesh, *element)) {
161 : ERROR("Invalid mesh: " << *mesh << " for Element " << *element);
162 : }
163 : const auto& my_block = domain.blocks()[element_id.block_id()];
164 : *element_map = ElementMap<dim, Frame::Grid>{element_id, my_block};
165 :
166 : if (my_block.is_time_dependent()) {
167 : *grid_to_inertial_map =
168 : my_block.moving_mesh_grid_to_inertial_map().get_clone();
169 : } else {
170 : *grid_to_inertial_map =
171 : ::domain::make_coordinate_map_base<Frame::Grid, Frame::Inertial>(
172 : ::domain::CoordinateMaps::Identity<dim>{});
173 : }
174 :
175 : for (const auto& [direction, neighbors] : element->neighbors()) {
176 : for (const auto& neighbor : neighbors) {
177 : const auto& neighbor_block = domain.blocks()[neighbor.block_id()];
178 : if (neighbors.are_conforming()) {
179 : const auto& neighbor_orientation = neighbors.orientation(neighbor);
180 : neighbor_mesh->emplace(
181 : DirectionalId{direction, neighbor},
182 : neighbor_orientation.inverse_map()(::domain::create_initial_mesh(
183 : initial_extents, neighbor_block, neighbor, i1_basis,
184 : i1_quadrature)));
185 : } else if (element->face_types().at(direction) ==
186 : ::domain::FaceType::SingleNonconforming) {
187 : // We do not insert neighbor meshes into neighbor_mesh for a direction
188 : // with domain::FaceType::MultipleNonconforming as this could
189 : // overflow the FixedHashMap size
190 : neighbor_mesh->emplace(
191 : DirectionalId{direction, neighbor},
192 : ::domain::create_initial_mesh(initial_extents, neighbor_block,
193 : neighbor, i1_basis, i1_quadrature));
194 : }
195 : }
196 : }
197 : }
198 : };
199 :
200 : /// \brief Initialize/update items related to coordinate maps after an AMR
201 : /// change
202 : template <size_t Dim>
203 1 : struct ProjectDomain : tt::ConformsTo<amr::protocols::Projector> {
204 0 : using return_tags = tmpl::list<::domain::Tags::ElementMap<Dim, Frame::Grid>,
205 : ::domain::CoordinateMaps::Tags::CoordinateMap<
206 : Dim, Frame::Grid, Frame::Inertial>>;
207 0 : using argument_tags =
208 : tmpl::list<::domain::Tags::Domain<Dim>, ::domain::Tags::Element<Dim>>;
209 :
210 0 : static void apply(
211 : const gsl::not_null<ElementMap<Dim, Frame::Grid>*> /*element_map*/,
212 : const gsl::not_null<std::unique_ptr<
213 : ::domain::CoordinateMapBase<Frame::Grid, Frame::Inertial, Dim>>*>
214 : /*grid_to_inertial_map*/,
215 : const ::Domain<Dim>& /*domain*/, const Element<Dim>& /*element*/,
216 : const std::pair<Mesh<Dim>, Element<Dim>>& /*old_mesh_and_element*/) {
217 : // Do not change anything for p-refinement
218 : }
219 :
220 : template <typename ParentOrChildrenItemsType>
221 0 : static void apply(
222 : const gsl::not_null<ElementMap<Dim, Frame::Grid>*> element_map,
223 : const gsl::not_null<std::unique_ptr<
224 : ::domain::CoordinateMapBase<Frame::Grid, Frame::Inertial, Dim>>*>
225 : grid_to_inertial_map,
226 : const ::Domain<Dim>& domain, const Element<Dim>& element,
227 : const ParentOrChildrenItemsType& /*parent_or_children_items*/) {
228 : const ElementId<Dim>& element_id = element.id();
229 : const auto& my_block = domain.blocks()[element_id.block_id()];
230 : *element_map = ElementMap<Dim, Frame::Grid>{element_id, my_block};
231 : if (my_block.is_time_dependent()) {
232 : *grid_to_inertial_map =
233 : my_block.moving_mesh_grid_to_inertial_map().get_clone();
234 : } else {
235 : *grid_to_inertial_map =
236 : ::domain::make_coordinate_map_base<Frame::Grid, Frame::Inertial>(
237 : ::domain::CoordinateMaps::Identity<Dim>{});
238 : }
239 : }
240 : };
241 : } // namespace evolution::dg::Initialization
|
