ERF_SolveWithEBMLMG.cpp File Reference

#include "ERF.H"
#include "ERF_EB.H"
#include "ERF_Utils.H"
#include "ERF_SolverUtils.H"
#include <AMReX_MLMG.H>
#include <AMReX_MLEBABecLap.H>

Include dependency graph for ERF_SolveWithEBMLMG.cpp:

Functions
template<typename T >
void	FillZeroAreaFaceFluxes (MultiFab &phi, Array< MultiFab, AMREX_SPACEDIM > &fluxes, const Geometry &geom, EBFArrayBoxFactory const &ebfact, T const &ebfact_u, T const &ebfact_v, T const &ebfact_w)
template<typename T >
void	solve_with_EB_mlmg (int lev, Vector< MultiFab > &rhs, Vector< MultiFab > &phi, Vector< Array< MultiFab, AMREX_SPACEDIM >> &fluxes, EBFArrayBoxFactory const &ebfact, T const &ebfact_u, T const &ebfact_v, T const &ebfact_w, const Geometry &geom, const Vector< amrex::IntVect > &ref_ratio, Array< std::string, 2 *AMREX_SPACEDIM > domain_bc_type, int mg_verbose, Real reltol, Real abstol)
template void	solve_with_EB_mlmg (int, Vector< MultiFab > &, Vector< MultiFab > &, Vector< Array< MultiFab, AMREX_SPACEDIM >> &, EBFArrayBoxFactory const &, EBFArrayBoxFactory const &, EBFArrayBoxFactory const &, EBFArrayBoxFactory const &, const Geometry &, const Vector< amrex::IntVect > &, Array< std::string, 2 *AMREX_SPACEDIM >, int, Real, Real)
template void	solve_with_EB_mlmg (int, Vector< MultiFab > &, Vector< MultiFab > &, Vector< Array< MultiFab, AMREX_SPACEDIM >> &, EBFArrayBoxFactory const &, eb_aux_ const &, eb_aux_ const &, eb_aux_ const &, const Geometry &, const Vector< amrex::IntVect > &, Array< std::string, 2 *AMREX_SPACEDIM >, int, Real, Real)

Function Documentation

FillZeroAreaFaceFluxes()

template<typename T >

void FillZeroAreaFaceFluxes	(	MultiFab &	phi,
		Array< MultiFab, AMREX_SPACEDIM > &	fluxes,
		const Geometry &	geom,
		EBFArrayBoxFactory const &	ebfact,
		T const &	ebfact_u,
		T const &	ebfact_v,
		T const &	ebfact_w
	)

Compute phi gradients where the area of the face is zero.

Template Parameters

T	EB factory or auxiliary EB data type for face-centered grids

Parameters

[in]	phi	Cell-centered solution used to compute gradients
[out]	fluxes	Face-centered gradient fluxes to fill
[in]	geom	Geometry used for inverse cell spacing
[in]	ebfact	Cell-centered embedded-boundary factory
[in]	ebfact_u	Embedded-boundary data on x-faces
[in]	ebfact_v	Embedded-boundary data on y-faces
[in]	ebfact_w	Embedded-boundary data on z-faces

Compute phi gradients where the area of the face is zero

At such a face one of the two adjacent cells is covered, so the gradient is extrapolated from the two nearest face gradients on the uncovered side, which reaches three cells away. Give phi three filled ghost cells to get that stencil everywhere; with fewer, faces near the edge of a box fall back to a first-order one-sided gradient, and the result then depends on the grid decomposition.

Template Parameters

T	EB factory or auxiliary EB data type for face-centered grids

Parameters

phi	Cell-centered solution used to compute gradients, with its ghost cells filled
fluxes	Face-centered gradient fluxes to fill
geom	Geometry used for inverse cell spacing
ebfact	Cell-centered embedded-boundary factory
ebfact_u	Embedded-boundary data on x-faces
ebfact_v	Embedded-boundary data on y-faces
ebfact_w	Embedded-boundary data on z-faces

