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

Generated by: LCOV version 1.14