Operator.cpp

Go to the documentation of this file.

 
#include <AMReX_MLNodeLinOp.H>
#include <AMReX_MLCellLinOp.H>
#include <AMReX_MLNodeLap_K.H>
#include <AMReX_MultiFabUtil.H>
#include "Util/Color.H"
#include "Set/Set.H"
#include "Operator.H"
 
using namespace amrex;
namespace Operator {
 
// constexpr amrex::IntVect AMREX_D_DECL(Operator<Grid::Node>::dx,Operator<Grid::Node>::dy,Operator<Grid::Node>::dz);
constexpr amrex::IntVect AMREX_D_DECL(Operator<Grid::Cell>::dx, Operator<Grid::Cell>::dy, Operator<Grid::Cell>::dz);
 
void Operator<Grid::Node>::Diagonal(bool recompute)
{
    BL_PROFILE(Color::FG::Yellow + "Operator::Diagonal()" + Color::Reset);
    if (!recompute && m_diagonal_computed) return;
    m_diagonal_computed = true;
 
    for (int amrlev = 0; amrlev < m_num_amr_levels; ++amrlev)
    {
        for (int mglev = 0; mglev < m_num_mg_levels[amrlev]; ++mglev)
        {
            Diagonal(amrlev, mglev, *m_diag[amrlev][mglev]);
        }
    }
}
 
void Operator<Grid::Node>::Diagonal(int amrlev, int mglev, amrex::MultiFab& diag)
{
    BL_PROFILE("Operator::Diagonal()");
    //Util::Message(INFO);
 
    int ncomp = diag.nComp();
    int nghost = 0;
 
    int sep = 2;
    int num = AMREX_D_TERM(sep, *sep, *sep);
    int cntr = 0;
 
    amrex::MultiFab x(m_diag[amrlev][mglev]->boxArray(), m_diag[amrlev][mglev]->DistributionMap(), ncomp, nghost);
    amrex::MultiFab Ax(m_diag[amrlev][mglev]->boxArray(), m_diag[amrlev][mglev]->DistributionMap(), ncomp, nghost);
 
    for (MFIter mfi(x, false); mfi.isValid(); ++mfi)
    {
        const Box& bx = mfi.validbox();
        amrex::FArrayBox& diagfab = diag[mfi];
        amrex::FArrayBox& xfab = x[mfi];
        amrex::FArrayBox& Axfab = Ax[mfi];
 
        diagfab.setVal<amrex::RunOn::Device>(0.0);
 
        for (int i = 0; i < num; i++)
        {
            for (int n = 0; n < ncomp; n++)
            {
                xfab.setVal<amrex::RunOn::Device>(0.0);
                Axfab.setVal<amrex::RunOn::Device>(0.0);
 
                //BL_PROFILE_VAR("Operator::Part1", part1); 
                AMREX_D_TERM(for (int m1 = bx.loVect()[0]; m1 <= bx.hiVect()[0]; m1++),
                    for (int m2 = bx.loVect()[1]; m2 <= bx.hiVect()[1]; m2++),
                        for (int m3 = bx.loVect()[2]; m3 <= bx.hiVect()[2]; m3++))
                {
                    amrex::IntVect m(AMREX_D_DECL(m1, m2, m3));
 
                    if (m1 % sep == i / sep && m2 % sep == i % sep) xfab(m, n) = 1.0;
                    else xfab(m, n) = 0.0;
                }
                //BL_PROFILE_VAR_STOP(part1);
 
                BL_PROFILE_VAR("Operator::Part2", part2);
                Util::Message(INFO, "Calling fapply...", cntr++);
                Fapply(amrlev, mglev, Ax, x);
                BL_PROFILE_VAR_STOP(part2);
 
                //BL_PROFILE_VAR("Operator::Part3", part3); 
                Axfab.mult<amrex::RunOn::Device>(xfab, n, n, 1);
                diagfab.plus<amrex::RunOn::Device>(Axfab, n, n, 1);
                //BL_PROFILE_VAR_STOP(part3);
            }
        }
    }
}
 
void Operator<Grid::Node>::Fsmooth(int amrlev, int mglev, amrex::MultiFab& x, const amrex::MultiFab& b) const
{
    BL_PROFILE("Operator::Fsmooth()");
 
    amrex::Box domain(m_geom[amrlev][mglev].growPeriodicDomain(1));
    domain.convert(amrex::IntVect::TheNodeVector());
 
    int ncomp = b.nComp();
    const bool relax_ghost_rows = relaxCoarseFineGhostRows();
    int nghost = relax_ghost_rows ? 2 : 0;
 
 
    amrex::MultiFab Ax(x.boxArray(), x.DistributionMap(), ncomp, nghost);
    amrex::MultiFab Dx(x.boxArray(), x.DistributionMap(), ncomp, nghost);
    amrex::MultiFab Rx(x.boxArray(), x.DistributionMap(), ncomp, nghost);
 
    if (!m_diagonal_computed) Util::Abort(INFO, "Operator::Diagonal() must be called before using Fsmooth");
 
    // This is a JACOBI iteration, not Gauss-Seidel.
    // So we need to do twice the number of iterations to get the same behavior as GS.
    for (int ctr = 0; ctr < 2; ctr++)
    {
        Fapply(amrlev, mglev, Ax, x); // find Ax
 
        amrex::MultiFab::Copy(Dx, x, 0, 0, ncomp, nghost); // Dx = x
        amrex::MultiFab::Multiply(Dx, *m_diag[amrlev][mglev], 0, 0, ncomp, nghost); // Dx *= diag  (Dx = x*diag)
 
        amrex::MultiFab::Copy(Rx, Ax, 0, 0, ncomp, nghost); // Rx = Ax
        amrex::MultiFab::Subtract(Rx, Dx, 0, 0, ncomp, nghost); // Rx -= Dx  (Rx = Ax - Dx)
 
        for (MFIter mfi(x, false); mfi.isValid(); ++mfi)
        {
            Box bx = mfi.grownnodaltilebox();
            
            const auto xfab = x.array(mfi);
            const auto bfab = b.const_array(mfi);
            const auto Rxfab = Rx.const_array(mfi);
            const auto diagfab = (*m_diag[amrlev][mglev]).const_array(mfi);
 
            if (!relax_ghost_rows)
            {
                const Box cbx = mfi.nodaltilebox() & domain;
                for (int n = 0; n < ncomp; ++n)
                {
                    amrex::LoopConcurrentOnCpu(cbx, [&] (int i, int j, int k)
                    {
                        xfab(i,j,k,n) = (1. - m_omega) * xfab(i,j,k, n)
                            + m_omega * (bfab(i,j,k, n) - Rxfab(i,j,k, n))
                            / diagfab(i,j,k,n);
                    });
                }
                continue;
            }
 
 
            for (int n = 0; n < ncomp; n++)
            {
                const Set::Scalar omega = m_omega;
                amrex::ParallelFor(bx, [=] AMREX_GPU_DEVICE(int i, int j, int k)
                {
 
                    // Skip ghost cells outside problem domain
                    if (!domain.contains(i,j,k))
                    {
                        //continue;
                    }
                    else if ( !bx.strictly_contains(i,j,k))
                    {
                        xfab(i, j, k, n) = 0.0;
                        //continue;
                    }
                    else
                    {
                        xfab(i,j,k,n) = (1. - omega) * xfab(i,j,k, n) + omega * (bfab(i,j,k, n) - Rxfab(i,j,k, n)) / diagfab(i,j,k,n);
                    }
                });
            }
        }
        amrex::Geometry geom = m_geom[amrlev][mglev];
        x.setMultiGhost(true);
        x.FillBoundary(geom.periodicity());
        nodalSync(amrlev, mglev, x);
    }
}
 
void Operator<Grid::Node>::normalize(int amrlev, int mglev, MultiFab& a_x) const
{
    if (!m_diagonal_computed)
        Util::Abort(INFO, "Operator::Diagonal() must be called before using normalize");
 
    a_x.divide(*m_diag[amrlev][mglev],0,getNComp(),relaxCoarseFineGhostRows() ? 2 : 0);
 
    a_x.setMultiGhost(true);
    a_x.FillBoundaryAndSync(Geom(amrlev,mglev).periodicity());
}
 
Operator<Grid::Node>::Operator(const Vector<Geometry>& a_geom,
    const Vector<BoxArray>& a_grids,
    const Vector<DistributionMapping>& a_dmap,
    const LPInfo& a_info,
    const Vector<FabFactory<FArrayBox> const*>& a_factory)
{
    BL_PROFILE("Operator::Operator()");
    Util::Message(INFO);
 
    if (!(a_grids[0].ixType() == amrex::IndexType::TheNodeType()))
        Util::Abort(INFO, "Operator must be defined using CELL CENTERED boxarrays.");
 
    define(a_geom, a_grids, a_dmap, a_info, a_factory);
}
 
Operator<Grid::Node>::~Operator()
{}
 
void Operator<Grid::Node>::define(const Vector<Geometry>& a_geom,
    const Vector<BoxArray>& a_grids,
    const Vector<DistributionMapping>& a_dmap,
    const LPInfo& a_info,
    const Vector<FabFactory<FArrayBox> const*>& a_factory)
{
    BL_PROFILE("Operator::~Operator()");
 
    // Make sure we're not trying to parallelize in vain.
    if (amrex::ParallelDescriptor::NProcs() > a_grids[0].size())
    {
        Util::Warning(INFO, "There are more processors than there are boxes in the amrlev=0 boxarray!!\n",
            "(NProcs = ", amrex::ParallelDescriptor::NProcs(), ", a_grids[0].size() = ", a_grids[0].size(), ")\n",
            "You should decrease max_grid_size or you will not get proper scaling!");
    }
 
    // This makes sure grids are node-centered;
    Vector<BoxArray> cc_grids = a_grids;
    for (auto& ba : cc_grids) {
        ba.enclosedCells();
    }
 
    MLNodeLinOp::define(a_geom, a_grids, a_dmap, a_info, a_factory);
 
    int nghost = 2;
    // Resize the multifab containing the operator diagonal
    m_diag.resize(m_num_amr_levels);
    for (int amrlev = 0; amrlev < m_num_amr_levels; ++amrlev)
    {
        m_diag[amrlev].resize(m_num_mg_levels[amrlev]);
 
        for (int mglev = 0; mglev < m_num_mg_levels[amrlev]; ++mglev)
        {
            m_diag[amrlev][mglev].reset(new MultiFab(amrex::convert(m_grids[amrlev][mglev], amrex::IntVect::TheNodeVector()),
                m_dmap[amrlev][mglev], getNComp(), nghost));
        }
    }
 
    // We need to instantiate the m_lobc objects.
    // WE DO NOT USE THEM - our BCs are implemented differently.
    // But they need to be the right size or the code will segfault.
    m_lobc.resize(getNComp(), { {AMREX_D_DECL(BCType::bogus,BCType::bogus,BCType::bogus)} });
    m_hibc.resize(getNComp(), { {AMREX_D_DECL(BCType::bogus,BCType::bogus,BCType::bogus)} });
}
 
void Operator<Grid::Node>::fixUpResidualMask(int amrlev, iMultiFab& resmsk)
{
    BL_PROFILE("Operator::fixUpResidualMask()");
 
    if (!m_masks_built) buildMasks();
 
    const iMultiFab& cfmask = *m_nd_fine_mask[amrlev];
 
#ifdef _OPENMP
#pragma omp parallel if (Gpu::notInLaunchRegion())
#endif
    for (MFIter mfi(resmsk, TilingIfNotGPU()); mfi.isValid(); ++mfi)
    {
        const Box& bx = mfi.tilebox();
        Array4<int> const& rmsk = resmsk.array(mfi);
        Array4<int const> const& fmsk = cfmask.const_array(mfi);
        AMREX_HOST_DEVICE_PARALLEL_FOR_3D(bx, i, j, k,
            {
                if (fmsk(i,j,k) == amrex::nodelap_detail::crse_fine_node) rmsk(i,j,k) = 1;
            });
    }
}
 
void Operator<Grid::Node>::prepareForSolve()
{
    BL_PROFILE("Operator::prepareForSolve()");
    MLNodeLinOp::prepareForSolve();
    buildMasks();
    averageDownCoeffs();
    Diagonal(true);
}
 
void Operator<Grid::Node>::restriction(int amrlev, int cmglev, MultiFab& crse, MultiFab& fine) const
{
    BL_PROFILE("Operator::restriction()");
 
    applyBC(amrlev, cmglev - 1, fine, BCMode::Homogeneous, StateMode::Solution);
 
    amrex::Box cdomain = m_geom[amrlev][cmglev].growPeriodicDomain(1);
    cdomain = cdomain.convert(amrex::IntVect::TheNodeVector());
 
    bool need_parallel_copy = !amrex::isMFIterSafe(crse, fine);
    MultiFab cfine;
    if (need_parallel_copy) {
        const BoxArray& ba = amrex::coarsen(fine.boxArray(), 2);
        cfine.define(ba, fine.DistributionMap(), fine.nComp(), fine.nGrow());
    }
 
    MultiFab* pcrse = (need_parallel_copy) ? &cfine : &crse;
 
    for (MFIter mfi(*pcrse, amrex::TilingIfNotGPU()); mfi.isValid(); ++mfi)
    {
        const Box& bx = mfi.grownnodaltilebox(-1,1) & cdomain;
 
        amrex::Array4<const amrex::Real> const& fdata = fine.array(mfi);
        amrex::Array4<amrex::Real> const& cdata = pcrse->array(mfi);
 
        const Dim3 lo = amrex::lbound(bx), hi = amrex::ubound(bx);
 
 
        for (int n = 0; n < crse.nComp(); n++)
        {
            // I,J,K == coarse coordinates
            // i,j,k == fine coordinates
            amrex::ParallelFor(bx, [=] AMREX_GPU_DEVICE(int I, int J, int K) {
                int i = 2 * I, j = 2 * J, k = 2 * K;
 
                if ((I == lo.x || I == hi.x) &&
                    (J == lo.y || J == hi.y) &&
                    (K == lo.z || K == hi.z)) // Corner
                {
                    cdata(I, J, K, n) = fdata(i, j, k, n);
                }
                else if ((J == lo.y || J == hi.y) &&
                    (K == lo.z || K == hi.z)) // X edge
                {
                    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);
                }
                else if ((K == lo.z || K == hi.z) &&
                    (I == lo.x || I == hi.x)) // Y edge
                {
                    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);
                }
                else if ((I == lo.x || I == hi.x) &&
                    (J == lo.y || J == hi.y)) // Z edge
                {
                    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);
                }
                else if (I == lo.x || I == hi.x) // X face
                {
                    cdata(I, J, K, n) =
                        (+fdata(i, j - 1, k - 1, n) + 2.0 * fdata(i, j, k - 1, n) + fdata(i, j + 1, k - 1, n)
                            + 2.0 * fdata(i, j - 1, k, n) + 4.0 * fdata(i, j, k, n) + 2.0 * fdata(i, j + 1, k, n)
                            + fdata(i, j - 1, k + 1, n) + 2.0 * fdata(i, j, k + 1, n) + fdata(i, j + 1, k + 1, n)) / 16.0;
                }
                else if (J == lo.y || J == hi.y) // Y face
                {
                    cdata(I, J, K, n) =
                        (+fdata(i - 1, j, k - 1, n) + 2.0 * fdata(i - 1, j, k, n) + fdata(i - 1, j, k + 1, n)
                            + 2.0 * fdata(i, j, k - 1, n) + 4.0 * fdata(i, j, k, n) + 2.0 * fdata(i, j, k + 1, n)
                            + fdata(i + 1, j, k - 1, n) + 2.0 * fdata(i + 1, j, k, n) + fdata(i + 1, j, k + 1, n)) / 16.0;
                }
                else if (K == lo.z || K == hi.z) // Z face
                {
                    cdata(I, J, K, n) =
                        (+fdata(i - 1, j - 1, k, n) + 2.0 * fdata(i, j - 1, k, n) + fdata(i + 1, j - 1, k, n)
                            + 2.0 * fdata(i - 1, j, k, n) + 4.0 * fdata(i, j, k, n) + 2.0 * fdata(i + 1, j, k, n)
                            + fdata(i - 1, j + 1, k, n) + 2.0 * fdata(i, j + 1, k, n) + fdata(i + 1, j + 1, k, n)) / 16.0;
                }
                else // Interior
                    cdata(I, J, K, 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) + 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) + fdata(i + 1, j + 1, k + 1, n)) / 64.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) +
                        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) +
                        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
                    +
                    (fdata(i - 1, j, k, n) + fdata(i, j - 1, k, n) + fdata(i, j, k - 1, n) +
                        fdata(i + 1, j, k, n) + fdata(i, j + 1, k, n) + fdata(i, j, k + 1, n)) / 16.0
                    +
                    fdata(i, j, k, n) / 8.0;
            });
        }
    }
 
    if (need_parallel_copy) {
        crse.ParallelCopy(cfine);
    }
 
    crse.setMultiGhost(true);
    crse.FillBoundary(Geom(amrlev,cmglev).periodicity());
    nodalSync(amrlev, cmglev, crse);
}
 
