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
void	FillZeroAreaFaceFluxes (MultiFab &phi, Array< MultiFab, AMREX_SPACEDIM > &fluxes, const Geometry &geom, EBFArrayBoxFactory const &ebfact, eb_aux_ const &ebfact_u, eb_aux_ const &ebfact_v, eb_aux_ const &ebfact_w)
void	solve_with_EB_mlmg (int lev, Vector< MultiFab > &rhs, Vector< MultiFab > &phi, Vector< Array< MultiFab, AMREX_SPACEDIM >> &fluxes, EBFArrayBoxFactory const &ebfact, eb_aux_ const &ebfact_u, eb_aux_ const &ebfact_v, eb_aux_ 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)

Function Documentation

FillZeroAreaFaceFluxes()

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

Compute phi gradients where the area of the face is zero

{
    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) == Real(0.0)) {
                    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) - 3.*p_arr(i-2,j,k) + 2.*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] * (3.*p_arr(i+1,j,k) - p_arr(i+2,j,k) - 2.*p_arr(i,j,k));
                        }
                    }
                }
            },
            // y-face
            [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
            {
                if (apy(i,j,k) == Real(0.0)) {
                    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) - 3.*p_arr(i,j-2,k) + 2.*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] * (3.*p_arr(i,j+1,k) - p_arr(i,j+2,k) - 2.*p_arr(i,j,k));
                        }
                    }
                }
            },
            // z-face
            [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
            {
                if (apz(i,j,k) == Real(0.0)) {
                    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) - 3.*p_arr(i,j,k-2) + 2.*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] * (3.*p_arr(i,j,k+1) - p_arr(i,j,k+2) - 2.*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()

void solve_with_EB_mlmg	(	int	lev,
		Vector< MultiFab > &	rhs,
		Vector< MultiFab > &	phi,
		Vector< Array< MultiFab, AMREX_SPACEDIM >> &	fluxes,
		EBFArrayBoxFactory const &	ebfact,
		eb_aux_ const &	ebfact_u,
		eb_aux_ const &	ebfact_v,
		eb_aux_ 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 0.

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

{
    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(0.0, 1.0);
 
    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(-1.);
    fluxes[0][1].mult(-1.);
    fluxes[0][2].mult(-1.);
}

Here is the call graph for this function: