Skip to content

Commit

Permalink
Add InitialPressureProjection to enforce hydrostatic equilibrium
Browse files Browse the repository at this point in the history 
in cases where that's not done through the background density and
pressure.

Right now, this is a runtime option that is off by default. It does
not check for any background state; it is up to the user to know
when to use it and to ensure that the background density is turned
off with incflo.ro_0 = 0
  • Loading branch information
@cgilet
cgilet committed Jun 7, 2024
1 parent 0a40f11 commit 8c1dda4
Show file tree
Hide file tree
Showing 4 changed files with 78 additions and 33 deletions.
9 changes: 8 additions & 1 deletion src/incflo.H
View file
Original file line number Diff line number Diff line change
Expand Up @@ -252,6 +252,11 @@ public:
amrex::Real time, amrex::Real scaling_factor, bool incremental);
void ApplyNodalProjection(amrex::Vector<amrex::MultiFab const*> density,
amrex::Real time, amrex::Real scaling_factor, bool incremental);
void ApplyNodalProjection(amrex::Vector<amrex::MultiFab const*> density,
amrex::Vector<amrex::MultiFab *> vel,
amrex::Vector<amrex::MultiFab *> divu_Source,
amrex::Real time, amrex::Real scaling_factor, bool incremental,
bool set_inflow_bc);

///////////////////////////////////////////////////////////////////////////
//
Expand Down Expand Up @@ -345,7 +350,8 @@ private:
amrex::Real m_dt_change_max = amrex::Real(1.1);

// Initial projection / iterations
bool m_do_initial_proj = true;
bool m_do_initial_proj = true;
bool m_do_initial_pressure_proj = false;
int m_initial_iterations = 3;

// Use Boussinesq approximation for buoyancy?
Expand Down Expand Up @@ -901,6 +907,7 @@ private:
void ReadIOParameters ();
void ResizeArrays (); // Resize arrays to fit (up to) max_level + 1 AMR levels
void InitialProjection ();
void InitialPressureProjection ();
void InitialIterations ();
#ifdef AMREX_USE_EB
void InitialRedistribution ();
Expand Down
4 changes: 4 additions & 0 deletions src/incflo.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -75,6 +75,10 @@ void incflo::InitData ()

if (m_do_initial_proj) {
InitialProjection();

if (m_do_initial_pressure_proj) {
InitialPressureProjection();
}
}

InitialIterations();
Expand Down
41 changes: 29 additions & 12 deletions src/projection/incflo_apply_nodal_projection.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -69,6 +69,34 @@ void incflo::ApplyNodalProjection (Vector<MultiFab const*> density,
}
}

Vector<MultiFab*> vel;
for (int lev = 0; lev <= finest_level; ++lev) {
vel.push_back(&(m_leveldata[lev]->velocity));
}

// Cell-centered divergence constraint source term, which is always zero for now
Vector<MultiFab* > Source(finest_level+1, nullptr);

bool set_inflow_bc = !proj_for_small_dt && !incremental;
ApplyNodalProjection(density, vel, Source, time, scaling_factor,
incremental, set_inflow_bc);

// Define "vel" to be U^{n+1} rather than (U^{n+1}-U^n)
if (proj_for_small_dt || incremental)
{
for (int lev = 0; lev <= finest_level; ++lev) {
MultiFab::Add(m_leveldata[lev]->velocity,
m_leveldata[lev]->velocity_o, 0, 0, AMREX_SPACEDIM, 0);
}
}
}

