From 9bf16a76ce93ad81968ea57e7a4b3c230fc8c7ab Mon Sep 17 00:00:00 2001 From: Ann Almgren Date: Sun, 27 Sep 2026 14:12:01 -0700 Subject: [PATCH] Fix issues #179, #181-#185, #188-#190 #179 NS_derive: the EB one-sided vorticity stencils in dermgvort read velocity two cells away without checking those cells are uncovered -- isConnected only certifies the immediate neighbor. Certify the intermediate cell before taking the quadratic form, drop to a two-point difference when only one cell is available, and leave the derivative at zero when both sides are disconnected (2D x/y and 3D x/y/z). #181 ViscBench: the loop over global grid indices dereferenced MultiFabs by global index and FillVar(FArrayBox*,...) fills on one rank only, so the tool was broken under MPI. Use the distributed FillVar(MultiFab&,...) overload, drive the exact-solution loop from MFIter, and print through amrex::Print() so the collective norm calls stay on every rank. #182 getForce: the getForceVerbose min/max loops indexed State and force by absolute state component, reading past the end of the velocity-only State FAB and the 0-indexed force FAB that scalar callers pass. Index each FAB by its own layout; the printed label keeps the absolute component. #183 main: the regrid-on-restart test lacked the max_step < 0 guard that the time loop below applies, so with only stop_time set it regridded and wrote a spurious plotfile/checkpoint before stepping anyway. #184 HIT fast force: the coarse hi bound was one node short whenever ihi mod ff_factor was in [1, ff_factor-2], so the interpolation read past ff_force. Use ihi/ff_factor + 1, which is exactly the node the kernel touches, and give ff_force an Elixir so its device memory is not freed while the asynchronous kernels are still reading it. #185 TurbulentForcing: nothing constrained the mode indices writing into the flat stack array tmp, so large nmodes or a high-aspect domain corrupted the stack; assert the mode counts fit array_size. The GPU copy out of tmp was asynchronous with no sync before the stack frame died; use the synchronous htod_memcpy. Also fixes the verbose block, which printed nxmodes in place of nymodes and nzmodes. #188 init_ConvectedVortex: cases meanFlowDir = +/-2 assigned vel_x = v_vort and vel_y = u_vort, giving a non-solenoidal, essentially irrotational field instead of a vortex. Un-swap them so only the mean flow term moves to the y component. Make "no mean flow" an explicit case 0 and abort on the host for out-of-range directions, which a bare default silently accepted. #189 MacProj: the known-edgestate mac_sync_compute overload passed a default-constructed MultiFab as the state, which ComputeAofs indexes unconditionally; pass an alias of Sync instead. mac_sync_compute built the divu constraint at prev_time only, omitting the 0.5*dt*dsdt half-time correction that velocity_advection and scalar_advection apply, so the sync reconstructed edge states with a different divu than the advance used. test_umac_periodic did host-side FArrayBox work on device-resident data; give it pinned host copies. #190 ScalMinMax: both limiters seeded the running maximum with numeric_limits::min(), the smallest positive normal rather than lowest(), so the upper clip bound was wrong for any all-negative stencil. With lowest() binding, skip covered cells so they keep the COVERED_VAL invariant. do_denminmax called ConservativeScalMinMax with Density as both the scalar and the normalizing density, clamping rho/rho == 1 against rhoo/rhoo == 1 -- a provable no-op; call ConvectiveScalMinMax on the raw density instead. The do_denminmax fix makes ns.do_denminmax=1 change results for the first time, so the eb_run2d and eb_run3d hotspot benchmarks need regolding. Co-Authored-By: Claude Opus 5 (1M context) --- Source/MacProj.cpp | 35 +++- Source/NS_derive.cpp | 167 +++++++++++++----- Source/NS_getForce.cpp | 19 +- Source/NavierStokesBase.cpp | 28 +-- Source/main.cpp | 3 +- Source/prob/prob_init.cpp | 19 +- Tutorials/HIT/NS_getForce.cpp | 20 +-- Tutorials/HIT/TurbulentForcing_def.H | 13 +- .../TaylorGreen/benchmarks/ViscBench.cpp | 34 ++-- 9 files changed, 232 insertions(+), 106 deletions(-) diff --git a/Source/MacProj.cpp b/Source/MacProj.cpp index 81b047e45..c69bb01a1 100644 --- a/Source/MacProj.cpp +++ b/Source/MacProj.cpp @@ -561,6 +561,14 @@ MacProj::mac_sync_compute (int level, forcing_term = std::make_unique(grids, dmap, num_state_comps, NavierStokesBase::nghost_force()); divu_fp.reset(ns_level.getDivCond(NavierStokesBase::nghost_force(),prev_time)); + // Get divu to time n+1/2, exactly as velocity_advection and + // scalar_advection do before predicting the edge states, so that the + // sync rebuilds the edge states the way the advance made them. + { + std::unique_ptr dsdt(ns_level.getDsdt(NavierStokesBase::nghost_force(),prev_time)); + MultiFab::Saxpy(*divu_fp, 0.5*dt, *dsdt, 0, 0, 1, NavierStokesBase::nghost_force()); + } + MultiFab& Gp = ns_level.get_old_data(Gradp_Type); visc_terms.setVal(0.0); // Initialize to make calls below safe @@ -770,12 +778,20 @@ MacProj::mac_sync_compute (int level, bool do_crse_add = false; bool do_fine_add = update_fluxreg; + // + // ComputeAofs always builds an Array4 from the state MultiFab, even though + // it does not use the values when the edge states are known. So it must be + // handed a defined MultiFab; an alias of Sync has the right BoxArray and + // DistributionMap. + // + MultiFab Smf(Sync, amrex::make_alias, Sync_indx, ncomp); + // // Compute the mac sync correction. // ns_level.ComputeAofs(Sync, /*Ssync_comp*/ Sync_indx, state_comp, ncomp, - /*State*/ MultiFab(), /*S_comp*/ int(),//not used when known_edgestates + /*State*/ Smf, /*S_comp*/ 0, //not used when known_edgestates /*forcing*/ nullptr, /*f_comp*/ int(), //not used when known_edgestates /*constraint divU*/ nullptr, //not used when known_edgestates fluxes, /*flux_comp*/ 0, @@ -988,6 +1004,11 @@ MacProj::test_umac_periodic (int level, Vector pshifts(27); std::vector< std::pair > isects; + // + // MultiFabCopyDescriptor and the FArrayBox comparison below both operate + // on the host, so give them host-accessible copies of u_mac. + // + Array h_umac; for (int dim = 0; dim < AMREX_SPACEDIM; dim++) { @@ -995,7 +1016,13 @@ MacProj::test_umac_periodic (int level, { Box eDomain = amrex::surroundingNodes(geom.Domain(),dim); - mfid[dim] = mfcd.RegisterMultiFab(&u_mac[dim]); + h_umac[dim].define(u_mac[dim].boxArray(), u_mac[dim].DistributionMap(), + u_mac[dim].nComp(), u_mac[dim].nGrowVect(), + MFInfo().SetArena(The_Pinned_Arena())); + amrex::dtoh_memcpy(h_umac[dim], u_mac[dim]); + Gpu::streamSynchronize(); + + mfid[dim] = mfcd.RegisterMultiFab(&h_umac[dim]); // How to combine pirm into one global pirm? // don't think std::vector::push_back() is thread safe @@ -1061,11 +1088,11 @@ MacProj::test_umac_periodic (int level, AMREX_ASSERT(pirm_i.m_srcBox.sameSize(pirm_i.m_dstBox)); AMREX_ASSERT(u_mac[dim].DistributionMap()[pirm_i.m_idx] == ParallelDescriptor::MyProc()); - diff.resize(pirm_i.m_srcBox, 1); + diff.resize(pirm_i.m_srcBox, 1, The_Pinned_Arena()); mfcd.FillFab(mfid[dim], pirm_i.m_fbid, diff); - diff.minus(u_mac[dim][pirm_i.m_idx],pirm_i.m_dstBox,diff.box(),0,0,1); + diff.minus(h_umac[dim][pirm_i.m_idx],pirm_i.m_dstBox,diff.box(),0,0,1); const Real max_norm = diff.norm(0); diff --git a/Source/NS_derive.cpp b/Source/NS_derive.cpp index ca60ed8d3..e7bf02f03 100644 --- a/Source/NS_derive.cpp +++ b/Source/NS_derive.cpp @@ -131,25 +131,51 @@ namespace derive_functions // Need to check if there are covered cells in neighbours -- // -- if so, use one-sided difference computation (but still quadratic) if (!flag_fab(i,j,k).isConnected( 1,0,0)) { - vx = - (c0 * dat_arr(i ,j,k,1) - + c1 * dat_arr(i-1,j,k,1) - + c2 * dat_arr(i-2,j,k,1)) * idx; + if (flag_fab(i,j,k).isConnected(-1,0,0)) { + if (flag_fab(i-1,j,k).isConnected(-1,0,0)) { + vx = - (c0 * dat_arr(i ,j,k,1) + + c1 * dat_arr(i-1,j,k,1) + + c2 * dat_arr(i-2,j,k,1)) * idx; + } else { + // Only one uncovered cell that way, drop to linear + vx = (dat_arr(i ,j,k,1) - dat_arr(i-1,j,k,1)) * idx; + } + } + // Covered on both sides, leave the derivative at zero } else if (!flag_fab(i,j,k).isConnected(-1,0,0)) { - vx = (c0 * dat_arr(i ,j,k,1) - + c1 * dat_arr(i+1,j,k,1) - + c2 * dat_arr(i+2,j,k,1)) * idx; + if (flag_fab(i+1,j,k).isConnected( 1,0,0)) { + vx = (c0 * dat_arr(i ,j,k,1) + + c1 * dat_arr(i+1,j,k,1) + + c2 * dat_arr(i+2,j,k,1)) * idx; + } else { + // Only one uncovered cell that way, drop to linear + vx = (dat_arr(i+1,j,k,1) - dat_arr(i ,j,k,1)) * idx; + } } else { vx = 0.5 * (dat_arr(i+1,j,k,1) - dat_arr(i-1,j,k,1)) * idx; } // Do the same in y-direction if (!flag_fab(i,j,k).isConnected( 0,1,0)) { - uy = - (c0 * dat_arr(i,j ,k,0) - + c1 * dat_arr(i,j-1,k,0) - + c2 * dat_arr(i,j-2,k,0)) * idy; + if (flag_fab(i,j,k).isConnected(0,-1,0)) { + if (flag_fab(i,j-1,k).isConnected(0,-1,0)) { + uy = - (c0 * dat_arr(i,j ,k,0) + + c1 * dat_arr(i,j-1,k,0) + + c2 * dat_arr(i,j-2,k,0)) * idy; + } else { + // Only one uncovered cell that way, drop to linear + uy = (dat_arr(i,j ,k,0) - dat_arr(i,j-1,k,0)) * idy; + } + } + // Covered on both sides, leave the derivative at zero } else if (!flag_fab(i,j,k).isConnected(0,-1,0)) { - uy = (c0 * dat_arr(i,j ,k,0) - + c1 * dat_arr(i,j+1,k,0) - + c2 * dat_arr(i,j+2,k,0)) * idy; + if (flag_fab(i,j+1,k).isConnected( 0,1,0)) { + uy = (c0 * dat_arr(i,j ,k,0) + + c1 * dat_arr(i,j+1,k,0) + + c2 * dat_arr(i,j+2,k,0)) * idy; + } else { + // Only one uncovered cell that way, drop to linear + uy = (dat_arr(i,j+1,k,0) - dat_arr(i,j ,k,0)) * idy; + } } else { uy = 0.5 * (dat_arr(i,j+1,k,0) - dat_arr(i,j-1,k,0)) * idy; } @@ -167,20 +193,35 @@ namespace derive_functions // -- if so, use one-sided difference computation (but still quadratic) if (!flag_fab(i,j,k).isConnected( 1,0,0)) { // Covered cell to the right, go fish left - vx = - (c0 * dat_arr(i ,j,k,1) - + c1 * dat_arr(i-1,j,k,1) - + c2 * dat_arr(i-2,j,k,1)) * idx; - wx = - (c0 * dat_arr(i ,j,k,2) - + c1 * dat_arr(i-1,j,k,2) - + c2 * dat_arr(i-2,j,k,2)) * idx; + if (flag_fab(i,j,k).isConnected(-1,0,0)) { + if (flag_fab(i-1,j,k).isConnected(-1,0,0)) { + vx = - (c0 * dat_arr(i ,j,k,1) + + c1 * dat_arr(i-1,j,k,1) + + c2 * dat_arr(i-2,j,k,1)) * idx; + wx = - (c0 * dat_arr(i ,j,k,2) + + c1 * dat_arr(i-1,j,k,2) + + c2 * dat_arr(i-2,j,k,2)) * idx; + } else { + // Only one uncovered cell that way, drop to linear + vx = (dat_arr(i ,j,k,1) - dat_arr(i-1,j,k,1)) * idx; + wx = (dat_arr(i ,j,k,2) - dat_arr(i-1,j,k,2)) * idx; + } + } + // Covered on both sides, leave the derivative at zero } else if (!flag_fab(i,j,k).isConnected(-1,0,0)) { // Covered cell to the left, go fish right - vx = (c0 * dat_arr(i ,j,k,1) - + c1 * dat_arr(i+1,j,k,1) - + c2 * dat_arr(i+2,j,k,1)) * idx; - wx = (c0 * dat_arr(i ,j,k,2) - + c1 * dat_arr(i+1,j,k,2) - + c2 * dat_arr(i+2,j,k,2)) * idx; + if (flag_fab(i+1,j,k).isConnected( 1,0,0)) { + vx = (c0 * dat_arr(i ,j,k,1) + + c1 * dat_arr(i+1,j,k,1) + + c2 * dat_arr(i+2,j,k,1)) * idx; + wx = (c0 * dat_arr(i ,j,k,2) + + c1 * dat_arr(i+1,j,k,2) + + c2 * dat_arr(i+2,j,k,2)) * idx; + } else { + // Only one uncovered cell that way, drop to linear + vx = (dat_arr(i+1,j,k,1) - dat_arr(i ,j,k,1)) * idx; + wx = (dat_arr(i+1,j,k,2) - dat_arr(i ,j,k,2)) * idx; + } } else { // No covered cells right or left, use standard stencil vx = 0.5 * (dat_arr(i+1,j,k,1) - dat_arr(i-1,j,k,1)) * idx; @@ -188,38 +229,68 @@ namespace derive_functions } // Do the same in y-direction if (!flag_fab(i,j,k).isConnected(0, 1,0)) { - uy = - (c0 * dat_arr(i,j ,k,0) - + c1 * dat_arr(i,j-1,k,0) - + c2 * dat_arr(i,j-2,k,0)) * idy; - wy = - (c0 * dat_arr(i,j ,k,2) - + c1 * dat_arr(i,j-1,k,2) - + c2 * dat_arr(i,j-2,k,2)) * idy; + if (flag_fab(i,j,k).isConnected(0,-1,0)) { + if (flag_fab(i,j-1,k).isConnected(0,-1,0)) { + uy = - (c0 * dat_arr(i,j ,k,0) + + c1 * dat_arr(i,j-1,k,0) + + c2 * dat_arr(i,j-2,k,0)) * idy; + wy = - (c0 * dat_arr(i,j ,k,2) + + c1 * dat_arr(i,j-1,k,2) + + c2 * dat_arr(i,j-2,k,2)) * idy; + } else { + // Only one uncovered cell that way, drop to linear + uy = (dat_arr(i,j ,k,0) - dat_arr(i,j-1,k,0)) * idy; + wy = (dat_arr(i,j ,k,2) - dat_arr(i,j-1,k,2)) * idy; + } + } + // Covered on both sides, leave the derivative at zero } else if (!flag_fab(i,j,k).isConnected(0,-1,0)) { - uy = (c0 * dat_arr(i,j ,k,0) - + c1 * dat_arr(i,j+1,k,0) - + c2 * dat_arr(i,j+2,k,0)) * idy; - wy = (c0 * dat_arr(i,j ,k,2) - + c1 * dat_arr(i,j+1,k,2) - + c2 * dat_arr(i,j+2,k,2)) * idy; + if (flag_fab(i,j+1,k).isConnected(0, 1,0)) { + uy = (c0 * dat_arr(i,j ,k,0) + + c1 * dat_arr(i,j+1,k,0) + + c2 * dat_arr(i,j+2,k,0)) * idy; + wy = (c0 * dat_arr(i,j ,k,2) + + c1 * dat_arr(i,j+1,k,2) + + c2 * dat_arr(i,j+2,k,2)) * idy; + } else { + // Only one uncovered cell that way, drop to linear + uy = (dat_arr(i,j+1,k,0) - dat_arr(i,j ,k,0)) * idy; + wy = (dat_arr(i,j+1,k,2) - dat_arr(i,j ,k,2)) * idy; + } } else { uy = 0.5 * (dat_arr(i,j+1,k,0) - dat_arr(i,j-1,k,0)) * idy; wy = 0.5 * (dat_arr(i,j+1,k,2) - dat_arr(i,j-1,k,2)) * idy; } // Do the same in z-direction if (!flag_fab(i,j,k).isConnected(0,0, 1)) { - uz = - (c0 * dat_arr(i,j,k ,0) - + c1 * dat_arr(i,j,k-1,0) - + c2 * dat_arr(i,j,k-2,0)) * idz; - vz = - (c0 * dat_arr(i,j,k ,1) - + c1 * dat_arr(i,j,k-1,1) - + c2 * dat_arr(i,j,k-2,1)) * idz; + if (flag_fab(i,j,k).isConnected(0,0,-1)) { + if (flag_fab(i,j,k-1).isConnected(0,0,-1)) { + uz = - (c0 * dat_arr(i,j,k ,0) + + c1 * dat_arr(i,j,k-1,0) + + c2 * dat_arr(i,j,k-2,0)) * idz; + vz = - (c0 * dat_arr(i,j,k ,1) + + c1 * dat_arr(i,j,k-1,1) + + c2 * dat_arr(i,j,k-2,1)) * idz; + } else { + // Only one uncovered cell that way, drop to linear + uz = (dat_arr(i,j,k ,0) - dat_arr(i,j,k-1,0)) * idz; + vz = (dat_arr(i,j,k ,1) - dat_arr(i,j,k-1,1)) * idz; + } + } + // Covered on both sides, leave the derivative at zero } else if (!flag_fab(i,j,k).isConnected(0,0,-1)) { - uz = (c0 * dat_arr(i,j,k ,0) - + c1 * dat_arr(i,j,k+1,0) - + c2 * dat_arr(i,j,k+2,0)) * idz; - vz = (c0 * dat_arr(i,j,k ,1) - + c1 * dat_arr(i,j,k+1,1) - + c2 * dat_arr(i,j,k+2,1)) * idz; + if (flag_fab(i,j,k+1).isConnected(0,0, 1)) { + uz = (c0 * dat_arr(i,j,k ,0) + + c1 * dat_arr(i,j,k+1,0) + + c2 * dat_arr(i,j,k+2,0)) * idz; + vz = (c0 * dat_arr(i,j,k ,1) + + c1 * dat_arr(i,j,k+1,1) + + c2 * dat_arr(i,j,k+2,1)) * idz; + } else { + // Only one uncovered cell that way, drop to linear + uz = (dat_arr(i,j,k+1,0) - dat_arr(i,j,k ,0)) * idz; + vz = (dat_arr(i,j,k+1,1) - dat_arr(i,j,k ,1)) * idz; + } } else { uz = 0.5 * (dat_arr(i,j,k+1,0) - dat_arr(i,j,k-1,0)) * idz; vz = 0.5 * (dat_arr(i,j,k+1,1) - dat_arr(i,j,k-1,1)) * idz; diff --git a/Source/NS_getForce.cpp b/Source/NS_getForce.cpp index e66bee13b..6d8cb0278 100644 --- a/Source/NS_getForce.cpp +++ b/Source/NS_getForce.cpp @@ -92,10 +92,13 @@ NavierStokesBase::getForce (FArrayBox& force, #endif // Compute min/max - for (int n=0; n(scomp+n) << " / " - << State.max(scomp+n) << '\n'; + // State may not hold all of the state components (scalar callers pass + // a velocity-only FAB), so index it by its own layout rather than by + // the absolute state component scomp+n. + for (int n=0; n(n) << " / " + << State.max(n) << '\n'; } for (int n=auxScomp; n(scomp+n) << " / " - << force.max(scomp+n) << '\n'; + << force.min(fcomp) << " / " + << force.max(fcomp) << '\n'; } amrex::Print() << "NavierStokesBase::getForce(): Leaving..." diff --git a/Source/NavierStokesBase.cpp b/Source/NavierStokesBase.cpp index 5096470cf..6d71f67b5 100644 --- a/Source/NavierStokesBase.cpp +++ b/Source/NavierStokesBase.cpp @@ -2774,17 +2774,15 @@ NavierStokesBase::scalar_advection_update (Real dt, // Must do FillPatch here instead of MF iterator because we need the // boundary values in the old data (especially at inflow) // - const int index_new_s = Density; - const int index_new_rho = Density; - const int index_old_s = index_new_s - Density; - const int index_old_rho = index_new_rho - Density; - FillPatchIterator S_fpi(*this,S_old,1,prev_time,State_Type,Density,1); MultiFab& Smf=S_fpi.get_mf(); - ConservativeScalMinMax(S_new, index_new_s, index_new_rho, - Smf, index_old_s, index_old_rho); - + // + // Clamp the new density to the min/max of the old-time density in + // the neighborhood. (The conservative s/rho form with s = rho is + // an identity, i.e. a no-op.) + // + ConvectiveScalMinMax(S_new, Density, Smf, 0); } ++sComp; @@ -4277,8 +4275,13 @@ NavierStokesBase::ConservativeScalMinMax ( amrex::MultiFab& Snew, const in amrex::ParallelFor(bx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept { +#ifdef AMREX_USE_EB + // Leave covered cells alone; their stencil holds no valid data. + if ( vfrac(i,j,k) == 0. ) { return; } +#endif + Real smn = std::numeric_limits::max(); - Real smx = std::numeric_limits::min(); + Real smx = std::numeric_limits::lowest(); #if (AMREX_SPACEDIM==3) int ks = -1; @@ -4335,8 +4338,13 @@ NavierStokesBase::ConvectiveScalMinMax ( amrex::MultiFab& Snew, const int amrex::ParallelFor(bx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept { +#ifdef AMREX_USE_EB + // Leave covered cells alone; their stencil holds no valid data. + if ( vfrac(i,j,k) == 0. ) { return; } +#endif + Real smn = std::numeric_limits::max(); - Real smx = std::numeric_limits::min(); + Real smx = std::numeric_limits::lowest(); #if (AMREX_SPACEDIM==3) int ks = -1; diff --git a/Source/main.cpp b/Source/main.cpp index 3967b3ec6..bcd5ae6bd 100644 --- a/Source/main.cpp +++ b/Source/main.cpp @@ -105,7 +105,8 @@ main (int argc, // if (Amr::RegridOnRestart()) { - if ( (amrptr->levelSteps(0) >= max_step ) || + if ( ( (max_step >= 0) && + (amrptr->levelSteps(0) >= max_step) ) || ( (stop_time >= 0.0) && (amrptr->cumTime() >= stop_time) ) ) { diff --git a/Source/prob/prob_init.cpp b/Source/prob/prob_init.cpp index d7f023df7..fe3ffb211 100644 --- a/Source/prob/prob_init.cpp +++ b/Source/prob/prob_init.cpp @@ -630,6 +630,14 @@ void NavierStokes::init_ConvectedVortex (Box const& vbx, { const auto domlo = amrex::lbound(domain); + // + // amrex::Abort() cannot be called from within the device lambda below, so + // validate the mean flow direction here. + // + if ( IC.meanFlowDir < -3 || IC.meanFlowDir > 3 ) { + amrex::Abort("\n init_ConvectedVortex: prob.meanFlowDir must be 0 (no mean flow) or +/-1, +/-2, +/-3\n in the inputs file."); + } + amrex::ParallelFor(vbx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept { AMREX_D_TERM(Real x = problo[0] + (i - domlo.x + 0.5)*dx[0];, @@ -661,13 +669,13 @@ void NavierStokes::init_ConvectedVortex (Box const& vbx, vel(i,j,k,2) = w_vort); break; case 2 : - AMREX_D_TERM(vel(i,j,k,0) = v_vort;, - vel(i,j,k,1) = IC.meanFlowMag + u_vort;, + AMREX_D_TERM(vel(i,j,k,0) = u_vort;, + vel(i,j,k,1) = IC.meanFlowMag + v_vort;, vel(i,j,k,2) = w_vort); break; case -2 : - AMREX_D_TERM(vel(i,j,k,0) = v_vort;, - vel(i,j,k,1) = -IC.meanFlowMag + u_vort;, + AMREX_D_TERM(vel(i,j,k,0) = u_vort;, + vel(i,j,k,1) = -IC.meanFlowMag + v_vort;, vel(i,j,k,2) = w_vort); break; case 3 : @@ -680,7 +688,8 @@ void NavierStokes::init_ConvectedVortex (Box const& vbx, vel(i,j,k,1) = -IC.meanFlowMag + v_vort;, vel(i,j,k,2) = w_vort); break; - default : // no mean flow, i.e. the vortex alone + case 0 : // no mean flow, i.e. the vortex alone + default : AMREX_D_TERM(vel(i,j,k,0) = u_vort;, vel(i,j,k,1) = v_vort;, vel(i,j,k,2) = w_vort); diff --git a/Tutorials/HIT/NS_getForce.cpp b/Tutorials/HIT/NS_getForce.cpp index b08f466a9..b25379a76 100644 --- a/Tutorials/HIT/NS_getForce.cpp +++ b/Tutorials/HIT/NS_getForce.cpp @@ -318,9 +318,10 @@ NavierStokesBase::getForce (FArrayBox& force, int ff_jlo = jlo/TurbulentForcing::ff_factor; int ff_klo = klo/TurbulentForcing::ff_factor; - int ff_ihi = (ihi+1)/TurbulentForcing::ff_factor; - int ff_jhi = (jhi+1)/TurbulentForcing::ff_factor; - int ff_khi = (khi+1)/TurbulentForcing::ff_factor; + // +1 so that the (ff_i+1) node used by the interpolation below always exists + int ff_ihi = ihi/TurbulentForcing::ff_factor + 1; + int ff_jhi = jhi/TurbulentForcing::ff_factor + 1; + int ff_khi = khi/TurbulentForcing::ff_factor + 1; // adjust for ghost cells if (ilo < (ff_ilo*TurbulentForcing::ff_factor)) { @@ -332,20 +333,13 @@ NavierStokesBase::getForce (FArrayBox& force, if (klo < (ff_klo*TurbulentForcing::ff_factor)) { ff_klo=ff_klo-1; } - if (ihi == (ff_ihi*TurbulentForcing::ff_factor)) { - ff_ihi=ff_ihi+1; - } - if (jhi == (ff_jhi*TurbulentForcing::ff_factor)) { - ff_jhi=ff_jhi+1; - } - if (khi == (ff_khi*TurbulentForcing::ff_factor)) { - ff_khi=ff_khi+1; - } // allocate coarse force array Box ffbx(IntVect(ff_ilo, ff_jlo, ff_klo), IntVect(ff_ihi, ff_jhi, ff_khi)); - // not sure if want elixir, gpu::sync, or async_arena here... FArrayBox ff_force(ffbx,AMREX_SPACEDIM); + // The Elixir defers freeing ff_force's device memory until the + // asynchronous kernels below that read it have completed. No-op on CPU. + Elixir ff_force_i = ff_force.elixir(); const auto& ffarr = ff_force.array(); // Construct node-based coarse forcing diff --git a/Tutorials/HIT/TurbulentForcing_def.H b/Tutorials/HIT/TurbulentForcing_def.H index d9a4325bd..429015f0b 100644 --- a/Tutorials/HIT/TurbulentForcing_def.H +++ b/Tutorials/HIT/TurbulentForcing_def.H @@ -105,13 +105,18 @@ TurbulentForcing::init_turbulent_forcing (const amrex::GpuArray