From 453ec3c9bcb1efb9a9161a2b0ffcc1d3a2459e75 Mon Sep 17 00:00:00 2001 From: Michael Kuhn Date: Tue, 11 Jun 2024 15:25:22 -0600 Subject: [PATCH] changes to include pressure sharpening and treat dx differently --- amr-wind/overset/OversetOps.cpp | 72 +++++--- amr-wind/overset/overset_ops_K.H | 47 +++++ amr-wind/overset/overset_ops_routines.H | 226 +++++++++++++++++++++--- 3 files changed, 297 insertions(+), 48 deletions(-) diff --git a/amr-wind/overset/OversetOps.cpp b/amr-wind/overset/OversetOps.cpp index b6b0112679..834c9ad7a1 100644 --- a/amr-wind/overset/OversetOps.cpp +++ b/amr-wind/overset/OversetOps.cpp @@ -199,19 +199,31 @@ 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) + // 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); + // Create scratch field for pressure source term + auto p_src = repo.create_scratch_field(1, 0, amr_wind::FieldLoc::NODE); + // Create scratch fields for normal vector and target vof field 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 + // Create scratch field for pressure gradient that has ghost cells + // 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 + 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(); @@ -219,7 +231,19 @@ void OversetOps::sharpen_nalu_data() // Populate approximate signed distance function overset_ops::populate_psi(levelset(lev), vof(lev), i_th, m_asdf_tiny); + + // Populate gp scratch field + gp(lev).setVal(0.0); // for external boundaries + amrex::MultiFab::Copy(gp(lev), gp_noghost(lev), 0, 0, 3, 0); // nonghost + gp(lev).FillBoundary(geom[lev].periodicity()); // internal + + // 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::ParallelDescriptor::ReduceRealMin(pvscale); // Convert levelset to vof to get target_vof m_mphase->levelset2vof(iblank_cell, *target_vof); @@ -270,13 +294,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_margin, - m_mphase->rho1(), m_mphase->rho2()); + (*target_vof)(lev), (*normal_vec)(lev), velocity(lev), gp(lev), + rho(lev), pvscale, m_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), vof(lev), + (*p_src)(lev), iblank_cell(lev)); // Measure convergence to determine if loop can stop if (cconv) { @@ -312,15 +337,19 @@ void OversetOps::sharpen_nalu_data() // Apply fluxes for (int lev = 0; lev < nlevels; ++lev) { + 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); + + // Update ghost cells + vof(lev).FillBoundary(geom[lev].periodicity()); + velocity(lev).FillBoundary(geom[lev].periodicity()); + gp(lev).FillBoundary(geom[lev].periodicity()); } - // Fillpatch for ghost cells - vof.fillpatch((*m_sim_ptr).time().current_time()); - velocity.fillpatch((*m_sim_ptr).time().current_time()); - // Update density (fillpatch built in) m_mphase->set_density_via_vof(); @@ -335,6 +364,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 diff --git a/amr-wind/overset/overset_ops_K.H b/amr-wind/overset/overset_ops_K.H index 89b0780c7b..db6c1d23d7 100644 --- a/amr-wind/overset/overset_ops_K.H +++ b/amr-wind/overset/overset_ops_K.H @@ -114,6 +114,53 @@ void AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE velocity_face( } } +void AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE gp_rho_face( + int i, + int j, + int k, + 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 + int ii = i; + int jj = j; + int kk = k; + ii += (dir == 0) ? -1 : 0; + jj += (dir == 1) ? -1 : 0; + kk += (dir == 2) ? -1 : 0; + + // Gradient of phi normal to interface + const amrex::Real gphi = (vof(i, j, k) - vof(ii, jj, kk)); + + // Get velocities on both sides + const amrex::Real u_ = gp(i, j, k, 0) / rho(i, j, k); + const amrex::Real v_ = gp(i, j, k, 1) / rho(i, j, k); + const amrex::Real w_ = gp(i, j, k, 2) / rho(i, j, k); + const amrex::Real u_nb = gp(ii, jj, kk, 0) / rho(ii, jj, kk); + const amrex::Real v_nb = gp(ii, jj, kk, 1) / rho(ii, jj, kk); + const amrex::Real w_nb = gp(ii, jj, kk, 2) / rho(ii, jj, kk); + // 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_; + } +} + } // 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 bb47c32e42..dcff14963e 100644 --- a/amr-wind/overset/overset_ops_routines.H +++ b/amr-wind/overset/overset_ops_routines.H @@ -123,15 +123,18 @@ void populate_sharpen_fluxes( amrex::MultiFab& mf_target_vof, amrex::MultiFab& mf_norm, amrex::MultiFab& mf_velocity, + amrex::MultiFab& mf_gp, + amrex::MultiFab& mf_density, + const amrex::Real Gamma, const amrex::Real margin, const amrex::Real rho1, const amrex::Real rho2) { for (amrex::MFIter mfi(mf_vof); mfi.isValid(); ++mfi) { const auto& vbx = mfi.validbox(); - const auto& xbx = amrex::surroundingNodes(vbx, 0); - const auto& ybx = amrex::surroundingNodes(vbx, 1); - const auto& zbx = amrex::surroundingNodes(vbx, 2); + const auto& xbxp1 = grow(amrex::surroundingNodes(vbx, 0), 1); + const auto& ybxp1 = grow(amrex::surroundingNodes(vbx, 1), 1); + const auto& zbxp1 = grow(amrex::surroundingNodes(vbx, 2), 1); const amrex::Array4& fx = mf_fx.array(mfi); const amrex::Array4& fy = mf_fy.array(mfi); const amrex::Array4& fz = mf_fz.array(mfi); @@ -141,13 +144,16 @@ void populate_sharpen_fluxes( const amrex::Array4& norm = mf_norm.const_array(mfi); const amrex::Array4& vel = mf_velocity.const_array(mfi); + const amrex::Array4& gp = mf_gp.const_array(mfi); + const amrex::Array4& rho = + mf_density.const_array(mfi); // Populate vof and density fluxes for each direction amrex::ParallelFor( - xbx, ybx, zbx, + xbxp1, ybxp1, zbxp1, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept { // vof flux amrex::Real flux = - alpha_flux(i, j, k, 0, margin, vof, tg_vof, norm); + Gamma * alpha_flux(i, j, k, 0, margin, vof, tg_vof, norm); fx(i, j, k, 0) = flux; // density flux flux *= (rho1 - rho2); @@ -158,10 +164,17 @@ void populate_sharpen_fluxes( fx(i, j, k, 2) = flux * uf; fx(i, j, k, 3) = flux * vf; fx(i, j, k, 4) = flux * wf; + gp_rho_face(i, j, k, 0, vof, gp, rho, uf, vf, wf); + fx(i, j, k, 5) = flux * uf; + fx(i, j, k, 6) = flux * vf; + fx(i, j, k, 7) = flux * wf; + // Turn "on" all flux faces, later modified in + // process_fluxes_calc_src + fx(i, j, k, 8) = 1.0; }, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept { amrex::Real flux = - alpha_flux(i, j, k, 1, margin, vof, tg_vof, norm); + Gamma * alpha_flux(i, j, k, 1, margin, vof, tg_vof, norm); fy(i, j, k, 0) = flux; flux *= (rho1 - rho2); fy(i, j, k, 1) = flux; @@ -170,10 +183,15 @@ void populate_sharpen_fluxes( fy(i, j, k, 2) = flux * uf; fy(i, j, k, 3) = flux * vf; fy(i, j, k, 4) = flux * wf; + gp_rho_face(i, j, k, 1, vof, gp, rho, uf, vf, wf); + fy(i, j, k, 5) = flux * uf; + fy(i, j, k, 6) = flux * vf; + fy(i, j, k, 7) = flux * wf; + fy(i, j, k, 8) = 1.0; }, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept { amrex::Real flux = - alpha_flux(i, j, k, 2, margin, vof, tg_vof, norm); + Gamma * alpha_flux(i, j, k, 2, margin, vof, tg_vof, norm); fz(i, j, k, 0) = flux; flux *= (rho1 - rho2); fz(i, j, k, 1) = flux; @@ -182,15 +200,23 @@ void populate_sharpen_fluxes( fz(i, j, k, 2) = flux * uf; fz(i, j, k, 3) = flux * vf; fz(i, j, k, 4) = flux * wf; + gp_rho_face(i, j, k, 2, vof, gp, rho, uf, vf, wf); + fz(i, j, k, 5) = flux * uf; + fz(i, j, k, 6) = flux * vf; + fz(i, j, k, 7) = flux * wf; + fz(i, j, k, 8) = 1.0; }); } } -// 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, + amrex::MultiFab& mf_vof, + amrex::MultiFab& mf_psource, amrex::iMultiFab& mf_iblank) { int ncompx = mf_fx.nComp(); @@ -198,16 +224,18 @@ void process_fluxes( int ncompz = mf_fz.nComp(); for (amrex::MFIter mfi(mf_iblank); mfi.isValid(); ++mfi) { const auto& vbx = mfi.validbox(); - const auto& xbx = amrex::surroundingNodes(vbx, 0); - const auto& ybx = amrex::surroundingNodes(vbx, 1); - const auto& zbx = amrex::surroundingNodes(vbx, 2); + const auto& xbxp1 = grow(amrex::surroundingNodes(vbx, 0), 1); + const auto& ybxp1 = grow(amrex::surroundingNodes(vbx, 1), 1); + const auto& zbxp1 = grow(amrex::surroundingNodes(vbx, 2), 1); const amrex::Array4& fx = mf_fx.array(mfi); const amrex::Array4& fy = mf_fy.array(mfi); const amrex::Array4& fz = mf_fz.array(mfi); + const amrex::Array4& vof = mf_vof.array(mfi); + const amrex::Array4& sp = mf_psource.array(mfi); const amrex::Array4& iblank = mf_iblank.const_array(mfi); // Only faces with iblank = -1 on both sides can have nonzero flux amrex::ParallelFor( - xbx, ybx, zbx, + xbxp1, ybxp1, zbxp1, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept { bool zero_all = (iblank(i - 1, j, k) + iblank(i, j, k) > -2); for (int n = 0; n < ncompx; ++n) { @@ -226,9 +254,123 @@ void process_fluxes( fz(i, j, k, n) *= zero_all ? 0. : 1.; } }); + // With knowledge of fluxes, compute pressure source term + amrex::Box const& nbx = mfi.nodaltilebox(); + constexpr amrex::Real tiny = + std::numeric_limits::epsilon(); + amrex::ParallelFor( + nbx, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept { + // Get components of outer product between f and gradp/rho + // from interpolating between values at neighboring faces + 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 + sp(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 calculate_pseudo_velocity_scale( + amrex::iMultiFab& mf_iblank, + amrex::GpuArray dx, + const amrex::Real pvmax) +{ + // Get minimum length scale from dx + const amrex::Real dx_min = std::min(std::min(dx[0], dx[1]), dx[2]); + // The dx of this level should be considered if iblank = -1 here + // Otherwise, set to max possible value for this mesh (level 0 values) + const amrex::Real pvscale = (mf_iblank.min(0) == -1) ? dx_min : pvmax; + return pvscale; +} + // 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, @@ -304,49 +446,71 @@ void apply_fluxes( amrex::MultiFab& mf_fx, amrex::MultiFab& mf_fy, amrex::MultiFab& mf_fz, + amrex::MultiFab& mf_psource, amrex::MultiFab& mf_vof, amrex::MultiFab& mf_dens, amrex::MultiFab& mf_vel, + amrex::MultiFab& mf_gp, + amrex::MultiFab& mf_pressure, + amrex::GpuArray dx, amrex::Real ptfac, amrex::Real vof_tol) { for (amrex::MFIter mfi(mf_vof); mfi.isValid(); ++mfi) { const auto& vbx = mfi.validbox(); + const auto& nbx = mfi.nodaltilebox(); const amrex::Array4& fx = mf_fx.array(mfi); const amrex::Array4& fy = mf_fy.array(mfi); const amrex::Array4& fz = mf_fz.array(mfi); const amrex::Array4& vof = mf_vof.array(mfi); const amrex::Array4& dens = mf_dens.array(mfi); const amrex::Array4& vel = mf_vel.array(mfi); + const amrex::Array4& gp = mf_gp.array(mfi); + const amrex::Array4& p = mf_pressure.array(mfi); + const amrex::Array4& sp = mf_psource.array(mfi); amrex::ParallelFor( vbx, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept { const amrex::Real olddens = dens(i, j, k); - vof(i, j, k) += ptfac * ((fx(i + 1, j, k, 0) - fx(i, j, k, 0)) + - (fy(i, j + 1, k, 0) - fy(i, j, k, 0)) + - (fz(i, j, k + 1, 0) - fz(i, j, k, 0))); + vof(i, j, k) += + ptfac * ((fx(i + 1, j, k, 0) - fx(i, j, k, 0)) / dx[0] + + (fy(i, j + 1, k, 0) - fy(i, j, k, 0)) / dx[1] + + (fz(i, j, k + 1, 0) - fz(i, j, k, 0)) / dx[2]); dens(i, j, k) += - ptfac * ((fx(i + 1, j, k, 1) - fx(i, j, k, 1)) + - (fy(i, j + 1, k, 1) - fy(i, j, k, 1)) + - (fz(i, j, k + 1, 1) - fz(i, j, k, 1))); + ptfac * ((fx(i + 1, j, k, 1) - fx(i, j, k, 1)) / dx[0] + + (fy(i, j + 1, k, 1) - fy(i, j, k, 1)) / dx[1] + + (fz(i, j, k + 1, 1) - fz(i, j, k, 1)) / dx[2]); vel(i, j, k, 0) = 1.0 / dens(i, j, k) * (olddens * vel(i, j, k, 0) + - ptfac * ((fx(i + 1, j, k, 2) - fx(i, j, k, 2)) + - (fy(i, j + 1, k, 2) - fy(i, j, k, 2)) + - (fz(i, j, k + 1, 2) - fz(i, j, k, 2)))); + ptfac * ((fx(i + 1, j, k, 2) - fx(i, j, k, 2)) / dx[0] + + (fy(i, j + 1, k, 2) - fy(i, j, k, 2)) / dx[1] + + (fz(i, j, k + 1, 2) - fz(i, j, k, 2)) / dx[2])); vel(i, j, k, 1) = 1.0 / dens(i, j, k) * (olddens * vel(i, j, k, 1) + - ptfac * ((fx(i + 1, j, k, 3) - fx(i, j, k, 3)) + - (fy(i, j + 1, k, 3) - fy(i, j, k, 3)) + - (fz(i, j, k + 1, 3) - fz(i, j, k, 3)))); + ptfac * ((fx(i + 1, j, k, 3) - fx(i, j, k, 3)) / dx[0] + + (fy(i, j + 1, k, 3) - fy(i, j, k, 3)) / dx[1] + + (fz(i, j, k + 1, 3) - fz(i, j, k, 3)) / dx[2])); vel(i, j, k, 2) = 1.0 / dens(i, j, k) * (olddens * vel(i, j, k, 2) + - ptfac * ((fx(i + 1, j, k, 4) - fx(i, j, k, 4)) + - (fy(i, j + 1, k, 4) - fy(i, j, k, 4)) + - (fz(i, j, k + 1, 4) - fz(i, j, k, 4)))); + ptfac * ((fx(i + 1, j, k, 4) - fx(i, j, k, 4)) / dx[0] + + (fy(i, j + 1, k, 4) - fy(i, j, k, 4)) / dx[1] + + (fz(i, j, k + 1, 4) - fz(i, j, k, 4)) / dx[2])); + + gp(i, j, k, 0) += + ptfac * ((fx(i + 1, j, k, 5) - fx(i, j, k, 5)) / dx[0] + + (fy(i, j + 1, k, 5) - fy(i, j, k, 5)) / dx[1] + + (fz(i, j, k + 1, 5) - fz(i, j, k, 5)) / dx[2]); + gp(i, j, k, 1) += + ptfac * ((fx(i + 1, j, k, 6) - fx(i, j, k, 6)) / dx[0] + + (fy(i, j + 1, k, 6) - fy(i, j, k, 6)) / dx[1] + + (fz(i, j, k + 1, 6) - fz(i, j, k, 6)) / dx[2]); + gp(i, j, k, 2) += + ptfac * ((fx(i + 1, j, k, 7) - fx(i, j, k, 7)) / dx[0] + + (fy(i, j + 1, k, 7) - fy(i, j, k, 7)) / dx[1] + + (fz(i, j, k + 1, 7) - fz(i, j, k, 7)) / dx[2]); // Ensure vof is bounded vof(i, j, k) = @@ -355,6 +519,12 @@ void apply_fluxes( : (vof(i, j, k) > 1. - vof_tol ? 1. : vof(i, j, k)); // Density bounds are enforced elsewhere }); + + // Apply pressure source/sink + amrex::ParallelFor( + nbx, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept { + p(i, j, k) += ptfac * sp(i, j, k); + }); } }