LCOV - code coverage report
Current view: top level - src/Operator - Operator.cpp (source / functions) Coverage Total Hit
Test: coverage_merged.info Lines: 60.4 % 462 279
Test Date: 2025-04-03 04:02:21 Functions: 52.8 % 36 19

            Line data    Source code
       1              : 
       2              : #include <AMReX_MLNodeLinOp.H>
       3              : #include <AMReX_MLCellLinOp.H>
       4              : #include <AMReX_MLNodeLap_K.H>
       5              : #include <AMReX_MultiFabUtil.H>
       6              : #include "Util/Color.H"
       7              : #include "Set/Set.H"
       8              : #include "Operator.H"
       9              : 
      10              : using namespace amrex;
      11              : namespace Operator {
      12              : 
      13              : // constexpr amrex::IntVect AMREX_D_DECL(Operator<Grid::Node>::dx,Operator<Grid::Node>::dy,Operator<Grid::Node>::dz);
      14              : constexpr amrex::IntVect AMREX_D_DECL(Operator<Grid::Cell>::dx, Operator<Grid::Cell>::dy, Operator<Grid::Cell>::dz);
      15              : 
      16          426 : void Operator<Grid::Node>::Diagonal(bool recompute)
      17              : {
      18              :     BL_PROFILE(Color::FG::Yellow + "Operator::Diagonal()" + Color::Reset);
      19          426 :     if (!recompute && m_diagonal_computed) return;
      20          426 :     m_diagonal_computed = true;
      21              : 
      22          898 :     for (int amrlev = 0; amrlev < m_num_amr_levels; ++amrlev)
      23              :     {
      24         1435 :         for (int mglev = 0; mglev < m_num_mg_levels[amrlev]; ++mglev)
      25              :         {
      26          963 :             Diagonal(amrlev, mglev, *m_diag[amrlev][mglev]);
      27              :         }
      28              :     }
      29              : }
      30              : 
      31            0 : void Operator<Grid::Node>::Diagonal(int amrlev, int mglev, amrex::MultiFab& diag)
      32              : {
      33              :     BL_PROFILE("Operator::Diagonal()");
      34              :     //Util::Message(INFO);
      35              : 
      36            0 :     int ncomp = diag.nComp();
      37            0 :     int nghost = 0;
      38              : 
      39            0 :     int sep = 2;
      40            0 :     int num = AMREX_D_TERM(sep, *sep, *sep);
      41            0 :     int cntr = 0;
      42              : 
      43            0 :     amrex::MultiFab x(m_diag[amrlev][mglev]->boxArray(), m_diag[amrlev][mglev]->DistributionMap(), ncomp, nghost);
      44            0 :     amrex::MultiFab Ax(m_diag[amrlev][mglev]->boxArray(), m_diag[amrlev][mglev]->DistributionMap(), ncomp, nghost);
      45              : 
      46            0 :     for (MFIter mfi(x, false); mfi.isValid(); ++mfi)
      47              :     {
      48            0 :         const Box& bx = mfi.validbox();
      49            0 :         amrex::FArrayBox& diagfab = diag[mfi];
      50            0 :         amrex::FArrayBox& xfab = x[mfi];
      51            0 :         amrex::FArrayBox& Axfab = Ax[mfi];
      52              : 
      53            0 :         diagfab.setVal(0.0);
      54              : 
      55            0 :         for (int i = 0; i < num; i++)
      56              :         {
      57            0 :             for (int n = 0; n < ncomp; n++)
      58              :             {
      59            0 :                 xfab.setVal(0.0);
      60            0 :                 Axfab.setVal(0.0);
      61              : 
      62              :                 //BL_PROFILE_VAR("Operator::Part1", part1); 
      63            0 :                 AMREX_D_TERM(for (int m1 = bx.loVect()[0]; m1 <= bx.hiVect()[0]; m1++),
      64              :                     for (int m2 = bx.loVect()[1]; m2 <= bx.hiVect()[1]; m2++),
      65              :                         for (int m3 = bx.loVect()[2]; m3 <= bx.hiVect()[2]; m3++))
      66              :                 {
      67            0 :                     amrex::IntVect m(AMREX_D_DECL(m1, m2, m3));
      68              : 
      69            0 :                     if (m1 % sep == i / sep && m2 % sep == i % sep) xfab(m, n) = 1.0;
      70            0 :                     else xfab(m, n) = 0.0;
      71              :                 }
      72              :                 //BL_PROFILE_VAR_STOP(part1);
      73              : 
      74              :                 BL_PROFILE_VAR("Operator::Part2", part2);
      75            0 :                 Util::Message(INFO, "Calling fapply...", cntr++);
      76            0 :                 Fapply(amrlev, mglev, Ax, x);
      77              :                 BL_PROFILE_VAR_STOP(part2);
      78              : 
      79              :                 //BL_PROFILE_VAR("Operator::Part3", part3); 
      80            0 :                 Axfab.mult(xfab, n, n, 1);
      81            0 :                 diagfab.plus(Axfab, n, n, 1);
      82              :                 //BL_PROFILE_VAR_STOP(part3);
      83              :             }
      84              :         }
      85            0 :     }
      86            0 : }
      87              : 
      88        39500 : void Operator<Grid::Node>::Fsmooth(int amrlev, int mglev, amrex::MultiFab& x, const amrex::MultiFab& b) const
      89              : {
      90              :     BL_PROFILE("Operator::Fsmooth()");
      91              : 
      92        39500 :     amrex::Box domain(m_geom[amrlev][mglev].Domain());
      93              : 
      94        39500 :     int ncomp = b.nComp();
      95        39500 :     int nghost = 2; //b.nGrow();
      96              : 
      97              : 
      98        39500 :     amrex::MultiFab Ax(x.boxArray(), x.DistributionMap(), ncomp, nghost);
      99        39500 :     amrex::MultiFab Dx(x.boxArray(), x.DistributionMap(), ncomp, nghost);
     100        39500 :     amrex::MultiFab Rx(x.boxArray(), x.DistributionMap(), ncomp, nghost);
     101              : 
     102        39500 :     if (!m_diagonal_computed) Util::Abort(INFO, "Operator::Diagonal() must be called before using Fsmooth");
     103              : 
     104              :     // This is a JACOBI iteration, not Gauss-Seidel.
     105              :     // So we need to do twice the number of iterations to get the same behavior as GS.
     106       118500 :     for (int ctr = 0; ctr < 2; ctr++)
     107              :     {
     108        79000 :         Fapply(amrlev, mglev, Ax, x); // find Ax
     109              : 
     110        79000 :         amrex::MultiFab::Copy(Dx, x, 0, 0, ncomp, nghost); // Dx = x
     111        79000 :         amrex::MultiFab::Multiply(Dx, *m_diag[amrlev][mglev], 0, 0, ncomp, nghost); // Dx *= diag  (Dx = x*diag)
     112              : 
     113        79000 :         amrex::MultiFab::Copy(Rx, Ax, 0, 0, ncomp, nghost); // Rx = Ax
     114        79000 :         amrex::MultiFab::Subtract(Rx, Dx, 0, 0, ncomp, nghost); // Rx -= Dx  (Rx = Ax - Dx)
     115              : 
     116       168176 :         for (MFIter mfi(x, false); mfi.isValid(); ++mfi)
     117              :         {
     118        89176 :             const Box& bx = mfi.validbox();
     119        89176 :             amrex::FArrayBox& xfab = x[mfi];
     120        89176 :             const amrex::FArrayBox& bfab = b[mfi];
     121              :             //const amrex::FArrayBox &Axfab   = Ax[mfi];
     122        89176 :             const amrex::FArrayBox& Rxfab = Rx[mfi];
     123        89176 :             const amrex::FArrayBox& diagfab = (*m_diag[amrlev][mglev])[mfi];
     124              : 
     125       305704 :             for (int n = 0; n < ncomp; n++)
     126              :             {
     127    121522704 :                 AMREX_D_TERM(for (int m1 = bx.loVect()[0] - 2; m1 <= bx.hiVect()[0] + 2; m1++),
     128              :                     for (int m2 = bx.loVect()[1] - 2; m2 <= bx.hiVect()[1] + 2; m2++),
     129              :                         for (int m3 = bx.loVect()[2] - 2; m3 <= bx.hiVect()[2] + 2; m3++))
     130              :                 {
     131              : 
     132    109470064 :                     amrex::IntVect m(AMREX_D_DECL(m1, m2, m3));
     133              : 
     134              :                     // Skip ghost cells outside problem domain
     135    249230272 :                     if (AMREX_D_TERM(m[0] < domain.loVect()[0], ||
     136              :                         m[1] < domain.loVect()[1], ||
     137     76887552 :                         m[2] < domain.loVect()[2])) continue;
     138    131859136 :                     if (AMREX_D_TERM(m[0] > domain.hiVect()[0] + 1, ||
     139              :                         m[1] > domain.hiVect()[1] + 1, ||
     140     27417792 :                         m[2] > domain.hiVect()[2] + 1)) continue;
     141              : 
     142    164702272 :                     if (AMREX_D_TERM(m[0] == bx.loVect()[0] - nghost || m[0] == bx.hiVect()[0] + nghost, ||
     143              :                         m[1] == bx.loVect()[1] - nghost || m[1] == bx.hiVect()[1] + nghost, ||
     144              :                         m[2] == bx.loVect()[2] - nghost || m[2] == bx.hiVect()[2] + nghost))
     145              :                     {
     146      2193664 :                         xfab(m, n) = 0.0;
     147      2193664 :                         continue;
     148              :                     }
     149              : 
     150              :                     //xfab(m,n) = xfab(m,n) + omega*(bfab(m,n) - Axfab(m,n))/diagfab(m,n);
     151    195495072 :                     xfab(m, n) = (1. - m_omega) * xfab(m, n) + m_omega * (bfab(m, n) - Rxfab(m, n)) / diagfab(m, n);
     152              :                 }
     153              :             }
     154        79000 :         }
     155              :     }
     156        39500 :     amrex::Geometry geom = m_geom[amrlev][mglev];
     157        39500 :     realFillBoundary(x, geom);
     158        39500 :     nodalSync(amrlev, mglev, x);
     159        39500 : }
     160              : 
     161       285958 : void Operator<Grid::Node>::normalize(int amrlev, int mglev, MultiFab& a_x) const
     162              : {
     163              :     BL_PROFILE("Operator::normalize()");
     164       285958 :     amrex::Box domain(m_geom[amrlev][mglev].Domain());
     165       285958 :     domain.convert(amrex::IntVect::TheNodeVector());
     166              : 
     167       285958 :     int ncomp = getNComp();
     168       285958 :     int nghost = 1; //x.nGrow();
     169              : 
     170       285958 :     if (!m_diagonal_computed)
     171            0 :         Util::Abort(INFO, "Operator::Diagonal() must be called before using normalize");
     172              : 
     173       571916 :     for (MFIter mfi(a_x, amrex::TilingIfNotGPU()); mfi.isValid(); ++mfi)
     174              :     {
     175              : 
     176       285958 :         Box bx = mfi.tilebox();
     177       285958 :         bx.grow(nghost);
     178       285958 :         bx = bx & domain;
     179              : 
     180       285958 :         amrex::Array4<amrex::Real> const& x = a_x.array(mfi);
     181       285958 :         amrex::Array4<const amrex::Real> const& diag = m_diag[amrlev][mglev]->array(mfi);
     182              : 
     183      1011917 :         for (int n = 0; n < ncomp; n++)
     184              :         {
     185       725959 :             amrex::ParallelFor(bx, [=] AMREX_GPU_DEVICE(int i, int j, int k) {
     186              : 
     187     29703906 :                 x(i, j, k) /= diag(i, j, k);
     188              : 
     189     14851953 :             });
     190              :         }
     191       285958 :     }
     192       285958 : }
     193              : 
     194            0 : Operator<Grid::Node>::Operator(const Vector<Geometry>& a_geom,
     195              :     const Vector<BoxArray>& a_grids,
     196              :     const Vector<DistributionMapping>& a_dmap,
     197              :     const LPInfo& a_info,
     198            0 :     const Vector<FabFactory<FArrayBox> const*>& a_factory)
     199              : {
     200              :     BL_PROFILE("Operator::Operator()");
     201            0 :     Util::Message(INFO);
     202              : 
     203            0 :     if (!(a_grids[0].ixType() == amrex::IndexType::TheNodeType()))
     204            0 :         Util::Abort(INFO, "Operator must be defined using CELL CENTERED boxarrays.");
     205              : 
     206            0 :     define(a_geom, a_grids, a_dmap, a_info, a_factory);
     207            0 : }
     208              : 
     209          426 : Operator<Grid::Node>::~Operator()
     210          426 : {}
     211              : 
     212          426 : void Operator<Grid::Node>::define(const Vector<Geometry>& a_geom,
     213              :     const Vector<BoxArray>& a_grids,
     214              :     const Vector<DistributionMapping>& a_dmap,
     215              :     const LPInfo& a_info,
     216              :     const Vector<FabFactory<FArrayBox> const*>& a_factory)
     217              : {
     218              :     BL_PROFILE("Operator::~Operator()");
     219              : 
     220              :     // Make sure we're not trying to parallelize in vain.
     221          426 :     if (amrex::ParallelDescriptor::NProcs() > a_grids[0].size())
     222              :     {
     223            0 :         Util::Warning(INFO, "There are more processors than there are boxes in the amrlev=0 boxarray!!\n",
     224            0 :             "(NProcs = ", amrex::ParallelDescriptor::NProcs(), ", a_grids[0].size() = ", a_grids[0].size(), ")\n",
     225              :             "You should decrease max_grid_size or you will not get proper scaling!");
     226              :     }
     227              : 
     228              :     // This makes sure grids are node-centered;
     229          426 :     Vector<BoxArray> cc_grids = a_grids;
     230          908 :     for (auto& ba : cc_grids) {
     231          482 :         ba.enclosedCells();
     232              :     }
     233              : 
     234          426 :     MLNodeLinOp::define(a_geom, a_grids, a_dmap, a_info, a_factory);
     235              : 
     236          426 :     int nghost = 2;
     237              :     // Resize the multifab containing the operator diagonal
     238          426 :     m_diag.resize(m_num_amr_levels);
     239          898 :     for (int amrlev = 0; amrlev < m_num_amr_levels; ++amrlev)
     240              :     {
     241          472 :         m_diag[amrlev].resize(m_num_mg_levels[amrlev]);
     242              : 
     243         1435 :         for (int mglev = 0; mglev < m_num_mg_levels[amrlev]; ++mglev)
     244              :         {
     245         2889 :             m_diag[amrlev][mglev].reset(new MultiFab(amrex::convert(m_grids[amrlev][mglev], amrex::IntVect::TheNodeVector()),
     246         2889 :                 m_dmap[amrlev][mglev], getNComp(), nghost));
     247              :         }
     248              :     }
     249              : 
     250              :     // We need to instantiate the m_lobc objects.
     251              :     // WE DO NOT USE THEM - our BCs are implemented differently.
     252              :     // But they need to be the right size or the code will segfault.
     253          426 :     m_lobc.resize(getNComp(), { {AMREX_D_DECL(BCType::bogus,BCType::bogus,BCType::bogus)} });
     254          426 :     m_hibc.resize(getNComp(), { {AMREX_D_DECL(BCType::bogus,BCType::bogus,BCType::bogus)} });
     255          426 : }
     256              : 
     257           46 : void Operator<Grid::Node>::fixUpResidualMask(int amrlev, iMultiFab& resmsk)
     258              : {
     259              :     BL_PROFILE("Operator::fixUpResidualMask()");
     260              : 
     261           46 :     if (!m_masks_built) buildMasks();
     262              : 
     263           46 :     const iMultiFab& cfmask = *m_nd_fine_mask[amrlev];
     264              : 
     265              : #ifdef _OPENMP
     266              : #pragma omp parallel if (Gpu::notInLaunchRegion())
     267              : #endif
     268          145 :     for (MFIter mfi(resmsk, TilingIfNotGPU()); mfi.isValid(); ++mfi)
     269              :     {
     270           99 :         const Box& bx = mfi.tilebox();
     271           99 :         Array4<int> const& rmsk = resmsk.array(mfi);
     272           99 :         Array4<int const> const& fmsk = cfmask.const_array(mfi);
     273        57103 :         AMREX_HOST_DEVICE_PARALLEL_FOR_3D(bx, i, j, k,
     274              :             {
     275              :                 if (fmsk(i,j,k) == amrex::nodelap_detail::crse_fine_node) rmsk(i,j,k) = 1;
     276              :             });
     277           46 :     }
     278           46 : }
     279              : 
     280          426 : void Operator<Grid::Node>::prepareForSolve()
     281              : {
     282              :     BL_PROFILE("Operator::prepareForSolve()");
     283          426 :     MLNodeLinOp::prepareForSolve();
     284          426 :     buildMasks();
     285          426 :     averageDownCoeffs();
     286          426 :     Diagonal(true);
     287          426 : }
     288              : 
     289         9253 : void Operator<Grid::Node>::restriction(int amrlev, int cmglev, MultiFab& crse, MultiFab& fine) const
     290              : {
     291              :     BL_PROFILE("Operator::restriction()");
     292              : 
     293         9253 :     applyBC(amrlev, cmglev - 1, fine, BCMode::Homogeneous, StateMode::Solution);
     294              : 
     295         9253 :     amrex::Box cdomain = m_geom[amrlev][cmglev].Domain();
     296         9253 :     cdomain.convert(amrex::IntVect::TheNodeVector());
     297              : 
     298         9253 :     bool need_parallel_copy = !amrex::isMFIterSafe(crse, fine);
     299         9253 :     MultiFab cfine;
     300         9253 :     if (need_parallel_copy) {
     301            0 :         const BoxArray& ba = amrex::coarsen(fine.boxArray(), 2);
     302            0 :         cfine.define(ba, fine.DistributionMap(), fine.nComp(), fine.nGrow());
     303            0 :     }
     304              : 
     305         9253 :     MultiFab* pcrse = (need_parallel_copy) ? &cfine : &crse;
     306              : 
     307        18506 :     for (MFIter mfi(*pcrse, amrex::TilingIfNotGPU()); mfi.isValid(); ++mfi)
     308              :     {
     309         9253 :         const Box& bx = mfi.tilebox();
     310              : 
     311         9253 :         amrex::Array4<const amrex::Real> const& fdata = fine.array(mfi);
     312         9253 :         amrex::Array4<amrex::Real> const& cdata = pcrse->array(mfi);
     313              : 
     314         9253 :         const Dim3 lo = amrex::lbound(cdomain), hi = amrex::ubound(cdomain);
     315              : 
     316              : 
     317        32529 :         for (int n = 0; n < crse.nComp(); n++)
     318              :         {
     319              :             // I,J,K == coarse coordinates
     320              :             // i,j,k == fine coordinates
     321        23276 :             amrex::ParallelFor(bx, [=] AMREX_GPU_DEVICE(int I, int J, int K) {
     322       613704 :                 int i = 2 * I, j = 2 * J, k = 2 * K;
     323              : 
     324       613704 :                 if ((I == lo.x || I == hi.x) &&
     325       328984 :                     (J == lo.y || J == hi.y) &&
     326       207584 :                     (K == lo.z || K == hi.z)) // Corner
     327              :                 {
     328       451032 :                     cdata(I, J, K, n) = fdata(i, j, k, n);
     329              :                 }
     330       463360 :                 else if ((J == lo.y || J == hi.y) &&
     331       182800 :                     (K == lo.z || K == hi.z)) // X edge
     332              :                 {
     333       484700 :                     cdata(I, J, K, n) = 0.25 * fdata(i - 1, j, k, n) + 0.5 * fdata(i, j, k, n) + 0.25 * fdata(i + 1, j, k, n);
     334              :                 }
     335       366420 :                 else if ((K == lo.z || K == hi.z) &&
     336       237630 :                     (I == lo.x || I == hi.x)) // Y edge
     337              :                 {
     338       463900 :                     cdata(I, J, K, n) = 0.25 * fdata(i, j - 1, k, n) + 0.5 * fdata(i, j, k, n) + 0.25 * fdata(i, j + 1, k, n);
     339              :                 }
     340       273640 :                 else if ((I == lo.x || I == hi.x) &&
     341        85860 :                     (J == lo.y || J == hi.y)) // Z edge
     342              :                 {
     343       286200 :                     cdata(I, J, K, n) = 0.25 * fdata(i, j, k - 1, n) + 0.5 * fdata(i, j, k, n) + 0.25 * fdata(i, j, k + 1, n);
     344              :                 }
     345       216400 :                 else if (I == lo.x || I == hi.x) // X face
     346              :                 {
     347        28620 :                     cdata(I, J, K, n) =
     348        85860 :                         (+fdata(i, j - 1, k - 1, n) + 2.0 * fdata(i, j, k - 1, n) + fdata(i, j + 1, k - 1, n)
     349        85860 :                             + 2.0 * fdata(i, j - 1, k, n) + 4.0 * fdata(i, j, k, n) + 2.0 * fdata(i, j + 1, k, n)
     350       114480 :                             + fdata(i, j - 1, k + 1, n) + 2.0 * fdata(i, j, k + 1, n) + fdata(i, j + 1, k + 1, n)) / 16.0;
     351              :                 }
     352       187780 :                 else if (J == lo.y || J == hi.y) // Y face
     353              :                 {
     354        28620 :                     cdata(I, J, K, n) =
     355        85860 :                         (+fdata(i - 1, j, k - 1, n) + 2.0 * fdata(i - 1, j, k, n) + fdata(i - 1, j, k + 1, n)
     356        85860 :                             + 2.0 * fdata(i, j, k - 1, n) + 4.0 * fdata(i, j, k, n) + 2.0 * fdata(i, j, k + 1, n)
     357       114480 :                             + fdata(i + 1, j, k - 1, n) + 2.0 * fdata(i + 1, j, k, n) + fdata(i + 1, j, k + 1, n)) / 16.0;
     358              :                 }
     359       159160 :                 else if (K == lo.z || K == hi.z) // Z face
     360              :                 {
     361       144850 :                     cdata(I, J, K, n) =
     362       434550 :                         (+fdata(i - 1, j - 1, k, n) + 2.0 * fdata(i, j - 1, k, n) + fdata(i + 1, j - 1, k, n)
     363       434550 :                             + 2.0 * fdata(i - 1, j, k, n) + 4.0 * fdata(i, j, k, n) + 2.0 * fdata(i + 1, j, k, n)
     364       579400 :                             + fdata(i - 1, j + 1, k, n) + 2.0 * fdata(i, j + 1, k, n) + fdata(i + 1, j + 1, k, n)) / 16.0;
     365              :                 }
     366              :                 else // Interior
     367        14310 :                     cdata(I, J, K, n) =
     368        57240 :                     (fdata(i - 1, j - 1, k - 1, n) + fdata(i - 1, j - 1, k + 1, n) + fdata(i - 1, j + 1, k - 1, n) + fdata(i - 1, j + 1, k + 1, n) +
     369        57240 :                         fdata(i + 1, j - 1, k - 1, n) + fdata(i + 1, j - 1, k + 1, n) + fdata(i + 1, j + 1, k - 1, n) + fdata(i + 1, j + 1, k + 1, n)) / 64.0
     370        14310 :                     +
     371        57240 :                     (fdata(i, j - 1, k - 1, n) + fdata(i, j - 1, k + 1, n) + fdata(i, j + 1, k - 1, n) + fdata(i, j + 1, k + 1, n) +
     372        57240 :                         fdata(i - 1, j, k - 1, n) + fdata(i + 1, j, k - 1, n) + fdata(i - 1, j, k + 1, n) + fdata(i + 1, j, k + 1, n) +
     373        57240 :                         fdata(i - 1, j - 1, k, n) + fdata(i - 1, j + 1, k, n) + fdata(i + 1, j - 1, k, n) + fdata(i + 1, j + 1, k, n)) / 32.0
     374        14310 :                     +
     375        42930 :                     (fdata(i - 1, j, k, n) + fdata(i, j - 1, k, n) + fdata(i, j, k - 1, n) +
     376        42930 :                         fdata(i + 1, j, k, n) + fdata(i, j + 1, k, n) + fdata(i, j, k + 1, n)) / 16.0
     377        14310 :                     +
     378        28620 :                     fdata(i, j, k, n) / 8.0;
     379       613704 :             });
     380              :         }
     381         9253 :     }
     382              : 
     383         9253 :     if (need_parallel_copy) {
     384            0 :         crse.ParallelCopy(cfine);
     385              :     }
     386              : 
     387         9253 :     amrex::Geometry geom = m_geom[amrlev][cmglev];
     388         9253 :     realFillBoundary(crse, geom);
     389         9253 :     nodalSync(amrlev, cmglev, crse);
     390         9253 : }
     391              : 
     392         9253 : void Operator<Grid::Node>::interpolation(int amrlev, int fmglev, MultiFab& fine, const MultiFab& crse) const
     393              : {
     394              :     BL_PROFILE("Operator::interpolation()");
     395              :     //int nghost = getNGrow();
     396        18506 :     amrex::Box fdomain = m_geom[amrlev][fmglev].Domain(); fdomain.convert(amrex::IntVect::TheNodeVector());
     397              : 
     398         9253 :     bool need_parallel_copy = !amrex::isMFIterSafe(crse, fine);
     399         9253 :     MultiFab cfine;
     400         9253 :     const MultiFab* cmf = &crse;
     401         9253 :     if (need_parallel_copy) {
     402            0 :         const BoxArray& ba = amrex::coarsen(fine.boxArray(), 2);
     403            0 :         cfine.define(ba, fine.DistributionMap(), crse.nComp(), crse.nGrow());
     404            0 :         cfine.ParallelCopy(crse);
     405            0 :         cmf = &cfine;
     406            0 :     }
     407              : 
     408        18506 :     for (MFIter mfi(fine, false); mfi.isValid(); ++mfi)
     409              :     {
     410         9253 :         const Box& fine_bx = mfi.validbox() & fdomain;
     411         9253 :         const Box& course_bx = amrex::coarsen(fine_bx, 2);
     412         9253 :         const Box& tmpbx = amrex::refine(course_bx, 2);
     413         9253 :         FArrayBox tmpfab;
     414         9253 :         tmpfab.resize(tmpbx, fine.nComp());
     415         9253 :         tmpfab.setVal(0.0);
     416         9253 :         const amrex::FArrayBox& crsefab = (*cmf)[mfi];
     417              : 
     418        18506 :         amrex::Array4<const amrex::Real> const& cdata = crsefab.const_array();
     419         9253 :         amrex::Array4<amrex::Real> const& fdata = tmpfab.array();
     420              : 
     421        32529 :         for (int n = 0; n < crse.nComp(); n++)
     422              :         {
     423              :             // I,J,K == coarse coordinates
     424              :             // i,j,k == fine coordinates
     425        23276 :             amrex::ParallelFor(fine_bx, [=] AMREX_GPU_DEVICE(int i, int j, int k) {
     426              : 
     427      2560084 :                 int I = i / 2, J = j / 2, K = k / 2;
     428              : 
     429      2560084 :                 if (i % 2 == 0 && j % 2 == 0 && k % 2 == 0) // Coincident
     430      1841112 :                     fdata(i, j, k, n) = cdata(I, J, K, n);
     431      1946380 :                 else if (j % 2 == 0 && k % 2 == 0) // X Edge
     432      1796848 :                     fdata(i, j, k, n) = 0.5 * (cdata(I, J, K, n) + cdata(I + 1, J, K, n));
     433      1497168 :                 else if (k % 2 == 0 && i % 2 == 0) // Y Edge
     434      1788528 :                     fdata(i, j, k, n) = 0.5 * (cdata(I, J, K, n) + cdata(I, J + 1, K, n));
     435      1050036 :                 else if (i % 2 == 0 && j % 2 == 0) // Z Edge
     436      1030320 :                     fdata(i, j, k, n) = 0.5 * (cdata(I, J, K, n) + cdata(I, J, K + 1, n));
     437       792456 :                 else if (i % 2 == 0) // X Face
     438       515160 :                     fdata(i, j, k, n) = 0.25 * (cdata(I, J, K, n) + cdata(I, J + 1, K, n) +
     439       515160 :                         cdata(I, J, K + 1, n) + cdata(I, J + 1, K + 1, n));
     440       620736 :                 else if (j % 2 == 0) // Y Face
     441       515160 :                     fdata(i, j, k, n) = 0.25 * (cdata(I, J, K, n) + cdata(I, J, K + 1, n) +
     442       515160 :                         cdata(I + 1, J, K, n) + cdata(I + 1, J, K + 1, n));
     443       449016 :                 else if (k % 2 == 0) // Z Face
     444      1003608 :                     fdata(i, j, k, n) = 0.25 * (cdata(I, J, K, n) + cdata(I + 1, J, K, n) +
     445      1003608 :                         cdata(I, J + 1, K, n) + cdata(I + 1, J + 1, K, n));
     446              :                 else // Center
     447       228960 :                     fdata(i, j, k, n) = 0.125 * (cdata(I, J, K, n) +
     448       343440 :                         cdata(I + 1, J, K, n) + cdata(I, J + 1, K, n) + cdata(I, J, K + 1, n) +
     449       343440 :                         cdata(I, J + 1, K + 1, n) + cdata(I + 1, J, K + 1, n) + cdata(I + 1, J + 1, K, n) +
     450       228960 :                         cdata(I + 1, J + 1, K + 1, n));
     451              : 
     452      2560084 :             });
     453              :         }
     454         9253 :         fine[mfi].plus(tmpfab, fine_bx, fine_bx, 0, 0, fine.nComp());
     455        18506 :     }
     456              : 
     457         9253 :     amrex::Geometry geom = m_geom[amrlev][fmglev];
     458         9253 :     realFillBoundary(fine, geom);
     459         9253 :     nodalSync(amrlev, fmglev, fine);
     460         9253 : }
     461              : 
     462           46 : void Operator<Grid::Node>::averageDownSolutionRHS(int camrlev, MultiFab& crse_sol, MultiFab& /*crse_rhs*/,
     463              :     const MultiFab& fine_sol, const MultiFab& /*fine_rhs*/)
     464              : {
     465              :     BL_PROFILE("Operator::averageDownSolutionRHS()");
     466           46 :     const auto& amrrr = AMRRefRatio(camrlev);
     467           46 :     amrex::average_down(fine_sol, crse_sol, 0, crse_sol.nComp(), amrrr);
     468              : 
     469           46 :     if (isSingular(0))
     470              :     {
     471            0 :         Util::Abort(INFO, "Singular operators not supported!");
     472              :     }
     473              : 
     474           46 : }
     475              : 
     476       441392 : void Operator<Grid::Node>::realFillBoundary(MultiFab& phi, const Geometry& geom)
     477              : {
     478              :     Util::RealFillBoundary(phi, geom);
     479       441392 : }
     480              : 
     481       346002 : void Operator<Grid::Node>::applyBC(int amrlev, int mglev, MultiFab& phi, BCMode/* bc_mode*/,
     482              :     amrex::MLLinOp::StateMode /**/, bool skip_fillboundary) const
     483              : {
     484              :     BL_PROFILE("Operator::applyBC()");
     485              : 
     486       346002 :     const Geometry& geom = m_geom[amrlev][mglev];
     487              : 
     488       346002 :     if (!skip_fillboundary) {
     489              : 
     490       346002 :         realFillBoundary(phi, geom);
     491              :     }
     492       346002 : }
     493              : 
     494              : const amrex::FArrayBox&
     495            0 : Operator<Grid::Node>::GetFab(const int num, const int amrlev, const int mglev, const amrex::MFIter& mfi) const
     496              : {
     497              :     BL_PROFILE("Operator::GetFab()");
     498            0 :     Util::Message(INFO);
     499            0 :     return m_a_coeffs[num][amrlev][mglev][mfi];
     500              : }
     501              : 
     502            0 : void Operator<Grid::Node>::RegisterNewFab(amrex::Vector<amrex::MultiFab>& input)
     503              : {
     504              :     BL_PROFILE("Operator::RegisterNewFab()");
     505            0 :     Util::Message(INFO);
     506              :     /// \todo assertions here
     507            0 :     m_a_coeffs.resize(m_a_coeffs.size() + 1);
     508            0 :     m_a_coeffs[m_num_a_fabs].resize(m_num_amr_levels);
     509            0 :     for (int amrlev = 0; amrlev < m_num_amr_levels; ++amrlev)
     510              :     {
     511            0 :         m_a_coeffs[m_num_a_fabs][amrlev].resize(m_num_mg_levels[amrlev]);
     512            0 :         for (int mglev = 0; mglev < m_num_mg_levels[amrlev]; ++mglev)
     513            0 :             m_a_coeffs[m_num_a_fabs][amrlev][mglev].define(m_grids[amrlev][mglev],
     514            0 :                 m_dmap[amrlev][mglev],
     515            0 :                 input[amrlev].nComp(),
     516            0 :                 input[amrlev].nGrow());
     517              : 
     518            0 :         amrex::MultiFab::Copy(m_a_coeffs[m_num_a_fabs][amrlev][0],
     519            0 :             input[amrlev], 0, 0,
     520            0 :             input[amrlev].nComp(),
     521            0 :             input[amrlev].nGrow());
     522              :     }
     523            0 :     m_num_a_fabs++;
     524            0 : }
     525              : 
     526              : 
     527            0 : void Operator<Grid::Node>::RegisterNewFab(amrex::Vector<std::unique_ptr<amrex::MultiFab> >& input)
     528              : {
     529              :     BL_PROFILE("Operator::RegisterNewFab()");
     530            0 :     Util::Message(INFO);
     531              :     /// \todo assertions here
     532            0 :     m_a_coeffs.resize(m_a_coeffs.size() + 1);
     533            0 :     m_a_coeffs[m_num_a_fabs].resize(m_num_amr_levels);
     534            0 :     for (int amrlev = 0; amrlev < m_num_amr_levels; ++amrlev)
     535              :     {
     536            0 :         m_a_coeffs[m_num_a_fabs][amrlev].resize(m_num_mg_levels[amrlev]);
     537            0 :         for (int mglev = 0; mglev < m_num_mg_levels[amrlev]; ++mglev)
     538            0 :             m_a_coeffs[m_num_a_fabs][amrlev][mglev].define(m_grids[amrlev][mglev],
     539            0 :                 m_dmap[amrlev][mglev],
     540            0 :                 input[amrlev]->nComp(),
     541            0 :                 input[amrlev]->nGrow());
     542              : 
     543            0 :         amrex::MultiFab::Copy(m_a_coeffs[m_num_a_fabs][amrlev][0],
     544            0 :             *input[amrlev], 0, 0,
     545            0 :             input[amrlev]->nComp(),
     546            0 :             input[amrlev]->nGrow());
     547              :     }
     548            0 :     m_num_a_fabs++;
     549            0 : }
     550              : 
     551          354 : void Operator<Grid::Node>::reflux(int crse_amrlev,
     552              :     MultiFab& res, const MultiFab& /*crse_sol*/, const MultiFab& /*crse_rhs*/,
     553              :     MultiFab& fine_res, MultiFab& /*fine_sol*/, const MultiFab& /*fine_rhs*/) const
     554              : {
     555              :     BL_PROFILE("Operator::Elastic::reflux()");
     556              : 
     557          354 :     int ncomp = AMREX_SPACEDIM;
     558              : 
     559          354 :     amrex::Box cdomain(m_geom[crse_amrlev][0].Domain());
     560          354 :     cdomain.convert(amrex::IntVect::TheNodeVector());
     561              : 
     562          354 :     const Geometry& cgeom = m_geom[crse_amrlev][0];
     563              : 
     564          354 :     const BoxArray& fba = fine_res.boxArray();
     565          354 :     const DistributionMapping& fdm = fine_res.DistributionMap();
     566              : 
     567          354 :     MultiFab fine_res_for_coarse(amrex::coarsen(fba, 2), fdm, ncomp, 2);
     568          354 :     fine_res_for_coarse.ParallelCopy(res, 0, 0, ncomp, 0, 0, cgeom.periodicity());
     569              : 
     570          354 :     applyBC(crse_amrlev + 1, 0, fine_res, BCMode::Inhomogeneous, StateMode::Solution);
     571              : 
     572              :     /// \todo Replace with Enum
     573              :     // const int coarse_coarse_node = 0;
     574          354 :     const int coarse_fine_node = 1;
     575          354 :     const int fine_fine_node = 2;
     576              : 
     577          354 :     amrex::iMultiFab nodemask(amrex::coarsen(fba, 2), fdm, 1, 2);
     578          354 :     nodemask.ParallelCopy(*m_nd_fine_mask[crse_amrlev], 0, 0, 1, 0, 0, cgeom.periodicity());
     579              : 
     580          708 :     amrex::iMultiFab cellmask(amrex::convert(amrex::coarsen(fba, 2), amrex::IntVect::TheCellVector()), fdm, 1, 2);
     581          354 :     cellmask.ParallelCopy(*m_cc_fine_mask[crse_amrlev], 0, 0, 1, 1, 1, cgeom.periodicity());
     582              : 
     583         2276 :     for (MFIter mfi(fine_res_for_coarse, false); mfi.isValid(); ++mfi)
     584              :     {
     585         1922 :         const Box& bx = mfi.validbox();
     586              : 
     587         1922 :         amrex::Array4<const int> const& nmask = nodemask.array(mfi);
     588              :         //amrex::Array4<const int> const& cmask = cellmask.array(mfi);
     589              : 
     590         1922 :         amrex::Array4<amrex::Real> const& cdata = fine_res_for_coarse.array(mfi);
     591         1922 :         amrex::Array4<const amrex::Real> const& fdata = fine_res.array(mfi);
     592              : 
     593         1922 :         const Dim3 lo = amrex::lbound(cdomain), hi = amrex::ubound(cdomain);
     594              : 
     595         5766 :         for (int n = 0; n < fine_res.nComp(); n++)
     596              :         {
     597              :             // I,J,K == coarse coordinates
     598              :             // i,j,k == fine coordinates
     599         3844 :             amrex::ParallelFor(bx, [=] AMREX_GPU_DEVICE(int I, int J, int K) {
     600       445594 :                 int i = I * 2, j = J * 2, k = K * 2;
     601              : 
     602       971292 :                 if (nmask(I, J, K) == fine_fine_node || nmask(I, J, K) == coarse_fine_node)
     603              :                 {
     604       445594 :                     if ((I == lo.x || I == hi.x) &&
     605            0 :                         (J == lo.y || J == hi.y) &&
     606            0 :                         (K == lo.z || K == hi.z)) // Corner
     607            0 :                         cdata(I, J, K, n) = fdata(i, j, k, n);
     608       445594 :                     else if ((J == lo.y || J == hi.y) &&
     609            0 :                         (K == lo.z || K == hi.z)) // X edge
     610            0 :                         cdata(I, J, K, n) = 0.25 * fdata(i - 1, j, k, n) + 0.5 * fdata(i, j, k, n) + 0.25 * fdata(i + 1, j, k, n);
     611       445594 :                     else if ((K == lo.z || K == hi.z) &&
     612       445594 :                         (I == lo.x || I == hi.x)) // Y edge
     613            0 :                         cdata(I, J, K, n) = 0.25 * fdata(i, j - 1, k, n) + 0.5 * fdata(i, j, k, n) + 0.25 * fdata(i, j + 1, k, n);
     614       445594 :                     else if ((I == lo.x || I == hi.x) &&
     615            0 :                         (J == lo.y || J == hi.y)) // Z edge
     616            0 :                         cdata(I, J, K, n) = 0.25 * fdata(i, j, k - 1, n) + 0.5 * fdata(i, j, k, n) + 0.25 * fdata(i, j, k + 1, n);
     617       445594 :                     else if (I == lo.x || I == hi.x) // X face
     618            0 :                         cdata(I, J, K, n) =
     619            0 :                         (+fdata(i, j - 1, k - 1, n) + 2.0 * fdata(i, j, k - 1, n) + fdata(i, j + 1, k - 1, n)
     620            0 :                             + 2.0 * fdata(i, j - 1, k, n) + 4.0 * fdata(i, j, k, n) + 2.0 * fdata(i, j + 1, k, n)
     621            0 :                             + fdata(i, j - 1, k + 1, n) + 2.0 * fdata(i, j, k + 1, n) + fdata(i, j + 1, k + 1, n)) / 16.0;
     622       445594 :                     else if (J == lo.y || J == hi.y) // Y face
     623            0 :                         cdata(I, J, K, n) =
     624            0 :                         (+fdata(i - 1, j, k - 1, n) + 2.0 * fdata(i - 1, j, k, n) + fdata(i - 1, j, k + 1, n)
     625            0 :                             + 2.0 * fdata(i, j, k - 1, n) + 4.0 * fdata(i, j, k, n) + 2.0 * fdata(i, j, k + 1, n)
     626            0 :                             + fdata(i + 1, j, k - 1, n) + 2.0 * fdata(i + 1, j, k, n) + fdata(i + 1, j, k + 1, n)) / 16.0;
     627       445594 :                     else if (K == lo.z || K == hi.z) // Z face
     628       445594 :                         cdata(I, J, K, n) =
     629      1336782 :                         (+fdata(i - 1, j - 1, k, n) + 2.0 * fdata(i, j - 1, k, n) + fdata(i + 1, j - 1, k, n)
     630      1336782 :                             + 2.0 * fdata(i - 1, j, k, n) + 4.0 * fdata(i, j, k, n) + 2.0 * fdata(i + 1, j, k, n)
     631      1782376 :                             + fdata(i - 1, j + 1, k, n) + 2.0 * fdata(i, j + 1, k, n) + fdata(i + 1, j + 1, k, n)) / 16.0;
     632              :                     else // Interior
     633            0 :                         cdata(I, J, K, n) =
     634            0 :                         (fdata(i - 1, j - 1, k - 1, n) + fdata(i - 1, j - 1, k + 1, n) + fdata(i - 1, j + 1, k - 1, n) + fdata(i - 1, j + 1, k + 1, n) +
     635            0 :                             fdata(i + 1, j - 1, k - 1, n) + fdata(i + 1, j - 1, k + 1, n) + fdata(i + 1, j + 1, k - 1, n) + fdata(i + 1, j + 1, k + 1, n)) / 64.0
     636            0 :                         +
     637            0 :                         (fdata(i, j - 1, k - 1, n) + fdata(i, j - 1, k + 1, n) + fdata(i, j + 1, k - 1, n) + fdata(i, j + 1, k + 1, n) +
     638            0 :                             fdata(i - 1, j, k - 1, n) + fdata(i + 1, j, k - 1, n) + fdata(i - 1, j, k + 1, n) + fdata(i + 1, j, k + 1, n) +
     639            0 :                             fdata(i - 1, j - 1, k, n) + fdata(i - 1, j + 1, k, n) + fdata(i + 1, j - 1, k, n) + fdata(i + 1, j + 1, k, n)) / 32.0
     640            0 :                         +
     641            0 :                         (fdata(i - 1, j, k, n) + fdata(i, j - 1, k, n) + fdata(i, j, k - 1, n) +
     642            0 :                             fdata(i + 1, j, k, n) + fdata(i, j + 1, k, n) + fdata(i, j, k + 1, n)) / 16.0
     643            0 :                         +
     644            0 :                         fdata(i, j, k, n) / 8.0;
     645              :                 }
     646              : 
     647       445594 :             });
     648              :         }
     649          354 :     }
     650              : 
     651              :     // Copy the fine residual restricted onto the coarse grid
     652              :     // into the final residual.
     653          354 :     res.ParallelCopy(fine_res_for_coarse, 0, 0, ncomp, 0, 0, cgeom.periodicity());
     654              : 
     655          354 :     const int mglev = 0;
     656              : 
     657              :     // Sync up ghost nodes
     658          354 :     amrex::Geometry geom = m_geom[crse_amrlev][mglev];
     659          354 :     realFillBoundary(res, geom);
     660          354 :     nodalSync(crse_amrlev, mglev, res);
     661          708 :     return;
     662          354 : }
     663              : 
     664              : void
     665         9618 : Operator<Grid::Node>::solutionResidual(int amrlev, MultiFab& resid, MultiFab& x, const MultiFab& b,
     666              :     const MultiFab* /*crse_bcdata*/)
     667              : {
     668         9618 :     const int mglev = 0;
     669         9618 :     const int ncomp = b.nComp();
     670         9618 :     apply(amrlev, mglev, resid, x, BCMode::Inhomogeneous, StateMode::Solution);
     671         9618 :     MultiFab::Xpay(resid, -1.0, b, 0, 0, ncomp, 2);
     672         9618 :     amrex::Geometry geom = m_geom[amrlev][mglev];
     673         9618 :     realFillBoundary(resid, geom);
     674         9618 : }
     675              : 
     676              : void
     677        27412 : Operator<Grid::Node>::correctionResidual(int amrlev, int mglev, MultiFab& resid, MultiFab& x, const MultiFab& b,
     678              :     BCMode /*bc_mode*/, const MultiFab* /*crse_bcdata*/)
     679              : {
     680        27412 :     resid.setVal(0.0);
     681        27412 :     apply(amrlev, mglev, resid, x, BCMode::Homogeneous, StateMode::Correction);
     682        27412 :     int ncomp = b.nComp();
     683        27412 :     MultiFab::Xpay(resid, -1.0, b, 0, 0, ncomp, resid.nGrow());
     684        27412 :     amrex::Geometry geom = m_geom[amrlev][mglev];
     685        27412 :     realFillBoundary(resid, geom);
     686        27412 : }
     687              : 
     688              : 
     689              : 
     690              : 
     691            0 : Operator<Grid::Cell>::Operator()
     692              : {
     693            0 :     m_ixtype = amrex::IntVect::TheCellVector();
     694            0 : }
     695              : 
     696              : void
     697            0 : Operator<Grid::Cell>::define(amrex::Vector<amrex::Geometry> a_geom,
     698              :     const amrex::Vector<amrex::BoxArray>& a_grids,
     699              :     const amrex::Vector<amrex::DistributionMapping>& a_dmap,
     700              :     BC::BC<Set::Scalar>& a_bc,
     701              :     const amrex::LPInfo& a_info,
     702              :     const amrex::Vector<amrex::FabFactory<amrex::FArrayBox> const*>& a_factory)
     703              : {
     704            0 :     m_bc = &a_bc;
     705              : 
     706            0 :     std::array<int, AMREX_SPACEDIM> is_periodic = m_bc->IsPeriodic();
     707              : 
     708            0 :     MLCellLinOp::define(a_geom, a_grids, a_dmap, a_info, a_factory);
     709              : 
     710            0 :     Util::Warning(INFO, "This section of code has not been tested.");
     711            0 :     for (int n = 0; n < getNComp(); n++)
     712              :     {
     713            0 :         m_lobc.push_back({ AMREX_D_DECL(is_periodic[0] ? amrex::LinOpBCType::Periodic : amrex::LinOpBCType::Dirichlet,
     714              :                         is_periodic[1] ? amrex::LinOpBCType::Periodic : amrex::LinOpBCType::Dirichlet,
     715              :                         is_periodic[2] ? amrex::LinOpBCType::Periodic : amrex::LinOpBCType::Dirichlet) });
     716            0 :         m_hibc.push_back({ AMREX_D_DECL(is_periodic[0] ? amrex::LinOpBCType::Periodic : amrex::LinOpBCType::Dirichlet,
     717              :                         is_periodic[1] ? amrex::LinOpBCType::Periodic : amrex::LinOpBCType::Dirichlet,
     718              :                         is_periodic[2] ? amrex::LinOpBCType::Periodic : amrex::LinOpBCType::Dirichlet) });
     719              :     }
     720              : 
     721            0 :     for (int ilev = 0; ilev < a_geom.size(); ++ilev)
     722            0 :         setLevelBC(ilev, nullptr);
     723              : 
     724            0 : }
     725              : 
     726              : 
     727              : void
     728            0 : Operator<Grid::Cell>::prepareForSolve()
     729              : {
     730            0 :     MLCellLinOp::prepareForSolve();
     731            0 : }
     732              : 
     733            0 : Operator<Grid::Cell>::BndryCondLoc::BndryCondLoc(const amrex::BoxArray& ba, const amrex::DistributionMapping& dm)
     734            0 :     : bcond(ba, dm),
     735            0 :     bcloc(ba, dm)
     736              : {
     737            0 : }
     738              : 
     739              : void
     740            0 : Operator<Grid::Cell>::BndryCondLoc::setLOBndryConds(const amrex::Geometry& /*geom*/, const amrex::Real* /*dx*/,
     741              :     const amrex::Array<BCType, AMREX_SPACEDIM>& /*lobc*/,
     742              :     const amrex::Array<BCType, AMREX_SPACEDIM>& /*hibc*/,
     743              :     int /*ratio*/, const amrex::RealVect& /*a_loc*/)
     744              : {
     745            0 :     Util::Warning(INFO, "This code has not been properlyt tested");
     746            0 : }
     747              : 
     748              : 
     749              : void
     750            0 : Operator<Grid::Cell>::averageDownCoeffs()
     751              : {
     752            0 :     for (int i = 0; i < m_num_a_fabs; i++)
     753              :     {
     754            0 :         for (int amrlev = m_num_amr_levels - 1; amrlev > 0; --amrlev)
     755              :         {
     756            0 :             auto& fine_a_coeffs = m_a_coeffs[i][amrlev];
     757            0 :             averageDownCoeffsSameAmrLevel(fine_a_coeffs);
     758              :         }
     759            0 :         averageDownCoeffsSameAmrLevel(m_a_coeffs[i][0]);
     760              :     }
     761            0 : }
     762              : 
     763              : void
     764            0 : Operator<Grid::Cell>::averageDownCoeffsSameAmrLevel(amrex::Vector<amrex::MultiFab>& a)
     765              : {
     766            0 :     int nmglevs = a.size();
     767            0 :     for (int mglev = 1; mglev < nmglevs; ++mglev)
     768              :     {
     769            0 :         amrex::average_down(a[mglev - 1], a[mglev], 0, a[0].nComp(), mg_coarsen_ratio);
     770              :     }
     771            0 : }
     772              : 
     773              : 
     774              : 
     775              : const amrex::FArrayBox&
     776            0 : Operator<Grid::Cell>::GetFab(const int num, const int amrlev, const int mglev, const amrex::MFIter& mfi) const
     777              : {
     778            0 :     return m_a_coeffs[num][amrlev][mglev][mfi];
     779              : }
     780              : 
     781              : 
     782              : void
     783            0 : Operator<Grid::Cell>::RegisterNewFab(amrex::Vector<amrex::MultiFab>& input)
     784              : {
     785              :     /// \todo assertions here
     786            0 :     m_a_coeffs.resize(m_a_coeffs.size() + 1);
     787            0 :     m_a_coeffs[m_num_a_fabs].resize(m_num_amr_levels);
     788            0 :     for (int amrlev = 0; amrlev < m_num_amr_levels; ++amrlev)
     789              :     {
     790            0 :         m_a_coeffs[m_num_a_fabs][amrlev].resize(m_num_mg_levels[amrlev]);
     791            0 :         for (int mglev = 0; mglev < m_num_mg_levels[amrlev]; ++mglev)
     792            0 :             m_a_coeffs[m_num_a_fabs][amrlev][mglev].define(m_grids[amrlev][mglev],
     793            0 :                 m_dmap[amrlev][mglev],
     794            0 :                 input[amrlev].nComp(),
     795            0 :                 input[amrlev].nGrow());
     796              : 
     797            0 :         amrex::MultiFab::Copy(m_a_coeffs[m_num_a_fabs][amrlev][0],
     798            0 :             input[amrlev], 0, 0,
     799            0 :             input[amrlev].nComp(),
     800            0 :             input[amrlev].nGrow());
     801              :     }
     802            0 :     m_num_a_fabs++;
     803            0 : }
     804              : 
     805              : 
     806              : void
     807            0 : Operator<Grid::Cell>::RegisterNewFab(amrex::Vector<std::unique_ptr<amrex::MultiFab> >& input)
     808              : {
     809              :     /// \todo assertions here
     810            0 :     m_a_coeffs.resize(m_a_coeffs.size() + 1);
     811            0 :     m_a_coeffs[m_num_a_fabs].resize(m_num_amr_levels);
     812            0 :     for (int amrlev = 0; amrlev < m_num_amr_levels; ++amrlev)
     813              :     {
     814            0 :         m_a_coeffs[m_num_a_fabs][amrlev].resize(m_num_mg_levels[amrlev]);
     815            0 :         for (int mglev = 0; mglev < m_num_mg_levels[amrlev]; ++mglev)
     816            0 :             m_a_coeffs[m_num_a_fabs][amrlev][mglev].define(m_grids[amrlev][mglev],
     817            0 :                 m_dmap[amrlev][mglev],
     818            0 :                 input[amrlev]->nComp(),
     819            0 :                 input[amrlev]->nGrow());
     820              : 
     821            0 :         amrex::MultiFab::Copy(m_a_coeffs[m_num_a_fabs][amrlev][0],
     822            0 :             *input[amrlev], 0, 0,
     823            0 :             input[amrlev]->nComp(),
     824            0 :             input[amrlev]->nGrow());
     825              :     }
     826            0 :     m_num_a_fabs++;
     827            0 : }
     828              : 
     829              : 
     830              : }
        

Generated by: LCOV version 2.0-1