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