LCOV - code coverage report
Current view: top level - src/IC - Trig.H (source / functions) Hit Total Coverage
Test: coverage_merged.info Lines: 33 62 53.2 %
Date: 2025-02-18 04:54:05 Functions: 4 8 50.0 %

          Line data    Source code
       1             : //
       2             : // This is an old-fashioned IC. 
       3             : // It is kept because it is still used with the Trig regression test, but 
       4             : // it is recommended to use IC::Expression instead.
       5             : //
       6             : 
       7             : #ifndef IC_TRIG_H_
       8             : #define IC_TRIG_H_
       9             : 
      10             : #include <complex>
      11             : 
      12             : #include "IO/ParmParse.H"
      13             : #include "AMReX_Vector.H"
      14             : #include "IC/IC.H"
      15             : #include "Util/Util.H"
      16             : #include "Set/Set.H"
      17             : 
      18             : namespace IC
      19             : {
      20             : /// \brief Initialize using a trigonometric series
      21             : class Trig : public IC
      22             : {
      23             : public:
      24             :     static constexpr const char* name = "trig";
      25             :     Trig (amrex::Vector<amrex::Geometry> &_geom,
      26             :         Set::Scalar _alpha = 1.0,
      27             :         AMREX_D_DECL(std::complex<int> _n1 = 0,
      28             :                 std::complex<int> _n2 = 0,
      29             :                 std::complex<int> _n3 = 0),
      30             :         int _dim = AMREX_SPACEDIM) :
      31             :         IC(_geom)
      32             :     {
      33             :         Define(_alpha,AMREX_D_DECL(_n1,_n2,_n3),_dim);
      34             :     }
      35           4 :     Trig(amrex::Vector<amrex::Geometry> &_geom,IO::ParmParse &pp, std::string name) : IC(_geom)
      36           4 :     {pp_queryclass(name,*this);}
      37             : 
      38             : 
      39             :     void Define(Set::Scalar _alpha = 1.0,
      40             :             AMREX_D_DECL(std::complex<int> _n1 = 0,
      41             :                 std::complex<int> _n2 = 0,
      42             :                 std::complex<int> _n3 = 0),
      43             :             int _dim = AMREX_SPACEDIM)
      44             :     {
      45             :         alpha = _alpha;
      46             :         AMREX_D_TERM(n1 = _n1;, n2 = _n2;, n3 = _n3;);
      47             :         dim = _dim;
      48             :         AMREX_D_DECL(phi1 = std::atan2(n1.imag(),n1.real()),
      49             :                 phi2 = std::atan2(n2.imag(),n2.real()),
      50             :                 phi3 = std::atan2(n3.imag(),n3.real()));
      51             :     }
      52             : 
      53             : 
      54           0 :     void Add(const int &lev, Set::Field<Set::Scalar> &a_field, Set::Scalar)
      55             :     {
      56           0 :         bool cellcentered = (a_field[0]->boxArray().ixType() == amrex::IndexType(amrex::IntVect::TheCellVector()));
      57             : 
      58           0 :         const Set::Scalar AMREX_D_DECL(L1 = geom[lev].ProbHi()[0] - geom[lev].ProbLo()[0],
      59             :                             L2 = geom[lev].ProbHi()[1] - geom[lev].ProbLo()[1],
      60             :                             L3 = geom[lev].ProbHi()[2] - geom[lev].ProbLo()[2]);
      61             : 
      62             : 
      63           0 :         for (amrex::MFIter mfi(*a_field[lev], amrex::TilingIfNotGPU()); mfi.isValid(); ++mfi)
      64             :         {
      65           0 :             amrex::Box bx;
      66           0 :             if (cellcentered) bx = mfi.growntilebox();
      67           0 :             else bx = mfi.grownnodaltilebox();
      68             : 
      69           0 :             amrex::Array4<Set::Scalar> const &field = a_field[lev]->array(mfi);
      70           0 :             amrex::ParallelFor(bx, [=] AMREX_GPU_DEVICE(int i, int j, int k)
      71             :             {
      72             :                 Set::Scalar AMREX_D_DECL(x1,x2,x3);
      73           0 :                 if (cellcentered)
      74             :                 {
      75           0 :                     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           0 :                 Set::Scalar trigfn = 1.0;
      87             :                 #if AMREX_SPACEDIM > 0
      88           0 :                 if (dim > 0)
      89           0 :                     trigfn *= (fabs(std::cos(phi1))*std::cos(n1.real()*Set::Constant::Pi*x1 / L1) +
      90           0 :                             fabs(std::sin(phi1))*std::sin(n1.imag()*Set::Constant::Pi*x1 / L1));
      91             :                 #endif
      92             :                 #if AMREX_SPACEDIM > 1
      93           0 :                 if (dim > 1)
      94           0 :                     trigfn *= (fabs(std::cos(phi2))*std::cos(n2.real()*Set::Constant::Pi*x2 / L2) +
      95           0 :                             fabs(std::sin(phi2))*std::sin(n2.imag()*Set::Constant::Pi*x2 / L2));
      96             :                 #endif
      97             :                 #if AMREX_SPACEDIM > 2
      98           0 :                 if (dim > 2)
      99           0 :                     trigfn *= (fabs(std::cos(phi3))*std::cos(n3.real()*Set::Constant::Pi*x3 / L3) +
     100           0 :                             fabs(std::sin(phi3))*std::sin(n3.imag()*Set::Constant::Pi*x3 / L3));
     101             :                 #endif
     102             :                         
     103           0 :                 field(i,j,k,comp) += alpha * trigfn;
     104             : 
     105           0 :             });
     106             :         }
     107           0 :     }
     108          14 :     void Add(const int &lev, Set::Field<Set::Vector> &a_field, Set::Scalar)
     109             :     {
     110          14 :         bool cellcentered = (a_field[0]->boxArray().ixType() == amrex::IndexType(amrex::IntVect::TheCellVector()));
     111             : 
     112          14 :         const Set::Scalar AMREX_D_DECL(L1 = geom[lev].ProbHi()[0] - geom[lev].ProbLo()[0],
     113             :                             L2 = geom[lev].ProbHi()[1] - geom[lev].ProbLo()[1],
     114             :                             L3 = geom[lev].ProbHi()[2] - geom[lev].ProbLo()[2]);
     115             : 
     116             : 
     117          28 :         for (amrex::MFIter mfi(*a_field[lev], amrex::TilingIfNotGPU()); mfi.isValid(); ++mfi)
     118             :         {
     119          14 :             amrex::Box bx;
     120          14 :             if (cellcentered) bx = mfi.growntilebox();
     121          14 :             else bx = mfi.grownnodaltilebox();
     122             : 
     123          14 :             amrex::Array4<Set::Vector> const &field = a_field[lev]->array(mfi);
     124          14 :             amrex::ParallelFor(bx, [=] AMREX_GPU_DEVICE(int i, int j, int k)
     125             :             {
     126             :                 Set::Scalar AMREX_D_DECL(x1,x2,x3);
     127       19166 :                 if (cellcentered)
     128             :                 {
     129           0 :                     AMREX_D_TERM(x1 = geom[lev].ProbLo()[0] + ((amrex::Real)(i) + 0.5) * geom[lev].CellSize()[0];,
     130             :                                 x2 = geom[lev].ProbLo()[1] + ((amrex::Real)(j) + 0.5) * geom[lev].CellSize()[1];,
     131             :                                 x3 = geom[lev].ProbLo()[2] + ((amrex::Real)(k) + 0.5) * geom[lev].CellSize()[2];);
     132             :                 }
     133             :                 else
     134             :                 {
     135       19166 :                     AMREX_D_TERM(x1 = geom[lev].ProbLo()[0] + (amrex::Real)(i) * geom[lev].CellSize()[0];,
     136             :                                 x2 = geom[lev].ProbLo()[1] + (amrex::Real)(j) * geom[lev].CellSize()[1];,
     137             :                                 x3 = geom[lev].ProbLo()[2] + (amrex::Real)(k) * geom[lev].CellSize()[2];);
     138             :                 }
     139             :                                    
     140       19166 :                 Set::Scalar trigfn = 1.0;
     141             :                 #if AMREX_SPACEDIM > 0
     142       19166 :                 if (dim > 0)
     143       19166 :                     trigfn *= (fabs(std::cos(phi1))*std::cos(n1.real()*Set::Constant::Pi*x1 / L1) +
     144       19166 :                             fabs(std::sin(phi1))*std::sin(n1.imag()*Set::Constant::Pi*x1 / L1));
     145             :                 #endif
     146             :                 #if AMREX_SPACEDIM > 1
     147       19166 :                 if (dim > 1)
     148       19166 :                     trigfn *= (fabs(std::cos(phi2))*std::cos(n2.real()*Set::Constant::Pi*x2 / L2) +
     149       19166 :                             fabs(std::sin(phi2))*std::sin(n2.imag()*Set::Constant::Pi*x2 / L2));
     150             :                 #endif
     151             :                 #if AMREX_SPACEDIM > 2
     152           0 :                 if (dim > 2)
     153           0 :                     trigfn *= (fabs(std::cos(phi3))*std::cos(n3.real()*Set::Constant::Pi*x3 / L3) +
     154           0 :                             fabs(std::sin(phi3))*std::sin(n3.imag()*Set::Constant::Pi*x3 / L3));
     155             :                 #endif
     156             :                         
     157       38332 :                 field(i,j,k)(comp) += alpha * trigfn;
     158             : 
     159       19166 :             });
     160             :         }
     161          14 :     }
     162             :     using IC::Add;
     163             : 
     164             : public:
     165           4 :     static void Parse(Trig & value, IO::ParmParse & pp)
     166             :     {
     167           8 :         std::vector<int> n_real(AMREX_SPACEDIM,0.0);
     168           4 :         pp_queryarr("nr",n_real); // Number of real (cosin) waves
     169           8 :         std::vector<int> n_imag(AMREX_SPACEDIM,0.0);
     170           4 :         pp_queryarr("ni",n_imag); // Number of imaginary (sin) waves
     171             : 
     172           4 :         AMREX_D_TERM(
     173             :             value.n1 = std::complex<int>(n_real[0],n_imag[0]);,
     174             :             value.n2 = std::complex<int>(n_real[1],n_imag[1]);,
     175             :             value.n3 = std::complex<int>(n_real[2],n_imag[2]););
     176             : 
     177           4 :         AMREX_D_DECL(
     178             :             value.phi1 = std::atan2(value.n1.imag(),value.n1.real()),
     179             :             value.phi2 = std::atan2(value.n2.imag(),value.n2.real()),
     180             :             value.phi3 = std::atan2(value.n3.imag(),value.n3.real()));
     181             : 
     182           4 :         pp_query("dim",value.dim); // Spatial dimension
     183           4 :         pp_query("alpha",value.alpha); // Multiplier
     184             : 
     185           4 :     }
     186             : 
     187             : private:
     188             :     int dim = AMREX_SPACEDIM;
     189             :     Set::Scalar alpha = 1.0;
     190             :     std::complex<int> AMREX_D_DECL(n1, n2, n3);
     191             :     Set::Scalar AMREX_D_DECL(phi1=0.0,phi2=0.0,phi3=0.0);
     192             : };
     193             : }
     194             : #endif

Generated by: LCOV version 1.14