void Operator<Grid::Node>::interpolation(int amrlev, int fmglev, MultiFab& fine, const MultiFab& crse) const
{
    BL_PROFILE("Operator::interpolation()");
    amrex::Box fdomain = m_geom[amrlev][fmglev].growPeriodicDomain(2);
    fdomain.convert(amrex::IntVect::TheNodeVector());
 
    bool need_parallel_copy = !amrex::isMFIterSafe(crse, fine);
    MultiFab cfine;
    const MultiFab* cmf = &crse;
    if (need_parallel_copy) {
        const BoxArray& ba = amrex::coarsen(fine.boxArray(), 2);
        cfine.define(ba, fine.DistributionMap(), crse.nComp(), crse.nGrow());
        cfine.ParallelCopy(crse);
        cmf = &cfine;
    }
 
    for (MFIter mfi(fine, false); mfi.isValid(); ++mfi)
    {
        Box fine_bx = mfi.validbox() & fdomain;
 
        const Box& course_bx = amrex::coarsen(fine_bx, 2);
        const Box& tmpbx = amrex::refine(course_bx, 2);
        FArrayBox tmpfab;
        tmpfab.resize(tmpbx, fine.nComp());
        tmpfab.setVal<amrex::RunOn::Device>(0.0);
        const amrex::FArrayBox& crsefab = (*cmf)[mfi];
 
        amrex::Array4<const amrex::Real> const& cdata = crsefab.const_array();
        amrex::Array4<amrex::Real> const& fdata = tmpfab.array();
 
        for (int n = 0; n < crse.nComp(); n++)
        {
            // I,J,K == coarse coordinates
            // i,j,k == fine coordinates
            amrex::ParallelFor(fine_bx, [=] AMREX_GPU_DEVICE(int i, int j, int k) {
 
                int I = i / 2, J = j / 2, K = k / 2;
 
                if (i % 2 == 0 && j % 2 == 0 && k % 2 == 0) // Coincident
                    fdata(i, j, k, n) = cdata(I, J, K, n);
                else if (j % 2 == 0 && k % 2 == 0) // X Edge
                    fdata(i, j, k, n) = 0.5 * (cdata(I, J, K, n) + cdata(I + 1, J, K, n));
                else if (k % 2 == 0 && i % 2 == 0) // Y Edge
                    fdata(i, j, k, n) = 0.5 * (cdata(I, J, K, n) + cdata(I, J + 1, K, n));
                else if (i % 2 == 0 && j % 2 == 0) // Z Edge
                    fdata(i, j, k, n) = 0.5 * (cdata(I, J, K, n) + cdata(I, J, K + 1, n));
                else if (i % 2 == 0) // X Face
                    fdata(i, j, k, n) = 0.25 * (cdata(I, J, K, n) + cdata(I, J + 1, K, n) +
                        cdata(I, J, K + 1, n) + cdata(I, J + 1, K + 1, n));
                else if (j % 2 == 0) // Y Face
                    fdata(i, j, k, n) = 0.25 * (cdata(I, J, K, n) + cdata(I, J, K + 1, n) +
                        cdata(I + 1, J, K, n) + cdata(I + 1, J, K + 1, n));
                else if (k % 2 == 0) // Z Face
                    fdata(i, j, k, n) = 0.25 * (cdata(I, J, K, n) + cdata(I + 1, J, K, n) +
                        cdata(I, J + 1, K, n) + cdata(I + 1, J + 1, K, n));
                else // Center
                    fdata(i, j, k, n) = 0.125 * (cdata(I, J, K, n) +
                        cdata(I + 1, J, K, n) + cdata(I, J + 1, K, n) + cdata(I, J, K + 1, n) +
                        cdata(I, J + 1, K + 1, n) + cdata(I + 1, J, K + 1, n) + cdata(I + 1, J + 1, K, n) +
                        cdata(I + 1, J + 1, K + 1, n));
 
            });
        }
        fine[mfi].plus<amrex::RunOn::Device>(tmpfab, fine_bx, fine_bx, 0, 0, fine.nComp());
    }
 
    fine.setMultiGhost(true);
    fine.FillBoundary(Geom(amrlev,fmglev).periodicity());
    nodalSync(amrlev, fmglev, fine);
}
 
