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            0 :     }
      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            0 :     }
      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            0 :     }
      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            0 :         }
     138              :     }
     139            0 : }
     140              : }