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