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 <vector>
12 :
13 : #include "DataStructures/Tensor/Tensor.hpp"
14 : #include "Domain/Tags.hpp"
15 : #include "Domain/Tags/FaceNormal.hpp"
16 : #include "Elliptic/BoundaryConditions/BoundaryCondition.hpp"
17 : #include "Elliptic/BoundaryConditions/BoundaryConditionType.hpp"
18 : #include "Elliptic/BoundaryConditions/Tags.hpp"
19 : #include "NumericalAlgorithms/DiscontinuousGalerkin/NormalDotFlux.hpp"
20 : #include "Options/String.hpp"
21 : #include "Parallel/Tags/Metavariables.hpp"
22 : #include "PointwiseFunctions/InitialDataUtilities/AnalyticSolution.hpp"
23 : #include "Utilities/CallWithDynamicType.hpp"
24 : #include "Utilities/ErrorHandling/Error.hpp"
25 : #include "Utilities/Gsl.hpp"
26 : #include "Utilities/Serialization/CharmPupable.hpp"
27 : #include "Utilities/TMPL.hpp"
28 : #include "Utilities/TaggedTuple.hpp"
29 :
30 1 : namespace elliptic::BoundaryConditions {
31 :
32 : /// \cond
33 : template <typename System, size_t Dim = System::volume_dim,
34 : typename FieldTags = typename System::primal_fields,
35 : typename FluxTags = typename System::primal_fluxes>
36 : struct AnalyticSolution;
37 : /// \endcond
38 :
39 : /*!
40 : * \brief Impose the analytic solution on the boundary. Works only if an
41 : * analytic solution exists.
42 : *
43 : * The analytic solution is retrieved from `::Tags::AnalyticSolutionsBase`. It
44 : * must hold solutions for both the `System::primal_fields` and the
45 : * `System::primal_fluxes`. The user can select to impose the analytic solution
46 : * as Dirichlet or Neumann boundary conditions for each field separately.
47 : * Dirichlet boundary conditions are imposed on the fields and Neumann boundary
48 : * conditions are imposed on the fluxes.
49 : *
50 : * See `elliptic::Actions::InitializeAnalyticSolutions` for an action that can
51 : * add the analytic solutions to the DataBox.
52 : */
53 : template <typename System, size_t Dim, typename... FieldTags,
54 : typename... FluxTags>
55 1 : class AnalyticSolution<System, Dim, tmpl::list<FieldTags...>,
56 : tmpl::list<FluxTags...>>
57 : : public BoundaryCondition<Dim> {
58 : private:
59 0 : using Base = BoundaryCondition<Dim>;
60 :
61 : public:
62 0 : struct Solution {
63 0 : using type = std::unique_ptr<elliptic::analytic_data::AnalyticSolution>;
64 0 : static constexpr Options::String help = {
65 : "The analytic solution to impose on the boundary"};
66 : };
67 :
68 0 : using options =
69 : tmpl::list<Solution,
70 : elliptic::OptionTags::BoundaryConditionType<FieldTags>...>;
71 0 : static constexpr Options::String help =
72 : "Boundary conditions from the analytic solution";
73 :
74 0 : AnalyticSolution() = default;
75 0 : AnalyticSolution(const AnalyticSolution& rhs) : Base(rhs) { *this = rhs; }
76 0 : AnalyticSolution& operator=(const AnalyticSolution& rhs) {
77 : if (rhs.solution_ != nullptr) {
78 : solution_ = rhs.solution_->get_clone();
79 : } else {
80 : solution_ = nullptr;
81 : }
82 : boundary_condition_types_ = rhs.boundary_condition_types_;
83 : return *this;
84 : }
85 0 : AnalyticSolution(AnalyticSolution&&) = default;
86 0 : AnalyticSolution& operator=(AnalyticSolution&&) = default;
87 0 : ~AnalyticSolution() = default;
88 :
89 : /// \cond
90 : explicit AnalyticSolution(CkMigrateMessage* m) : Base(m) {}
91 : using PUP::able::register_constructor;
92 : WRAPPED_PUPable_decl_template(AnalyticSolution);
93 : /// \endcond
94 :
95 : /// Select which `elliptic::BoundaryConditionType` to apply for each field
96 1 : explicit AnalyticSolution(
97 : std::unique_ptr<elliptic::analytic_data::AnalyticSolution> solution,
98 : // This pack expansion repeats the type `elliptic::BoundaryConditionType`
99 : // for each system field
100 : const typename elliptic::OptionTags::BoundaryConditionType<
101 : FieldTags>::type... boundary_condition_types)
102 : : solution_(std::move(solution)),
103 : boundary_condition_types_{boundary_condition_types...} {}
104 :
105 0 : std::unique_ptr<domain::BoundaryConditions::BoundaryCondition> get_clone()
106 : const override {
107 : return std::make_unique<AnalyticSolution>(*this);
108 : }
109 :
110 0 : std::vector<elliptic::BoundaryConditionType> boundary_condition_types()
111 : const override {
112 : std::vector<elliptic::BoundaryConditionType> result{};
113 : const auto collect = [&result](
114 : const auto tag_v,
115 : const elliptic::BoundaryConditionType bc_type) {
116 : using tag = std::decay_t<decltype(tag_v)>;
117 : for (size_t i = 0; i < tag::type::size(); ++i) {
118 : result.push_back(bc_type);
119 : }
120 : };
121 : EXPAND_PACK_LEFT_TO_RIGHT(collect(
122 : FieldTags{}, get<elliptic::Tags::BoundaryConditionType<FieldTags>>(
123 : boundary_condition_types_)));
124 : return result;
125 : }
126 :
127 0 : using argument_tags =
128 : tmpl::list<Parallel::Tags::Metavariables,
129 : domain::Tags::Coordinates<Dim, Frame::Inertial>,
130 : ::Tags::Normalized<domain::Tags::UnnormalizedFaceNormal<
131 : Dim, Frame::Inertial>>>;
132 0 : using volume_tags = tmpl::list<Parallel::Tags::Metavariables>;
133 :
134 : template <typename Metavariables>
135 0 : void apply(const gsl::not_null<typename FieldTags::type*>... fields,
136 : const gsl::not_null<typename FieldTags::type*>... n_dot_fluxes,
137 : const TensorMetafunctions::prepend_spatial_index<
138 : typename FieldTags::type, Dim, UpLo::Lo,
139 : Frame::Inertial>&... /*deriv_fields*/,
140 : const Metavariables& /*meta*/,
141 : const tnsr::I<DataVector, Dim>& face_inertial_coords,
142 : const tnsr::i<DataVector, Dim>& face_normal) const {
143 : // Retrieve variables for both Dirichlet and Neumann conditions, then decide
144 : // which to impose. We could also retrieve either the field for the flux for
145 : // each field individually based on the selection, but that would incur the
146 : // overhead of calling into the analytic solution multiple times and
147 : // possibly computing temporary quantities multiple times. This performance
148 : // consideration is probably irrelevant because the boundary conditions are
149 : // only evaluated once at the beginning of the solve.
150 : using analytic_tags = tmpl::list<FieldTags..., FluxTags...>;
151 : using factory_classes =
152 : typename Metavariables::factory_creation::factory_classes;
153 : const auto solution_vars = call_with_dynamic_type<
154 : tuples::tagged_tuple_from_typelist<analytic_tags>,
155 : tmpl::at<factory_classes, elliptic::analytic_data::AnalyticSolution>>(
156 : solution_.get(), [&face_inertial_coords](const auto* const derived) {
157 : return derived->variables(face_inertial_coords, analytic_tags{});
158 : });
159 : const auto impose_boundary_condition = [this, &solution_vars, &face_normal](
160 : auto field_tag_v,
161 : auto flux_tag_v,
162 : const auto field,
163 : const auto n_dot_flux) {
164 : using field_tag = decltype(field_tag_v);
165 : using flux_tag = decltype(flux_tag_v);
166 : switch (get<elliptic::Tags::BoundaryConditionType<field_tag>>(
167 : boundary_condition_types_)) {
168 : case elliptic::BoundaryConditionType::Dirichlet:
169 : *field = get<field_tag>(solution_vars);
170 : break;
171 : case elliptic::BoundaryConditionType::Neumann:
172 : normal_dot_flux(n_dot_flux, face_normal,
173 : get<flux_tag>(solution_vars));
174 : break;
175 : default:
176 : ERROR("Unsupported boundary condition type: "
177 : << get<elliptic::Tags::BoundaryConditionType<field_tag>>(
178 : boundary_condition_types_));
179 : }
180 : };
181 : EXPAND_PACK_LEFT_TO_RIGHT(impose_boundary_condition(FieldTags{}, FluxTags{},
182 : fields, n_dot_fluxes));
183 : }
184 :
185 0 : using argument_tags_linearized = tmpl::list<>;
186 0 : using volume_tags_linearized = tmpl::list<>;
187 :
188 0 : void apply_linearized(
189 : const gsl::not_null<typename FieldTags::type*>... fields,
190 : const gsl::not_null<typename FieldTags::type*>... n_dot_fluxes,
191 : const TensorMetafunctions::prepend_spatial_index<
192 : typename FieldTags::type, Dim, UpLo::Lo,
193 : Frame::Inertial>&... /*deriv_fields*/) const {
194 : const auto impose_boundary_condition = [this](auto field_tag_v,
195 : const auto field,
196 : const auto n_dot_flux) {
197 : using field_tag = decltype(field_tag_v);
198 : switch (get<elliptic::Tags::BoundaryConditionType<field_tag>>(
199 : boundary_condition_types_)) {
200 : case elliptic::BoundaryConditionType::Dirichlet:
201 : for (auto& field_component : *field) {
202 : field_component = 0.;
203 : }
204 : break;
205 : case elliptic::BoundaryConditionType::Neumann:
206 : for (auto& n_dot_flux_component : *n_dot_flux) {
207 : n_dot_flux_component = 0.;
208 : }
209 : break;
210 : default:
211 : ERROR("Unsupported boundary condition type: "
212 : << get<elliptic::Tags::BoundaryConditionType<field_tag>>(
213 : boundary_condition_types_));
214 : }
215 : };
216 : EXPAND_PACK_LEFT_TO_RIGHT(
217 : impose_boundary_condition(FieldTags{}, fields, n_dot_fluxes));
218 : }
219 :
220 : // NOLINTNEXTLINE
221 0 : void pup(PUP::er& p) override;
222 :
223 : private:
224 0 : std::unique_ptr<elliptic::analytic_data::AnalyticSolution> solution_{nullptr};
225 : tuples::TaggedTuple<elliptic::Tags::BoundaryConditionType<FieldTags>...>
226 0 : boundary_condition_types_{};
227 : };
228 :
229 : template <typename System, size_t Dim, typename... FieldTags,
230 : typename... FluxTags>
231 : void AnalyticSolution<System, Dim, tmpl::list<FieldTags...>,
232 : tmpl::list<FluxTags...>>::pup(PUP::er& p) {
233 : Base::pup(p);
234 : p | solution_;
235 : p | boundary_condition_types_;
236 : }
237 :
238 : /// \cond
239 : template <typename System, size_t Dim, typename... FieldTags,
240 : typename... FluxTags>
241 : PUP::able::PUP_ID AnalyticSolution<System, Dim, tmpl::list<FieldTags...>,
242 : tmpl::list<FluxTags...>>::my_PUP_ID =
243 : 0; // NOLINT
244 : /// \endcond
245 :
246 : } // namespace elliptic::BoundaryConditions
|
