diff --git a/amr-wind/overset/CMakeLists.txt b/amr-wind/overset/CMakeLists.txt index 757b857a1b..8874741902 100644 --- a/amr-wind/overset/CMakeLists.txt +++ b/amr-wind/overset/CMakeLists.txt @@ -2,4 +2,5 @@ target_sources(${amr_wind_lib_name} PRIVATE TiogaInterface.cpp OversetOps.cpp + overset_ops_routines.cpp ) diff --git a/amr-wind/overset/OversetOps.H b/amr-wind/overset/OversetOps.H index 0658019d4d..ae935202eb 100644 --- a/amr-wind/overset/OversetOps.H +++ b/amr-wind/overset/OversetOps.H @@ -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}; diff --git a/amr-wind/overset/OversetOps.cpp b/amr-wind/overset/OversetOps.cpp index 3c03ec7498..361a97014d 100644 --- a/amr-wind/overset/OversetOps.cpp +++ b/amr-wind/overset/OversetOps.cpp @@ -22,6 +22,7 @@ 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); @@ -29,12 +30,25 @@ void OversetOps::initialize(CFDSim& sim) 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(); + + // 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); @@ -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(); @@ -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( @@ -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 @@ -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); @@ -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) { @@ -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) @@ -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(); @@ -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 @@ -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(); @@ -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 @@ -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()); } } diff --git a/amr-wind/overset/overset_ops_K.H b/amr-wind/overset/overset_ops_K.H index 6b03c6a63e..143731ad36 100644 --- a/amr-wind/overset/overset_ops_K.H +++ b/amr-wind/overset/overset_ops_K.H @@ -1,6 +1,7 @@ #ifndef OVERSET_OPS_K_H_ #define OVERSET_OPS_K_H_ +#include "amr-wind/core/vs/vector_space.H" #include namespace amr_wind::overset_ops { @@ -103,6 +104,164 @@ void AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE velocity_face( } } +void AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE gp_rho_face( + const int i, + const int j, + const int k, + const int dir, + amrex::Array4 const& vof, + amrex::Array4 const& gp, + amrex::Array4 const& rho, + amrex::Real& uface, + amrex::Real& vface, + amrex::Real& wface) +{ + // Set up neighbor indices + const amrex::IntVect iv{i, j, k}; + const amrex::IntVect dv{(int)(dir == 0), (int)(dir == 1), (int)(dir == 2)}; + const amrex::IntVect ivm = iv - dv; + + // Gradient of phi normal to interface + const amrex::Real gphi = (vof(iv) - vof(ivm)); + + // Get velocities on both sides + const amrex::Real u_ = gp(iv, 0) / rho(iv); + const amrex::Real v_ = gp(iv, 1) / rho(iv); + const amrex::Real w_ = gp(iv, 2) / rho(iv); + const amrex::Real u_nb = gp(ivm, 0) / rho(ivm); + const amrex::Real v_nb = gp(ivm, 1) / rho(ivm); + const amrex::Real w_nb = gp(ivm, 2) / rho(ivm); + // Average value when gphi = 0 + uface = 0.5 * (u_ + u_nb); + vface = 0.5 * (v_ + v_nb); + wface = 0.5 * (w_ + w_nb); + // Use simple upwinding + if (gphi > 0.0) { + uface = u_nb; + vface = v_nb; + wface = w_nb; + } + if (gphi < 0.0) { + uface = u_; + vface = v_; + wface = w_; + } +} + +vs::Tensor AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE gp_flux_tensor( + const int i, + const int j, + const int k, + amrex::Array4 const& fx, + amrex::Array4 const& fy, + amrex::Array4 const& fz, + const amrex::Real tiny) +{ + // Averaging depends on iblank array, encoded in 8th component + const amrex::Real avg_fx = fx(i, j, k, 8) + fx(i, j, k - 1, 8) + + fx(i, j - 1, k, 8) + fx(i, j - 1, k - 1, 8) + + tiny; + const amrex::Real avg_fy = fy(i, j, k, 8) + fy(i - 1, j, k, 8) + + fy(i, j, k - 1, 8) + fy(i - 1, j, k - 1, 8) + + tiny; + const amrex::Real avg_fz = fz(i, j, k, 8) + fz(i - 1, j, k, 8) + + fz(i, j - 1, k, 8) + fz(i - 1, j - 1, k, 8) + + tiny; + + auto f_otimes_gradp = vs::Tensor::I(); + // Average fluxes (faces) to nodes where pressure exists + f_otimes_gradp.xx() = (fx(i, j, k, 5) + fx(i, j, k - 1, 5) + + fx(i, j - 1, k, 5) + fx(i, j - 1, k - 1, 5)) / + avg_fx; + f_otimes_gradp.xy() = (fx(i, j, k, 6) + fx(i, j, k - 1, 6) + + fx(i, j - 1, k, 6) + fx(i, j - 1, k - 1, 6)) / + avg_fx; + f_otimes_gradp.xz() = (fx(i, j, k, 7) + fx(i, j, k - 1, 7) + + fx(i, j - 1, k, 7) + fx(i, j - 1, k - 1, 7)) / + avg_fx; + f_otimes_gradp.yx() = (fy(i, j, k, 5) + fy(i - 1, j, k, 5) + + fy(i, j, k - 1, 5) + fy(i - 1, j, k - 1, 5)) / + avg_fy; + f_otimes_gradp.yy() = (fy(i, j, k, 6) + fy(i - 1, j, k, 6) + + fy(i, j, k - 1, 6) + fy(i - 1, j, k - 1, 6)) / + avg_fy; + f_otimes_gradp.yz() = (fy(i, j, k, 7) + fy(i - 1, j, k, 7) + + fy(i, j, k - 1, 7) + fy(i - 1, j, k - 1, 7)) / + avg_fy; + f_otimes_gradp.zx() = (fz(i, j, k, 5) + fz(i - 1, j, k, 5) + + fz(i, j - 1, k, 5) + fz(i - 1, j - 1, k, 5)) / + avg_fz; + f_otimes_gradp.zy() = (fz(i, j, k, 6) + fz(i - 1, j, k, 6) + + fz(i, j - 1, k, 6) + fz(i - 1, j - 1, k, 6)) / + avg_fz; + f_otimes_gradp.zz() = (fz(i, j, k, 7) + fz(i - 1, j, k, 7) + + fz(i, j - 1, k, 7) + fz(i - 1, j - 1, k, 7)) / + avg_fz; + + return f_otimes_gradp; +} + +vs::Tensor AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE normal_reinit_tensor( + const int i, + const int j, + const int k, + amrex::Array4 const& fx, + amrex::Array4 const& fy, + amrex::Array4 const& fz, + amrex::Array4 const& vof, + const amrex::Real tiny) +{ + // Averaging depends on iblank array, encoded in 8th component + const amrex::Real avg_fx = fx(i, j, k, 8) + fx(i, j, k - 1, 8) + + fx(i, j - 1, k, 8) + fx(i, j - 1, k - 1, 8) + + tiny; + const amrex::Real avg_fy = fy(i, j, k, 8) + fy(i - 1, j, k, 8) + + fy(i, j, k - 1, 8) + fy(i - 1, j, k - 1, 8) + + tiny; + const amrex::Real avg_fz = fz(i, j, k, 8) + fz(i - 1, j, k, 8) + + fz(i, j - 1, k, 8) + fz(i - 1, j - 1, k, 8) + + tiny; + + auto n_zeta = vs::Vector::one(); + // Get components of normal in each position + n_zeta.x() = + ((vof(i, j, k) - vof(i - 1, j, k)) * fx(i, j, k, 8) + + (vof(i, j - 1, k) - vof(i - 1, j - 1, k)) * fx(i, j - 1, k, 8) + + (vof(i, j, k - 1) - vof(i - 1, j, k - 1)) * fx(i, j, k - 1, 8) + + (vof(i, j - 1, k - 1) - vof(i - 1, j - 1, k - 1)) * + fx(i, j - 1, k - 1, 8)) / + avg_fx; + n_zeta.y() = + ((vof(i, j, k) - vof(i, j - 1, k)) * fy(i, j, k, 8) + + (vof(i - 1, j, k) - vof(i - 1, j - 1, k)) * fy(i - 1, j, k, 8) + + (vof(i, j, k - 1) - vof(i, j - 1, k - 1)) * fy(i, j, k - 1, 8) + + (vof(i - 1, j, k - 1) - vof(i - 1, j - 1, k - 1)) * + fy(i - 1, j, k - 1, 8)) / + avg_fy; + n_zeta.z() = + ((vof(i, j, k) - vof(i, j, k - 1)) * fz(i, j, k, 8) + + (vof(i - 1, j, k) - vof(i - 1, j, k - 1)) * fz(i - 1, j, k, 8) + + (vof(i, j - 1, k) - vof(i, j - 1, k - 1)) * fz(i, j - 1, k, 8) + + (vof(i - 1, j - 1, k) - vof(i - 1, j - 1, k - 1)) * + fz(i - 1, j - 1, k, 8)) / + avg_fz; + + n_zeta.normalize(); + + // To do outer product between vectors, set up tensors + auto n_tensor = vs::Tensor::I(); + auto n_tensor_T = vs::Tensor::I(); + n_tensor.rows(n_zeta, n_zeta, n_zeta); + n_tensor_T.cols(n_zeta, n_zeta, n_zeta); + // Multiply tensors elementwise + auto n_otimes_n = vs::Tensor::I(); + n_otimes_n.rows( + n_tensor.x() * n_tensor_T.x(), n_tensor.y() * n_tensor_T.y(), + n_tensor.z() * n_tensor_T.z()); + + return n_otimes_n; +} + } // namespace amr_wind::overset_ops #endif \ No newline at end of file diff --git a/amr-wind/overset/overset_ops_routines.H b/amr-wind/overset/overset_ops_routines.H index 9e289f8ac2..de9afbdff8 100644 --- a/amr-wind/overset/overset_ops_routines.H +++ b/amr-wind/overset/overset_ops_routines.H @@ -1,6 +1,8 @@ #ifndef OVERSET_OPS_ROUTINES_H_ #define OVERSET_OPS_ROUTINES_H_ +#include "AMReX_iMultiFab.H" +#include "AMReX_MultiFab.H" #include "amr-wind/equation_systems/vof/volume_fractions.H" #include "amr-wind/overset/overset_ops_K.H" @@ -11,90 +13,22 @@ void populate_psi( amrex::MultiFab& mf_psi, const amrex::MultiFab& mf_vof, const amrex::Real i_th, - const amrex::Real asdf_tiny) -{ - const auto& psi = mf_psi.arrays(); - const auto& vof = mf_vof.const_arrays(); - amrex::ParallelFor( - mf_psi, mf_psi.n_grow, - [=] AMREX_GPU_DEVICE(int nbx, int i, int j, int k) noexcept { - psi[nbx](i, j, k) = asdf(vof[nbx](i, j, k), i_th, asdf_tiny); - }); -} + const amrex::Real asdf_tiny); // Modify a vof field to not have values that barely differ from 0 or 1 -void process_vof(amrex::MultiFab& mf_vof, const amrex::Real vof_tol) -{ - const auto& vof = mf_vof.arrays(); - amrex::ParallelFor( - mf_vof, mf_vof.n_grow, - [=] AMREX_GPU_DEVICE(int nbx, int i, int j, int k) noexcept { - vof[nbx](i, j, k) = - vof[nbx](i, j, k) < vof_tol - ? 0.0 - : (vof[nbx](i, j, k) > 1. - vof_tol ? 1. - : vof[nbx](i, j, k)); - }); -} +void process_vof(amrex::MultiFab& mf_vof, const amrex::Real vof_tol); // Combine overset target vof field with current non-overset vof field void harmonize_vof( amrex::MultiFab& mf_vof_target, const amrex::MultiFab& mf_vof_original, - const amrex::iMultiFab& mf_iblank) -{ - const auto& tg_vof = mf_vof_target.arrays(); - const auto& og_vof = mf_vof_original.const_arrays(); - const auto& iblank = mf_iblank.const_arrays(); - amrex::ParallelFor( - mf_vof_target, - [=] AMREX_GPU_DEVICE(int nbx, int i, int j, int k) noexcept { - // Replace amr-wind vof values with originals - if (iblank[nbx](i, j, k) != -1) { - tg_vof[nbx](i, j, k) = og_vof[nbx](i, j, k); - } - }); -} + const amrex::iMultiFab& mf_iblank); // Populate normal vector with special treatment of overset boundary void populate_normal_vector( amrex::MultiFab& mf_normvec, const amrex::MultiFab& mf_vof, - const amrex::iMultiFab& mf_iblank) -{ - const auto& normvec = mf_normvec.arrays(); - const auto& vof = mf_vof.const_arrays(); - const auto& iblank = mf_iblank.const_arrays(); - // Calculate gradients in each direction with centered diff - 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); - // Normalize normal - mmag = std::sqrt(mx * mx + my * my + mz * mz + 1e-20); - // Save normal - normvec[nbx](i, j, k, 0) = mx / mmag; - normvec[nbx](i, j, k, 1) = my / mmag; - normvec[nbx](i, j, k, 2) = mz / mmag; - }); -} + const amrex::iMultiFab& mf_iblank); // Calculate fluxes for reinitialization over entire domain without concern for // overset bdy @@ -106,281 +40,57 @@ void populate_sharpen_fluxes( const amrex::MultiFab& mf_target_vof, const amrex::MultiFab& mf_norm, const amrex::MultiFab& mf_velocity, + const amrex::MultiFab& mf_gp, + const amrex::MultiFab& mf_density, + const amrex::Real Gamma, const amrex::Real margin, const amrex::Real rho1, - const amrex::Real rho2) -{ - const auto& fx = mf_fx.arrays(); - const auto& fy = mf_fy.arrays(); - const auto& fz = mf_fz.arrays(); - const auto& vof = mf_vof.const_arrays(); - const auto& tg_vof = mf_target_vof.const_arrays(); - const auto& norm = mf_norm.const_arrays(); - const auto& vel = mf_velocity.const_arrays(); - amrex::ParallelFor( - mf_fx, mf_fx.n_grow, - [=] AMREX_GPU_DEVICE(int nbx, int i, int j, int k) noexcept { - // vof flux - amrex::Real flux = alpha_flux( - i, j, k, 0, margin, vof[nbx], tg_vof[nbx], norm[nbx]); - fx[nbx](i, j, k, 0) = flux; - // density flux - flux *= (rho1 - rho2); - fx[nbx](i, j, k, 1) = flux; - // momentum fluxes (dens flux * face vel) - amrex::Real uf, vf, wf; - velocity_face(i, j, k, 0, vof[nbx], vel[nbx], uf, vf, wf); - fx[nbx](i, j, k, 2) = flux * uf; - fx[nbx](i, j, k, 3) = flux * vf; - fx[nbx](i, j, k, 4) = flux * wf; - }); - amrex::ParallelFor( - mf_fy, mf_fy.n_grow, - [=] AMREX_GPU_DEVICE(int nbx, int i, int j, int k) noexcept { - amrex::Real flux = alpha_flux( - i, j, k, 1, margin, vof[nbx], tg_vof[nbx], norm[nbx]); - fy[nbx](i, j, k, 0) = flux; - flux *= (rho1 - rho2); - fy[nbx](i, j, k, 1) = flux; - amrex::Real uf, vf, wf; - velocity_face(i, j, k, 1, vof[nbx], vel[nbx], uf, vf, wf); - fy[nbx](i, j, k, 2) = flux * uf; - fy[nbx](i, j, k, 3) = flux * vf; - fy[nbx](i, j, k, 4) = flux * wf; - }); - amrex::ParallelFor( - mf_fz, mf_fz.n_grow, - [=] AMREX_GPU_DEVICE(int nbx, int i, int j, int k) noexcept { - amrex::Real flux = alpha_flux( - i, j, k, 2, margin, vof[nbx], tg_vof[nbx], norm[nbx]); - fz[nbx](i, j, k, 0) = flux; - flux *= (rho1 - rho2); - fz[nbx](i, j, k, 1) = flux; - amrex::Real uf, vf, wf; - velocity_face(i, j, k, 2, vof[nbx], vel[nbx], uf, vf, wf); - fz[nbx](i, j, k, 2) = flux * uf; - fz[nbx](i, j, k, 3) = flux * vf; - fz[nbx](i, j, k, 4) = flux * wf; - }); -} + const amrex::Real rho2); -// Process reinitialization fluxes - zero non-internal to overset region -void process_fluxes( +// Process reinitialization fluxes - zero non-internal to overset region; +// also calculate pressure source / sink term as a function of fluxes +void process_fluxes_calc_src( amrex::MultiFab& mf_fx, amrex::MultiFab& mf_fy, amrex::MultiFab& mf_fz, - const amrex::iMultiFab& mf_iblank) -{ - const auto& fx = mf_fx.arrays(); - const auto& fy = mf_fy.arrays(); - const auto& fz = mf_fz.arrays(); - const auto& iblank = mf_iblank.const_arrays(); - amrex::ParallelFor( - mf_fx, mf_fx.n_grow, mf_fx.n_comp, - [=] AMREX_GPU_DEVICE(int nbx, int i, int j, int k, int n) noexcept { - bool zero_all = - (iblank[nbx](i - 1, j, k) + iblank[nbx](i, j, k) > -2); - fx[nbx](i, j, k, n) *= zero_all ? 0. : 1.; - }); - amrex::ParallelFor( - mf_fy, mf_fy.n_grow, mf_fy.n_comp, - [=] AMREX_GPU_DEVICE(int nbx, int i, int j, int k, int n) noexcept { - bool zero_all = - (iblank[nbx](i, j - 1, k) + iblank[nbx](i, j, k) > -2); - fy[nbx](i, j, k, n) *= zero_all ? 0. : 1.; - }); - amrex::ParallelFor( - mf_fz, mf_fz.n_grow, mf_fz.n_comp, - [=] AMREX_GPU_DEVICE(int nbx, int i, int j, int k, int n) noexcept { - bool zero_all = - (iblank[nbx](i, j, k - 1) + iblank[nbx](i, j, k) > -2); - fz[nbx](i, j, k, n) *= zero_all ? 0. : 1.; - }); -} + amrex::MultiFab& mf_psource, + const amrex::MultiFab& mf_vof, + const amrex::iMultiFab& mf_iblank); + +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 amrex::Real calculate_pseudo_dt_flux( - amrex::MultiFab& mf_fx, - amrex::MultiFab& mf_fy, - amrex::MultiFab& mf_fz, - amrex::MultiFab& mf_vof, - amrex::Real tol) -{ - // 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) / 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_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) / 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_fab = amrex::min(pdt_fab, pdt_lim); - }); - return pdt_fab; - }); - const amrex::Real pdt_fz = amrex::ReduceMin( - mf_fz, mf_vof, 0, - [=] AMREX_GPU_HOST_DEVICE( - amrex::Box const& bx, amrex::Array4 const& fz, - 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 (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); - } 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_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; -} + const amrex::MultiFab& mf_fx, + const amrex::MultiFab& mf_fy, + const amrex::MultiFab& mf_fz, + const amrex::MultiFab& mf_vof, + const amrex::Real tol); // Apply reinitialization fluxes to modify fields void apply_fluxes( - amrex::MultiFab& mf_fx, - amrex::MultiFab& mf_fy, - amrex::MultiFab& mf_fz, + const amrex::MultiFab& mf_fx, + const amrex::MultiFab& mf_fy, + const amrex::MultiFab& mf_fz, + const amrex::MultiFab& mf_psource, amrex::MultiFab& mf_vof, amrex::MultiFab& mf_dens, amrex::MultiFab& mf_vel, - amrex::Real ptfac, - amrex::Real vof_tol) -{ - const auto& fx = mf_fx.const_arrays(); - const auto& fy = mf_fy.const_arrays(); - const auto& fz = mf_fz.const_arrays(); - const auto& vof = mf_vof.arrays(); - const auto& dens = mf_dens.arrays(); - const auto& vel = mf_vel.arrays(); - - amrex::ParallelFor( - mf_vof, [=] AMREX_GPU_DEVICE(int nbx, int i, int j, int k) noexcept { - const amrex::Real olddens = dens[nbx](i, j, k); - vof[nbx](i, j, k) += - ptfac * ((fx[nbx](i + 1, j, k, 0) - fx[nbx](i, j, k, 0)) + - (fy[nbx](i, j + 1, k, 0) - fy[nbx](i, j, k, 0)) + - (fz[nbx](i, j, k + 1, 0) - fz[nbx](i, j, k, 0))); - dens[nbx](i, j, k) += - ptfac * ((fx[nbx](i + 1, j, k, 1) - fx[nbx](i, j, k, 1)) + - (fy[nbx](i, j + 1, k, 1) - fy[nbx](i, j, k, 1)) + - (fz[nbx](i, j, k + 1, 1) - fz[nbx](i, j, k, 1))); - vel[nbx](i, j, k, 0) = - 1.0 / dens[nbx](i, j, k) * - (olddens * vel[nbx](i, j, k, 0) + - ptfac * ((fx[nbx](i + 1, j, k, 2) - fx[nbx](i, j, k, 2)) + - (fy[nbx](i, j + 1, k, 2) - fy[nbx](i, j, k, 2)) + - (fz[nbx](i, j, k + 1, 2) - fz[nbx](i, j, k, 2)))); - vel[nbx](i, j, k, 1) = - 1.0 / dens[nbx](i, j, k) * - (olddens * vel[nbx](i, j, k, 1) + - ptfac * ((fx[nbx](i + 1, j, k, 3) - fx[nbx](i, j, k, 3)) + - (fy[nbx](i, j + 1, k, 3) - fy[nbx](i, j, k, 3)) + - (fz[nbx](i, j, k + 1, 3) - fz[nbx](i, j, k, 3)))); - vel[nbx](i, j, k, 2) = - 1.0 / dens[nbx](i, j, k) * - (olddens * vel[nbx](i, j, k, 2) + - ptfac * ((fx[nbx](i + 1, j, k, 4) - fx[nbx](i, j, k, 4)) + - (fy[nbx](i, j + 1, k, 4) - fy[nbx](i, j, k, 4)) + - (fz[nbx](i, j, k + 1, 4) - fz[nbx](i, j, k, 4)))); - - // Ensure vof is bounded - vof[nbx](i, j, k) = - vof[nbx](i, j, k) < vof_tol - ? 0.0 - : (vof[nbx](i, j, k) > 1. - vof_tol ? 1. - : vof[nbx](i, j, k)); - // Density bounds are enforced elsewhere - }); -} + amrex::MultiFab& mf_gp, + amrex::MultiFab& mf_pressure, + const amrex::GpuArray dx, + const amrex::Real ptfac, + const amrex::Real vof_tol); // Get the size of the smallest VOF flux to quantify convergence amrex::Real measure_convergence( - amrex::MultiFab& mf_fx, amrex::MultiFab& mf_fy, amrex::MultiFab& mf_fz) -{ - // Get the maximum flux magnitude, but just for vof fluxes - const amrex::Real err_fx = amrex::ReduceMax( - mf_fx, 0, - [=] AMREX_GPU_HOST_DEVICE( - amrex::Box const& bx, - amrex::Array4 const& fx) -> amrex::Real { - amrex::Real err_fab = -1.0; - amrex::Loop(bx, [=, &err_fab](int i, int j, int k) noexcept { - err_fab = amrex::max(err_fab, std::abs(fx(i, j, k, 0))); - }); - return err_fab; - }); - const amrex::Real err_fy = amrex::ReduceMax( - mf_fy, 0, - [=] AMREX_GPU_HOST_DEVICE( - amrex::Box const& bx, - amrex::Array4 const& fy) -> amrex::Real { - amrex::Real err_fab = -1.0; - amrex::Loop(bx, [=, &err_fab](int i, int j, int k) noexcept { - err_fab = amrex::max(err_fab, std::abs(fy(i, j, k, 0))); - }); - return err_fab; - }); - const amrex::Real err_fz = amrex::ReduceMax( - mf_fz, 0, - [=] AMREX_GPU_HOST_DEVICE( - amrex::Box const& bx, - amrex::Array4 const& fz) -> amrex::Real { - amrex::Real err_fab = -1.0; - amrex::Loop(bx, [=, &err_fab](int i, int j, int k) noexcept { - err_fab = amrex::max(err_fab, std::abs(fz(i, j, k, 0))); - }); - return err_fab; - }); - const amrex::Real err = amrex::max(err_fx, amrex::max(err_fy, err_fz)); - return err; -} + amrex::MultiFab& mf_fx, amrex::MultiFab& mf_fy, amrex::MultiFab& mf_fz); // Set levelset field to another quantity to view in plotfile for debugging -void equate_field(amrex::MultiFab& mf_dest, const amrex::MultiFab& mf_src) -{ - 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)); - }); -} +void equate_field(amrex::MultiFab& mf_dest, const amrex::MultiFab& mf_src); // Replace pressure gradient with hydrostatic field in overset regions void replace_gradp_hydrostatic( @@ -389,64 +99,20 @@ void replace_gradp_hydrostatic( const amrex::MultiFab& mf_refdens, const amrex::iMultiFab& mf_iblank, const amrex::Real grav_z, - const bool is_pptb) -{ - const auto& gp = mf_gp.arrays(); - const auto& rho = mf_density.const_arrays(); - const auto& rho0 = mf_refdens.const_arrays(); - const auto& iblank = mf_iblank.const_arrays(); - amrex::ParallelFor( - mf_gp, mf_gp.n_grow, - [=] AMREX_GPU_DEVICE(int nbx, int i, int j, int k) noexcept { - if (iblank[nbx](i, j, k) == -1) { - const amrex::Real dfac = - is_pptb ? rho[nbx](i, j, k) - rho0[nbx](i, j, k) - : rho[nbx](i, j, k); - gp[nbx](i, j, k, 0) = 0.; - gp[nbx](i, j, k, 1) = 0.; - gp[nbx](i, j, k, 2) = dfac * grav_z; - } - }); -} + const bool is_pptb); // Swap pressure gradient values in overset region void replace_gradp( amrex::MultiFab& mf_gp, const amrex::MultiFab& mf_gp0, - const amrex::iMultiFab& mf_iblank) -{ - const auto gp = mf_gp.arrays(); - const auto gp0 = mf_gp0.const_arrays(); - const auto iblank = mf_iblank.const_arrays(); - amrex::ParallelFor( - mf_gp, mf_gp.n_grow, - [=] AMREX_GPU_DEVICE(int nbx, int i, int j, int k) noexcept { - if (iblank[nbx](i, j, k) == -1) { - gp[nbx](i, j, k, 0) = gp0[nbx](i, j, k, 0); - gp[nbx](i, j, k, 1) = gp0[nbx](i, j, k, 1); - gp[nbx](i, j, k, 2) = gp0[nbx](i, j, k, 2); - } - }); -} + const amrex::iMultiFab& mf_iblank); // Apply pressure gradient to velocity field void apply_pressure_gradient( amrex::MultiFab& mf_vel, const amrex::MultiFab& mf_density, const amrex::MultiFab& mf_gp, - const amrex::Real scaling_factor) -{ - const auto& vel = mf_vel.arrays(); - const auto& rho = mf_density.const_arrays(); - const auto& gp = mf_gp.const_arrays(); - amrex::ParallelFor( - mf_vel, [=] AMREX_GPU_DEVICE(int nbx, int i, int j, int k) noexcept { - const amrex::Real soverrho = scaling_factor / rho[nbx](i, j, k); - vel[nbx](i, j, k, 0) -= gp[nbx](i, j, k, 0) * soverrho; - vel[nbx](i, j, k, 1) -= gp[nbx](i, j, k, 1) * soverrho; - vel[nbx](i, j, k, 2) -= gp[nbx](i, j, k, 2) * soverrho; - }); -} + const amrex::Real scaling_factor); } // namespace amr_wind::overset_ops diff --git a/amr-wind/overset/overset_ops_routines.cpp b/amr-wind/overset/overset_ops_routines.cpp new file mode 100644 index 0000000000..57becc49ec --- /dev/null +++ b/amr-wind/overset/overset_ops_routines.cpp @@ -0,0 +1,531 @@ +#include "amr-wind/overset/overset_ops_routines.H" + +namespace amr_wind::overset_ops { +// Populate approximate signed distance function using vof field +void populate_psi( + amrex::MultiFab& mf_psi, + const amrex::MultiFab& mf_vof, + const amrex::Real i_th, + const amrex::Real asdf_tiny) +{ + const auto& psi = mf_psi.arrays(); + const auto& vof = mf_vof.const_arrays(); + amrex::ParallelFor( + mf_psi, mf_psi.n_grow, + [=] AMREX_GPU_DEVICE(int nbx, int i, int j, int k) noexcept { + psi[nbx](i, j, k) = asdf(vof[nbx](i, j, k), i_th, asdf_tiny); + }); +} + +// Modify a vof field to not have values that barely differ from 0 or 1 +void process_vof(amrex::MultiFab& mf_vof, const amrex::Real vof_tol) +{ + const auto& vof = mf_vof.arrays(); + amrex::ParallelFor( + mf_vof, mf_vof.n_grow, + [=] AMREX_GPU_DEVICE(int nbx, int i, int j, int k) noexcept { + vof[nbx](i, j, k) = + vof[nbx](i, j, k) < vof_tol + ? 0.0 + : (vof[nbx](i, j, k) > 1. - vof_tol ? 1. + : vof[nbx](i, j, k)); + }); +} + +// Combine overset target vof field with current non-overset vof field +void harmonize_vof( + amrex::MultiFab& mf_vof_target, + const amrex::MultiFab& mf_vof_original, + const amrex::iMultiFab& mf_iblank) +{ + const auto& tg_vof = mf_vof_target.arrays(); + const auto& og_vof = mf_vof_original.const_arrays(); + const auto& iblank = mf_iblank.const_arrays(); + amrex::ParallelFor( + mf_vof_target, + [=] AMREX_GPU_DEVICE(int nbx, int i, int j, int k) noexcept { + // Replace amr-wind vof values with originals + if (iblank[nbx](i, j, k) != -1) { + tg_vof[nbx](i, j, k) = og_vof[nbx](i, j, k); + } + }); +} + +// Populate normal vector with special treatment of overset boundary +void populate_normal_vector( + amrex::MultiFab& mf_normvec, + const amrex::MultiFab& mf_vof, + const amrex::iMultiFab& mf_iblank) +{ + const auto& normvec = mf_normvec.arrays(); + const auto& vof = mf_vof.const_arrays(); + const auto& iblank = mf_iblank.const_arrays(); + // Calculate gradients in each direction with centered diff + 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); + // Normalize normal + mmag = std::sqrt(mx * mx + my * my + mz * mz + 1e-20); + // Save normal + normvec[nbx](i, j, k, 0) = mx / mmag; + normvec[nbx](i, j, k, 1) = my / mmag; + normvec[nbx](i, j, k, 2) = mz / mmag; + }); +} + +// Calculate fluxes for reinitialization over entire domain without concern for +// overset bdy +void populate_sharpen_fluxes( + amrex::MultiFab& mf_fx, + amrex::MultiFab& mf_fy, + amrex::MultiFab& mf_fz, + const amrex::MultiFab& mf_vof, + const amrex::MultiFab& mf_target_vof, + const amrex::MultiFab& mf_norm, + const amrex::MultiFab& mf_velocity, + const amrex::MultiFab& mf_gp, + const amrex::MultiFab& mf_density, + const amrex::Real Gamma, + const amrex::Real margin, + const amrex::Real rho1, + const amrex::Real rho2) +{ + const auto& fx = mf_fx.arrays(); + const auto& fy = mf_fy.arrays(); + const auto& fz = mf_fz.arrays(); + const auto& vof = mf_vof.const_arrays(); + const auto& tg_vof = mf_target_vof.const_arrays(); + const auto& norm = mf_norm.const_arrays(); + const auto& vel = mf_velocity.const_arrays(); + const auto& gp = mf_gp.const_arrays(); + const auto& rho = mf_density.const_arrays(); + amrex::ParallelFor( + mf_fx, mf_fx.n_grow, + [=] AMREX_GPU_DEVICE(int nbx, int i, int j, int k) noexcept { + // vof flux + amrex::Real flux = Gamma * alpha_flux( + i, j, k, 0, margin, vof[nbx], + tg_vof[nbx], norm[nbx]); + fx[nbx](i, j, k, 0) = flux; + // density flux + flux *= (rho1 - rho2); + fx[nbx](i, j, k, 1) = flux; + // momentum fluxes (dens flux * face vel) + amrex::Real uf, vf, wf; + velocity_face(i, j, k, 0, vof[nbx], vel[nbx], uf, vf, wf); + fx[nbx](i, j, k, 2) = flux * uf; + fx[nbx](i, j, k, 3) = flux * vf; + fx[nbx](i, j, k, 4) = flux * wf; + // pressure gradient fluxes + gp_rho_face(i, j, k, 0, vof[nbx], gp[nbx], rho[nbx], uf, vf, wf); + fx[nbx](i, j, k, 5) = flux * uf; + fx[nbx](i, j, k, 6) = flux * vf; + fx[nbx](i, j, k, 7) = flux * wf; + // Turn "on" all flux faces, later modified in + // process_fluxes_calc_src + fx[nbx](i, j, k, 8) = 1.0; + }); + amrex::ParallelFor( + mf_fy, mf_fy.n_grow, + [=] AMREX_GPU_DEVICE(int nbx, int i, int j, int k) noexcept { + amrex::Real flux = Gamma * alpha_flux( + i, j, k, 1, margin, vof[nbx], + tg_vof[nbx], norm[nbx]); + fy[nbx](i, j, k, 0) = flux; + flux *= (rho1 - rho2); + fy[nbx](i, j, k, 1) = flux; + amrex::Real uf, vf, wf; + velocity_face(i, j, k, 1, vof[nbx], vel[nbx], uf, vf, wf); + fy[nbx](i, j, k, 2) = flux * uf; + fy[nbx](i, j, k, 3) = flux * vf; + fy[nbx](i, j, k, 4) = flux * wf; + gp_rho_face(i, j, k, 1, vof[nbx], gp[nbx], rho[nbx], uf, vf, wf); + fy[nbx](i, j, k, 5) = flux * uf; + fy[nbx](i, j, k, 6) = flux * vf; + fy[nbx](i, j, k, 7) = flux * wf; + fy[nbx](i, j, k, 8) = 1.0; + }); + amrex::ParallelFor( + mf_fz, mf_fz.n_grow, + [=] AMREX_GPU_DEVICE(int nbx, int i, int j, int k) noexcept { + amrex::Real flux = Gamma * alpha_flux( + i, j, k, 2, margin, vof[nbx], + tg_vof[nbx], norm[nbx]); + fz[nbx](i, j, k, 0) = flux; + flux *= (rho1 - rho2); + fz[nbx](i, j, k, 1) = flux; + amrex::Real uf, vf, wf; + velocity_face(i, j, k, 2, vof[nbx], vel[nbx], uf, vf, wf); + fz[nbx](i, j, k, 2) = flux * uf; + fz[nbx](i, j, k, 3) = flux * vf; + fz[nbx](i, j, k, 4) = flux * wf; + gp_rho_face(i, j, k, 2, vof[nbx], gp[nbx], rho[nbx], uf, vf, wf); + fz[nbx](i, j, k, 5) = flux * uf; + fz[nbx](i, j, k, 6) = flux * vf; + fz[nbx](i, j, k, 7) = flux * wf; + fz[nbx](i, j, k, 8) = 1.0; + }); +} + +// Process reinitialization fluxes - zero non-internal to overset region; +// also calculate pressure source / sink term as a function of fluxes +void process_fluxes_calc_src( + amrex::MultiFab& mf_fx, + amrex::MultiFab& mf_fy, + amrex::MultiFab& mf_fz, + amrex::MultiFab& mf_psource, + const amrex::MultiFab& mf_vof, + const amrex::iMultiFab& mf_iblank) +{ + 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(); + constexpr amrex::Real tiny = std::numeric_limits::epsilon(); + // Zero fluxes based on iblank array + amrex::ParallelFor( + mf_fx, mf_fx.n_grow, mf_fx.n_comp, + [=] AMREX_GPU_DEVICE(int nbx, int i, int j, int k, int n) noexcept { + const bool zero_all = + (iblank[nbx](i - 1, j, k) + iblank[nbx](i, j, k) > -2); + fx[nbx](i, j, k, n) *= zero_all ? 0. : 1.; + }); + amrex::ParallelFor( + mf_fy, mf_fy.n_grow, mf_fy.n_comp, + [=] AMREX_GPU_DEVICE(int nbx, int i, int j, int k, int n) noexcept { + const bool zero_all = + (iblank[nbx](i, j - 1, k) + iblank[nbx](i, j, k) > -2); + fy[nbx](i, j, k, n) *= zero_all ? 0. : 1.; + }); + amrex::ParallelFor( + mf_fz, mf_fz.n_grow, mf_fz.n_comp, + [=] AMREX_GPU_DEVICE(int nbx, int i, int j, int k, int n) noexcept { + const bool zero_all = + (iblank[nbx](i, j, k - 1) + iblank[nbx](i, j, k) > -2); + fz[nbx](i, j, k, n) *= zero_all ? 0. : 1.; + }); + // With knowledge of fluxes, compute pressure source term + amrex::ParallelFor( + 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); + }); +} + +amrex::Real calculate_pseudo_velocity_scale( + const amrex::iMultiFab& mf_iblank, + const amrex::GpuArray dx, + const amrex::Real pvmax) +{ + // The dx of this level should be considered if iblank = -1 here + // Otherwise, set to max possible value for this mesh (level 0 values) + return (mf_iblank.min(0) == -1) ? std::min(std::min(dx[0], dx[1]), dx[2]) + : pvmax; +} + +// Calculate a type of CFL by measuring how much % VOF is being removed per cell +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::Real tol) +{ + // 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) / 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_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) / 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_fab = amrex::min(pdt_fab, pdt_lim); + }); + return pdt_fab; + }); + const amrex::Real pdt_fz = amrex::ReduceMin( + mf_fz, mf_vof, 0, + [=] AMREX_GPU_HOST_DEVICE( + amrex::Box const& bx, amrex::Array4 const& fz, + 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 (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); + } 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_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; +} + +// Apply reinitialization fluxes to modify fields +void apply_fluxes( + const amrex::MultiFab& mf_fx, + const amrex::MultiFab& mf_fy, + const amrex::MultiFab& mf_fz, + const amrex::MultiFab& mf_psource, + amrex::MultiFab& mf_vof, + amrex::MultiFab& mf_dens, + amrex::MultiFab& mf_vel, + amrex::MultiFab& mf_gp, + amrex::MultiFab& mf_pressure, + const amrex::GpuArray dx, + const amrex::Real ptfac, + const amrex::Real vof_tol) +{ + const auto& fx = mf_fx.const_arrays(); + const auto& fy = mf_fy.const_arrays(); + const auto& fz = mf_fz.const_arrays(); + const auto& sp = mf_psource.const_arrays(); + const auto& vof = mf_vof.arrays(); + const auto& dens = mf_dens.arrays(); + const auto& vel = mf_vel.arrays(); + const auto& gp = mf_gp.arrays(); + const auto& p = mf_pressure.arrays(); + + amrex::ParallelFor( + mf_vof, [=] AMREX_GPU_DEVICE(int nbx, int i, int j, int k) noexcept { + const amrex::Real olddens = dens[nbx](i, j, k); + vof[nbx](i, j, k) += + ptfac * + ((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]); + dens[nbx](i, j, k) += + ptfac * + ((fx[nbx](i + 1, j, k, 1) - fx[nbx](i, j, k, 1)) / dx[0] + + (fy[nbx](i, j + 1, k, 1) - fy[nbx](i, j, k, 1)) / dx[1] + + (fz[nbx](i, j, k + 1, 1) - fz[nbx](i, j, k, 1)) / dx[2]); + vel[nbx](i, j, k, 0) = + 1.0 / dens[nbx](i, j, k) * + (olddens * vel[nbx](i, j, k, 0) + + ptfac * + ((fx[nbx](i + 1, j, k, 2) - fx[nbx](i, j, k, 2)) / dx[0] + + (fy[nbx](i, j + 1, k, 2) - fy[nbx](i, j, k, 2)) / dx[1] + + (fz[nbx](i, j, k + 1, 2) - fz[nbx](i, j, k, 2)) / dx[2])); + vel[nbx](i, j, k, 1) = + 1.0 / dens[nbx](i, j, k) * + (olddens * vel[nbx](i, j, k, 1) + + ptfac * + ((fx[nbx](i + 1, j, k, 3) - fx[nbx](i, j, k, 3)) / dx[0] + + (fy[nbx](i, j + 1, k, 3) - fy[nbx](i, j, k, 3)) / dx[1] + + (fz[nbx](i, j, k + 1, 3) - fz[nbx](i, j, k, 3)) / dx[2])); + vel[nbx](i, j, k, 2) = + 1.0 / dens[nbx](i, j, k) * + (olddens * vel[nbx](i, j, k, 2) + + ptfac * + ((fx[nbx](i + 1, j, k, 4) - fx[nbx](i, j, k, 4)) / dx[0] + + (fy[nbx](i, j + 1, k, 4) - fy[nbx](i, j, k, 4)) / dx[1] + + (fz[nbx](i, j, k + 1, 4) - fz[nbx](i, j, k, 4)) / dx[2])); + gp[nbx](i, j, k, 0) += + ptfac * + ((fx[nbx](i + 1, j, k, 5) - fx[nbx](i, j, k, 5)) / dx[0] + + (fy[nbx](i, j + 1, k, 5) - fy[nbx](i, j, k, 5)) / dx[1] + + (fz[nbx](i, j, k + 1, 5) - fz[nbx](i, j, k, 5)) / dx[2]); + gp[nbx](i, j, k, 1) += + ptfac * + ((fx[nbx](i + 1, j, k, 6) - fx[nbx](i, j, k, 6)) / dx[0] + + (fy[nbx](i, j + 1, k, 6) - fy[nbx](i, j, k, 6)) / dx[1] + + (fz[nbx](i, j, k + 1, 6) - fz[nbx](i, j, k, 6)) / dx[2]); + gp[nbx](i, j, k, 2) += + ptfac * + ((fx[nbx](i + 1, j, k, 7) - fx[nbx](i, j, k, 7)) / dx[0] + + (fy[nbx](i, j + 1, k, 7) - fy[nbx](i, j, k, 7)) / dx[1] + + (fz[nbx](i, j, k + 1, 7) - fz[nbx](i, j, k, 7)) / dx[2]); + + // Ensure vof is bounded + vof[nbx](i, j, k) = + vof[nbx](i, j, k) < vof_tol + ? 0.0 + : (vof[nbx](i, j, k) > 1. - vof_tol ? 1. + : vof[nbx](i, j, k)); + // Density bounds are enforced elsewhere + }); + amrex::ParallelFor( + mf_pressure, + [=] AMREX_GPU_DEVICE(int nbx, int i, int j, int k) noexcept { + p[nbx](i, j, k) += ptfac * sp[nbx](i, j, k); + }); +} + +// Get the size of the smallest VOF flux to quantify convergence +amrex::Real measure_convergence( + amrex::MultiFab& mf_fx, amrex::MultiFab& mf_fy, amrex::MultiFab& mf_fz) +{ + // Get the maximum flux magnitude, but just for vof fluxes + const amrex::Real err_fx = amrex::ReduceMax( + mf_fx, 0, + [=] AMREX_GPU_HOST_DEVICE( + amrex::Box const& bx, + amrex::Array4 const& fx) -> amrex::Real { + amrex::Real err_fab = -1.0; + amrex::Loop(bx, [=, &err_fab](int i, int j, int k) noexcept { + err_fab = amrex::max(err_fab, std::abs(fx(i, j, k, 0))); + }); + return err_fab; + }); + const amrex::Real err_fy = amrex::ReduceMax( + mf_fy, 0, + [=] AMREX_GPU_HOST_DEVICE( + amrex::Box const& bx, + amrex::Array4 const& fy) -> amrex::Real { + amrex::Real err_fab = -1.0; + amrex::Loop(bx, [=, &err_fab](int i, int j, int k) noexcept { + err_fab = amrex::max(err_fab, std::abs(fy(i, j, k, 0))); + }); + return err_fab; + }); + const amrex::Real err_fz = amrex::ReduceMax( + mf_fz, 0, + [=] AMREX_GPU_HOST_DEVICE( + amrex::Box const& bx, + amrex::Array4 const& fz) -> amrex::Real { + amrex::Real err_fab = -1.0; + amrex::Loop(bx, [=, &err_fab](int i, int j, int k) noexcept { + err_fab = amrex::max(err_fab, std::abs(fz(i, j, k, 0))); + }); + return err_fab; + }); + const amrex::Real err = amrex::max(err_fx, amrex::max(err_fy, err_fz)); + return err; +} + +// Set levelset field to another quantity to view in plotfile for debugging +void equate_field(amrex::MultiFab& mf_dest, const amrex::MultiFab& mf_src) +{ + 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)); + }); +} + +// Replace pressure gradient with hydrostatic field in overset regions +void replace_gradp_hydrostatic( + amrex::MultiFab& mf_gp, + const amrex::MultiFab& mf_density, + const amrex::MultiFab& mf_refdens, + const amrex::iMultiFab& mf_iblank, + const amrex::Real grav_z, + const bool is_pptb) +{ + const auto& gp = mf_gp.arrays(); + const auto& rho = mf_density.const_arrays(); + const auto& rho0 = mf_refdens.const_arrays(); + const auto& iblank = mf_iblank.const_arrays(); + amrex::ParallelFor( + mf_gp, mf_gp.n_grow, + [=] AMREX_GPU_DEVICE(int nbx, int i, int j, int k) noexcept { + if (iblank[nbx](i, j, k) == -1) { + const amrex::Real dfac = + is_pptb ? rho[nbx](i, j, k) - rho0[nbx](i, j, k) + : rho[nbx](i, j, k); + gp[nbx](i, j, k, 0) = 0.; + gp[nbx](i, j, k, 1) = 0.; + gp[nbx](i, j, k, 2) = dfac * grav_z; + } + }); +} + +// Swap pressure gradient values in overset region +void replace_gradp( + amrex::MultiFab& mf_gp, + const amrex::MultiFab& mf_gp0, + const amrex::iMultiFab& mf_iblank) +{ + const auto gp = mf_gp.arrays(); + const auto gp0 = mf_gp0.const_arrays(); + const auto iblank = mf_iblank.const_arrays(); + amrex::ParallelFor( + mf_gp, mf_gp.n_grow, + [=] AMREX_GPU_DEVICE(int nbx, int i, int j, int k) noexcept { + if (iblank[nbx](i, j, k) == -1) { + gp[nbx](i, j, k, 0) = gp0[nbx](i, j, k, 0); + gp[nbx](i, j, k, 1) = gp0[nbx](i, j, k, 1); + gp[nbx](i, j, k, 2) = gp0[nbx](i, j, k, 2); + } + }); +} + +// Apply pressure gradient to velocity field +void apply_pressure_gradient( + amrex::MultiFab& mf_vel, + const amrex::MultiFab& mf_density, + const amrex::MultiFab& mf_gp, + const amrex::Real scaling_factor) +{ + const auto& vel = mf_vel.arrays(); + const auto& rho = mf_density.const_arrays(); + const auto& gp = mf_gp.const_arrays(); + amrex::ParallelFor( + mf_vel, [=] AMREX_GPU_DEVICE(int nbx, int i, int j, int k) noexcept { + const amrex::Real soverrho = scaling_factor / rho[nbx](i, j, k); + vel[nbx](i, j, k, 0) -= gp[nbx](i, j, k, 0) * soverrho; + vel[nbx](i, j, k, 1) -= gp[nbx](i, j, k, 1) * soverrho; + vel[nbx](i, j, k, 2) -= gp[nbx](i, j, k, 2) * soverrho; + }); +} + +} // namespace amr_wind::overset_ops \ No newline at end of file diff --git a/amr-wind/physics/multiphase/MultiPhase.H b/amr-wind/physics/multiphase/MultiPhase.H index 578884d87b..0dcde410dc 100644 --- a/amr-wind/physics/multiphase/MultiPhase.H +++ b/amr-wind/physics/multiphase/MultiPhase.H @@ -71,6 +71,13 @@ public: amrex::Real water_level() const { return m_water_level0; } + bool perturb_pressure() const { return m_use_perturb_pressure; } + + bool reconstruct_true_pressure() const + { + return m_reconstruct_true_pressure; + } + private: const CFDSim& m_sim; diff --git a/unit_tests/multiphase/test_vof_overset_ops.cpp b/unit_tests/multiphase/test_vof_overset_ops.cpp index 82b80083fe..fa27e9d538 100644 --- a/unit_tests/multiphase/test_vof_overset_ops.cpp +++ b/unit_tests/multiphase/test_vof_overset_ops.cpp @@ -1,7 +1,7 @@ #include "aw_test_utils/MeshTest.H" #include "aw_test_utils/iter_tools.H" #include "aw_test_utils/test_utils.H" -#include "amr-wind/overset/overset_ops_K.H" +#include "amr-wind/overset/overset_ops_routines.H" namespace amr_wind_tests { @@ -141,6 +141,51 @@ void init_velocity_etc( }); } +void init_gp_rho_etc( + amr_wind::Field& gp, + amr_wind::Field& rho, + amr_wind::Field& vof, + const int& dir, + const amrex::Real& rho1, + const amrex::Real& rho2) +{ + run_algorithm(gp, [&](const int lev, const amrex::MFIter& mfi) { + auto gp_arr = gp(lev).array(mfi); + auto rho_arr = rho(lev).array(mfi); + auto vof_arr = vof(lev).array(mfi); + const auto& bx = mfi.validbox(); + amrex::ParallelFor( + grow(bx, 2), [=] AMREX_GPU_DEVICE(int i, int j, int k) { + gp_arr(i, j, k) = 0.0; + vof_arr(i, j, k) = 0.0; + const int idx = (dir == 0 ? i : (dir == 1 ? j : k)); + if (idx == -1) { + vof_arr(i, j, k) = 0.41; + gp_arr(i, j, k, 0) = 1.0; + gp_arr(i, j, k, 1) = 2.0; + gp_arr(i, j, k, 2) = 3.0; + } else if (idx == 0) { + vof_arr(i, j, k) = 0.55; + gp_arr(i, j, k, 0) = 1.5; + gp_arr(i, j, k, 1) = 2.5; + gp_arr(i, j, k, 2) = 3.5; + } else if (idx == 1) { + vof_arr(i, j, k) = 0.2; + gp_arr(i, j, k, 0) = 2.0; + gp_arr(i, j, k, 1) = 3.0; + gp_arr(i, j, k, 2) = 4.0; + } else if (idx == 2) { + vof_arr(i, j, k) = 0.2; + gp_arr(i, j, k, 0) = 2.5; + gp_arr(i, j, k, 1) = 3.5; + gp_arr(i, j, k, 2) = 4.5; + } + rho_arr(i, j, k) = + rho1 * vof_arr(i, j, k) + rho2 * (1. - vof_arr(i, j, k)); + }); + }); +} + void calc_alpha_flux( amr_wind::Field& flux, amr_wind::Field& vof, @@ -184,6 +229,157 @@ void calc_velocity_face( }); } +void calc_gp_rho_face( + amr_wind::Field& flux, + amr_wind::Field& gp, + amr_wind::Field& rho, + amr_wind::Field& vof, + const int& dir, + const int& ioff = 0, + bool random_fflag = false) +{ + run_algorithm(flux, [&](const int lev, const amrex::MFIter& mfi) { + auto f_arr = flux(lev).array(mfi); + auto vof_arr = vof(lev).array(mfi); + auto gp_arr = gp(lev).array(mfi); + auto rho_arr = rho(lev).array(mfi); + const auto& bx = mfi.validbox(); + amrex::ParallelFor( + grow(bx, flux.num_grow()), + [=] AMREX_GPU_DEVICE(int i, int j, int k) { + amrex::Real u_f, v_f, w_f; + amrex::Real fac = 1.0; + if (random_fflag) { + fac = ((i + j + k) % 8 == 0) ? 0. : fac; + f_arr(i, j, k, 3 + ioff) = fac; + } + amr_wind::overset_ops::gp_rho_face( + i, j, k, dir, vof_arr, gp_arr, rho_arr, u_f, v_f, w_f); + f_arr(i, j, k, 0 + ioff) = u_f * fac; + f_arr(i, j, k, 1 + ioff) = v_f * fac; + f_arr(i, j, k, 2 + ioff) = w_f * fac; + }); + }); +} + +void calc_psrc( + amr_wind::Field& f_fx, + amr_wind::Field& f_fy, + amr_wind::Field& f_fz, + amr_wind::Field& f_vof, + amr_wind::Field& f_psrc, + amr_wind::Field& f_psrc_manual) +{ + constexpr amrex::Real tiny = std::numeric_limits::epsilon(); + run_algorithm(f_vof, [&](const int lev, const amrex::MFIter& mfi) { + auto fx = f_fx(lev).const_array(mfi); + auto fy = f_fy(lev).const_array(mfi); + auto fz = f_fz(lev).const_array(mfi); + auto vof = f_vof(lev).const_array(mfi); + auto psrc = f_psrc(lev).array(mfi); + auto psmn = f_psrc_manual(lev).array(mfi); + const auto& bx = mfi.validbox(); + amrex::ParallelFor(bx, [=] AMREX_GPU_DEVICE(int i, int j, int k) { + psrc(i, j, k) = amr_wind::overset_ops::gp_flux_tensor( + i, j, k, fx, fy, fz, tiny) && + amr_wind::overset_ops::normal_reinit_tensor( + i, j, k, fx, fy, fz, vof, tiny); + + const amrex::Real fxgpx_flux = + (fx(i, j, k, 5) + fx(i, j, k - 1, 5) + fx(i, j - 1, k, 5) + + fx(i, j - 1, k - 1, 5)) / + (fx(i, j, k, 8) + fx(i, j, k - 1, 8) + fx(i, j - 1, k, 8) + + fx(i, j - 1, k - 1, 8) + tiny); + const amrex::Real fxgpy_flux = + (fx(i, j, k, 6) + fx(i, j, k - 1, 6) + fx(i, j - 1, k, 6) + + fx(i, j - 1, k - 1, 6)) / + (fx(i, j, k, 8) + fx(i, j, k - 1, 8) + fx(i, j - 1, k, 8) + + fx(i, j - 1, k - 1, 8) + tiny); + const amrex::Real fxgpz_flux = + (fx(i, j, k, 7) + fx(i, j, k - 1, 7) + fx(i, j - 1, k, 7) + + fx(i, j - 1, k - 1, 7)) / + (fx(i, j, k, 8) + fx(i, j, k - 1, 8) + fx(i, j - 1, k, 8) + + fx(i, j - 1, k - 1, 8) + tiny); + + const amrex::Real fygpx_flux = + (fy(i, j, k, 5) + fy(i - 1, j, k, 5) + fy(i, j, k - 1, 5) + + fy(i - 1, j, k - 1, 5)) / + (fy(i, j, k, 8) + fy(i - 1, j, k, 8) + fy(i, j, k - 1, 8) + + fy(i - 1, j, k - 1, 8) + tiny); + const amrex::Real fygpy_flux = + (fy(i, j, k, 6) + fy(i - 1, j, k, 6) + fy(i, j, k - 1, 6) + + fy(i - 1, j, k - 1, 6)) / + (fy(i, j, k, 8) + fy(i - 1, j, k, 8) + fy(i, j, k - 1, 8) + + fy(i - 1, j, k - 1, 8) + tiny); + const amrex::Real fygpz_flux = + (fy(i, j, k, 7) + fy(i - 1, j, k, 7) + fy(i, j, k - 1, 7) + + fy(i - 1, j, k - 1, 7)) / + (fy(i, j, k, 8) + fy(i - 1, j, k, 8) + fy(i, j, k - 1, 8) + + fy(i - 1, j, k - 1, 8) + tiny); + + const amrex::Real fzgpx_flux = + (fz(i, j, k, 5) + fz(i - 1, j, k, 5) + fz(i, j - 1, k, 5) + + fz(i - 1, j - 1, k, 5)) / + (fz(i, j, k, 8) + fz(i - 1, j, k, 8) + fz(i, j - 1, k, 8) + + fz(i - 1, j - 1, k, 8) + tiny); + const amrex::Real fzgpy_flux = + (fz(i, j, k, 6) + fz(i - 1, j, k, 6) + fz(i, j - 1, k, 6) + + fz(i - 1, j - 1, k, 6)) / + (fz(i, j, k, 8) + fz(i - 1, j, k, 8) + fz(i, j - 1, k, 8) + + fz(i - 1, j - 1, k, 8) + tiny); + const amrex::Real fzgpz_flux = + (fz(i, j, k, 7) + fz(i - 1, j, k, 7) + fz(i, j - 1, k, 7) + + fz(i - 1, j - 1, k, 7)) / + (fz(i, j, k, 8) + fz(i - 1, j, k, 8) + fz(i, j - 1, k, 8) + + fz(i - 1, j - 1, k, 8) + tiny); + + // Calculate interface normal at pressure node + amrex::Real n0 = + ((vof(i, j, k) - vof(i - 1, j, k)) * fx(i, j, k, 8) + + (vof(i, j - 1, k) - vof(i - 1, j - 1, k)) * + fx(i, j - 1, k, 8) + + (vof(i, j, k - 1) - vof(i - 1, j, k - 1)) * + fx(i, j, k - 1, 8) + + (vof(i, j - 1, k - 1) - vof(i - 1, j - 1, k - 1)) * + fx(i, j - 1, k - 1, 8)) / + (fx(i, j, k, 8) + fx(i, j, k - 1, 8) + fx(i, j - 1, k, 8) + + fx(i, j - 1, k - 1, 8) + tiny); + amrex::Real n1 = + ((vof(i, j, k) - vof(i, j - 1, k)) * fy(i, j, k, 8) + + (vof(i - 1, j, k) - vof(i - 1, j - 1, k)) * + fy(i - 1, j, k, 8) + + (vof(i, j, k - 1) - vof(i, j - 1, k - 1)) * + fy(i, j, k - 1, 8) + + (vof(i - 1, j, k - 1) - vof(i - 1, j - 1, k - 1)) * + fy(i - 1, j, k - 1, 8)) / + (fy(i, j, k, 8) + fy(i, j, k - 1, 8) + fy(i - 1, j, k, 8) + + fy(i - 1, j, k - 1, 8) + tiny); + amrex::Real n2 = + ((vof(i, j, k) - vof(i, j, k - 1)) * fz(i, j, k, 8) + + (vof(i - 1, j, k) - vof(i - 1, j, k - 1)) * + fz(i - 1, j, k, 8) + + (vof(i, j - 1, k) - vof(i, j - 1, k - 1)) * + fz(i, j - 1, k, 8) + + (vof(i - 1, j - 1, k) - vof(i - 1, j - 1, k - 1)) * + fz(i - 1, j - 1, k, 8)) / + (fz(i, j, k, 8) + fz(i, j - 1, k, 8) + fz(i - 1, j, k, 8) + + fz(i - 1, j - 1, k, 8) + tiny); + // Normalize vector + amrex::Real nmag = std::sqrt(n0 * n0 + n1 * n1 + n2 * n2) + tiny; + n0 /= nmag; + n1 /= nmag; + n2 /= nmag; + + // Perform double contraction + psmn(i, j, k) = fxgpx_flux * n0 * n0 + fxgpy_flux * n0 * n1 + + fxgpz_flux * n0 * n2 + fygpx_flux * n1 * n0 + + fygpy_flux * n1 * n1 + fygpz_flux * n1 * n2 + + fzgpx_flux * n2 * n0 + fzgpy_flux * n2 * n1 + + fzgpz_flux * n2 * n2; + }); + }); +} + amrex::Real check_alpha_flux_impl(amr_wind::Field& flux, const int& dir) { amrex::Real error_total = 0; @@ -295,6 +491,92 @@ amrex::Real check_velocity_face_impl(amr_wind::Field& flux, const int& dir) } return error_total; } + +amrex::Real check_gp_rho_face_impl( + amr_wind::Field& flux, + const int& dir, + const amrex::Real& rho1, + const amrex::Real& rho2) +{ + amrex::Real error_total = 0; + + for (int lev = 0; lev < flux.repo().num_active_levels(); ++lev) { + + error_total += amrex::ReduceSum( + flux(lev), 0, + [=] AMREX_GPU_HOST_DEVICE( + amrex::Box const& bx, + amrex::Array4 const& f_arr) -> amrex::Real { + amrex::Real error = 0; + + amrex::Loop(bx, [=, &error](int i, int j, int k) noexcept { + const int idx = (dir == 0 ? i : (dir == 1 ? j : k)); + if (idx == 0) { + // gphi > 0, uwpind from the "left" + const amrex::Real rho_answer = + rho1 * 0.41 + rho2 * (1. - 0.41); + amrex::Real flux_answer = 1.0 / rho_answer; + error += std::abs(f_arr(i, j, k, 0) - flux_answer); + flux_answer = 2.0 / rho_answer; + error += std::abs(f_arr(i, j, k, 1) - flux_answer); + flux_answer = 3.0 / rho_answer; + error += std::abs(f_arr(i, j, k, 2) - flux_answer); + } else if (idx == 1) { + // gphi < 0, upwind from the "right" + const amrex::Real rho_answer = + rho1 * 0.2 + rho2 * (1. - 0.2); + amrex::Real flux_answer = 2.0 / rho_answer; + error += std::abs(f_arr(i, j, k, 0) - flux_answer); + flux_answer = 3.0 / rho_answer; + error += std::abs(f_arr(i, j, k, 1) - flux_answer); + flux_answer = 4.0 / rho_answer; + error += std::abs(f_arr(i, j, k, 2) - flux_answer); + } else if (idx == 2) { + // gphi = 0, average both sides + const amrex::Real rho_r = + rho1 * 0.2 + rho2 * (1. - 0.2); + const amrex::Real rho_l = + rho1 * 0.2 + rho2 * (1. - 0.2); + amrex::Real flux_answer = + 0.5 * (2.5 / rho_r + 2.0 / rho_l); + error += std::abs(f_arr(i, j, k, 0) - flux_answer); + flux_answer = 0.5 * (3.5 / rho_r + 3.0 / rho_l); + error += std::abs(f_arr(i, j, k, 1) - flux_answer); + flux_answer = 0.5 * (4.5 / rho_r + 4.0 / rho_l); + error += std::abs(f_arr(i, j, k, 2) - flux_answer); + } + }); + + return error; + }); + } + return error_total; +} + +amrex::Real +check_psrc_manual_impl(amr_wind::Field& psrc, amr_wind::Field& psrc_manual) +{ + amrex::Real error_total = 0; + + for (int lev = 0; lev < psrc.repo().num_active_levels(); ++lev) { + + error_total += amrex::ReduceSum( + psrc(lev), psrc_manual(lev), 0, + [=] AMREX_GPU_HOST_DEVICE( + amrex::Box const& bx, + amrex::Array4 const& f_arr, + amrex::Array4 const& fm_arr) -> amrex::Real { + amrex::Real error = 0; + + amrex::Loop(bx, [=, &error](int i, int j, int k) noexcept { + error += std::abs(f_arr(i, j, k) - fm_arr(i, j, k)); + }); + + return error; + }); + } + return error_total; +} } // namespace TEST_F(VOFOversetOps, alpha_flux) @@ -403,4 +685,206 @@ TEST_F(VOFOversetOps, velocity_face) EXPECT_NEAR(error_total, 0.0, 1e-15); } +TEST_F(VOFOversetOps, gp_rho_face) +{ + populate_parameters(); + initialize_mesh(); + + const amrex::Real rho_liq = 1000.; + const amrex::Real rho_gas = 1.; + + auto& repo = sim().repo(); + const int nghost = 3; + auto& gp = repo.declare_field("gp", 3, nghost); + auto& rho = repo.declare_field("density", 1, nghost); + auto& vof = repo.declare_field("vof", 1, nghost); + + // Create flux fields + auto& flux_x = + repo.declare_field("flux_x", 3, 0, 1, amr_wind::FieldLoc::XFACE); + auto& flux_y = + repo.declare_field("flux_y", 3, 0, 1, amr_wind::FieldLoc::YFACE); + auto& flux_z = + repo.declare_field("flux_z", 3, 0, 1, amr_wind::FieldLoc::ZFACE); + + // -- Variations in x direction -- // + int dir = 0; + // Initialize gp, rho, and vof + init_gp_rho_etc(gp, rho, vof, dir, rho_liq, rho_gas); + // Populate flux field + calc_gp_rho_face(flux_x, gp, rho, vof, dir); + // Check results + amrex::Real error_total = + check_gp_rho_face_impl(flux_x, dir, rho_liq, rho_gas); + amrex::ParallelDescriptor::ReduceRealSum(error_total); + EXPECT_NEAR(error_total, 0.0, 1e-15); + + // -- Variations in y direction -- // + dir = 1; + // Initialize gp, rho, and vof + init_gp_rho_etc(gp, rho, vof, dir, rho_liq, rho_gas); + // Populate flux field + calc_gp_rho_face(flux_y, gp, rho, vof, dir); + // Check results + error_total = check_gp_rho_face_impl(flux_y, dir, rho_liq, rho_gas); + amrex::ParallelDescriptor::ReduceRealSum(error_total); + EXPECT_NEAR(error_total, 0.0, 1e-15); + + // -- Variations in z direction -- // + dir = 2; + // Initialize gp, rho, and vof + init_gp_rho_etc(gp, rho, vof, dir, rho_liq, rho_gas); + // Populate flux field + calc_gp_rho_face(flux_z, gp, rho, vof, dir); + // Check results + error_total = check_gp_rho_face_impl(flux_z, dir, rho_liq, rho_gas); + amrex::ParallelDescriptor::ReduceRealSum(error_total); + EXPECT_NEAR(error_total, 0.0, 1e-15); +} + +TEST_F(VOFOversetOps, pseudo_vscale_dt) +{ + populate_parameters(); + initialize_mesh(); + + auto& repo = sim().repo(); + const int nghost = 3; + auto& vof = repo.declare_field("vof", 1, nghost); + 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); + 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 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]); + + // Create flux fields + auto& flux_x = + repo.declare_field("flux_x", 1, 0, 1, amr_wind::FieldLoc::XFACE); + auto& flux_y = + repo.declare_field("flux_y", 1, 0, 1, amr_wind::FieldLoc::YFACE); + auto& flux_z = + repo.declare_field("flux_z", 1, 0, 1, amr_wind::FieldLoc::ZFACE); + + // -- Variations in x direction -- // + int dir = 0; + // Initialize vof and other fields + init_vof_etc(vof, tg_vof, norm, dir); + // Populate flux field + calc_alpha_flux(flux_x, vof, tg_vof, norm, dir, margin); + // Calculate pseudo dt + amrex::Real ptfac = 100.0; + 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); + ptfac = amrex::min(ptfac, ptfac_lev); + } + amrex::ParallelDescriptor::ReduceRealMin(ptfac); + EXPECT_DOUBLE_EQ(ptfac, pdt_answer); + // Zero flux for subsequent tests + flux_x.setVal(0.0); + + // -- Variations in y direction -- // + dir = 1; + // Initialize vof and other fields + init_vof_etc(vof, tg_vof, norm, dir); + // Populate flux field + calc_alpha_flux(flux_y, vof, tg_vof, norm, dir, margin); + // Calculate pseudo dt + ptfac = 100.0; + 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); + ptfac = amrex::min(ptfac, ptfac_lev); + } + amrex::ParallelDescriptor::ReduceRealMin(ptfac); + EXPECT_DOUBLE_EQ(ptfac, pdt_answer); + // Zero flux for subsequent test + flux_y.setVal(0.0); + + // -- Variations in z direction -- // + dir = 2; + // Initialize vof and other fields + init_vof_etc(vof, tg_vof, norm, dir); + // Populate flux field + calc_alpha_flux(flux_z, vof, tg_vof, norm, dir, margin); + // Calculate pseudo dt + ptfac = 100.0; + 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); + 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; + 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); +} + +TEST_F(VOFOversetOps, psource_manual) +{ + populate_parameters(); + initialize_mesh(); + + const amrex::Real rho_liq = 1000.; + const amrex::Real rho_gas = 1.; + + auto& repo = sim().repo(); + const int nghost = 3; + auto& gp = repo.declare_field("gp", 3, nghost); + auto& rho = repo.declare_field("density", 1, nghost); + auto& vof = repo.declare_field("vof", 1, nghost); + + auto& psrc = repo.declare_field("psrc", 1, 0, 1, amr_wind::FieldLoc::NODE); + auto& psmn = + repo.declare_field("psrc_manual", 1, 0, 1, amr_wind::FieldLoc::NODE); + + // Create flux fields + auto& flux_x = + repo.declare_field("flux_x", 9, 1, 1, amr_wind::FieldLoc::XFACE); + auto& flux_y = + repo.declare_field("flux_y", 9, 1, 1, amr_wind::FieldLoc::YFACE); + auto& flux_z = + repo.declare_field("flux_z", 9, 1, 1, amr_wind::FieldLoc::ZFACE); + + // Initialize gp_rho fluxes in every direction + int dir = 0; + init_gp_rho_etc(gp, rho, vof, dir, rho_liq, rho_gas); + calc_gp_rho_face(flux_x, gp, rho, vof, dir, 5, true); + dir = 1; + init_gp_rho_etc(gp, rho, vof, dir, rho_liq, rho_gas); + calc_gp_rho_face(flux_y, gp, rho, vof, dir, 5, true); + dir = 2; + init_gp_rho_etc(gp, rho, vof, dir, rho_liq, rho_gas); + calc_gp_rho_face(flux_z, gp, rho, vof, dir, 5, true); + + // Calculate psrc with vector/tensor ops and manually + calc_psrc(flux_x, flux_y, flux_z, vof, psrc, psmn); + + // Check difference + amrex::Real error_total = check_psrc_manual_impl(psrc, psmn); + amrex::ParallelDescriptor::ReduceRealSum(error_total); + EXPECT_NEAR(error_total, 0.0, 1e-10); +} + } // namespace amr_wind_tests \ No newline at end of file