diff --git a/amr-wind/equation_systems/icns/icns_advection.H b/amr-wind/equation_systems/icns/icns_advection.H index 7e9eb9d5dc..78e41e645b 100644 --- a/amr-wind/equation_systems/icns/icns_advection.H +++ b/amr-wind/equation_systems/icns/icns_advection.H @@ -61,6 +61,7 @@ private: #endif MLMGOptions m_options; bool m_has_overset{false}; + bool m_disable_mphase_overset_ops{false}; bool m_need_init{true}; bool m_variable_density{false}; bool m_mesh_mapping{false}; diff --git a/amr-wind/equation_systems/icns/icns_advection.cpp b/amr-wind/equation_systems/icns/icns_advection.cpp index 597dbb767b..5b312662d6 100644 --- a/amr-wind/equation_systems/icns/icns_advection.cpp +++ b/amr-wind/equation_systems/icns/icns_advection.cpp @@ -62,6 +62,11 @@ MacProjOp::MacProjOp( pp.query("use_fft", m_use_fft); } #endif + + // Always off for single-phase flows, this can turn off multiphase overset + // operations for all cases + amrex::ParmParse pp_o("Overset"); + pp_o.query("disable_mphase_ops", m_disable_mphase_overset_ops); } void MacProjOp::enforce_inout_solvability( @@ -78,7 +83,7 @@ void MacProjOp::enforce_inout_solvability( void MacProjOp::init_projector(const MacProjOp::FaceFabPtrVec& beta) noexcept { // Prepare masking for projection - if (m_has_overset) { + if (m_has_overset && !m_disable_mphase_overset_ops) { amr_wind::overset_ops::prepare_mask_cell_for_mac(m_repo); } @@ -107,7 +112,7 @@ void MacProjOp::init_projector(const MacProjOp::FaceFabPtrVec& beta) noexcept void MacProjOp::init_projector(const amrex::Real beta) noexcept { // Prepare masking for projection - if (m_has_overset) { + if (m_has_overset && !m_disable_mphase_overset_ops) { amr_wind::overset_ops::prepare_mask_cell_for_mac(m_repo); } @@ -365,7 +370,7 @@ void MacProjOp::operator()(const FieldState fstate, const amrex::Real dt) } // Prepare masking for remainder of algorithm - if (m_has_overset) { + if (m_has_overset && !m_disable_mphase_overset_ops) { amr_wind::overset_ops::revert_mask_cell_after_mac(m_repo); } diff --git a/amr-wind/equation_systems/vof/volume_fractions.H b/amr-wind/equation_systems/vof/volume_fractions.H index 7582ce3489..a1535815da 100644 --- a/amr-wind/equation_systems/vof/volume_fractions.H +++ b/amr-wind/equation_systems/vof/volume_fractions.H @@ -5,6 +5,16 @@ #include #include "amr-wind/utilities/constants.H" +namespace { + +AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE int +idx(const int i_off, const int j_off, const int k_off) +{ + return (k_off + 1) + 3 * (j_off + 1) + 9 * (i_off + 1); +} + +} // namespace + namespace amr_wind::multiphase { /** Young's finite-difference gradient scheme. @@ -70,56 +80,106 @@ youngs_finite_difference_normal_neumann( const int i, const int j, const int k, - const int ibdy, - const int jbdy, - const int kbdy, + amrex::Array4 const& iblank, amrex::Array4 const& volfrac, amrex::Real& mx, amrex::Real& my, amrex::Real& mz) noexcept { - // Do neumann condition via indices - const int im1 = ibdy == -1 ? i : i - 1; - const int jm1 = jbdy == -1 ? j : j - 1; - const int km1 = kbdy == -1 ? k : k - 1; - const int ip1 = ibdy == +1 ? i : i + 1; - const int jp1 = jbdy == +1 ? j : j + 1; - const int kp1 = kbdy == +1 ? k : k + 1; - - amrex::Real mm1 = volfrac(im1, jm1, km1) + volfrac(im1, jm1, kp1) + - volfrac(im1, jp1, km1) + volfrac(im1, jp1, kp1) + - 2.0 * (volfrac(im1, jm1, k) + volfrac(im1, jp1, k) + - volfrac(im1, j, km1) + volfrac(im1, j, kp1)) + - 4.0 * volfrac(im1, j, k); - amrex::Real mm2 = volfrac(ip1, jm1, km1) + volfrac(ip1, jm1, kp1) + - volfrac(ip1, jp1, km1) + volfrac(ip1, jp1, kp1) + - 2.0 * (volfrac(ip1, jm1, k) + volfrac(ip1, jp1, k) + - volfrac(ip1, j, km1) + volfrac(ip1, j, kp1)) + - 4.0 * volfrac(ip1, j, k); + // Do neumann condition using iblank array + // Cells that are not iblank == -1 get calculated with extrapolation using + // the average of all valid (relying only on iblank == -1 cells) gradients + const int ibl_c = iblank(i, j, k); + amrex::GpuArray dvf{0., 0., 0., 0., 0., 0.}; + const amrex::IntVect iv{i, j, k}; + for (int ddir = 0; ddir < 3; ++ddir) { + int ct_0 = 0; + int ct_1 = 0; + const int odir0 = (ddir + 1) % 3; + const int odir1 = (ddir + 2) % 3; + for (int n0 = -1; n0 <= 1; ++n0) { + for (int n1 = -1; n1 <= 1; ++n1) { + amrex::IntVect iv_c = iv; + iv_c[odir0] += n0; + iv_c[odir1] += n1; + amrex::IntVect iv_p = iv_c; + amrex::IntVect iv_m = iv_c; + iv_p[ddir] += 1; + iv_m[ddir] += -1; + if (iblank(iv_m) == iblank(iv_c) && iblank(iv_m) == ibl_c) { + dvf[ddir] += volfrac(iv_c) - volfrac(iv_m); + ct_0++; + } + if (iblank(iv_p) == iblank(iv_c) && iblank(iv_p) == ibl_c) { + dvf[ddir + 3] += volfrac(iv_p) - volfrac(iv_c); + ct_1++; + } + } + } + dvf[ddir] /= (amrex::Real)ct_0 + constants::EPS; + dvf[ddir + 3] /= (amrex::Real)ct_1 + constants::EPS; + } + + amrex::GpuArray vf; + vf.fill(0.); + for (int n0 = -1; n0 <= 1; ++n0) { + for (int n1 = -1; n1 <= 1; ++n1) { + for (int n2 = -1; n2 <= 1; ++n2) { + amrex::IntVect iv_l = iv; + iv_l[0] += n0; + iv_l[1] += n1; + iv_l[2] += n2; + vf[idx(n0, n1, n2)] = + iblank(iv_l) == ibl_c + ? volfrac(iv_l) + : std::max( + 0., + std::min( + 1., volfrac(iv) + + (amrex::Real)n0 * + dvf[0 + 3 * std::max(0, n0)] + + (amrex::Real)n1 * + dvf[1 + 3 * std::max(0, n1)] + + (amrex::Real)n2 * + dvf[2 + 3 * std::max(0, n2)])); + } + } + } + + amrex::Real mm1 = vf[idx(-1, -1, -1)] + vf[idx(-1, -1, 1)] + + vf[idx(-1, 1, -1)] + vf[idx(-1, 1, 1)] + + 2.0 * (vf[idx(-1, -1, 0)] + vf[idx(-1, 1, 0)] + + vf[idx(-1, 0, -1)] + vf[idx(-1, 0, 1)]) + + 4.0 * vf[idx(-1, 0, 0)]; + amrex::Real mm2 = vf[idx(1, -1, -1)] + vf[idx(1, -1, 1)] + + vf[idx(1, 1, -1)] + vf[idx(1, 1, 1)] + + 2.0 * (vf[idx(1, -1, 0)] + vf[idx(1, 1, 0)] + + vf[idx(1, 0, -1)] + vf[idx(1, 0, 1)]) + + 4.0 * vf[idx(1, 0, 0)]; mx = mm1 - mm2; - mm1 = volfrac(im1, jm1, km1) + volfrac(im1, jm1, kp1) + - volfrac(ip1, jm1, km1) + volfrac(ip1, jm1, kp1) + - 2.0 * (volfrac(im1, jm1, k) + volfrac(ip1, jm1, k) + - volfrac(i, jm1, km1) + volfrac(i, jm1, kp1)) + - 4.0 * volfrac(i, jm1, k); - mm2 = volfrac(im1, jp1, km1) + volfrac(im1, jp1, kp1) + - volfrac(ip1, jp1, km1) + volfrac(ip1, jp1, kp1) + - 2.0 * (volfrac(im1, jp1, k) + volfrac(ip1, jp1, k) + - volfrac(i, jp1, km1) + volfrac(i, jp1, kp1)) + - 4.0 * volfrac(i, jp1, k); + mm1 = vf[idx(-1, -1, -1)] + vf[idx(-1, -1, 1)] + vf[idx(1, -1, -1)] + + vf[idx(1, -1, 1)] + + 2.0 * (vf[idx(-1, -1, 0)] + vf[idx(1, -1, 0)] + vf[idx(0, -1, -1)] + + vf[idx(0, -1, 1)]) + + 4.0 * vf[idx(0, -1, 0)]; + mm2 = vf[idx(-1, 1, -1)] + vf[idx(-1, 1, 1)] + vf[idx(1, 1, -1)] + + vf[idx(1, 1, 1)] + + 2.0 * (vf[idx(-1, 1, 0)] + vf[idx(1, 1, 0)] + vf[idx(0, 1, -1)] + + vf[idx(0, 1, 1)]) + + 4.0 * vf[idx(0, 1, 0)]; my = mm1 - mm2; - mm1 = volfrac(im1, jm1, km1) + volfrac(im1, jp1, km1) + - volfrac(ip1, jm1, km1) + volfrac(ip1, jp1, km1) + - 2.0 * (volfrac(im1, j, km1) + volfrac(ip1, j, km1) + - volfrac(i, jm1, km1) + volfrac(i, jp1, km1)) + - 4.0 * volfrac(i, j, km1); - mm2 = volfrac(im1, jm1, kp1) + volfrac(im1, jp1, kp1) + - volfrac(ip1, jm1, kp1) + volfrac(ip1, jp1, kp1) + - 2.0 * (volfrac(im1, j, kp1) + volfrac(ip1, j, kp1) + - volfrac(i, jm1, kp1) + volfrac(i, jp1, kp1)) + - 4.0 * volfrac(i, j, kp1); + mm1 = vf[idx(-1, -1, -1)] + vf[idx(-1, 1, -1)] + vf[idx(1, -1, -1)] + + vf[idx(1, 1, -1)] + + 2.0 * (vf[idx(-1, 0, -1)] + vf[idx(1, 0, -1)] + vf[idx(0, -1, -1)] + + vf[idx(0, 1, -1)]) + + 4.0 * vf[idx(0, 0, -1)]; + mm2 = vf[idx(-1, -1, 1)] + vf[idx(-1, 1, 1)] + vf[idx(1, -1, 1)] + + vf[idx(1, 1, 1)] + + 2.0 * (vf[idx(-1, 0, 1)] + vf[idx(1, 0, 1)] + vf[idx(0, -1, 1)] + + vf[idx(0, 1, 1)]) + + 4.0 * vf[idx(0, 0, 1)]; mz = mm1 - mm2; } diff --git a/amr-wind/overset/OversetOps.H b/amr-wind/overset/OversetOps.H index 81df229a0d..81c2a4481c 100644 --- a/amr-wind/overset/OversetOps.H +++ b/amr-wind/overset/OversetOps.H @@ -27,6 +27,10 @@ private: // Check for multiphase sim bool m_vof_exists{false}; + // Option to turn off all multiphase-specific treatments (for testing and + // comparison) + bool m_disable_mphase_ops{false}; + // Coupling options bool m_replace_gradp_postsolve{true}; @@ -40,6 +44,8 @@ private: amrex::Real m_relative_length_scale = 1.5; amrex::Real m_upw_margin = 0.1; amrex::Real m_target_cutoff = 0.0; // proc_tgvof_tol + // Maximum pseudoCFL + amrex::Real m_pCFL = 1.0; // Tolerance for VOF-related bound checks const amrex::Real m_vof_tol = 1e-12; diff --git a/amr-wind/overset/OversetOps.cpp b/amr-wind/overset/OversetOps.cpp index 1f8d510bfa..1515f27fab 100644 --- a/amr-wind/overset/OversetOps.cpp +++ b/amr-wind/overset/OversetOps.cpp @@ -22,13 +22,17 @@ void OversetOps::initialize(CFDSim& sim) pp.query("reinit_rlscale", m_relative_length_scale); pp.query("reinit_upw_margin", m_upw_margin); pp.query("reinit_target_cutoff", m_target_cutoff); + pp.query("reinit_pseudo_CFL", m_pCFL); // Queries for coupling options pp.query("replace_gradp_postsolve", m_replace_gradp_postsolve); pp.query("verbose", m_verbose); + // Option to turn it all off + pp.query("disable_mphase_ops", m_disable_mphase_ops); + m_vof_exists = m_sim_ptr->repo().field_exists("vof"); - if (m_vof_exists) { + if (m_vof_exists && !m_disable_mphase_ops) { m_mphase = &m_sim_ptr->physics_manager().get(); // Check combination of pressure options @@ -59,7 +63,7 @@ void OversetOps::pre_advance_work() // Update pressure gradient using updated overset pressure field update_gradp(); - if (m_vof_exists) { + if (m_vof_exists && !m_disable_mphase_ops) { // Reinitialize fields sharpen_nalu_data(); if (m_mphase->perturb_pressure()) { @@ -148,7 +152,7 @@ void OversetOps::update_gradp() void OversetOps::post_advance_work() { // Replace and reapply pressure gradient if requested - if (m_vof_exists && m_replace_gradp_postsolve) { + if (m_vof_exists && m_replace_gradp_postsolve && !m_disable_mphase_ops) { replace_masked_gradp(); } } @@ -161,29 +165,37 @@ void OversetOps::parameter_output() const { // Print the details if (m_verbose > 0 && m_vof_exists) { - // Important parameters - amrex::Print() << "\nOverset Coupling Parameters: \n" - << "---- Replace overset pres grad: " - << m_replace_gradp_postsolve << std::endl; - 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 - << "---- Convergence tolerance: " << m_convg_tol - << std::endl - << "---- Relative length scale: " - << m_relative_length_scale << std::endl - << "---- Upwinding VOF margin : " << m_upw_margin - << std::endl; - if (m_verbose > 1) { - // Less important or less used parameters - amrex::Print() << "---- Calc. conv. interval : " - << m_calc_convg_interval << std::endl - << "---- Target field cutoff : " << m_target_cutoff + if (m_disable_mphase_ops) { + amrex::Print() + << "\nOverset Coupling: multiphase operations are disabled\n"; + } else { + // Important parameters + amrex::Print() << "\nOverset Coupling Parameters: \n" + << "---- Replace overset pres grad: " + << m_replace_gradp_postsolve << std::endl; + 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 + << "---- Convergence tolerance: " << m_convg_tol + << std::endl + << "---- Relative length scale: " + << m_relative_length_scale << std::endl + << "---- Upwinding VOF margin : " << m_upw_margin + << std::endl + << "---- Pseudo CFL: " << m_pCFL << std::endl; + if (m_verbose > 1) { + // Less important or less used parameters + amrex::Print() + << "---- Calc. conv. interval : " << m_calc_convg_interval + << std::endl + << "---- Target field cutoff : " << m_target_cutoff + << std::endl; + } } amrex::Print() << std::endl; } @@ -195,8 +207,9 @@ void OversetOps::sharpen_nalu_data() const auto nlevels = repo.num_active_levels(); const auto geom = m_sim_ptr->mesh().Geom(); - // Get blanking for cells + // Get blanking for cells and nodes const auto& iblank_cell = repo.get_int_field("iblank_cell"); + const auto& iblank_node = repo.get_int_field("iblank_node"); // Get fields that will be modified auto& vof = repo.get_field("vof"); @@ -214,6 +227,7 @@ void OversetOps::sharpen_nalu_data() 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]); + auto delta_vof = repo.create_scratch_field(1); // Sharpening fluxes (at faces) have 1 ghost, requiring fields to have >= 2 auto gp_scr = repo.create_scratch_field(3, 2); @@ -227,6 +241,9 @@ void OversetOps::sharpen_nalu_data() // Prep things that do not change with iterations for (int lev = 0; lev < nlevels; ++lev) { + // Ensure vof is bounded after interpolation from nalu-wind + overset_ops::process_vof(vof(lev), 0.); + // Thickness used here is user parameter, whatever works best const auto dx = (geom[lev]).CellSizeArray(); const amrex::Real i_th = @@ -253,17 +270,24 @@ void OversetOps::sharpen_nalu_data() // Process target vof for tiny margins from single-phase for (int lev = 0; lev < nlevels; ++lev) { - // A tolerance of 0 should do nothing overset_ops::process_vof((*target_vof)(lev), m_target_cutoff); } amrex::Gpu::streamSynchronize(); - // Replace vof with original values in amr domain + amrex::Real target_err0 = 0.; for (int lev = 0; lev < nlevels; ++lev) { + // Replace vof with original values in amr domain overset_ops::harmonize_vof( (*target_vof)(lev), vof(lev), iblank_cell(lev)); + // Get initial error from target field + const amrex::Real target_err_lev = + overset_ops::measure_target_convergence( + (*target_vof)(lev), vof(lev)); + target_err0 += target_err_lev; } amrex::Gpu::streamSynchronize(); + amrex::ParallelDescriptor::ReduceRealSum(target_err0); + amrex::Real target_err_last = target_err0; // Put fluxes in vector for averaging down during iterations amrex::Vector> fluxes( @@ -275,18 +299,16 @@ void OversetOps::sharpen_nalu_data() } // Pseudo-time loop - amrex::Real err = 100.0 * m_convg_tol; + amrex::Real err = target_err0; + amrex::Real target_err = err; int n = 0; + // Will skip if there is little difference between target and current while (n < m_n_iterations && err > m_convg_tol) { ++n; bool calc_convg = n % m_calc_convg_interval == 0; // Zero error if being calculated this step err = calc_convg ? 0.0 : err; - - // Maximum possible value of pseudo time factor (dtau) - amrex::Real ptfac = 1.0; - // Maximum pseudoCFL, 0.5 seems to work well - const amrex::Real pCFL = 0.5; + target_err = calc_convg ? 0.0 : target_err; for (int lev = 0; lev < nlevels; ++lev) { // Populate normal vector @@ -303,13 +325,13 @@ void OversetOps::sharpen_nalu_data() // Process fluxes overset_ops::process_fluxes_calc_src( (*flux_x)(lev), (*flux_y)(lev), (*flux_z)(lev), (*p_src)(lev), - vof(lev), iblank_cell(lev)); + vof(lev), iblank_cell(lev), iblank_node(lev)); // 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( + overset_ops::measure_flux_convergence( (*flux_x)(lev), (*flux_y)(lev), (*flux_z)(lev)) / pvscale; err = amrex::max(err, err_lev); @@ -324,21 +346,38 @@ void OversetOps::sharpen_nalu_data() repo.mesh().refRatio(lev - 1), geom[lev - 1]); } + if (calc_convg) { + target_err = 0.; + } + // Get pseudo dt (dtau) + amrex::Real ptfac = 1.0; for (int lev = 0; lev < nlevels; ++lev) { + amrex::iMultiFab level_mask; + if (lev < nlevels - 1) { + level_mask = makeFineMask( + repo.mesh().boxArray(lev), repo.mesh().DistributionMap(lev), + repo.mesh().boxArray(lev + 1), repo.mesh().refRatio(lev), 1, + 0); + } else { + level_mask.define( + repo.mesh().boxArray(lev), repo.mesh().DistributionMap(lev), + 1, 0, amrex::MFInfo()); + level_mask.setVal(1); + } const auto dx = (geom[lev]).CellSizeArray(); + (*delta_vof)(lev).setVal(0.); // 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), dx, - m_convg_tol); + (*flux_x)(lev), (*flux_y)(lev), (*flux_z)(lev), vof(lev), + (*delta_vof)(lev), iblank_cell(lev), level_mask, dx); ptfac = amrex::min(ptfac, ptfac_lev); } amrex::Gpu::streamSynchronize(); amrex::ParallelDescriptor::ReduceRealMin(ptfac); - - // Conform pseudo dt (dtau) to pseudo CFL - ptfac = pCFL * ptfac; + ptfac *= m_pCFL; + ptfac = amrex::max(0., ptfac); // Apply fluxes for (int lev = 0; lev < nlevels; ++lev) { @@ -352,6 +391,13 @@ void OversetOps::sharpen_nalu_data() vof(lev).FillBoundary(geom[lev].periodicity()); velocity(lev).FillBoundary(geom[lev].periodicity()); gp(lev).FillBoundary(geom[lev].periodicity()); + + if (calc_convg) { + const amrex::Real target_err_lev = + overset_ops::measure_target_convergence( + (*target_vof)(lev), vof(lev)); + target_err += target_err_lev; + } } amrex::Gpu::streamSynchronize(); @@ -361,12 +407,39 @@ void OversetOps::sharpen_nalu_data() // Ensure that err is same across processors if (calc_convg) { amrex::ParallelDescriptor::ReduceRealMax(err); + amrex::ParallelDescriptor::ReduceRealSum(target_err); } if (m_verbose > 0) { - amrex::Print() << "OversetOps: sharpen step " << n << " conv. err " - << err << " tol " << m_convg_tol << std::endl; + amrex::Print() << "OversetOps: sharpen step " << std::setw(2) << n + << " max vof flux " << std::scientific + << std::setprecision(4) << err << " targ_err " + << target_err << " /init " + << target_err / target_err0 << " p-dt " << ptfac + << std::endl; + } + if (target_err > target_err_last * (1.0 + constants::LOOSE_TOL)) { + amrex::Print() << "OversetOps: WARNING, target error increased " + << target_err_last << " -> " << target_err + << std::endl; } + target_err_last = target_err; + } + + // Finish with consistency across levels + for (int lev = nlevels - 1; lev > 0; --lev) { + amrex::average_down( + velocity(lev), velocity(lev - 1), 0, AMREX_SPACEDIM, + repo.mesh().refRatio(lev - 1)); + velocity(lev - 1).FillBoundary(geom[lev - 1].periodicity()); + amrex::average_down( + vof(lev), vof(lev - 1), 0, 1, repo.mesh().refRatio(lev - 1)); + vof(lev - 1).FillBoundary(geom[lev - 1].periodicity()); + amrex::average_down( + rho(lev), rho(lev - 1), 0, 1, repo.mesh().refRatio(lev - 1)); + rho(lev - 1).FillBoundary(geom[lev - 1].periodicity()); + amrex::average_down_nodal( + p(lev), p(lev - 1), repo.mesh().refRatio(lev - 1)); } // Fillpatch for pressure to make sure pressure stencil has all points @@ -374,9 +447,10 @@ void OversetOps::sharpen_nalu_data() // 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 + // evolves) with the difference b/w post-sharpening vof and target + field_ops::lincomb(*target_vof, 1., *target_vof, 0, -1., vof, 0, 0, 1, 0); for (int lev = 0; lev < nlevels; ++lev) { - overset_ops::equate_field(levelset(lev), velocity(lev)); + overset_ops::equate_field(levelset(lev), (*target_vof)(lev), 0); } amrex::Gpu::streamSynchronize(); } diff --git a/amr-wind/overset/overset_ops_K.H b/amr-wind/overset/overset_ops_K.H index 4c416ef3a9..3ad4ac6785 100644 --- a/amr-wind/overset/overset_ops_K.H +++ b/amr-wind/overset/overset_ops_K.H @@ -58,6 +58,9 @@ amrex::Real AMREX_GPU_DEVICE AMREX_FORCE_INLINE alpha_flux( : dphi_eval; } } + // Limit vof fluxes based on neighbors + dphi_eval = std::min(vof(iv), std::max(-vof(ivm), dphi_eval)); + dphi_eval = std::min(1. - vof(ivm), std::max(-(1. - vof(iv)), dphi_eval)); return dphi_eval; } diff --git a/amr-wind/overset/overset_ops_routines.H b/amr-wind/overset/overset_ops_routines.H index 1aa47ae248..8afa060af5 100644 --- a/amr-wind/overset/overset_ops_routines.H +++ b/amr-wind/overset/overset_ops_routines.H @@ -65,21 +65,24 @@ void process_fluxes_calc_src( amrex::MultiFab& mf_fz, amrex::MultiFab& mf_psource, const amrex::MultiFab& mf_vof, - const amrex::iMultiFab& mf_iblank); + const amrex::iMultiFab& mf_iblank_cell, + const amrex::iMultiFab& mf_iblank_node); amrex::Real calculate_pseudo_velocity_scale( const amrex::iMultiFab& mf_iblank, const amrex::GpuArray dx, const amrex::Real pvmax); -// Calculate a type of CFL by measuring how much % VOF is being removed per cell +// Calculate a type of CFL using fluxes and bounds of vof amrex::Real calculate_pseudo_dt_flux( const amrex::MultiFab& mf_fx, const amrex::MultiFab& mf_fy, const amrex::MultiFab& mf_fz, const amrex::MultiFab& mf_vof, - const amrex::GpuArray& dx, - const amrex::Real tol); + amrex::MultiFab& mf_dvof, + const amrex::iMultiFab& mf_iblank, + const amrex::iMultiFab& mf_lmask, + const amrex::GpuArray& dx); // Apply reinitialization fluxes to modify fields void apply_fluxes( @@ -97,11 +100,16 @@ void apply_fluxes( const amrex::Real vof_tol); // Get the size of the smallest VOF flux to quantify convergence -amrex::Real measure_convergence( +amrex::Real measure_flux_convergence( amrex::MultiFab& mf_fx, amrex::MultiFab& mf_fy, amrex::MultiFab& mf_fz); +// Get the difference between the current vof and the target vof +amrex::Real measure_target_convergence( + amrex::MultiFab& mf_vof_target, amrex::MultiFab& mf_vof); + // Set levelset field to another quantity to view in plotfile for debugging -void equate_field(amrex::MultiFab& mf_dest, const amrex::MultiFab& mf_src); +void equate_field( + amrex::MultiFab& mf_dest, const amrex::MultiFab& mf_src, const int icomp); // Swap pressure gradient values in overset region void replace_gradp( diff --git a/amr-wind/overset/overset_ops_routines.cpp b/amr-wind/overset/overset_ops_routines.cpp index 6e1b3b58e4..cfa67a8454 100644 --- a/amr-wind/overset/overset_ops_routines.cpp +++ b/amr-wind/overset/overset_ops_routines.cpp @@ -233,24 +233,10 @@ void populate_normal_vector( amrex::ParallelFor( mf_normvec, mf_normvec.n_grow - amrex::IntVect(1), [=] AMREX_GPU_DEVICE(int nbx, int i, int j, int k) noexcept { - // Neumann condition across nalu bdy - int ibdy = - (iblank[nbx](i, j, k) != iblank[nbx](i - 1, j, k)) ? -1 : 0; - int jbdy = - (iblank[nbx](i, j, k) != iblank[nbx](i, j - 1, k)) ? -1 : 0; - int kbdy = - (iblank[nbx](i, j, k) != iblank[nbx](i, j, k - 1)) ? -1 : 0; - // no cell should be isolated such that -1 and 1 are needed - ibdy = - (iblank[nbx](i, j, k) != iblank[nbx](i + 1, j, k)) ? +1 : ibdy; - jbdy = - (iblank[nbx](i, j, k) != iblank[nbx](i, j + 1, k)) ? +1 : jbdy; - kbdy = - (iblank[nbx](i, j, k) != iblank[nbx](i, j, k + 1)) ? +1 : kbdy; // Calculate normal amrex::Real mx, my, mz, mmag; multiphase::youngs_finite_difference_normal_neumann( - i, j, k, ibdy, jbdy, kbdy, vof[nbx], mx, my, mz); + i, j, k, iblank[nbx], vof[nbx], mx, my, mz); // Normalize normal mmag = std::sqrt(mx * mx + my * my + mz * mz + 1e-20); // Save normal @@ -362,14 +348,16 @@ void process_fluxes_calc_src( amrex::MultiFab& mf_fz, amrex::MultiFab& mf_psource, const amrex::MultiFab& mf_vof, - const amrex::iMultiFab& mf_iblank) + const amrex::iMultiFab& mf_iblank_cell, + const amrex::iMultiFab& mf_iblank_node) { const auto& fx = mf_fx.arrays(); const auto& fy = mf_fy.arrays(); const auto& fz = mf_fz.arrays(); const auto& sp = mf_psource.arrays(); const auto& vof = mf_vof.const_arrays(); - const auto& iblank = mf_iblank.const_arrays(); + const auto& iblank = mf_iblank_cell.const_arrays(); + const auto& ibl_nd = mf_iblank_node.const_arrays(); constexpr amrex::Real tiny = std::numeric_limits::epsilon(); // Zero fluxes based on iblank array amrex::ParallelFor( @@ -398,9 +386,13 @@ void process_fluxes_calc_src( mf_psource, [=] AMREX_GPU_DEVICE(int nbx, int i, int j, int k) noexcept { sp[nbx](i, j, k) = - gp_flux_tensor(i, j, k, fx[nbx], fy[nbx], fz[nbx], tiny) && - normal_reinit_tensor( - i, j, k, fx[nbx], fy[nbx], fz[nbx], vof[nbx], tiny); + ibl_nd[nbx](i, j, k) == -1 + ? gp_flux_tensor( + i, j, k, fx[nbx], fy[nbx], fz[nbx], tiny) && + normal_reinit_tensor( + i, j, k, fx[nbx], fy[nbx], fz[nbx], vof[nbx], + tiny) + : 0.; }); } @@ -415,74 +407,63 @@ amrex::Real calculate_pseudo_velocity_scale( : pvmax; } -// Calculate a type of CFL by measuring how much % VOF is being removed per cell +// Calculate a type of CFL using bounds of VOF amrex::Real calculate_pseudo_dt_flux( const amrex::MultiFab& mf_fx, const amrex::MultiFab& mf_fy, const amrex::MultiFab& mf_fz, const amrex::MultiFab& mf_vof, - const amrex::GpuArray& dx, - const amrex::Real tol) + amrex::MultiFab& mf_dvof, + const amrex::iMultiFab& mf_iblank, + const amrex::iMultiFab& mf_lmask, + const amrex::GpuArray& dx) { - // Get the maximum flux magnitude, but just for vof fluxes - const amrex::Real pdt_fx = amrex::ReduceMin( - mf_fx, mf_vof, 0, - [=] AMREX_GPU_HOST_DEVICE( - amrex::Box const& bx, amrex::Array4 const& fx, - amrex::Array4 const& vof) -> amrex::Real { - amrex::Real pdt_fab = 1.0; - amrex::Loop(bx, [=, &pdt_fab](int i, int j, int k) noexcept { - 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) * 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) * dx[0] / -fx(i, j, k, 0); - } - pdt_fab = amrex::min(pdt_fab, pdt_lim); - }); - return pdt_fab; - }); - const amrex::Real pdt_fy = amrex::ReduceMin( - mf_fy, mf_vof, 0, - [=] AMREX_GPU_HOST_DEVICE( - amrex::Box const& bx, amrex::Array4 const& fy, - amrex::Array4 const& vof) -> amrex::Real { - amrex::Real pdt_fab = 1.0; - amrex::Loop(bx, [=, &pdt_fab](int i, int j, int k) noexcept { - 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) * 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) * dx[1] / -fy(i, j, k, 0); - } - pdt_fab = amrex::min(pdt_fab, pdt_lim); - }); - return pdt_fab; + const auto& fx = mf_fx.const_arrays(); + const auto& fy = mf_fy.const_arrays(); + const auto& fz = mf_fz.const_arrays(); + const auto& ibc = mf_iblank.const_arrays(); + const auto& lmsk = mf_lmask.const_arrays(); + auto dvof_arrs = mf_dvof.arrays(); + + // Apply VOF fluxes fully + amrex::ParallelFor( + mf_vof, [=] AMREX_GPU_DEVICE(int nbx, int i, int j, int k) noexcept { + const amrex::Real ibfac = ibc[nbx](i, j, k) == -1 ? 1.0 : 0.0; + const amrex::Real lmfac = lmsk[nbx](i, j, k) == 1 ? 1.0 : 0.0; + dvof_arrs[nbx](i, j, k) = + ibfac * lmfac * + ((fx[nbx](i + 1, j, k, 0) - fx[nbx](i, j, k, 0)) / dx[0] + + (fy[nbx](i, j + 1, k, 0) - fy[nbx](i, j, k, 0)) / dx[1] + + (fz[nbx](i, j, k + 1, 0) - fz[nbx](i, j, k, 0)) / dx[2]); }); - const amrex::Real pdt_fz = amrex::ReduceMin( - mf_fz, mf_vof, 0, + // Get get pseudo-dt based on overshoots, undershoots + const amrex::Real pdt = amrex::ReduceMin( + mf_vof, mf_dvof, mf_iblank, 0, [=] AMREX_GPU_HOST_DEVICE( - amrex::Box const& bx, amrex::Array4 const& fz, - amrex::Array4 const& vof) -> amrex::Real { + amrex::Box const& bx, amrex::Array4 const& vof, + amrex::Array4 const& dvof, + amrex::Array4 const& iblank) -> amrex::Real { amrex::Real pdt_fab = 1.0; amrex::Loop(bx, [=, &pdt_fab](int i, int j, int k) noexcept { 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) * 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) * dx[2] / -fz(i, j, k, 0); + if (iblank(i, j, k) == -1) { + pdt_lim = vof(i, j, k) + dvof(i, j, k) < -constants::EPS + ? amrex::min( + -vof(i, j, k) / + (dvof(i, j, k) - constants::EPS), + pdt_lim) + : pdt_lim; + pdt_lim = vof(i, j, k) + dvof(i, j, k) > 1. + constants::EPS + ? amrex::min( + (1. - vof(i, j, k)) / + (dvof(i, j, k) + constants::EPS), + pdt_lim) + : pdt_lim; } pdt_fab = amrex::min(pdt_fab, pdt_lim); }); return pdt_fab; }); - const amrex::Real pdt = amrex::min(pdt_fx, amrex::min(pdt_fy, pdt_fz)); return pdt; } @@ -576,8 +557,8 @@ void apply_fluxes( }); } -// Get the size of the smallest VOF flux to quantify convergence -amrex::Real measure_convergence( +// Get the size of the largest VOF flux to quantify convergence +amrex::Real measure_flux_convergence( amrex::MultiFab& mf_fx, amrex::MultiFab& mf_fy, amrex::MultiFab& mf_fz) { // Get the maximum flux magnitude, but just for vof fluxes @@ -618,17 +599,33 @@ amrex::Real measure_convergence( return err; } +// Get the difference between current vof and target +amrex::Real measure_target_convergence( + amrex::MultiFab& mf_vof_target, amrex::MultiFab& mf_vof) +{ + // Get the maximum flux magnitude, but just for vof fluxes + return amrex::ReduceSum( + mf_vof, mf_vof_target, 0, + [=] AMREX_GPU_HOST_DEVICE( + amrex::Box const& bx, amrex::Array4 const& vof, + amrex::Array4 const& targ) -> amrex::Real { + amrex::Real err_fab = 0.0; + amrex::Loop(bx, [=, &err_fab](int i, int j, int k) noexcept { + err_fab += std::abs(vof(i, j, k) - targ(i, j, k)); + }); + return err_fab; + }); +} + // Set levelset field to another quantity to view in plotfile for debugging -void equate_field(amrex::MultiFab& mf_dest, const amrex::MultiFab& mf_src) +void equate_field( + amrex::MultiFab& mf_dest, const amrex::MultiFab& mf_src, const int icomp) { const auto& dest = mf_dest.arrays(); const auto& src = mf_src.const_arrays(); amrex::ParallelFor( mf_dest, [=] AMREX_GPU_DEVICE(int nbx, int i, int j, int k) noexcept { - dest[nbx](i, j, k) = std::sqrt( - src[nbx](i, j, k, 0) * src[nbx](i, j, k, 0) + - src[nbx](i, j, k, 1) * src[nbx](i, j, k, 1) + - src[nbx](i, j, k, 2) * src[nbx](i, j, k, 2)); + dest[nbx](i, j, k) = src[nbx](i, j, k, icomp); }); } diff --git a/amr-wind/physics/multiphase/MultiPhase.cpp b/amr-wind/physics/multiphase/MultiPhase.cpp index ee908b5d35..9eb1781927 100644 --- a/amr-wind/physics/multiphase/MultiPhase.cpp +++ b/amr-wind/physics/multiphase/MultiPhase.cpp @@ -554,36 +554,9 @@ void MultiPhase::levelset2vof( amrex::ParallelFor( levelset, [=] AMREX_GPU_DEVICE(int nbx, int i, int j, int k) noexcept { - // Neumann of levelset across iblank boundaries - int ibdy = - (iblank_arrs[nbx](i, j, k) != iblank_arrs[nbx](i - 1, j, k)) - ? -1 - : 0; - int jbdy = - (iblank_arrs[nbx](i, j, k) != iblank_arrs[nbx](i, j - 1, k)) - ? -1 - : 0; - int kbdy = - (iblank_arrs[nbx](i, j, k) != iblank_arrs[nbx](i, j, k - 1)) - ? -1 - : 0; - // no cell should be isolated such that -1 and 1 are - // needed - ibdy = - (iblank_arrs[nbx](i, j, k) != iblank_arrs[nbx](i + 1, j, k)) - ? +1 - : ibdy; - jbdy = - (iblank_arrs[nbx](i, j, k) != iblank_arrs[nbx](i, j + 1, k)) - ? +1 - : jbdy; - kbdy = - (iblank_arrs[nbx](i, j, k) != iblank_arrs[nbx](i, j, k + 1)) - ? +1 - : kbdy; amrex::Real mx, my, mz; multiphase::youngs_finite_difference_normal_neumann( - i, j, k, ibdy, jbdy, kbdy, phi_arrs[nbx], mx, my, mz); + i, j, k, iblank_arrs[nbx], phi_arrs[nbx], mx, my, mz); mx = std::abs(mx / 32.); my = std::abs(my / 32.); mz = std::abs(mz / 32.); diff --git a/amr-wind/physics/multiphase/SloshingTank.cpp b/amr-wind/physics/multiphase/SloshingTank.cpp index b11b497af8..d4b0a816b1 100644 --- a/amr-wind/physics/multiphase/SloshingTank.cpp +++ b/amr-wind/physics/multiphase/SloshingTank.cpp @@ -77,9 +77,9 @@ void SloshingTank::initialize_fields(int level, const amrex::Geometry& geom) -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 - amrex::Real ih_g = + const amrex::Real ih_g = amrex::max(0.0, amrex::min(probhi[2] - z0, probhi[2] - z)); - amrex::Real ih_l = + 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; diff --git a/unit_tests/multiphase/test_vof_overset_ops.cpp b/unit_tests/multiphase/test_vof_overset_ops.cpp index 04b21b703f..2786559e85 100644 --- a/unit_tests/multiphase/test_vof_overset_ops.cpp +++ b/unit_tests/multiphase/test_vof_overset_ops.cpp @@ -924,14 +924,17 @@ TEST_F(VOFOversetOps, pseudo_vscale_dt) auto& repo = sim().repo(); const int nghost = 3; auto& vof = repo.declare_field("vof", 1, nghost); + auto& dvof = repo.declare_field("dvof", 1, 0); auto& tg_vof = repo.declare_field("target_vof", 1, nghost); auto& norm = repo.declare_field("int_normal", 3, nghost); auto& iblank = repo.declare_int_field("iblank_cell", 1, nghost); + auto& level_mask = repo.declare_int_field("level_mask", 1, nghost); iblank.setVal(-1); + level_mask.setVal(1); constexpr amrex::Real margin = 0.1; - constexpr amrex::Real convg_tol = 1e-8; - // With vof and target_vof arrays, max vof removed is 50%, doubling pdt - constexpr amrex::Real pdt_answer = 2.0; + // With vof and target_vof arrays, resulting alpha flux means dt is limited + // in cell of index 2, with initial of of 0.2, fluxes of -(0.1 + 0.04) / dx + constexpr amrex::Real pdt_answer = (0.2 / (0.1 + 0.04) * 8) / 8.; // 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 = @@ -970,8 +973,8 @@ 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), dx_lev0, - convg_tol) / + flux_x(lev), flux_y(lev), flux_z(lev), vof(lev), dvof(lev), + iblank(lev), level_mask(lev), dx_lev0) / pvscale; ptfac = amrex::min(ptfac, ptfac_lev); } @@ -991,8 +994,8 @@ 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), dx_lev0, - convg_tol) / + flux_x(lev), flux_y(lev), flux_z(lev), vof(lev), dvof(lev), + iblank(lev), level_mask(lev), dx_lev0) / pvscale; ptfac = amrex::min(ptfac, ptfac_lev); } @@ -1012,8 +1015,8 @@ 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), dx_lev0, - convg_tol) / + flux_x(lev), flux_y(lev), flux_z(lev), vof(lev), dvof(lev), + iblank(lev), level_mask(lev), dx_lev0) / pvscale; ptfac = amrex::min(ptfac, ptfac_lev); } @@ -1025,8 +1028,8 @@ 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), dx_lev0, - convg_tol) / + flux_x(lev), flux_y(lev), flux_z(lev), vof(lev), dvof(lev), + iblank(lev), level_mask(lev), dx_lev0) / pvscale; ptfac = amrex::min(ptfac, ptfac_lev); } diff --git a/unit_tests/multiphase/test_vof_plic.cpp b/unit_tests/multiphase/test_vof_plic.cpp index 8d12e54632..b46f45fd60 100644 --- a/unit_tests/multiphase/test_vof_plic.cpp +++ b/unit_tests/multiphase/test_vof_plic.cpp @@ -245,40 +245,26 @@ amrex::Real normal_vector_neumann_test_impl( amrex::Real error = 0.0; amrex::Loop(bx, [=, &error](int i, int j, int k) noexcept { - int ibdy = - (iblank(i, j, k) != iblank(i - 1, j, k)) ? -1 : 0; - int jbdy = - (iblank(i, j, k) != iblank(i, j - 1, k)) ? -1 : 0; - int kbdy = - (iblank(i, j, k) != iblank(i, j, k - 1)) ? -1 : 0; - ibdy = (iblank(i, j, k) != iblank(i + 1, j, k)) ? +1 : ibdy; - jbdy = (iblank(i, j, k) != iblank(i, j + 1, k)) ? +1 : jbdy; - kbdy = (iblank(i, j, k) != iblank(i, j, k + 1)) ? +1 : kbdy; amrex::Real mxn, myn, mzn; amr_wind::multiphase:: youngs_finite_difference_normal_neumann( - i, j, k, ibdy, jbdy, kbdy, vof_arr, mxn, myn, mzn); + i, j, k, iblank, vof_arr, mxn, myn, mzn); amrex::Real mx, my, mz; amr_wind::multiphase::youngs_finite_difference_normal( i, j, k, vof_arr, mx, my, mz); - // Use L1 norm, check against non-neumann implementation - // Slope across overset boundary should be different - constexpr amrex::Real slp_tol = 1e-8; - if (ibdy != 0) { - // x slope should be different - error += std::abs(mx - mxn) > slp_tol ? 0. : 1.0; - } - if (jbdy != 0) { - // y slope should be different - error += std::abs(my - myn) > slp_tol ? 0. : 1.0; - } - if (kbdy != 0) { - // z slope should be different - error += std::abs(mz - mzn) > slp_tol ? 0. : 1.0; + // Look for iblank == 1 in surrounding cells (norm stencil) + amrex::Real ibl_sum = 0.; + for (int n0 = -1; n0 <= 1; ++n0) { + for (int n1 = -1; n1 <= 1; ++n1) { + for (int n2 = -1; n2 <= 1; ++n2) { + ibl_sum += iblank(i + n0, j + n1, k + n2); + } + } } - // Slope should otherwise be the same - if (ibdy == 0 && jbdy == 0 && kbdy == 0) { + + // Assume there will be no error when there are enough cells + if (ibl_sum < -10) { error += std::abs(mx - mxn); error += std::abs(my - myn); error += std::abs(mz - mzn); @@ -308,19 +294,10 @@ amrex::Real normal_vector_neumann_test_impl( amrex::Real error = 0.0; amrex::Loop(bx, [=, &error](int i, int j, int k) noexcept { - int ibdy = - (iblank(i, j, k) != iblank(i - 1, j, k)) ? -1 : 0; - int jbdy = - (iblank(i, j, k) != iblank(i, j - 1, k)) ? -1 : 0; - int kbdy = - (iblank(i, j, k) != iblank(i, j, k - 1)) ? -1 : 0; - ibdy = (iblank(i, j, k) != iblank(i + 1, j, k)) ? +1 : ibdy; - jbdy = (iblank(i, j, k) != iblank(i, j + 1, k)) ? +1 : jbdy; - kbdy = (iblank(i, j, k) != iblank(i, j, k + 1)) ? +1 : kbdy; amrex::Real mxn, myn, mzn; amr_wind::multiphase:: youngs_finite_difference_normal_neumann( - i, j, k, ibdy, jbdy, kbdy, vof_arr, mxn, myn, mzn); + i, j, k, iblank, vof_arr, mxn, myn, mzn); // Use L1 norm, check for 0 if (iblank(i, j, k) == 1) {