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/GhostDataForReconstruction.hpp"
28 : #include "Evolution/DgSubcell/Tags/Mesh.hpp"
29 : #include "Evolution/DiscontinuousGalerkin/NormalVectorTags.hpp"
30 : #include "Evolution/Systems/GrMhd/GhValenciaDivClean/BoundaryConditions/BoundaryCondition.hpp"
31 : #include "Evolution/Systems/GrMhd/GhValenciaDivClean/BoundaryConditions/DirichletAnalytic.hpp"
32 : #include "Evolution/Systems/GrMhd/GhValenciaDivClean/BoundaryConditions/Factory.hpp"
33 : #include "Evolution/Systems/GrMhd/GhValenciaDivClean/FiniteDifference/Reconstructor.hpp"
34 : #include "Evolution/Systems/GrMhd/GhValenciaDivClean/System.hpp"
35 : #include "Evolution/Systems/GrMhd/GhValenciaDivClean/Tags.hpp"
36 : #include "Evolution/VariableFixing/FixToAtmosphere.hpp"
37 : #include "Evolution/VariableFixing/Tags.hpp"
38 : #include "NumericalAlgorithms/Spectral/Mesh.hpp"
39 : #include "PointwiseFunctions/GeneralRelativity/Tags.hpp"
40 : #include "PointwiseFunctions/Hydro/Tags.hpp"
41 : #include "Utilities/CallWithDynamicType.hpp"
42 : #include "Utilities/ErrorHandling/Assert.hpp"
43 : #include "Utilities/Gsl.hpp"
44 : #include "Utilities/PrettyType.hpp"
45 : #include "Utilities/TMPL.hpp"
46 :
47 1 : namespace grmhd::GhValenciaDivClean::fd {
48 : /*!
49 : * \brief Computes finite difference ghost data for external boundary
50 : * conditions.
51 : *
52 : * If the element is at the external boundary, computes FD ghost data with a
53 : * given boundary condition and stores it into neighbor data with {direction,
54 : * ElementId::external_boundary_id()} as the mortar_id key.
55 : *
56 : * \note Subcell needs to be enabled for boundary elements. Otherwise this
57 : * function would be never called.
58 : *
59 : */
60 : template <typename System>
61 1 : struct BoundaryConditionGhostData {
62 : template <typename DbTagsList>
63 0 : static void apply(
64 : gsl::not_null<db::DataBox<DbTagsList>*> box, const Element<3>& element,
65 : const Reconstructor<System>& reconstructor);
66 :
67 : private:
68 : template <typename FdBoundaryConditionHelper, typename DbTagsList,
69 : typename... FdBoundaryConditionArgsTags>
70 : // A helper function for calling fd_ghost() of BoundaryCondition subclasses
71 0 : static void apply_subcell_boundary_condition_impl(
72 : FdBoundaryConditionHelper& fd_boundary_condition_helper,
73 : const gsl::not_null<db::DataBox<DbTagsList>*>& box,
74 : tmpl::list<FdBoundaryConditionArgsTags...>) {
75 : return fd_boundary_condition_helper(
76 : db::get<FdBoundaryConditionArgsTags>(*box)...);
77 : }
78 : };
79 :
80 : template <typename System>
81 : template <typename DbTagsList>
82 0 : void BoundaryConditionGhostData<System>::apply(
83 : const gsl::not_null<db::DataBox<DbTagsList>*> box,
84 : const Element<3>& element,
85 : const Reconstructor<System>& reconstructor) {
86 : const auto& external_boundary_condition =
87 : db::get<domain::Tags::ExternalBoundaryConditions<3>>(*box).at(
88 : element.id().block_id());
89 :
90 : // Check if the element is on the external boundary. If not, the caller is
91 : // doing something wrong (e.g. trying to compute FD ghost data with boundary
92 : // conditions at an element which is not on the external boundary).
93 : ASSERT(not element.external_boundaries().empty(),
94 : "The element (ID : " << element.id()
95 : << ") is not on external boundaries");
96 :
97 : const Mesh<3> subcell_mesh =
98 : db::get<evolution::dg::subcell::Tags::Mesh<3>>(*box);
99 :
100 : const size_t ghost_zone_size{reconstructor.ghost_zone_size()};
101 :
102 : // Tags and tags list for FD reconstruction
103 : using RestMassDensity = hydro::Tags::RestMassDensity<DataVector>;
104 : using ElectronFraction = hydro::Tags::ElectronFraction<DataVector>;
105 : using Temperature = hydro::Tags::Temperature<DataVector>;
106 : using LorentzFactorTimesSpatialVelocity =
107 : hydro::Tags::LorentzFactorTimesSpatialVelocity<DataVector, 3>;
108 : using MagneticField = hydro::Tags::MagneticField<DataVector, 3>;
109 : using DivergenceCleaningField =
110 : hydro::Tags::DivergenceCleaningField<DataVector>;
111 : using SpacetimeMetric = gr::Tags::SpacetimeMetric<DataVector, 3>;
112 : using Pi = gh::Tags::Pi<DataVector, 3>;
113 : using Phi = gh::Tags::Phi<DataVector, 3>;
114 :
115 : using reconstruction_tags = GhValenciaDivClean::Tags::
116 : primitive_grmhd_and_spacetime_reconstruction_tags;
117 : using NeighborVariables = Variables<reconstruction_tags>;
118 : constexpr size_t number_of_tensor_components =
119 : NeighborVariables::number_of_independent_components;
120 :
121 : for (const auto& direction : element.external_boundaries()) {
122 : const auto& boundary_condition_at_direction =
123 : *external_boundary_condition.at(direction);
124 :
125 : const size_t num_face_pts{
126 : subcell_mesh.extents().slice_away(direction.dimension()).product()};
127 :
128 : // Allocate a vector to store the computed FD ghost data and assign a
129 : // non-owning Variables on it.
130 : auto& all_ghost_data = db::get_mutable_reference<
131 : evolution::dg::subcell::Tags::GhostDataForReconstruction<3>>(box);
132 : // Put the computed ghost data into neighbor data with {direction,
133 : // ElementId::external_boundary_id()} as the mortar_id key
134 : const DirectionalId<3> mortar_id{direction,
135 : ElementId<3>::external_boundary_id()};
136 : all_ghost_data[mortar_id] = evolution::dg::subcell::GhostData{1};
137 : DataVector& boundary_ghost_data =
138 : all_ghost_data.at(mortar_id).neighbor_ghost_data_for_reconstruction();
139 : boundary_ghost_data.destructive_resize(number_of_tensor_components *
140 : ghost_zone_size * num_face_pts);
141 : Variables<reconstruction_tags> ghost_data_vars{boundary_ghost_data.data(),
142 : boundary_ghost_data.size()};
143 :
144 : // We don't need to care about boundary ghost data when using the periodic
145 : // condition, so exclude it from the type list
146 : using factory_classes =
147 : typename std::decay_t<decltype(db::get<Parallel::Tags::Metavariables>(
148 : *box))>::factory_creation::factory_classes;
149 : using derived_boundary_conditions_for_subcell = tmpl::remove_if<
150 : tmpl::at<factory_classes, typename System::
151 : boundary_conditions_base>,
152 : tmpl::or_<
153 : std::is_base_of<domain::BoundaryConditions::MarkAsPeriodic,
154 : tmpl::_1>,
155 : std::is_base_of<domain::BoundaryConditions::MarkAsNone, tmpl::_1>>>;
156 :
157 : // Now apply subcell boundary conditions
158 : call_with_dynamic_type<void, derived_boundary_conditions_for_subcell>(
159 : &boundary_condition_at_direction,
160 : [&box, &direction, &ghost_data_vars](const auto* boundary_condition) {
161 : using BoundaryCondition = std::decay_t<decltype(*boundary_condition)>;
162 : using bcondition_interior_evolved_vars_tags =
163 : typename BoundaryCondition::fd_interior_evolved_variables_tags;
164 : using bcondition_interior_temporary_tags =
165 : typename BoundaryCondition::fd_interior_temporary_tags;
166 : using bcondition_interior_primitive_vars_tags =
167 : typename BoundaryCondition::fd_interior_primitive_variables_tags;
168 : using bcondition_gridless_tags =
169 : typename BoundaryCondition::fd_gridless_tags;
170 :
171 : using bcondition_interior_tags =
172 : tmpl::append<bcondition_interior_evolved_vars_tags,
173 : bcondition_interior_temporary_tags,
174 : bcondition_interior_primitive_vars_tags,
175 : bcondition_gridless_tags>;
176 :
177 : if constexpr (BoundaryCondition::bc_type ==
178 : evolution::BoundaryConditions::Type::Ghost or
179 : BoundaryCondition::bc_type ==
180 : evolution::BoundaryConditions::Type::
181 : GhostAndTimeDerivative) {
182 : const auto apply_fd_ghost =
183 : [&boundary_condition, &direction,
184 : &ghost_data_vars](const auto&... boundary_ghost_data_args) {
185 : (*boundary_condition)
186 : .fd_ghost(
187 : make_not_null(&get<SpacetimeMetric>(ghost_data_vars)),
188 : make_not_null(&get<Pi>(ghost_data_vars)),
189 : make_not_null(&get<Phi>(ghost_data_vars)),
190 : make_not_null(&get<RestMassDensity>(ghost_data_vars)),
191 : make_not_null(
192 : &get<ElectronFraction>(ghost_data_vars)),
193 : make_not_null(&get<Temperature>(ghost_data_vars)),
194 : make_not_null(&get<LorentzFactorTimesSpatialVelocity>(
195 : ghost_data_vars)),
196 : make_not_null(&get<MagneticField>(ghost_data_vars)),
197 : make_not_null(
198 : &get<DivergenceCleaningField>(ghost_data_vars)),
199 : direction, boundary_ghost_data_args...);
200 : };
201 : apply_subcell_boundary_condition_impl(apply_fd_ghost, box,
202 : bcondition_interior_tags{});
203 : } else {
204 : ERROR("Unsupported boundary condition "
205 : << pretty_type::short_name<BoundaryCondition>()
206 : << " when using finite-difference");
207 : }
208 : });
209 : if (dynamic_cast<const BoundaryConditions::DirichletAnalytic<System>*>(
210 : &boundary_condition_at_direction) != nullptr and
211 : reconstructor.reconstruct_rho_times_temperature()) {
212 : // If we reconstruct rho*T we end up having to divide by rho to compute
213 : // T. In some cases, like when evolving a TOV star with an analytic
214 : // boundary condition, the boundary condition sets rho=0. While
215 : // unphysical in general, this is how we have implemented the
216 : // solutions. We deal with this by applying our atmosphere treatment to
217 : // the reconstructed stated.
218 :
219 : const auto& atmosphere_fixer =
220 : db::get<::Tags::VariableFixer<VariableFixing::FixToAtmosphere<3>>>(
221 : *box);
222 : const auto& equation_of_state =
223 : db::get<hydro::Tags::GrmhdEquationOfState>(*box);
224 :
225 : tnsr::ii<DataVector, 3, Frame::Inertial> spatial_metric{};
226 : for (size_t i = 0; i < 3; ++i) {
227 : for (size_t j = i; j < 3; ++j) {
228 : spatial_metric.get(i, j).set_data_ref(
229 : &get<SpacetimeMetric>(ghost_data_vars).get(i + 1, j + 1));
230 : }
231 : }
232 :
233 : Variables<tmpl::list<hydro::Tags::SpatialVelocity<DataVector, 3>,
234 : hydro::Tags::LorentzFactor<DataVector>,
235 : hydro::Tags::SpecificInternalEnergy<DataVector>,
236 : hydro::Tags::Pressure<DataVector>>>
237 : temp_hydro_vars{ghost_data_vars.number_of_grid_points()};
238 :
239 : const auto& lorentz_factor_times_spatial_velocity =
240 : get<hydro::Tags::LorentzFactorTimesSpatialVelocity<DataVector, 3>>(
241 : ghost_data_vars);
242 : auto& lorentz_factor =
243 : get<hydro::Tags::LorentzFactor<DataVector>>(temp_hydro_vars);
244 : get(lorentz_factor) = 0.0;
245 : for (size_t i = 0; i < 3; ++i) {
246 : get(lorentz_factor) +=
247 : spatial_metric.get(i, i) *
248 : square(lorentz_factor_times_spatial_velocity.get(i));
249 : for (size_t j = i + 1; j < 3; ++j) {
250 : get(lorentz_factor) += 2.0 * spatial_metric.get(i, j) *
251 : lorentz_factor_times_spatial_velocity.get(i) *
252 : lorentz_factor_times_spatial_velocity.get(j);
253 : }
254 : }
255 : get(lorentz_factor) = sqrt(1.0 + get(lorentz_factor));
256 : auto& spatial_velocity =
257 : get<hydro::Tags::SpatialVelocity<DataVector, 3>>(temp_hydro_vars) =
258 : lorentz_factor_times_spatial_velocity;
259 : for (size_t i = 0; i < 3; ++i) {
260 : spatial_velocity.get(i) /= get(lorentz_factor);
261 : }
262 :
263 : get<hydro::Tags::SpecificInternalEnergy<DataVector>>(temp_hydro_vars) =
264 : equation_of_state
265 : .specific_internal_energy_from_density_and_temperature(
266 : get<RestMassDensity>(ghost_data_vars),
267 : get<Temperature>(ghost_data_vars),
268 : get<ElectronFraction>(ghost_data_vars));
269 : get<hydro::Tags::Pressure<DataVector>>(temp_hydro_vars) =
270 : equation_of_state.pressure_from_density_and_temperature(
271 : get<RestMassDensity>(ghost_data_vars),
272 : get<Temperature>(ghost_data_vars),
273 : get<ElectronFraction>(ghost_data_vars));
274 :
275 : atmosphere_fixer(
276 : make_not_null(&get<RestMassDensity>(ghost_data_vars)),
277 : make_not_null(&get<hydro::Tags::SpecificInternalEnergy<DataVector>>(
278 : temp_hydro_vars)),
279 : make_not_null(&get<hydro::Tags::SpatialVelocity<DataVector, 3>>(
280 : temp_hydro_vars)),
281 : make_not_null(
282 : &get<hydro::Tags::LorentzFactor<DataVector>>(temp_hydro_vars)),
283 : make_not_null(
284 : &get<hydro::Tags::Pressure<DataVector>>(temp_hydro_vars)),
285 : make_not_null(&get<Temperature>(ghost_data_vars)),
286 : get<ElectronFraction>(ghost_data_vars), spatial_metric,
287 : equation_of_state);
288 : }
289 : } // for (direction : external boundaries)
290 : }
291 : } // namespace grmhd::GhValenciaDivClean::fd
|
