Line data    Source code

       1            0 : // Distributed under the MIT License.
       2              : // See LICENSE.txt for details.
       3              : 
       4              : #pragma once
       5              : 
       6              : #include <array>
       7              : #include <cstddef>
       8              : #include <optional>
       9              : #include <unordered_map>
      10              : #include <utility>
      11              : #include <vector>
      12              : 
      13              : #include "Domain/Block.hpp"
      14              : #include "Domain/Creators/Tags/Domain.hpp"
      15              : #include "Domain/Creators/Tags/InitialExtents.hpp"
      16              : #include "Domain/Creators/Tags/InitialRefinementLevels.hpp"
      17              : #include "Domain/Domain.hpp"
      18              : #include "Domain/ElementMap.hpp"
      19              : #include "Domain/ExcisionSphere.hpp"
      20              : #include "Domain/InterfaceLogicalCoordinates.hpp"
      21              : #include "Domain/Structure/CreateInitialMesh.hpp"
      22              : #include "Domain/Structure/ElementId.hpp"
      23              : #include "Domain/Structure/InitialElementIds.hpp"
      24              : #include "Evolution/DiscontinuousGalerkin/Initialization/QuadratureTag.hpp"
      25              : #include "Evolution/Systems/CurvedScalarWave/Worldtube/Tags.hpp"
      26              : #include "NumericalAlgorithms/Spectral/Basis.hpp"
      27              : #include "NumericalAlgorithms/Spectral/Mesh.hpp"
      28              : #include "NumericalAlgorithms/Spectral/Quadrature.hpp"
      29              : #include "Utilities/Gsl.hpp"
      30              : #include "Utilities/TMPL.hpp"
      31              : 
      32              : namespace CurvedScalarWave::Worldtube::Initialization {
      33              : 
      34              : /*!
      35              :  * \brief Initializes a map of the grid coordinates centered on the worldtube of
      36              :  * all faces that abut the worldtube with corresponding ElementIds.
      37              :  *
      38              :  * \details The worldtube singleton computes an internal solution in the grid
      39              :  * frame and uses this map to evaluate compute it at the grid coordinates of
      40              :  * each element face abutting the worldtube each time step. This data is sent to
      41              :  * the corresponding element where it is used to apply pointwise boundary
      42              :  * conditions. The `ElementFacesGridCoordinates` holds a map of all the element
      43              :  * ids abutting the worldtube with the corresonding grid coordinates.
      44              :  *
      45              :  * \warning This currently assumes that initial domain remains the same and
      46              :  * there is no AMR. To support this, the worldtube could send the
      47              :  * coefficients of its internal solution to each element which can evaluate
      48              :  * it on their current grid in the boundary conditions.
      49              :  *
      50              :  *  DataBox changes:
      51              :  * - Adds: nothing
      52              :  * - Removes: nothing
      53              :  * - Modifies: Tags::ElementFacesCoordinatesMap<Dim>
      54              :  */
      55              : template <size_t Dim>
      56            1 : struct InitializeElementFacesGridCoordinates {
      57            0 :   using return_tags = tmpl::list<Tags::ElementFacesGridCoordinates<Dim>>;
      58            0 :   using simple_tags = tmpl::list<Tags::ElementFacesGridCoordinates<Dim>>;
      59            0 :   using compute_tags = tmpl::list<>;
      60            0 :   using simple_tags_from_options =
      61              :       tmpl::list<::domain::Tags::InitialExtents<Dim>,
      62              :                  ::domain::Tags::InitialRefinementLevels<Dim>,
      63              :                  evolution::dg::Tags::Quadrature>;
      64            0 :   using const_global_cache_tags = tmpl::list<>;
      65            0 :   using mutable_global_cache_tags = tmpl::list<>;
      66            0 :   using argument_tags = tmpl::flatten<
      67              :       tmpl::list<simple_tags_from_options, ::domain::Tags::Domain<Dim>,
      68              :                  Tags::ExcisionSphere<Dim>>>;
      69            0 :   static void apply(
      70              :       const gsl::not_null<std::unordered_map<
      71              :           ElementId<Dim>, tnsr::I<DataVector, Dim, Frame::Grid>>*>
      72              :           element_faces_grid_coords,
      73              :       const std::vector<std::array<size_t, Dim>>& initial_extents,
      74              :       const std::vector<std::array<size_t, Dim>>& initial_refinement,
      75              :       const Spectral::Quadrature& i1_quadrature, const Domain<Dim>& domain,
      76              :       const ::ExcisionSphere<Dim>& excision_sphere) {
      77              :     const Spectral::Basis i1_basis{Spectral::Basis::Legendre};
      78              :     const auto& blocks = domain.blocks();
      79              :     const auto& worldtube_grid_coords = excision_sphere.center();
      80              :     const auto& neighboring_blocks = excision_sphere.abutting_directions();
      81              :     for (const auto& [block_id, _] : neighboring_blocks) {
      82              :       const auto element_ids =
      83              :           initial_element_ids(block_id, initial_refinement.at(block_id));
      84              :       for (const auto element_id : element_ids) {
      85              :         const auto direction = excision_sphere.abutting_direction(element_id);
      86              :         if (direction.has_value()) {
      87              :           const auto& current_block = blocks.at(block_id);
      88              :           const auto mesh = ::domain::create_initial_mesh(
      89              :               initial_extents, current_block, element_id, i1_basis,
      90              :               i1_quadrature);
      91              :           const auto face_mesh = mesh.slice_away(direction.value().dimension());
      92              :           const ElementMap<Dim, Frame::Grid> element_map{
      93              :               element_id,
      94              :               current_block.is_time_dependent()
      95              :                   ? current_block.moving_mesh_logical_to_grid_map().get_clone()
      96              :                   : current_block.stationary_map().get_to_grid_frame()};
      97              :           const auto face_logical_coords =
      98              :               interface_logical_coordinates(face_mesh, direction.value());
      99              :           auto faces_grid_coords = element_map(face_logical_coords);
     100              :           for (size_t i = 0; i < 3; ++i) {
     101              :             faces_grid_coords.get(i) -= worldtube_grid_coords.get(i);
     102              :           }
     103              :           element_faces_grid_coords->operator[](element_id) =
     104              :               std::move(faces_grid_coords);
     105              :         }
     106              :       }
     107              :     }
     108              :   }
     109              : };
     110              : }  // namespace CurvedScalarWave::Worldtube::Initialization