Skip to content

Commit

Permalink
vof operations changes: (#1037)
Browse files Browse the repository at this point in the history 
* deprecate vof-only sharpening
* add replace mask
* make VOF advection into an operator applied from n to n+1
  • Loading branch information
@mbkuhn
mbkuhn authored Jun 5, 2024
1 parent 02feaf5 commit eb2df80
Show file tree
Hide file tree
Showing 4 changed files with 58 additions and 190 deletions.
48 changes: 30 additions & 18 deletions amr-wind/equation_systems/vof/vof_advection.H
View file
Original file line number Diff line number Diff line change
Expand Up @@ -27,6 +27,7 @@ struct AdvectionOp<VOF, fvm::Godunov>
{
amrex::ParmParse pp_multiphase("VOF");
pp_multiphase.query("remove_debris", m_rm_debris);
pp_multiphase.query("replace_masked", m_replace_mask);

// Setup density factor arrays for multiplying velocity flux
fields_in.repo.declare_face_normal_field(
Expand All @@ -47,13 +48,25 @@ struct AdvectionOp<VOF, fvm::Godunov>

auto& repo = fields.repo;
const auto& geom = repo.mesh().Geom();
const int nlevels = repo.num_active_levels();

auto& aa_x = repo.get_field("advalpha_x");
auto& aa_y = repo.get_field("advalpha_y");
auto& aa_z = repo.get_field("advalpha_z");

// cppcheck-suppress constVariable
auto& dof_field = fields.field;
// Old and new states
auto& dof_old = fields.field.state(amr_wind::FieldState::Old);
auto& dof_new = fields.field;
// Working state of vof is nph, to keep others untouched during step
auto& dof_field = fields.field.state(amr_wind::FieldState::NPH);

// Initialize as old state values
for (int lev = 0; lev < repo.num_active_levels(); ++lev) {
amrex::MultiFab::Copy(
dof_field(lev), dof_old(lev), 0, 0, dof_field.num_comp(),
dof_field.num_grow());
}

//
// Advect volume using Implicit Eulerian Sweeping method with PLIC
// reconstruction
Expand All @@ -75,8 +88,6 @@ struct AdvectionOp<VOF, fvm::Godunov>
isweep = 1;
}

const int nlevels = repo.num_active_levels();

amrex::Vector<amrex::Array<amrex::MultiFab*, AMREX_SPACEDIM>> fluxes(
nlevels);
amrex::Vector<amrex::Array<amrex::MultiFab*, AMREX_SPACEDIM>> advas(
Expand All @@ -90,20 +101,6 @@ struct AdvectionOp<VOF, fvm::Godunov>
advas[lev][2] = &aa_z(lev);
}

// Sharpen acquired vof field if hybrid solver is being used
if (repo.int_field_exists("iblank_cell")) {
auto& f_iblank = repo.get_int_field("iblank_cell");
amr_wind::multiphase::sharpen_acquired_vof(
nlevels, f_iblank, dof_field);
// Update old vof so that old density can be updated
// cppcheck-suppress constVariableReference
auto& old_dof_field = dof_field.state(amr_wind::FieldState::Old);
for (int lev = 0; lev < repo.num_active_levels(); ++lev) {
amrex::MultiFab::Copy(
old_dof_field(lev), dof_field(lev), 0, 0,
dof_field.num_comp(), dof_field.num_grow());
}
}
// Split advection step 1, with cmask calculation
multiphase::split_advection_step(
isweep, 0, nlevels, dof_field, fluxes, (*fluxC), advas, u_mac,
Expand All @@ -116,6 +113,20 @@ struct AdvectionOp<VOF, fvm::Godunov>
multiphase::split_advection_step(
isweep, 2, nlevels, dof_field, fluxes, (*fluxC), advas, u_mac,
v_mac, w_mac, dof_field.bc_type(), geom, dt, m_rm_debris);

// Replace masked cells using overset
if (repo.int_field_exists("iblank_cell") && m_replace_mask) {
auto& f_iblank = repo.get_int_field("iblank_cell");
amr_wind::multiphase::replace_masked_vof(
nlevels, f_iblank, dof_field, dof_new);
}

// Copy working version of vof to new state
for (int lev = 0; lev < repo.num_active_levels(); ++lev) {
amrex::MultiFab::Copy(
dof_new(lev), dof_field(lev), 0, 0, dof_field.num_comp(),
dof_field.num_grow());
}
}

PDEFields& fields;
Expand All @@ -124,6 +135,7 @@ struct AdvectionOp<VOF, fvm::Godunov>
Field& w_mac;
int isweep = 0;
bool m_rm_debris{true};
bool m_replace_mask{true};
// Lagrangian transport is deprecated, only Eulerian is supported
};

Expand Down
26 changes: 9 additions & 17 deletions amr-wind/equation_systems/vof/vof_hybridsolver_ops.H
View file
Original file line number Diff line number Diff line change
Expand Up @@ -7,39 +7,31 @@

namespace amr_wind::multiphase {

AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE amrex::Real sharpen_kernel(
int i,
int j,
int k,
amrex::Array4<amrex::Real const> const& volfrac) noexcept
{
return std::max(
0.0,
std::min(
1.0,
0.5 + (volfrac(i, j, k) < 0.5 ? -1.0 : 1.0) *
std::pow(std::abs(volfrac(i, j, k) - 0.5), 1.0 / 3.0)));
}

static void sharpen_acquired_vof(
const int nlevels, amr_wind::IntField& f_iblank, amr_wind::Field& f_vof)
static void replace_masked_vof(
const int nlevels,
amr_wind::IntField& f_iblank,
amr_wind::Field& f_vof,
amr_wind::Field& f_vof_new)
{
// Sharpen data from nalu-wind (in iblank regions)
for (int lev = 0; lev < nlevels; ++lev) {
auto& iblank = f_iblank(lev);
auto& vof = f_vof(lev);
const auto& vof_new = f_vof_new(lev);

for (amrex::MFIter mfi(iblank); mfi.isValid(); ++mfi) {
const auto& gbx = mfi.growntilebox();
const amrex::Array4<const int>& native_flag =
iblank.const_array(mfi);
const amrex::Array4<amrex::Real>& volfrac = vof.array(mfi);
const amrex::Array4<const amrex::Real>& vfmasked =
vof_new.const_array(mfi);
amrex::ParallelFor(
gbx, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept {
// In iblanked regions, sharpen VOF and limit it
volfrac(i, j, k) = (native_flag(i, j, k) > 0)
? volfrac(i, j, k)
: sharpen_kernel(i, j, k, volfrac);
: vfmasked(i, j, k);
});
}
}
Expand Down
6 changes: 6 additions & 0 deletions unit_tests/multiphase/test_vof_cons.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -222,6 +222,12 @@ class VOFConsTest : public MeshTest
seqn.initialize();

for (int n = 0; n < niter; ++n) {
// Copy new to old to prep for advection
for (int lev = 0; lev < repo.num_active_levels(); ++lev) {
amrex::MultiFab::Copy(
vof.state(amr_wind::FieldState::Old)(lev), vof(lev), 0, 0,
vof.num_comp(), vof.num_grow());
}
// Perform VOF solve
seqn.compute_advection_term(amr_wind::FieldState::Old);
seqn.post_solve_actions();
Expand Down
168 changes: 13 additions & 155 deletions unit_tests/multiphase/test_vof_tools.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -201,7 +201,7 @@ amrex::Real interface_band_test_impl(amr_wind::Field& vof)
return error_total;
}

amrex::Real sharpen_test_impl(amr_wind::Field& vof)
amrex::Real initvof_test_impl(amr_wind::Field& vof)
{
amrex::Real error_total = 0;

Expand All @@ -217,29 +217,18 @@ amrex::Real sharpen_test_impl(amr_wind::Field& vof)

amrex::Loop(bx, [=, &error](int i, int j, int k) noexcept {
// Initial VOF distribution is 0, 0.3, 0.7, or 1.0
// Expected answer based on implementation of
// amr_wind::multiphase::sharpen_kernel
amrex::Real vof_answer = 0.0;
if (i + j + k > 5) {
// because VOF < 0.5
const amrex::Real sign = -1.0;
const amrex::Real delta = std::abs(0.3 - 0.5);
vof_answer = 0.5 + sign * std::pow(delta, 1.0 / 3.0);
vof_answer = 0.3;
}
if (i + j + k > 10) {
// because VOF > 0.5
const amrex::Real sign = 1.0;
const amrex::Real delta = std::abs(0.7 - 0.5);
vof_answer = 0.5 + sign * std::pow(delta, 1.0 / 3.0);
vof_answer = 0.7;
}
// Set up a liquid cell
if (i == 0 && j == 0 && k == 0) {
vof_answer = 1.0;
}

// Limit answer to VOF bounds
vof_answer = std::max(0.0, std::min(1.0, vof_answer));

// Difference between actual and expected
error += std::abs(vof_arr(i, j, k) - vof_answer);
});
Expand All @@ -250,58 +239,6 @@ amrex::Real sharpen_test_impl(amr_wind::Field& vof)
return error_total;
}

amrex::Real sharpen_test_density_impl(
amr_wind::Field& dens, const amrex::Real rho1, const amrex::Real rho2)
{
amrex::Real error_total = 0;

for (int lev = 0; lev < dens.repo().num_active_levels(); ++lev) {

error_total += amrex::ReduceSum(
dens(lev), 0,
[=] AMREX_GPU_HOST_DEVICE(
amrex::Box const& bx,
amrex::Array4<amrex::Real const> const& dens_arr)
-> amrex::Real {
amrex::Real error = 0;

amrex::Loop(bx, [=, &error](int i, int j, int k) noexcept {
// Initial VOF distribution is 0, 0.3, 0.7, or 1.0
// Expected answer based on implementation of
// amr_wind::multiphase::sharpen_kernel
amrex::Real vof_answer = 0.0;
if (i + j + k > 5) {
// because VOF < 0.5
const amrex::Real sign = -1.0;
const amrex::Real delta = std::abs(0.3 - 0.5);
vof_answer = 0.5 + sign * std::pow(delta, 1.0 / 3.0);
}
if (i + j + k > 10) {
// because VOF > 0.5
const amrex::Real sign = 1.0;
const amrex::Real delta = std::abs(0.7 - 0.5);
vof_answer = 0.5 + sign * std::pow(delta, 1.0 / 3.0);
}
// Set up a liquid cell
if (i == 0 && j == 0 && k == 0) {
vof_answer = 1.0;
}

// Limit answer to VOF bounds
vof_answer = std::max(0.0, std::min(1.0, vof_answer));

// Difference between actual and expected
const amrex::Real dens_answer =
rho1 * vof_answer + rho2 * (1.0 - vof_answer);
error += std::abs(dens_arr(i, j, k) - dens_answer);
});

return error;
});
}
return error_total;
}

} // namespace

TEST_F(VOFToolTest, interface_band)
Expand Down Expand Up @@ -362,7 +299,7 @@ TEST_F(VOFToolTest, levelset_to_vof)
EXPECT_NEAR(error_total, 0.0, 0.016);
}

TEST_F(VOFToolTest, sharpen_acquired_vof)
TEST_F(VOFToolTest, replace_masked_vof)
{

populate_parameters();
Expand All @@ -372,100 +309,21 @@ TEST_F(VOFToolTest, sharpen_acquired_vof)
const int ncomp = 1;
const int nghost = 3;
const int nghost_int = 1;
auto& vof = repo.declare_field("vof", ncomp, nghost);
auto& vof_new = repo.declare_field("vof", ncomp, nghost);
auto& vof_mod = repo.declare_field("vof_mod", ncomp, nghost);
auto& iblank = repo.declare_int_field("iblank_cell", ncomp, nghost_int);

// Use as if entire domain is from nalu
iblank.setVal(-1);
// Initialize and sharpen vof
init_vof(vof);
amr_wind::multiphase::sharpen_acquired_vof(1, iblank, vof);

// Check results
amrex::Real error_total = sharpen_test_impl(vof);
amrex::ParallelDescriptor::ReduceRealSum(error_total);
EXPECT_NEAR(error_total, 0.0, 1e-15);
}

TEST_F(VOFToolTest, sharpen_replace_old_density)
{

const amrex::Real rho1 = 1000.0;
const amrex::Real rho2 = 1.0;
populate_parameters();
{
amrex::ParmParse pp("geometry");
amrex::Vector<int> periodic{{1, 1, 1}};
pp.addarr("is_periodic", periodic);
}
{
amrex::ParmParse pp("MultiPhase");
pp.add("density_fluid1", rho1);
pp.add("density_fluid2", rho2);
}
{
amrex::ParmParse pp("incflo");
amrex::Vector<std::string> physics{"MultiPhase"};
pp.addarr("physics", physics);
pp.add("use_godunov", (int)1);
}

initialize_mesh();

auto& repo = sim().repo();

// PDE manager, for access to VOF PDE later
auto& pde_mgr = sim().pde_manager();
// Setup of icns provides MAC velocities
pde_mgr.register_icns();

// Initialize physics for the sake of MultiPhase routines
sim().init_physics();

// Initialize volume fraction field
auto& vof = repo.get_field("vof");

// Set up iblank field
const int ncomp = 1;
const int nghost_int = 1;
auto& iblank = repo.declare_int_field("iblank_cell", ncomp, nghost_int);
// Use as if entire domain is from nalu
iblank.setVal(-1);

// Initialize vof
init_vof(vof);
init_vof(vof_new);
// Initialize other field (working copy of vof)
vof_mod.setVal(100.0);
// Replace masked vof values with new ones
amr_wind::multiphase::replace_masked_vof(1, iblank, vof_mod, vof_new);

// Initialize advective velocities with something
auto& umac = repo.get_field("u_mac");
auto& vmac = repo.get_field("v_mac");
auto& wmac = repo.get_field("w_mac");
umac.setVal(1.0);
vmac.setVal(1.0);
wmac.setVal(1.0);

// Get vof equation handle and perform init
auto& seqn = pde_mgr(
amr_wind::pde::VOF::pde_name() + "-" +
amr_wind::fvm::Godunov::scheme_name());
seqn.initialize();

// Do advection step (which results in sharpening and in copying of
// sharpened vof to old vof)
seqn.compute_advection_term(amr_wind::FieldState::Old);

// Check that old vof is equal to the sharpened solution
amrex::Real error_total =
sharpen_test_impl(vof.state(amr_wind::FieldState::Old));
amrex::ParallelDescriptor::ReduceRealSum(error_total);
EXPECT_NEAR(error_total, 0.0, 1e-15);

// Do post solve (which recalculates old density with old vof)
seqn.post_solve_actions();

// Check old density
auto& dens = repo.get_field("density");
error_total = sharpen_test_density_impl(
dens.state(amr_wind::FieldState::Old), rho1, rho2);
// Check results
amrex::Real error_total = initvof_test_impl(vof_mod);
amrex::ParallelDescriptor::ReduceRealSum(error_total);
EXPECT_NEAR(error_total, 0.0, 1e-15);
}
Expand Down

0 comments on commit eb2df80

Please sign in to comment.