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           16 :     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           36 :         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           36 :         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           36 :     }
     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            0 :         }
     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