DirichletAnalytic.hpp
1 // 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 
13 #include "DataStructures/DataVector.hpp"
14 #include "DataStructures/Tags/TempTensor.hpp"
17 #include "Evolution/BoundaryConditions/Type.hpp"
18 #include "Evolution/Systems/RadiationTransport/M1Grey/BoundaryConditions/BoundaryCondition.hpp"
19 #include "Evolution/Systems/RadiationTransport/M1Grey/Fluxes.hpp"
21 #include "Evolution/Systems/RadiationTransport/M1Grey/Tags.hpp"
22 #include "Options/Options.hpp"
24 #include "PointwiseFunctions/AnalyticData/Tags.hpp"
25 #include "PointwiseFunctions/AnalyticSolutions/AnalyticSolution.hpp"
26 #include "PointwiseFunctions/GeneralRelativity/Tags.hpp"
27 #include "PointwiseFunctions/Hydro/Tags.hpp"
28 #include "Time/Tags.hpp"
29 #include "Utilities/Gsl.hpp"
30 #include "Utilities/TMPL.hpp"
31 
32 /// \cond
33 namespace domain::Tags {
34 template <size_t Dim, typename Frame>
35 struct Coordinates;
36 } // namespace domain::Tags
37 /// \endcond
38 
40 
41 /// \cond
42 template <typename NeutrinoSpeciesList>
43 class DirichletAnalytic;
44 /// \endcond
45 
46 /*!
47  * \brief Sets Dirichlet boundary conditions using the analytic solution or
48  * analytic data.
49  */
50 template <typename... NeutrinoSpecies>
51 class DirichletAnalytic<tmpl::list<NeutrinoSpecies...>> final
52  : public BoundaryCondition<tmpl::list<NeutrinoSpecies...>> {
53  public:
54  using options = tmpl::list<>;
55  static constexpr Options::String help{
56  "DirichletAnalytic boundary conditions using either analytic solution or "
57  "analytic data."};
58 
59  DirichletAnalytic() = default;
60  DirichletAnalytic(DirichletAnalytic&&) noexcept = default;
61  DirichletAnalytic& operator=(DirichletAnalytic&&) noexcept = default;
62  DirichletAnalytic(const DirichletAnalytic&) = default;
63  DirichletAnalytic& operator=(const DirichletAnalytic&) = default;
64  ~DirichletAnalytic() override = default;
65 
66  explicit DirichletAnalytic(CkMigrateMessage* msg) noexcept
67  : BoundaryCondition<tmpl::list<NeutrinoSpecies...>>(msg) {}
68 
71 
72  auto get_clone() const noexcept -> std::unique_ptr<
74  return std::make_unique<DirichletAnalytic>(*this);
75  }
76 
77  static constexpr evolution::BoundaryConditions::Type bc_type =
78  evolution::BoundaryConditions::Type::Ghost;
79 
80  void pup(PUP::er& p) override {
81  BoundaryCondition<tmpl::list<NeutrinoSpecies...>>::pup(p);
82  }
83 
84  using dg_interior_evolved_variables_tags = tmpl::list<>;
85  using dg_interior_temporary_tags =
86  tmpl::list<domain::Tags::Coordinates<3, Frame::Inertial>>;
87  using dg_interior_primitive_variables_tags = tmpl::list<>;
88  using dg_gridless_tags =
89  tmpl::list<::Tags::Time, ::Tags::AnalyticSolutionOrData>;
90 
91  template <typename AnalyticSolutionOrData>
93  const gsl::not_null<typename Tags::TildeE<
94  Frame::Inertial, NeutrinoSpecies>::type*>... tilde_e,
95  const gsl::not_null<typename Tags::TildeS<
96  Frame::Inertial, NeutrinoSpecies>::type*>... tilde_s,
97 
98  const gsl::not_null<typename ::Tags::Flux<
100  Frame::Inertial>::type*>... flux_tilde_e,
101  const gsl::not_null<typename ::Tags::Flux<
103  Frame::Inertial>::type*>... flux_tilde_s,
104 
105  const gsl::not_null<tnsr::II<DataVector, 3, Frame::Inertial>*>
106  inv_spatial_metric,
107 
108  const std::optional<
109  tnsr::I<DataVector, 3, Frame::Inertial>>& /*face_mesh_velocity*/,
110  const tnsr::i<DataVector, 3, Frame::Inertial>& /*normal_covector*/,
111  const tnsr::I<DataVector, 3, Frame::Inertial>& /*normal_vector*/,
112  const tnsr::I<DataVector, 3, Frame::Inertial>& coords, const double time,
113  const AnalyticSolutionOrData& analytic_solution_or_data) const noexcept {
114  auto boundary_values = [&analytic_solution_or_data, &coords,
115  &time]() noexcept {
116  if constexpr (std::is_base_of_v<MarkAsAnalyticSolution,
117  AnalyticSolutionOrData>) {
118  return analytic_solution_or_data.variables(
119  coords, time,
127 
128  } else {
129  (void)time;
130  return analytic_solution_or_data.variables(
131  coords,
139  }
140  }();
141 
142  *inv_spatial_metric =
143  get<gr::Tags::InverseSpatialMetric<3>>(boundary_values);
144 
145  // Allocate the temporary tensors outside the loop over species
146  Variables<tmpl::list<::Tags::TempScalar<0>, ::Tags::TempII<0, 3>,
149  buffer((*inv_spatial_metric)[0].size());
150 
151  const auto assign_boundary_values_for_neutrino_species =
152  [&boundary_values, &inv_spatial_metric, &buffer](
153  const gsl::not_null<Scalar<DataVector>*> local_tilde_e,
155  local_tilde_s,
157  local_flux_tilde_e,
159  local_flux_tilde_s,
160  auto species_v) noexcept {
161  using species = decltype(species_v);
162  using tilde_e_tag = Tags::TildeE<Frame::Inertial, species>;
163  using tilde_s_tag = Tags::TildeS<Frame::Inertial, species>;
164  *local_tilde_e = get<tilde_e_tag>(boundary_values);
165  *local_tilde_s = get<tilde_s_tag>(boundary_values);
166 
167  // Compute pressure tensor tilde_p from the M1Closure
168  const auto& fluid_velocity =
169  get<hydro::Tags::SpatialVelocity<DataVector, 3>>(boundary_values);
170  const auto& fluid_lorentz_factor =
171  get<hydro::Tags::LorentzFactor<DataVector>>(boundary_values);
172  const auto& lapse = get<gr::Tags::Lapse<>>(boundary_values);
173  const auto& shift = get<gr::Tags::Shift<3>>(boundary_values);
174  const auto& spatial_metric =
175  get<gr::Tags::SpatialMetric<3>>(boundary_values);
176  auto& closure_factor = get<::Tags::TempScalar<0>>(buffer);
177  // The M1Closure reads in values from `closure_factor` as an initial
178  // guess. We need to specify some value (else code fails in DEBUG
179  // mode when the temp Variables are initialized with NaNs)... for now
180  // we use 0 because it's easy, but improvements may be possible.
181  get(closure_factor) = 0.;
182  auto& pressure_tensor = get<::Tags::TempII<0, 3>>(buffer);
183  auto& comoving_energy_density = get<::Tags::TempScalar<1>>(buffer);
184  auto& comoving_momentum_density_normal =
185  get<::Tags::TempScalar<2>>(buffer);
186  auto& comoving_momentum_density_spatial =
187  get<::Tags::Tempi<0, 3>>(buffer);
188  detail::compute_closure_impl(
189  make_not_null(&closure_factor), make_not_null(&pressure_tensor),
190  make_not_null(&comoving_energy_density),
191  make_not_null(&comoving_momentum_density_normal),
192  make_not_null(&comoving_momentum_density_spatial), *local_tilde_e,
193  *local_tilde_s, fluid_velocity, fluid_lorentz_factor,
194  spatial_metric, *inv_spatial_metric);
195  // Store det of metric in (otherwise unused) comoving_energy_density
196  get(comoving_energy_density) = get(determinant(spatial_metric));
197  for (auto& component : pressure_tensor) {
198  component *= get(comoving_energy_density);
199  }
200  const auto& tilde_p = pressure_tensor;
201 
202  auto& tilde_s_M = get<::Tags::TempI<0, 3>>(buffer);
203  detail::compute_fluxes_impl(local_flux_tilde_e, local_flux_tilde_s,
204  &tilde_s_M, *local_tilde_e,
205  *local_tilde_s, tilde_p, lapse, shift,
206  spatial_metric, *inv_spatial_metric);
207  return '0';
208  };
209 
210  expand_pack(assign_boundary_values_for_neutrino_species(
211  tilde_e, tilde_s, flux_tilde_e, flux_tilde_s, NeutrinoSpecies{})...);
212  return {};
213  }
214 };
215 
216 /// \cond
217 template <typename... NeutrinoSpecies>
218 // NOLINTNEXTLINE
219 PUP::able::PUP_ID DirichletAnalytic<tmpl::list<NeutrinoSpecies...>>::my_PUP_ID =
220  0;
221 /// \endcond
222 
223 } // namespace RadiationTransport::M1Grey::BoundaryConditions
expand_pack
constexpr void expand_pack(Ts &&...) noexcept
Allows zero-cost unordered expansion of a parameter.
Definition: TMPL.hpp:585
CharmPupable.hpp
domain::BoundaryConditions::BoundaryCondition
Base class from which all system-specific base classes must inherit.
Definition: BoundaryCondition.hpp:18
Frame::Inertial
Definition: IndexType.hpp:44
gr::Tags::SpatialMetric
Definition: Tags.hpp:26
M1Closure.hpp
Options.hpp
db::get
const auto & get(const DataBox< TagList > &box) noexcept
Retrieve the item with tag Tag from the DataBox.
Definition: DataBox.hpp:791
RadiationTransport::M1Grey::BoundaryConditions::BoundaryCondition
The base class off of which all boundary conditions must inherit.
Definition: BoundaryCondition.hpp:22
RadiationTransport::M1Grey::Tags::TildeE
The densitized energy density of neutrinos of a given species .
Definition: Tags.hpp:31
determinant
void determinant(const gsl::not_null< Scalar< T > * > det_tensor, const Tensor< T, Symm, index_list< Index0, Index1 >> &tensor) noexcept
Computes the determinant of a rank-2 Tensor tensor.
Definition: Determinant.hpp:139
WRAPPED_PUPable_decl_base_template
#define WRAPPED_PUPable_decl_base_template(baseClassName, className)
Mark derived template classes as serializable.
Definition: CharmPupable.hpp:32
RadiationTransport::M1Grey::Tags::TildeS
The densitized momentum density of neutrinos of a given species .
Definition: Tags.hpp:41
gr::lapse
Scalar< DataType > lapse(const tnsr::I< DataType, SpatialDim, Frame > &shift, const tnsr::aa< DataType, SpatialDim, Frame > &spacetime_metric) noexcept
Compute lapse from shift and spacetime metric.
hydro::Tags::LorentzFactor
The Lorentz factor , where is the spatial velocity of the fluid.
Definition: Tags.hpp:64
memory
gr::Tags::Shift
Definition: Tags.hpp:48
RadiationTransport::M1Grey::BoundaryConditions
Boundary conditions for the M1Grey radiation transport system.
Definition: BoundaryCondition.hpp:13
Variables.hpp
gr::shift
tnsr::I< DataType, SpatialDim, Frame > shift(const tnsr::aa< DataType, SpatialDim, Frame > &spacetime_metric, const tnsr::II< DataType, SpatialDim, Frame > &inverse_spatial_metric) noexcept
Compute shift from spacetime metric and inverse spatial metric.
Tags::TempTensor
Definition: TempTensor.hpp:21
Gsl.hpp
Options::String
const char *const String
The string used in option structs.
Definition: Options.hpp:32
hydro::Tags::SpatialVelocity
The spatial velocity of the fluid, where . Here is the spatial part of the 4-velocity of the fluid,...
Definition: Tags.hpp:142
Tensor.hpp
gr::spatial_metric
tnsr::ii< DataType, SpatialDim, Frame > spatial_metric(const tnsr::aa< DataType, SpatialDim, Frame > &spacetime_metric) noexcept
Compute spatial metric from spacetime metric.
optional
make_not_null
gsl::not_null< T * > make_not_null(T *ptr) noexcept
Construct a not_null from a pointer. Often this will be done as an implicit conversion,...
Definition: Gsl.hpp:880
std::unique_ptr
Prefixes.hpp
type_traits
gr::Tags::Lapse
Definition: Tags.hpp:52
TMPL.hpp
domain::Tags
Tags for the domain.
Definition: FaceNormal.hpp:107
gr::Tags::InverseSpatialMetric
Inverse of the spatial metric.
Definition: Tags.hpp:33
gsl::not_null
Require a pointer to not be a nullptr
Definition: ReadSpecPiecewisePolynomial.hpp:13
string