void Operator<Grid::Node>::averageDownSolutionRHS(int camrlev, MultiFab& crse_sol, MultiFab& /*crse_rhs*/,
    const MultiFab& fine_sol, const MultiFab& /*fine_rhs*/)
{
    BL_PROFILE("Operator::averageDownSolutionRHS()");
    const auto& amrrr = AMRRefRatio(camrlev);
    amrex::average_down(fine_sol, crse_sol, 0, crse_sol.nComp(), amrrr);
 
    if (isSingular(0))
    {
        Util::Abort(INFO, "Singular operators not supported!");
    }
 
}
 
void Operator<Grid::Node>::interpolationAmr(int famrlev, MultiFab& fine,
    const MultiFab& crse, IntVect const& nghost) const
{
    BL_PROFILE("Operator::interpolationAmr()");
    if (!useQuadraticAmrInterpolation())
    {
        amrex::MLNodeLinOp::interpolationAmr(famrlev, fine, crse, nghost);
        return;
    }
    Util::Assert(INFO, TEST(AMRRefRatio(famrlev - 1) == 2));
    const int ncomp = getNComp();
 
    for (MFIter mfi(fine, false); mfi.isValid(); ++mfi)
    {
        Box fbx = mfi.tilebox();
        const Box valid = mfi.validbox();
        fbx.grow(nghost);
        const Dim3 vlo = amrex::lbound(valid), vhi = amrex::ubound(valid);
        Array4<Real> const& ffab = fine.array(mfi);
        Array4<Real const> const& cfab = crse.const_array(mfi);
 
        amrex::LoopConcurrentOnCpu(fbx, ncomp,
            [=] (int i, int j, int k, int n)
        {
            int ci[3][3] = {};
            Real cw[3][3] = {};
            int nc[3] = {1, 1, 1};
            cw[0][0] = cw[1][0] = cw[2][0] = 1.0;
 
            const int fi[3] = {i, j, k};
            const int flo[3] = {vlo.x, vlo.y, vlo.z};
            const int fhi[3] = {vhi.x, vhi.y, vhi.z};
            for (int d = 0; d < AMREX_SPACEDIM; ++d)
            {
                const int q = fi[d] >= 0 ? fi[d] / 2 : (fi[d] - 1) / 2;
                if (fi[d] % 2 == 0)
                {
                    ci[d][0] = q;
                }
                else if (fi[d] < flo[d])
                {
                    nc[d] = 3;
                    ci[d][0] = q;     cw[d][0] =  3.0 / 8.0;
                    ci[d][1] = q + 1; cw[d][1] =  3.0 / 4.0;
                    ci[d][2] = q + 2; cw[d][2] = -1.0 / 8.0;
                }
                else if (fi[d] > fhi[d])
                {
                    nc[d] = 3;
                    ci[d][0] = q - 1; cw[d][0] = -1.0 / 8.0;
                    ci[d][1] = q;     cw[d][1] =  3.0 / 4.0;
                    ci[d][2] = q + 1; cw[d][2] =  3.0 / 8.0;
                }
                else
                {
                    nc[d] = 2;
                    ci[d][0] = q;     cw[d][0] = 0.5;
                    ci[d][1] = q + 1; cw[d][1] = 0.5;
                }
            }
 
            Real value = 0.0;
            for (int a = 0; a < nc[0]; ++a)
                for (int b = 0; b < nc[1]; ++b)
                    for (int c = 0; c < nc[2]; ++c)
                        value += cw[0][a] * cw[1][b] * cw[2][c]
                            * cfab(ci[0][a], ci[1][b], ci[2][c], n);
            ffab(i, j, k, n) = value;
        });
    }
}
 
