From 9411fd47c5d73ff80bbf65180565669b0a62b694 Mon Sep 17 00:00:00 2001 From: Weiqun Zhang Date: Tue, 11 Oct 2022 15:42:00 -0700 Subject: [PATCH] Fix MLMG::getGradSolution & getFluxes for inhomogeneous Neumann and Robin BC Because of the way how inhomogeneous and Robin BC are handled, we must add the inhomogeneous fluxes back, otherwise they would be zero at those boundaries. --- Src/LinearSolvers/MLMG/AMReX_MLCellABecLap.H | 4 + .../MLMG/AMReX_MLCellABecLap.cpp | 111 ++++++++++++++++++ Src/LinearSolvers/MLMG/AMReX_MLCellLinOp.H | 5 + Src/LinearSolvers/MLMG/AMReX_MLCellLinOp.cpp | 2 + 4 files changed, 122 insertions(+) diff --git a/Src/LinearSolvers/MLMG/AMReX_MLCellABecLap.H b/Src/LinearSolvers/MLMG/AMReX_MLCellABecLap.H index 8849a2be292..0cc6456b7c8 100644 --- a/Src/LinearSolvers/MLMG/AMReX_MLCellABecLap.H +++ b/Src/LinearSolvers/MLMG/AMReX_MLCellABecLap.H @@ -61,6 +61,10 @@ public: virtual void applyInhomogNeumannTerm (int amrlev, Any& rhs) const final override; + virtual void addInhomogNeumannFlux ( + int amrlev, const Array& grad, + MultiFab const& sol, bool mult_bcoef) const final override; + virtual void applyOverset (int amlev, Any& rhs) const override; #if defined(AMREX_USE_HYPRE) && (AMREX_SPACEDIM > 1) diff --git a/Src/LinearSolvers/MLMG/AMReX_MLCellABecLap.cpp b/Src/LinearSolvers/MLMG/AMReX_MLCellABecLap.cpp index af094d89406..db57162c21f 100644 --- a/Src/LinearSolvers/MLMG/AMReX_MLCellABecLap.cpp +++ b/Src/LinearSolvers/MLMG/AMReX_MLCellABecLap.cpp @@ -189,6 +189,7 @@ MLCellABecLap::getFluxes (const Vector >& a_flux a_flux[alev][idim]->mult(betainv); } } + addInhomogNeumannFlux(alev, a_flux[alev], *a_sol[alev], true); } } @@ -416,6 +417,116 @@ MLCellABecLap::applyInhomogNeumannTerm (int amrlev, Any& a_rhs) const } } +void +MLCellABecLap::addInhomogNeumannFlux ( + int amrlev, const Array& grad, MultiFab const& sol, + bool mult_bcoef) const +{ + /* + * if (mult_bcoef == true) + * grad is -bceof*grad phi + * else + * grad is grad phi + */ + Real fac = mult_bcoef ? Real(-1.0) : Real(1.0); + + bool has_inhomog_neumann = hasInhomogNeumannBC(); + bool has_robin = hasRobinBC(); + + if (!has_inhomog_neumann && !has_robin) return; + + int ncomp = getNComp(); + const int mglev = 0; + + const auto dxinv = m_geom[amrlev][mglev].InvCellSize(); + const Box domain = m_geom[amrlev][mglev].growPeriodicDomain(1); + + Array bcoef = {AMREX_D_DECL(nullptr,nullptr,nullptr)}; + if (mult_bcoef) { + bcoef = getBCoeffs(amrlev,mglev); + } + + const auto& bndry = *m_bndry_sol[amrlev]; + + MFItInfo mfi_info; + if (Gpu::notInLaunchRegion()) mfi_info.SetDynamic(true); + +#ifdef AMREX_USE_OMP +#pragma omp parallel if (Gpu::notInLaunchRegion()) +#endif + for (MFIter mfi(sol, mfi_info); mfi.isValid(); ++mfi) + { + Box const& vbx = mfi.validbox(); + for (OrientationIter orit; orit.isValid(); ++orit) { + const Orientation ori = orit(); + const int idim = ori.coordDir(); + const Box& ccb = amrex::adjCell(vbx, ori); + const Dim3 os = IntVect::TheDimensionVector(idim).dim3(); + const Real dxi = dxinv[idim]; + if (! domain.contains(ccb)) { + for (int icomp = 0; icomp < ncomp; ++icomp) { + auto const& phi = sol.const_array(mfi,icomp); + auto const bv = bndry.bndryValues(ori).multiFab().const_array(mfi,icomp); + auto const bc = bcoef[idim] ? bcoef[idim]->const_array(mfi,icomp) + : Array4{}; + auto const& f = grad[idim]->array(mfi,icomp); + if (ori.isLow()) { + if (m_lobc_orig[icomp][idim] == + LinOpBCType::inhomogNeumann) { + AMREX_HOST_DEVICE_FOR_3D(ccb, i, j, k, + { + int ii = i+os.x; + int jj = j+os.y; + int kk = k+os.z; + Real b = bc ? bc(ii,jj,kk) : Real(1.0); + f(ii,jj,kk) = fac*b*bv(i,j,k); + }); + } else if (m_lobc_orig[icomp][idim] == + LinOpBCType::Robin) { + Array4 const& rbc = (*m_robin_bcval[amrlev])[mfi].const_array(icomp*3); + AMREX_HOST_DEVICE_FOR_3D(ccb, i, j, k, + { + int ii = i+os.x; + int jj = j+os.y; + int kk = k+os.z; + Real tmp = Real(1.0) / + (rbc(i,j,k,1)*dxi + rbc(i,j,k,0)*Real(0.5)); + Real RA = rbc(i,j,k,2) * tmp; + Real RB = (rbc(i,j,k,1)*dxi - rbc(i,j,k,0)*Real(0.5)) * tmp; + Real b = bc ? bc(ii,jj,kk) : Real(1.0); + f(ii,jj,kk) = fac*b*dxi*((Real(1.0)-RB)*phi(ii,jj,kk)-RA); + }); + } + } else { + if (m_hibc_orig[icomp][idim] == + LinOpBCType::inhomogNeumann) { + AMREX_HOST_DEVICE_FOR_3D(ccb, i, j, k, + { + Real b = bc ? bc(i,j,k) : Real(1.0); + f(i,j,k) = fac*b*bv(i,j,k); + }); + } else if (m_hibc_orig[icomp][idim] == + LinOpBCType::Robin) { + Array4 const& rbc = (*m_robin_bcval[amrlev])[mfi].const_array(icomp*3); + AMREX_HOST_DEVICE_FOR_3D(ccb, i, j, k, + { + Real tmp = Real(1.0) / + (rbc(i,j,k,1)*dxi + rbc(i,j,k,0)*Real(0.5)); + Real RA = rbc(i,j,k,2) * tmp; + Real RB = (rbc(i,j,k,1)*dxi - rbc(i,j,k,0)*Real(0.5)) * tmp; + Real b = bc ? bc(i,j,k) : Real(1.0); + f(i,j,k) = fac*b*dxi*(RA+(RB-Real(1.0))* + phi(i-os.x,j-os.y,k-os.z)); + }); + } + } + } + } + } + } +} + + void MLCellABecLap::applyOverset (int amrlev, Any& a_rhs) const { diff --git a/Src/LinearSolvers/MLMG/AMReX_MLCellLinOp.H b/Src/LinearSolvers/MLMG/AMReX_MLCellLinOp.H index 457f7565df3..9a6bb222113 100644 --- a/Src/LinearSolvers/MLMG/AMReX_MLCellLinOp.H +++ b/Src/LinearSolvers/MLMG/AMReX_MLCellLinOp.H @@ -164,6 +164,11 @@ public: virtual void AnyAverageDownAndSync (Vector& sol) const override; + virtual void addInhomogNeumannFlux (int /*amrlev*/, + const Array& /*grad*/, + MultiFab const& /*sol*/, + bool /*mult_bcoef*/) const {} + struct BCTL { BoundCond type; Real location; diff --git a/Src/LinearSolvers/MLMG/AMReX_MLCellLinOp.cpp b/Src/LinearSolvers/MLMG/AMReX_MLCellLinOp.cpp index e4c9cef953f..5c8edcbb1a6 100644 --- a/Src/LinearSolvers/MLMG/AMReX_MLCellLinOp.cpp +++ b/Src/LinearSolvers/MLMG/AMReX_MLCellLinOp.cpp @@ -938,6 +938,8 @@ MLCellLinOp::compGrad (int amrlev, const Array& grad, }); #endif } + + addInhomogNeumannFlux(amrlev, grad, sol, false); } void