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 <pup.h>
       8             : #include <string>
       9             : #include <vector>
      10             : 
      11             : #include "DataStructures/Tensor/TypeAliases.hpp"
      12             : #include "Elliptic/BoundaryConditions/BoundaryCondition.hpp"
      13             : #include "Elliptic/BoundaryConditions/BoundaryConditionType.hpp"
      14             : #include "Options/Context.hpp"
      15             : #include "Options/String.hpp"
      16             : #include "Utilities/Gsl.hpp"
      17             : #include "Utilities/Serialization/CharmPupable.hpp"
      18             : #include "Utilities/TMPL.hpp"
      19             : 
      20             : /// \cond
      21             : class DataVector;
      22             : /// \endcond
      23             : 
      24             : namespace Poisson::BoundaryConditions {
      25             : 
      26             : /// Impose Robin boundary conditions \f$a u + b n_i \nabla^i u = c\f$. The
      27             : /// boundary condition is imposed as Neumann-type (i.e. on \f$n_i \nabla^i u\f$)
      28             : /// if \f$|b| > 0\f$ and as Dirichlet-type (i.e. on \f$u\f$) if \f$b = 0\f$.
      29             : template <size_t Dim>
      30           1 : class Robin : public elliptic::BoundaryConditions::BoundaryCondition<Dim> {
      31             :  private:
      32           0 :   using Base = elliptic::BoundaryConditions::BoundaryCondition<Dim>;
      33             : 
      34             :  public:
      35           0 :   static constexpr Options::String help =
      36             :       "Robin boundary conditions a * u + b * n_i grad(u)^i = c. The boundary "
      37             :       "condition is imposed as Neumann-type (i.e. on n_i grad(u)^i) if abs(b) "
      38             :       "> 0 and as Dirichlet-type (i.e. on u) if b = 0.";
      39             : 
      40           0 :   struct DirichletWeight {
      41           0 :     using type = double;
      42           0 :     static constexpr Options::String help = "The parameter 'a'";
      43             :   };
      44             : 
      45           0 :   struct NeumannWeight {
      46           0 :     using type = double;
      47           0 :     static constexpr Options::String help = "The parameter 'b'";
      48             :   };
      49             : 
      50           0 :   struct Constant {
      51           0 :     using type = double;
      52           0 :     static constexpr Options::String help = "The parameter 'c'";
      53             :   };
      54             : 
      55           0 :   using options = tmpl::list<DirichletWeight, NeumannWeight, Constant>;
      56             : 
      57           0 :   Robin() = default;
      58           0 :   Robin(const Robin&) = default;
      59           0 :   Robin& operator=(const Robin&) = default;
      60           0 :   Robin(Robin&&) = default;
      61           0 :   Robin& operator=(Robin&&) = default;
      62           0 :   ~Robin() = default;
      63             : 
      64             :   /// \cond
      65             :   explicit Robin(CkMigrateMessage* m) : Base(m) {}
      66             :   using PUP::able::register_constructor;
      67             :   WRAPPED_PUPable_decl_template(Robin);
      68             :   /// \endcond
      69             : 
      70           0 :   std::unique_ptr<domain::BoundaryConditions::BoundaryCondition> get_clone()
      71             :       const override {
      72             :     return std::make_unique<Robin>(*this);
      73             :   }
      74             : 
      75           0 :   Robin(double dirichlet_weight, double neumann_weight, double constant,
      76             :         const Options::Context& context = {});
      77             : 
      78           0 :   double dirichlet_weight() const { return dirichlet_weight_; }
      79           0 :   double neumann_weight() const { return neumann_weight_; }
      80           0 :   double constant() const { return constant_; }
      81             : 
      82           0 :   std::vector<elliptic::BoundaryConditionType> boundary_condition_types()
      83             :       const override {
      84             :     return {1, neumann_weight_ == 0.
      85             :                    ? elliptic::BoundaryConditionType::Dirichlet
      86             :                    : elliptic::BoundaryConditionType::Neumann};
      87             :   }
      88             : 
      89           0 :   using argument_tags = tmpl::list<>;
      90           0 :   using volume_tags = tmpl::list<>;
      91             : 
      92           0 :   void apply(gsl::not_null<Scalar<DataVector>*> field,
      93             :              gsl::not_null<Scalar<DataVector>*> n_dot_field_gradient,
      94             :              const tnsr::i<DataVector, Dim>& deriv_field) const;
      95             : 
      96           0 :   using argument_tags_linearized = tmpl::list<>;
      97           0 :   using volume_tags_linearized = tmpl::list<>;
      98             : 
      99           0 :   void apply_linearized(
     100             :       gsl::not_null<Scalar<DataVector>*> field_correction,
     101             :       gsl::not_null<Scalar<DataVector>*> n_dot_field_gradient_correction,
     102             :       const tnsr::i<DataVector, Dim>& deriv_field_correction) const;
     103             : 
     104           0 :   void pup(PUP::er& p) override;
     105             : 
     106             :  private:
     107           0 :   double dirichlet_weight_ = std::numeric_limits<double>::signaling_NaN();
     108           0 :   double neumann_weight_ = std::numeric_limits<double>::signaling_NaN();
     109           0 :   double constant_ = std::numeric_limits<double>::signaling_NaN();
     110             : };
     111             : 
     112             : template <size_t Dim>
     113           0 : bool operator==(const Robin<Dim>& lhs, const Robin<Dim>& rhs);
     114             : 
     115             : template <size_t Dim>
     116           0 : bool operator!=(const Robin<Dim>& lhs, const Robin<Dim>& rhs);
     117             : 
     118             : }  // namespace Poisson::BoundaryConditions