diff --git a/Src/LinearSolvers/Projections/AMReX_MacProjector.H b/Src/LinearSolvers/Projections/AMReX_MacProjector.H index 5bc9c5a0f58..3070d5e59ed 100644 --- a/Src/LinearSolvers/Projections/AMReX_MacProjector.H +++ b/Src/LinearSolvers/Projections/AMReX_MacProjector.H @@ -95,6 +95,8 @@ private: void setOptions (); + void averageDownVelocity (); + std::unique_ptr m_abeclap; #ifdef AMREX_USE_EB std::unique_ptr m_eb_abeclap; @@ -107,6 +109,7 @@ private: Vector > m_umac; Vector m_rhs; Vector m_phi; + Vector m_divu; Vector > m_fluxes; Vector m_geom; diff --git a/Src/LinearSolvers/Projections/AMReX_MacProjector.cpp b/Src/LinearSolvers/Projections/AMReX_MacProjector.cpp index 88a451116da..dc69b69ec52 100644 --- a/Src/LinearSolvers/Projections/AMReX_MacProjector.cpp +++ b/Src/LinearSolvers/Projections/AMReX_MacProjector.cpp @@ -35,6 +35,7 @@ MacProjector::MacProjector (const Vector >& a_um m_rhs.resize(nlevs); m_phi.resize(nlevs); m_fluxes.resize(nlevs); + m_divu.resize(nlevs); #ifdef AMREX_USE_EB bool has_eb = a_umac[0][0]->hasEBFabFactory(); @@ -91,8 +92,14 @@ MacProjector::MacProjector (const Vector >& a_um } for (int ilev = 0, N = a_divu.size(); ilev < N; ++ilev) { - if (a_divu[ilev]) { - MultiFab::Copy(m_rhs[ilev], *a_divu[ilev], 0, 0, 1, 0); + if (a_divu[ilev]) + { +#ifdef AMREX_USE_EB + m_divu[ilev].define(ba[ilev],dm[ilev],1,0,MFInfo(),a_umac[ilev][0]->Factory()); +#else + m_divu[ilev].define(ba[ilev],dm[ilev],1,0); +#endif + MultiFab::Copy(m_divu[ilev], *a_divu[ilev], 0, 0, 1, 0); } } @@ -129,6 +136,8 @@ MacProjector::project (Real reltol, Real atol) { const int nlevs = m_rhs.size(); + averageDownVelocity(); + for (int ilev = 0; ilev < nlevs; ++ilev) { Array u; @@ -146,12 +155,20 @@ MacProjector::project (Real reltol, Real atol) m_umac[ilev][idim]->FillBoundary(m_geom[ilev].periodicity()); } } - + EB_computeDivergence(divu, u, m_geom[ilev], (m_umac_loc == MLMG::Location::FaceCentroid)); #else computeDivergence(divu, u, m_geom[ilev]); #endif - MultiFab::Subtract(m_rhs[ilev], divu, 0, 0, 1, 0); + + // Setup RHS as (m_divu - divu) where m_divu is a user-provided source term + MultiFab::Copy(m_rhs[ilev], divu, 0, 0, 1, 0); + m_rhs[ilev].mult(-1.0); + + if (m_divu[ilev].ok()) + { + MultiFab::Add(m_rhs[ilev],m_divu[ilev],0,0,1,0); + } } m_mlmg->solve(amrex::GetVecOfPtrs(m_phi), amrex::GetVecOfConstPtrs(m_rhs), reltol, atol); @@ -166,6 +183,8 @@ MacProjector::project (Real reltol, Real atol) #endif } } + + averageDownVelocity(); } void @@ -173,6 +192,8 @@ MacProjector::project (const Vector& phi_inout, Real reltol, Real ato { const int nlevs = m_rhs.size(); + averageDownVelocity(); + for (int ilev = 0; ilev < nlevs; ++ilev) { Array u; @@ -194,7 +215,14 @@ MacProjector::project (const Vector& phi_inout, Real reltol, Real ato #else computeDivergence(divu, u, m_geom[ilev]); #endif - MultiFab::Subtract(m_rhs[ilev], divu, 0, 0, 1, 0); + // Setup RHS as (m_divu - divu) where m_divu is a user-provided source term + MultiFab::Copy(m_rhs[ilev], divu, 0, 0, 1, 0); + m_rhs[ilev].mult(-1.0); + + if (m_divu[ilev].ok()) + { + MultiFab::Add(m_rhs[ilev],m_divu[ilev],0,0,1,0); + } MultiFab::Copy(m_phi[ilev], *phi_inout[ilev], 0, 0, 1, 0); } @@ -212,6 +240,9 @@ MacProjector::project (const Vector& phi_inout, Real reltol, Real ato #endif } } + + averageDownVelocity(); + for (int ilev = 0; ilev < nlevs; ++ilev) { MultiFab::Copy(*phi_inout[ilev], m_phi[ilev], 0, 0, 1, 0); @@ -294,4 +325,27 @@ MacProjector::setOptions () } } +void +MacProjector::averageDownVelocity () +{ + int finest_level = m_umac.size() - 1; + + + for (int lev = finest_level; lev > 0; --lev) + { + + IntVect rr = m_geom[lev].Domain().size() / m_geom[lev-1].Domain().size(); + +#ifdef AMREX_USE_EB + EB_average_down_faces(GetArrOfConstPtrs(m_umac[lev]), + m_umac[lev-1], + rr, m_geom[lev-1]); +#else + average_down_faces(GetArrOfConstPtrs(m_umac[lev]), + m_umac[lev-1], + rr, m_geom[lev-1]); +#endif + } +} + }