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  static std::string name() noexcept { return "DirichletAnalytic"; }
59 
60  DirichletAnalytic() = default;
61  DirichletAnalytic(DirichletAnalytic&&) noexcept = default;
62  DirichletAnalytic& operator=(DirichletAnalytic&&) noexcept = default;
63  DirichletAnalytic(const DirichletAnalytic&) = default;
64  DirichletAnalytic& operator=(const DirichletAnalytic&) = default;
65  ~DirichletAnalytic() override = default;
66 
67  explicit DirichletAnalytic(CkMigrateMessage* msg) noexcept
68  : BoundaryCondition<tmpl::list<NeutrinoSpecies...>>(msg) {}
69 
72 
73  auto get_clone() const noexcept -> std::unique_ptr<
75  return std::make_unique<DirichletAnalytic>(*this);
76  }
77 
78  static constexpr evolution::BoundaryConditions::Type bc_type =
79  evolution::BoundaryConditions::Type::Ghost;
80 
81  void pup(PUP::er& p) override {
82  BoundaryCondition<tmpl::list<NeutrinoSpecies...>>::pup(p);
83  }
84 
85  using dg_interior_evolved_variables_tags = tmpl::list<>;
86  using dg_interior_temporary_tags =
87  tmpl::list<domain::Tags::Coordinates<3, Frame::Inertial>>;
88  using dg_interior_primitive_variables_tags = tmpl::list<>;
89  using dg_gridless_tags =
90  tmpl::list<::Tags::Time, ::Tags::AnalyticSolutionOrData>;
91 
92  template <typename AnalyticSolutionOrData>
94  const gsl::not_null<typename Tags::TildeE<
95  Frame::Inertial, NeutrinoSpecies>::type*>... tilde_e,
96  const gsl::not_null<typename Tags::TildeS<
97  Frame::Inertial, NeutrinoSpecies>::type*>... tilde_s,
98 
99  const gsl::not_null<typename ::Tags::Flux<
101  Frame::Inertial>::type*>... flux_tilde_e,
102  const gsl::not_null<typename ::Tags::Flux<
104  Frame::Inertial>::type*>... flux_tilde_s,
105 
106  const gsl::not_null<tnsr::II<DataVector, 3, Frame::Inertial>*>
107  inv_spatial_metric,
108 
109  const std::optional<
110  tnsr::I<DataVector, 3, Frame::Inertial>>& /*face_mesh_velocity*/,
111  const tnsr::i<DataVector, 3, Frame::Inertial>& /*normal_covector*/,
112  const tnsr::I<DataVector, 3, Frame::Inertial>& /*normal_vector*/,
113  const tnsr::I<DataVector, 3, Frame::Inertial>& coords, const double time,
114  const AnalyticSolutionOrData& analytic_solution_or_data) const noexcept {
115  auto boundary_values = [&analytic_solution_or_data, &coords,
116  &time]() noexcept {
117  if constexpr (std::is_base_of_v<MarkAsAnalyticSolution,
118  AnalyticSolutionOrData>) {
119  return analytic_solution_or_data.variables(
120  coords, time,
128 
129  } else {
130  (void)time;
131  return analytic_solution_or_data.variables(
132  coords,
140  }
141  }();
142 
143  *inv_spatial_metric =
144  get<gr::Tags::InverseSpatialMetric<3>>(boundary_values);
145 
146  // Allocate the temporary tensors outside the loop over species
147  Variables<tmpl::list<::Tags::TempScalar<0>, ::Tags::TempII<0, 3>,
150  buffer((*inv_spatial_metric)[0].size());
151 
152  const auto assign_boundary_values_for_neutrino_species =
153  [&boundary_values, &inv_spatial_metric, &buffer](
154  const gsl::not_null<Scalar<DataVector>*> local_tilde_e,
156  local_tilde_s,
158  local_flux_tilde_e,
160  local_flux_tilde_s,
161  auto species_v) noexcept {
162  using species = decltype(species_v);
163  using tilde_e_tag = Tags::TildeE<Frame::Inertial, species>;
164  using tilde_s_tag = Tags::TildeS<Frame::Inertial, species>;
165  *local_tilde_e = get<tilde_e_tag>(boundary_values);
166  *local_tilde_s = get<tilde_s_tag>(boundary_values);
167 
168  // Compute pressure tensor tilde_p from the M1Closure
169  const auto& fluid_velocity =
170  get<hydro::Tags::SpatialVelocity<DataVector, 3>>(boundary_values);
171  const auto& fluid_lorentz_factor =
172  get<hydro::Tags::LorentzFactor<DataVector>>(boundary_values);
173  const auto& lapse = get<gr::Tags::Lapse<>>(boundary_values);
174  const auto& shift = get<gr::Tags::Shift<3>>(boundary_values);
175  const auto& spatial_metric =
176  get<gr::Tags::SpatialMetric<3>>(boundary_values);
177  auto& closure_factor = get<::Tags::TempScalar<0>>(buffer);
178  // The M1Closure reads in values from `closure_factor` as an initial
179  // guess. We need to specify some value (else code fails in DEBUG
180  // mode when the temp Variables are initialized with NaNs)... for now
181  // we use 0 because it's easy, but improvements may be possible.
182  get(closure_factor) = 0.;
183  auto& pressure_tensor = get<::Tags::TempII<0, 3>>(buffer);
184  auto& comoving_energy_density = get<::Tags::TempScalar<1>>(buffer);
185  auto& comoving_momentum_density_normal =
186  get<::Tags::TempScalar<2>>(buffer);
187  auto& comoving_momentum_density_spatial =
188  get<::Tags::Tempi<0, 3>>(buffer);
189  detail::compute_closure_impl(
190  make_not_null(&closure_factor), make_not_null(&pressure_tensor),
191  make_not_null(&comoving_energy_density),
192  make_not_null(&comoving_momentum_density_normal),
193  make_not_null(&comoving_momentum_density_spatial), *local_tilde_e,
194  *local_tilde_s, fluid_velocity, fluid_lorentz_factor,
195  spatial_metric, *inv_spatial_metric);
196  // Store det of metric in (otherwise unused) comoving_energy_density
197  get(comoving_energy_density) = get(determinant(spatial_metric));
198  for (auto& component : pressure_tensor) {
199  component *= get(comoving_energy_density);
200  }
201  const auto& tilde_p = pressure_tensor;
202 
203  auto& tilde_s_M = get<::Tags::TempI<0, 3>>(buffer);
204  detail::compute_fluxes_impl(local_flux_tilde_e, local_flux_tilde_s,
205  &tilde_s_M, *local_tilde_e,
206  *local_tilde_s, tilde_p, lapse, shift,
207  spatial_metric, *inv_spatial_metric);
208  return '0';
209  };
210 
211  expand_pack(assign_boundary_values_for_neutrino_species(
212  tilde_e, tilde_s, flux_tilde_e, flux_tilde_s, NeutrinoSpecies{})...);
213  return {};
214  }
215 };
216 
217 /// \cond
218 template <typename... NeutrinoSpecies>
219 // NOLINTNEXTLINE
220 PUP::able::PUP_ID DirichletAnalytic<tmpl::list<NeutrinoSpecies...>>::my_PUP_ID =
221  0;
222 /// \endcond
223 
224 } // namespace RadiationTransport::M1Grey::BoundaryConditions
expand_pack
constexpr void expand_pack(Ts &&...) noexcept
Allows zero-cost unordered expansion of a parameter.
Definition: TMPL.hpp:549
std::string
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
get
constexpr Tag::type & get(Variables< TagList > &v) noexcept
Return Tag::type pointing into the contiguous array.
Definition: Variables.hpp:660
M1Closure.hpp
Options.hpp
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: ReadSpecThirdOrderPiecewisePolynomial.hpp:13
string