AnalyticSolution.hpp
1 // Distributed under the MIT License.
2 // See LICENSE.txt for details.
3 
4 #pragma once
5 
6 #include <cstddef>
7 #include <ostream>
8 #include <pup.h>
9 #include <string>
10 
11 #include "DataStructures/Tensor/EagerMath/Magnitude.hpp"
12 #include "DataStructures/Tensor/Slice.hpp"
15 #include "Domain/Structure/IndexToSliceAt.hpp"
16 #include "Domain/Tags.hpp"
17 #include "Elliptic/BoundaryConditions/BoundaryCondition.hpp"
18 #include "Elliptic/BoundaryConditions/BoundaryConditionType.hpp"
19 #include "Elliptic/BoundaryConditions/Tags.hpp"
20 #include "NumericalAlgorithms/DiscontinuousGalerkin/NormalDotFlux.hpp"
21 #include "Options/Options.hpp"
23 #include "PointwiseFunctions/AnalyticSolutions/Tags.hpp"
25 #include "Utilities/Gsl.hpp"
26 #include "Utilities/TMPL.hpp"
27 #include "Utilities/TaggedTuple.hpp"
28 #include "Utilities/TypeTraits/IsA.hpp"
29 
31 
32 /// \cond
33 template <typename System, size_t Dim, typename FieldTags, typename FluxTags,
34  typename Registrars>
35 struct AnalyticSolution;
36 
37 namespace Registrars {
38 template <typename System, size_t Dim = System::volume_dim,
39  typename FieldTags = typename System::primal_fields,
40  typename FluxTags = typename System::primal_fluxes>
41 struct AnalyticSolution {
42  template <typename Registrars>
43  using f = BoundaryConditions::AnalyticSolution<System, Dim, FieldTags,
44  FluxTags, Registrars>;
45 };
46 } // namespace Registrars
47 
48 template <typename System, size_t Dim = System::volume_dim,
49  typename FieldTags = typename System::primal_fields,
50  typename FluxTags = typename System::primal_fluxes,
51  typename Registrars =
52  tmpl::list<BoundaryConditions::Registrars::AnalyticSolution<
53  System, Dim, FieldTags, FluxTags>>>
54 struct AnalyticSolution;
55 /// \endcond
56 
57 /*!
58  * \brief Impose the analytic solution on the boundary. Works only if an
59  * analytic solution exists.
60  *
61  * The analytic solution is retrieved from `::Tags::AnalyticSolutionsBase`. It
62  * must hold solutions for both the `System::primal_fields` and the
63  * `System::primal_fluxes`. The user can select to impose the analytic solution
64  * as Dirichlet or Neumann boundary conditions for each field separately.
65  * Dirichlet boundary conditions are imposed on the fields and Neumann boundary
66  * conditions are imposed on the fluxes.
67  *
68  * See `elliptic::Actions::InitializeAnalyticSolutions` for an action that can
69  * add the analytic solutions to the DataBox.
70  */
71 template <typename System, size_t Dim, typename... FieldTags,
72  typename... FluxTags, typename Registrars>
73 class AnalyticSolution<System, Dim, tmpl::list<FieldTags...>,
74  tmpl::list<FluxTags...>, Registrars>
75  : public BoundaryCondition<Dim, Registrars> {
76  private:
78 
79  public:
80  using options =
81  tmpl::list<elliptic::OptionTags::BoundaryConditionType<FieldTags>...>;
82  static constexpr Options::String help =
83  "Boundary conditions from the analytic solution";
84 
85  AnalyticSolution() = default;
86  AnalyticSolution(const AnalyticSolution&) noexcept = default;
87  AnalyticSolution& operator=(const AnalyticSolution&) noexcept = default;
88  AnalyticSolution(AnalyticSolution&&) noexcept = default;
89  AnalyticSolution& operator=(AnalyticSolution&&) noexcept = default;
90  ~AnalyticSolution() noexcept = default;
91 
92  /// \cond
93  explicit AnalyticSolution(CkMigrateMessage* m) noexcept : Base(m) {}
94  using PUP::able::register_constructor;
95  WRAPPED_PUPable_decl_template(AnalyticSolution);
96  /// \endcond
97 
98  /// Select which `elliptic::BoundaryConditionType` to apply for each field
99  explicit AnalyticSolution(
100  // This pack expansion repeats the type `elliptic::BoundaryConditionType`
101  // for each system field
103  FieldTags>::type... boundary_condition_types) noexcept
104  : boundary_condition_types_{boundary_condition_types...} {}
105 
107  const noexcept override {
108  return std::make_unique<AnalyticSolution>(*this);
109  }
110 
111  const auto& boundary_condition_types() const noexcept {
112  return boundary_condition_types_;
113  }
114 
115  using argument_tags =
116  tmpl::list<::Tags::AnalyticSolutionsBase, domain::Tags::Mesh<Dim>,
119  Dim, Frame::Inertial>>>;
120  using volume_tags =
121  tmpl::list<::Tags::AnalyticSolutionsBase, domain::Tags::Mesh<Dim>>;
122 
123  template <typename OptionalAnalyticSolutions>
124  void apply(const gsl::not_null<typename FieldTags::type*>... fields,
125  const gsl::not_null<typename FieldTags::type*>... n_dot_fluxes,
126  const OptionalAnalyticSolutions& optional_analytic_solutions,
127  const Mesh<Dim>& volume_mesh, const Direction<Dim>& direction,
128  const tnsr::i<DataVector, Dim>& face_normal) const noexcept {
129  const auto& analytic_solutions = [&optional_analytic_solutions]() noexcept
130  -> const auto& {
131  if constexpr (tt::is_a_v<std::optional, OptionalAnalyticSolutions>) {
132  if (not optional_analytic_solutions.has_value()) {
134  "Trying to impose boundary conditions from an analytic solution, "
135  "but no analytic solution is available. You probably selected "
136  "the 'AnalyticSolution' boundary condition but chose to solve a "
137  "problem that has no analytic solution. If this is the case, you "
138  "should probably select a different boundary condition.");
139  }
140  return *optional_analytic_solutions;
141  } else {
142  return optional_analytic_solutions;
143  }
144  }
145  ();
146  const size_t slice_index =
147  index_to_slice_at(volume_mesh.extents(), direction);
148  const auto impose_boundary_condition = [this, &analytic_solutions,
149  &volume_mesh, &direction,
150  &slice_index, &face_normal](
151  auto field_tag_v,
152  auto flux_tag_v,
153  const auto field,
154  const auto n_dot_flux) noexcept {
155  using field_tag = decltype(field_tag_v);
156  using flux_tag = decltype(flux_tag_v);
158  boundary_condition_types())) {
161  field, get<::Tags::Analytic<field_tag>>(analytic_solutions),
162  volume_mesh.extents(), direction.dimension(), slice_index);
163  break;
165  normal_dot_flux(
166  n_dot_flux, face_normal,
167  data_on_slice(get<::Tags::Analytic<flux_tag>>(analytic_solutions),
168  volume_mesh.extents(), direction.dimension(),
169  slice_index));
170  break;
171  default:
172  ERROR("Unsupported boundary condition type: "
174  boundary_condition_types()));
175  }
176  };
177  EXPAND_PACK_LEFT_TO_RIGHT(impose_boundary_condition(FieldTags{}, FluxTags{},
178  fields, n_dot_fluxes));
179  }
180 
181  using argument_tags_linearized = tmpl::list<>;
182  using volume_tags_linearized = tmpl::list<>;
183 
184  void apply_linearized(
186  const gsl::not_null<typename FieldTags::type*>... n_dot_fluxes)
187  const noexcept {
188  const auto impose_boundary_condition = [this](
189  auto field_tag_v,
190  const auto field,
191  const auto n_dot_flux) noexcept {
192  using field_tag = decltype(field_tag_v);
194  boundary_condition_types())) {
196  for (auto& field_component : *field) {
197  field_component = 0.;
198  }
199  break;
201  for (auto& n_dot_flux_component : *n_dot_flux) {
202  n_dot_flux_component = 0.;
203  }
204  break;
205  default:
206  ERROR("Unsupported boundary condition type: "
208  boundary_condition_types()));
209  }
210  };
212  impose_boundary_condition(FieldTags{}, fields, n_dot_fluxes));
213  }
214 
215  // NOLINTNEXTLINE
216  void pup(PUP::er& p) noexcept override;
217 
218  private:
219  friend bool operator==(const AnalyticSolution& lhs,
220  const AnalyticSolution& rhs) noexcept {
221  return lhs.boundary_condition_types_ == rhs.boundary_condition_types_;
222  }
223 
224  friend bool operator!=(const AnalyticSolution& lhs,
225  const AnalyticSolution& rhs) noexcept {
226  return not(lhs == rhs);
227  }
228 
230  boundary_condition_types_{};
231 };
232 
233 template <typename System, size_t Dim, typename... FieldTags,
234  typename... FluxTags, typename Registrars>
235 void AnalyticSolution<System, Dim, tmpl::list<FieldTags...>,
236  tmpl::list<FluxTags...>,
237  Registrars>::pup(PUP::er& p) noexcept {
238  Base::pup(p);
239  p | boundary_condition_types_;
240 }
241 
242 /// \cond
243 template <typename System, size_t Dim, typename... FieldTags,
244  typename... FluxTags, typename Registrars>
245 PUP::able::PUP_ID AnalyticSolution<System, Dim, tmpl::list<FieldTags...>,
246  tmpl::list<FluxTags...>,
247  Registrars>::my_PUP_ID = 0; // NOLINT
248 /// \endcond
249 
250 } // namespace elliptic::BoundaryConditions
domain::Tags::UnnormalizedFaceNormal
Definition: FaceNormal.hpp:112
std::apply
T apply(T... args)
CharmPupable.hpp
EXPAND_PACK_LEFT_TO_RIGHT
#define EXPAND_PACK_LEFT_TO_RIGHT(...)
Expand a parameter pack evaluating the terms from left to right.
Definition: TMPL.hpp:565
Frame::Inertial
Definition: IndexType.hpp:44
elliptic::Tags::BoundaryConditionType
The elliptic::BoundaryConditionType to impose on the variable represented by Tag, e....
Definition: Tags.hpp:33
index_to_slice_at
size_t index_to_slice_at(const Index< Dim > &extents, const Direction< Dim > &direction, const size_t offset=0) noexcept
Finds the index in the perpendicular dimension of an element boundary.
Definition: IndexToSliceAt.hpp:18
std::rel_ops::operator!=
T operator!=(T... args)
Options.hpp
Tags.hpp
Error.hpp
elliptic::BoundaryConditionType::Dirichlet
@ Dirichlet
Dirichlet boundary conditions like .
elliptic::BoundaryConditionType::Neumann
@ Neumann
Neumann boundary conditions like , where is the normal to the domain boundary.
db::get
const auto & get(const DataBox< TagList > &box) noexcept
Retrieve the item with tag Tag from the DataBox.
Definition: DataBox.hpp:791
Direction< Dim >
ERROR
#define ERROR(m)
prints an error message to the standard error stream and aborts the program.
Definition: Error.hpp:37
Tags::Analytic
Prefix indicating the analytic solution value for a quantity.
Definition: Tags.hpp:88
WRAPPED_PUPable_decl_template
#define WRAPPED_PUPable_decl_template(className)
Mark derived classes as serializable.
Definition: CharmPupable.hpp:22
Tags::Normalized
Definition: Magnitude.hpp:137
db::data_on_slice
Variables< tmpl::list< TagsToSlice... > > data_on_slice(const db::DataBox< TagsList > &box, const Index< VolumeDim > &element_extents, const size_t sliced_dim, const size_t fixed_index, tmpl::list< TagsToSlice... >) noexcept
Slices volume Tensors from a DataBox into a Variables
Definition: DataOnSlice.hpp:33
cstddef
ERROR_NO_TRACE
#define ERROR_NO_TRACE(m)
Same as ERROR but does not print a backtrace. Intended to be used for user errors,...
Definition: Error.hpp:72
elliptic::BoundaryConditions::BoundaryCondition
Base class for boundary conditions for elliptic systems.
Definition: BoundaryCondition.hpp:91
tuples::TaggedTuple
An associative container that is indexed by structs.
Definition: TaggedTuple.hpp:271
elliptic::BoundaryConditions
Boundary conditions for elliptic systems.
Definition: AnalyticSolution.hpp:30
Mesh
Holds the number of grid points, basis, and quadrature in each direction of the computational grid.
Definition: Mesh.hpp:49
elliptic::BoundaryConditions::AnalyticSolution< System, Dim, tmpl::list< FieldTags... >, tmpl::list< FluxTags... >, Registrars >::AnalyticSolution
AnalyticSolution(const typename elliptic::OptionTags::BoundaryConditionType< FieldTags >::type... boundary_condition_types) noexcept
Select which elliptic::BoundaryConditionType to apply for each field.
Definition: AnalyticSolution.hpp:99
Gsl.hpp
Options::String
const char *const String
The string used in option structs.
Definition: Options.hpp:32
domain::Tags::Direction
Definition: Tags.hpp:383
Tensor.hpp
Direction.hpp
ostream
std::unique_ptr< domain::BoundaryConditions::BoundaryCondition >
elliptic::OptionTags::BoundaryConditionType
Definition: Tags.hpp:19
TMPL.hpp
gsl::not_null
Require a pointer to not be a nullptr
Definition: ReadSpecPiecewisePolynomial.hpp:13
string