Line data    Source code

       1             : #ifndef MODEL_INTERFACE_GB_READ_H
       2             : #define MODEL_INTERFACE_GB_READ_H
       3             : 
       4             : #include <iostream>
       5             : #include <fstream>
       6             : #include <vector>
       7             : 
       8             : #include "AMReX.H"
       9             : #include "GB.H"
      10             : 
      11             : #include "Numeric/Interpolator/Linear.H"
      12             : 
      13             : #define PI 3.14159265 
      14             : 
      15             : namespace Model
      16             : {
      17             : namespace Interface
      18             : {
      19             : namespace GB
      20             : {
      21             : /// Reads the data from a file and computes energies and its derivates
      22             : ///
      23             : class Read : public GB
      24             : {
      25             : public:
      26             :     static constexpr const char* name = "read";
      27             : 
      28             : 
      29             :     /// \brief Read in data
      30             :     ///
      31             :     /// Reads the data from a file and abort if it is not possible to open the file or if the range of thetas do not give a range between 0 and 2pi. It also computes the derivatives of the energy and stores them in vectors that will be used in W, DW, and DDW
      32             :     ///
      33             :     /// \f[ \int_0^1x^2dx = \frac{1}{3} \f]
      34             :     ///
      35             :     Read() {}
      36             :     Read(IO::ParmParse &pp) {pp_queryclass(*this);}
      37           0 :     Read(IO::ParmParse &pp, std::string name) {pp_queryclass(name,*this);}
      38             :     Read(std::string filename)
      39             :     {
      40             :         Define(filename);
      41             :     }
      42             : 
      43           0 :     void Define(std::string filename)
      44             :     {
      45           0 :         m_w   = Numeric::Interpolator::Linear<Set::Scalar>::Read(filename,0);
      46           0 :         m_dw  = Numeric::Interpolator::Linear<Set::Scalar>::Read(filename,1);
      47           0 :         m_ddw = Numeric::Interpolator::Linear<Set::Scalar>::Read(filename,2);
      48           0 :         return;
      49             :         std::ifstream input;
      50             :         input.open(filename);
      51             :         std::string line;
      52             :         std::vector<Set::Scalar> theta, thetasmall, w, dw, ddw;
      53             :         while(std::getline(input,line))
      54             :         {
      55             :             std::vector<std::string> dat = Util::String::Split(line);
      56             :             theta.push_back(std::stof(dat[0]));
      57             :             w.push_back(std::stof(dat[1]));
      58             :         }
      59             :         for (unsigned int i = 1; i < theta.size()-1; i++)
      60             :         {
      61             :             thetasmall.push_back(theta[i]);
      62             :             dw.push_back((w[i+1] - w[i-1]) / (theta[i+1] - theta[i-1]));
      63             :             ddw.push_back((w[i+1] - 2.0*w[i] + w[i-1]) / ((theta[i+1]-theta[i]) * (theta[i]-theta[i-1])));
      64             :         }
      65             :         m_w.define(w,theta);
      66             :         m_dw.define(dw,thetasmall);
      67             :         m_ddw.define(ddw,thetasmall);
      68             : 
      69             :         for (Set::Scalar t = -0.001; t < 2*Set::Constant::Pi+1.0; t+=0.001)
      70             :         {
      71             :             Set::Scalar w = m_w(t), dw = m_dw(t), ddw = m_ddw(t);
      72             :             if (std::isnan(w) || std::isnan(dw) || std::isnan(ddw) ||
      73             :                 std::isinf(w) || std::isinf(dw) || std::isinf(ddw)) 
      74             :                 Util::Abort(INFO,"Error in GB Read: t=",t," w=",w," dw=",dw," ddw=",ddw);
      75             :         }
      76             :     };
      77           0 :     Set::Scalar W(const Set::Scalar theta) const
      78             :     {
      79           0 :         return m_w(theta);
      80             :     };
      81           0 :     Set::Scalar DW(const Set::Scalar theta) const
      82             :     {
      83           0 :         return m_dw(theta);
      84             :     };
      85           0 :     Set::Scalar DDW(const Set::Scalar theta) const
      86             :     {
      87           0 :         return m_ddw(theta);
      88             :     };
      89             : 
      90             : private:
      91             :     Numeric::Interpolator::Linear<Set::Scalar> m_w, m_dw, m_ddw;
      92             :       
      93             : public:
      94           0 :     static void Parse(Read & value, IO::ParmParse & pp)
      95             :     {
      96           0 :         std::string filename;
      97           0 :         pp_query_file("filename",filename); // Filename containing GB data
      98           0 :         value.Define(filename);
      99           0 :     }
     100             :       
     101             : };
     102             : }
     103             : }
     104             : }
     105             : #endif