Line data Source code
1 0 : // Distributed under the MIT License.
2 : // See LICENSE.txt for details.
3 :
4 : #pragma once
5 :
6 : #include <cstddef>
7 : #include <memory>
8 : #include <ostream>
9 : #include <pup.h>
10 : #include <string>
11 : #include <type_traits>
12 : #include <utility>
13 : #include <vector>
14 :
15 : #include "DataStructures/DataBox/MetavariablesTag.hpp"
16 : #include "DataStructures/TaggedTuple.hpp"
17 : #include "DataStructures/Tensor/Tensor.hpp"
18 : #include "Domain/Tags.hpp"
19 : #include "Domain/Tags/FaceNormal.hpp"
20 : #include "Elliptic/BoundaryConditions/BoundaryCondition.hpp"
21 : #include "Elliptic/BoundaryConditions/BoundaryConditionType.hpp"
22 : #include "Elliptic/BoundaryConditions/Tags.hpp"
23 : #include "NumericalAlgorithms/DiscontinuousGalerkin/NormalDotFlux.hpp"
24 : #include "Options/String.hpp"
25 : #include "PointwiseFunctions/InitialDataUtilities/InitialGuess.hpp"
26 : #include "Utilities/CallWithDynamicType.hpp"
27 : #include "Utilities/ErrorHandling/Error.hpp"
28 : #include "Utilities/Gsl.hpp"
29 : #include "Utilities/Serialization/CharmPupable.hpp"
30 : #include "Utilities/Serialization/Serialize.hpp"
31 : #include "Utilities/TMPL.hpp"
32 :
33 1 : namespace elliptic::BoundaryConditions {
34 : namespace detail {
35 :
36 : template <typename Solution, size_t Dim, typename Tag, typename = std::void_t<>>
37 : struct has_boundary_variables : std::false_type {};
38 :
39 : template <typename Solution, size_t Dim, typename Tag>
40 : struct has_boundary_variables<
41 : Solution, Dim, Tag,
42 : std::void_t<decltype(std::declval<const Solution&>().variables(
43 : std::declval<const tnsr::I<DataVector, Dim, Frame::Inertial>&>(),
44 : tmpl::list<Tag>{}))>> : std::true_type {};
45 :
46 : template <typename Solution, size_t Dim, typename Tag>
47 : constexpr bool has_boundary_variables_v =
48 : has_boundary_variables<Solution, Dim, Tag>::value;
49 :
50 : } // namespace detail
51 :
52 : /// \cond
53 : template <typename System, size_t Dim = System::volume_dim,
54 : typename FieldTags = typename System::primal_fields,
55 : typename FluxTags = typename System::primal_fluxes>
56 : struct AnalyticSolution;
57 : /// \endcond
58 :
59 : /*!
60 : * \brief Impose the analytic solution on the boundary.
61 : *
62 : * The user can select to impose the analytic solution as Dirichlet or
63 : * Neumann boundary conditions for each field separately. Dirichlet
64 : * boundary conditions are imposed on the fields and Neumann boundary
65 : * conditions are imposed on the fluxes.
66 : */
67 : template <typename System, size_t Dim, typename... FieldTags,
68 : typename... FluxTags>
69 1 : class AnalyticSolution<System, Dim, tmpl::list<FieldTags...>,
70 : tmpl::list<FluxTags...>>
71 : : public BoundaryCondition<Dim> {
72 : private:
73 0 : using Base = BoundaryCondition<Dim>;
74 :
75 : public:
76 0 : struct Solution {
77 0 : using type = std::unique_ptr<elliptic::analytic_data::InitialGuess>;
78 0 : static constexpr Options::String help = {
79 : "The analytic data to impose on the boundary"};
80 : };
81 :
82 0 : using options =
83 : tmpl::list<Solution,
84 : elliptic::OptionTags::BoundaryConditionType<FieldTags>...>;
85 0 : static constexpr Options::String help =
86 : "Boundary conditions from the analytic solution";
87 :
88 0 : AnalyticSolution() = default;
89 0 : AnalyticSolution(const AnalyticSolution& rhs) : Base(rhs) { *this = rhs; }
90 0 : AnalyticSolution& operator=(const AnalyticSolution& rhs) {
91 : if (rhs.solution_ != nullptr) {
92 : solution_ = serialize_and_deserialize<
93 : std::unique_ptr<elliptic::analytic_data::InitialGuess>>(
94 : rhs.solution_);
95 : } else {
96 : solution_ = nullptr;
97 : }
98 : boundary_condition_types_ = rhs.boundary_condition_types_;
99 : return *this;
100 : }
101 0 : AnalyticSolution(AnalyticSolution&&) = default;
102 0 : AnalyticSolution& operator=(AnalyticSolution&&) = default;
103 0 : ~AnalyticSolution() = default;
104 :
105 : /// \cond
106 : explicit AnalyticSolution(CkMigrateMessage* m) : Base(m) {}
107 : using PUP::able::register_constructor;
108 : WRAPPED_PUPable_decl_template(AnalyticSolution);
109 : /// \endcond
110 :
111 : /// Select which `elliptic::BoundaryConditionType` to apply for each field
112 1 : explicit AnalyticSolution(
113 : std::unique_ptr<elliptic::analytic_data::InitialGuess> solution,
114 : // This pack expansion repeats the type `elliptic::BoundaryConditionType`
115 : // for each system field
116 : const typename elliptic::OptionTags::BoundaryConditionType<
117 : FieldTags>::type... boundary_condition_types)
118 : : solution_(std::move(solution)),
119 : boundary_condition_types_{boundary_condition_types...} {}
120 :
121 0 : std::unique_ptr<domain::BoundaryConditions::BoundaryCondition> get_clone()
122 : const override {
123 : return std::make_unique<AnalyticSolution>(*this);
124 : }
125 :
126 0 : std::vector<elliptic::BoundaryConditionType> boundary_condition_types()
127 : const override {
128 : std::vector<elliptic::BoundaryConditionType> result{};
129 : const auto collect = [&result](
130 : const auto tag_v,
131 : const elliptic::BoundaryConditionType bc_type) {
132 : using tag = std::decay_t<decltype(tag_v)>;
133 : for (size_t i = 0; i < tag::type::size(); ++i) {
134 : result.push_back(bc_type);
135 : }
136 : };
137 : EXPAND_PACK_LEFT_TO_RIGHT(collect(
138 : FieldTags{}, get<elliptic::Tags::BoundaryConditionType<FieldTags>>(
139 : boundary_condition_types_)));
140 : return result;
141 : }
142 :
143 0 : using argument_tags =
144 : tmpl::list<Parallel::Tags::Metavariables,
145 : domain::Tags::Coordinates<Dim, Frame::Inertial>,
146 : ::Tags::Normalized<domain::Tags::UnnormalizedFaceNormal<
147 : Dim, Frame::Inertial>>>;
148 0 : using volume_tags = tmpl::list<Parallel::Tags::Metavariables>;
149 :
150 : template <typename Metavariables>
151 0 : void apply(const gsl::not_null<typename FieldTags::type*>... fields,
152 : const gsl::not_null<typename FieldTags::type*>... n_dot_fluxes,
153 : const TensorMetafunctions::prepend_spatial_index<
154 : typename FieldTags::type, Dim, UpLo::Lo,
155 : Frame::Inertial>&... /*deriv_fields*/,
156 : const Metavariables& /*meta*/,
157 : const tnsr::I<DataVector, Dim>& face_inertial_coords,
158 : const tnsr::i<DataVector, Dim>& face_normal) const {
159 : using factory_classes =
160 : typename Metavariables::factory_creation::factory_classes;
161 : call_with_dynamic_type<
162 : void, tmpl::at<factory_classes, elliptic::analytic_data::InitialGuess>>(
163 : solution_.get(), [this, &face_inertial_coords, &face_normal, &fields...,
164 : &n_dot_fluxes...](const auto* const derived) {
165 : const auto impose_boundary_condition = [this, &face_inertial_coords,
166 : &face_normal, derived](
167 : auto field_tag_v,
168 : auto flux_tag_v,
169 : const auto field,
170 : const auto n_dot_flux) {
171 : using field_tag = std::decay_t<decltype(field_tag_v)>;
172 : using flux_tag = std::decay_t<decltype(flux_tag_v)>;
173 : using derived_type = std::decay_t<decltype(*derived)>;
174 : switch (get<elliptic::Tags::BoundaryConditionType<field_tag>>(
175 : boundary_condition_types_)) {
176 : case elliptic::BoundaryConditionType::Dirichlet: {
177 : if constexpr (detail::has_boundary_variables_v<
178 : derived_type, Dim, field_tag>) {
179 : const auto solution_vars = derived->variables(
180 : face_inertial_coords, tmpl::list<field_tag>{});
181 : *field = get<field_tag>(solution_vars);
182 : } else {
183 : ERROR(
184 : "The analytic data does not provide the field required "
185 : "for this Dirichlet boundary condition.");
186 : }
187 : break;
188 : }
189 : case elliptic::BoundaryConditionType::Neumann: {
190 : if constexpr (detail::has_boundary_variables_v<derived_type,
191 : Dim, flux_tag>) {
192 : const auto solution_vars = derived->variables(
193 : face_inertial_coords, tmpl::list<flux_tag>{});
194 : normal_dot_flux(n_dot_flux, face_normal,
195 : get<flux_tag>(solution_vars));
196 : } else {
197 : ERROR(
198 : "The analytic data does not provide the flux required "
199 : "for this Neumann boundary condition.");
200 : }
201 : break;
202 : }
203 : default:
204 : ERROR("Unsupported boundary condition type: "
205 : << get<elliptic::Tags::BoundaryConditionType<field_tag>>(
206 : boundary_condition_types_));
207 : }
208 : };
209 : EXPAND_PACK_LEFT_TO_RIGHT(impose_boundary_condition(
210 : FieldTags{}, FluxTags{}, fields, n_dot_fluxes));
211 : });
212 : }
213 :
214 0 : using argument_tags_linearized = tmpl::list<>;
215 0 : using volume_tags_linearized = tmpl::list<>;
216 :
217 0 : void apply_linearized(
218 : const gsl::not_null<typename FieldTags::type*>... fields,
219 : const gsl::not_null<typename FieldTags::type*>... n_dot_fluxes,
220 : const TensorMetafunctions::prepend_spatial_index<
221 : typename FieldTags::type, Dim, UpLo::Lo,
222 : Frame::Inertial>&... /*deriv_fields*/) const {
223 : const auto impose_boundary_condition =
224 : [this](auto field_tag_v, const auto field, const auto n_dot_flux) {
225 : using field_tag = decltype(field_tag_v);
226 : switch (get<elliptic::Tags::BoundaryConditionType<field_tag>>(
227 : boundary_condition_types_)) {
228 : case elliptic::BoundaryConditionType::Dirichlet:
229 : for (auto& field_component : *field) {
230 : field_component = 0.;
231 : }
232 : break;
233 : case elliptic::BoundaryConditionType::Neumann:
234 : for (auto& n_dot_flux_component : *n_dot_flux) {
235 : n_dot_flux_component = 0.;
236 : }
237 : break;
238 : default:
239 : ERROR("Unsupported boundary condition type: "
240 : << get<elliptic::Tags::BoundaryConditionType<field_tag>>(
241 : boundary_condition_types_));
242 : }
243 : };
244 : EXPAND_PACK_LEFT_TO_RIGHT(
245 : impose_boundary_condition(FieldTags{}, fields, n_dot_fluxes));
246 : }
247 :
248 : // NOLINTNEXTLINE
249 0 : void pup(PUP::er& p) override;
250 :
251 : private:
252 0 : std::unique_ptr<elliptic::analytic_data::InitialGuess> solution_{nullptr};
253 : tuples::TaggedTuple<elliptic::Tags::BoundaryConditionType<FieldTags>...>
254 0 : boundary_condition_types_{};
255 : };
256 :
257 : template <typename System, size_t Dim, typename... FieldTags,
258 : typename... FluxTags>
259 : void AnalyticSolution<System, Dim, tmpl::list<FieldTags...>,
260 : tmpl::list<FluxTags...>>::pup(PUP::er& p) {
261 : Base::pup(p);
262 : p | solution_;
263 : p | boundary_condition_types_;
264 : }
265 :
266 : /// \cond
267 : template <typename System, size_t Dim, typename... FieldTags,
268 : typename... FluxTags>
269 : PUP::able::PUP_ID AnalyticSolution<System, Dim, tmpl::list<FieldTags...>,
270 : tmpl::list<FluxTags...>>::my_PUP_ID =
271 : 0; // NOLINT
272 : /// \endcond
273 :
274 : } // namespace elliptic::BoundaryConditions
|
