Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Jump to bottom

More multiphase overset improvements #1182

Merged
merged 9 commits into from
Aug 5, 2024
Merged

More multiphase overset improvements #1182

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
9 commits
Select commit Hold shift + click to select a range
7c89771
change how fluxes are averaged to avoid segfault
mbkuhn Aug 2, 2024
af897bf
remove pressure initialization from sloshing tank
mbkuhn Aug 2, 2024
0132b5c
add p0 to p in post_init
mbkuhn Aug 5, 2024
c939701
Revert "add p0 to p in post_init"
mbkuhn Aug 5, 2024
e4bddf6
correct approach on first timestep
mbkuhn Aug 5, 2024
481b8b3
Merge branch 'main' into mphase_stokescase_fixes
mbkuhn Aug 5, 2024
d26b95b
removing unused function declarations
mbkuhn Aug 5, 2024
7ea3a66
better use of m_sim_ptr
mbkuhn Aug 5, 2024
0baf252
proper use of refRatio
mbkuhn Aug 5, 2024
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
6 changes: 0 additions & 6 deletions amr-wind/incflo.H
View file
Open in desktop
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
Open in desktop
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
Open in desktop
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
Open in desktop
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