void incflo::ApplyNodalProjection (Vector<MultiFab const*> density,
Vector<MultiFab *> vel,
Vector<MultiFab *> /*divu_Source*/, // only incompressible for now
Real time, Real scaling_factor, bool incremental,
bool set_inflow_bc)
{
Vector<amrex::MultiFab> sigma(finest_level+1);
if (!m_constant_density)
{
Expand Down Expand Up @@ -97,7 +125,6 @@ void incflo::ApplyNodalProjection (Vector<MultiFab const*> density,
auto bclo = get_nodal_projection_bc(Orientation::low);
auto bchi = get_nodal_projection_bc(Orientation::high);

Vector<MultiFab*> vel;
for (int lev = 0; lev <= finest_level; ++lev) {
#ifdef AMREX_USE_EB
if (m_eb_flow.enabled) {
Expand All @@ -106,9 +133,8 @@ void incflo::ApplyNodalProjection (Vector<MultiFab const*> density,
set_eb_tracer(lev, time, *get_tracer_eb()[lev], 1);
}
#endif
vel.push_back(&(m_leveldata[lev]->velocity));
vel[lev]->setBndry(0.0);
if (!proj_for_small_dt && !incremental) {
if (set_inflow_bc) {
// Only the inflow boundary gets set here
IntVect nghost(1);
amrex::Vector<amrex::BCRec> inflow_bcr;
Expand Down Expand Up @@ -167,15 +193,6 @@ void incflo::ApplyNodalProjection (Vector<MultiFab const*> density,

nodal_projector->project(m_nodal_mg_rtol, m_nodal_mg_atol);

// Define "vel" to be U^{n+1} rather than (U^{n+1}-U^n)
if (proj_for_small_dt || incremental)
{
for (int lev = 0; lev <= finest_level; ++lev) {
MultiFab::Add(m_leveldata[lev]->velocity,
m_leveldata[lev]->velocity_o, 0, 0, AMREX_SPACEDIM, 0);
}
}

// Get phi and fluxes
auto phi = nodal_projector->getPhi();
auto gradphi = nodal_projector->getGradPhi();
Expand Down
57 changes: 37 additions & 20 deletions src/setup/init.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -38,6 +38,7 @@ void incflo::ReadParameters ()
pp.query("steady_state_tol", m_steady_state_tol);
pp.query("initial_iterations", m_initial_iterations);
pp.query("do_initial_proj", m_do_initial_proj);
pp.query("do_initial_pressure_proj", m_do_initial_pressure_proj);

pp.query("fixed_dt", m_fixed_dt);
pp.query("cfl", m_cfl);
Expand Down Expand Up @@ -342,26 +343,6 @@ void incflo::InitialProjection()
{
BL_PROFILE("incflo::InitialProjection()");

// *************************************************************************************
// Allocate space for the temporary MAC velocities
// *************************************************************************************
Vector<MultiFab> u_mac_tmp(finest_level+1), v_mac_tmp(finest_level+1), w_mac_tmp(finest_level+1);
int ngmac = nghost_mac();

for (int lev = 0; lev <= finest_level; ++lev) {
AMREX_D_TERM(u_mac_tmp[lev].define(amrex::convert(grids[lev],IntVect::TheDimensionVector(0)), dmap[lev],
1, ngmac, MFInfo(), Factory(lev));,
v_mac_tmp[lev].define(amrex::convert(grids[lev],IntVect::TheDimensionVector(1)), dmap[lev],
1, ngmac, MFInfo(), Factory(lev));,
w_mac_tmp[lev].define(amrex::convert(grids[lev],IntVect::TheDimensionVector(2)), dmap[lev],
1, ngmac, MFInfo(), Factory(lev)););
if (ngmac > 0) {
AMREX_D_TERM(u_mac_tmp[lev].setBndry(0.0);,
v_mac_tmp[lev].setBndry(0.0);,
w_mac_tmp[lev].setBndry(0.0););
}
}

Real dummy_dt = 1.0;
bool incremental = false;
for (int lev = 0; lev <= finest_level; lev++)
Expand All @@ -378,6 +359,42 @@ void incflo::InitialProjection()
}
}

// Project to enforce hydrostatic equilibrium
void incflo::InitialPressureProjection()
{
BL_PROFILE("incflo::InitialPressureProjection()");

if (m_verbose > 0) { Print() << " Initial pressure projection \n"; }

Real dummy_dt = 1.0;
int nGhost = 1;

// fixme??? are density ghosts fill already???
//I think we only need to worry about this if doing outflow bcs...
for (int lev = 0; lev <= finest_level; lev++)
{
m_leveldata[lev]->density.FillBoundary(geom[lev].periodicity());
}

// Set the velocity to the gravity field
Vector<MultiFab> vel(finest_level + 1);
for (int lev = 0; lev <= finest_level; ++lev) {
vel[lev].define(grids[lev], dmap[lev], AMREX_SPACEDIM, nGhost,
MFInfo(), *m_factory[lev]);
for (int idim = 0; idim < AMREX_SPACEDIM; ++idim) {
vel[lev].setVal(m_gravity[idim], idim, 1, 1);
}
}

// Cell-centered divergence condition source term
// Always zero this here
Vector<MultiFab*> Source(finest_level+1, nullptr);

ApplyNodalProjection(get_density_new_const(), GetVecOfPtrs(vel), Source,
m_cur_time, dummy_dt, false /*incremental*/,
true /*set_inflow_bc*/);
}

#ifdef AMREX_USE_EB
void
incflo::InitialRedistribution ()
Expand Down

0 comments on commit 8c1dda4

Please sign in to comment.