Line data    Source code

       1             : #ifndef MODEL_INTERFACE_GB_SH_H
       2             : #define MODEL_INTERFACE_GB_SH_H
       3             : 
       4             : #include <iostream>
       5             : #include <fstream>
       6             : 
       7             : #include "AMReX.H"
       8             : #include "GB.H"
       9             : #include "Set/Set.H"
      10             : #include "Util/Util.H"
      11             : 
      12             : #define PI 3.14159265 
      13             : 
      14             : namespace Model
      15             : {
      16             : namespace Interface
      17             : {
      18             : namespace GB
      19             : {
      20             : ///
      21             : /// A 2D interface model class. Under construction
      22             : ///
      23             : /// Here are some useful trig identities for reference:
      24             : ///
      25             : ///     \f{align}{\sin(2 \arcsin x) &= 2 x \sqrt{1-x^2} &
      26             : ///               \cos(2 \arcsin x) &= 1 - 2 x^2 \\ %
      27             : ///               \sin(2 \arccos x) &= 2 x \sqrt{1-x^2} &
      28             : ///               \cos(2 \arccos x) &= 2 x^2 - 1 \\ %
      29             : ///               \sin(2 \arctan x) &= \frac{2 x}{1 + x^2}  &
      30             : ///               \cos(2 \arctan x) &= \frac{1-x^2}{1+x^2}    \f}
      31             : /// Specializing for tangent:
      32             : ///     \f{align}{\sin(2 \arctan y/x) &= \frac{2 x y}{x^2 + y^2}     &
      33             : ///               \cos(2 \arctan y/x) &= \frac{x^2-y^2}{x^2+y^2}  \f}
      34             : ///
      35             : class SH : public GB
      36             : {
      37             : public:
      38             :     static constexpr const char* name = "sh";
      39             : 
      40             :     SH(){};
      41             :     SH(IO::ParmParse &pp){pp_queryclass(*this);};
      42           0 :     SH(IO::ParmParse &pp,std::string name){pp_queryclass(name,*this);};
      43       12222 :     SH(const amrex::Real a_theta0, const amrex::Real a_phi0, const amrex::Real a_sigma0, const amrex::Real a_sigma1)
      44       12222 :     {
      45       12222 :         Define(a_theta0,a_phi0,a_sigma0,a_sigma1);
      46       12222 :     };
      47       12222 :     void Define(const amrex::Real a_theta0, const amrex::Real a_phi0, const amrex::Real a_sigma0, const amrex::Real a_sigma1)
      48             :     {
      49       12222 :         theta0 = a_theta0;
      50       12222 :         phi0   = a_phi0;
      51       12222 :         sigma0 = a_sigma0;
      52       12222 :         sigma1 = a_sigma1;
      53       12222 :     };
      54             :     void Randomize()
      55             :     {
      56             :         theta0 = Util::Random()*Set::Constant::Pi;
      57             :         sigma0 = Util::Random();
      58             :         sigma1 = Util::Random();
      59             :     };
      60           0 :     Set::Scalar W(const Set::Scalar a_theta, const Set::Scalar a_phi) const
      61             :     {
      62           0 :         return sigma0 + sigma1*(1.0 - (cos(2*(a_phi))*cos(2*(a_phi)) * sin(2*(a_theta))*sin(2*(a_theta))));
      63             :     };
      64             :     std::array<Set::Scalar,2> DW(const Set::Scalar a_theta, const Set::Scalar a_phi) const
      65             :     {
      66             :         return {- sigma1 * 4.0 * cos(2*a_phi) * cos(2*a_phi) * sin(2*a_theta) * cos(2*a_theta),
      67             :                 + sigma1 * 4.0 * sin(2*a_phi) * cos(2*a_phi) * sin(2*a_theta) * sin(2*a_theta)} ;
      68             :     };
      69             : 
      70           0 :     virtual Set::Scalar W(const Set::Vector &a_n) const override
      71             :     {
      72           0 :         Set::Vector n = a_n / a_n.lpNorm<2>();
      73           0 :         Set::Scalar theta = std::acos(n(0));
      74           0 :         Set::Scalar phi   = std::atan2(n(2),n(1));
      75           0 :         return W(theta,phi);
      76             :     };
      77           0 :     virtual Set::Scalar DW(const Set::Vector &a_n, const Set::Vector &t_n) const override
      78             :     {
      79           0 :         const Set::Scalar alpha = 1E-4;
      80           0 :         return (W(a_n + alpha*t_n) - W(a_n - alpha*t_n)) / 2.0 / alpha;
      81             :     };
      82           0 :     virtual Set::Scalar DDW(const Set::Vector &a_n, const Set::Vector &t_n) const override
      83             :     {
      84           0 :         const Set::Scalar alpha = 1E-4;
      85           0 :         return (W(a_n + alpha*t_n) - 2.0*W(a_n) + W(a_n - alpha*t_n)) / alpha / alpha;
      86             :     };
      87             : 
      88             : 
      89           0 :     Set::Scalar W(const Set::Scalar ) const override   {return NAN;};
      90           0 :     Set::Scalar DW(const Set::Scalar ) const override  {return NAN;};
      91           0 :     Set::Scalar DDW(const Set::Scalar ) const override {return NAN;};
      92             :   
      93             : private:
      94             :     Set::Scalar theta0 = NAN, phi0 = NAN, sigma0 = NAN, sigma1 = NAN;
      95             :     
      96             : public:
      97           0 :     static void Parse(SH & value, amrex::ParmParse & pp)
      98             :     {
      99           0 :         pp_query("theta0",value.theta0); // Theta offset (degrees)
     100           0 :         value.theta0 *= 0.01745329251; // convert degrees into radians
     101           0 :         pp_query("phi0",value.phi0); // Phi offset (radians)
     102           0 :         value.phi0 *= 0.01745329251;   // convert degrees into radians
     103           0 :         pp_query("sigma0",value.sigma0); // Minimum energy value
     104           0 :         pp_query("sigma1",value.sigma1); // Energy multiplier
     105             : 
     106           0 :         std::string reg_str = "wilhelm";
     107           0 :         pp_query("regularization",reg_str); // Type of regularization to use: {wilhelm,k23}
     108           0 :         if (reg_str == "wilhelm") value.regularization = Regularization::Wilhelm;
     109           0 :         else if (reg_str == "k23") value.regularization = Regularization::K23;
     110           0 :         else Util::Abort(INFO,"Invalid regularization string ",reg_str);
     111           0 :     }
     112             : 
     113             : };
     114             : }
     115             : }
     116             : }
     117             : #endif