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 <limits>
       9             : #include <memory>
      10             : 
      11             : #include "DataStructures/DataBox/PrefixHelpers.hpp"
      12             : #include "DataStructures/DataBox/Prefixes.hpp"
      13             : #include "DataStructures/DataVector.hpp"
      14             : #include "DataStructures/Tensor/TypeAliases.hpp"
      15             : #include "Domain/Structure/DirectionalIdMap.hpp"
      16             : #include "Domain/Tags.hpp"
      17             : #include "Evolution/DgSubcell/GhostData.hpp"
      18             : #include "Evolution/DgSubcell/Tags/GhostDataForReconstruction.hpp"
      19             : #include "Evolution/DgSubcell/Tags/Mesh.hpp"
      20             : #include "Evolution/Systems/NewtonianEuler/FiniteDifference/Reconstructor.hpp"
      21             : #include "Evolution/Systems/NewtonianEuler/Tags.hpp"
      22             : #include "PointwiseFunctions/Hydro/EquationsOfState/EquationOfState.hpp"
      23             : #include "PointwiseFunctions/Hydro/Tags.hpp"
      24             : 
      25           1 : namespace NewtonianEuler::fd {
      26             : /*!
      27             :  * \brief Adaptive-order WENO reconstruction hybridizing orders 5 and 3. See
      28             :  * ::fd::reconstruction::aoweno_53() for details.
      29             :  */
      30             : template <size_t Dim>
      31           1 : class AoWeno53Prim : public Reconstructor<Dim> {
      32             :  private:
      33             :   // Conservative vars tags
      34           0 :   using MassDensityCons = NewtonianEuler::Tags::MassDensityCons;
      35           0 :   using EnergyDensity = NewtonianEuler::Tags::EnergyDensity;
      36           0 :   using MomentumDensity = NewtonianEuler::Tags::MomentumDensity<Dim>;
      37             : 
      38             :   // Primitive vars tags
      39           0 :   using MassDensity = hydro::Tags::RestMassDensity<DataVector>;
      40           0 :   using Velocity = hydro::Tags::SpatialVelocity<DataVector, Dim>;
      41           0 :   using SpecificInternalEnergy =
      42             :       hydro::Tags::SpecificInternalEnergy<DataVector>;
      43           0 :   using Pressure = hydro::Tags::Pressure<DataVector>;
      44             : 
      45           0 :   using prims_tags =
      46             :       tmpl::list<MassDensity, Velocity, SpecificInternalEnergy, Pressure>;
      47           0 :   using cons_tags = tmpl::list<MassDensityCons, MomentumDensity, EnergyDensity>;
      48           0 :   using flux_tags = db::wrap_tags_in<::Tags::Flux, cons_tags, tmpl::size_t<Dim>,
      49             :                                      Frame::Inertial>;
      50           0 :   using prim_tags_for_reconstruction =
      51             :       tmpl::list<MassDensity, Velocity, Pressure>;
      52             : 
      53             :  public:
      54           0 :   struct GammaHi {
      55           0 :     using type = double;
      56           0 :     static constexpr Options::String help = {
      57             :         "The linear weight for the 5th-order stencil."};
      58             :   };
      59           0 :   struct GammaLo {
      60           0 :     using type = double;
      61           0 :     static constexpr Options::String help = {
      62             :         "The linear weight for the central 3rd-order stencil."};
      63             :   };
      64           0 :   struct Epsilon {
      65           0 :     using type = double;
      66           0 :     static constexpr Options::String help = {
      67             :         "The parameter added to the oscillation indicators to avoid division "
      68             :         "by zero"};
      69             :   };
      70           0 :   struct NonlinearWeightExponent {
      71           0 :     using type = size_t;
      72           0 :     static constexpr Options::String help = {
      73             :         "The exponent q to which the oscillation indicators are raised"};
      74             :   };
      75             : 
      76           0 :   using options =
      77             :       tmpl::list<GammaHi, GammaLo, Epsilon, NonlinearWeightExponent>;
      78           0 :   static constexpr Options::String help{
      79             :       "Adaptive-order WENO reconstruction hybridizing orders 5 and 3 using "
      80             :       "primitive variables."};
      81             : 
      82           0 :   AoWeno53Prim() = default;
      83           0 :   AoWeno53Prim(AoWeno53Prim&&) = default;
      84           0 :   AoWeno53Prim& operator=(AoWeno53Prim&&) = default;
      85           0 :   AoWeno53Prim(const AoWeno53Prim&) = default;
      86           0 :   AoWeno53Prim& operator=(const AoWeno53Prim&) = default;
      87           0 :   ~AoWeno53Prim() override = default;
      88             : 
      89           0 :   AoWeno53Prim(double gamma_hi, double gamma_lo, double epsilon,
      90             :                size_t nonlinear_weight_exponent);
      91             : 
      92           0 :   explicit AoWeno53Prim(CkMigrateMessage* msg);
      93             : 
      94           0 :   WRAPPED_PUPable_decl_base_template(Reconstructor<Dim>, AoWeno53Prim);
      95             : 
      96           0 :   auto get_clone() const -> std::unique_ptr<Reconstructor<Dim>> override;
      97             : 
      98           0 :   void pup(PUP::er& p) override;
      99             : 
     100           0 :   size_t ghost_zone_size() const override { return 3; }
     101             : 
     102           0 :   using reconstruction_argument_tags =
     103             :       tmpl::list<::Tags::Variables<prims_tags>,
     104             :                  hydro::Tags::EquationOfState<false, 2>,
     105             :                  domain::Tags::Element<Dim>,
     106             :                  evolution::dg::subcell::Tags::GhostDataForReconstruction<Dim>,
     107             :                  evolution::dg::subcell::Tags::Mesh<Dim>>;
     108             : 
     109             :   template <typename TagsList>
     110           0 :   void reconstruct(
     111             :       gsl::not_null<std::array<Variables<TagsList>, Dim>*> vars_on_lower_face,
     112             :       gsl::not_null<std::array<Variables<TagsList>, Dim>*> vars_on_upper_face,
     113             :       const Variables<prims_tags>& volume_prims,
     114             :       const EquationsOfState::EquationOfState<false, 2>& eos,
     115             :       const Element<Dim>& element,
     116             :       const DirectionalIdMap<Dim, evolution::dg::subcell::GhostData>&
     117             :           ghost_data,
     118             :       const Mesh<Dim>& subcell_mesh) const;
     119             : 
     120             :   /// Called by an element doing DG when the neighbor is doing subcell.
     121             :   template <typename TagsList>
     122           1 :   void reconstruct_fd_neighbor(
     123             :       gsl::not_null<Variables<TagsList>*> vars_on_face,
     124             :       const Variables<prims_tags>& subcell_volume_prims,
     125             :       const EquationsOfState::EquationOfState<false, 2>& eos,
     126             :       const Element<Dim>& element,
     127             :       const DirectionalIdMap<Dim, evolution::dg::subcell::GhostData>&
     128             :           ghost_data,
     129             :       const Mesh<Dim>& subcell_mesh,
     130             :       const Direction<Dim> direction_to_reconstruct) const;
     131             : 
     132             :  private:
     133             :   template <size_t LocalDim>
     134             :   // NOLINTNEXTLINE(readability-redundant-declaration)
     135           0 :   friend bool operator==(const AoWeno53Prim<LocalDim>& lhs,
     136             :                          const AoWeno53Prim<LocalDim>& rhs);
     137             : 
     138           0 :   double gamma_hi_ = std::numeric_limits<double>::signaling_NaN();
     139           0 :   double gamma_lo_ = std::numeric_limits<double>::signaling_NaN();
     140           0 :   double epsilon_ = std::numeric_limits<double>::signaling_NaN();
     141           0 :   size_t nonlinear_weight_exponent_ = 0;
     142             : 
     143           0 :   void (*reconstruct_)(gsl::not_null<std::array<gsl::span<double>, Dim>*>,
     144             :                        gsl::not_null<std::array<gsl::span<double>, Dim>*>,
     145             :                        const gsl::span<const double>&,
     146             :                        const DirectionMap<Dim, gsl::span<const double>>&,
     147             :                        const Index<Dim>&, size_t, double, double, double);
     148           0 :   void (*reconstruct_lower_neighbor_)(gsl::not_null<DataVector*>,
     149             :                                       const DataVector&, const DataVector&,
     150             :                                       const Index<Dim>&, const Index<Dim>&,
     151             :                                       const Direction<Dim>&, const double&,
     152             :                                       const double&, const double&);
     153           0 :   void (*reconstruct_upper_neighbor_)(gsl::not_null<DataVector*>,
     154             :                                       const DataVector&, const DataVector&,
     155             :                                       const Index<Dim>&, const Index<Dim>&,
     156             :                                       const Direction<Dim>&, const double&,
     157             :                                       const double&, const double&);
     158             : };
     159             : 
     160             : template <size_t Dim>
     161           0 : bool operator!=(const AoWeno53Prim<Dim>& lhs, const AoWeno53Prim<Dim>& rhs) {
     162             :   return not(lhs == rhs);
     163             : }
     164             : }  // namespace NewtonianEuler::fd