Line data    Source code

       1           0 : // Distributed under the MIT License.
       2             : // See LICENSE.txt for details.
       3             : 
       4             : #pragma once
       5             : 
       6             : #include <boost/math/special_functions/spherical_harmonic.hpp>
       7             : #include <cstddef>
       8             : #include <optional>
       9             : #include <vector>
      10             : 
      11             : #include "DataStructures/DataBox/Tag.hpp"
      12             : #include "DataStructures/DataVector.hpp"
      13             : #include "DataStructures/Tags/TempTensor.hpp"
      14             : #include "DataStructures/Tensor/EagerMath/DotProduct.hpp"
      15             : #include "DataStructures/Tensor/Slice.hpp"
      16             : #include "DataStructures/Variables.hpp"
      17             : #include "Domain/AreaElement.hpp"
      18             : #include "Domain/ExcisionSphere.hpp"
      19             : #include "Domain/Structure/Direction.hpp"
      20             : #include "Domain/Structure/IndexToSliceAt.hpp"
      21             : #include "Domain/Tags.hpp"
      22             : #include "Domain/TagsTimeDependent.hpp"
      23             : #include "Evolution/Systems/CurvedScalarWave/System.hpp"
      24             : #include "Evolution/Systems/CurvedScalarWave/Tags.hpp"
      25             : #include "Evolution/Systems/CurvedScalarWave/Worldtube/Inboxes.hpp"
      26             : #include "Evolution/Systems/CurvedScalarWave/Worldtube/SingletonChare.hpp"
      27             : #include "Evolution/Systems/CurvedScalarWave/Worldtube/Tags.hpp"
      28             : #include "NumericalAlgorithms/LinearOperators/DefiniteIntegral.hpp"
      29             : #include "NumericalAlgorithms/LinearOperators/PartialDerivatives.hpp"
      30             : #include "NumericalAlgorithms/SphericalHarmonics/RealSphericalHarmonics.hpp"
      31             : #include "Parallel/AlgorithmExecution.hpp"
      32             : #include "Parallel/GlobalCache.hpp"
      33             : #include "Parallel/Invoke.hpp"
      34             : #include "PointwiseFunctions/GeneralRelativity/Tags.hpp"
      35             : #include "Utilities/ErrorHandling/Assert.hpp"
      36             : #include "Utilities/Gsl.hpp"
      37             : #include "Utilities/TMPL.hpp"
      38             : 
      39             : /// \cond
      40             : namespace Tags {
      41             : struct TimeStepId;
      42             : }  // namespace Tags
      43             : /// \endcond
      44             : 
      45             : namespace CurvedScalarWave::Worldtube::Actions {
      46             : /*!
      47             :  * \brief Projects the regular field \f$\Psi^R\f$ and its time derivative
      48             :  * \f$\partial_t \Psi^R\f$ onto real spherical harmonics and sends the result to
      49             :  * the worldtube.
      50             :  *
      51             :  * \details The regular field is obtained by subtracting the singular/puncture
      52             :  * field from the numerical DG field.
      53             :  * All spherical harmonics are computed for \f$l <= n\f$, where \f$n\f$ is the
      54             :  * worldtube expansion order. The projection is done by integrating over the DG
      55             :  * grid of the element face using \ref definite_integral with the euclidean area
      56             :  * element. The worldtube adds up all integrals from the different elements to
      57             :  * obtain the integral over the entire sphere.
      58             :  *
      59             :  * DataBox:
      60             :  * - Uses:
      61             :  *    - `tags_to_slice_on_face`
      62             :  *    - `Worldtube::Tags::ExpansionOrder`
      63             :  *    - `Worldtube::Tags::FaceCoordinates<Dim, Frame::Grid, true>`
      64             :  *    - `Worldtube::Tags::PunctureField`
      65             :  *    - `Worldtube::Tags::ExcisionSphere`
      66             :  *    - `Tags::TimeStepId`
      67             :  */
      68           1 : struct SendToWorldtube {
      69           0 :   static constexpr size_t Dim = 3;
      70           0 :   using tags_to_send = tmpl::list<CurvedScalarWave::Tags::Psi,
      71             :                                   ::Tags::dt<CurvedScalarWave::Tags::Psi>>;
      72           0 :   using tags_to_slice_to_face = tmpl::list<
      73             :       CurvedScalarWave::Tags::Psi, CurvedScalarWave::Tags::Pi,
      74             :       CurvedScalarWave::Tags::Phi<Dim>, gr::Tags::Shift<DataVector, Dim>,
      75             :       gr::Tags::Lapse<DataVector>,
      76             :       domain::Tags::InverseJacobian<Dim, Frame::ElementLogical, Frame::Grid>>;
      77             : 
      78           0 :   using inbox_tags = tmpl::list<Worldtube::Tags::SphericalHarmonicsInbox<Dim>>;
      79           0 :   using simple_tags = tmpl::list<Tags::RegularFieldAdvectiveTerm<Dim>>;
      80             : 
      81             :   template <typename DbTagsList, typename... InboxTags, typename Metavariables,
      82             :             typename ArrayIndex, typename ActionList,
      83             :             typename ParallelComponent>
      84           0 :   static Parallel::iterable_action_return_t apply(
      85             :       db::DataBox<DbTagsList>& box,
      86             :       tuples::TaggedTuple<InboxTags...>& /*inboxes*/,
      87             :       Parallel::GlobalCache<Metavariables>& cache,
      88             :       const ArrayIndex& /*array_index*/, const ActionList /*meta*/,
      89             :       const ParallelComponent* const /*meta*/) {
      90             :     const auto& puncture_field =
      91             :         db::get<Worldtube::Tags::PunctureField<Dim>>(box);
      92             :     if (puncture_field.has_value()) {
      93             :       const auto& element_id = db::get<domain::Tags::Element<Dim>>(box).id();
      94             :       const auto& excision_sphere = db::get<Tags::ExcisionSphere<Dim>>(box);
      95             :       const auto direction = excision_sphere.abutting_direction(element_id);
      96             :       ASSERT(direction.has_value(), "Should be abutting");
      97             :       const auto& mesh = db::get<domain::Tags::Mesh<Dim>>(box);
      98             :       const auto face_mesh = mesh.slice_away(direction->dimension());
      99             :       const size_t face_size = face_mesh.number_of_grid_points();
     100             : 
     101             :       Variables<tmpl::push_back<tags_to_slice_to_face, ::Tags::TempScalar<0>>>
     102             :           vars_on_face(face_size);
     103             : 
     104             :       tmpl::for_each<tags_to_slice_to_face>(
     105             :           [&box, &vars_on_face, &mesh, &direction](auto tag_to_slice_v) {
     106             :             using tag_to_slice = typename decltype(tag_to_slice_v)::type;
     107             :             data_on_slice(make_not_null(&get<tag_to_slice>(vars_on_face)),
     108             :                           db::get<tag_to_slice>(box), mesh.extents(),
     109             :                           direction.value().dimension(),
     110             :                           index_to_slice_at(mesh.extents(), direction.value()));
     111             :           });
     112             :       const auto& face_lapse = get<gr::Tags::Lapse<DataVector>>(vars_on_face);
     113             :       const auto& face_shift =
     114             :           get<gr::Tags::Shift<DataVector, Dim>>(vars_on_face);
     115             :       auto& face_inv_jacobian =
     116             :           get<domain::Tags::InverseJacobian<Dim, Frame::ElementLogical,
     117             :                                             Frame::Grid>>(vars_on_face);
     118             :       const auto& face_psi = get<CurvedScalarWave::Tags::Psi>(vars_on_face);
     119             :       const auto& face_pi = get<CurvedScalarWave::Tags::Pi>(vars_on_face);
     120             :       const auto& face_phi =
     121             :           get<CurvedScalarWave::Tags::Phi<Dim>>(vars_on_face);
     122             :       // re-use allocations
     123             :       Scalar<DataVector>& area_element =
     124             :           get<::Tags::TempScalar<0>>(vars_on_face);
     125             :       euclidean_area_element(make_not_null(&area_element), face_inv_jacobian,
     126             :                              direction.value());
     127             :       // re-use allocations
     128             :       DataVector& psi_regular_times_det = get<0, 0>(face_inv_jacobian);
     129             :       DataVector& dt_psi_regular_times_det = get<0, 1>(face_inv_jacobian);
     130             :       // the regular field is the full numerical field minus the puncture field
     131             :       psi_regular_times_det =
     132             :           get(face_psi) -
     133             :           get(get<CurvedScalarWave::Tags::Psi>(puncture_field.value()));
     134             : 
     135             :       // transform Pi to dt Psi. This is equivalent to the evolution equation
     136             :       // for dt Psi but we need to calculate it again because boundary
     137             :       // corrections from the DG scheme are already applied to the stored value.
     138             :       dt_psi_regular_times_det =
     139             :           -get(face_lapse) * get(face_pi) +
     140             :           get(dot_product(face_shift, face_phi)) -
     141             :           get(get<::Tags::dt<CurvedScalarWave::Tags::Psi>>(
     142             :               puncture_field.value()));
     143             : 
     144             :       const auto& mesh_velocity = db::get<domain::Tags::MeshVelocity<Dim>>(box);
     145             :       ASSERT(mesh_velocity.has_value(),
     146             :              "Expected a moving grid for worldrube evolution.");
     147             :       // is an optional so we can't put the tag into variables. The shift is not
     148             :       // used at this point so we use the allocation for the mesh_velocity to
     149             :       // save memory.
     150             :       auto& mesh_velocity_on_face =
     151             :           get<gr::Tags::Shift<DataVector, Dim>>(vars_on_face);
     152             :       data_on_slice(make_not_null(&mesh_velocity_on_face),
     153             :                     mesh_velocity.value(), mesh.extents(),
     154             :                     direction.value().dimension(),
     155             :                     index_to_slice_at(mesh.extents(), direction.value()));
     156             :       db::mutate<Tags::RegularFieldAdvectiveTerm<Dim>>(
     157             :           [&face_phi, &mesh_velocity_on_face,
     158             :            &di_psi_puncture =
     159             :                get<::Tags::deriv<CurvedScalarWave::Tags::Psi, tmpl::size_t<3>,
     160             :                                  Frame::Inertial>>(puncture_field.value())](
     161             :               const gsl::not_null<Scalar<DataVector>*> regular_advective_term) {
     162             :             tenex::evaluate<>(regular_advective_term,
     163             :                               (face_phi(ti::i) - di_psi_puncture(ti::i)) *
     164             :                                   mesh_velocity_on_face(ti::I));
     165             :           },
     166             :           make_not_null(&box));
     167             :       // The time derivative is transformed into the grid frame using the
     168             :       // advective term which comes from the transformation of the time
     169             :       // derivative due to the moving mesh.
     170             :       dt_psi_regular_times_det +=
     171             :           get(get<Tags::RegularFieldAdvectiveTerm<Dim>>(box));
     172             : 
     173             :       psi_regular_times_det *= get(area_element);
     174             :       dt_psi_regular_times_det *= get(area_element);
     175             :       const auto& centered_face_coords =
     176             :           db::get<Tags::FaceCoordinates<Dim, Frame::Grid, true>>(box);
     177             :       ASSERT(centered_face_coords.has_value(),
     178             :              "Should be an abutting element here, but face coords are not "
     179             :              "calculated!");
     180             :       const auto& x = get<0>(centered_face_coords.value());
     181             :       const auto& y = get<1>(centered_face_coords.value());
     182             :       const auto& z = get<2>(centered_face_coords.value());
     183             : 
     184             :       const size_t order = db::get<Worldtube::Tags::ExpansionOrder>(box);
     185             :       const size_t num_modes = (order + 1) * (order + 1);
     186             :       Variables<tags_to_send> Ylm_coefs(num_modes);
     187             :       // re-use allocations
     188             :       auto& theta = get<0, 2>(face_inv_jacobian);
     189             :       auto& phi = get<1, 0>(face_inv_jacobian);
     190             :       auto& spherical_harmonic = get<1, 1>(face_inv_jacobian);
     191             : 
     192             :       theta = atan2(hypot(x, y), z);
     193             :       phi = atan2(y, x);
     194             :       size_t index = 0;
     195             :       // project onto spherical harmonics
     196             :       for (size_t l = 0; l <= order; ++l) {
     197             :         for (int m = -l; m <= static_cast<int>(l); ++m, ++index) {
     198             :           spherical_harmonic = ylm::real_spherical_harmonic(theta, phi, l, m);
     199             :           get(get<CurvedScalarWave::Tags::Psi>(Ylm_coefs)).at(index) =
     200             :               definite_integral(psi_regular_times_det * spherical_harmonic,
     201             :                                 face_mesh);
     202             :           get(get<::Tags::dt<CurvedScalarWave::Tags::Psi>>(Ylm_coefs))
     203             :               .at(index) = definite_integral(
     204             :               dt_psi_regular_times_det * spherical_harmonic, face_mesh);
     205             :         }
     206             :       }
     207             :       ASSERT(index == num_modes,
     208             :              "Internal indexing error. "
     209             :                  << num_modes << " modes should have been calculated but "
     210             :                  << index << " modes were computed.");
     211             : 
     212             :       auto& worldtube_component = Parallel::get_parallel_component<
     213             :           Worldtube::WorldtubeSingleton<Metavariables>>(cache);
     214             :       Parallel::receive_data<Worldtube::Tags::SphericalHarmonicsInbox<Dim>>(
     215             :           worldtube_component, db::get<::Tags::TimeStepId>(box),
     216             :           std::make_pair(element_id, std::move(Ylm_coefs)));
     217             :     }
     218             :     return {Parallel::AlgorithmExecution::Continue, std::nullopt};
     219             :   }
     220             : };
     221             : }  // namespace CurvedScalarWave::Worldtube::Actions