Skip to content

Commit

Permalink
Multiphase Hybrid - FY24 Q2 milestone (sharpening the pressure field) (…
Browse files Browse the repository at this point in the history 
…#1099)
  • Loading branch information
@mbkuhn
mbkuhn authored Jul 22, 2024
1 parent 5ebb2ab commit fe73b70
Show file tree
Hide file tree
Showing 8 changed files with 1,320 additions and 412 deletions.
1 change: 1 addition & 0 deletions amr-wind/overset/CMakeLists.txt
View file
Original file line number Diff line number Diff line change
Expand Up @@ -2,4 +2,5 @@ target_sources(${amr_wind_lib_name}
PRIVATE
TiogaInterface.cpp
OversetOps.cpp
overset_ops_routines.cpp
)
8 changes: 3 additions & 5 deletions amr-wind/overset/OversetOps.H
View file
Original file line number Diff line number Diff line change
Expand Up @@ -21,17 +21,15 @@ private:
// Functions called within public functions
void parameter_output() const;
void sharpen_nalu_data();
void form_perturb_pressure();
void set_hydrostatic_gradp();
void replace_masked_gradp();

// Check for multiphase sim
bool m_vof_exists{false};
// Check for perturbational pressure
// (will be removed soon)
bool m_perturb_p{false};

// This is the only option for now
const bool m_use_hydrostatic_gradp{true};
// To avoid using sharpened pressure, use hydrostatic pressure
bool m_use_hydrostatic_gradp{false};

// Coupling options
bool m_disable_nodal_proj{false};
Expand Down
126 changes: 94 additions & 32 deletions amr-wind/overset/OversetOps.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -22,19 +22,33 @@ void OversetOps::initialize(CFDSim& sim)
pp.query("reinit_target_cutoff", m_target_cutoff);

// Queries for coupling options
pp.query("use_hydrostatic_gradp", m_use_hydrostatic_gradp);
pp.query("replace_gradp_postsolve", m_replace_gradp_postsolve);
// OversetOps does not control these coupling options, merely reports them
pp.query("disable_coupled_nodal_proj", m_disable_nodal_proj);
pp.query("disable_coupled_mac_proj", m_disable_mac_proj);

pp.query("verbose", m_verbose);

amrex::ParmParse pp_icns("ICNS");
pp_icns.query("use_perturb_pressure", m_perturb_p);

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

// Check combination of pressure options
if (m_mphase->perturb_pressure() &&
!m_mphase->reconstruct_true_pressure()) {
amrex::Abort(
"OversetOps: perturbational pressure is turned on, but true "
"pressure reconstruction is turned off. This approach will be "
"incorrect when coupling with Nalu-Wind.");
}
if (!m_mphase->perturb_pressure() &&
m_mphase->reconstruct_true_pressure()) {
amrex::Print()
<< "WARNING (OversetOps): true pressure reconstruction is "
"turned on, but it will remain inactive because "
"perturbational pressure is turned off.\n";
}
}
if (m_replace_gradp_postsolve) {
m_gp_copy = &(*m_sim_ptr).repo().declare_field("gp_copy", 3);
Expand All @@ -45,7 +59,6 @@ void OversetOps::initialize(CFDSim& sim)

void OversetOps::pre_advance_work()
{
// Pressure gradient not updated for current multiphase approach
if (!(m_vof_exists && m_use_hydrostatic_gradp)) {
// Update pressure gradient using updated overset pressure field
update_gradp();
Expand All @@ -57,12 +70,19 @@ void OversetOps::pre_advance_work()
if (m_use_hydrostatic_gradp) {
// 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) {
auto& gp = (*m_sim_ptr).repo().get_field("gp");
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(
Expand Down Expand Up @@ -154,6 +174,11 @@ void OversetOps::parameter_output() const
<< "---- Replace overset pres grad: "
<< m_replace_gradp_postsolve << std::endl;
if (m_vof_exists) {
amrex::Print() << "---- Perturbational pressure : "
<< m_mphase->perturb_pressure() << std::endl
<< "---- Reconstruct true pressure: "
<< m_mphase->reconstruct_true_pressure()
<< std::endl;
amrex::Print() << "Overset Reinitialization Parameters:\n"
<< "---- Maximum iterations : " << m_n_iterations
<< std::endl
Expand Down Expand Up @@ -190,29 +215,50 @@ void OversetOps::sharpen_nalu_data()
auto& levelset = repo.get_field("levelset");
auto& rho = repo.get_field("density");
auto& velocity = repo.get_field("velocity");
auto& gp_noghost = repo.get_field("gp");
auto& p = repo.get_field("p");

// Create scratch fields for fluxes
// 5 components are vof, density, and 3 of velocity
auto flux_x = repo.create_scratch_field(5, 0, amr_wind::FieldLoc::XFACE);
auto flux_y = repo.create_scratch_field(5, 0, amr_wind::FieldLoc::YFACE);
auto flux_z = repo.create_scratch_field(5, 0, amr_wind::FieldLoc::ZFACE);
// Create scratch field for approximate signed distance function and grad
// (components 0-2 are gradient, 3 is asdf)
auto normal_vec = repo.create_scratch_field(3, vof.num_grow()[0] - 1);
// 9 components are vof, density, 3 of velocity, 3 of gp, and psource flag
auto flux_x = repo.create_scratch_field(9, 1, amr_wind::FieldLoc::XFACE);
auto flux_y = repo.create_scratch_field(9, 1, amr_wind::FieldLoc::YFACE);
auto flux_z = repo.create_scratch_field(9, 1, amr_wind::FieldLoc::ZFACE);

auto p_src = repo.create_scratch_field(1, 0, amr_wind::FieldLoc::NODE);
auto normal_vec = repo.create_scratch_field(3, vof.num_grow()[0] - 1);
auto target_vof = repo.create_scratch_field(1, vof.num_grow()[0]);

// Create levelset field
// Sharpening fluxes (at faces) have 1 ghost, requiring fields to have >= 2
auto gp_scr = repo.create_scratch_field(3, 2);
auto& gp = *gp_scr;

// Give initial max possible value of pseudo-velocity scale
const auto dx_lev0 = (geom[0]).CellSizeArray();
const amrex::Real max_pvscale =
std::min(std::min(dx_lev0[0], dx_lev0[1]), dx_lev0[2]);
amrex::Real pvscale = max_pvscale;

// Prep things that do not change with iterations
for (int lev = 0; lev < nlevels; ++lev) {
// Thickness used here is user parameter, whatever works best
auto dx = (geom[lev]).CellSizeArray();
const auto dx = (geom[lev]).CellSizeArray();
const amrex::Real i_th =
m_relative_length_scale * std::cbrt(dx[0] * dx[1] * dx[2]);

// Populate approximate signed distance function
overset_ops::populate_psi(levelset(lev), vof(lev), i_th, m_asdf_tiny);

gp(lev).setVal(0.0);
amrex::MultiFab::Copy(gp(lev), gp_noghost(lev), 0, 0, 3, 0);
gp(lev).FillBoundary(geom[lev].periodicity());

// Get pseudo-velocity scale, proportional to smallest dx in iblank
const amrex::Real pvscale_lev =
overset_ops::calculate_pseudo_velocity_scale(
iblank_cell(lev), dx, max_pvscale);
pvscale = std::min(pvscale, pvscale_lev);
}
amrex::Gpu::synchronize();
amrex::ParallelDescriptor::ReduceRealMin(pvscale);

// Convert levelset to vof to get target_vof
m_mphase->levelset2vof(iblank_cell, *target_vof);
Expand Down Expand Up @@ -262,13 +308,14 @@ void OversetOps::sharpen_nalu_data()
// Sharpening fluxes for vof, density, and momentum
overset_ops::populate_sharpen_fluxes(
(*flux_x)(lev), (*flux_y)(lev), (*flux_z)(lev), vof(lev),
(*target_vof)(lev), (*normal_vec)(lev), velocity(lev),
m_upw_margin, m_mphase->rho1(), m_mphase->rho2());
(*target_vof)(lev), (*normal_vec)(lev), velocity(lev), gp(lev),
rho(lev), pvscale, m_upw_margin, m_mphase->rho1(),
m_mphase->rho2());

// Process fluxes
overset_ops::process_fluxes(
(*flux_x)(lev), (*flux_y)(lev), (*flux_z)(lev),
iblank_cell(lev));
overset_ops::process_fluxes_calc_src(
(*flux_x)(lev), (*flux_y)(lev), (*flux_z)(lev), (*p_src)(lev),
vof(lev), iblank_cell(lev));

// Measure convergence to determine if loop can stop
if (calc_convg) {
Expand All @@ -282,11 +329,9 @@ void OversetOps::sharpen_nalu_data()

// Average down fluxes across levels for consistency
for (int lev = nlevels - 1; lev > 0; --lev) {
amrex::IntVect rr =
geom[lev].Domain().size() / geom[lev - 1].Domain().size();
amrex::average_down_faces(
GetArrOfConstPtrs(fluxes[lev]), fluxes[lev - 1], rr,
geom[lev - 1]);
GetArrOfConstPtrs(fluxes[lev]), fluxes[lev - 1],
repo.mesh().refRatio(lev), geom[lev - 1]);
}

// Get pseudo dt (dtau)
Expand All @@ -305,15 +350,19 @@ void OversetOps::sharpen_nalu_data()

// Apply fluxes
for (int lev = 0; lev < nlevels; ++lev) {
const auto dx = (geom[lev]).CellSizeArray();

overset_ops::apply_fluxes(
(*flux_x)(lev), (*flux_y)(lev), (*flux_z)(lev), vof(lev),
rho(lev), velocity(lev), ptfac, m_vof_tol);
(*flux_x)(lev), (*flux_y)(lev), (*flux_z)(lev), (*p_src)(lev),
vof(lev), rho(lev), velocity(lev), gp(lev), p(lev), dx, ptfac,
m_vof_tol);

vof(lev).FillBoundary(geom[lev].periodicity());
velocity(lev).FillBoundary(geom[lev].periodicity());
gp(lev).FillBoundary(geom[lev].periodicity());
}
amrex::Gpu::synchronize();

// Fillpatch for ghost cells
vof.fillpatch((*m_sim_ptr).time().current_time());
velocity.fillpatch((*m_sim_ptr).time().current_time());
amrex::Gpu::synchronize();

// Update density (fillpatch built in)
m_mphase->set_density_via_vof();
Expand All @@ -329,6 +378,9 @@ void OversetOps::sharpen_nalu_data()
}
}

// Fillpatch for pressure to make sure pressure stencil has all points
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
// evolves) with the post-sharpening velocity magnitude
Expand All @@ -338,6 +390,16 @@ void OversetOps::sharpen_nalu_data()
amrex::Gpu::synchronize();
}

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++) {
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();
Expand All @@ -351,7 +413,7 @@ void OversetOps::set_hydrostatic_gradp()
Field* rho0{nullptr};
auto& rho = repo.get_field("density");
auto& gp = repo.get_field("gp");
if (m_perturb_p) {
if (m_mphase->perturb_pressure()) {
rho0 = &((*m_sim_ptr).repo().get_field("reference_density"));
} else {
// Point to existing field, won't be used
Expand All @@ -362,7 +424,7 @@ void OversetOps::set_hydrostatic_gradp()
for (int lev = 0; lev < nlevels; ++lev) {
overset_ops::replace_gradp_hydrostatic(
gp(lev), rho(lev), (*rho0)(lev), iblank_cell(lev),
m_mphase->gravity()[2], m_perturb_p);
m_mphase->gravity()[2], m_mphase->perturb_pressure());
}
}

Expand Down
Loading

0 comments on commit fe73b70

Please sign in to comment.