Line data    Source code

       1             : #ifndef IC_TRIG2_H_
       2             : #define IC_TRIG2_H_
       3             : 
       4             : #include <complex>
       5             : 
       6             : #include "AMReX_Vector.H"
       7             : #include "IC/IC.H"
       8             : #include "Util/Util.H"
       9             : #include "Set/Set.H"
      10             : 
      11             : namespace IC
      12             : {
      13             : /// \brief Initialize using a trigonometric series
      14             : class Trig2 : public IC
      15             : {
      16             : public:
      17          38 :     Trig2 (amrex::Vector<amrex::Geometry> &_geom,
      18             :         Set::Scalar _alpha = 1.0,
      19             :         AMREX_D_DECL(Set::Scalar _theta1 = 0,
      20             :                 Set::Scalar _theta2 = 0,
      21             :                 Set::Scalar _theta3 = 0),
      22             :         AMREX_D_DECL(int _n1 = 0,
      23             :                 int _n2 = 0,
      24             :                 int _n3 = 0),
      25          38 :         int _dim = AMREX_SPACEDIM) :
      26          38 :         IC(_geom)
      27             :     {
      28          38 :         Define(_alpha,
      29             :                 AMREX_D_DECL(_theta1,_theta2,_theta3),
      30             :                 AMREX_D_DECL(_n1,_n2,_n3),
      31             :                 _dim);
      32          38 :     }
      33             : 
      34             : 
      35          38 :     void Define(Set::Scalar _alpha = 1.0,
      36             :             AMREX_D_DECL(Set::Scalar _theta1 = 0,
      37             :                 Set::Scalar _theta2 = 0,
      38             :                 Set::Scalar _theta3 = 0),
      39             :             AMREX_D_DECL(int _n1 = 0,
      40             :                 int _n2 = 0,
      41             :                 int _n3 = 0),
      42             :             int _dim = AMREX_SPACEDIM)
      43             :     {
      44          38 :         alpha = _alpha;
      45          38 :         AMREX_D_TERM(theta1 = _theta1;, theta2 = _theta2;, theta3 = _theta3;);
      46          38 :         AMREX_D_TERM(n1 = _n1;, n2 = _n2;, n3 = _n3;);
      47          38 :         dim = _dim;
      48          38 :     }
      49             : 
      50             : 
      51          38 :     virtual void Add(const int &lev, Set::Field<Set::Scalar> &field, Set::Scalar)
      52             :     {
      53          38 :         bool cellcentered = (field[0]->boxArray().ixType() == amrex::IndexType(amrex::IntVect::TheCellVector()));
      54             : 
      55          38 :         const std::complex<Set::Scalar> I(0.0,1.0);
      56             : 
      57          38 :         const amrex::Real AMREX_D_DECL(L1 = geom[lev].ProbHi()[0] - geom[lev].ProbLo()[0],
      58             :                             L2 = geom[lev].ProbHi()[1] - geom[lev].ProbLo()[1],
      59             :                             L3 = geom[lev].ProbHi()[2] - geom[lev].ProbLo()[2]);
      60         848 :         for (amrex::MFIter mfi(*field[lev],true); mfi.isValid(); ++mfi)
      61             :         {
      62         810 :             const amrex::Box& box = mfi.validbox();
      63         810 :             amrex::BaseFab<amrex::Real> &field_box = (*field[lev])[mfi];
      64             : 
      65     2564458 :             AMREX_D_TERM(for (int i = box.loVect()[0] - field[lev]->nGrow(); i<=box.hiVect()[0] + field[lev]->nGrow(); i++),
      66             :                     for (int j = box.loVect()[1] - field[lev]->nGrow(); j<=box.hiVect()[1] + field[lev]->nGrow(); j++),
      67             :                     for (int k = box.loVect()[2] - field[lev]->nGrow(); k<=box.hiVect()[2] + field[lev]->nGrow(); k++))
      68             :             {
      69             : 
      70             :                 Set::Scalar AMREX_D_DECL(x1,x2,x3);
      71     2393550 :                 if (cellcentered)
      72             :                 {
      73     2393550 :                     AMREX_D_TERM(x1 = geom[lev].ProbLo()[0] + ((amrex::Real)(i) + 0.5) * geom[lev].CellSize()[0];,
      74             :                             x2 = geom[lev].ProbLo()[1] + ((amrex::Real)(j) + 0.5) * geom[lev].CellSize()[1];,
      75             :                             x3 = geom[lev].ProbLo()[2] + ((amrex::Real)(k) + 0.5) * geom[lev].CellSize()[2];);
      76             :                 }
      77             :                 else
      78             :                 {
      79           0 :                     AMREX_D_TERM(x1 = geom[lev].ProbLo()[0] + ((amrex::Real)(i)) * geom[lev].CellSize()[0];,
      80             :                             x2 = geom[lev].ProbLo()[1] + ((amrex::Real)(j)) * geom[lev].CellSize()[1];,
      81             :                             x3 = geom[lev].ProbLo()[2] + ((amrex::Real)(k)) * geom[lev].CellSize()[2];);
      82             :                 }
      83             : 
      84             :                 
      85     2393550 :                 Set::Scalar trigfn = 1.0;
      86             : #if AMREX_SPACEDIM > 0
      87     2393550 :                 if (dim > 0) trigfn *= std::exp(I*(theta1 + ((amrex::Real)n1)*x1*Set::Constant::Pi / L1)).real();
      88             : 
      89             : #endif
      90             : #if AMREX_SPACEDIM > 1
      91     2393550 :                 if (dim > 1) trigfn *= std::exp(I*(theta2 + ((amrex::Real)n2)*x2*Set::Constant::Pi / L2)).real();
      92             : #endif
      93             : #if AMREX_SPACEDIM > 2
      94     2369250 :                 if (dim > 2) trigfn *= std::exp(I*(theta3 + ((amrex::Real)n3)*x3*Set::Constant::Pi / L3)).real();
      95             : #endif
      96     4787100 :                 field_box(amrex::IntVect(AMREX_D_DECL(i,j,k)),comp) += alpha * trigfn;
      97             : 
      98             :             }
      99             :         }
     100          38 :     }
     101             : 
     102             : private:
     103             :     int dim = AMREX_SPACEDIM;
     104             :     Set::Scalar alpha;
     105             :     int AMREX_D_DECL(n1, n2, n3);
     106             :     Set::Scalar AMREX_D_DECL(theta1, theta2, theta3);
     107             : };
     108             : }
     109             : #endif