void Operator<Grid::Node>::realFillBoundary(MultiFab& phi, const Geometry& geom)
{
    Util::RealFillBoundary(phi, geom);
}
 
void Operator<Grid::Node>::applyBC(int amrlev, int mglev, MultiFab& phi, BCMode/* bc_mode*/,
    amrex::MLLinOp::StateMode /**/, bool skip_fillboundary) const
{
    BL_PROFILE("Operator::applyBC()");
 
    const Geometry& geom = m_geom[amrlev][mglev];
 
    if (!skip_fillboundary) {
        //phi.FillBoundary(geom.periodicity());
        //phi.setMultiGhost(true);
        phi.FillBoundaryAndSync(geom.periodicity());
    }
}
 
const amrex::FArrayBox&
Operator<Grid::Node>::GetFab(const int num, const int amrlev, const int mglev, const amrex::MFIter& mfi) const
{
    BL_PROFILE("Operator::GetFab()");
    Util::Message(INFO);
    return m_a_coeffs[num][amrlev][mglev][mfi];
}
 
void Operator<Grid::Node>::RegisterNewFab(amrex::Vector<amrex::MultiFab>& input)
{
    BL_PROFILE("Operator::RegisterNewFab()");
    Util::Message(INFO);
    /// \todo assertions here
    m_a_coeffs.resize(m_a_coeffs.size() + 1);
    m_a_coeffs[m_num_a_fabs].resize(m_num_amr_levels);
    for (int amrlev = 0; amrlev < m_num_amr_levels; ++amrlev)
    {
        m_a_coeffs[m_num_a_fabs][amrlev].resize(m_num_mg_levels[amrlev]);
        for (int mglev = 0; mglev < m_num_mg_levels[amrlev]; ++mglev)
            m_a_coeffs[m_num_a_fabs][amrlev][mglev].define(m_grids[amrlev][mglev],
                m_dmap[amrlev][mglev],
                input[amrlev].nComp(),
                input[amrlev].nGrow());
 
        amrex::MultiFab::Copy(m_a_coeffs[m_num_a_fabs][amrlev][0],
            input[amrlev], 0, 0,
            input[amrlev].nComp(),
            input[amrlev].nGrow());
    }
    m_num_a_fabs++;
}
 
 
void Operator<Grid::Node>::RegisterNewFab(amrex::Vector<std::unique_ptr<amrex::MultiFab> >& input)
{
    BL_PROFILE("Operator::RegisterNewFab()");
    Util::Message(INFO);
    /// \todo assertions here
    m_a_coeffs.resize(m_a_coeffs.size() + 1);
    m_a_coeffs[m_num_a_fabs].resize(m_num_amr_levels);
    for (int amrlev = 0; amrlev < m_num_amr_levels; ++amrlev)
    {
        m_a_coeffs[m_num_a_fabs][amrlev].resize(m_num_mg_levels[amrlev]);
        for (int mglev = 0; mglev < m_num_mg_levels[amrlev]; ++mglev)
            m_a_coeffs[m_num_a_fabs][amrlev][mglev].define(m_grids[amrlev][mglev],
                m_dmap[amrlev][mglev],
                input[amrlev]->nComp(),
                input[amrlev]->nGrow());
 
        amrex::MultiFab::Copy(m_a_coeffs[m_num_a_fabs][amrlev][0],
            *input[amrlev], 0, 0,
            input[amrlev]->nComp(),
            input[amrlev]->nGrow());
    }
    m_num_a_fabs++;
}
 
