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 <optional>
8 : #include <type_traits>
9 : #include <unordered_set>
10 : #include <utility>
11 :
12 : #include "DataStructures/DataBox/DataBox.hpp"
13 : #include "DataStructures/DataBox/MetavariablesTag.hpp"
14 : #include "DataStructures/DataVector.hpp"
15 : #include "DataStructures/Variables.hpp"
16 : #include "Domain/BoundaryConditions/BoundaryCondition.hpp"
17 : #include "Domain/BoundaryConditions/None.hpp"
18 : #include "Domain/BoundaryConditions/Periodic.hpp"
19 : #include "Domain/Creators/Tags/ExternalBoundaryConditions.hpp"
20 : #include "Domain/Domain.hpp"
21 : #include "Domain/Structure/Direction.hpp"
22 : #include "Domain/Structure/Element.hpp"
23 : #include "Domain/Structure/ElementId.hpp"
24 : #include "Domain/Tags.hpp"
25 : #include "Domain/TagsTimeDependent.hpp"
26 : #include "Evolution/BoundaryConditions/Type.hpp"
27 : #include "Evolution/DgSubcell/Tags/CellCenteredFlux.hpp"
28 : #include "Evolution/DgSubcell/Tags/GhostDataForReconstruction.hpp"
29 : #include "Evolution/DgSubcell/Tags/Mesh.hpp"
30 : #include "Evolution/DiscontinuousGalerkin/NormalVectorTags.hpp"
31 : #include "Evolution/Systems/GrMhd/ValenciaDivClean/BoundaryConditions/BoundaryCondition.hpp"
32 : #include "Evolution/Systems/GrMhd/ValenciaDivClean/BoundaryConditions/Factory.hpp"
33 : #include "Evolution/Systems/GrMhd/ValenciaDivClean/FiniteDifference/Reconstructor.hpp"
34 : #include "Evolution/Systems/GrMhd/ValenciaDivClean/System.hpp"
35 : #include "Evolution/VariableFixing/FixToAtmosphere.hpp"
36 : #include "Evolution/VariableFixing/Tags.hpp"
37 : #include "NumericalAlgorithms/Spectral/Mesh.hpp"
38 : #include "PointwiseFunctions/GeneralRelativity/Tags.hpp"
39 : #include "PointwiseFunctions/Hydro/Tags.hpp"
40 : #include "Utilities/CallWithDynamicType.hpp"
41 : #include "Utilities/ErrorHandling/Assert.hpp"
42 : #include "Utilities/Gsl.hpp"
43 : #include "Utilities/PrettyType.hpp"
44 : #include "Utilities/TMPL.hpp"
45 :
46 1 : namespace grmhd::ValenciaDivClean::fd {
47 : /*!
48 : * \brief Computes finite difference ghost data for external boundary
49 : * conditions.
50 : *
51 : * If the element is at the external boundary, computes FD ghost data with a
52 : * given boundary condition and stores it into neighbor data with {direction,
53 : * ElementId::external_boundary_id()} as the mortar_id key.
54 : *
55 : * \note Subcell needs to be enabled for boundary elements. Otherwise this
56 : * function would be never called.
57 : *
58 : */
59 1 : struct BoundaryConditionGhostData {
60 : template <typename DbTagsList>
61 0 : static void apply(gsl::not_null<db::DataBox<DbTagsList>*> box,
62 : const Element<3>& element,
63 : const Reconstructor& reconstructor);
64 :
65 : private:
66 : template <typename FdBoundaryConditionHelper, typename DbTagsList,
67 : typename... FdBoundaryConditionArgsTags>
68 : // A helper function for calling fd_ghost() of BoundaryCondition subclasses
69 0 : static void apply_subcell_boundary_condition_impl(
70 : FdBoundaryConditionHelper& fd_boundary_condition_helper,
71 : const gsl::not_null<db::DataBox<DbTagsList>*>& box,
72 : tmpl::list<FdBoundaryConditionArgsTags...>) {
73 : return fd_boundary_condition_helper(
74 : db::get<FdBoundaryConditionArgsTags>(*box)...);
75 : }
76 : };
77 :
78 : template <typename DbTagsList>
79 0 : void BoundaryConditionGhostData::apply(
80 : const gsl::not_null<db::DataBox<DbTagsList>*> box,
81 : const Element<3>& element, const Reconstructor& reconstructor) {
82 : const auto& external_boundary_condition =
83 : db::get<domain::Tags::ExternalBoundaryConditions<3>>(*box).at(
84 : element.id().block_id());
85 :
86 : // Check if the element is on the external boundary. If not, the caller is
87 : // doing something wrong (e.g. trying to compute FD ghost data with boundary
88 : // conditions at an element which is not on the external boundary).
89 : ASSERT(not element.external_boundaries().empty(),
90 : "The element (ID : " << element.id()
91 : << ") is not on external boundaries");
92 :
93 : const Mesh<3> subcell_mesh =
94 : db::get<evolution::dg::subcell::Tags::Mesh<3>>(*box);
95 :
96 : const size_t ghost_zone_size{reconstructor.ghost_zone_size()};
97 :
98 : // Tags and tags list for FD reconstruction
99 : using RestMassDensity = hydro::Tags::RestMassDensity<DataVector>;
100 : using ElectronFraction = hydro::Tags::ElectronFraction<DataVector>;
101 : using Temperature = hydro::Tags::Temperature<DataVector>;
102 : using LorentzFactorTimesSpatialVelocity =
103 : hydro::Tags::LorentzFactorTimesSpatialVelocity<DataVector, 3>;
104 : using MagneticField = hydro::Tags::MagneticField<DataVector, 3>;
105 : using DivergenceCleaningField =
106 : hydro::Tags::DivergenceCleaningField<DataVector>;
107 :
108 : using prims_for_reconstruction =
109 : tmpl::list<RestMassDensity, ElectronFraction, Temperature,
110 : LorentzFactorTimesSpatialVelocity, MagneticField,
111 : DivergenceCleaningField>;
112 :
113 : size_t num_prims_tensor_components = 0;
114 : tmpl::for_each<prims_for_reconstruction>([&num_prims_tensor_components](
115 : auto tag) {
116 : num_prims_tensor_components += tmpl::type_from<decltype(tag)>::type::size();
117 : });
118 :
119 : using flux_variables =
120 : typename grmhd::ValenciaDivClean::System::flux_variables;
121 : const bool compute_cell_centered_flux =
122 : db::get<
123 : evolution::dg::subcell::Tags::CellCenteredFlux<flux_variables, 3>>(
124 : *box)
125 : .has_value();
126 :
127 : for (const auto& direction : element.external_boundaries()) {
128 : const auto& boundary_condition_at_direction =
129 : *external_boundary_condition.at(direction);
130 :
131 : const size_t num_face_pts{
132 : subcell_mesh.extents().slice_away(direction.dimension()).product()};
133 :
134 : // Allocate a vector to store the computed FD ghost data and assign a
135 : // non-owning Variables on it.
136 : using FluxVars =
137 : Variables<db::wrap_tags_in<::Tags::Flux, flux_variables,
138 : tmpl::size_t<3>, Frame::Inertial>>;
139 : const size_t prims_size =
140 : num_prims_tensor_components * ghost_zone_size * num_face_pts;
141 : const size_t fluxes_size =
142 : (compute_cell_centered_flux ? FluxVars::number_of_independent_components
143 : : 0) *
144 : ghost_zone_size * num_face_pts;
145 :
146 : auto& all_ghost_data = db::get_mutable_reference<
147 : evolution::dg::subcell::Tags::GhostDataForReconstruction<3>>(box);
148 : // Put the computed ghost data into neighbor data with {direction,
149 : // ElementId::external_boundary_id()} as the mortar_id key
150 : const DirectionalId<3> mortar_id{direction,
151 : ElementId<3>::external_boundary_id()};
152 :
153 : all_ghost_data[mortar_id] = evolution::dg::subcell::GhostData{1};
154 : DataVector& boundary_ghost_data =
155 : all_ghost_data.at(mortar_id).neighbor_ghost_data_for_reconstruction();
156 : boundary_ghost_data.destructive_resize(prims_size + fluxes_size);
157 : Variables<prims_for_reconstruction> ghost_data_vars{
158 : boundary_ghost_data.data(), prims_size};
159 : std::optional<FluxVars> cell_centered_ghost_fluxes{};
160 : if (compute_cell_centered_flux) {
161 : cell_centered_ghost_fluxes = FluxVars{};
162 : cell_centered_ghost_fluxes.value().set_data_ref(
163 : std::next(boundary_ghost_data.data(),
164 : static_cast<std::ptrdiff_t>(prims_size)),
165 : fluxes_size);
166 : }
167 :
168 : // We don't need to care about boundary ghost data when using the periodic
169 : // condition, so exclude it from the type list
170 : using factory_classes =
171 : typename std::decay_t<decltype(db::get<Parallel::Tags::Metavariables>(
172 : *box))>::factory_creation::factory_classes;
173 : using derived_boundary_conditions_for_subcell = tmpl::remove_if<
174 : tmpl::at<factory_classes, typename System::boundary_conditions_base>,
175 : tmpl::or_<
176 : std::is_base_of<domain::BoundaryConditions::MarkAsPeriodic,
177 : tmpl::_1>,
178 : std::is_base_of<domain::BoundaryConditions::MarkAsNone, tmpl::_1>>>;
179 :
180 : // Now apply subcell boundary conditions
181 : call_with_dynamic_type<void, derived_boundary_conditions_for_subcell>(
182 : &boundary_condition_at_direction,
183 : [&box, &cell_centered_ghost_fluxes, &direction,
184 : &ghost_data_vars](const auto* boundary_condition) {
185 : using BoundaryCondition = std::decay_t<decltype(*boundary_condition)>;
186 : using bcondition_interior_evolved_vars_tags =
187 : typename BoundaryCondition::fd_interior_evolved_variables_tags;
188 : using bcondition_interior_temporary_tags =
189 : typename BoundaryCondition::fd_interior_temporary_tags;
190 : using bcondition_interior_primitive_vars_tags =
191 : typename BoundaryCondition::fd_interior_primitive_variables_tags;
192 : using bcondition_gridless_tags =
193 : typename BoundaryCondition::fd_gridless_tags;
194 :
195 : using bcondition_interior_tags =
196 : tmpl::append<bcondition_interior_evolved_vars_tags,
197 : bcondition_interior_temporary_tags,
198 : bcondition_interior_primitive_vars_tags,
199 : bcondition_gridless_tags>;
200 :
201 : if constexpr (BoundaryCondition::bc_type ==
202 : evolution::BoundaryConditions::Type::Ghost) {
203 : const auto apply_fd_ghost =
204 : [&boundary_condition, &cell_centered_ghost_fluxes, &direction,
205 : &ghost_data_vars](const auto&... boundary_ghost_data_args) {
206 : (*boundary_condition)
207 : .fd_ghost(
208 : make_not_null(&get<RestMassDensity>(ghost_data_vars)),
209 : make_not_null(
210 : &get<ElectronFraction>(ghost_data_vars)),
211 : make_not_null(&get<Temperature>(ghost_data_vars)),
212 : make_not_null(&get<LorentzFactorTimesSpatialVelocity>(
213 : ghost_data_vars)),
214 : make_not_null(&get<MagneticField>(ghost_data_vars)),
215 : make_not_null(
216 : &get<DivergenceCleaningField>(ghost_data_vars)),
217 : make_not_null(&cell_centered_ghost_fluxes), direction,
218 : boundary_ghost_data_args...);
219 : };
220 : apply_subcell_boundary_condition_impl(apply_fd_ghost, box,
221 : bcondition_interior_tags{});
222 : } else if constexpr (BoundaryCondition::bc_type ==
223 : evolution::BoundaryConditions::Type::
224 : DemandOutgoingCharSpeeds) {
225 : // This boundary condition only checks if all the characteristic
226 : // speeds are directed outward.
227 : const auto& volume_mesh_velocity =
228 : db::get<domain::Tags::MeshVelocity<3, Frame::Inertial>>(*box);
229 : if (volume_mesh_velocity.has_value()) {
230 : ERROR("Subcell currently does not support moving mesh");
231 : }
232 :
233 : std::optional<tnsr::I<DataVector, 3>> face_mesh_velocity{};
234 :
235 : const auto& normal_covector_and_magnitude =
236 : db::get<evolution::dg::Tags::NormalCovectorAndMagnitude<3>>(
237 : *box);
238 : const auto outward_directed_normal_covector =
239 : get<evolution::dg::Tags::NormalCovector<3>>(
240 : normal_covector_and_magnitude.at(direction).value());
241 :
242 : const auto apply_fd_demand_outgoing_char_speeds =
243 : [&boundary_condition, &cell_centered_ghost_fluxes, &direction,
244 : &face_mesh_velocity, &ghost_data_vars,
245 : &outward_directed_normal_covector](
246 : const auto&... boundary_ghost_data_args) {
247 : return (*boundary_condition)
248 : .fd_demand_outgoing_char_speeds(
249 : make_not_null(&get<RestMassDensity>(ghost_data_vars)),
250 : make_not_null(
251 : &get<ElectronFraction>(ghost_data_vars)),
252 : make_not_null(&get<Temperature>(ghost_data_vars)),
253 : make_not_null(&get<LorentzFactorTimesSpatialVelocity>(
254 : ghost_data_vars)),
255 : make_not_null(&get<MagneticField>(ghost_data_vars)),
256 : make_not_null(
257 : &get<DivergenceCleaningField>(ghost_data_vars)),
258 : make_not_null(&cell_centered_ghost_fluxes), direction,
259 : face_mesh_velocity, outward_directed_normal_covector,
260 : boundary_ghost_data_args...);
261 : };
262 : apply_subcell_boundary_condition_impl(
263 : apply_fd_demand_outgoing_char_speeds, box,
264 : bcondition_interior_tags{});
265 :
266 : return;
267 : } else {
268 : ERROR("Unsupported boundary condition "
269 : << pretty_type::short_name<BoundaryCondition>()
270 : << " when using finite-difference");
271 : }
272 : });
273 : if (dynamic_cast<const BoundaryConditions::DirichletAnalytic*>(
274 : &boundary_condition_at_direction) != nullptr and
275 : reconstructor.reconstruct_rho_times_temperature()) {
276 : // If we reconstruct rho*T we end up having to divide by rho to compute
277 : // T. In some cases, like when evolving a TOV star with an analytic
278 : // boundary condition, the boundary condition sets rho=0. While
279 : // unphysical in general, this is how we have implemented the
280 : // solutions. We deal with this by applying our atmosphere treatment to
281 : // the reconstructed state.
282 :
283 : const auto& atmosphere_fixer =
284 : db::get<::Tags::VariableFixer<VariableFixing::FixToAtmosphere<3>>>(
285 : *box);
286 : const auto& equation_of_state =
287 : db::get<hydro::Tags::GrmhdEquationOfState>(*box);
288 :
289 : Variables<tmpl::list<gr::Tags::SpatialMetric<DataVector, 3>,
290 : hydro::Tags::SpatialVelocity<DataVector, 3>,
291 : hydro::Tags::LorentzFactor<DataVector>,
292 : hydro::Tags::SpecificInternalEnergy<DataVector>,
293 : hydro::Tags::Pressure<DataVector>>>
294 : temp_hydro_vars{ghost_data_vars.number_of_grid_points()};
295 :
296 : tnsr::ii<DataVector, 3, Frame::Inertial>& spatial_metric =
297 : get<gr::Tags::SpatialMetric<DataVector, 3>>(temp_hydro_vars);
298 : for (size_t i = 0; i < 3; ++i) {
299 : for (size_t j = i; j < 3; ++j) {
300 : spatial_metric.get(i, j) = (i == j ? 1.0 : 0.0);
301 : }
302 : }
303 :
304 : const auto& lorentz_factor_times_spatial_velocity =
305 : get<hydro::Tags::LorentzFactorTimesSpatialVelocity<DataVector, 3>>(
306 : ghost_data_vars);
307 : auto& lorentz_factor =
308 : get<hydro::Tags::LorentzFactor<DataVector>>(temp_hydro_vars);
309 : get(lorentz_factor) = 0.0;
310 : for (size_t i = 0; i < 3; ++i) {
311 : get(lorentz_factor) +=
312 : spatial_metric.get(i, i) *
313 : square(lorentz_factor_times_spatial_velocity.get(i));
314 : for (size_t j = i + 1; j < 3; ++j) {
315 : get(lorentz_factor) += 2.0 * spatial_metric.get(i, j) *
316 : lorentz_factor_times_spatial_velocity.get(i) *
317 : lorentz_factor_times_spatial_velocity.get(j);
318 : }
319 : }
320 : get(lorentz_factor) = sqrt(1.0 + get(lorentz_factor));
321 : auto& spatial_velocity =
322 : get<hydro::Tags::SpatialVelocity<DataVector, 3>>(temp_hydro_vars) =
323 : lorentz_factor_times_spatial_velocity;
324 : for (size_t i = 0; i < 3; ++i) {
325 : spatial_velocity.get(i) /= get(lorentz_factor);
326 : }
327 :
328 : get<hydro::Tags::SpecificInternalEnergy<DataVector>>(temp_hydro_vars) =
329 : equation_of_state
330 : .specific_internal_energy_from_density_and_temperature(
331 : get<RestMassDensity>(ghost_data_vars),
332 : get<Temperature>(ghost_data_vars),
333 : get<ElectronFraction>(ghost_data_vars));
334 : get<hydro::Tags::Pressure<DataVector>>(temp_hydro_vars) =
335 : equation_of_state.pressure_from_density_and_temperature(
336 : get<RestMassDensity>(ghost_data_vars),
337 : get<Temperature>(ghost_data_vars),
338 : get<ElectronFraction>(ghost_data_vars));
339 :
340 : atmosphere_fixer(
341 : make_not_null(&get<RestMassDensity>(ghost_data_vars)),
342 : make_not_null(&get<hydro::Tags::SpecificInternalEnergy<DataVector>>(
343 : temp_hydro_vars)),
344 : make_not_null(&get<hydro::Tags::SpatialVelocity<DataVector, 3>>(
345 : temp_hydro_vars)),
346 : make_not_null(
347 : &get<hydro::Tags::LorentzFactor<DataVector>>(temp_hydro_vars)),
348 : make_not_null(
349 : &get<hydro::Tags::Pressure<DataVector>>(temp_hydro_vars)),
350 : make_not_null(&get<Temperature>(ghost_data_vars)),
351 : get<ElectronFraction>(ghost_data_vars), spatial_metric,
352 : equation_of_state);
353 : }
354 : }
355 : }
356 : } // namespace grmhd::ValenciaDivClean::fd
|
