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