SpECTRE Documentation Coverage Report
Current view: top level - Elliptic/DiscontinuousGalerkin/SubdomainOperator - InitializeSubdomain.hpp Hit Total Coverage
Commit: 802361e1b4aed23d7a052290519255a28e9b74f4 Lines: 2 21 9.5 %
Date: 2021-11-24 08:38:40
Legend: Lines: hit not hit

          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 <tuple>
       9             : #include <type_traits>
      10             : #include <unordered_map>
      11             : #include <utility>
      12             : #include <vector>
      13             : 
      14             : #include "DataStructures/DataBox/DataBox.hpp"
      15             : #include "DataStructures/DataBox/PrefixHelpers.hpp"
      16             : #include "DataStructures/DataVector.hpp"
      17             : #include "DataStructures/SliceVariables.hpp"
      18             : #include "DataStructures/Tensor/EagerMath/Magnitude.hpp"
      19             : #include "DataStructures/Tensor/Tensor.hpp"
      20             : #include "DataStructures/Variables.hpp"
      21             : #include "DataStructures/VariablesTag.hpp"
      22             : #include "Domain/ElementMap.hpp"
      23             : #include "Domain/FaceNormal.hpp"
      24             : #include "Domain/LogicalCoordinates.hpp"
      25             : #include "Domain/Structure/CreateInitialMesh.hpp"
      26             : #include "Domain/Structure/Direction.hpp"
      27             : #include "Domain/Structure/DirectionMap.hpp"
      28             : #include "Domain/Structure/Element.hpp"
      29             : #include "Domain/Structure/ElementId.hpp"
      30             : #include "Domain/Structure/IndexToSliceAt.hpp"
      31             : #include "Domain/Tags.hpp"
      32             : #include "Domain/Tags/FaceNormal.hpp"
      33             : #include "Domain/Tags/Faces.hpp"
      34             : #include "Elliptic/DiscontinuousGalerkin/Initialization.hpp"
      35             : #include "Elliptic/DiscontinuousGalerkin/SubdomainOperator/Tags.hpp"
      36             : #include "Elliptic/Utilities/ApplyAt.hpp"
      37             : #include "NumericalAlgorithms/DiscontinuousGalerkin/MortarHelpers.hpp"
      38             : #include "NumericalAlgorithms/DiscontinuousGalerkin/Tags.hpp"
      39             : #include "NumericalAlgorithms/Spectral/Mesh.hpp"
      40             : #include "NumericalAlgorithms/Spectral/Spectral.hpp"
      41             : #include "ParallelAlgorithms/LinearSolver/Schwarz/OverlapHelpers.hpp"
      42             : #include "ParallelAlgorithms/LinearSolver/Schwarz/Tags.hpp"
      43             : #include "Utilities/Gsl.hpp"
      44             : #include "Utilities/MakeArray.hpp"
      45             : #include "Utilities/Requires.hpp"
      46             : #include "Utilities/TMPL.hpp"
      47             : 
      48             : /// \cond
      49             : namespace Parallel {
      50             : template <typename Metavariables>
      51             : struct GlobalCache;
      52             : }  // namespace Parallel
      53             : namespace tuples {
      54             : template <typename... Tags>
      55             : struct TaggedTuple;
      56             : }  // namespace tuples
      57             : /// \endcond
      58             : 
      59             : /// Actions related to the DG subdomain operator
      60           1 : namespace elliptic::dg::subdomain_operator::Actions {
      61             : 
      62             : namespace detail {
      63             : // Initialize the geometry of a neighbor into which an overlap extends
      64             : template <size_t Dim>
      65             : struct InitializeOverlapGeometry {
      66             :   using return_tags =
      67             :       tmpl::list<elliptic::dg::subdomain_operator::Tags::ExtrudingExtent,
      68             :                  elliptic::dg::subdomain_operator::Tags::NeighborMortars<
      69             :                      domain::Tags::Mesh<Dim>, Dim>,
      70             :                  elliptic::dg::subdomain_operator::Tags::NeighborMortars<
      71             :                      domain::Tags::UnnormalizedFaceNormalMagnitude<Dim>, Dim>,
      72             :                  elliptic::dg::subdomain_operator::Tags::NeighborMortars<
      73             :                      domain::Tags::Mesh<Dim - 1>, Dim>,
      74             :                  elliptic::dg::subdomain_operator::Tags::NeighborMortars<
      75             :                      ::Tags::MortarSize<Dim - 1>, Dim>>;
      76             :   using argument_tags =
      77             :       tmpl::list<domain::Tags::Element<Dim>, domain::Tags::Mesh<Dim>>;
      78             :   void operator()(
      79             :       gsl::not_null<size_t*> extruding_extent,
      80             :       gsl::not_null<::dg::MortarMap<Dim, Mesh<Dim>>*> neighbor_meshes,
      81             :       gsl::not_null<::dg::MortarMap<Dim, Scalar<DataVector>>*>
      82             :           neighbor_face_normal_magnitudes,
      83             :       gsl::not_null<::dg::MortarMap<Dim, Mesh<Dim - 1>>*>
      84             :           neighbor_mortar_meshes,
      85             :       gsl::not_null<::dg::MortarMap<Dim, ::dg::MortarSize<Dim - 1>>*>
      86             :           neighbor_mortar_sizes,
      87             :       const Element<Dim>& element, const Mesh<Dim>& mesh,
      88             :       const std::vector<std::array<size_t, Dim>>& initial_extents,
      89             :       const ElementId<Dim>& element_id, const Direction<Dim>& overlap_direction,
      90             :       const size_t max_overlap) const;
      91             : };
      92             : }  // namespace detail
      93             : 
      94             : /*!
      95             :  * \brief Initialize the geometry for the DG subdomain operator
      96             :  *
      97             :  * Initializes tags that define the geometry of overlap regions with neighboring
      98             :  * elements. The data needs to be updated if the geometry of neighboring
      99             :  * elements changes.
     100             :  *
     101             :  * Note that the geometry depends on the system and on the choice of background
     102             :  * through the normalization of face normals, which involves a metric.
     103             :  *
     104             :  * DataBox:
     105             :  * - Uses:
     106             :  *   - `BackgroundTag`
     107             :  *   - `domain::Tags::Element<Dim>`
     108             :  *   - `domain::Tags::InitialExtents<Dim>`
     109             :  *   - `domain::Tags::InitialRefinementLevels<Dim>`
     110             :  *   - `domain::Tags::Domain<Dim>`
     111             :  *   - `LinearSolver::Schwarz::Tags::MaxOverlap<OptionsGroup>`
     112             :  * - Adds: Many tags prefixed with `LinearSolver::Schwarz::Tags::Overlaps`. See
     113             :  *   `elliptic::dg::Actions::InitializeDomain` and
     114             :  *   `elliptic::dg::Actions::initialize_operator` for a complete list.
     115             :  */
     116             : template <typename System, typename BackgroundTag, typename OptionsGroup>
     117           1 : struct InitializeSubdomain {
     118             :  private:
     119           0 :   static constexpr size_t Dim = System::volume_dim;
     120           0 :   static constexpr bool is_curved =
     121             :       not std::is_same_v<typename System::inv_metric_tag, void>;
     122           0 :   static constexpr bool has_background_fields =
     123             :       not std::is_same_v<typename System::background_fields, tmpl::list<>>;
     124             : 
     125           0 :   using InitializeGeometry = elliptic::dg::InitializeGeometry<Dim>;
     126           0 :   using InitializeOverlapGeometry = detail::InitializeOverlapGeometry<Dim>;
     127           0 :   using InitializeFacesAndMortars =
     128             :       elliptic::dg::InitializeFacesAndMortars<Dim>;
     129           0 :   using NormalizeFaceNormal =
     130             :       elliptic::dg::NormalizeFaceNormal<Dim, typename System::inv_metric_tag>;
     131             : 
     132             :   template <typename Tag>
     133           0 :   using overlaps_tag =
     134             :       LinearSolver::Schwarz::Tags::Overlaps<Tag, Dim, OptionsGroup>;
     135             :   // Only slice those background fields to internal boundaries that are
     136             :   // necessary for the DG operator, i.e. the background fields in the
     137             :   // System::fluxes_computer::argument_tags
     138           0 :   using fluxes_non_background_args =
     139             :       tmpl::list_difference<typename System::fluxes_computer::argument_tags,
     140             :                             typename System::background_fields>;
     141           0 :   using background_fields_internal =
     142             :       tmpl::list_difference<typename System::fluxes_computer::argument_tags,
     143             :                             fluxes_non_background_args>;
     144             :   // Slice all background fields to external boundaries for use in boundary
     145             :   // conditions
     146           0 :   using background_fields_external = typename System::background_fields;
     147             : 
     148             :  public:
     149           0 :   using initialization_tags =
     150             :       tmpl::list<domain::Tags::InitialExtents<Dim>,
     151             :                  domain::Tags::InitialRefinementLevels<Dim>>;
     152           0 :   using const_global_cache_tags =
     153             :       tmpl::list<LinearSolver::Schwarz::Tags::MaxOverlap<OptionsGroup>>;
     154           0 :   using simple_tags = db::wrap_tags_in<
     155             :       overlaps_tag,
     156             :       tmpl::append<
     157             :           typename InitializeGeometry::return_tags,
     158             :           typename InitializeFacesAndMortars::return_tags,
     159             :           typename InitializeOverlapGeometry::return_tags,
     160             :           tmpl::conditional_t<
     161             :               has_background_fields,
     162             :               tmpl::list<::Tags::Variables<typename System::background_fields>>,
     163             :               tmpl::list<>>,
     164             :           domain::make_faces_tags<Dim, typename System::background_fields>>>;
     165           0 :   using compute_tags = tmpl::list<>;
     166             : 
     167             :   template <typename DataBox, typename... InboxTags, typename Metavariables,
     168             :             typename ActionList, typename ParallelComponent>
     169           0 :   static std::tuple<DataBox&&> apply(
     170             :       DataBox& box, const tuples::TaggedTuple<InboxTags...>& /*inboxes*/,
     171             :       const Parallel::GlobalCache<Metavariables>& /*cache*/,
     172             :       const ElementId<Dim>& /*element_id*/, const ActionList /*meta*/,
     173             :       const ParallelComponent* const /*meta*/) {
     174             :     if constexpr (tmpl::all<initialization_tags,
     175             :                             tmpl::bind<db::tag_is_retrievable, tmpl::_1,
     176             :                                        tmpl::pin<DataBox>>>::value) {
     177             :       const auto& element = db::get<domain::Tags::Element<Dim>>(box);
     178             :       for (const auto& [direction, neighbors] : element.neighbors()) {
     179             :         const auto& orientation = neighbors.orientation();
     180             :         const auto direction_from_neighbor = orientation(direction.opposite());
     181             :         for (const auto& neighbor_id : neighbors) {
     182             :           const LinearSolver::Schwarz::OverlapId<Dim> overlap_id{direction,
     183             :                                                                  neighbor_id};
     184             :           // Initialize background-agnostic geometry on overlaps
     185             :           elliptic::util::mutate_apply_at<
     186             :               db::wrap_tags_in<overlaps_tag,
     187             :                                typename InitializeGeometry::return_tags>,
     188             :               typename InitializeGeometry::argument_tags,
     189             :               typename InitializeGeometry::argument_tags>(
     190             :               InitializeGeometry{}, make_not_null(&box), overlap_id,
     191             :               neighbor_id);
     192             :           // Initialize faces and mortars on overlaps
     193             :           elliptic::util::mutate_apply_at<
     194             :               db::wrap_tags_in<overlaps_tag,
     195             :                                typename InitializeFacesAndMortars::return_tags>,
     196             :               db::wrap_tags_in<
     197             :                   overlaps_tag,
     198             :                   typename InitializeFacesAndMortars::argument_tags>,
     199             :               tmpl::list<>>(InitializeFacesAndMortars{}, make_not_null(&box),
     200             :                             overlap_id,
     201             :                             db::get<domain::Tags::InitialExtents<Dim>>(box));
     202             :           // Initialize subdomain-specific tags on overlaps
     203             :           elliptic::util::mutate_apply_at<
     204             :               db::wrap_tags_in<overlaps_tag,
     205             :                                typename InitializeOverlapGeometry::return_tags>,
     206             :               db::wrap_tags_in<
     207             :                   overlaps_tag,
     208             :                   typename InitializeOverlapGeometry::argument_tags>,
     209             :               tmpl::list<>>(
     210             :               InitializeOverlapGeometry{}, make_not_null(&box), overlap_id,
     211             :               db::get<domain::Tags::InitialExtents<Dim>>(box), neighbor_id,
     212             :               direction_from_neighbor,
     213             :               db::get<LinearSolver::Schwarz::Tags::MaxOverlap<OptionsGroup>>(
     214             :                   box));
     215             :           // Background fields
     216             :           if constexpr (has_background_fields) {
     217             :             initialize_background_fields(make_not_null(&box), overlap_id);
     218             :           }
     219             :           // Normalize face normals
     220             :           normalize_face_normals(make_not_null(&box), overlap_id);
     221             :         }  // neighbors in direction
     222             :       }    // directions
     223             :     } else {
     224             :       ERROR(
     225             :           "Dependencies not fulfilled. Did you forget to terminate the phase "
     226             :           "after removing options?");
     227             :     }
     228             :     return {std::move(box)};
     229             :   }
     230             : 
     231             :  private:
     232             :   template <typename DbTagsList>
     233           0 :   static void initialize_background_fields(
     234             :       const gsl::not_null<db::DataBox<DbTagsList>*> box,
     235             :       const LinearSolver::Schwarz::OverlapId<Dim>& overlap_id) {
     236             :     const auto& background = db::get<BackgroundTag>(*box);
     237             :     DirectionMap<Dim, Variables<typename System::background_fields>>
     238             :         face_background_fields{};
     239             :     elliptic::util::mutate_apply_at<
     240             :         tmpl::list<overlaps_tag<
     241             :             ::Tags::Variables<typename System::background_fields>>>,
     242             :         db::wrap_tags_in<
     243             :             overlaps_tag,
     244             :             tmpl::list<domain::Tags::Coordinates<Dim, Frame::Inertial>,
     245             :                        domain::Tags::Mesh<Dim>,
     246             :                        domain::Tags::InverseJacobian<Dim, Frame::ElementLogical,
     247             :                                                      Frame::Inertial>>>,
     248             :         tmpl::list<>>(
     249             :         [&background, &face_background_fields](
     250             :             const gsl::not_null<Variables<typename System::background_fields>*>
     251             :                 background_fields,
     252             :             const tnsr::I<DataVector, Dim>& inertial_coords,
     253             :             const Mesh<Dim>& mesh,
     254             :             const InverseJacobian<DataVector, Dim, Frame::ElementLogical,
     255             :                                   Frame::Inertial>& inv_jacobian) {
     256             :           *background_fields = variables_from_tagged_tuple(
     257             :               background.variables(inertial_coords, mesh, inv_jacobian,
     258             :                                    typename System::background_fields{}));
     259             :           for (const auto& direction : Direction<Dim>::all_directions()) {
     260             :             // Slice the background fields to the face instead of evaluating
     261             :             // them on the face coords to avoid re-computing them, and because
     262             :             // this is also what the DG operator currently does. The result is
     263             :             // the same on Gauss-Lobatto grids, but may need adjusting when
     264             :             // adding support for Gauss grids.
     265             :             data_on_slice(make_not_null(&face_background_fields[direction]),
     266             :                           *background_fields, mesh.extents(),
     267             :                           direction.dimension(),
     268             :                           index_to_slice_at(mesh.extents(), direction));
     269             :           }
     270             :         },
     271             :         box, overlap_id);
     272             :     // Move face background fields into DataBox
     273             :     const auto mutate_assign_face_background_field =
     274             :         [&box, &overlap_id, &face_background_fields](
     275             :             auto tag_v, const Direction<Dim>& direction) {
     276             :           using tag = tmpl::type_from<std::decay_t<decltype(tag_v)>>;
     277             :           db::mutate<overlaps_tag<domain::Tags::Faces<Dim, tag>>>(
     278             :               box, [&face_background_fields, &overlap_id,
     279             :                     &direction](const auto stored_value) {
     280             :                 (*stored_value)[overlap_id][direction] =
     281             :                     get<tag>(face_background_fields.at(direction));
     282             :               });
     283             :         };
     284             :     const auto& element =
     285             :         db::get<overlaps_tag<domain::Tags::Element<Dim>>>(*box).at(overlap_id);
     286             :     for (const auto& direction : element.internal_boundaries()) {
     287             :       tmpl::for_each<background_fields_internal>(
     288             :           [&mutate_assign_face_background_field, &direction](auto tag_v) {
     289             :             mutate_assign_face_background_field(tag_v, direction);
     290             :           });
     291             :     }
     292             :     for (const auto& direction : element.external_boundaries()) {
     293             :       tmpl::for_each<background_fields_external>(
     294             :           [&mutate_assign_face_background_field, &direction](auto tag_v) {
     295             :             mutate_assign_face_background_field(tag_v, direction);
     296             :           });
     297             :     }
     298             :   }
     299             : 
     300             :   template <typename DbTagsList>
     301           0 :   static void normalize_face_normals(
     302             :       const gsl::not_null<db::DataBox<DbTagsList>*> box,
     303             :       const LinearSolver::Schwarz::OverlapId<Dim>& overlap_id) {
     304             :     // Faces of the overlapped element (internal and external)
     305             :     const auto& element =
     306             :         db::get<overlaps_tag<domain::Tags::Element<Dim>>>(*box).at(overlap_id);
     307             :     for (const auto& direction : Direction<Dim>::all_directions()) {
     308             :       elliptic::util::mutate_apply_at<
     309             :           db::wrap_tags_in<overlaps_tag,
     310             :                            domain::make_faces_tags<
     311             :                                Dim, typename NormalizeFaceNormal::return_tags>>,
     312             :           db::wrap_tags_in<
     313             :               overlaps_tag,
     314             :               domain::make_faces_tags<
     315             :                   Dim, typename NormalizeFaceNormal::argument_tags>>,
     316             :           tmpl::list<>>(NormalizeFaceNormal{}, box,
     317             :                         std::make_tuple(overlap_id, direction));
     318             :     }
     319             :     // Faces on the other side of the overlapped element's mortars
     320             :     const auto& domain = db::get<domain::Tags::Domain<Dim>>(*box);
     321             :     const auto& neighbor_meshes = db::get<overlaps_tag<
     322             :         elliptic::dg::subdomain_operator::Tags::NeighborMortars<
     323             :             domain::Tags::Mesh<Dim>, Dim>>>(*box)
     324             :                                       .at(overlap_id);
     325             :     for (const auto& [direction, neighbors] : element.neighbors()) {
     326             :       const auto& orientation = neighbors.orientation();
     327             :       const auto direction_from_neighbor = orientation(direction.opposite());
     328             :       for (const auto& neighbor_id : neighbors) {
     329             :         const ::dg::MortarId<Dim> mortar_id{direction, neighbor_id};
     330             :         const auto neighbor_face_mesh =
     331             :             neighbor_meshes.at(mortar_id).slice_away(
     332             :                 direction_from_neighbor.dimension());
     333             :         const auto& neighbor_block = domain.blocks()[neighbor_id.block_id()];
     334             :         ElementMap<Dim, Frame::Inertial> neighbor_element_map{
     335             :             neighbor_id, neighbor_block.stationary_map().get_clone()};
     336             :         const auto neighbor_face_normal = unnormalized_face_normal(
     337             :             neighbor_face_mesh, neighbor_element_map, direction_from_neighbor);
     338             :         using neighbor_face_normal_magnitudes_tag = overlaps_tag<
     339             :             elliptic::dg::subdomain_operator::Tags::NeighborMortars<
     340             :                 domain::Tags::UnnormalizedFaceNormalMagnitude<Dim>, Dim>>;
     341             :         if constexpr (is_curved) {
     342             :           const auto& background = db::get<BackgroundTag>(*box);
     343             :           const auto neighbor_face_inertial_coords =
     344             :               neighbor_element_map(interface_logical_coordinates(
     345             :                   neighbor_face_mesh, direction_from_neighbor));
     346             :           const auto inv_metric_on_face =
     347             :               get<typename System::inv_metric_tag>(background.variables(
     348             :                   neighbor_face_inertial_coords,
     349             :                   tmpl::list<typename System::inv_metric_tag>{}));
     350             :           elliptic::util::mutate_apply_at<
     351             :               tmpl::list<neighbor_face_normal_magnitudes_tag>, tmpl::list<>,
     352             :               tmpl::list<>>(
     353             :               [&neighbor_face_normal,
     354             :                &inv_metric_on_face](const auto neighbor_face_normal_magnitude) {
     355             :                 magnitude(neighbor_face_normal_magnitude, neighbor_face_normal,
     356             :                           inv_metric_on_face);
     357             :               },
     358             :               box, std::make_tuple(overlap_id, mortar_id));
     359             :         } else {
     360             :           elliptic::util::mutate_apply_at<
     361             :               tmpl::list<neighbor_face_normal_magnitudes_tag>, tmpl::list<>,
     362             :               tmpl::list<>>(
     363             :               [&neighbor_face_normal](
     364             :                   const auto neighbor_face_normal_magnitude) {
     365             :                 magnitude(neighbor_face_normal_magnitude, neighbor_face_normal);
     366             :               },
     367             :               box, std::make_tuple(overlap_id, mortar_id));
     368             :         }
     369             :       }  // neighbors
     370             :     }    // internal directions
     371             :   }
     372             : };
     373             : 
     374             : }  // namespace elliptic::dg::subdomain_operator::Actions

Generated by: LCOV version 1.14