From e8b5656e44a5633bcad93ad38e57239c4dced207 Mon Sep 17 00:00:00 2001 From: Michael B Kuhn <31661049+mbkuhn@users.noreply.github.com> Date: Fri, 2 Aug 2024 10:15:32 -0600 Subject: [PATCH] Multiphase overset -- correcting some things (#1175) * 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 --- amr-wind/overset/OversetOps.cpp | 9 ++- amr-wind/overset/overset_ops_K.H | 8 +-- amr-wind/overset/overset_ops_routines.H | 1 + amr-wind/overset/overset_ops_routines.cpp | 13 ++-- .../multiphase/test_vof_overset_ops.cpp | 69 ++++++++++++------- 5 files changed, 62 insertions(+), 38 deletions(-) diff --git a/amr-wind/overset/OversetOps.cpp b/amr-wind/overset/OversetOps.cpp index 920c45f8e7..ca20436d5e 100644 --- a/amr-wind/overset/OversetOps.cpp +++ b/amr-wind/overset/OversetOps.cpp @@ -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); } } @@ -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); } diff --git a/amr-wind/overset/overset_ops_K.H b/amr-wind/overset/overset_ops_K.H index bd8780d3f0..3f5fd3ae68 100644 --- a/amr-wind/overset/overset_ops_K.H +++ b/amr-wind/overset/overset_ops_K.H @@ -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); diff --git a/amr-wind/overset/overset_ops_routines.H b/amr-wind/overset/overset_ops_routines.H index de9afbdff8..b14043b26a 100644 --- a/amr-wind/overset/overset_ops_routines.H +++ b/amr-wind/overset/overset_ops_routines.H @@ -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& dx, const amrex::Real tol); // Apply reinitialization fluxes to modify fields diff --git a/amr-wind/overset/overset_ops_routines.cpp b/amr-wind/overset/overset_ops_routines.cpp index 57becc49ec..f9ae2ea1ed 100644 --- a/amr-wind/overset/overset_ops_routines.cpp +++ b/amr-wind/overset/overset_ops_routines.cpp @@ -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& dx, const amrex::Real tol) { // Get the maximum flux magnitude, but just for vof fluxes @@ -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); }); @@ -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); }); @@ -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); }); diff --git a/unit_tests/multiphase/test_vof_overset_ops.cpp b/unit_tests/multiphase/test_vof_overset_ops.cpp index fa27e9d538..876fb4c576 100644 --- a/unit_tests/multiphase/test_vof_overset_ops.cpp +++ b/unit_tests/multiphase/test_vof_overset_ops.cpp @@ -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); } }); @@ -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); @@ -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); @@ -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); @@ -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)