void Operator<Grid::Node>::reflux(int crse_amrlev,
    MultiFab& res, const MultiFab& /*crse_sol*/, const MultiFab& /*crse_rhs*/,
    MultiFab& fine_res, MultiFab& /*fine_sol*/, const MultiFab& /*fine_rhs*/) const
{
    BL_PROFILE("Operator::Elastic::reflux()");
 
    int ncomp = AMREX_SPACEDIM;
 
    amrex::Box cdomain(m_geom[crse_amrlev][0].growPeriodicDomain(2));
    cdomain.convert(amrex::IntVect::TheNodeVector());
 
    const Geometry& cgeom = m_geom[crse_amrlev][0];
 
    const BoxArray& fba = fine_res.boxArray();
    const DistributionMapping& fdm = fine_res.DistributionMap();
 
    MultiFab fine_res_for_coarse(amrex::coarsen(fba, 2), fdm, ncomp, 2);
    fine_res_for_coarse.ParallelCopy(res, 0, 0, ncomp, 0, 0, cgeom.periodicity());
 
    applyBC(crse_amrlev + 1, 0, fine_res, BCMode::Inhomogeneous, StateMode::Solution);
 
    /// \todo Replace with Enum
    // const int coarse_coarse_node = 0;
    const int coarse_fine_node = 1;
    const int fine_fine_node = 2;
 
    amrex::iMultiFab nodemask(amrex::coarsen(fba, 2), fdm, 1, 2);
    nodemask.ParallelCopy(*m_nd_fine_mask[crse_amrlev], 0, 0, 1, 0, 0, cgeom.periodicity());
 
    for (MFIter mfi(fine_res_for_coarse, false); mfi.isValid(); ++mfi)
    {
        const Box& bx = mfi.grownnodaltilebox(-1,1) & cdomain;
 
        amrex::Array4<const int> const& nmask = nodemask.array(mfi);
        //amrex::Array4<const int> const& cmask = cellmask.array(mfi);
 
        amrex::Array4<amrex::Real> const& cdata = fine_res_for_coarse.array(mfi);
        amrex::Array4<const amrex::Real> const& fdata = fine_res.array(mfi);
 
        const Dim3 lo = amrex::lbound(bx), hi = amrex::ubound(bx);
 
        for (int n = 0; n < fine_res.nComp(); n++)
        {
            // I,J,K == coarse coordinates
            // i,j,k == fine coordinates
            amrex::ParallelFor(bx, [=] AMREX_GPU_DEVICE(int I, int J, int K) {
                int i = I * 2, j = J * 2, k = K * 2;
 
                if (nmask(I, J, K) == fine_fine_node || nmask(I, J, K) == coarse_fine_node)
                {
                    if ((I == lo.x || I == hi.x) &&
                        (J == lo.y || J == hi.y) &&
                        (K == lo.z || K == hi.z)) // Corner
                        cdata(I, J, K, n) = fdata(i, j, k, n);
                    else if ((J == lo.y || J == hi.y) &&
                        (K == lo.z || K == hi.z)) // X edge
                        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);
                    else if ((K == lo.z || K == hi.z) &&
                        (I == lo.x || I == hi.x)) // Y edge
                        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);
                    else if ((I == lo.x || I == hi.x) &&
                        (J == lo.y || J == hi.y)) // Z edge
                        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);
                    else if (I == lo.x || I == hi.x) // X face
                        cdata(I, J, K, n) =
                        (+fdata(i, j - 1, k - 1, n) + 2.0 * fdata(i, j, k - 1, n) + fdata(i, j + 1, k - 1, n)
                            + 2.0 * fdata(i, j - 1, k, n) + 4.0 * fdata(i, j, k, n) + 2.0 * fdata(i, j + 1, k, n)
                            + fdata(i, j - 1, k + 1, n) + 2.0 * fdata(i, j, k + 1, n) + fdata(i, j + 1, k + 1, n)) / 16.0;
                    else if (J == lo.y || J == hi.y) // Y face
                        cdata(I, J, K, n) =
                        (+fdata(i - 1, j, k - 1, n) + 2.0 * fdata(i - 1, j, k, n) + fdata(i - 1, j, k + 1, n)
                            + 2.0 * fdata(i, j, k - 1, n) + 4.0 * fdata(i, j, k, n) + 2.0 * fdata(i, j, k + 1, n)
                            + fdata(i + 1, j, k - 1, n) + 2.0 * fdata(i + 1, j, k, n) + fdata(i + 1, j, k + 1, n)) / 16.0;
                    else if (K == lo.z || K == hi.z) // Z face
                        cdata(I, J, K, n) =
                        (+fdata(i - 1, j - 1, k, n) + 2.0 * fdata(i, j - 1, k, n) + fdata(i + 1, j - 1, k, n)
                            + 2.0 * fdata(i - 1, j, k, n) + 4.0 * fdata(i, j, k, n) + 2.0 * fdata(i + 1, j, k, n)
                            + fdata(i - 1, j + 1, k, n) + 2.0 * fdata(i, j + 1, k, n) + fdata(i + 1, j + 1, k, n)) / 16.0;
                    else // Interior
                        cdata(I, J, K, 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) + 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) + fdata(i + 1, j + 1, k + 1, n)) / 64.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) +
                            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) +
                            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
                        +
                        (fdata(i - 1, j, k, n) + fdata(i, j - 1, k, n) + fdata(i, j, k - 1, n) +
                            fdata(i + 1, j, k, n) + fdata(i, j + 1, k, n) + fdata(i, j, k + 1, n)) / 16.0
                        +
                        fdata(i, j, k, n) / 8.0;
                }
 
            });
        }
    }
 
    // Copy the fine residual restricted onto the coarse grid
    // into the final residual.
    res.ParallelCopy(fine_res_for_coarse, 0, 0, ncomp, 0, 0, cgeom.periodicity());
 
    // Sync up ghost nodes
    res.setMultiGhost(true);
    res.FillBoundaryAndSync(Geom(crse_amrlev).periodicity());
    nodalSync(crse_amrlev, 0, res);
 
    return;
}
 
