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