Skip to content

Commit

Permalink
More multiphase overset improvements (#1182)
Browse files Browse the repository at this point in the history 
* proper use of refRatio
* remove pressure initialization from sloshing tank
* correct approach on first timestep
* removing unused function declarations
  • Loading branch information
@mbkuhn
mbkuhn authored Aug 5, 2024
1 parent 3e9194c commit 2fbd345
Show file tree
Hide file tree
Showing 4 changed files with 33 additions and 81 deletions.
6 changes: 0 additions & 6 deletions amr-wind/incflo.H
View file
Original file line number Diff line number Diff line change
Expand Up @@ -143,11 +143,6 @@ public:
amrex::Real scaling_factor,
bool incremental);

void UpdateGradP(
amrex::Vector<amrex::MultiFab const*> density,
amrex::Real time,
amrex::Real scaling_factor);

//! Initialize Physics instances as well as PDEs (include turbulence models)
void init_physics_and_pde();

Expand Down Expand Up @@ -204,7 +199,6 @@ private:
//
///////////////////////////////////////////////////////////////////////////

void set_background_pressure();
void ReadParameters();
void InitialProjection();
void InitialIterations();
Expand Down
61 changes: 32 additions & 29 deletions amr-wind/overset/OversetOps.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -30,9 +30,9 @@ void OversetOps::initialize(CFDSim& sim)

pp.query("verbose", m_verbose);

m_vof_exists = (*m_sim_ptr).repo().field_exists("vof");
m_vof_exists = m_sim_ptr->repo().field_exists("vof");
if (m_vof_exists) {
m_mphase = &(*m_sim_ptr).physics_manager().get<MultiPhase>();
m_mphase = &m_sim_ptr->physics_manager().get<MultiPhase>();

// Check combination of pressure options
if (m_mphase->perturb_pressure() &&
Expand All @@ -51,7 +51,7 @@ void OversetOps::initialize(CFDSim& sim)
}
}
if (m_replace_gradp_postsolve) {
m_gp_copy = &(*m_sim_ptr).repo().declare_field("gp_copy", 3);
m_gp_copy = &m_sim_ptr->repo().declare_field("gp_copy", 3);
}

parameter_output();
Expand All @@ -60,6 +60,14 @@ void OversetOps::initialize(CFDSim& sim)
void OversetOps::pre_advance_work()
{
if (!(m_vof_exists && m_use_hydrostatic_gradp)) {
if (m_mphase != nullptr) {
// Avoid modifying pressure upon initialization, assume pressure = 0
if (m_mphase->perturb_pressure() &&
m_sim_ptr->time().current_time() > 0.0) {
// Modify to be consistent with internal source terms
form_perturb_pressure();
}
}
// Update pressure gradient using updated overset pressure field
update_gradp();
}
Expand All @@ -71,20 +79,15 @@ void OversetOps::pre_advance_work()
// Use hydrostatic pressure gradient
set_hydrostatic_gradp();
} else {
if (m_mphase->perturb_pressure()) {
// Modify to be consistent with internal source terms
form_perturb_pressure();
}
// Update pressure gradient using sharpened pressure field
update_gradp();
}
}

// If pressure gradient will be replaced, store current pressure gradient
if (m_replace_gradp_postsolve) {
const auto& gp = (*m_sim_ptr).repo().get_field("gp");
for (int lev = 0; lev < (*m_sim_ptr).repo().num_active_levels();
++lev) {
const auto& gp = m_sim_ptr->repo().get_field("gp");
for (int lev = 0; lev < m_sim_ptr->repo().num_active_levels(); ++lev) {
amrex::MultiFab::Copy(
(*m_gp_copy)(lev), gp(lev), 0, 0, gp(lev).nComp(),
(m_gp_copy)->num_grow());
Expand All @@ -97,30 +100,30 @@ void OversetOps::update_gradp()
BL_PROFILE("amr-wind::OversetOps::update_gradp");

// Get relevant fields
auto& grad_p = (*m_sim_ptr).repo().get_field("gp");
auto& pressure = (*m_sim_ptr).repo().get_field("p");
auto& velocity = (*m_sim_ptr).repo().get_field("velocity");
auto& grad_p = m_sim_ptr->repo().get_field("gp");
auto& pressure = m_sim_ptr->repo().get_field("p");
auto& velocity = m_sim_ptr->repo().get_field("velocity");

// Set up projection object
std::unique_ptr<Hydro::NodalProjector> nodal_projector;
// Boundary conditions from field, periodicity
auto bclo = nodal_projection::get_projection_bc(
amrex::Orientation::low, pressure,
(*m_sim_ptr).mesh().Geom(0).isPeriodic());
m_sim_ptr->mesh().Geom(0).isPeriodic());
auto bchi = nodal_projection::get_projection_bc(
amrex::Orientation::high, pressure,
(*m_sim_ptr).mesh().Geom(0).isPeriodic());
m_sim_ptr->mesh().Geom(0).isPeriodic());
// Velocity multifab is needed for proper initialization, but only the size
// matters for the purpose of calculating gradp, the values do not matter
int finest_level = (*m_sim_ptr).mesh().finestLevel();
int finest_level = m_sim_ptr->mesh().finestLevel();
Vector<MultiFab*> vel;
for (int lev = 0; lev <= finest_level; ++lev) {
vel.push_back(&(velocity(lev)));
}
amr_wind::MLMGOptions options("nodal_proj");
// Create nodal projector with unity scaling factor for simplicity
nodal_projector = std::make_unique<Hydro::NodalProjector>(
vel, 1.0, (*m_sim_ptr).mesh().Geom(0, finest_level), options.lpinfo());
vel, 1.0, m_sim_ptr->mesh().Geom(0, finest_level), options.lpinfo());
// Set MLMG and NodalProjector options
options(*nodal_projector);
nodal_projector->setDomainBC(bclo, bchi);
Expand Down Expand Up @@ -203,9 +206,9 @@ void OversetOps::parameter_output() const

