SpECTRE Documentation Coverage Report
Current view: top level - Elliptic/DiscontinuousGalerkin - Initialization.hpp Hit Total Coverage
Commit: d0fc80462417e83e5cddfa1b9901bb4a9b6af4d6 Lines: 3 22 13.6 %
Date: 2024-03-29 00:33:31
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 <type_traits>
       9             : #include <unordered_map>
      10             : #include <utility>
      11             : #include <vector>
      12             : 
      13             : #include "DataStructures/SliceVariables.hpp"
      14             : #include "DataStructures/Tensor/EagerMath/Magnitude.hpp"
      15             : #include "DataStructures/Tensor/EagerMath/RaiseOrLowerIndex.hpp"
      16             : #include "DataStructures/Tensor/Tensor.hpp"
      17             : #include "DataStructures/Variables.hpp"
      18             : #include "DataStructures/VariablesTag.hpp"
      19             : #include "Domain/Creators/Tags/Domain.hpp"
      20             : #include "Domain/Creators/Tags/InitialExtents.hpp"
      21             : #include "Domain/Creators/Tags/InitialRefinementLevels.hpp"
      22             : #include "Domain/Domain.hpp"
      23             : #include "Domain/ElementMap.hpp"
      24             : #include "Domain/FaceNormal.hpp"
      25             : #include "Domain/FunctionsOfTime/FunctionOfTime.hpp"
      26             : #include "Domain/FunctionsOfTime/Tags.hpp"
      27             : #include "Domain/InterfaceLogicalCoordinates.hpp"
      28             : #include "Domain/Structure/CreateInitialMesh.hpp"
      29             : #include "Domain/Structure/Direction.hpp"
      30             : #include "Domain/Structure/DirectionMap.hpp"
      31             : #include "Domain/Structure/DirectionalId.hpp"
      32             : #include "Domain/Structure/Element.hpp"
      33             : #include "Domain/Structure/IndexToSliceAt.hpp"
      34             : #include "Domain/Tags.hpp"
      35             : #include "Domain/Tags/FaceNormal.hpp"
      36             : #include "Domain/Tags/Faces.hpp"
      37             : #include "Domain/Tags/NeighborMesh.hpp"
      38             : #include "Domain/Tags/SurfaceJacobian.hpp"
      39             : #include "Elliptic/DiscontinuousGalerkin/Penalty.hpp"
      40             : #include "Elliptic/DiscontinuousGalerkin/Tags.hpp"
      41             : #include "NumericalAlgorithms/DiscontinuousGalerkin/MortarHelpers.hpp"
      42             : #include "NumericalAlgorithms/DiscontinuousGalerkin/ProjectToBoundary.hpp"
      43             : #include "NumericalAlgorithms/DiscontinuousGalerkin/Tags.hpp"
      44             : #include "NumericalAlgorithms/LinearOperators/PartialDerivatives.hpp"
      45             : #include "NumericalAlgorithms/Spectral/Mesh.hpp"
      46             : #include "NumericalAlgorithms/Spectral/Projection.hpp"
      47             : #include "NumericalAlgorithms/Spectral/Quadrature.hpp"
      48             : #include "Parallel/Tags/Metavariables.hpp"
      49             : #include "ParallelAlgorithms/Amr/Protocols/Projector.hpp"
      50             : #include "Utilities/CallWithDynamicType.hpp"
      51             : #include "Utilities/ErrorHandling/Assert.hpp"
      52             : #include "Utilities/Gsl.hpp"
      53             : #include "Utilities/ProtocolHelpers.hpp"
      54             : #include "Utilities/TMPL.hpp"
      55             : #include "Utilities/TaggedTuple.hpp"
      56             : 
      57             : namespace elliptic::dg {
      58             : 
      59             : namespace detail {
      60             : template <size_t Dim>
      61             : void initialize_coords_and_jacobians(
      62             :     gsl::not_null<tnsr::I<DataVector, Dim, Frame::ElementLogical>*>
      63             :         logical_coords,
      64             :     gsl::not_null<tnsr::I<DataVector, Dim, Frame::Inertial>*> inertial_coords,
      65             :     gsl::not_null<InverseJacobian<DataVector, Dim, Frame::ElementLogical,
      66             :                                   Frame::Inertial>*>
      67             :         inv_jacobian,
      68             :     gsl::not_null<Scalar<DataVector>*> det_inv_jacobian,
      69             :     gsl::not_null<Scalar<DataVector>*> det_jacobian,
      70             :     gsl::not_null<InverseJacobian<DataVector, Dim, Frame::ElementLogical,
      71             :                                   Frame::Inertial>*>
      72             :         det_times_inv_jacobian,
      73             :     const Mesh<Dim>& mesh, const ElementMap<Dim, Frame::Inertial>& element_map,
      74             :     const domain::FunctionsOfTimeMap& functions_of_time);
      75             : }  // namespace detail
      76             : 
      77             : /*!
      78             :  * \brief Initialize the background-independent geometry for the elliptic DG
      79             :  * operator
      80             :  *
      81             :  * ## Geometric aliasing
      82             :  *
      83             :  * The geometric quantities such as Jacobians are evaluated on the DG grid.
      84             :  * Since we know them analytically, we could also evaluate them on a
      85             :  * higher-order grid or with a stronger quadrature (Gauss instead of
      86             :  * Gauss-Lobatto) to combat geometric aliasing. See discussions in
      87             :  * \cite Vincent2019qpd and \cite Fischer2021voj .
      88             :  */
      89             : template <size_t Dim>
      90           1 : struct InitializeGeometry {
      91           0 :   using return_tags = tmpl::list<
      92             :       domain::Tags::Mesh<Dim>, domain::Tags::Element<Dim>,
      93             :       domain::Tags::NeighborMesh<Dim>, domain::Tags::ElementMap<Dim>,
      94             :       domain::Tags::Coordinates<Dim, Frame::ElementLogical>,
      95             :       domain::Tags::Coordinates<Dim, Frame::Inertial>,
      96             :       domain::Tags::InverseJacobian<Dim, Frame::ElementLogical,
      97             :                                     Frame::Inertial>,
      98             :       domain::Tags::DetInvJacobian<Frame::ElementLogical, Frame::Inertial>,
      99             :       domain::Tags::DetJacobian<Frame::ElementLogical, Frame::Inertial>,
     100             :       domain::Tags::DetTimesInvJacobian<Dim, Frame::ElementLogical,
     101             :                                         Frame::Inertial>>;
     102           0 :   using argument_tags =
     103             :       tmpl::list<domain::Tags::InitialExtents<Dim>,
     104             :                  domain::Tags::InitialRefinementLevels<Dim>,
     105             :                  domain::Tags::Domain<Dim>, domain::Tags::FunctionsOfTime,
     106             :                  elliptic::dg::Tags::Quadrature>;
     107           0 :   using volume_tags = argument_tags;
     108           0 :   static void apply(
     109             :       gsl::not_null<Mesh<Dim>*> mesh, gsl::not_null<Element<Dim>*> element,
     110             :       gsl::not_null<DirectionalIdMap<Dim, Mesh<Dim>>*> neighbor_meshes,
     111             :       gsl::not_null<ElementMap<Dim, Frame::Inertial>*> element_map,
     112             :       gsl::not_null<tnsr::I<DataVector, Dim, Frame::ElementLogical>*>
     113             :           logical_coords,
     114             :       gsl::not_null<tnsr::I<DataVector, Dim, Frame::Inertial>*> inertial_coords,
     115             :       gsl::not_null<InverseJacobian<DataVector, Dim, Frame::ElementLogical,
     116             :                                     Frame::Inertial>*>
     117             :           inv_jacobian,
     118             :       gsl::not_null<Scalar<DataVector>*> det_inv_jacobian,
     119             :       gsl::not_null<Scalar<DataVector>*> det_jacobian,
     120             :       gsl::not_null<InverseJacobian<DataVector, Dim, Frame::ElementLogical,
     121             :                                     Frame::Inertial>*>
     122             :           det_times_inv_jacobian,
     123             :       const std::vector<std::array<size_t, Dim>>& initial_extents,
     124             :       const std::vector<std::array<size_t, Dim>>& initial_refinement,
     125             :       const Domain<Dim>& domain,
     126             :       const domain::FunctionsOfTimeMap& functions_of_time,
     127             :       Spectral::Quadrature quadrature, const ElementId<Dim>& element_id);
     128             : };
     129             : 
     130             : template <size_t Dim>
     131           0 : struct ProjectGeometry : tt::ConformsTo<::amr::protocols::Projector> {
     132           0 :   using return_tags = tmpl::list<
     133             :       domain::Tags::ElementMap<Dim>,
     134             :       domain::Tags::Coordinates<Dim, Frame::ElementLogical>,
     135             :       domain::Tags::Coordinates<Dim, Frame::Inertial>,
     136             :       domain::Tags::InverseJacobian<Dim, Frame::ElementLogical,
     137             :                                     Frame::Inertial>,
     138             :       domain::Tags::DetInvJacobian<Frame::ElementLogical, Frame::Inertial>,
     139             :       domain::Tags::DetJacobian<Frame::ElementLogical, Frame::Inertial>,
     140             :       domain::Tags::DetTimesInvJacobian<Dim, Frame::ElementLogical,
     141             :                                         Frame::Inertial>>;
     142           0 :   using argument_tags =
     143             :       tmpl::list<domain::Tags::Mesh<Dim>, domain::Tags::Element<Dim>,
     144             :                  domain::Tags::Domain<Dim>, domain::Tags::FunctionsOfTime>;
     145           0 :   using volume_tags =
     146             :       tmpl::list<domain::Tags::Domain<Dim>, domain::Tags::FunctionsOfTime>;
     147             : 
     148             :   // p-refinement
     149           0 :   static void apply(
     150             :       const gsl::not_null<ElementMap<Dim, Frame::Inertial>*> element_map,
     151             :       const gsl::not_null<tnsr::I<DataVector, Dim, Frame::ElementLogical>*>
     152             :           logical_coords,
     153             :       const gsl::not_null<tnsr::I<DataVector, Dim, Frame::Inertial>*>
     154             :           inertial_coords,
     155             :       const gsl::not_null<InverseJacobian<
     156             :           DataVector, Dim, Frame::ElementLogical, Frame::Inertial>*>
     157             :           inv_jacobian,
     158             :       const gsl::not_null<Scalar<DataVector>*> det_inv_jacobian,
     159             :       const gsl::not_null<Scalar<DataVector>*> det_jacobian,
     160             :       const gsl::not_null<InverseJacobian<
     161             :           DataVector, Dim, Frame::ElementLogical, Frame::Inertial>*>
     162             :           det_times_inv_jacobian,
     163             :       const Mesh<Dim>& mesh, const Element<Dim>& /*element*/,
     164             :       const Domain<Dim>& /*domain*/,
     165             :       const domain::FunctionsOfTimeMap& functions_of_time,
     166             :       const std::pair<Mesh<Dim>, Element<Dim>>& old_mesh_and_element) {
     167             :     if (mesh == old_mesh_and_element.first) {
     168             :       // Neighbors may have changed, but this element hasn't. Nothing to do.
     169             :       return;
     170             :     }
     171             :     // The element map doesn't change with p-refinement
     172             :     detail::initialize_coords_and_jacobians(
     173             :         logical_coords, inertial_coords, inv_jacobian, det_inv_jacobian,
     174             :         det_jacobian, det_times_inv_jacobian, mesh, *element_map,
     175             :         functions_of_time);
     176             :   }
     177             : 
     178             :   // h-refinement
     179             :   template <typename... ParentOrChildrenItemsType>
     180           0 :   static void apply(
     181             :       const gsl::not_null<ElementMap<Dim, Frame::Inertial>*> element_map,
     182             :       const gsl::not_null<tnsr::I<DataVector, Dim, Frame::ElementLogical>*>
     183             :           logical_coords,
     184             :       const gsl::not_null<tnsr::I<DataVector, Dim, Frame::Inertial>*>
     185             :           inertial_coords,
     186             :       const gsl::not_null<InverseJacobian<
     187             :           DataVector, Dim, Frame::ElementLogical, Frame::Inertial>*>
     188             :           inv_jacobian,
     189             :       const gsl::not_null<Scalar<DataVector>*> det_inv_jacobian,
     190             :       const gsl::not_null<Scalar<DataVector>*> det_jacobian,
     191             :       const gsl::not_null<InverseJacobian<
     192             :           DataVector, Dim, Frame::ElementLogical, Frame::Inertial>*>
     193             :           det_times_inv_jacobian,
     194             :       const Mesh<Dim>& mesh, const Element<Dim>& element,
     195             :       const Domain<Dim>& domain,
     196             :       const domain::FunctionsOfTimeMap& functions_of_time,
     197             :       const ParentOrChildrenItemsType&... /*parent_or_children_items*/) {
     198             :     const auto& element_id = element.id();
     199             :     const auto& block = domain.blocks()[element_id.block_id()];
     200             :     *element_map = ElementMap<Dim, Frame::Inertial>{element_id, block};
     201             :     detail::initialize_coords_and_jacobians(
     202             :         logical_coords, inertial_coords, inv_jacobian, det_inv_jacobian,
     203             :         det_jacobian, det_times_inv_jacobian, mesh, *element_map,
     204             :         functions_of_time);
     205             :   }
     206             : };
     207             : 
     208             : namespace detail {
     209             : // Compute derivative of the Jacobian numerically
     210             : template <size_t Dim>
     211             : void deriv_unnormalized_face_normals_impl(
     212             :     gsl::not_null<DirectionMap<Dim, tnsr::ij<DataVector, Dim>>*>
     213             :         deriv_unnormalized_face_normals,
     214             :     const Mesh<Dim>& mesh, const Element<Dim>& element,
     215             :     const InverseJacobian<DataVector, Dim, Frame::ElementLogical,
     216             :                           Frame::Inertial>& inv_jacobian);
     217             : 
     218             : // Get element-logical coordinates of the mortar collocation points
     219             : template <size_t Dim>
     220             : tnsr::I<DataVector, Dim, Frame::ElementLogical> mortar_logical_coordinates(
     221             :     const Mesh<Dim - 1>& mortar_mesh,
     222             :     const ::dg::MortarSize<Dim - 1>& mortar_size,
     223             :     const Direction<Dim>& direction);
     224             : 
     225             : template <size_t Dim>
     226             : void mortar_jacobian(
     227             :     gsl::not_null<Scalar<DataVector>*> mortar_jacobian,
     228             :     gsl::not_null<Scalar<DataVector>*> perpendicular_element_size,
     229             :     const Mesh<Dim - 1>& mortar_mesh,
     230             :     const ::dg::MortarSize<Dim - 1>& mortar_size,
     231             :     const Direction<Dim>& direction,
     232             :     const tnsr::I<DataVector, Dim, Frame::ElementLogical>&
     233             :         mortar_logical_coords,
     234             :     const std::optional<tnsr::II<DataVector, Dim>>& inv_metric_on_mortar,
     235             :     const ElementMap<Dim, Frame::Inertial>& element_map,
     236             :     const domain::FunctionsOfTimeMap& functions_of_time);
     237             : }  // namespace detail
     238             : 
     239             : /// Initialize the geometry on faces and mortars for the elliptic DG operator
     240             : ///
     241             : /// To normalize face normals this function needs the inverse background metric.
     242             : /// Pass the tag representing the inverse background metric to the
     243             : /// `InvMetricTag` template parameter, and the tag representing the analytic
     244             : /// background from which it can be retrieved to the `BackgroundTag` template
     245             : /// parameter. Set `InvMetricTag` and `BackgroundTag` to `void` to normalize
     246             : /// face normals with the Euclidean magnitude.
     247             : ///
     248             : /// Mortar Jacobians are added only on nonconforming internal element
     249             : /// boundaries, i.e., when `Spectral::needs_projection()` is true.
     250             : ///
     251             : /// The `::Tags::deriv<domain::Tags::UnnormalizedFaceNormal<Dim>>` is only added
     252             : /// on external boundaries, for use by boundary conditions.
     253             : template <size_t Dim, typename InvMetricTag, typename BackgroundTag>
     254           1 : struct InitializeFacesAndMortars : tt::ConformsTo<::amr::protocols::Projector> {
     255           0 :   using return_tags = tmpl::append<
     256             :       domain::make_faces_tags<
     257             :           Dim,
     258             :           tmpl::list<domain::Tags::Direction<Dim>,
     259             :                      domain::Tags::Coordinates<Dim, Frame::Inertial>,
     260             :                      domain::Tags::FaceNormal<Dim>,
     261             :                      domain::Tags::FaceNormalVector<Dim>,
     262             :                      domain::Tags::UnnormalizedFaceNormalMagnitude<Dim>,
     263             :                      domain::Tags::DetSurfaceJacobian<Frame::ElementLogical,
     264             :                                                       Frame::Inertial>,
     265             :                      // This is the volume inverse Jacobian on the face grid
     266             :                      // points, multiplied by the determinant of the _face_
     267             :                      // Jacobian (the tag above)
     268             :                      domain::Tags::DetTimesInvJacobian<
     269             :                          Dim, Frame::ElementLogical, Frame::Inertial>,
     270             :                      // Possible optimization: The derivative of the face normal
     271             :                      // could be omitted for some systems, but its memory usage
     272             :                      // is probably insignificant since it's only added on
     273             :                      // external boundaries.
     274             :                      ::Tags::deriv<domain::Tags::UnnormalizedFaceNormal<Dim>,
     275             :                                    tmpl::size_t<Dim>, Frame::Inertial>>>,
     276             :       tmpl::list<::Tags::Mortars<domain::Tags::Mesh<Dim - 1>, Dim>,
     277             :                  ::Tags::Mortars<::Tags::MortarSize<Dim - 1>, Dim>,
     278             :                  ::Tags::Mortars<domain::Tags::DetSurfaceJacobian<
     279             :                                      Frame::ElementLogical, Frame::Inertial>,
     280             :                                  Dim>,
     281             :                  ::Tags::Mortars<elliptic::dg::Tags::PenaltyFactor, Dim>>>;
     282           0 :   using argument_tags = tmpl::append<
     283             :       tmpl::list<domain::Tags::Mesh<Dim>, domain::Tags::Element<Dim>,
     284             :                  domain::Tags::NeighborMesh<Dim>, domain::Tags::ElementMap<Dim>,
     285             :                  domain::Tags::InverseJacobian<Dim, Frame::ElementLogical,
     286             :                                                Frame::Inertial>,
     287             :                  domain::Tags::Domain<Dim>, domain::Tags::FunctionsOfTime,
     288             :                  elliptic::dg::Tags::PenaltyParameter>,
     289             :       tmpl::conditional_t<
     290             :           std::is_same_v<BackgroundTag, void>, tmpl::list<>,
     291             :           tmpl::list<BackgroundTag, Parallel::Tags::Metavariables>>>;
     292           0 :   using volume_tags = tmpl::append<
     293             :       tmpl::list<domain::Tags::Domain<Dim>, domain::Tags::FunctionsOfTime,
     294             :                  elliptic::dg::Tags::PenaltyParameter>,
     295             :       tmpl::conditional_t<
     296             :           std::is_same_v<BackgroundTag, void>, tmpl::list<>,
     297             :           tmpl::list<BackgroundTag, Parallel::Tags::Metavariables>>>;
     298             :   template <typename... AmrData>
     299           0 :   static void apply(
     300             :       const gsl::not_null<DirectionMap<Dim, Direction<Dim>>*> face_directions,
     301             :       const gsl::not_null<DirectionMap<Dim, tnsr::I<DataVector, Dim>>*>
     302             :           faces_inertial_coords,
     303             :       const gsl::not_null<DirectionMap<Dim, tnsr::i<DataVector, Dim>>*>
     304             :           face_normals,
     305             :       const gsl::not_null<DirectionMap<Dim, tnsr::I<DataVector, Dim>>*>
     306             :           face_normal_vectors,
     307             :       const gsl::not_null<DirectionMap<Dim, Scalar<DataVector>>*>
     308             :           face_normal_magnitudes,
     309             :       const gsl::not_null<DirectionMap<Dim, Scalar<DataVector>>*>
     310             :           face_jacobians,
     311             :       const gsl::not_null<DirectionMap<
     312             :           Dim, InverseJacobian<DataVector, Dim, Frame::ElementLogical,
     313             :                                Frame::Inertial>>*>
     314             :           face_jacobian_times_inv_jacobian,
     315             :       const gsl::not_null<DirectionMap<Dim, tnsr::ij<DataVector, Dim>>*>
     316             :           deriv_unnormalized_face_normals,
     317             :       const gsl::not_null<::dg::MortarMap<Dim, Mesh<Dim - 1>>*> mortar_meshes,
     318             :       const gsl::not_null<::dg::MortarMap<Dim, ::dg::MortarSize<Dim - 1>>*>
     319             :           mortar_sizes,
     320             :       const gsl::not_null<::dg::MortarMap<Dim, Scalar<DataVector>>*>
     321             :           mortar_jacobians,
     322             :       const gsl::not_null<::dg::MortarMap<Dim, Scalar<DataVector>>*>
     323             :           penalty_factors,
     324             :       const Mesh<Dim>& mesh, const Element<Dim>& element,
     325             :       const DirectionalIdMap<Dim, Mesh<Dim>>& neighbor_meshes,
     326             :       const ElementMap<Dim, Frame::Inertial>& element_map,
     327             :       const InverseJacobian<DataVector, Dim, Frame::ElementLogical,
     328             :                             Frame::Inertial>& inv_jacobian,
     329             :       const Domain<Dim>& domain,
     330             :       const domain::FunctionsOfTimeMap& functions_of_time,
     331             :       const double penalty_parameter, const AmrData&... amr_data) {
     332             :     apply(face_directions, faces_inertial_coords, face_normals,
     333             :           face_normal_vectors, face_normal_magnitudes, face_jacobians,
     334             :           face_jacobian_times_inv_jacobian, deriv_unnormalized_face_normals,
     335             :           mortar_meshes, mortar_sizes, mortar_jacobians, penalty_factors, mesh,
     336             :           element, neighbor_meshes, element_map, inv_jacobian, domain,
     337             :           functions_of_time, penalty_parameter, nullptr, nullptr, amr_data...);
     338             :   }
     339             :   template <typename Background, typename Metavariables, typename... AmrData>
     340           0 :   static void apply(
     341             :       const gsl::not_null<DirectionMap<Dim, Direction<Dim>>*> face_directions,
     342             :       const gsl::not_null<DirectionMap<Dim, tnsr::I<DataVector, Dim>>*>
     343             :           faces_inertial_coords,
     344             :       const gsl::not_null<DirectionMap<Dim, tnsr::i<DataVector, Dim>>*>
     345             :           face_normals,
     346             :       const gsl::not_null<DirectionMap<Dim, tnsr::I<DataVector, Dim>>*>
     347             :           face_normal_vectors,
     348             :       const gsl::not_null<DirectionMap<Dim, Scalar<DataVector>>*>
     349             :           face_normal_magnitudes,
     350             :       const gsl::not_null<DirectionMap<Dim, Scalar<DataVector>>*>
     351             :           face_jacobians,
     352             :       const gsl::not_null<DirectionMap<
     353             :           Dim, InverseJacobian<DataVector, Dim, Frame::ElementLogical,
     354             :                                Frame::Inertial>>*>
     355             :           face_jacobian_times_inv_jacobian,
     356             :       const gsl::not_null<DirectionMap<Dim, tnsr::ij<DataVector, Dim>>*>
     357             :           deriv_unnormalized_face_normals,
     358             :       const gsl::not_null<::dg::MortarMap<Dim, Mesh<Dim - 1>>*> mortar_meshes,
     359             :       const gsl::not_null<::dg::MortarMap<Dim, ::dg::MortarSize<Dim - 1>>*>
     360             :           mortar_sizes,
     361             :       const gsl::not_null<::dg::MortarMap<Dim, Scalar<DataVector>>*>
     362             :           mortar_jacobians,
     363             :       const gsl::not_null<::dg::MortarMap<Dim, Scalar<DataVector>>*>
     364             :           penalty_factors,
     365             :       const Mesh<Dim>& mesh, const Element<Dim>& element,
     366             :       const DirectionalIdMap<Dim, Mesh<Dim>>& neighbor_meshes,
     367             :       const ElementMap<Dim, Frame::Inertial>& element_map,
     368             :       const InverseJacobian<DataVector, Dim, Frame::ElementLogical,
     369             :                             Frame::Inertial>& inv_jacobian,
     370             :       const Domain<Dim>& domain,
     371             :       const domain::FunctionsOfTimeMap& functions_of_time,
     372             :       const double penalty_parameter, const Background& background,
     373             :       const Metavariables& /*meta*/, const AmrData&... /*amr_data*/) {
     374             :     static_assert(std::is_same_v<InvMetricTag, void> or
     375             :                       not(std::is_same_v<Background, std::nullptr_t>),
     376             :                   "Supply an analytic background from which the 'InvMetricTag' "
     377             :                   "can be retrieved");
     378             :     [[maybe_unused]] const auto get_inv_metric =
     379             :         [&background]([[maybe_unused]] const tnsr::I<DataVector, Dim>&
     380             :                           local_inertial_coords)
     381             :         -> std::optional<tnsr::II<DataVector, Dim>> {
     382             :       if constexpr (not std::is_same_v<InvMetricTag, void>) {
     383             :         using factory_classes = typename std::decay_t<
     384             :             Metavariables>::factory_creation::factory_classes;
     385             :         return call_with_dynamic_type<tnsr::II<DataVector, Dim>,
     386             :                                       tmpl::at<factory_classes, Background>>(
     387             :             &background, [&local_inertial_coords](const auto* const derived) {
     388             :               return get<InvMetricTag>(derived->variables(
     389             :                   local_inertial_coords, tmpl::list<InvMetricTag>{}));
     390             :             });
     391             :       } else {
     392             :         (void)background;
     393             :         return std::nullopt;
     394             :       }
     395             :     };
     396             :     ASSERT(std::equal(mesh.quadrature().begin() + 1, mesh.quadrature().end(),
     397             :                       mesh.quadrature().begin()),
     398             :            "This function is implemented assuming the quadrature is isotropic");
     399             :     // Faces
     400             :     for (const auto& direction : Direction<Dim>::all_directions()) {
     401             :       const auto face_mesh = mesh.slice_away(direction.dimension());
     402             :       (*face_directions)[direction] = direction;
     403             :       // Possible optimization: Not all systems need the coordinates on internal
     404             :       // faces.
     405             :       const auto face_logical_coords =
     406             :           interface_logical_coordinates(face_mesh, direction);
     407             :       auto& face_inertial_coords = (*faces_inertial_coords)[direction];
     408             :       face_inertial_coords =
     409             :           element_map(face_logical_coords, 0., functions_of_time);
     410             :       auto& face_normal = (*face_normals)[direction];
     411             :       auto& face_normal_vector = (*face_normal_vectors)[direction];
     412             :       auto& face_normal_magnitude = (*face_normal_magnitudes)[direction];
     413             :       // Buffer the inv Jacobian on the face here, then multiply by the face
     414             :       // Jacobian below
     415             :       auto& inv_jacobian_on_face =
     416             :           (*face_jacobian_times_inv_jacobian)[direction];
     417             :       inv_jacobian_on_face =
     418             :           element_map.inv_jacobian(face_logical_coords, 0., functions_of_time);
     419             :       unnormalized_face_normal(make_not_null(&face_normal), face_mesh,
     420             :                                inv_jacobian_on_face, direction);
     421             :       if constexpr (std::is_same_v<InvMetricTag, void>) {
     422             :         magnitude(make_not_null(&face_normal_magnitude), face_normal);
     423             :         for (size_t d = 0; d < Dim; ++d) {
     424             :           face_normal.get(d) /= get(face_normal_magnitude);
     425             :           face_normal_vector.get(d) = face_normal.get(d);
     426             :         }
     427             :       } else {
     428             :         const auto inv_metric_on_face = *get_inv_metric(face_inertial_coords);
     429             :         magnitude(make_not_null(&face_normal_magnitude), face_normal,
     430             :                   inv_metric_on_face);
     431             :         for (size_t d = 0; d < Dim; ++d) {
     432             :           face_normal.get(d) /= get(face_normal_magnitude);
     433             :         }
     434             :         raise_or_lower_index(make_not_null(&face_normal_vector), face_normal,
     435             :                              inv_metric_on_face);
     436             :       }
     437             :       auto& face_jacobian = (*face_jacobians)[direction];
     438             :       get(face_jacobian) =
     439             :           get(face_normal_magnitude) / get(determinant(inv_jacobian_on_face));
     440             :       for (auto& component : inv_jacobian_on_face) {
     441             :         component *= get(face_jacobian);
     442             :       }
     443             :     }
     444             :     // Compute the Jacobian derivative numerically, because our coordinate maps
     445             :     // currently don't provide it analytically.
     446             :     detail::deriv_unnormalized_face_normals_impl(
     447             :         deriv_unnormalized_face_normals, mesh, element, inv_jacobian);
     448             :     // Mortars (internal directions)
     449             :     mortar_meshes->clear();
     450             :     mortar_sizes->clear();
     451             :     mortar_jacobians->clear();
     452             :     penalty_factors->clear();
     453             :     const auto& element_id = element.id();
     454             :     for (const auto& [direction, neighbors] : element.neighbors()) {
     455             :       const auto face_mesh = mesh.slice_away(direction.dimension());
     456             :       const auto& orientation = neighbors.orientation();
     457             :       for (const auto& neighbor_id : neighbors) {
     458             :         const ::dg::MortarId<Dim> mortar_id{direction, neighbor_id};
     459             :         const auto& neighbor_mesh = neighbor_meshes.at(mortar_id);
     460             :         mortar_meshes->emplace(
     461             :             mortar_id, ::dg::mortar_mesh(
     462             :                            face_mesh, orientation.inverse_map()(neighbor_mesh)
     463             :                                           .slice_away(direction.dimension())));
     464             :         mortar_sizes->emplace(
     465             :             mortar_id, ::dg::mortar_size(element_id, neighbor_id,
     466             :                                          direction.dimension(), orientation));
     467             :         // Mortar Jacobian
     468             :         const auto& mortar_mesh = mortar_meshes->at(mortar_id);
     469             :         const auto& mortar_size = mortar_sizes->at(mortar_id);
     470             :         const auto mortar_logical_coords = detail::mortar_logical_coordinates(
     471             :             mortar_mesh, mortar_size, direction);
     472             :         const auto mortar_inertial_coords =
     473             :             element_map(mortar_logical_coords, 0., functions_of_time);
     474             :         Scalar<DataVector> perpendicular_element_size{};
     475             :         if (Spectral::needs_projection(face_mesh, mortar_mesh, mortar_size)) {
     476             :           auto& mortar_jacobian = (*mortar_jacobians)[mortar_id];
     477             :           detail::mortar_jacobian(make_not_null(&mortar_jacobian),
     478             :                                   make_not_null(&perpendicular_element_size),
     479             :                                   mortar_mesh, mortar_size, direction,
     480             :                                   mortar_logical_coords,
     481             :                                   get_inv_metric(mortar_inertial_coords),
     482             :                                   element_map, functions_of_time);
     483             :         } else {
     484             :           // Mortar is identical to face, and we have computed the face normal
     485             :           // magnitude already above
     486             :           get(perpendicular_element_size) =
     487             :               2. / get(face_normal_magnitudes->at(direction));
     488             :         }
     489             :         // Penalty factor
     490             :         // The penalty factor (like all quantities on mortars) must agree when
     491             :         // calculated on both sides of the mortar. So we switch perspective to
     492             :         // the neighbor here.
     493             :         const auto direction_in_neighbor = orientation(direction).opposite();
     494             :         const auto reoriented_mortar_mesh = ::dg::mortar_mesh(
     495             :             orientation(mesh).slice_away(direction_in_neighbor.dimension()),
     496             :             neighbor_mesh.slice_away(direction_in_neighbor.dimension()));
     497             :         const auto mortar_size_in_neighbor = ::dg::mortar_size(
     498             :             neighbor_id, element_id, direction_in_neighbor.dimension(),
     499             :             orientation.inverse_map());
     500             :         const auto mortar_logical_coords_in_neighbor =
     501             :             detail::mortar_logical_coordinates(reoriented_mortar_mesh,
     502             :                                                mortar_size_in_neighbor,
     503             :                                                direction_in_neighbor);
     504             :         const ElementMap<Dim, Frame::Inertial> neighbor_element_map{
     505             :             neighbor_id, domain.blocks()[neighbor_id.block_id()]};
     506             :         const auto reoriented_mortar_inertial_coords = neighbor_element_map(
     507             :             mortar_logical_coords_in_neighbor, 0., functions_of_time);
     508             :         Scalar<DataVector> buffer{};
     509             :         Scalar<DataVector> reoriented_neighbor_element_size{};
     510             :         detail::mortar_jacobian(
     511             :             make_not_null(&buffer),
     512             :             make_not_null(&reoriented_neighbor_element_size),
     513             :             reoriented_mortar_mesh, mortar_size_in_neighbor,
     514             :             direction_in_neighbor, mortar_logical_coords_in_neighbor,
     515             :             get_inv_metric(reoriented_mortar_inertial_coords),
     516             :             neighbor_element_map, functions_of_time);
     517             :         // Orient the result back to the perspective of this element
     518             :         const auto neighbor_element_size = orient_variables_on_slice(
     519             :             get(reoriented_neighbor_element_size),
     520             :             reoriented_mortar_mesh.extents(), direction_in_neighbor.dimension(),
     521             :             orientation.inverse_map());
     522             :         penalty_factors->emplace(
     523             :             mortar_id, elliptic::dg::penalty(
     524             :                            blaze::min(get(perpendicular_element_size),
     525             :                                       neighbor_element_size),
     526             :                            std::max(mesh.extents(direction.dimension()),
     527             :                                     neighbor_mesh.extents(
     528             :                                         direction_in_neighbor.dimension())),
     529             :                            penalty_parameter));
     530             :       }  // neighbors
     531             :     }    // internal directions
     532             :     // Mortars (external directions)
     533             :     for (const auto& direction : element.external_boundaries()) {
     534             :       const auto face_mesh = mesh.slice_away(direction.dimension());
     535             :       const auto mortar_id =
     536             :           DirectionalId<Dim>{direction, ElementId<Dim>::external_boundary_id()};
     537             :       mortar_meshes->emplace(mortar_id, face_mesh);
     538             :       mortar_sizes->emplace(mortar_id,
     539             :                             make_array<Dim - 1>(Spectral::MortarSize::Full));
     540             :       penalty_factors->emplace(
     541             :           mortar_id,
     542             :           elliptic::dg::penalty(2. / get(face_normal_magnitudes->at(direction)),
     543             :                                 mesh.extents(direction.dimension()),
     544             :                                 penalty_parameter));
     545             :     }  // external directions
     546             :   }
     547             : };
     548             : 
     549             : /// Initialize background quantities for the elliptic DG operator, possibly
     550             : /// including the metric necessary for normalizing face normals
     551             : template <size_t Dim, typename BackgroundFields, typename BackgroundTag>
     552           1 : struct InitializeBackground : tt::ConformsTo<::amr::protocols::Projector> {
     553           0 :   using return_tags =
     554             :       tmpl::list<::Tags::Variables<BackgroundFields>,
     555             :                  domain::Tags::Faces<Dim, ::Tags::Variables<BackgroundFields>>>;
     556           0 :   using argument_tags =
     557             :       tmpl::list<domain::Tags::Coordinates<Dim, Frame::Inertial>,
     558             :                  domain::Tags::Mesh<Dim>,
     559             :                  domain::Tags::InverseJacobian<Dim, Frame::ElementLogical,
     560             :                                                Frame::Inertial>,
     561             :                  BackgroundTag, Parallel::Tags::Metavariables>;
     562             : 
     563             :   template <typename BackgroundBase, typename Metavariables,
     564             :             typename... AmrData>
     565           0 :   static void apply(
     566             :       const gsl::not_null<Variables<BackgroundFields>*> background_fields,
     567             :       const gsl::not_null<DirectionMap<Dim, Variables<BackgroundFields>>*>
     568             :           face_background_fields,
     569             :       const tnsr::I<DataVector, Dim>& inertial_coords, const Mesh<Dim>& mesh,
     570             :       const InverseJacobian<DataVector, Dim, Frame::ElementLogical,
     571             :                             Frame::Inertial>& inv_jacobian,
     572             :       const BackgroundBase& background, const Metavariables& /*meta*/,
     573             :       const AmrData&... amr_data) {
     574             :     if constexpr (sizeof...(AmrData) == 1) {
     575             :       if constexpr (std::is_same_v<AmrData...,
     576             :                                    std::pair<Mesh<Dim>, Element<Dim>>>) {
     577             :         if (((mesh == amr_data.first) and ...)) {
     578             :           // This element hasn't changed during AMR. Nothing to do.
     579             :           return;
     580             :         }
     581             :       }
     582             :     }
     583             : 
     584             :     // Background fields in the volume
     585             :     using factory_classes =
     586             :         typename std::decay_t<Metavariables>::factory_creation::factory_classes;
     587             :     *background_fields =
     588             :         call_with_dynamic_type<Variables<BackgroundFields>,
     589             :                                tmpl::at<factory_classes, BackgroundBase>>(
     590             :             &background, [&inertial_coords, &mesh,
     591             :                           &inv_jacobian](const auto* const derived) {
     592             :               return variables_from_tagged_tuple(derived->variables(
     593             :                   inertial_coords, mesh, inv_jacobian, BackgroundFields{}));
     594             :             });
     595             :     // Background fields on faces
     596             :     for (const auto& direction : Direction<Dim>::all_directions()) {
     597             :       // Possible optimization: Only the background fields in the
     598             :       // System::fluxes_computer::argument_tags are needed on internal faces.
     599             :       // On Gauss grids we could evaluate the background directly on the faces
     600             :       // instead of projecting the data. However, we need to take derivatives of
     601             :       // the background fields, so we could evaluate them on a Gauss-Lobatto
     602             :       // grid in the volume. We could even evaluate the background fields on a
     603             :       // higher-order grid and project down to get more accurate derivatives if
     604             :       // needed.
     605             :       (*face_background_fields)[direction].initialize(
     606             :           mesh.slice_away(direction.dimension()).number_of_grid_points());
     607             :       ::dg::project_contiguous_data_to_boundary(
     608             :           make_not_null(&(*face_background_fields)[direction]),
     609             :           *background_fields, mesh, direction);
     610             :     }
     611             :   }
     612             : };
     613             : 
     614             : }  // namespace elliptic::dg

Generated by: LCOV version 1.14