From 83d124d9bdc9125409d3702772af503bd5fbe3fb Mon Sep 17 00:00:00 2001 From: Candace Gilet Date: Thu, 20 Jun 2024 15:23:07 -0400 Subject: [PATCH] Add InitialPressureProjection to enforce hydrostatic equilibrium (#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 --- src/incflo.H | 9 ++- src/incflo.cpp | 4 ++ .../incflo_apply_nodal_projection.cpp | 41 +++++++++---- src/setup/init.cpp | 57 ++++++++++++------- test_no_eb_2d/benchmark.rayleigh_taylor | 3 + test_no_eb_3d/benchmark.rayleigh_taylor | 3 + 6 files changed, 84 insertions(+), 33 deletions(-) diff --git a/src/incflo.H b/src/incflo.H index 4f7babd2..8863c4d5 100644 --- a/src/incflo.H +++ b/src/incflo.H @@ -262,6 +262,11 @@ public: amrex::Real time, amrex::Real scaling_factor, bool incremental); void ApplyNodalProjection(amrex::Vector density, amrex::Real time, amrex::Real scaling_factor, bool incremental); + void ApplyNodalProjection(amrex::Vector const& density, + amrex::Vector vel, + amrex::Vector const& divu_Source, + amrex::Real time, amrex::Real scaling_factor, bool incremental, + bool set_inflow_bc); /////////////////////////////////////////////////////////////////////////// // @@ -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? @@ -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 (); diff --git a/src/incflo.cpp b/src/incflo.cpp index 1bf9a392..33823686 100644 --- a/src/incflo.cpp +++ b/src/incflo.cpp @@ -75,6 +75,10 @@ void incflo::InitData () if (m_do_initial_proj) { InitialProjection(); + + if (m_do_initial_pressure_proj) { + InitialPressureProjection(); + } } InitialIterations(); diff --git a/src/projection/incflo_apply_nodal_projection.cpp b/src/projection/incflo_apply_nodal_projection.cpp index cb061be5..1fb8ba93 100644 --- a/src/projection/incflo_apply_nodal_projection.cpp +++ b/src/projection/incflo_apply_nodal_projection.cpp @@ -69,6 +69,34 @@ void incflo::ApplyNodalProjection (Vector density, } } + Vector 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 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 const& density, + Vector vel, + Vector const& /*divu_Source*/, // only incompressible for now + Real time, Real scaling_factor, bool incremental, + bool set_inflow_bc) +{ Vector sigma(finest_level+1); if (!m_constant_density) { @@ -97,7 +125,6 @@ void incflo::ApplyNodalProjection (Vector density, auto bclo = get_nodal_projection_bc(Orientation::low); auto bchi = get_nodal_projection_bc(Orientation::high); - Vector vel; for (int lev = 0; lev <= finest_level; ++lev) { #ifdef AMREX_USE_EB if (m_eb_flow.enabled) { @@ -106,9 +133,8 @@ void incflo::ApplyNodalProjection (Vector 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 inflow_bcr; @@ -167,15 +193,6 @@ void incflo::ApplyNodalProjection (Vector 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(); diff --git a/src/setup/init.cpp b/src/setup/init.cpp index f0510c23..519477ca 100644 --- a/src/setup/init.cpp +++ b/src/setup/init.cpp @@ -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); @@ -345,26 +346,6 @@ void incflo::InitialProjection() { BL_PROFILE("incflo::InitialProjection()"); - // ************************************************************************************* - // Allocate space for the temporary MAC velocities - // ************************************************************************************* - Vector 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++) @@ -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 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 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 () diff --git a/test_no_eb_2d/benchmark.rayleigh_taylor b/test_no_eb_2d/benchmark.rayleigh_taylor index 6a35e829..5e79b3a7 100644 --- a/test_no_eb_2d/benchmark.rayleigh_taylor +++ b/test_no_eb_2d/benchmark.rayleigh_taylor @@ -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 diff --git a/test_no_eb_3d/benchmark.rayleigh_taylor b/test_no_eb_3d/benchmark.rayleigh_taylor index 6670613c..d465621e 100644 --- a/test_no_eb_3d/benchmark.rayleigh_taylor +++ b/test_no_eb_3d/benchmark.rayleigh_taylor @@ -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