void OversetOps::sharpen_nalu_data()
{
const auto& repo = (*m_sim_ptr).repo();
const auto& repo = m_sim_ptr->repo();
const auto nlevels = repo.num_active_levels();
const auto geom = (*m_sim_ptr).mesh().Geom();
const auto geom = m_sim_ptr->mesh().Geom();

// Get blanking for cells
const auto& iblank_cell = repo.get_int_field("iblank_cell");
Expand Down Expand Up @@ -333,7 +336,7 @@ void OversetOps::sharpen_nalu_data()
for (int lev = nlevels - 1; lev > 0; --lev) {
amrex::average_down_faces(
GetArrOfConstPtrs(fluxes[lev]), fluxes[lev - 1],
repo.mesh().refRatio(lev), geom[lev - 1]);
repo.mesh().refRatio(lev - 1), geom[lev - 1]);
}

// Get pseudo dt (dtau)
Expand Down Expand Up @@ -382,7 +385,7 @@ void OversetOps::sharpen_nalu_data()
}

// Fillpatch for pressure to make sure pressure stencil has all points
p.fillpatch((*m_sim_ptr).time().current_time());
p.fillpatch(m_sim_ptr->time().current_time());

// Purely for debugging via visualization, should be removed later
// Currently set up to overwrite the levelset field (not used as time
Expand All @@ -395,19 +398,19 @@ void OversetOps::sharpen_nalu_data()

