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