Skip to content

Commit

Permalink
Multiphase overset -- correcting some things (#1175)
Browse files Browse the repository at this point in the history 
* changes to pseudo_dt to use dx correctly
* measure convergence nondimensionally
* fix known norm usage sign error
* updating unit test to go with other changes
  • Loading branch information
@mbkuhn
mbkuhn authored Aug 2, 2024
1 parent 0c610ba commit e8b5656
Show file tree
Hide file tree
Showing 5 changed files with 62 additions and 38 deletions.
9 changes: 6 additions & 3 deletions amr-wind/overset/OversetOps.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -320,8 +320,10 @@ void OversetOps::sharpen_nalu_data()
// Measure convergence to determine if loop can stop
if (calc_convg) {
// Update error at specified interval of steps
const amrex::Real err_lev = overset_ops::measure_convergence(
(*flux_x)(lev), (*flux_y)(lev), (*flux_z)(lev));
const amrex::Real err_lev =
overset_ops::measure_convergence(
(*flux_x)(lev), (*flux_y)(lev), (*flux_z)(lev)) /
pvscale;
err = amrex::max(err, err_lev);
}
}
Expand All @@ -336,10 +338,11 @@ void OversetOps::sharpen_nalu_data()

// Get pseudo dt (dtau)
for (int lev = 0; lev < nlevels; ++lev) {
const auto dx = (geom[lev]).CellSizeArray();
// Compare vof fluxes to vof in source cells
// Convergence tolerance determines what size of fluxes matter
const amrex::Real ptfac_lev = overset_ops::calculate_pseudo_dt_flux(
(*flux_x)(lev), (*flux_y)(lev), (*flux_z)(lev), vof(lev),
(*flux_x)(lev), (*flux_y)(lev), (*flux_z)(lev), vof(lev), dx,
m_convg_tol);
ptfac = amrex::min(ptfac, ptfac_lev);
}
Expand Down
8 changes: 4 additions & 4 deletions amr-wind/overset/overset_ops_K.H
View file
Original file line number Diff line number Diff line change
Expand Up @@ -32,12 +32,12 @@ amrex::Real AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE alpha_flux(
// Gradient of phi normal to interface
const amrex::Real gphi = (vof(iv) - vof(ivm));
// Normal vector in each cell (already normalized)
const amrex::Real norm_ = normal(iv, dir);
const amrex::Real norm_nb = normal(ivm, dir);
const amrex::Real norm_ = std::abs(normal(iv, dir));
const amrex::Real norm_nb = std::abs(normal(ivm, dir));

// Determine which delta_phi (and multiply by normal)
// The sign depends on side of flux face (like upwinding)
const amrex::Real dphi_ = (tg_vof(iv) - vof(iv)) * (-norm_);
// The sign depends on side of flux face
const amrex::Real dphi_ = (vof(iv) - tg_vof(iv)) * norm_;
const amrex::Real dphi_nb = (tg_vof(ivm) - vof(ivm)) * norm_nb;
// Average value used across the interface
amrex::Real dphi_eval = 0.5 * (dphi_ + dphi_nb);
Expand Down
1 change: 1 addition & 0 deletions amr-wind/overset/overset_ops_routines.H
View file
Original file line number Diff line number Diff line change
Expand Up @@ -68,6 +68,7 @@ amrex::Real calculate_pseudo_dt_flux(
const amrex::MultiFab& mf_fy,
const amrex::MultiFab& mf_fz,
const amrex::MultiFab& mf_vof,
const amrex::GpuArray<amrex::Real, AMREX_SPACEDIM>& dx,
const amrex::Real tol);

// Apply reinitialization fluxes to modify fields
Expand Down
13 changes: 7 additions & 6 deletions amr-wind/overset/overset_ops_routines.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -252,6 +252,7 @@ amrex::Real calculate_pseudo_dt_flux(
const amrex::MultiFab& mf_fy,
const amrex::MultiFab& mf_fz,
const amrex::MultiFab& mf_vof,
const amrex::GpuArray<amrex::Real, AMREX_SPACEDIM>& dx,
const amrex::Real tol)
{
// Get the maximum flux magnitude, but just for vof fluxes
Expand All @@ -265,10 +266,10 @@ amrex::Real calculate_pseudo_dt_flux(
amrex::Real pdt_lim = 1.0;
if (fx(i, j, k, 0) > tol && vof(i, j, k) > tol) {
// VOF is removed from cell i
pdt_lim = vof(i, j, k) / fx(i, j, k, 0);
pdt_lim = vof(i, j, k) * dx[0] / fx(i, j, k, 0);
} else if (fx(i, j, k, 0) < -tol && vof(i - 1, j, k) > tol) {
// VOF is removed from cell i-1
pdt_lim = vof(i - 1, j, k) / -fx(i, j, k, 0);
pdt_lim = vof(i - 1, j, k) * dx[0] / -fx(i, j, k, 0);
}
pdt_fab = amrex::min(pdt_fab, pdt_lim);
});
Expand All @@ -284,10 +285,10 @@ amrex::Real calculate_pseudo_dt_flux(
amrex::Real pdt_lim = 1.0;
if (fy(i, j, k, 0) > tol && vof(i, j, k) > tol) {
// VOF is removed from cell j
pdt_lim = vof(i, j, k) / fy(i, j, k, 0);
pdt_lim = vof(i, j, k) * dx[1] / fy(i, j, k, 0);
} else if (fy(i, j, k, 0) < -tol && vof(i, j - 1, k) > tol) {
// VOF is removed from cell j-1
pdt_lim = vof(i, j - 1, k) / -fy(i, j, k, 0);
pdt_lim = vof(i, j - 1, k) * dx[1] / -fy(i, j, k, 0);
}
pdt_fab = amrex::min(pdt_fab, pdt_lim);
});
Expand All @@ -303,10 +304,10 @@ amrex::Real calculate_pseudo_dt_flux(
amrex::Real pdt_lim = 1.0;
if (fz(i, j, k, 0) > tol && vof(i, j, k) > tol) {
// VOF is removed from cell k
pdt_lim = vof(i, j, k) / fz(i, j, k, 0);
pdt_lim = vof(i, j, k) * dx[2] / fz(i, j, k, 0);
} else if (fz(i, j, k, 0) < -tol && vof(i, j, k - 1) > tol) {
// VOF is removed from cell k-1
pdt_lim = vof(i, j, k - 1) / -fz(i, j, k, 0);
pdt_lim = vof(i, j, k - 1) * dx[2] / -fz(i, j, k, 0);
}
pdt_fab = amrex::min(pdt_fab, pdt_lim);
});
Expand Down
69 changes: 44 additions & 25 deletions unit_tests/multiphase/test_vof_overset_ops.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -399,41 +399,41 @@ amrex::Real check_alpha_flux_impl(amr_wind::Field& flux, const int& dir)
if (idx == 0) {
// Both within margin, will average
const amrex::Real flux_answer =
0.5 * ((0.5 - 0.41) * -1.0 + (0.5 - 0.55) * -1.0);
0.5 * ((0.5 - 0.41) * 1.0 + (0.5 - 0.55) * -1.0);
error += std::abs(f_arr(i, j, k) - flux_answer);
} else if (idx == 1) {
// Current is low, neighbor within margin
const amrex::Real flux_answer = (0.2 - 0.1) * 1.0;
const amrex::Real flux_answer = (0.2 - 0.1) * -1.0;
error += std::abs(f_arr(i, j, k) - flux_answer);
} else if (idx == 2) {
// Low on both sides, gphi > 0
const amrex::Real flux_answer = (0.2 - 0.1) * -1.0;
const amrex::Real flux_answer = (0.2 - 0.1) * 1.0;
error += std::abs(f_arr(i, j, k) - flux_answer);
} else if (idx == 3) {
// Opposite sides of margin, no conditional fits
const amrex::Real flux_answer =
0.5 * ((0.8 - 0.7) * 1.0 + (0.22 - 0.2) * 1.0);
0.5 * ((0.8 - 0.7) * -1.0 + (0.22 - 0.2) * 1.0);
error += std::abs(f_arr(i, j, k) - flux_answer);
} else if (idx == 4) {
// High on both sides, gphi < 0
const amrex::Real flux_answer = (0.8 - 0.7) * -1.0;
const amrex::Real flux_answer = (0.8 - 0.7) * 1.0;
error += std::abs(f_arr(i, j, k) - flux_answer);
} else if (idx == 5) {
// High on both sides, gphi > 0
const amrex::Real flux_answer = (0.75 - 0.8) * 1.0;
const amrex::Real flux_answer = (0.75 - 0.8) * -1.0;
error += std::abs(f_arr(i, j, k) - flux_answer);
} else if (idx == 6) {
// Current is within margin, neighbor is high
const amrex::Real flux_answer = (0.75 - 0.8) * -1.0;
const amrex::Real flux_answer = (0.75 - 0.8) * 1.0;
error += std::abs(f_arr(i, j, k) - flux_answer);
} else if (idx == 7) {
// Current is low, neighbor is within margin
const amrex::Real flux_answer = (0.2 - 0.3) * 1.0;
const amrex::Real flux_answer = (0.2 - 0.3) * -1.0;
error += std::abs(f_arr(i, j, k) - flux_answer);
} else if (idx == 8) {
// Opposite sides of margin, no conditional fits
const amrex::Real flux_answer =
0.5 * ((0.6 - 0.7) * -1.0 + (0.2 - 0.3) * -1.0);
0.5 * ((0.6 - 0.7) * -1.0 + (0.2 - 0.3) * 1.0);
error += std::abs(f_arr(i, j, k) - flux_answer);
}
});
Expand Down Expand Up @@ -756,13 +756,27 @@ TEST_F(VOFOversetOps, pseudo_vscale_dt)
iblank.setVal(-1);
constexpr amrex::Real margin = 0.1;
constexpr amrex::Real convg_tol = 1e-8;
// With simple fluxes, pseudo dt should be 100%
constexpr amrex::Real pdt_answer = 1.0;
// With vof and target_vof arrays, max vof removed is 50%, doubling pdt
constexpr amrex::Real pdt_answer = 2.0;
// With a single level, pseudo velocity scale should be dx of lev 0
const auto dx_lev0 = repo.mesh().Geom(0).CellSizeArray();
const amrex::Real pvs_answer =
std::min(std::min(dx_lev0[0], dx_lev0[1]), dx_lev0[2]);

