Skip to content

Commit

Permalink
MacProjector: allow for re-use of the object and enhance multi-level …
Browse files Browse the repository at this point in the history
…algorithm. (AMReX-Codes#1225)

* MacProjector: build RHS each time project() is called

* MacProjector: average down velocities before and after projection
  • Loading branch information
mic84 authored Aug 7, 2020
1 parent b9ed579 commit 6de33d5
Show file tree
Hide file tree
Showing 2 changed files with 62 additions and 5 deletions.
3 changes: 3 additions & 0 deletions Src/LinearSolvers/Projections/AMReX_MacProjector.H
Original file line number Diff line number Diff line change
Expand Up @@ -95,6 +95,8 @@ private:

void setOptions ();

void averageDownVelocity ();

std::unique_ptr<MLABecLaplacian> m_abeclap;
#ifdef AMREX_USE_EB
std::unique_ptr<MLEBABecLap> m_eb_abeclap;
Expand All @@ -107,6 +109,7 @@ private:
Vector<Array<MultiFab*,AMREX_SPACEDIM> > m_umac;
Vector<MultiFab> m_rhs;
Vector<MultiFab> m_phi;
Vector<MultiFab> m_divu;
Vector<Array<MultiFab,AMREX_SPACEDIM> > m_fluxes;

Vector<Geometry> m_geom;
Expand Down
64 changes: 59 additions & 5 deletions Src/LinearSolvers/Projections/AMReX_MacProjector.cpp
Original file line number Diff line number Diff line change
Expand Up @@ -35,6 +35,7 @@ MacProjector::MacProjector (const Vector<Array<MultiFab*,AMREX_SPACEDIM> >& 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();
Expand Down Expand Up @@ -91,8 +92,14 @@ MacProjector::MacProjector (const Vector<Array<MultiFab*,AMREX_SPACEDIM> >& 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);
}
}

Expand Down Expand Up @@ -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<MultiFab const*, AMREX_SPACEDIM> u;
Expand All @@ -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);
Expand All @@ -166,13 +183,17 @@ MacProjector::project (Real reltol, Real atol)
#endif
}
}

averageDownVelocity();
}

void
MacProjector::project (const Vector<MultiFab*>& phi_inout, Real reltol, Real atol)
{
const int nlevs = m_rhs.size();

averageDownVelocity();

for (int ilev = 0; ilev < nlevs; ++ilev)
{
Array<MultiFab const*, AMREX_SPACEDIM> u;
Expand All @@ -194,7 +215,14 @@ MacProjector::project (const Vector<MultiFab*>& 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);
}
Expand All @@ -212,6 +240,9 @@ MacProjector::project (const Vector<MultiFab*>& 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);
Expand Down Expand Up @@ -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
}
}

}

0 comments on commit 6de33d5

Please sign in to comment.