void OversetOps::form_perturb_pressure()
{
auto& pressure = (*m_sim_ptr).repo().get_field("p");
const auto& p0 = (*m_sim_ptr).repo().get_field("reference_pressure");
for (int lev = 0; lev < (*m_sim_ptr).repo().num_active_levels(); lev++) {
auto& pressure = m_sim_ptr->repo().get_field("p");
const auto& p0 = m_sim_ptr->repo().get_field("reference_pressure");
for (int lev = 0; lev < m_sim_ptr->repo().num_active_levels(); lev++) {
amrex::MultiFab::Subtract(
pressure(lev), p0(lev), 0, 0, 1, pressure.num_grow()[0]);
}
}

void OversetOps::set_hydrostatic_gradp()
{
const auto& repo = (*m_sim_ptr).repo();
const auto& repo = m_sim_ptr->repo();
const auto nlevels = repo.num_active_levels();
const auto geom = (*m_sim_ptr).mesh().Geom();
const auto geom = m_sim_ptr->mesh().Geom();

// Get blanking for cells
const auto& iblank_cell = repo.get_int_field("iblank_cell");
Expand All @@ -417,7 +420,7 @@ void OversetOps::set_hydrostatic_gradp()
auto& rho = repo.get_field("density");
auto& gp = repo.get_field("gp");
if (m_mphase->perturb_pressure()) {
rho0 = &((*m_sim_ptr).repo().get_field("reference_density"));
rho0 = &(m_sim_ptr->repo().get_field("reference_density"));
} else {
// Point to existing field, won't be used
rho0 = &rho;
Expand All @@ -433,11 +436,11 @@ void OversetOps::set_hydrostatic_gradp()

void OversetOps::replace_masked_gradp()
{
const auto& repo = (*m_sim_ptr).repo();
const auto& repo = m_sim_ptr->repo();
const auto nlevels = repo.num_active_levels();

// Get timestep
const amrex::Real dt = (*m_sim_ptr).time().delta_t();
const amrex::Real dt = m_sim_ptr->time().delta_t();

// Get blanking for cells
const auto& iblank_cell = repo.get_int_field("iblank_cell");
Expand Down
9 changes: 0 additions & 9 deletions amr-wind/physics/multiphase/SloshingTank.H
View file
Original file line number Diff line number Diff line change
Expand Up @@ -35,7 +35,6 @@ public:
private:
Field& m_velocity;
Field& m_levelset;
Field& m_pressure;

//! Initial free surface amplitude magnitude
amrex::Real m_amplitude{0.1};
Expand All @@ -45,14 +44,6 @@ private:

//! Initial zero-level free-surface water depth
amrex::Real m_waterlevel{0.0};

//! Stuff to get from MultiPhase physics
amrex::Real m_rho1{std::numeric_limits<amrex::Real>::quiet_NaN()};
amrex::Real m_rho2{std::numeric_limits<amrex::Real>::quiet_NaN()};
amrex::Vector<amrex::Real> m_gravity{
std::numeric_limits<amrex::Real>::quiet_NaN(),
std::numeric_limits<amrex::Real>::quiet_NaN(),
std::numeric_limits<amrex::Real>::quiet_NaN()};
};

} // namespace amr_wind
Expand Down
38 changes: 1 addition & 37 deletions amr-wind/physics/multiphase/SloshingTank.cpp
View file
Original file line number Diff line number Diff line change
@@ -1,5 +1,4 @@
#include "amr-wind/physics/multiphase/SloshingTank.H"
#include "amr-wind/physics/multiphase/MultiPhase.H"
#include "amr-wind/CFDSim.H"
#include "AMReX_ParmParse.H"
#include "amr-wind/fvm/gradient.H"
Expand All @@ -10,30 +9,22 @@ namespace amr_wind {
SloshingTank::SloshingTank(CFDSim& sim)
: m_velocity(sim.repo().get_field("velocity"))
, m_levelset(sim.repo().get_field("levelset"))
, m_pressure(sim.repo().get_field("p"))
{
amrex::ParmParse pp(identifier());
pp.query("amplitude", m_amplitude);
pp.query("peak_enhance", m_kappa);
pp.query("water_level", m_waterlevel);

const auto& mphase = sim.physics_manager().get<MultiPhase>();
m_rho1 = mphase.rho1();
m_rho2 = mphase.rho2();
m_gravity = mphase.gravity();
}

/** Initialize the velocity and levelset fields at the beginning of the
* simulation. Initializing pressure has little effect for pure AMR-Wind
* simulations, but it improves initialization in the overset solver.
* simulation.
*/
void SloshingTank::initialize_fields(int level, const amrex::Geometry& geom)
{
auto& velocity = m_velocity(level);
velocity.setVal(0.0);

auto& levelset = m_levelset(level);
auto& pressure = m_pressure(level);
const auto& dx = geom.CellSizeArray();
const auto& problo = geom.ProbLoArray();
const auto& probhi = geom.ProbHiArray();
Expand All @@ -42,14 +33,10 @@ void SloshingTank::initialize_fields(int level, const amrex::Geometry& geom)
const amrex::Real water_level = m_waterlevel;
const amrex::Real Lx = probhi[0] - problo[0];
const amrex::Real Ly = probhi[1] - problo[1];
const amrex::Real rho1 = m_rho1;
const amrex::Real rho2 = m_rho2;
const amrex::Real grav_z = m_gravity[2];

for (amrex::MFIter mfi(levelset); mfi.isValid(); ++mfi) {
const auto& vbx = mfi.validbox();
auto phi = levelset.array(mfi);
auto p = pressure.array(mfi);

amrex::ParallelFor(
vbx, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept {
Expand All @@ -63,29 +50,6 @@ void SloshingTank::initialize_fields(int level, const amrex::Geometry& geom)
std::pow(y - problo[1] - 0.5 * Ly, 2)));
phi(i, j, k) = z0 - z;
});
amrex::Box const& nbx = mfi.grownnodaltilebox();
amrex::ParallelFor(
nbx, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept {
// For pressure nodes, no offset
const amrex::Real x = problo[0] + i * dx[0];
const amrex::Real y = problo[1] + j * dx[1];
const amrex::Real z = problo[2] + k * dx[2];
const amrex::Real z0 =
water_level +
Amp * std::exp(
-kappa * (std::pow(x - problo[0] - 0.5 * Lx, 2) +
std::pow(y - problo[1] - 0.5 * Ly, 2)));
// Integrated (top-down in z) phase heights to pressure node
const amrex::Real ih_g =
amrex::max(0.0, amrex::min(probhi[2] - z0, probhi[2] - z));
const amrex::Real ih_l =
amrex::max(0.0, amrex::min(z0 - z, z0 - problo[2]));
// Integrated rho at pressure node
const amrex::Real irho = rho1 * ih_l + rho2 * ih_g;

// Add term to reference pressure
p(i, j, k) = -irho * grav_z;
});
}
}

Expand Down

0 comments on commit 2fbd345

Please sign in to comment.