diff --git a/Source/Diffusion.H b/Source/Diffusion.H index 78c2a8f5..b44157b1 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 b9c8fc92..13250590 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 3c48a0f4..3f647492 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 80c3282a..e6be18a5 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 4081d8d2..2edc390e 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 0ac31d60..bc4ce8a7 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 6d71f67b..f1920d7b 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 3879791e..78f6d1f6 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 e96d1944..79090b28 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 193e3ffa..3a2ae1de 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 640584d4..b642d130 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