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/DataVector.hpp"
14 : #include "DataStructures/Variables.hpp"
15 : #include "Domain/BoundaryConditions/BoundaryCondition.hpp"
16 : #include "Domain/BoundaryConditions/None.hpp"
17 : #include "Domain/BoundaryConditions/Periodic.hpp"
18 : #include "Domain/Creators/Tags/ExternalBoundaryConditions.hpp"
19 : #include "Domain/Domain.hpp"
20 : #include "Domain/Structure/Direction.hpp"
21 : #include "Domain/Structure/Element.hpp"
22 : #include "Domain/Structure/ElementId.hpp"
23 : #include "Domain/Tags.hpp"
24 : #include "Domain/TagsTimeDependent.hpp"
25 : #include "Evolution/BoundaryConditions/Type.hpp"
26 : #include "Evolution/DgSubcell/Tags/CellCenteredFlux.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/ValenciaDivClean/BoundaryConditions/BoundaryCondition.hpp"
31 : #include "Evolution/Systems/GrMhd/ValenciaDivClean/BoundaryConditions/Factory.hpp"
32 : #include "Evolution/Systems/GrMhd/ValenciaDivClean/FiniteDifference/Reconstructor.hpp"
33 : #include "Evolution/Systems/GrMhd/ValenciaDivClean/System.hpp"
34 : #include "NumericalAlgorithms/Spectral/Mesh.hpp"
35 : #include "Parallel/Tags/Metavariables.hpp"
36 : #include "PointwiseFunctions/Hydro/Tags.hpp"
37 : #include "Utilities/CallWithDynamicType.hpp"
38 : #include "Utilities/ErrorHandling/Assert.hpp"
39 : #include "Utilities/Gsl.hpp"
40 : #include "Utilities/PrettyType.hpp"
41 : #include "Utilities/TMPL.hpp"
42 :
43 1 : namespace grmhd::ValenciaDivClean::fd {
44 : /*!
45 : * \brief Computes finite difference ghost data for external boundary
46 : * conditions.
47 : *
48 : * If the element is at the external boundary, computes FD ghost data with a
49 : * given boundary condition and stores it into neighbor data with {direction,
50 : * ElementId::external_boundary_id()} as the mortar_id key.
51 : *
52 : * \note Subcell needs to be enabled for boundary elements. Otherwise this
53 : * function would be never called.
54 : *
55 : */
56 1 : struct BoundaryConditionGhostData {
57 : template <typename DbTagsList>
58 0 : static void apply(const gsl::not_null<db::DataBox<DbTagsList>*> box,
59 : const Element<3>& element,
60 : const Reconstructor& reconstructor);
61 :
62 : private:
63 : template <typename FdBoundaryConditionHelper, typename DbTagsList,
64 : typename... FdBoundaryConditionArgsTags>
65 : // A helper function for calling fd_ghost() of BoundaryCondition subclasses
66 0 : static void apply_subcell_boundary_condition_impl(
67 : FdBoundaryConditionHelper& fd_boundary_condition_helper,
68 : const gsl::not_null<db::DataBox<DbTagsList>*>& box,
69 : tmpl::list<FdBoundaryConditionArgsTags...>) {
70 : return fd_boundary_condition_helper(
71 : db::get<FdBoundaryConditionArgsTags>(*box)...);
72 : }
73 : };
74 :
75 : template <typename DbTagsList>
76 : void BoundaryConditionGhostData::apply(
77 : const gsl::not_null<db::DataBox<DbTagsList>*> box,
78 : const Element<3>& element, const Reconstructor& reconstructor) {
79 : const auto& external_boundary_condition =
80 : db::get<domain::Tags::ExternalBoundaryConditions<3>>(*box).at(
81 : element.id().block_id());
82 :
83 : // Check if the element is on the external boundary. If not, the caller is
84 : // doing something wrong (e.g. trying to compute FD ghost data with boundary
85 : // conditions at an element which is not on the external boundary).
86 : ASSERT(not element.external_boundaries().empty(),
87 : "The element (ID : " << element.id()
88 : << ") is not on external boundaries");
89 :
90 : const Mesh<3> subcell_mesh =
91 : db::get<evolution::dg::subcell::Tags::Mesh<3>>(*box);
92 :
93 : const size_t ghost_zone_size{reconstructor.ghost_zone_size()};
94 :
95 : // Tags and tags list for FD reconstruction
96 : using RestMassDensity = hydro::Tags::RestMassDensity<DataVector>;
97 : using ElectronFraction = hydro::Tags::ElectronFraction<DataVector>;
98 : using Temperature = hydro::Tags::Temperature<DataVector>;
99 : using LorentzFactorTimesSpatialVelocity =
100 : hydro::Tags::LorentzFactorTimesSpatialVelocity<DataVector, 3>;
101 : using MagneticField = hydro::Tags::MagneticField<DataVector, 3>;
102 : using DivergenceCleaningField =
103 : hydro::Tags::DivergenceCleaningField<DataVector>;
104 :
105 : using prims_for_reconstruction =
106 : tmpl::list<RestMassDensity, ElectronFraction, Temperature,
107 : LorentzFactorTimesSpatialVelocity, MagneticField,
108 : DivergenceCleaningField>;
109 :
110 : size_t num_prims_tensor_components = 0;
111 : tmpl::for_each<prims_for_reconstruction>([&num_prims_tensor_components](
112 : auto tag) {
113 : num_prims_tensor_components += tmpl::type_from<decltype(tag)>::type::size();
114 : });
115 :
116 : using flux_variables =
117 : typename grmhd::ValenciaDivClean::System::flux_variables;
118 : const bool compute_cell_centered_flux =
119 : db::get<
120 : evolution::dg::subcell::Tags::CellCenteredFlux<flux_variables, 3>>(
121 : *box)
122 : .has_value();
123 :
124 : for (const auto& direction : element.external_boundaries()) {
125 : const auto& boundary_condition_at_direction =
126 : *external_boundary_condition.at(direction);
127 :
128 : const size_t num_face_pts{
129 : subcell_mesh.extents().slice_away(direction.dimension()).product()};
130 :
131 : // Allocate a vector to store the computed FD ghost data and assign a
132 : // non-owning Variables on it.
133 : using FluxVars =
134 : Variables<db::wrap_tags_in<::Tags::Flux, flux_variables,
135 : tmpl::size_t<3>, Frame::Inertial>>;
136 : const size_t prims_size =
137 : num_prims_tensor_components * ghost_zone_size * num_face_pts;
138 : const size_t fluxes_size =
139 : (compute_cell_centered_flux ? FluxVars::number_of_independent_components
140 : : 0) *
141 : ghost_zone_size * num_face_pts;
142 :
143 : auto& all_ghost_data = db::get_mutable_reference<
144 : evolution::dg::subcell::Tags::GhostDataForReconstruction<3>>(box);
145 : // Put the computed ghost data into neighbor data with {direction,
146 : // ElementId::external_boundary_id()} as the mortar_id key
147 : const DirectionalId<3> mortar_id{direction,
148 : ElementId<3>::external_boundary_id()};
149 :
150 : all_ghost_data[mortar_id] = evolution::dg::subcell::GhostData{1};
151 : DataVector& boundary_ghost_data =
152 : all_ghost_data.at(mortar_id).neighbor_ghost_data_for_reconstruction();
153 : boundary_ghost_data.destructive_resize(prims_size + fluxes_size);
154 : Variables<prims_for_reconstruction> ghost_data_vars{
155 : boundary_ghost_data.data(), prims_size};
156 : std::optional<FluxVars> cell_centered_ghost_fluxes{};
157 : if (compute_cell_centered_flux) {
158 : cell_centered_ghost_fluxes = FluxVars{};
159 : cell_centered_ghost_fluxes.value().set_data_ref(
160 : std::next(boundary_ghost_data.data(),
161 : static_cast<std::ptrdiff_t>(prims_size)),
162 : fluxes_size);
163 : }
164 :
165 : // We don't need to care about boundary ghost data when using the periodic
166 : // condition, so exclude it from the type list
167 : using factory_classes =
168 : typename std::decay_t<decltype(db::get<Parallel::Tags::Metavariables>(
169 : *box))>::factory_creation::factory_classes;
170 : using derived_boundary_conditions_for_subcell = tmpl::remove_if<
171 : tmpl::at<factory_classes, typename System::boundary_conditions_base>,
172 : tmpl::or_<
173 : std::is_base_of<domain::BoundaryConditions::MarkAsPeriodic,
174 : tmpl::_1>,
175 : std::is_base_of<domain::BoundaryConditions::MarkAsNone, tmpl::_1>>>;
176 :
177 : // Now apply subcell boundary conditions
178 : call_with_dynamic_type<void, derived_boundary_conditions_for_subcell>(
179 : &boundary_condition_at_direction,
180 : [&box, &cell_centered_ghost_fluxes, &direction,
181 : &ghost_data_vars](const auto* boundary_condition) {
182 : using BoundaryCondition = std::decay_t<decltype(*boundary_condition)>;
183 : using bcondition_interior_evolved_vars_tags =
184 : typename BoundaryCondition::fd_interior_evolved_variables_tags;
185 : using bcondition_interior_temporary_tags =
186 : typename BoundaryCondition::fd_interior_temporary_tags;
187 : using bcondition_interior_primitive_vars_tags =
188 : typename BoundaryCondition::fd_interior_primitive_variables_tags;
189 : using bcondition_gridless_tags =
190 : typename BoundaryCondition::fd_gridless_tags;
191 :
192 : using bcondition_interior_tags =
193 : tmpl::append<bcondition_interior_evolved_vars_tags,
194 : bcondition_interior_temporary_tags,
195 : bcondition_interior_primitive_vars_tags,
196 : bcondition_gridless_tags>;
197 :
198 : if constexpr (BoundaryCondition::bc_type ==
199 : evolution::BoundaryConditions::Type::Ghost) {
200 : const auto apply_fd_ghost =
201 : [&boundary_condition, &cell_centered_ghost_fluxes, &direction,
202 : &ghost_data_vars](const auto&... boundary_ghost_data_args) {
203 : (*boundary_condition)
204 : .fd_ghost(
205 : make_not_null(&get<RestMassDensity>(ghost_data_vars)),
206 : make_not_null(
207 : &get<ElectronFraction>(ghost_data_vars)),
208 : make_not_null(&get<Temperature>(ghost_data_vars)),
209 : make_not_null(&get<LorentzFactorTimesSpatialVelocity>(
210 : ghost_data_vars)),
211 : make_not_null(&get<MagneticField>(ghost_data_vars)),
212 : make_not_null(
213 : &get<DivergenceCleaningField>(ghost_data_vars)),
214 : make_not_null(&cell_centered_ghost_fluxes), direction,
215 : boundary_ghost_data_args...);
216 : };
217 : apply_subcell_boundary_condition_impl(apply_fd_ghost, box,
218 : bcondition_interior_tags{});
219 : } else if constexpr (BoundaryCondition::bc_type ==
220 : evolution::BoundaryConditions::Type::
221 : DemandOutgoingCharSpeeds) {
222 : // This boundary condition only checks if all the characteristic
223 : // speeds are directed outward.
224 : const auto& volume_mesh_velocity =
225 : db::get<domain::Tags::MeshVelocity<3, Frame::Inertial>>(*box);
226 : if (volume_mesh_velocity.has_value()) {
227 : ERROR("Subcell currently does not support moving mesh");
228 : }
229 :
230 : std::optional<tnsr::I<DataVector, 3>> face_mesh_velocity{};
231 :
232 : const auto& normal_covector_and_magnitude =
233 : db::get<evolution::dg::Tags::NormalCovectorAndMagnitude<3>>(
234 : *box);
235 : const auto outward_directed_normal_covector =
236 : get<evolution::dg::Tags::NormalCovector<3>>(
237 : normal_covector_and_magnitude.at(direction).value());
238 :
239 : const auto apply_fd_demand_outgoing_char_speeds =
240 : [&boundary_condition, &cell_centered_ghost_fluxes, &direction,
241 : &face_mesh_velocity, &ghost_data_vars,
242 : &outward_directed_normal_covector](
243 : const auto&... boundary_ghost_data_args) {
244 : return (*boundary_condition)
245 : .fd_demand_outgoing_char_speeds(
246 : make_not_null(&get<RestMassDensity>(ghost_data_vars)),
247 : make_not_null(
248 : &get<ElectronFraction>(ghost_data_vars)),
249 : make_not_null(&get<Temperature>(ghost_data_vars)),
250 : make_not_null(&get<LorentzFactorTimesSpatialVelocity>(
251 : ghost_data_vars)),
252 : make_not_null(&get<MagneticField>(ghost_data_vars)),
253 : make_not_null(
254 : &get<DivergenceCleaningField>(ghost_data_vars)),
255 : make_not_null(&cell_centered_ghost_fluxes), direction,
256 : face_mesh_velocity, outward_directed_normal_covector,
257 : boundary_ghost_data_args...);
258 : };
259 : apply_subcell_boundary_condition_impl(
260 : apply_fd_demand_outgoing_char_speeds, box,
261 : bcondition_interior_tags{});
262 :
263 : return;
264 : } else {
265 : ERROR("Unsupported boundary condition "
266 : << pretty_type::short_name<BoundaryCondition>()
267 : << " when using finite-difference");
268 : }
269 : });
270 : }
271 : }
272 : } // namespace grmhd::ValenciaDivClean::fd
|
