Skip to content

Commit

Permalink
Add InitialPressureProjection to enforce hydrostatic equilibrium (AMR…
Browse files Browse the repository at this point in the history 
…eX-Fluids#102)

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

---------

Co-authored-by: Ann Almgren <asalmgren@lbl.gov>
  • Loading branch information
@cgilet @asalmgren
cgilet and asalmgren authored Jun 20, 2024
1 parent a101f71 commit 83d124d
Show file tree
Hide file tree
Showing 6 changed files with 84 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 @@ -262,6 +262,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*> const& density,
amrex::Vector<amrex::MultiFab *> vel,
amrex::Vector<amrex::MultiFab *> const& divu_Source,
amrex::Real time, amrex::Real scaling_factor, bool incremental,
bool set_inflow_bc);

///////////////////////////////////////////////////////////////////////////
//
Expand Down Expand Up @@ -355,7 +360,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 @@ -914,6 +920,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*> const& density,
Vector<MultiFab *> vel,
Vector<MultiFab *> const& /*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 @@ -345,26 +346,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 @@ -381,6 +362,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
3 changes: 3 additions & 0 deletions test_no_eb_2d/benchmark.rayleigh_taylor
View file
Original file line number Diff line number Diff line change
Expand Up @@ -37,6 +37,9 @@ amr.plt_vort = 1
incflo.advection_type = "Godunov"

incflo.constant_density = false
incflo.ro_0 = 0 # Turn off background state
incflo.do_initial_pressure_proj = 1 # to enforce hydrostatic equilibrium

incflo.advect_tracer = true

incflo.diffusion_type = 1 # 0 = Explicit, 1 = Crank-Nicolson, 2 = Implicit
Expand Down
3 changes: 3 additions & 0 deletions test_no_eb_3d/benchmark.rayleigh_taylor
View file
Original file line number Diff line number Diff line change
Expand Up @@ -37,6 +37,9 @@ amr.plt_vort = 1
incflo.advection_type = "Godunov"

incflo.constant_density = false
incflo.ro_0 = 0 # Turn off background state
incflo.do_initial_pressure_proj = 1 # to enforce hydrostatic equilibrium

incflo.advect_tracer = true

incflo.diffusion_type = 1 # 0 = Explicit, 1 = Crank-Nicolson, 2 = Implicit
Expand Down

0 comments on commit 83d124d

Please sign in to comment.