// Pseudo-velocity scale, should be the smallest dx in iblank region
const auto& iblank_cell = repo.get_int_field("iblank_cell");
const amrex::Real max_pvscale = 100.;
amrex::Real pvscale = max_pvscale;
for (int lev = 0; lev < repo.num_active_levels(); ++lev) {
const amrex::Real pvscale_lev =
amr_wind::overset_ops::calculate_pseudo_velocity_scale(
iblank_cell(lev), repo.mesh().Geom(lev).CellSizeArray(),
max_pvscale);
pvscale = std::min(pvscale, pvscale_lev);
}
amrex::ParallelDescriptor::ReduceRealMin(pvscale);
EXPECT_DOUBLE_EQ(pvscale, pvs_answer);

// Create flux fields
auto& flux_x =
repo.declare_field("flux_x", 1, 0, 1, amr_wind::FieldLoc::XFACE);
Expand All @@ -782,7 +796,9 @@ TEST_F(VOFOversetOps, pseudo_vscale_dt)
for (int lev = 0; lev < repo.num_active_levels(); ++lev) {
const amrex::Real ptfac_lev =
amr_wind::overset_ops::calculate_pseudo_dt_flux(
flux_x(lev), flux_y(lev), flux_z(lev), vof(lev), convg_tol);
flux_x(lev), flux_y(lev), flux_z(lev), vof(lev), dx_lev0,
convg_tol) /
pvscale;
ptfac = amrex::min(ptfac, ptfac_lev);
}
amrex::ParallelDescriptor::ReduceRealMin(ptfac);
Expand All @@ -801,7 +817,9 @@ TEST_F(VOFOversetOps, pseudo_vscale_dt)
for (int lev = 0; lev < repo.num_active_levels(); ++lev) {
const amrex::Real ptfac_lev =
amr_wind::overset_ops::calculate_pseudo_dt_flux(
flux_x(lev), flux_y(lev), flux_z(lev), vof(lev), convg_tol);
flux_x(lev), flux_y(lev), flux_z(lev), vof(lev), dx_lev0,
convg_tol) /
pvscale;
ptfac = amrex::min(ptfac, ptfac_lev);
}
amrex::ParallelDescriptor::ReduceRealMin(ptfac);
Expand All @@ -820,25 +838,26 @@ TEST_F(VOFOversetOps, pseudo_vscale_dt)
for (int lev = 0; lev < repo.num_active_levels(); ++lev) {
const amrex::Real ptfac_lev =
amr_wind::overset_ops::calculate_pseudo_dt_flux(
flux_x(lev), flux_y(lev), flux_z(lev), vof(lev), convg_tol);
flux_x(lev), flux_y(lev), flux_z(lev), vof(lev), dx_lev0,
convg_tol) /
pvscale;
ptfac = amrex::min(ptfac, ptfac_lev);
}
amrex::ParallelDescriptor::ReduceRealMin(ptfac);
EXPECT_DOUBLE_EQ(ptfac, pdt_answer);