void
Operator<Grid::Node>::solutionResidual(int amrlev, MultiFab& resid, MultiFab& x, const MultiFab& b,
    const MultiFab* /*crse_bcdata*/)
{
    const int mglev = 0;
    const int ncomp = b.nComp();
    apply(amrlev, mglev, resid, x, BCMode::Inhomogeneous, StateMode::Solution);
    MultiFab::Xpay(resid, -1.0, b, 0, 0, ncomp, 2);
    resid.setMultiGhost(true);
    resid.FillBoundaryAndSync(Geom(amrlev).periodicity());
}
 
void
Operator<Grid::Node>::correctionResidual(int amrlev, int mglev, MultiFab& resid, MultiFab& x, const MultiFab& b,
    BCMode /*bc_mode*/, const MultiFab* /*crse_bcdata*/)
{
    resid.setVal(0.0);
    apply(amrlev, mglev, resid, x, BCMode::Homogeneous, StateMode::Correction);
    int ncomp = b.nComp();
    MultiFab::Xpay(resid, -1.0, b, 0, 0, ncomp, resid.nGrow());
    resid.setMultiGhost(true);
    resid.FillBoundaryAndSync(Geom(amrlev).periodicity());
}
 
 
 
 
Operator<Grid::Cell>::Operator()
{
    m_ixtype = amrex::IntVect::TheCellVector();
}
 
