Line data    Source code

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