{
    BL_PROFILE("ERF::FillZeroAreaFaceFluxes()");
 
    const GpuArray<Real, AMREX_SPACEDIM> dxInv = geom.InvCellSizeArray();
 
    for (MFIter mfi(phi,TileNoZ()); mfi.isValid(); ++mfi)
    {
        const Box& tbx = mfi.tilebox();
        const Box& xbx = mfi.nodaltilebox(0);
        const Box& ybx = mfi.nodaltilebox(1);
        const Box& zbx = mfi.nodaltilebox(2);
 
        // EBCellFlagFab const& cflag_fab = (ebfact.get_const_factory())->getMultiEBCellFlagFab()[mfi];
        EBCellFlagFab const& cflag_fab = (ebfact.getMultiEBCellFlagFab())[mfi];
        Array4<const EBCellFlag> cflag = cflag_fab.const_array();
 
        if (cflag_fab.getType(tbx) == FabType::singlevalued)
        {
            Array4<const Real> apx = ebfact.getAreaFrac()[0]->const_array(mfi);
            Array4<const Real> apy = ebfact.getAreaFrac()[1]->const_array(mfi);
            Array4<const Real> apz = ebfact.getAreaFrac()[2]->const_array(mfi);
 
            Array4<const EBCellFlag> u_cflag = ebfact_u.getMultiEBCellFlagFab()[mfi].const_array();
            Array4<const EBCellFlag> v_cflag = ebfact_v.getMultiEBCellFlagFab()[mfi].const_array();
            Array4<const EBCellFlag> w_cflag = ebfact_w.getMultiEBCellFlagFab()[mfi].const_array();
 
            Array4<Real const> const& p_arr = phi.const_array(mfi);
            Array4<Real> const& fx = fluxes[0].array(mfi);
            Array4<Real> const& fy = fluxes[1].array(mfi);
            Array4<Real> const& fz = fluxes[2].array(mfi);
 
            // The region of phi that we may read: this FAB's own cells, less any that lie
            // outside the domain in a non-periodic direction, since the FillBoundary done
            // before this call leaves those unfilled. The three-cell stencils below reach
            // three cells past the tile at a face on the box edge, so they must be checked
            // against this rather than assumed to fit (issue #3699)
            const Box rbx = phi[mfi].box() & geom.growPeriodicDomain(phi.nGrowVect());
            const auto rlo = lbound(rbx);
            const auto rhi = ubound(rbx);
 
            ParallelFor(xbx, ybx, zbx,
            // x-face
            [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
            {
                if (apx(i,j,k) == zero) {
                    if (!u_cflag(i,j,k).isCovered()) {
                        if (cflag(i,j,k).isCovered() && !cflag(i-1,j,k).isCovered()) {
                            // Extrapolate from the uncovered side, at i-1 and below
                            if (i-3 >= rlo.x) {
                                fx(i,j,k) = dxInv[0] * (p_arr(i-3,j,k) - three*p_arr(i-2,j,k) + two*p_arr(i-1,j,k));
                            } else if (i-2 >= rlo.x) {
                                fx(i,j,k) = dxInv[0] * (p_arr(i-1,j,k) - p_arr(i-2,j,k));
                            }
                        } else if (cflag(i-1,j,k).isCovered() && !cflag(i,j,k).isCovered()) {
                            // Extrapolate from the uncovered side, at i and above
                            if (i+2 <= rhi.x) {
                                fx(i,j,k) = dxInv[0] * (three*p_arr(i+1,j,k) - p_arr(i+2,j,k) - two*p_arr(i,j,k));
                            } else if (i+1 <= rhi.x) {
                                fx(i,j,k) = dxInv[0] * (p_arr(i+1,j,k) - p_arr(i,j,k));
                            }
                        }
                    }
                }
            },
            // y-face
            [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
            {
                if (apy(i,j,k) == zero) {
                    if (!v_cflag(i,j,k).isCovered()) {
                        if (cflag(i,j,k).isCovered() && !cflag(i,j-1,k).isCovered()) {
                            // Extrapolate from the uncovered side, at j-1 and below
                            if (j-3 >= rlo.y) {
                                fy(i,j,k) = dxInv[1] * (p_arr(i,j-3,k) - three*p_arr(i,j-2,k) + two*p_arr(i,j-1,k));
                            } else if (j-2 >= rlo.y) {
                                fy(i,j,k) = dxInv[1] * (p_arr(i,j-1,k) - p_arr(i,j-2,k));
                            }
                        } else if (cflag(i,j-1,k).isCovered() && !cflag(i,j,k).isCovered()) {
                            // Extrapolate from the uncovered side, at j and above
                            if (j+2 <= rhi.y) {
                                fy(i,j,k) = dxInv[1] * (three*p_arr(i,j+1,k) - p_arr(i,j+2,k) - two*p_arr(i,j,k));
                            } else if (j+1 <= rhi.y) {
                                fy(i,j,k) = dxInv[1] * (p_arr(i,j+1,k) - p_arr(i,j,k));
                            }
                        }
                    }
                }
            },
            // z-face
            [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
            {
                if (apz(i,j,k) == zero) {
                    if (!w_cflag(i,j,k).isCovered()) {
                        if (cflag(i,j,k).isCovered() && !cflag(i,j,k-1).isCovered()) {
                            // Extrapolate from the uncovered side, at k-1 and below
                            if (k-3 >= rlo.z) {
                                fz(i,j,k) = dxInv[2] * (p_arr(i,j,k-3) - three*p_arr(i,j,k-2) + two*p_arr(i,j,k-1));
                            } else if (k-2 >= rlo.z) {
                                fz(i,j,k) = dxInv[2] * (p_arr(i,j,k-1) - p_arr(i,j,k-2));
                            }
                        } else if (cflag(i,j,k-1).isCovered() && !cflag(i,j,k).isCovered()) {
                            // Extrapolate from the uncovered side, at k and above
                            if (k+2 <= rhi.z) {
                                fz(i,j,k) = dxInv[2] * (three*p_arr(i,j,k+1) - p_arr(i,j,k+2) - two*p_arr(i,j,k));
                            } else if (k+1 <= rhi.z) {
                                fz(i,j,k) = dxInv[2] * (p_arr(i,j,k+1) - p_arr(i,j,k));
                            }
                        }
                    }
                }
            });
        } // single-valued
    } // mfi
}

Referenced by solve_with_EB_mlmg().

Here is the call graph for this function:

Here is the caller graph for this function:

solve_with_EB_mlmg() [1/3]

template<typename T >

void solve_with_EB_mlmg	(	int	lev,
		Vector< MultiFab > &	rhs,
		Vector< MultiFab > &	phi,
		Vector< Array< MultiFab, AMREX_SPACEDIM >> &	fluxes,
		EBFArrayBoxFactory const &	ebfact,
		T const &	ebfact_u,
		T const &	ebfact_v,
		T const &	ebfact_w,
		const Geometry &	geom,
		const Vector< amrex::IntVect > &	ref_ratio,
		Array< std::string, 2 *AMREX_SPACEDIM >	domain_bc_type,
		int	mg_verbose,
		Real	reltol,
		Real	abstol
	)

Solve the Poisson equation using EB_enabled MLMG Note that the level may or may not be level zero

Important: we solve on the whole level even if there are disjoint regions

Template Parameters

T	EB factory or auxiliary EB data type for face-centered grids

Parameters

lev	Level index for the solve
rhs	Right-hand side MultiFab vector
phi	Solution MultiFab vector to fill
fluxes	Face-centered gradient fluxes to fill
ebfact	Cell-centered embedded-boundary factory
ebfact_u	Embedded-boundary data on x-faces
ebfact_v	Embedded-boundary data on y-faces
ebfact_w	Embedded-boundary data on z-faces
geom	Geometry for the solve level
ref_ratio	Coarse-fine refinement ratios
domain_bc_type	Domain boundary-condition names
mg_verbose	MLMG verbosity level
reltol	Relative solver tolerance
abstol	Absolute solver tolerance

{
    BL_PROFILE("ERF::solve_with_EB_mlmg()");
 
    auto const dom_lo = lbound(geom.Domain());
    auto const dom_hi = ubound(geom.Domain());
 
    LPInfo info;
    // Allow a hidden direction if the domain is one cell wide in any lateral direction
    if (dom_lo.x == dom_hi.x) {
        info.setHiddenDirection(0);
    } else if (dom_lo.y == dom_hi.y) {
        info.setHiddenDirection(1);
    }
 
    // Make sure the solver only sees the levels over which we are solving
    Vector<BoxArray>            ba_tmp;   ba_tmp.push_back(rhs[0].boxArray());
    Vector<DistributionMapping> dm_tmp;   dm_tmp.push_back(rhs[0].DistributionMap());
    Vector<Geometry>          geom_tmp; geom_tmp.push_back(geom);
 
    auto bclo = get_lo_projection_bc(geom,domain_bc_type);
    auto bchi = get_hi_projection_bc(geom,domain_bc_type);
 
    amrex::Print() << "BCLO " << bclo[0] << " " << bclo[1] << " " << bclo[2] << std::endl;
    amrex::Print() << "BCHI " << bchi[0] << " " << bchi[1] << " " << bchi[2] << std::endl;
 
    // ****************************************************************************
    // Multigrid solve
    // ****************************************************************************
 
    MLEBABecLap mleb (geom_tmp, ba_tmp, dm_tmp, info, {&ebfact});
 
    mleb.setMaxOrder(2);
    mleb.setDomainBC(bclo, bchi);
    if (lev > 0) {
        mleb.setCoarseFineBC(nullptr, ref_ratio[lev-1], LinOpBCType::Neumann);
    }
    mleb.setLevelBC(0, nullptr);
 
    //
    // This sets A = 0, B = 1 so that
    // the operator A alpha - b del dot beta grad to b
    // becomes  - del dot beta grad
    //
    mleb.setScalars(zero, one);
 
    Array<MultiFab,AMREX_SPACEDIM> bcoef;
    for (int idim = 0; idim < AMREX_SPACEDIM; ++idim) {
        bcoef[idim].define(convert(ba_tmp[0],IntVect::TheDimensionVector(idim)),
                            dm_tmp[0], 1, 0, MFInfo(), ebfact);
        bcoef[idim].setVal(-1.0);
    }
    mleb.setBCoeffs(0, amrex::GetArrOfConstPtrs(bcoef));
 
    MLMG mlmg(mleb);
 
    int max_iter = 100;
    mlmg.setMaxIter(max_iter);
    mlmg.setVerbose(mg_verbose);
    mlmg.setBottomVerbose(0);
 
    mlmg.solve(GetVecOfPtrs(phi), GetVecOfConstPtrs(rhs), reltol, abstol);
 
    mlmg.getFluxes(GetVecOfArrOfPtrs(fluxes));
 
    // Add the flux values (gradient phi) to the face-centered cell with zero area fraction
    // (mlmg.getFluxes does not fill gradient phi at faces with zero area fraction)
    //
    // That routine extrapolates the gradient onto a zero-area face from three cells on the
    // uncovered side, and at a face on the edge of a box the stencil reaches three cells
    // past the box -- well beyond the single ghost cell phi carries, and beyond whatever
    // the solve happened to leave in it (issue #3699). Do the extrapolation out of a
    // three-ghost-cell scratch copy instead of widening phi itself, which would cost MLMG
    // its alias fast path for the solution. Ghost cells outside the domain in a
    // non-periodic direction are still unfilled; FillZeroAreaFaceFluxes detects that and
    // drops to a lower order one-sided form there.
    MultiFab phi_ext(phi[0].boxArray(), phi[0].DistributionMap(), 1, 3);
    MultiFab::Copy(phi_ext, phi[0], 0, 0, 1, 0);
    phi_ext.FillBoundary(geom.periodicity());
 
    FillZeroAreaFaceFluxes(phi_ext, fluxes[0], geom, ebfact, ebfact_u, ebfact_v, ebfact_w);
 
    // ImposeBCsOnPhi(lev,phi[0], geom[lev].Domain());
 
    //
    // This arises because we solve MINUS del dot beta grad phi = div (rho u)
    //
    fluxes[0][0].mult(-one);
    fluxes[0][1].mult(-one);
    fluxes[0][2].mult(-one);
}

Here is the call graph for this function:

solve_with_EB_mlmg() [2/3]

template void solve_with_EB_mlmg	(	int	,
		Vector< MultiFab > &	,
		Vector< MultiFab > &	,
		Vector< Array< MultiFab, AMREX_SPACEDIM >> &	,
		EBFArrayBoxFactory const &	,
		eb_aux_ const &	,
		eb_aux_ const &	,
		eb_aux_ const &	,
		const Geometry &	,
		const Vector< amrex::IntVect > &	,
		Array< std::string, 2 *AMREX_SPACEDIM >	,
		int	,
		Real	,
		Real
	)

solve_with_EB_mlmg() [3/3]

template void solve_with_EB_mlmg	(	int	,
		Vector< MultiFab > &	,
		Vector< MultiFab > &	,
		Vector< Array< MultiFab, AMREX_SPACEDIM >> &	,
		EBFArrayBoxFactory const &	,
		EBFArrayBoxFactory const &	,
		EBFArrayBoxFactory const &	,
		EBFArrayBoxFactory const &	,
		const Geometry &	,
		const Vector< amrex::IntVect > &	,
		Array< std::string, 2 *AMREX_SPACEDIM >	,
		int	,
		Real	,
		Real
	)