void
Operator<Grid::Cell>::define(amrex::Vector<amrex::Geometry> a_geom,
    const amrex::Vector<amrex::BoxArray>& a_grids,
    const amrex::Vector<amrex::DistributionMapping>& a_dmap,
    BC::BC<Set::Scalar>& a_bc,
    const amrex::LPInfo& a_info,
    const amrex::Vector<amrex::FabFactory<amrex::FArrayBox> const*>& a_factory)
{
    m_bc = &a_bc;
 
    std::array<int, AMREX_SPACEDIM> is_periodic = m_bc->IsPeriodic();
 
    MLCellLinOp::define(a_geom, a_grids, a_dmap, a_info, a_factory);
 
    Util::Warning(INFO, "This section of code has not been tested.");
    for (int n = 0; n < getNComp(); n++)
    {
        m_lobc.push_back({ AMREX_D_DECL(is_periodic[0] ? amrex::LinOpBCType::Periodic : amrex::LinOpBCType::Dirichlet,
                        is_periodic[1] ? amrex::LinOpBCType::Periodic : amrex::LinOpBCType::Dirichlet,
                        is_periodic[2] ? amrex::LinOpBCType::Periodic : amrex::LinOpBCType::Dirichlet) });
        m_hibc.push_back({ AMREX_D_DECL(is_periodic[0] ? amrex::LinOpBCType::Periodic : amrex::LinOpBCType::Dirichlet,
                        is_periodic[1] ? amrex::LinOpBCType::Periodic : amrex::LinOpBCType::Dirichlet,
                        is_periodic[2] ? amrex::LinOpBCType::Periodic : amrex::LinOpBCType::Dirichlet) });
    }
 
    for (int ilev = 0; ilev < a_geom.size(); ++ilev)
        setLevelBC(ilev, nullptr);
 
}
 
 
void
Operator<Grid::Cell>::prepareForSolve()
{
    MLCellLinOp::prepareForSolve();
}
 
