Line data    Source code

       1           0 : // Distributed under the MIT License.
       2             : // See LICENSE.txt for details.
       3             : 
       4             : #pragma once
       5             : 
       6             : #include <memory>
       7             : #include <optional>
       8             : #include <pup.h>
       9             : #include <string>
      10             : #include <type_traits>
      11             : 
      12             : #include "DataStructures/DataBox/Prefixes.hpp"
      13             : #include "DataStructures/DataVector.hpp"
      14             : #include "DataStructures/Tags/TempTensor.hpp"
      15             : #include "DataStructures/Tensor/Tensor.hpp"
      16             : #include "DataStructures/Variables.hpp"
      17             : #include "Evolution/BoundaryConditions/Type.hpp"
      18             : #include "Evolution/Systems/RadiationTransport/M1Grey/BoundaryConditions/BoundaryCondition.hpp"
      19             : #include "Evolution/Systems/RadiationTransport/M1Grey/Fluxes.hpp"
      20             : #include "Evolution/Systems/RadiationTransport/M1Grey/M1Closure.hpp"
      21             : #include "Evolution/Systems/RadiationTransport/M1Grey/Tags.hpp"
      22             : #include "Options/String.hpp"
      23             : #include "PointwiseFunctions/AnalyticData/Tags.hpp"
      24             : #include "PointwiseFunctions/AnalyticSolutions/AnalyticSolution.hpp"
      25             : #include "PointwiseFunctions/GeneralRelativity/Tags.hpp"
      26             : #include "PointwiseFunctions/Hydro/Tags.hpp"
      27             : #include "Utilities/Gsl.hpp"
      28             : #include "Utilities/Serialization/CharmPupable.hpp"
      29             : #include "Utilities/TMPL.hpp"
      30             : 
      31             : /// \cond
      32             : namespace Tags {
      33             : struct Time;
      34             : }  // namespace Tags
      35             : namespace domain::Tags {
      36             : template <size_t Dim, typename Frame>
      37             : struct Coordinates;
      38             : }  // namespace domain::Tags
      39             : /// \endcond
      40             : 
      41             : namespace RadiationTransport::M1Grey::BoundaryConditions {
      42             : 
      43             : /// \cond
      44             : template <typename NeutrinoSpeciesList>
      45             : class DirichletAnalytic;
      46             : /// \endcond
      47             : 
      48             : /*!
      49             :  * \brief Sets Dirichlet boundary conditions using the analytic solution or
      50             :  * analytic data.
      51             :  */
      52             : template <typename... NeutrinoSpecies>
      53           1 : class DirichletAnalytic<tmpl::list<NeutrinoSpecies...>> final
      54             :     : public BoundaryCondition<tmpl::list<NeutrinoSpecies...>> {
      55             :  public:
      56           0 :   using options = tmpl::list<>;
      57           0 :   static constexpr Options::String help{
      58             :       "DirichletAnalytic boundary conditions using either analytic solution or "
      59             :       "analytic data."};
      60             : 
      61           0 :   DirichletAnalytic() = default;
      62           0 :   DirichletAnalytic(DirichletAnalytic&&) = default;
      63           0 :   DirichletAnalytic& operator=(DirichletAnalytic&&) = default;
      64           0 :   DirichletAnalytic(const DirichletAnalytic&) = default;
      65           0 :   DirichletAnalytic& operator=(const DirichletAnalytic&) = default;
      66           0 :   ~DirichletAnalytic() override = default;
      67             : 
      68           0 :   explicit DirichletAnalytic(CkMigrateMessage* msg)
      69             :       : BoundaryCondition<tmpl::list<NeutrinoSpecies...>>(msg) {}
      70             : 
      71           0 :   WRAPPED_PUPable_decl_base_template(
      72             :       domain::BoundaryConditions::BoundaryCondition, DirichletAnalytic);
      73             : 
      74           0 :   auto get_clone() const -> std::unique_ptr<
      75             :       domain::BoundaryConditions::BoundaryCondition> override {
      76             :     return std::make_unique<DirichletAnalytic>(*this);
      77             :   }
      78             : 
      79           0 :   static constexpr evolution::BoundaryConditions::Type bc_type =
      80             :       evolution::BoundaryConditions::Type::Ghost;
      81             : 
      82           0 :   void pup(PUP::er& p) override {
      83             :     BoundaryCondition<tmpl::list<NeutrinoSpecies...>>::pup(p);
      84             :   }
      85             : 
      86           0 :   using dg_interior_evolved_variables_tags = tmpl::list<>;
      87           0 :   using dg_interior_temporary_tags =
      88             :       tmpl::list<domain::Tags::Coordinates<3, Frame::Inertial>>;
      89           0 :   using dg_interior_primitive_variables_tags = tmpl::list<>;
      90           0 :   using dg_gridless_tags =
      91             :       tmpl::list<::Tags::Time, ::Tags::AnalyticSolutionOrData>;
      92             : 
      93             :   template <typename AnalyticSolutionOrData>
      94           0 :   std::optional<std::string> dg_ghost(
      95             :       const gsl::not_null<typename Tags::TildeE<
      96             :           Frame::Inertial, NeutrinoSpecies>::type*>... tilde_e,
      97             :       const gsl::not_null<typename Tags::TildeS<
      98             :           Frame::Inertial, NeutrinoSpecies>::type*>... tilde_s,
      99             : 
     100             :       const gsl::not_null<typename ::Tags::Flux<
     101             :           Tags::TildeE<Frame::Inertial, NeutrinoSpecies>, tmpl::size_t<3>,
     102             :           Frame::Inertial>::type*>... flux_tilde_e,
     103             :       const gsl::not_null<typename ::Tags::Flux<
     104             :           Tags::TildeS<Frame::Inertial, NeutrinoSpecies>, tmpl::size_t<3>,
     105             :           Frame::Inertial>::type*>... flux_tilde_s,
     106             : 
     107             :       const gsl::not_null<tnsr::II<DataVector, 3, Frame::Inertial>*>
     108             :           inv_spatial_metric,
     109             : 
     110             :       const std::optional<
     111             :           tnsr::I<DataVector, 3, Frame::Inertial>>& /*face_mesh_velocity*/,
     112             :       const tnsr::i<DataVector, 3, Frame::Inertial>& /*normal_covector*/,
     113             :       const tnsr::I<DataVector, 3, Frame::Inertial>& /*normal_vector*/,
     114             :       const tnsr::I<DataVector, 3, Frame::Inertial>& coords, const double time,
     115             :       const AnalyticSolutionOrData& analytic_solution_or_data) const {
     116             :     auto boundary_values = [&analytic_solution_or_data, &coords, &time]() {
     117             :       if constexpr (is_analytic_solution_v<AnalyticSolutionOrData>) {
     118             :         return analytic_solution_or_data.variables(
     119             :             coords, time,
     120             :             tmpl::list<Tags::TildeE<Frame::Inertial, NeutrinoSpecies>...,
     121             :                        Tags::TildeS<Frame::Inertial, NeutrinoSpecies>...,
     122             :                        hydro::Tags::LorentzFactor<DataVector>,
     123             :                        hydro::Tags::SpatialVelocity<DataVector, 3>,
     124             :                        gr::Tags::Lapse<DataVector>,
     125             :                        gr::Tags::Shift<DataVector, 3>,
     126             :                        gr::Tags::SpatialMetric<DataVector, 3>,
     127             :                        gr::Tags::InverseSpatialMetric<DataVector, 3>>{});
     128             : 
     129             :       } else {
     130             :         (void)time;
     131             :         return analytic_solution_or_data.variables(
     132             :             coords,
     133             :             tmpl::list<Tags::TildeE<Frame::Inertial, NeutrinoSpecies>...,
     134             :                        Tags::TildeS<Frame::Inertial, NeutrinoSpecies>...,
     135             :                        hydro::Tags::LorentzFactor<DataVector>,
     136             :                        hydro::Tags::SpatialVelocity<DataVector, 3>,
     137             :                        gr::Tags::Lapse<DataVector>,
     138             :                        gr::Tags::Shift<DataVector, 3>,
     139             :                        gr::Tags::SpatialMetric<DataVector, 3>,
     140             :                        gr::Tags::InverseSpatialMetric<DataVector, 3>>{});
     141             :       }
     142             :     }();
     143             : 
     144             :     *inv_spatial_metric =
     145             :         get<gr::Tags::InverseSpatialMetric<DataVector, 3>>(boundary_values);
     146             : 
     147             :     // Allocate the temporary tensors outside the loop over species
     148             :     Variables<tmpl::list<::Tags::TempScalar<0>, ::Tags::TempII<0, 3>,
     149             :                          ::Tags::TempScalar<1>, ::Tags::TempScalar<2>,
     150             :                          ::Tags::Tempi<0, 3>, ::Tags::TempI<0, 3>>>
     151             :         buffer((*inv_spatial_metric)[0].size());
     152             : 
     153             :     const auto assign_boundary_values_for_neutrino_species =
     154             :         [&boundary_values, &inv_spatial_metric, &buffer](
     155             :             const gsl::not_null<Scalar<DataVector>*> local_tilde_e,
     156             :             const gsl::not_null<tnsr::i<DataVector, 3, Frame::Inertial>*>
     157             :                 local_tilde_s,
     158             :             const gsl::not_null<tnsr::I<DataVector, 3, Frame::Inertial>*>
     159             :                 local_flux_tilde_e,
     160             :             const gsl::not_null<tnsr::Ij<DataVector, 3, Frame::Inertial>*>
     161             :                 local_flux_tilde_s,
     162             :             auto species_v) {
     163             :           using species = decltype(species_v);
     164             :           using tilde_e_tag = Tags::TildeE<Frame::Inertial, species>;
     165             :           using tilde_s_tag = Tags::TildeS<Frame::Inertial, species>;
     166             :           *local_tilde_e = get<tilde_e_tag>(boundary_values);
     167             :           *local_tilde_s = get<tilde_s_tag>(boundary_values);
     168             : 
     169             :           // Compute pressure tensor tilde_p from the M1Closure
     170             :           const auto& fluid_velocity =
     171             :               get<hydro::Tags::SpatialVelocity<DataVector, 3>>(boundary_values);
     172             :           const auto& fluid_lorentz_factor =
     173             :               get<hydro::Tags::LorentzFactor<DataVector>>(boundary_values);
     174             :           const auto& lapse = get<gr::Tags::Lapse<DataVector>>(boundary_values);
     175             :           const auto& shift =
     176             :               get<gr::Tags::Shift<DataVector, 3>>(boundary_values);
     177             :           const auto& spatial_metric =
     178             :               get<gr::Tags::SpatialMetric<DataVector, 3>>(boundary_values);
     179             :           auto& closure_factor = get<::Tags::TempScalar<0>>(buffer);
     180             :           // The M1Closure reads in values from `closure_factor` as an initial
     181             :           // guess. We need to specify some value (else code fails in DEBUG
     182             :           // mode when the temp Variables are initialized with NaNs)... for now
     183             :           // we use 0 because it's easy, but improvements may be possible.
     184             :           get(closure_factor) = 0.;
     185             :           auto& pressure_tensor = get<::Tags::TempII<0, 3>>(buffer);
     186             :           auto& comoving_energy_density = get<::Tags::TempScalar<1>>(buffer);
     187             :           auto& comoving_momentum_density_normal =
     188             :               get<::Tags::TempScalar<2>>(buffer);
     189             :           auto& comoving_momentum_density_spatial =
     190             :               get<::Tags::Tempi<0, 3>>(buffer);
     191             :           detail::compute_closure_impl(
     192             :               make_not_null(&closure_factor), make_not_null(&pressure_tensor),
     193             :               make_not_null(&comoving_energy_density),
     194             :               make_not_null(&comoving_momentum_density_normal),
     195             :               make_not_null(&comoving_momentum_density_spatial), *local_tilde_e,
     196             :               *local_tilde_s, fluid_velocity, fluid_lorentz_factor,
     197             :               spatial_metric, *inv_spatial_metric);
     198             :           // Store det of metric in (otherwise unused) comoving_energy_density
     199             :           get(comoving_energy_density) = get(determinant(spatial_metric));
     200             :           for (auto& component : pressure_tensor) {
     201             :             component *= get(comoving_energy_density);
     202             :           }
     203             :           const auto& tilde_p = pressure_tensor;
     204             : 
     205             :           auto& tilde_s_M = get<::Tags::TempI<0, 3>>(buffer);
     206             :           detail::compute_fluxes_impl(local_flux_tilde_e, local_flux_tilde_s,
     207             :                                       &tilde_s_M, *local_tilde_e,
     208             :                                       *local_tilde_s, tilde_p, lapse, shift,
     209             :                                       spatial_metric, *inv_spatial_metric);
     210             :           return '0';
     211             :         };
     212             : 
     213             :     expand_pack(assign_boundary_values_for_neutrino_species(
     214             :         tilde_e, tilde_s, flux_tilde_e, flux_tilde_s, NeutrinoSpecies{})...);
     215             :     return {};
     216             :   }
     217             : };
     218             : 
     219             : /// \cond
     220             : template <typename... NeutrinoSpecies>
     221             : // NOLINTNEXTLINE
     222             : PUP::able::PUP_ID DirichletAnalytic<tmpl::list<NeutrinoSpecies...>>::my_PUP_ID =
     223             :     0;
     224             : /// \endcond
     225             : 
     226             : }  // namespace RadiationTransport::M1Grey::BoundaryConditions