Line data    Source code

       1             : #ifndef TEST_NUMERIC_STENCIL
       2             : #define TEST_NUMERIC_STENCIL
       3             : 
       4             : #include <AMReX.H>
       5             : 
       6             : #include "Set/Set.H"
       7             : #include "IC/Trig2.H"
       8             : #include "Numeric/Stencil.H"
       9             : 
      10             : namespace Test
      11             : {
      12             : /// Tests for the Numeric namespace classes
      13             : namespace Numeric
      14             : {
      15             : class Stencil
      16             : {
      17             : 
      18             : public:
      19           8 :     Stencil() {};
      20           2 :     ~Stencil() {};
      21             : 
      22           2 :     void Define(amrex::IntVect _ncells)
      23             :     {
      24           2 :         ncells = _ncells;
      25           2 :         amrex::RealBox rb({AMREX_D_DECL(0.,0.,0.)}, {AMREX_D_DECL(L,L,L)});
      26           2 :         amrex::Geometry::Setup(&rb, 0);
      27             :         
      28           4 :         amrex::Box domain(amrex::IntVect{AMREX_D_DECL(0,0,0)}, ncells,
      29           2 :                 amrex::IntVect::TheCellVector());
      30           2 :         geom.resize(1);
      31           2 :         grids.resize(1);
      32           2 :         dmap.resize(1);
      33           2 :         phi.resize(1);
      34           2 :         DphiNumeric.resize(1);
      35           2 :         DphiExact.resize(1);
      36             : 
      37           2 :         geom[0].define(domain);
      38           2 :         grids[0].define(domain);
      39           4 :         grids[0].maxSize(ncells[0]/2);
      40           2 :         dmap[0].define(grids[0]);
      41             : 
      42           2 :         phi        .Define(0,grids[0],dmap[0],1,nghost);
      43           2 :         DphiNumeric.Define(0,grids[0],dmap[0],1,nghost);
      44           2 :         DphiExact  .Define(0,grids[0],dmap[0],1,nghost);
      45             : 
      46           2 :         phi[0]->setVal(0.0);
      47             :         // Initial function is
      48             :         //   cos(2*pi*x/L) * cos(2*pi*y/L) * cos(2*pi*z / L)
      49           2 :         IC::Trig2 ic(geom,1.0,AMREX_D_DECL(0.0,0.0,0.0),AMREX_D_DECL(1,1,1));
      50           2 :         ic.Add(0,phi,0.0);
      51           2 :     }
      52           2 :     void Define(int _ncells)
      53             :     {
      54           2 :         Define(amrex::IntVect{AMREX_D_DECL(_ncells,_ncells,_ncells)});
      55           2 :     }
      56             : 
      57             :     template<int i,int j,int k>
      58          36 :     bool Derivative(int verbose,std::string plotfile = "")
      59             :     {
      60          36 :         const Set::Scalar tolerance = 1E-2;
      61             :         
      62          36 :         DphiExact[0]->setVal(0.0);
      63          36 :         DphiNumeric[0]->setVal(0.0);
      64             : 
      65             :         //
      66             :         // Compute exact derivative
      67             :         //
      68          36 :         Set::Scalar AMREX_D_DECL(theta1 = 0.0, theta2 = 0.0, theta3 = 0.0);
      69          36 :         Set::Scalar fac = 1.0;
      70         143 :         AMREX_D_TERM( for (int p = 0; p < i; p++) theta1 += 0.5*Set::Constant::Pi;,
      71             :                 for (int q = 0; q < j; q++) theta2 += 0.5*Set::Constant::Pi;,
      72             :                 for (int r = 0; r < k; r++) theta3 += 0.5*Set::Constant::Pi; );
      73             :                   
      74         143 :         AMREX_D_TERM( for (int p = 0; p < i; p++) fac *= 1.0*Set::Constant::Pi / L;,
      75             :                 for (int q = 0; q < j; q++) fac *= 1.0*Set::Constant::Pi / L;,
      76             :                 for (int r = 0; r < k; r++) fac *= 1.0*Set::Constant::Pi / L;);
      77          72 :         IC::Trig2 ic(geom,fac,AMREX_D_DECL(theta1,theta2,theta3),AMREX_D_DECL(1,1,1));
      78          36 :         ic.Add(0,DphiExact,0.0);
      79             : 
      80             : 
      81             :         //
      82             :         // Compute numeric derivative
      83             :         //
      84          36 :         const amrex::Real* DX = geom[0].CellSize();
      85         810 :         for ( amrex::MFIter mfi(*phi[0],amrex::TilingIfNotGPU()); mfi.isValid(); ++mfi )
      86             :         {
      87         774 :             const amrex::Box& bx = mfi.tilebox();
      88         774 :             amrex::Array4<const amrex::Real> const& Phi = phi[0]->array(mfi);
      89         774 :             amrex::Array4<amrex::Real> const& DPhi    = DphiNumeric[0]->array(mfi);
      90     1821582 :             amrex::ParallelFor (bx,[=] AMREX_GPU_DEVICE(int p, int q, int r){ DPhi(p,q,r) = ::Numeric::Stencil<Set::Scalar,i,j,k>::D(Phi,p,q,r,0,DX); });
      91             :         }
      92             : 
      93         108 :         amrex::MultiFab diff(grids[0],dmap[0],ncomp,nghost);
      94          36 :         amrex::MultiFab::Copy    (diff,*DphiExact[0]  ,0,0,ncomp,nghost);
      95          36 :         amrex::MultiFab::Subtract(diff,*DphiNumeric[0],0,0,ncomp,nghost);
      96             : 
      97          36 :         Set::Scalar error = diff.norm0(0,0,false) / DphiExact[0]->norm0(0,0,false);
      98             : 
      99          36 :         if (verbose) Util::Message(INFO,"Infinity norm of error = ", 100*error,"%");
     100             : 
     101          36 :         if (plotfile != "")
     102           0 :             WritePlotFile(plotfile);
     103             : 
     104          36 :         if (error < tolerance) return 0;
     105           0 :         else return 1;
     106             :     }
     107             : 
     108             : 
     109           0 :     void WritePlotFile(std::string plotfile)
     110             :     {
     111           0 :         amrex::Vector<amrex::MultiFab> plotmf(1);
     112           0 :         plotmf[0].define(grids[0], dmap[0], 3, 0);
     113           0 :         for (amrex::MFIter mfi(plotmf[0],false); mfi.isValid(); ++mfi)
     114             :         {
     115           0 :             amrex::Box bx = mfi.validbox();
     116           0 :             plotmf[0][mfi].copy((*phi[0])[mfi]        ,bx, 0, bx, 0, 1); 
     117           0 :             plotmf[0][mfi].copy((*DphiExact[0])[mfi]  ,bx, 0, bx, 1, 1); 
     118           0 :             plotmf[0][mfi].copy((*DphiNumeric[0])[mfi],bx, 0, bx, 2, 1); 
     119             :         }
     120           0 :         amrex::WriteSingleLevelPlotfile(plotfile,plotmf[0],varnames,geom[0],0.0,1);
     121           0 :     }        
     122             : 
     123             : 
     124             : private:
     125             :     amrex::IntVect ncells;
     126             :     const Set::Scalar L = 1.0;
     127             :     const int ncomp = 1, nghost = 2;
     128             :     amrex::Vector<amrex::Geometry> geom;
     129             :     amrex::Vector<amrex::BoxArray> grids;
     130             :     amrex::Vector<amrex::DistributionMapping> dmap;
     131             :     
     132             :     Set::Field<Set::Scalar> phi;
     133             :     Set::Field<Set::Scalar> DphiExact;
     134             :     Set::Field<Set::Scalar> DphiNumeric;
     135             :     amrex::Vector<std::string> varnames = {"phi","dphi_exact","dphi_numeric"};
     136             : };
     137             : }
     138             : }
     139             : 
     140             : 
     141             : 
     142             : #endif