// Pseudo-velocity scale, should be the smallest dx in iblank region
const auto& iblank_cell = repo.get_int_field("iblank_cell");
const amrex::Real max_pvscale = 100.;
amrex::Real pvscale = max_pvscale;
// Redo z direction with max being 1, as in overall algorithm
ptfac = 1.0;
for (int lev = 0; lev < repo.num_active_levels(); ++lev) {
const amrex::Real pvscale_lev =
amr_wind::overset_ops::calculate_pseudo_velocity_scale(
iblank_cell(lev), repo.mesh().Geom(lev).CellSizeArray(),
max_pvscale);
pvscale = std::min(pvscale, pvscale_lev);
const amrex::Real ptfac_lev =
amr_wind::overset_ops::calculate_pseudo_dt_flux(
flux_x(lev), flux_y(lev), flux_z(lev), vof(lev), dx_lev0,
convg_tol) /
pvscale;
ptfac = amrex::min(ptfac, ptfac_lev);
}
amrex::ParallelDescriptor::ReduceRealMin(pvscale);
EXPECT_DOUBLE_EQ(pvscale, pvs_answer);
amrex::ParallelDescriptor::ReduceRealMin(ptfac);
EXPECT_DOUBLE_EQ(ptfac, 1.0);
}

TEST_F(VOFOversetOps, psource_manual)
Expand Down

0 comments on commit e8b5656

Please sign in to comment.