Line data    Source code

       1             : #include <AMReX_MultiFabUtil.H>
       2             : #include <AMReX_REAL.H>
       3             : #include <AMReX_MLCGSolver.H>
       4             : #include "Set/Set.H"
       5             : 
       6             : #include <AMReX_ArrayLim.H>
       7             : 
       8             : #include "Set/Set.H"
       9             : #include "Diagonal.H"
      10             : 
      11             : 
      12             : 
      13             : namespace Operator
      14             : {
      15           0 : Diagonal::Diagonal (const Vector<Geometry>& a_geom,
      16             :     const Vector<BoxArray>& a_grids,
      17             :     const Vector<DistributionMapping>& a_dmap,
      18           0 :     const LPInfo& a_info)
      19           0 : {define(a_geom, a_grids, a_dmap, a_info);}
      20           0 : Diagonal::~Diagonal () {}
      21             : 
      22             : 
      23             : void
      24           0 : Diagonal::Fapply (int amrlev,
      25             :         int mglev,
      26             :         MultiFab& f,
      27             :         const MultiFab& u
      28             :         ) const
      29             : {
      30             :     BL_PROFILE("Operator::Diagonal::Diagonal::Fapply()");
      31           0 :     std::cout << "in fapply, amrlev = " << amrlev << std::endl; 
      32             : 
      33           0 :     amrex::Box domain(m_geom[amrlev][mglev].Domain());
      34             : 
      35           0 :     for (MFIter mfi(f, true); mfi.isValid(); ++mfi)
      36             :     {
      37           0 :         const amrex::FArrayBox &ufab    = u[mfi];
      38           0 :         if(ufab.contains_inf()) Util::Abort(INFO, "Inf in ufab [before update]");
      39           0 :         if(ufab.contains_nan()) Util::Abort(INFO, "Nan in ufab [before update]");
      40             :     }
      41             : 
      42           0 :     for (MFIter mfi(f, true); mfi.isValid(); ++mfi)
      43             :     {
      44           0 :         const Box& bx = mfi.tilebox();
      45           0 :         const amrex::FArrayBox &ufab    = u[mfi];
      46           0 :         amrex::FArrayBox       &ffab    = f[mfi];
      47             : 
      48           0 :         AMREX_D_TERM(for (int m1 = bx.loVect()[0]; m1<=bx.hiVect()[0]; m1++),
      49             :                 for (int m2 = bx.loVect()[1]; m2<=bx.hiVect()[1]; m2++),
      50             :                 for (int m3 = bx.loVect()[2]; m3<=bx.hiVect()[2]; m3++))
      51             :         {
      52           0 :             amrex::IntVect m(AMREX_D_DECL(m1,m2,m3));
      53           0 :             for (int i=0; i<AMREX_SPACEDIM; i++)
      54             :             {
      55           0 :                 if (AMREX_D_TERM(m1 == domain.loVect()[0],
      56             :                         || m2 == domain.loVect()[1],
      57           0 :                         || m3 == domain.loVect()[2]) ||
      58           0 :                     AMREX_D_TERM(m1 == domain.hiVect()[0]+1,
      59             :                         || m2 == domain.hiVect()[1]+1,
      60             :                         || m3 == domain.hiVect()[2]+1))
      61             :                 {
      62           0 :                     ffab(m,i) = 2*ufab(m,i);
      63             :                 }
      64           0 :                 else if (AMREX_D_TERM(m1 == bx.loVect()[0], 
      65             :                             || m2 == bx.loVect()[1],
      66           0 :                             || m3 == bx.loVect()[2]) || 
      67           0 :                     AMREX_D_TERM(m1 == bx.hiVect()[0],
      68             :                             || m2 == bx.hiVect()[1],
      69             :                             || m3 == bx.hiVect()[2]))
      70           0 :                     continue;
      71             :                 else 
      72           0 :                     ffab(m,i) = 2*ufab(m,i);
      73             :             }
      74             :         }
      75             :     }
      76           0 :     for (MFIter mfi(f, true); mfi.isValid(); ++mfi)
      77             :     {
      78           0 :         amrex::FArrayBox       &ffab    = f[mfi];
      79           0 :         if(ffab.contains_inf()) Util::Abort(INFO, "Inf in ffab [after update]");
      80           0 :         if(ffab.contains_nan()) Util::Abort(INFO, "Nan in ffab [after update]");
      81             :     }
      82           0 : }
      83             : 
      84             : 
      85             : void
      86           0 : Diagonal::Fsmooth (int amrlev,
      87             :         int mglev,  
      88             :         MultiFab& u, 
      89             :         const MultiFab& rhs
      90             :         ) const
      91             : {
      92           0 :     for (int redblack = 0; redblack < 2; redblack++)
      93             :     {
      94             :         BL_PROFILE("Operator::Diagonal::Diagonal::Fsmooth()");
      95             : 
      96           0 :         amrex::Box domain(m_geom[amrlev][mglev].Domain());
      97             : 
      98           0 :         for (MFIter mfi(u,MFItInfo().EnableTiling().SetDynamic(true));
      99           0 :             mfi.isValid(); ++mfi)
     100             :         {
     101           0 :             const Box&       bx     = mfi.tilebox();
     102           0 :             FArrayBox&       ufab    = u[mfi];
     103           0 :             const FArrayBox& rhsfab  = rhs[mfi];
     104             : 
     105           0 :             AMREX_D_TERM(for (int m1 = bx.loVect()[0]; m1<=bx.hiVect()[0]; m1++),
     106             :                     for (int m2 = bx.loVect()[1]; m2<=bx.hiVect()[1]; m2++),
     107             :                     for (int m3 = bx.loVect()[2]; m3<=bx.hiVect()[2]; m3++))
     108             :             {
     109           0 :                 if ((AMREX_D_TERM(m1, + m2, + m3))%2 == redblack) continue;
     110           0 :                 amrex::IntVect m(AMREX_D_DECL(m1,m2,m3));
     111             : 
     112           0 :                 for (int i=0; i<AMREX_SPACEDIM; i++)
     113             :                 {
     114           0 :                     for (int k=0; k<AMREX_SPACEDIM; k++)
     115             :                     {
     116           0 :                         if (AMREX_D_TERM(m1 == domain.loVect()[0],
     117             :                                 || m2 == domain.loVect()[1],
     118           0 :                                 || m3 == domain.loVect()[2]) ||
     119           0 :                             AMREX_D_TERM(m1 == domain.hiVect()[0]+1,
     120             :                                 || m2 == domain.hiVect()[1]+1,
     121             :                                 || m3 == domain.hiVect()[2]+1))
     122             :                         {
     123           0 :                             ufab(m,k) = rhsfab(m,k) / 2.0;
     124             :                         }
     125           0 :                         else if (AMREX_D_TERM(m1 == bx.loVect()[0],
     126             :                                     || m2 == bx.loVect()[1],
     127           0 :                                     || m3 == bx.loVect()[2]) || 
     128           0 :                             AMREX_D_TERM(m1 == bx.hiVect()[0], 
     129             :                                     || m2 == bx.hiVect()[1],
     130             :                                     || m3 == bx.hiVect()[2]))
     131           0 :                             continue;
     132             :                         else 
     133           0 :                             ufab(m,k) = rhsfab(m,k) / 2.0;
     134             :                     }
     135             :                 }
     136             :             }
     137             :         }
     138             :     }
     139           0 : }
     140             : }