Operator<Grid::Cell>::BndryCondLoc::BndryCondLoc(const amrex::BoxArray& ba, const amrex::DistributionMapping& dm)
    : bcond(ba, dm),
    bcloc(ba, dm)
{
}
 
void
Operator<Grid::Cell>::BndryCondLoc::setLOBndryConds(const amrex::Geometry& /*geom*/, const amrex::Real* /*dx*/,
    const amrex::Array<BCType, AMREX_SPACEDIM>& /*lobc*/,
    const amrex::Array<BCType, AMREX_SPACEDIM>& /*hibc*/,
    int /*ratio*/, const amrex::RealVect& /*a_loc*/)
{
    Util::Warning(INFO, "This code has not been properlyt tested");
}
 
 
void
Operator<Grid::Cell>::averageDownCoeffs()
{
    for (int i = 0; i < m_num_a_fabs; i++)
    {
        for (int amrlev = m_num_amr_levels - 1; amrlev > 0; --amrlev)
        {
            auto& fine_a_coeffs = m_a_coeffs[i][amrlev];
            averageDownCoeffsSameAmrLevel(fine_a_coeffs);
        }
        averageDownCoeffsSameAmrLevel(m_a_coeffs[i][0]);
    }
}
 
void
Operator<Grid::Cell>::averageDownCoeffsSameAmrLevel(amrex::Vector<amrex::MultiFab>& a)
{
    int nmglevs = a.size();
    for (int mglev = 1; mglev < nmglevs; ++mglev)
    {
        amrex::average_down(a[mglev - 1], a[mglev], 0, a[0].nComp(), mg_coarsen_ratio);
    }
}
 
 
 
const amrex::FArrayBox&
Operator<Grid::Cell>::GetFab(const int num, const int amrlev, const int mglev, const amrex::MFIter& mfi) const
{
    return m_a_coeffs[num][amrlev][mglev][mfi];
}
 
 
void
Operator<Grid::Cell>::RegisterNewFab(amrex::Vector<amrex::MultiFab>& input)
{
    /// \todo assertions here
    m_a_coeffs.resize(m_a_coeffs.size() + 1);
    m_a_coeffs[m_num_a_fabs].resize(m_num_amr_levels);
    for (int amrlev = 0; amrlev < m_num_amr_levels; ++amrlev)
    {
        m_a_coeffs[m_num_a_fabs][amrlev].resize(m_num_mg_levels[amrlev]);
        for (int mglev = 0; mglev < m_num_mg_levels[amrlev]; ++mglev)
            m_a_coeffs[m_num_a_fabs][amrlev][mglev].define(m_grids[amrlev][mglev],
                m_dmap[amrlev][mglev],
                input[amrlev].nComp(),
                input[amrlev].nGrow());
 
        amrex::MultiFab::Copy(m_a_coeffs[m_num_a_fabs][amrlev][0],
            input[amrlev], 0, 0,
            input[amrlev].nComp(),
            input[amrlev].nGrow());
    }
    m_num_a_fabs++;
}
 
 
void
Operator<Grid::Cell>::RegisterNewFab(amrex::Vector<std::unique_ptr<amrex::MultiFab> >& input)
{
    /// \todo assertions here
    m_a_coeffs.resize(m_a_coeffs.size() + 1);
    m_a_coeffs[m_num_a_fabs].resize(m_num_amr_levels);
    for (int amrlev = 0; amrlev < m_num_amr_levels; ++amrlev)
    {
        m_a_coeffs[m_num_a_fabs][amrlev].resize(m_num_mg_levels[amrlev]);
        for (int mglev = 0; mglev < m_num_mg_levels[amrlev]; ++mglev)
            m_a_coeffs[m_num_a_fabs][amrlev][mglev].define(m_grids[amrlev][mglev],
                m_dmap[amrlev][mglev],
                input[amrlev]->nComp(),
                input[amrlev]->nGrow());
 
        amrex::MultiFab::Copy(m_a_coeffs[m_num_a_fabs][amrlev][0],
            *input[amrlev], 0, 0,
            input[amrlev]->nComp(),
            input[amrlev]->nGrow());
    }
    m_num_a_fabs++;
}
 
 
}