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