SpECTRE Documentation Coverage Report
Current view: top level - Elliptic/BoundaryConditions - AnalyticSolution.hpp Hit Total Coverage
Commit: 817e13c5144619b701c7cd870655d8dbf94ab8ce Lines: 3 27 11.1 %
Date: 2024-07-19 22:17:05
Legend: Lines: hit not hit

          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

Generated by: LCOV version 1.14