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

phi	Cell-centered solution used to compute gradients
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);
 
            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()) {
                            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 (cflag(i-1,j,k).isCovered() && !cflag(i,j,k).isCovered()) {
                            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));
                        }
                    }
                }
            },
            // 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()) {
                            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 (cflag(i,j-1,k).isCovered() && !cflag(i,j,k).isCovered()) {
                            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));
                        }
                    }
                }
            },
            // 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()) {
                            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 (cflag(i,j,k-1).isCovered() && !cflag(i,j,k).isCovered()) {
                            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));
                        }
                    }
                }
            });
        } // 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)
    FillZeroAreaFaceFluxes(phi[0], 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
	)