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 <utility>
9 :
10 : #include "DataStructures/DataBox/DataBox.hpp"
11 : #include "DataStructures/DataBox/MetavariablesTag.hpp"
12 : #include "DataStructures/DataBox/PrefixHelpers.hpp"
13 : #include "DataStructures/DataBox/Prefixes.hpp"
14 : #include "DataStructures/TaggedTuple.hpp"
15 : #include "DataStructures/Variables.hpp"
16 : #include "DataStructures/VariablesTag.hpp"
17 : #include "Domain/BlockLogicalCoordinates.hpp"
18 : #include "Domain/Creators/Tags/Domain.hpp"
19 : #include "Domain/Domain.hpp"
20 : #include "Domain/Structure/ElementId.hpp"
21 : #include "Domain/Tags.hpp"
22 : #include "Elliptic/DiscontinuousGalerkin/Tags.hpp"
23 : #include "Elliptic/Systems/SelfForce/GeneralRelativity/AnalyticData/CircularOrbit.hpp"
24 : #include "Elliptic/Systems/SelfForce/GeneralRelativity/Equations.hpp"
25 : #include "Elliptic/Systems/SelfForce/GeneralRelativity/Tags.hpp"
26 : #include "NumericalAlgorithms/DiscontinuousGalerkin/ApplyMassMatrix.hpp"
27 : #include "NumericalAlgorithms/DiscontinuousGalerkin/MortarHelpers.hpp"
28 : #include "NumericalAlgorithms/DiscontinuousGalerkin/NormalDotFlux.hpp"
29 : #include "NumericalAlgorithms/DiscontinuousGalerkin/Tags.hpp"
30 : #include "Parallel/AlgorithmExecution.hpp"
31 : #include "ParallelAlgorithms/Amr/Protocols/Projector.hpp"
32 : #include "Utilities/CallWithDynamicType.hpp"
33 : #include "Utilities/MakeWithValue.hpp"
34 : #include "Utilities/TMPL.hpp"
35 :
36 : /// \cond
37 : namespace Parallel {
38 : template <typename Metavariables>
39 : struct GlobalCache;
40 : } // namespace Parallel
41 : /// \endcond
42 :
43 0 : namespace GrSelfForce::Actions {
44 :
45 : /*!
46 : * \brief Initialize the effective source and singular field for the
47 : * gravitational self-force problem in a regularized region around the puncture.
48 : * Replaces 'elliptic::Actions::InitializeFixedSources' in the standard elliptic
49 : * initialization.
50 : *
51 : * \details The regularized region is chosen as the full block containing the
52 : * puncture. In this region, the `Tags::FieldIsRegularized` is set to true and
53 : * the effective source and the singular field are set by the
54 : * GrSelfForce::AnalyticData::CircularOrbit class. Outside this region, the
55 : * fixed source and singular field are set to zero.
56 : */
57 : template <typename System, typename BackgroundTag>
58 1 : struct InitializeEffectiveSource : tt::ConformsTo<::amr::protocols::Projector> {
59 : private:
60 0 : static constexpr size_t Dim = System::volume_dim;
61 0 : using fixed_sources_tag = ::Tags::Variables<
62 : db::wrap_tags_in<::Tags::FixedSource, typename System::primal_fields>>;
63 0 : using singular_vars_tag = ::Tags::Variables<tmpl::list<
64 : Tags::SingularField,
65 : ::Tags::deriv<Tags::SingularField, tmpl::size_t<Dim>, Frame::Inertial>>>;
66 0 : using singular_vars_on_mortars_tag =
67 : ::Tags::Variables<tmpl::list<Tags::SingularField,
68 : ::Tags::NormalDotFlux<Tags::SingularField>>>;
69 0 : using analytic_tags_list = tmpl::push_back<
70 : typename fixed_sources_tag::tags_list, Tags::SingularField,
71 : ::Tags::deriv<Tags::SingularField, tmpl::size_t<Dim>, Frame::Inertial>,
72 : Tags::BoyerLindquistRadius>;
73 :
74 : public: // Iterable action
75 0 : using const_global_cache_tags =
76 : tmpl::list<elliptic::dg::Tags::Massive, BackgroundTag,
77 : Tags::NullSlicingBlocks<Dim>>;
78 0 : using simple_tags =
79 : tmpl::list<fixed_sources_tag, singular_vars_tag,
80 : ::Tags::Mortars<singular_vars_on_mortars_tag, Dim>,
81 : Tags::BoyerLindquistRadius, Tags::FieldIsRegularized,
82 : ::Tags::Mortars<Tags::FieldIsRegularized, Dim>>;
83 0 : using compute_tags = tmpl::list<>;
84 :
85 : template <typename DbTagsList, typename... InboxTags, typename Metavariables,
86 : typename ActionList, typename ParallelComponent>
87 0 : static Parallel::iterable_action_return_t apply(
88 : db::DataBox<DbTagsList>& box,
89 : const tuples::TaggedTuple<InboxTags...>& /*inboxes*/,
90 : const Parallel::GlobalCache<Metavariables>& /*cache*/,
91 : const ElementId<Dim>& /*array_index*/, const ActionList /*meta*/,
92 : const ParallelComponent* const /*meta*/) {
93 : db::mutate_apply<InitializeEffectiveSource>(make_not_null(&box));
94 : return {Parallel::AlgorithmExecution::Continue, std::nullopt};
95 : }
96 :
97 : public: // DataBox mutator, amr::protocols::Projector
98 0 : using return_tags = simple_tags;
99 0 : using argument_tags = tmpl::list<
100 : domain::Tags::Coordinates<Dim, Frame::Inertial>,
101 : domain::Tags::Domain<Dim>, domain::Tags::Element<Dim>,
102 : ::Tags::Mortars<domain::Tags::Coordinates<Dim, Frame::Inertial>, Dim>,
103 : BackgroundTag, elliptic::dg::Tags::Massive, domain::Tags::Mesh<Dim>,
104 : domain::Tags::DetInvJacobian<Frame::ElementLogical, Frame::Inertial>,
105 : Parallel::Tags::Metavariables>;
106 :
107 : template <typename Background, typename Metavariables, typename... AmrData>
108 0 : static void apply(
109 : const gsl::not_null<typename fixed_sources_tag::type*> fixed_sources,
110 : const gsl::not_null<typename singular_vars_tag::type*> singular_vars,
111 : const gsl::not_null<
112 : DirectionalIdMap<Dim, typename singular_vars_on_mortars_tag::type>*>
113 : singular_vars_on_mortars,
114 : const gsl::not_null<Scalar<DataVector>*> bl_radius,
115 : const gsl::not_null<bool*> field_is_regularized,
116 : const gsl::not_null<DirectionalIdMap<Dim, bool>*>
117 : neighbors_field_is_regularized,
118 : const tnsr::I<DataVector, Dim>& inertial_coords,
119 : const Domain<Dim>& domain, const Element<Dim>& element,
120 : const DirectionalIdMap<Dim, tnsr::I<DataVector, Dim>>&
121 : all_mortar_inertial_coords,
122 : const Background& background, const bool massive, const Mesh<Dim>& mesh,
123 : const Scalar<DataVector>& det_inv_jacobian, const Metavariables& /*meta*/,
124 : const AmrData&... /*amr_data*/) {
125 : const auto& circular_orbit =
126 : dynamic_cast<const GrSelfForce::AnalyticData::CircularOrbit&>(
127 : background);
128 :
129 : // Determine the regularized region
130 : const tnsr::I<double, Dim> puncture_pos =
131 : circular_orbit.puncture_position();
132 : const auto puncture_in_element =
133 : [&puncture_pos, &domain](const ElementId<Dim>& element_id) -> bool {
134 : // Regularize full block that contains the puncture
135 : const auto& block = domain.blocks()[element_id.block_id()];
136 : const auto block_logical_coords =
137 : block_logical_coordinates_single_point(puncture_pos, block);
138 : return block_logical_coords.has_value();
139 : };
140 : *field_is_regularized = puncture_in_element(element.id());
141 : neighbors_field_is_regularized->clear();
142 : for (const auto& [direction, neighbors] : element.neighbors()) {
143 : for (const auto& neighbor_id : neighbors) {
144 : const dg::MortarId<Dim> mortar_id{direction, neighbor_id};
145 : neighbors_field_is_regularized->emplace(
146 : mortar_id, puncture_in_element(neighbor_id));
147 : }
148 : }
149 :
150 : // Set the effective source if solving for the regular field
151 : const size_t num_points = mesh.number_of_grid_points();
152 : if (*field_is_regularized) {
153 : const auto vars =
154 : circular_orbit.variables(inertial_coords, analytic_tags_list{});
155 : fixed_sources->initialize(num_points);
156 : singular_vars->initialize(num_points);
157 : get<::Tags::FixedSource<Tags::MMode>>(*fixed_sources) =
158 : get<::Tags::FixedSource<Tags::MMode>>(vars);
159 : get<Tags::SingularField>(*singular_vars) = get<Tags::SingularField>(vars);
160 : get<::Tags::deriv<Tags::SingularField, tmpl::size_t<Dim>,
161 : Frame::Inertial>>(*singular_vars) =
162 : get<::Tags::deriv<Tags::SingularField, tmpl::size_t<Dim>,
163 : Frame::Inertial>>(vars);
164 : *bl_radius = get<Tags::BoyerLindquistRadius>(vars);
165 : // Apply DG mass matrix to the fixed sources if the DG operator is massive
166 : if (massive) {
167 : *fixed_sources /= get(det_inv_jacobian);
168 : ::dg::apply_mass_matrix(fixed_sources, mesh);
169 : }
170 : } else {
171 : *fixed_sources =
172 : Variables<typename fixed_sources_tag::tags_list>{num_points, 0.};
173 : *singular_vars =
174 : Variables<typename singular_vars_tag::tags_list>{num_points, 0.};
175 : *bl_radius = Scalar<DataVector>{num_points, 0.};
176 : }
177 :
178 : // Set the singular field and flux on mortars, which is needed to transform
179 : // between regularized and full fields
180 : singular_vars_on_mortars->clear();
181 : for (const auto& [mortar_id, mortar_inertial_coords] :
182 : all_mortar_inertial_coords) {
183 : if (*field_is_regularized ==
184 : neighbors_field_is_regularized->at(mortar_id)) {
185 : // Both elements solve for the same field. No transformation needed.
186 : continue;
187 : }
188 : const auto vars_on_mortar = circular_orbit.variables(
189 : mortar_inertial_coords, analytic_tags_list{});
190 : const auto background_on_mortar = circular_orbit.variables(
191 : mortar_inertial_coords,
192 : tmpl::list<Tags::Alpha, Tags::Beta, Tags::GammaRstar,
193 : Tags::GammaTheta>{});
194 : auto& singular_vars_on_mortar = (*singular_vars_on_mortars)[mortar_id];
195 : singular_vars_on_mortar.initialize(
196 : mortar_inertial_coords.begin()->size());
197 : get<Tags::SingularField>(singular_vars_on_mortar) =
198 : get<Tags::SingularField>(vars_on_mortar);
199 : const auto& deriv_singular_field_on_mortar =
200 : get<::Tags::deriv<Tags::SingularField, tmpl::size_t<Dim>,
201 : Frame::Inertial>>(vars_on_mortar);
202 : const auto& alpha_on_mortar = get<Tags::Alpha>(background_on_mortar);
203 : FluxTensorType singular_field_flux_on_mortar{};
204 : GrSelfForce::fluxes(make_not_null(&singular_field_flux_on_mortar),
205 : alpha_on_mortar, deriv_singular_field_on_mortar);
206 : // Assuming mortar normal is just (1, 0) or (0, 1). This is true for the
207 : // rectangular domains used in the self-force problem (AlignedLattice in
208 : // r, theta coordinates).
209 : tnsr::i<DataVector, Dim> mortar_normal{
210 : mortar_inertial_coords.begin()->size(), 0.};
211 : mortar_normal.get(mortar_id.direction().dimension()) =
212 : mortar_id.direction().sign();
213 : normal_dot_flux(
214 : make_not_null(&get<::Tags::NormalDotFlux<Tags::SingularField>>(
215 : singular_vars_on_mortar)),
216 : mortar_normal, singular_field_flux_on_mortar);
217 : }
218 : }
219 : };
220 :
221 : } // namespace GrSelfForce::Actions
|
