From c806e0cae88194e340042b2bfd9a24a802b20152 Mon Sep 17 00:00:00 2001 From: Ann Almgren Date: Sun, 27 Sep 2026 16:55:15 -0700 Subject: [PATCH] Fix issues #186, #191-#193, #195-#200 Source/Diffusion.{H,cpp}, Source/NavierStokes.cpp (#200): - diffuse_tensor_Vsync solved with shear viscosity hard-coded to 1.0, ignoring the beta it was passed. Restore its use; thread a betaCC parameter through diffuse_Vsync so MLEBTensorOp::setEBShearViscosity gets the old-time cell-centered viscosity in EB builds. - diffuse_tensor_Vsync scaled only component 0 of the RHS by rhsscale while setScalars scaled the operator for all components. - diffuse_tensor_Vsync multiplied the rho_flag==3 RHS by old-time density while acoef and the caller both use new-time density. - diffuse_tensor_velocity never multiplied its RHS by rhsscale. - diffuse_scalar never called setCoarseFineBC on a level>0 solve with no coarse data, so MLMG assumed a refinement ratio of 2. Source/Projection.cpp (#199): - The RZ branch of scaleVar zeroed velocity ghosts outside the domain in the axial direction, destroying inflow values the nodal divu stencil reads. Scale them by radius instead, as the pre-port Fortran radmpyvel did, and mirror the change in rescaleVar. - computeRhoG's 3D y-hi face mixed density rows in rho_ii, both in the main loop and in the x-lo ext_dir edge branch; the latter also read one cell outside the rho FAB. Source/NavierStokesBase.cpp (#198): - Restart with ns.gradp_in_checkpoint=0 called computeGradP on a Gradp_Type StateData that AmrLevel::restart had left undefined. Define it first, using the midpoint of Press's own interval so that Gradp's curTime()/prevTime() match Press's exactly. Source/NavierStokesBase.cpp (#197): - post_timestep_particle passed an undefined MultiFab to Timestamp whenever particles.timestamp_indices was not set. - post_timestep_particle and ParticleDerive built MultiFabs on level lev's BoxArray with the *current* level's FabFactory. - ParticleDerive("total_particle_count") accumulated fine into coarse with a host BoxIterator loop over device data; it is now a ParallelFor with Gpu::Atomic::AddNoRet, which also removes an OpenMP race between tiles. Source/NavierStokesBase.cpp (#196): - The cut-cell CrseAdd/FineAdd overloads multiply by the EB area fraction, but AMReX-Hydro already area-weighted the fluxes, so cut c/f faces were refluxed with ap^2. Hand those overloads un-weighted copies. - The StateRedist "state" for the mac_sync was copied from Vsync/Ssync without filling ghost cells, making the sync correction near box boundaries grid-decomposition dependent. - use_wts_in_divnc was sitting in ApplyMLRedistribution's fac_for_deltaR slot, so dm_as_fine was scaled by +dt even in the mac_sync, where the fluxes go to FineAdd with -dt. Source/NavierStokesBase.{H,cpp}, Source/NS_average.cpp (#195, #186): - time_avg/time_avg_fluct/dt_avg were sized only in post_init and post_restart, so a level created by a later regrid indexed them out of bounds. Add grow_avg_vectors() and call it from post_regrid and from both init() overloads. - The brand-new-level init() never initialized Average_Type. FillCoarsePatch it, as the init(AmrLevel&) twin does, and seed the level's accumulators from the coarser level so that the interpolated integral is normalized consistently. - checkPoint wrote the single shared TimeAverage file once per level with trunc, so only the finest level's data survived, and dt_avg was never checkpointed. Level 0 now writes one (time_avg, time_avg_fluct, dt_avg) triple per level, and post_restart reads forward to its own level's triple. - time_average dereferenced Average_Type old data that is never allocated on a fresh start with ns.init_iter=0. Source/SyncRegister.cpp (#193): - The bndry_mask threshold was SPACEDIM^SPACEDIM-0.5 (26.5 in 3D) but the accumulated count maxes out at 2^SPACEDIM, so in 3D no node was ever masked out of the sync-projection RHS. 2D is unchanged; 3D multilevel answers move. - outflow_dirs was one element short of the number of directions the loop below can record. Source/NS_LES.cpp (#192): - The Smagorinsky branch added each velocity-gradient component to itself rather than to its transpose, so mu_t was built from |grad u| instead of |S| and was overpredicted in every rotational flow. Pair each component with its transpose. Util/ConvertCheckpoint/ConvertCheckpointGrids.cpp (#191): - nsets_save was a 1-element Vector written at indices 0..ndesc-1 and resized only at lev==1. Drop the global and derive nsets from the data pointers, the way StateData::checkPoint does. - ConvertData unconditionally dereferenced old_data, which is null when a state was checkpointed with nsets==1. - The trailing state types were given zero ghost cells, but cell_cons_interp's slope stencil reads one cell outside the coarse fab; and the source ghosts were never filled, so periodic-boundary ghosts fed the setVal(10.) sentinel into the slopes. Util/ConvertCheckpoint/Make.package: - Drop the reference to AMReX_FABUTIL_$(DIM)D.F, which no longer exists in AMReX; without this the utility does not build at all. Note that regression baselines move for 3D multilevel runs (#193), any viscous multilevel run (#200), EB cut-cell reflux (#196), RZ with axial inflow (#199) and Smagorinsky LES (#192). Co-Authored-By: Claude Opus 5 (1M context) --- Source/Diffusion.H | 2 + Source/Diffusion.cpp | 51 ++-- Source/NS_LES.cpp | 8 +- Source/NS_average.cpp | 10 + Source/NavierStokes.cpp | 4 +- Source/NavierStokesBase.H | 6 + Source/NavierStokesBase.cpp | 241 +++++++++++++++--- Source/Projection.cpp | 22 +- Source/SyncRegister.cpp | 9 +- .../ConvertCheckpointGrids.cpp | 76 +++--- Util/ConvertCheckpoint/Make.package | 1 - 11 files changed, 333 insertions(+), 97 deletions(-) diff --git a/Source/Diffusion.H b/Source/Diffusion.H index 78c2a8f5b..b44157b11 100644 --- a/Source/Diffusion.H +++ b/Source/Diffusion.H @@ -118,6 +118,7 @@ public: const amrex::MultiFab& rho_half, int rho_flag, const amrex::MultiFab* const* beta, + const amrex::MultiFab* betaCC, int betaComp = 0, bool update_fluxreg = true); @@ -127,6 +128,7 @@ public: const amrex::MultiFab& rho_half, int rho_flag, const amrex::MultiFab* const* beta, + const amrex::MultiFab* betaCC, int betaComp, bool update_fluxreg); diff --git a/Source/Diffusion.cpp b/Source/Diffusion.cpp index b9c8fc92b..13250590f 100644 --- a/Source/Diffusion.cpp +++ b/Source/Diffusion.cpp @@ -394,6 +394,14 @@ Diffusion::diffuse_scalar (const Vector& S_old, } opn.setCoarseFineBC(Solnc.get(), cratio[0]); } + else if (level > 0) { + // + // No coarse data supplied, but this is a level > 0 solve, so + // we still must tell MLMG the true refinement ratio for the + // homogeneous Dirichlet coarse-fine boundary condition. + // + opn.setCoarseFineBC(nullptr, cratio[0]); + } MultiFab::Copy(Soln,*S_old[0],sigma,0,nComp,ng); if (rho_flag == 2) { #ifdef AMREX_USE_OMP @@ -517,6 +525,9 @@ Diffusion::diffuse_scalar (const Vector& S_old, } opnp1.setCoarseFineBC(Solnc.get(), cratio[0]); } + else if (level > 0) { + opnp1.setCoarseFineBC(nullptr, cratio[0]); + } #ifdef AMREX_USE_OMP #pragma omp parallel if (Gpu::notInLaunchRegion()) #endif @@ -847,7 +858,6 @@ Diffusion::diffuse_tensor_velocity (Real dt, // const Real tol_abs = 0.0; // cribbing from scalar const Real tol_rel = visc_tol; - const Real tol_abs = get_scaled_abs_tol(Rhs, visc_tol); LPInfo info; info.setAgglomeration(agglomeration); @@ -889,10 +899,10 @@ Diffusion::diffuse_tensor_velocity (Real dt, tensorop.setLevelBC(0, &Soln); } + Real rhsscale = 1.0; { MultiFab acoef; std::pair scalars; - Real rhsscale = 1.0; const MultiFab& rho = (rho_flag == 1) ? rho_half : navier_stokes->get_new_data(State_Type); const int rho_comp = (rho_flag == 1) ? 0 : Density; computeAlpha(acoef, scalars, a, b, @@ -901,6 +911,12 @@ Diffusion::diffuse_tensor_velocity (Real dt, tensorop.setScalars(scalars.first, scalars.second); tensorop.setACoeffs(0, acoef); } + // + // computeAlpha scaled the operator scalars by rhsscale; the RHS must be + // scaled to match (cf. diffuse_scalar and diffuse_Ssync). + // + Rhs.mult(rhsscale,0,AMREX_SPACEDIM); + const Real tol_abs = get_scaled_abs_tol(Rhs, visc_tol); #ifdef AMREX_USE_EB setViscosity(tensorop, betanp1, betaComp, *betanp1CC); @@ -966,6 +982,7 @@ Diffusion::diffuse_Vsync (MultiFab& Vsync, const MultiFab& rho_half, int rho_flag, const MultiFab* const* beta, + const MultiFab* betaCC, int betaComp, bool update_fluxreg) { @@ -979,7 +996,7 @@ Diffusion::diffuse_Vsync (MultiFab& Vsync, AMREX_ASSERT(beta[d]->min(0,0) >= 0.0); #endif - diffuse_tensor_Vsync(Vsync,dt,be_cn_theta,rho_half,rho_flag,beta,betaComp,update_fluxreg); + diffuse_tensor_Vsync(Vsync,dt,be_cn_theta,rho_half,rho_flag,beta,betaCC,betaComp,update_fluxreg); // // applyBC has put "incorrect" values in the ghost cells // outside external Dirichlet boundaries. Reset these to zero @@ -1017,8 +1034,9 @@ Diffusion::diffuse_tensor_Vsync (MultiFab& Vsync, Real be_cn_theta, const MultiFab& rho_half, int rho_flag, - const MultiFab* const* /*beta*/, - int /*betaComp*/, + const MultiFab* const* beta, + const MultiFab* betaCC, + int betaComp, bool update_fluxreg) { AMREX_ASSERT(rho_flag == 1 || rho_flag == 3); @@ -1045,7 +1063,9 @@ Diffusion::diffuse_tensor_Vsync (MultiFab& Vsync, { const Box& bx = mfi.tilebox(); auto const& rhs = Rhs.array(mfi); - auto const& rho = (rho_flag == 1) ? rho_half.array(mfi) : navier_stokes->get_old_data(State_Type).array(mfi,Density); + // NOTE: must use the new-time density here to match the acoef built + // below and mac_sync's normalization of Vsync by rho^{n+1}. + auto const& rho = (rho_flag == 1) ? rho_half.array(mfi) : navier_stokes->get_new_data(State_Type).array(mfi,Density); amrex::ParallelFor(bx, [rhs, rho] AMREX_GPU_DEVICE (int i, int j, int k) noexcept @@ -1113,22 +1133,13 @@ Diffusion::diffuse_tensor_Vsync (MultiFab& Vsync, tensorop.setACoeffs(0, acoef); } - { - FluxBoxes fb_bcoef; - MultiFab** face_bcoef = nullptr; - face_bcoef = fb_bcoef.define(navier_stokes); - for (int dir=0; dirsetVal(1.0); - } - #ifdef AMREX_USE_EB - MultiFab bcoefCC(grids,dmap,1,0,MFInfo(),navier_stokes->Factory()); - bcoefCC.setVal(1.0); - setViscosity(tensorop, face_bcoef, 0, bcoefCC); + AMREX_ALWAYS_ASSERT(betaCC != nullptr); + setViscosity(tensorop, beta, betaComp, *betaCC); #else - setViscosity(tensorop, face_bcoef, 0); + amrex::ignore_unused(betaCC); + setViscosity(tensorop, beta, betaComp); #endif - } MLMG mlmg(tensorop); if (max_iter > 0) { @@ -1144,7 +1155,7 @@ Diffusion::diffuse_tensor_Vsync (MultiFab& Vsync, mlmg.setMaxFmgIter(max_fmg_iter); mlmg.setVerbose(verbose); - Rhs.mult(rhsscale,0,1); + Rhs.mult(rhsscale,0,AMREX_SPACEDIM); mlmg.setFinalFillBC(true); mlmg.solve({&Soln}, {&Rhs}, tol_rel, tol_abs); diff --git a/Source/NS_LES.cpp b/Source/NS_LES.cpp index 3c48a0f42..3f647492c 100644 --- a/Source/NS_LES.cpp +++ b/Source/NS_LES.cpp @@ -125,7 +125,13 @@ NavierStokesBase::calc_mut_LES(MultiFab* mu_LES[AMREX_SPACEDIM], const Real time Real smag = 0; for (int i_symij = 0; i_symij < dim_fluxes; ++i_symij) { - Real symij = src(i,j,k,i_symij) + src(i,j,k,i_symij); + // compVelGrad stores du_m/dx_n in component AMREX_SPACEDIM*n+m, + // so the transpose of flat index c is + // (c%AMREX_SPACEDIM)*AMREX_SPACEDIM + c/AMREX_SPACEDIM. Pairing + // each component with its transpose gives symij = 2*S_mn, so + // smag below is 2*S:S and mu_t = (Cs*dx)^2*sqrt(2 S_ij S_ij). + int i_symji = (i_symij%AMREX_SPACEDIM)*AMREX_SPACEDIM + i_symij/AMREX_SPACEDIM; + Real symij = src(i,j,k,i_symij) + src(i,j,k,i_symji); smag += symij * symij; } diff --git a/Source/NS_average.cpp b/Source/NS_average.cpp index 80c3282a6..e6be18a5a 100644 --- a/Source/NS_average.cpp +++ b/Source/NS_average.cpp @@ -26,6 +26,16 @@ NavierStokesBase::time_average(amrex::Real& a_time_avg, amrex::Real& a_time_av { MultiFab& Sstate = get_new_data(State_Type); MultiFab& Savg = get_new_data(Average_Type); + // + // Old data may not have been allocated yet, e.g. when this is called from + // post_init and no advance has taken place (ns.init_iter <= 0). + // Zero is the right seed: the kernel below copies S_avg back into + // S_avg_old on every sample, so both time levels agree by construction. + // + if (! state[Average_Type].hasOldData()) { + state[Average_Type].allocOldData(); + state[Average_Type].oldData().setVal(0.); + } MultiFab& Savg_old = get_old_data(Average_Type); #ifdef _OPENMP diff --git a/Source/NavierStokes.cpp b/Source/NavierStokes.cpp index 4081d8d27..2edc390e3 100644 --- a/Source/NavierStokes.cpp +++ b/Source/NavierStokes.cpp @@ -1561,7 +1561,9 @@ NavierStokes::mac_sync () loc_viscn = fb_viscn.define(this); getViscosity(loc_viscn, viscTime); - diffusion->diffuse_Vsync(Vsync,dt,be_cn_theta,Rh,rho_flag,loc_viscn,0); + // viscn_cc holds the cell-centered viscosity at prevTime; needed by + // MLEBTensorOp::setEBShearViscosity in EB builds. + diffusion->diffuse_Vsync(Vsync,dt,be_cn_theta,Rh,rho_flag,loc_viscn,viscn_cc,0); } FluxBoxes fb_SC; diff --git a/Source/NavierStokesBase.H b/Source/NavierStokesBase.H index 0ac31d60e..bc4ce8a71 100644 --- a/Source/NavierStokesBase.H +++ b/Source/NavierStokesBase.H @@ -787,6 +787,12 @@ public: static amrex::Vector time_avg; static amrex::Vector time_avg_fluct; static amrex::Vector dt_avg; + // + // Grow time_avg/time_avg_fluct/dt_avg so that every level up to + // a_finest_level has an entry. New entries are zero-filled; existing + // accumulations are preserved. + // + static void grow_avg_vectors (int a_finest_level); static int avg_interval; static int compute_fluctuations; // diff --git a/Source/NavierStokesBase.cpp b/Source/NavierStokesBase.cpp index 6d71f67b5..f1920d7bf 100644 --- a/Source/NavierStokesBase.cpp +++ b/Source/NavierStokesBase.cpp @@ -852,6 +852,25 @@ NavierStokesBase::calc_dsdt (Real /*time*/, } } +// +// Grow the on-the-fly averaging accumulators so that every level up to +// a_finest_level has an entry. A regrid can create levels that did not exist +// when these were sized in post_init/post_restart. resize() zero-fills the new +// entries and leaves the existing accumulations alone. +// +void +NavierStokesBase::grow_avg_vectors (int a_finest_level) +{ + const int nlev = a_finest_level + 1; + + if (static_cast(NavierStokesBase::time_avg.size()) < nlev) + { + NavierStokesBase::time_avg.resize(nlev,0.); + NavierStokesBase::time_avg_fluct.resize(nlev,0.); + NavierStokesBase::dt_avg.resize(nlev,0.); + } +} + void NavierStokesBase::checkPoint (const std::string& dir, std::ostream& os, @@ -860,7 +879,11 @@ NavierStokesBase::checkPoint (const std::string& dir, { AmrLevel::checkPoint(dir, os, how, dump_old); - if (avg_interval > 0) + // + // There is a single TimeAverage file per checkpoint, so only level 0 + // writes it -- and it writes the data for every level. + // + if (avg_interval > 0 && level == 0) { VisMF::IO_Buffer io_buffer(VisMF::IO_Buffer_Size); @@ -882,8 +905,15 @@ NavierStokesBase::checkPoint (const std::string& dir, // write out title line TImeAverageFile << "Writing time_average to checkpoint\n"; - TImeAverageFile << NavierStokesBase::time_avg[level] << "\n"; - TImeAverageFile << NavierStokesBase::time_avg_fluct[level] << "\n"; + // + // One (time_avg, time_avg_fluct, dt_avg) triple per level. + // + for (int lev = 0; lev <= parent->finestLevel(); lev++) + { + TImeAverageFile << NavierStokesBase::time_avg[lev] << "\n"; + TImeAverageFile << NavierStokesBase::time_avg_fluct[lev] << "\n"; + TImeAverageFile << NavierStokesBase::dt_avg[lev] << "\n"; + } } } @@ -1732,6 +1762,7 @@ NavierStokesBase::init (AmrLevel &old) FillPatch(old,Gp_new,Gp_new.nGrow(),cur_pres_time,Gradp_Type,0,AMREX_SPACEDIM); if (avg_interval > 0){ + grow_avg_vectors(level); MultiFab& Save_new = get_new_data(Average_Type); FillPatch(old,Save_new,0,cur_time,Average_Type,0,AMREX_SPACEDIM*2); } @@ -1795,6 +1826,23 @@ NavierStokesBase::init () FillCoarsePatch(P_new,0,cur_pres_time,Press_Type,0,1); FillCoarsePatch(Gp_new,0,cur_pres_time,Gradp_Type,0,AMREX_SPACEDIM,Gp_new.nGrow()); // + // Get the best coarse time-average data. post_regrid has not run yet, so + // make room for this level here; the interpolated average is the integral + // over the *coarser* level's averaging window, so the accumulators must be + // seeded from the coarser level too or der_vel_avg will mis-normalize it. + // + if (avg_interval > 0) + { + grow_avg_vectors(level); + + MultiFab& Save_new = get_new_data(Average_Type); + FillCoarsePatch(Save_new,0,cur_time,Average_Type,0,AMREX_SPACEDIM*2); + + NavierStokesBase::time_avg[level] = NavierStokesBase::time_avg[level-1]; + NavierStokesBase::time_avg_fluct[level] = NavierStokesBase::time_avg_fluct[level-1]; + NavierStokesBase::dt_avg[level] = NavierStokesBase::dt_avg[level-1]; + } + // // Get best coarse divU and dSdt data. // if (have_divu) @@ -2445,6 +2493,14 @@ void NavierStokesBase::post_regrid (int lbase, int /*new_finest*/) { + // + // A regrid may have created levels that did not exist when the on-the-fly + // averaging data were sized in post_init/post_restart. + // + if (avg_interval > 0) { + grow_avg_vectors(parent->finestLevel()); + } + #ifdef AMREX_PARTICLES if (NSPC && level == lbase) { @@ -2514,9 +2570,19 @@ NavierStokesBase::post_restart () // read in title line std::getline(isp, line); - isp >> NavierStokesBase::time_avg[level]; - isp >> NavierStokesBase::time_avg_fluct[level]; - NavierStokesBase::dt_avg[level] = 0; + // + // The file holds one (time_avg, time_avg_fluct, dt_avg) triple per + // level, so read forward to this level's triple. Note that a + // checkpoint written by an older version of the code holds only + // level 0's (time_avg, time_avg_fluct) pair; the failed extractions + // then leave zeros behind, which is the best we can do. + // + for (int lev = 0; lev <= level; lev++) + { + isp >> NavierStokesBase::time_avg[level]; + isp >> NavierStokesBase::time_avg_fluct[level]; + isp >> NavierStokesBase::dt_avg[level]; + } } } @@ -2717,10 +2783,31 @@ NavierStokesBase::restart (Amr& papa, <<'\n'; // - // Compute GradP from the Pressure + // AmrLevel::restart skipped state[Gradp_Type].restart() because Gradp + // was not in the checkpoint, so the StateData is still undefined. + // Define it here (mirroring the Average_Type recovery in post_restart) + // before computing GradP from the Pressure. + // + Real cur_time = state[Press_Type].curTime(); + Real prev_time = state[Press_Type].prevTime(); + Real dt_gp = cur_time - prev_time; + + // + // Gradp_Type is registered StateDescriptor::Interval, so StateData::define + // builds new_time = [t,t+dt] and old_time = [t-dt,t]. Passing the midpoint + // of Press's own interval makes Gradp's curTime()/prevTime() match Press's + // exactly, so both get_data() lookups inside computeGradP resolve. + // + state[Gradp_Type].define(geom.Domain(), grids, dmap, desc_lst[Gradp_Type], + 0.5*(prev_time+cur_time), dt_gp, Factory()); + + computeGradP(cur_time); + + // + // Now allocate the old data and fill it. // - computeGradP(state[Press_Type].curTime()); - computeGradP(state[Press_Type].prevTime()); + state[Gradp_Type].allocOldData(); + computeGradP(prev_time); } define_workspace(); @@ -3916,7 +4003,9 @@ NavierStokesBase::post_timestep_particle (int crse_iteration) if (tindices.size() > 0) { - tmf.define(S_new.boxArray(), S_new.DistributionMap(), tindices.size(), ng, MFInfo(), Factory()); + // NOTE: must use the factory of the level being timestamped, + // not this level's. + tmf.define(S_new.boxArray(), S_new.DistributionMap(), tindices.size(), ng, MFInfo(), amr_level.Factory()); if (n > 0) { @@ -3944,6 +4033,15 @@ NavierStokesBase::post_timestep_particle (int crse_iteration) timestamp_add_extras(lev, curr_time, tmf); } } + else + { + // + // Timestamp indexes tmf's BoxArray/DistributionMap even + // when tindices is empty, so tmf must always be defined. + // + tmf.define(S_new.boxArray(), S_new.DistributionMap(), 1, ng, + MFInfo(), amr_level.Factory()); + } NSPC->Timestamp(basename, tmf, lev, curr_time, tindices); } @@ -4004,7 +4102,9 @@ NavierStokesBase::ParticleDerive (const std::string& name, { BoxArray ba = parent->boxArray(lev); - MultiFab temp_dat(ba,parent->DistributionMap(lev),1,0,MFInfo(),Factory()); + // NOTE: must use lev's factory, not this level's. + MultiFab temp_dat(ba,parent->DistributionMap(lev),1,0,MFInfo(), + parent->getLevel(lev).Factory()); trr *= parent->refRatio(lev-1); @@ -4018,22 +4118,28 @@ NavierStokesBase::ParticleDerive (const std::string& name, NSPC->Increment(temp_dat,lev); #ifdef _OPENMP -#pragma omp parallel +#pragma omp parallel if (Gpu::notInLaunchRegion()) #endif - for (MFIter mfi(temp_dat,true); mfi.isValid(); ++mfi) + for (MFIter mfi(temp_dat,TilingIfNotGPU()); mfi.isValid(); ++mfi) { - const FArrayBox& ffab = temp_dat[mfi]; - FArrayBox& cfab = ctemp_dat[mfi]; - const Box& fbx = mfi.tilebox(); - - AMREX_ASSERT(cfab.box() == amrex::coarsen(fbx,trr)); - - for (IntVect p = fbx.smallEnd(); p <= fbx.bigEnd(); fbx.next(p)) + const Box& fbx = mfi.tilebox(); + auto const& ffab = temp_dat.const_array(mfi); + auto const& cfab = ctemp_dat.array(mfi); + const IntVect ratio = trr; + + // NOTE: Increment() fills temp_dat on the device, so this + // accumulation must run there too. The atomic also + // removes the OpenMP race between tiles mapping to + // the same coarse cell. + amrex::ParallelFor(fbx, [ffab, cfab, ratio] + AMREX_GPU_DEVICE (int i, int j, int k) noexcept { - const Real val = ffab(p); - if (val > 0) - cfab(amrex::coarsen(p,trr)) += val; - } + const Real val = ffab(i,j,k); + if (val > 0) { + const auto cp = amrex::coarsen(Dim3{i,j,k}, ratio); + Gpu::Atomic::AddNoRet(&cfab(cp.x,cp.y,cp.z), val); + } + }); } temp_dat.clear(); @@ -4874,6 +4980,13 @@ NavierStokesBase::ComputeAofs ( MultiFab& advc, int a_comp, // Advection term "A rstate_tmp.define(S.boxArray(),S.DistributionMap(),ncomp,S.nGrow(), MFInfo(),ebfact); MultiFab::Copy(rstate_tmp,advc,a_comp,0,ncomp,S.nGrow()); + // + // Vsync/Ssync only hold valid-region data (reflux writes valid cells + // only), so fill the grid-overlap ghost cells before using this as the + // state redistribution "state". Otherwise the sync correction near + // every box boundary depends on the grid decomposition. + // + rstate_tmp.FillBoundary(geom.periodicity()); } MultiFab const* rstate = (is_sync && redistribution_type == "StateRedist") ? &rstate_tmp : &S; @@ -4930,6 +5043,14 @@ NavierStokesBase::ComputeAofs ( MultiFab& advc, int a_comp, // Advection term "A FArrayBox fab_drho_as_crse(Box::TheUnitBox(),ncomp); IArrayBox fab_rrflag_as_crse(Box::TheUnitBox()); + // This is a hack-y way of testing whether this ComputeAofs call + // came from the mac_sync (do_crse_add = false) + // or from the regular advance (do_crse_add = true). When the call + // comes from the mac_sync, the multiplier in FineAdd needs to have + // the opposite sign -- and so does the re-redistributed mass that + // the redistribution puts into dm_as_fine/drho_as_crse. + Real sync_factor = do_crse_add ? 1.0 : -1.0; + if (flagfab.getType(grow(bx,4)) != FabType::regular) { AMREX_D_TERM( auto apx = ebfact.getAreaFrac()[0]->const_array(mfi);, @@ -4998,7 +5119,9 @@ NavierStokesBase::ComputeAofs ( MultiFab& advc, int a_comp, // Advection term "A geom, dt, redistribution_type, as_crse, p_drho_as_crse->array(), p_rrflag_as_crse->array(), as_fine, dm_as_fine.array(), coarse_fine_mask->const_array(mfi), - level_mask_notcovered, use_wts_in_divnc); + level_mask_notcovered, + /*fac_for_deltaR*/ sync_factor, + use_wts_in_divnc); } else { bool use_wts_in_divnc = true; ApplyRedistribution( bx, ncomp, redist_arr, update_arr, @@ -5043,6 +5166,63 @@ NavierStokesBase::ComputeAofs ( MultiFab& advc, int a_comp, // Advection term "A FArrayBox fy_fr_fab(fy_fab,amrex::make_alias,flux_comp,ncomp);, FArrayBox fz_fr_fab(fz_fab,amrex::make_alias,flux_comp,ncomp);); + // + // The cut-cell overloads of CrseAdd/FineAdd multiply the fluxes by + // the EB area fraction themselves, but the fluxes computed by + // AMReX-Hydro are already area-fraction weighted. Un-weight them + // here so cut faces are not refluxed with ap^2. + // + AMREX_D_TERM(FArrayBox fx_unwtd_fab;, + FArrayBox fy_unwtd_fab;, + FArrayBox fz_unwtd_fab;); + + const bool need_unwtd_flux = + ( do_reflux && + ( (do_crse_add && level < parent->finestLevel()) || + (do_fine_add && level > 0) ) && + flagfab.getType(amrex::grow(bx,1)) != FabType::regular ); + + if ( need_unwtd_flux ) + { + AMREX_D_TERM(const FArrayBox& apx_fab = (*areafrac[0])[mfi];, + const FArrayBox& apy_fab = (*areafrac[1])[mfi];, + const FArrayBox& apz_fab = (*areafrac[2])[mfi];); + + AMREX_D_TERM(fx_unwtd_fab.resize(fx_fr_fab.box(),ncomp,The_Async_Arena());, + fy_unwtd_fab.resize(fy_fr_fab.box(),ncomp,The_Async_Arena());, + fz_unwtd_fab.resize(fz_fr_fab.box(),ncomp,The_Async_Arena());); + AMREX_D_TERM(fx_unwtd_fab.template setVal(0.0);, + fy_unwtd_fab.template setVal(0.0);, + fz_unwtd_fab.template setVal(0.0);); + + AMREX_D_TERM(auto const& fxu = fx_unwtd_fab.array();, + auto const& fyu = fy_unwtd_fab.array();, + auto const& fzu = fz_unwtd_fab.array();); + AMREX_D_TERM(auto const& fxw = fx_fr_fab.const_array();, + auto const& fyw = fy_fr_fab.const_array();, + auto const& fzw = fz_fr_fab.const_array();); + AMREX_D_TERM(auto const& apx_a = apx_fab.const_array();, + auto const& apy_a = apy_fab.const_array();, + auto const& apz_a = apz_fab.const_array();); + + AMREX_D_TERM(const Box& uxbx = fx_unwtd_fab.box() & apx_fab.box();, + const Box& uybx = fy_unwtd_fab.box() & apy_fab.box();, + const Box& uzbx = fz_unwtd_fab.box() & apz_fab.box();); + + amrex::ParallelFor(uxbx, ncomp, [fxu,fxw,apx_a] + AMREX_GPU_DEVICE (int i, int j, int k, int n) noexcept + { fxu(i,j,k,n) = (apx_a(i,j,k) > 0.) ? fxw(i,j,k,n)/apx_a(i,j,k) : 0.; }); + + amrex::ParallelFor(uybx, ncomp, [fyu,fyw,apy_a] + AMREX_GPU_DEVICE (int i, int j, int k, int n) noexcept + { fyu(i,j,k,n) = (apy_a(i,j,k) > 0.) ? fyw(i,j,k,n)/apy_a(i,j,k) : 0.; }); +#if (AMREX_SPACEDIM == 3) + amrex::ParallelFor(uzbx, ncomp, [fzu,fzw,apz_a] + AMREX_GPU_DEVICE (int i, int j, int k, int n) noexcept + { fzu(i,j,k,n) = (apz_a(i,j,k) > 0.) ? fzw(i,j,k,n)/apz_a(i,j,k) : 0.; }); +#endif + } + // Now update the flux registers (inside test on AMREX_USE_EB) if ( do_reflux && do_crse_add && (level < parent->finestLevel()) ) { if (flagfab.getType(amrex::grow(bx,1)) == FabType::regular) @@ -5053,20 +5233,13 @@ NavierStokesBase::ComputeAofs ( MultiFab& advc, int a_comp, // Advection term "A } else if (flagfab.getType(bx) != FabType::covered ) { getAdvFluxReg(level + 1).CrseAdd(mfi, - {AMREX_D_DECL(&fx_fr_fab,&fy_fr_fab,&fz_fr_fab)}, + {AMREX_D_DECL(&fx_unwtd_fab,&fy_unwtd_fab,&fz_unwtd_fab)}, dxDp, dt, (*volfrac)[mfi], {AMREX_D_DECL(&(*areafrac[0])[mfi], &(*areafrac[1])[mfi], &(*areafrac[2])[mfi])}, 0, state_indx, ncomp, amrex::RunOn::Device); } } // do_reflux && level < finest_level - // This is a hack-y way of testing whether this ComputeAofs call - // came from the mac_sync (do_crse_add = false) - // or from the regular advance (do_crse_add = true). When the call - // comes from the mac_sync, the multiplier in FineAdd needs to have - // the opposite sign - Real sync_factor = do_crse_add ? 1.0 : -1.0; - if ( do_reflux && do_fine_add && (level > 0)) { if (flagfab.getType(amrex::grow(bx,1)) == FabType::regular) { @@ -5075,7 +5248,7 @@ NavierStokesBase::ComputeAofs ( MultiFab& advc, int a_comp, // Advection term "A dxDp, sync_factor*dt, 0, state_indx, ncomp, amrex::RunOn::Device); } else if (flagfab.getType(bx) != FabType::covered ) { advflux_reg->FineAdd(mfi, - {AMREX_D_DECL(&fx_fr_fab,&fy_fr_fab,&fz_fr_fab)}, + {AMREX_D_DECL(&fx_unwtd_fab,&fy_unwtd_fab,&fz_unwtd_fab)}, dxDp, sync_factor*dt, (*volfrac)[mfi], {AMREX_D_DECL(&(*areafrac[0])[mfi], &(*areafrac[1])[mfi], &(*areafrac[2])[mfi])}, dm_as_fine, 0, state_indx, ncomp, amrex::RunOn::Device); diff --git a/Source/Projection.cpp b/Source/Projection.cpp index 3879791e1..78f6d1f6b 100644 --- a/Source/Projection.cpp +++ b/Source/Projection.cpp @@ -1292,8 +1292,12 @@ Projection::scaleVar (MultiFab* sig, amrex::ParallelFor(bx, AMREX_SPACEDIM, [=] AMREX_GPU_DEVICE (int i, int j, int k, int n) noexcept { - if ( i >= domlox && i <= domhix && - j >= domloy && j <= domhiy) + // NOTE: cells outside the domain in the axial (j) direction are + // multiplied by the radius rather than zeroed, because the divu + // stencil in the nodal solver includes them and they might hold + // inflow values. set_boundary_velocity() later zeroes the ghost + // cells that really must be zero. + if ( i >= domlox && i <= domhix ) { velarr(i,j,k,n) = (static_cast(i)+ Real(0.5))*dxr*velarr(i,j,k,n); } @@ -1411,11 +1415,17 @@ Projection::rescaleVar (MultiFab* sig, amrex::ParallelFor(bx, AMREX_SPACEDIM, [=] AMREX_GPU_DEVICE (int i, int j, int k, int n) noexcept { - if ( i >= domlox && i <= domhix && - j >= domloy && j <= domhiy) + // Mirrors scaleVar: the axial out-of-domain ghosts were scaled + // there, so undo the scaling here rather than trashing them. + if ( i >= domlox && i <= domhix ) { velarr(i,j,k,n) /= (static_cast(i)+ Real(0.5))*dxr; } + else if ( n == 0 && i > domhix ) + { + // high-r inflow, scaled in scaleVar + velarr(i,j,k,n) /= (static_cast(i)+ Real(0.5))*dxr; + } else { // set vals outside the domain @@ -2305,7 +2315,7 @@ Projection::computeRhoG(FArrayBox* rhoFab, for (int k = hi.z-1; k >= lo.z; k--) { rho_i = 0.5 * (rho(i,j-1,k) + rho(i-1,j-1,k)); - rho_ii = 0.5 * (rho(i,j-1,k) + rho(i-1,j-2,k)); + rho_ii = 0.5 * (rho(i,j-2,k) + rho(i-1,j-2,k)); add_rhog(rho_i, rho_ii, rhog, phi(i,j,k)); } } @@ -2320,7 +2330,7 @@ Projection::computeRhoG(FArrayBox* rhoFab, if ( has_extdir_lo ) { for (int k = hi.z-1; k >= lo.z; k--) { rho_i = rho(i-1,j-1,k); - rho_ii = rho(i-2,j-1,k); + rho_ii = rho(i-1,j-2,k); add_rhog(rho_i, rho_ii, rhog, phi(i,j,k)); } } else if ( has_hoextrap_lo ) { diff --git a/Source/SyncRegister.cpp b/Source/SyncRegister.cpp index e96d19441..79090b28d 100644 --- a/Source/SyncRegister.cpp +++ b/Source/SyncRegister.cpp @@ -64,7 +64,9 @@ SyncRegister::InitRHS (MultiFab& rhs, const Geometry& geom, const BCRec& phys_bc const int* phys_lo = phys_bc.lo(); const int* phys_hi = phys_bc.hi(); - int outflow_dirs[AMREX_SPACEDIM-1]={-1}; + // One entry per coordinate direction: the loop below appends for any axis + // with an outflow face on either side. + int outflow_dirs[AMREX_SPACEDIM]={AMREX_D_DECL(-1,-1,-1)}; int nOutflow = 0; for (int dir = 0; dir < AMREX_SPACEDIM; dir++) { @@ -262,7 +264,10 @@ SyncRegister::InitRHS (MultiFab& rhs, const Geometry& geom, const BCRec& phys_bc FArrayBox& fab = fs[fsi]; const Box& bx = fab.box(); auto const& mask = fab.array(); - const Real maxcount = AMREX_D_TERM(AMREX_SPACEDIM,*AMREX_SPACEDIM,*AMREX_SPACEDIM) - 0.5; + // The sum built above counts the 2^AMREX_SPACEDIM cells touching a + // node, so the threshold is 2^AMREX_SPACEDIM - 0.5 (3.5 in 2D, + // 7.5 in 3D), as in the Fortran convertmask this replaced. + const Real maxcount = AMREX_D_TERM(2.,*2.,*2.) - 0.5; amrex::ParallelFor(bx, [mask,maxcount] AMREX_GPU_DEVICE (int i, int j, int k) noexcept { diff --git a/Util/ConvertCheckpoint/ConvertCheckpointGrids.cpp b/Util/ConvertCheckpoint/ConvertCheckpointGrids.cpp index 193e3ffae..3a2ae1deb 100644 --- a/Util/ConvertCheckpoint/ConvertCheckpointGrids.cpp +++ b/Util/ConvertCheckpoint/ConvertCheckpointGrids.cpp @@ -65,8 +65,6 @@ Real avg_time_fluct; bool TimeAverageFile_exist = false; int flag_eb = 0; -Vector nsets_save(1); - VisMF::How how = VisMF::OneFilePerCPU; // --------------------------------------------------------------- @@ -234,8 +232,6 @@ static void ReadCheckpointFile(const std::string& fileName) { is >> fakeAmr_src.level_count[i]; } - int ndesc_save; - // READ LEVEL DATA for(int lev(0); lev <= fakeAmr_src.finest_level; ++lev) { @@ -274,12 +270,6 @@ static void ReadCheckpointFile(const std::string& fileName) { is >> nstate; int ndesc = nstate; - // This should be the same at all levels - ndesc_save = ndesc; - - // ndesc depends on which descriptor so we store a value for each - if (lev == 1) nsets_save.resize(ndesc_save); - falRef.state.resize(ndesc); falRef.new_state.resize(ndesc); @@ -296,8 +286,6 @@ static void ReadCheckpointFile(const std::string& fileName) { int nsets; is >> nsets; - nsets_save[ii] = nsets; - falRef.state[ii].old_data = 0; falRef.state[ii].new_data = 0; @@ -386,7 +374,7 @@ static void WriteCheckpointFile(const std::string& inFileName, const std::string HeaderFile.rdbuf()->pubsetbuf(io_buffer.dataPtr(), io_buffer.size()); - int old_prec(0), i, ndesc_save; + int old_prec(0), i; if(ParallelDescriptor::IOProcessor()) { // Only the IOProcessor() writes to the header file. @@ -430,7 +418,6 @@ static void WriteCheckpointFile(const std::string& inFileName, const std::string std::ostream &os = HeaderFile; FakeAmrLevel &falRef = fakeAmr_trgt.fakeAmrLevels[lev]; int ndesc = falRef.state.size(); - ndesc_save = ndesc; // Build directory to hold the MultiFabs in the StateData at this level. char buf[64]; @@ -479,10 +466,17 @@ static void WriteCheckpointFile(const std::string& inFileName, const std::string const std::string name(PathNameInHeader); const std::string fullpathname(FullPathName); + // + // Derive the number of data sets from the pointers themselves, the + // way StateData::checkPoint does. This covers levels (and state + // types) whose nsets differ, and single-level checkpoints. + // bool dump_old(true); if(dump_old == true && falRef.state[i].old_data == 0) { dump_old = false; } + const int nsets_loc = (falRef.state[i].new_data == 0) ? 0 + : (dump_old ? 2 : 1); if(ParallelDescriptor::IOProcessor()) { // The relative name gets written to the Header file. @@ -501,7 +495,7 @@ static void WriteCheckpointFile(const std::string& inFileName, const std::string << falRef.state[i].new_time.start << '\n' << falRef.state[i].new_time.stop << '\n'; - if (nsets_save[i] > 0) { + if (nsets_loc > 0) { if(dump_old) { os << 2 << '\n' << mf_name_new << '\n' << mf_name_old << '\n'; } else { @@ -513,14 +507,14 @@ static void WriteCheckpointFile(const std::string& inFileName, const std::string } - if (nsets_save[i] > 0) { + if (nsets_loc > 0) { BL_ASSERT(falRef.state[i].new_data); std::string mf_fullpath_new = fullpathname; mf_fullpath_new += NewSuffix; VisMF::Write(*(falRef.state[i].new_data),mf_fullpath_new,how); } - if (nsets_save[i] > 1) { + if (nsets_loc > 1) { BL_ASSERT(dump_old); BL_ASSERT(falRef.state[i].old_data); std::string mf_fullpath_old = fullpathname; @@ -654,21 +648,22 @@ static void ConvertData() { if(n == 1) ngrow_loc = 1; } - // We don't have the same number of ghost-cells for each data type - // Warning, this should be adapted for EB - if (falRef_src.state.size() == 4 && n == falRef_src.state.size()-1){ - ngrow_loc = 0; // For this case, we just have Average_Type and no Divu_Type and Dsdt_type - } - else if (falRef_src.state.size() == 5 && n == falRef_src.state.size()-1){ - ngrow_loc = 0; // For this case, we have Divu_Type and Dsdt_type, no Average_Type - } - else if (falRef_src.state.size() == 6 && n >= falRef_src.state.size()-2){ - ngrow_loc = 0; // Here we have both Average_Type and Divu and Dsdt types - } - + // + // NOTE: every cell-centered state needs at least one ghost cell of + // workspace here: CellConservativeLinear::CoarseBox grows + // coarsen(fine_box) by one, and the slope kernels read i+/-1. + // The copies below use src_nghost=0, so the source's own ghost + // count does not matter. + // // Assuming that OldState and NewState have the same number of components - int ncomps = (falRef_src.state[n].old_data)->nComp(); + int ncomps = (falRef_src.state[n].new_data)->nComp(); + + // + // old_data is null when the checkpoint stores only one data set + // (nsets == 1), e.g. a checkpoint written right after initialization. + // + bool has_old = (falRef_src.state[n].old_data != 0); BoxArray new_grids_state = falRef_trgt.state[n].grids; BoxArray save_grids_state = falRef_trgt.state[n].grids; @@ -694,7 +689,9 @@ static void ConvertData() { OldData_src -> setVal(10.); NewData_src -> copy(*(falRef_src.state[n].new_data),0,0,ncomps,0,ngrow_loc); - OldData_src -> copy(*(falRef_src.state[n].old_data),0,0,ncomps,0,ngrow_loc); + if (has_old) { + OldData_src -> copy(*(falRef_src.state[n].old_data),0,0,ncomps,0,ngrow_loc); + } MultiFab * NewData_trgt = new MultiFab(new_grids_state,dm_trgt,ncomps,ngrow_loc); MultiFab * OldData_trgt = new MultiFab(new_grids_state,dm_trgt,ncomps,ngrow_loc); @@ -712,6 +709,17 @@ static void ConvertData() { const Geometry& fgeom = falRef_trgt.geom; const Geometry& cgeom = falRef_src.geom; + // + // The copies above are non-periodic, so ghost cells at periodic + // domain boundaries still hold the setVal(10.) sentinel. Fill them + // from valid data before computing the interpolation slopes. + // NOTE: ghosts at non-periodic physical boundaries still hold the + // sentinel -- the state descriptors (and hence the real BCs) + // are never restored by this tool. + // + NewData_src->FillBoundary(cgeom.periodicity()); + OldData_src->FillBoundary(cgeom.periodicity()); + for (MFIter mfi(*NewData_trgt); mfi.isValid(); ++mfi) { FArrayBox& ffab = (*NewData_trgt)[mfi]; @@ -749,7 +757,11 @@ static void ConvertData() { } falRef_trgt.state[n].new_data = NewData_trgt; - falRef_trgt.state[n].old_data = OldData_trgt; + // + // Leave old_data null when the source had none, so that the writer + // emits a single data set (nsets == 1) for this state. + // + falRef_trgt.state[n].old_data = has_old ? OldData_trgt : 0; } } diff --git a/Util/ConvertCheckpoint/Make.package b/Util/ConvertCheckpoint/Make.package index 640584d46..b642d1309 100644 --- a/Util/ConvertCheckpoint/Make.package +++ b/Util/ConvertCheckpoint/Make.package @@ -1,5 +1,4 @@ CEXE_headers += AMReX_DataServices.H AMReX_AmrData.H AMReX_XYPlotDataList.H AMReX_AmrvisConstants.H CEXE_sources += AMReX_DataServices.cpp AMReX_AmrData.cpp AMReX_XYPlotDataList.cpp -FEXE_sources += AMReX_FABUTIL_$(DIM)D.F