Line data Source code
1 0 : // Distributed under the MIT License.
2 : // See LICENSE.txt for details.
3 :
4 : #pragma once
5 :
6 : #include <cstddef>
7 : #include <memory>
8 : #include <string>
9 : #include <unordered_map>
10 :
11 : #include "DataStructures/DataBox/Protocols/Mutator.hpp"
12 : #include "DataStructures/DataVector.hpp"
13 : #include "DataStructures/Tensor/Tensor.hpp"
14 : #include "DataStructures/Variables.hpp"
15 : #include "Domain/Block.hpp"
16 : #include "Domain/CoordinateMaps/CoordinateMap.hpp"
17 : #include "Domain/CoordinateMaps/Tags.hpp"
18 : #include "Domain/Creators/Tags/Domain.hpp"
19 : #include "Domain/Domain.hpp"
20 : #include "Domain/ElementMap.hpp"
21 : #include "Domain/FunctionsOfTime/FunctionOfTime.hpp"
22 : #include "Domain/FunctionsOfTime/Tags.hpp"
23 : #include "Domain/Structure/Element.hpp"
24 : #include "Domain/Structure/ElementId.hpp"
25 : #include "Domain/Tags.hpp"
26 : #include "Evolution/DgSubcell/Tags/Coordinates.hpp"
27 : #include "Evolution/DgSubcell/Tags/DidRollback.hpp"
28 : #include "Evolution/DgSubcell/Tags/Inactive.hpp"
29 : #include "Evolution/DgSubcell/Tags/Mesh.hpp"
30 : #include "Evolution/DgSubcell/Tags/OnSubcellFaces.hpp"
31 : #include "Evolution/Initialization/InitialData.hpp"
32 : #include "NumericalAlgorithms/Spectral/LogicalCoordinates.hpp"
33 : #include "NumericalAlgorithms/Spectral/Mesh.hpp"
34 : #include "NumericalAlgorithms/Spectral/Quadrature.hpp"
35 : #include "PointwiseFunctions/AnalyticData/Tags.hpp"
36 : #include "PointwiseFunctions/InitialDataUtilities/InitialData.hpp"
37 : #include "PointwiseFunctions/InitialDataUtilities/Tags/InitialData.hpp"
38 : #include "Time/Tags/Time.hpp"
39 : #include "Utilities/CallWithDynamicType.hpp"
40 : #include "Utilities/ErrorHandling/Assert.hpp"
41 : #include "Utilities/Gsl.hpp"
42 : #include "Utilities/ProtocolHelpers.hpp"
43 : #include "Utilities/TMPL.hpp"
44 :
45 : namespace evolution::dg::subcell {
46 :
47 : /*!
48 : * \brief Allocate or assign background general relativity quantities on
49 : * cell-centered and face-centered FD grid points, for evolution systems run on
50 : * a curved spacetime without solving Einstein equations (e.g. ValenciaDivclean,
51 : * ForceFree),
52 : *
53 : * \warning This mutator assumes that the GR analytic data or solution
54 : * specifying background spacetime metric is time-independent.
55 : *
56 : */
57 : template <typename System, typename Metavariables, bool ComputeOnlyOnRollback>
58 1 : struct BackgroundGrVars : tt::ConformsTo<db::protocols::Mutator> {
59 0 : static constexpr size_t volume_dim = System::volume_dim;
60 :
61 0 : using gr_vars_tag = typename System::spacetime_variables_tag;
62 0 : using inactive_gr_vars_tag =
63 : evolution::dg::subcell::Tags::Inactive<gr_vars_tag>;
64 0 : using subcell_faces_gr_tag = evolution::dg::subcell::Tags::OnSubcellFaces<
65 : typename System::flux_spacetime_variables_tag, volume_dim>;
66 :
67 0 : using GrVars = typename gr_vars_tag::type;
68 0 : using InactiveGrVars = typename inactive_gr_vars_tag::type;
69 0 : using SubcellFaceGrVars = typename subcell_faces_gr_tag::type;
70 :
71 0 : using argument_tags = tmpl::list<
72 : ::Tags::Time, domain::Tags::FunctionsOfTime,
73 : domain::Tags::Domain<volume_dim>, domain::Tags::Element<volume_dim>,
74 : domain::Tags::ElementMap<volume_dim, Frame::Grid>,
75 : domain::CoordinateMaps::Tags::CoordinateMap<volume_dim, Frame::Grid,
76 : Frame::Inertial>,
77 : evolution::dg::subcell::Tags::Mesh<volume_dim>,
78 : evolution::dg::subcell::Tags::Coordinates<3, Frame::Inertial>,
79 : subcell::Tags::DidRollback, evolution::initial_data::Tags::InitialData>;
80 :
81 0 : using return_tags =
82 : tmpl::list<gr_vars_tag, inactive_gr_vars_tag, subcell_faces_gr_tag>;
83 :
84 : template <typename T>
85 0 : static void apply(
86 : const gsl::not_null<GrVars*> active_gr_vars,
87 : const gsl::not_null<InactiveGrVars*> inactive_gr_vars,
88 : const gsl::not_null<SubcellFaceGrVars*> subcell_face_gr_vars,
89 : const double time,
90 : const std::unordered_map<
91 : std::string,
92 : std::unique_ptr<::domain::FunctionsOfTime::FunctionOfTime>>&
93 : functions_of_time,
94 : const Domain<volume_dim>& domain, const Element<volume_dim>& element,
95 : const ElementMap<volume_dim, Frame::Grid>& logical_to_grid_map,
96 : const domain::CoordinateMapBase<Frame::Grid, Frame::Inertial, volume_dim>&
97 : grid_to_inertial_map,
98 : const Mesh<volume_dim>& subcell_mesh,
99 : const tnsr::I<DataVector, volume_dim, Frame::Inertial>&
100 : subcell_inertial_coords,
101 : const bool did_rollback, const T& solution_or_data) {
102 : const size_t num_subcell_pts = subcell_mesh.number_of_grid_points();
103 :
104 : if (gsl::at(*subcell_face_gr_vars, 0).number_of_grid_points() != 0) {
105 : // Evolution phase
106 :
107 : // Check if the mesh is actually moving i.e. block coordinate map is
108 : // time-dependent. If not, we can skip the evaluation of GR variables
109 : // since they may stay with their values assigned at the initialization
110 : // phase.
111 : const auto& element_id = element.id();
112 : const size_t block_id = element_id.block_id();
113 : const Block<volume_dim>& block = domain.blocks()[block_id];
114 :
115 : if (block.is_time_dependent()) {
116 : if (did_rollback or not ComputeOnlyOnRollback) {
117 : if (did_rollback) {
118 : // Right after rollback, subcell GR vars are stored in the
119 : // `inactive` one.
120 : ASSERT(inactive_gr_vars->number_of_grid_points() == num_subcell_pts,
121 : "The size of subcell GR variables ("
122 : << inactive_gr_vars->number_of_grid_points()
123 : << ") is not equal to the number of FD grid points ("
124 : << subcell_mesh.number_of_grid_points() << ").");
125 :
126 : cell_centered_impl(inactive_gr_vars, time, subcell_inertial_coords,
127 : solution_or_data);
128 :
129 : } else {
130 : // In this case the element didn't rollback but started from FD.
131 : // Therefore subcell GR vars are in the `active` one.
132 : ASSERT(active_gr_vars->number_of_grid_points() == num_subcell_pts,
133 : "The size of subcell GR variables ("
134 : << active_gr_vars->number_of_grid_points()
135 : << ") is not equal to the number of FD grid points ("
136 : << subcell_mesh.number_of_grid_points() << ").");
137 :
138 : cell_centered_impl(active_gr_vars, time, subcell_inertial_coords,
139 : solution_or_data);
140 : }
141 : if constexpr (not std::is_same_v<
142 : typename SubcellFaceGrVars::value_type::tags_list,
143 : tmpl::list<>>){
144 : face_centered_impl(subcell_face_gr_vars, time, functions_of_time,
145 : logical_to_grid_map, grid_to_inertial_map,
146 : subcell_mesh, solution_or_data);
147 : }
148 : }
149 : }
150 :
151 : } else {
152 : // Initialization phase
153 : (*inactive_gr_vars).initialize(num_subcell_pts);
154 :
155 : ASSERT(Mesh<volume_dim>(subcell_mesh.extents(0), subcell_mesh.basis(0),
156 : subcell_mesh.quadrature(0)) == subcell_mesh,
157 : "The subcell mesh must have isotropic basis, quadrature. and "
158 : "extents but got "
159 : << subcell_mesh);
160 : if constexpr (not std::is_same_v<
161 : typename SubcellFaceGrVars::value_type::tags_list,
162 : tmpl::list<>>){
163 : const size_t num_face_centered_mesh_grid_pts =
164 : (subcell_mesh.extents(0) + 1) * subcell_mesh.extents(1) *
165 : subcell_mesh.extents(2);
166 : for (size_t d = 0; d < volume_dim; ++d) {
167 : gsl::at(*subcell_face_gr_vars, d)
168 : .initialize(num_face_centered_mesh_grid_pts);
169 : }
170 : face_centered_impl(subcell_face_gr_vars, time, functions_of_time,
171 : logical_to_grid_map, grid_to_inertial_map,
172 : subcell_mesh, solution_or_data);
173 : }
174 :
175 : cell_centered_impl(inactive_gr_vars, time, subcell_inertial_coords,
176 : solution_or_data);
177 : }
178 : }
179 :
180 : private:
181 : template <typename Vars, typename T>
182 0 : static void cell_centered_impl(
183 : const gsl::not_null<Vars*> background_gr_vars, const double time,
184 : const tnsr::I<DataVector, volume_dim, Frame::Inertial>& inertial_coords,
185 : const T& solution_or_data) {
186 : GrVars temp{background_gr_vars->data(), background_gr_vars->size()};
187 :
188 : using derived_classes =
189 : tmpl::at<typename Metavariables::factory_creation::factory_classes,
190 : evolution::initial_data::InitialData>;
191 : call_with_dynamic_type<void, derived_classes>(
192 : &solution_or_data, [&temp, &inertial_coords,
193 : &time](const auto* const solution_or_data_ptr) {
194 : temp.assign_subset(evolution::Initialization::initial_data(
195 : *solution_or_data_ptr, inertial_coords, time,
196 : typename GrVars::tags_list{}));
197 : });
198 : }
199 :
200 : template <typename T>
201 0 : static void face_centered_impl(
202 : const gsl::not_null<SubcellFaceGrVars*> face_centered_gr_vars,
203 : const double time,
204 : const std::unordered_map<
205 : std::string,
206 : std::unique_ptr<::domain::FunctionsOfTime::FunctionOfTime>>&
207 : functions_of_time,
208 : const ElementMap<volume_dim, Frame::Grid>& logical_to_grid_map,
209 : const domain::CoordinateMapBase<Frame::Grid, Frame::Inertial, volume_dim>&
210 : grid_to_inertial_map,
211 : const Mesh<volume_dim>& subcell_mesh, const T& solution_or_data) {
212 : ASSERT(Mesh<volume_dim>(subcell_mesh.extents(0), subcell_mesh.basis(0),
213 : subcell_mesh.quadrature(0)) == subcell_mesh,
214 : "The subcell mesh must have isotropic basis, quadrature. and "
215 : "extents but got "
216 : << subcell_mesh);
217 :
218 : for (size_t dim = 0; dim < volume_dim; ++dim) {
219 : const auto basis = make_array<volume_dim>(subcell_mesh.basis(0));
220 : auto quadrature = make_array<volume_dim>(subcell_mesh.quadrature(0));
221 : auto extents = make_array<volume_dim>(subcell_mesh.extents(0));
222 :
223 : gsl::at(extents, dim) = subcell_mesh.extents(0) + 1;
224 : gsl::at(quadrature, dim) = Spectral::Quadrature::FaceCentered;
225 :
226 : const Mesh<volume_dim> face_centered_mesh{extents, basis, quadrature};
227 : const auto face_centered_logical_coords =
228 : logical_coordinates(face_centered_mesh);
229 : const auto face_centered_inertial_coords = grid_to_inertial_map(
230 : logical_to_grid_map(face_centered_logical_coords), time,
231 : functions_of_time);
232 :
233 : using derived_classes =
234 : tmpl::at<typename Metavariables::factory_creation::factory_classes,
235 : evolution::initial_data::InitialData>;
236 : call_with_dynamic_type<void, derived_classes>(
237 : &solution_or_data,
238 : [&face_centered_gr_vars, &face_centered_inertial_coords, &dim,
239 : &time](const auto* const solution_or_data_ptr) {
240 : gsl::at(*face_centered_gr_vars, dim)
241 : .assign_subset(evolution::Initialization::initial_data(
242 : *solution_or_data_ptr, face_centered_inertial_coords, time,
243 : typename SubcellFaceGrVars::value_type::tags_list{}));
244 : });
245 : }
246 : }
247 : };
248 :
249 : } // namespace evolution::dg::subcell
|
