diff --git a/Source/BoundaryConditions/ERF_FillPatch.cpp b/Source/BoundaryConditions/ERF_FillPatch.cpp index 4a98a01ba..6836eba2e 100644 --- a/Source/BoundaryConditions/ERF_FillPatch.cpp +++ b/Source/BoundaryConditions/ERF_FillPatch.cpp @@ -43,20 +43,9 @@ ERF::FillPatch (int lev, Real time, FPr_c[lev-1].FillSet(*mfs_vel[Vars::cons], time, null_bc, domain_bcs_type); } if (cf_set_width >= 0 && !cons_only) { - // - // This is an optimization since we won't need more than one ghost - // cell of momentum in the integrator if not using NumDiff - // - //IntVect ngu = (solverChoice.use_NumDiff) ? IntVect(1,1,1) : mfs_vel[Vars::xvel]->nGrowVect(); - //IntVect ngv = (solverChoice.use_NumDiff) ? IntVect(1,1,1) : mfs_vel[Vars::yvel]->nGrowVect(); - //IntVect ngw = (solverChoice.use_NumDiff) ? IntVect(1,1,0) : mfs_vel[Vars::zvel]->nGrowVect(); - IntVect ngu = IntVect::TheZeroVector(); - IntVect ngv = IntVect::TheZeroVector(); - IntVect ngw = IntVect::TheZeroVector(); - - VelocityToMomentum(*mfs_vel[Vars::xvel], ngu, - *mfs_vel[Vars::yvel], ngv, - *mfs_vel[Vars::zvel], ngw, + VelocityToMomentum(*mfs_vel[Vars::xvel], IntVect{0}, + *mfs_vel[Vars::yvel], IntVect{0}, + *mfs_vel[Vars::zvel], IntVect{0}, *mfs_vel[Vars::cons], *mfs_mom[IntVars::xmom], *mfs_mom[IntVars::ymom], @@ -545,16 +534,12 @@ ERF::FillIntermediatePatch (int lev, Real time, // We always come in to this call with momenta so we need to leave with momenta! // We need to make sure we convert back on all ghost cells/faces because this is // how velocity from fine-fine copies (as well as physical and interpolated bcs) will be filled - if (!cons_only) { - IntVect ngu = mfs_vel[Vars::xvel]->nGrowVect(); - IntVect ngv = mfs_vel[Vars::yvel]->nGrowVect(); - IntVect ngw = mfs_vel[Vars::zvel]->nGrowVect(); - - if (!solverChoice.use_NumDiff) { - ngu = IntVect(1,1,1); - ngv = IntVect(1,1,1); - ngw = IntVect(1,1,1); - } + if (!cons_only) + { + IntVect ngu = (!solverChoice.use_NumDiff) ? IntVect(1,1,1) : mfs_vel[Vars::xvel]->nGrowVect(); + IntVect ngv = (!solverChoice.use_NumDiff) ? IntVect(1,1,1) : mfs_vel[Vars::yvel]->nGrowVect(); + IntVect ngw = (!solverChoice.use_NumDiff) ? IntVect(1,1,0) : mfs_vel[Vars::zvel]->nGrowVect(); + VelocityToMomentum(*mfs_vel[Vars::xvel], ngu, *mfs_vel[Vars::yvel], ngv, *mfs_vel[Vars::zvel], ngw, @@ -600,13 +585,9 @@ ERF::FillCoarsePatch (int lev, Real time) // Convert velocity to momentum at the COARSE level // ************************************************ // - IntVect ngu = IntVect(0,0,0); - IntVect ngv = IntVect(0,0,0); - IntVect ngw = IntVect(0,0,0); - - VelocityToMomentum(vars_new[lev-1][Vars::xvel], ngu, - vars_new[lev-1][Vars::yvel], ngv, - vars_new[lev-1][Vars::zvel], ngw, + VelocityToMomentum(vars_new[lev-1][Vars::xvel], IntVect{0}, + vars_new[lev-1][Vars::yvel], IntVect{0}, + vars_new[lev-1][Vars::zvel], IntVect{0}, vars_new[lev-1][Vars::cons], rU_new[lev-1], rV_new[lev-1], @@ -640,7 +621,7 @@ ERF::FillCoarsePatch (int lev, Real time) // with InterpFromCoarseLevel which ASSUMES that all ghost cells have already been filled // ************************************************************************************************ // - InterpFromCoarseLevel(rU_new[lev], ngu, IntVect(0,0,0), rU_new[lev-1], 0, 0, 1, + InterpFromCoarseLevel(rU_new[lev], IntVect{0}, IntVect{0}, rU_new[lev-1], 0, 0, 1, geom[lev-1], geom[lev], refRatio(lev-1), mapper_f, domain_bcs_type, BCVars::xvel_bc); @@ -650,7 +631,7 @@ ERF::FillCoarsePatch (int lev, Real time) // with InterpFromCoarseLevel which ASSUMES that all ghost cells have already been filled // ************************************************************************************************ // - InterpFromCoarseLevel(rV_new[lev], ngv, IntVect(0,0,0), rV_new[lev-1], 0, 0, 1, + InterpFromCoarseLevel(rV_new[lev], IntVect{0}, IntVect{0}, rV_new[lev-1], 0, 0, 1, geom[lev-1], geom[lev], refRatio(lev-1), mapper_f, domain_bcs_type, BCVars::yvel_bc); @@ -658,7 +639,7 @@ ERF::FillCoarsePatch (int lev, Real time) // Interpolate z-momentum from coarse to fine level // with InterpFromCoarseLevel which ASSUMES that all ghost cells have already been filled // ************************************************************************************************ - InterpFromCoarseLevel(rW_new[lev], ngw, IntVect(0,0,0), rW_new[lev-1], 0, 0, 1, + InterpFromCoarseLevel(rW_new[lev], IntVect{0}, IntVect{0}, rW_new[lev-1], 0, 0, 1, geom[lev-1], geom[lev], refRatio(lev-1), mapper_f, domain_bcs_type, BCVars::zvel_bc); // diff --git a/Source/Diffusion/ERF_ComputeStrain_N.cpp b/Source/Diffusion/ERF_ComputeStrain_N.cpp index 996851723..6fbafdd5c 100644 --- a/Source/Diffusion/ERF_ComputeStrain_N.cpp +++ b/Source/Diffusion/ERF_ComputeStrain_N.cpp @@ -208,8 +208,8 @@ ComputeStrain_N (Box bxcc, Box tbxxy, Box tbxxz, Box tbxyz, Box domain, //*********************************************************************************** // Cell centered strains ParallelFor(bxcc, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept { - tau11(i,j,k) = (u(i+1, j , k )/mf_u(i+1,j,0) - u(i, j, k)/mf_u(i,j,0))*dxInv[0]*mf_u(i,j,0)*mf_u(i,j,0); - tau22(i,j,k) = (v(i , j+1, k )/mf_v(i,j+1,0) - v(i, j, k)/mf_v(i,j,0))*dxInv[1]*mf_v(i,j,0)*mf_v(i,j,0); + tau11(i,j,k) = (u(i+1, j , k )/mf_u(i+1,j,0) - u(i, j, k)/mf_u(i,j,0))*dxInv[0]*mf_m(i,j,0)*mf_m(i,j,0); + tau22(i,j,k) = (v(i , j+1, k )/mf_v(i,j+1,0) - v(i, j, k)/mf_v(i,j,0))*dxInv[1]*mf_m(i,j,0)*mf_m(i,j,0); tau33(i,j,k) = (w(i , j , k+1) - w(i, j, k))*dxInv[2]; }); diff --git a/Source/Diffusion/ERF_ComputeStrain_T.cpp b/Source/Diffusion/ERF_ComputeStrain_T.cpp index 450d4317d..377a6a782 100644 --- a/Source/Diffusion/ERF_ComputeStrain_T.cpp +++ b/Source/Diffusion/ERF_ComputeStrain_T.cpp @@ -25,6 +25,7 @@ using namespace amrex; * @param[in] z_nd nodal array of physical z heights * @param[in] bc_ptr container with boundary condition types * @param[in] dxInv inverse cell size array + * @param[in] mf_m map factor at cell center * @param[in] mf_u map factor at x-face * @param[in] mf_v map factor at y-face */ @@ -38,7 +39,7 @@ ComputeStrain_T (Box bxcc, Box tbxxy, Box tbxxz, Box tbxyz, Box domain, const Array4& z_nd, const Array4& detJ, const BCRec* bc_ptr, const GpuArray& dxInv, - const Array4& /*mf_m*/, + const Array4& mf_m, const Array4& mf_u, const Array4& mf_v) { @@ -452,9 +453,9 @@ ComputeStrain_T (Box bxcc, Box tbxxy, Box tbxxz, Box tbxyz, Box domain, met_h_zeta = detJ(i,j,k); tau11(i,j,k) = ( (u(i+1, j, k)/mf_u(i+1,j,0) - u(i, j, k)/mf_u(i,j,0))*dxInv[0] - - (met_h_xi/met_h_zeta)*GradUz ) * mf_u(i,j,0)*mf_u(i,j,0); + - (met_h_xi/met_h_zeta)*GradUz ) * mf_m(i,j,0)*mf_m(i,j,0); tau22(i,j,k) = ( (v(i, j+1, k)/mf_v(i,j+1,0) - v(i, j, k)/mf_v(i,j,0))*dxInv[1] - - (met_h_eta/met_h_zeta)*GradVz ) * mf_v(i,j,0)*mf_v(i,j,0); + - (met_h_eta/met_h_zeta)*GradVz ) * mf_m(i,j,0)*mf_m(i,j,0); tau33(i,j,k) = (w(i, j, k+1) - w(i, j, k))*dxInv[2]/met_h_zeta; }); diff --git a/Source/ERF_make_new_arrays.cpp b/Source/ERF_make_new_arrays.cpp index ed9f652bd..68630c53a 100644 --- a/Source/ERF_make_new_arrays.cpp +++ b/Source/ERF_make_new_arrays.cpp @@ -147,6 +147,9 @@ ERF::init_stuff (int lev, const BoxArray& ba, const DistributionMapping& dm, rW_new[lev].define(convert(ba, IntVect(0,0,1)), dm, 1, ngrow_vels); // We do this here just so they won't be undefined in the initial FillPatch + rU_old[lev].setVal(1.2e21); + rV_old[lev].setVal(3.4e22); + rW_old[lev].setVal(5.6e23); rU_new[lev].setVal(1.2e21); rV_new[lev].setVal(3.4e22); rW_new[lev].setVal(5.6e23); diff --git a/Source/ERF_make_new_level.cpp b/Source/ERF_make_new_level.cpp index 2de79276f..52d59ba27 100644 --- a/Source/ERF_make_new_level.cpp +++ b/Source/ERF_make_new_level.cpp @@ -26,7 +26,7 @@ void ERF::MakeNewLevelFromScratch (int lev, Real time, const BoxArray& ba_in, BoxArray ba; DistributionMapping dm; Box domain(Geom(0).Domain()); - if (lev == 0 && + if (lev == 0 && restart_chkfile.empty() && (max_grid_size[0][0] >= domain.length(0)) && (max_grid_size[0][1] >= domain.length(1)) && ba_in.size() != ParallelDescriptor::NProcs()) diff --git a/Source/TimeIntegration/ERF_advance_dycore.cpp b/Source/TimeIntegration/ERF_advance_dycore.cpp index 2872007e8..4761c8646 100644 --- a/Source/TimeIntegration/ERF_advance_dycore.cpp +++ b/Source/TimeIntegration/ERF_advance_dycore.cpp @@ -252,9 +252,10 @@ void ERF::advance_dycore(int level, // This is an optimization since we won't need more than one ghost // cell of momentum in the integrator if not using NumDiff // - IntVect ngu = (solverChoice.use_NumDiff) ? IntVect(1,1,1) : xvel_old.nGrowVect(); - IntVect ngv = (solverChoice.use_NumDiff) ? IntVect(1,1,1) : yvel_old.nGrowVect(); - IntVect ngw = (solverChoice.use_NumDiff) ? IntVect(1,1,0) : zvel_old.nGrowVect(); + IntVect ngu = (!solverChoice.use_NumDiff) ? IntVect(1,1,1) : xvel_old.nGrowVect(); + IntVect ngv = (!solverChoice.use_NumDiff) ? IntVect(1,1,1) : yvel_old.nGrowVect(); + IntVect ngw = (!solverChoice.use_NumDiff) ? IntVect(1,1,0) : zvel_old.nGrowVect(); + VelocityToMomentum(xvel_old, ngu, yvel_old, ngv, zvel_old, ngw, density, state_old[IntVars::xmom], state_old[IntVars::ymom], diff --git a/Source/Utils/ERF_VelocityToMomentum.cpp b/Source/Utils/ERF_VelocityToMomentum.cpp index 271356276..213572e41 100644 --- a/Source/Utils/ERF_VelocityToMomentum.cpp +++ b/Source/Utils/ERF_VelocityToMomentum.cpp @@ -54,19 +54,9 @@ void VelocityToMomentum (const MultiFab& xvel_in, tby = mfi.tilebox(IntVect(0,1,0),yvel_ngrow); tbz = mfi.tilebox(IntVect(0,0,1),zvel_ngrow); -#if 0 - if (l_use_ndiff) { - tbx = mfi.tilebox(IntVect(1,0,0),xvel_ngrow); - tby = mfi.tilebox(IntVect(0,1,0),yvel_ngrow); - tbz = mfi.tilebox(IntVect(0,0,1),zvel_ngrow); - } else { - tbx = mfi.tilebox(IntVect(1,0,0),IntVect(1,1,1)); - if (tbx.smallEnd(2) < 0) tbx.setSmall(2,0); - tby = mfi.tilebox(IntVect(0,1,0),IntVect(1,1,1)); - if (tby.smallEnd(2) < 0) tby.setSmall(2,0); - tbz = mfi.tilebox(IntVect(0,0,1),IntVect(1,1,0)); - } -#endif + // Don't actually try to fill w above or below the domain + if (tbz.smallEnd(2) < domain.smallEnd(2)) tbz.setSmall(2,domain.smallEnd(2)); + if (tbz.bigEnd(2) > domain.bigEnd(2)+1) tbz.setBig(2,domain.bigEnd(2)+1); // Conserved/state variables on cell centers -- we use this for density const Array4& dens_arr = density.array(mfi); @@ -81,6 +71,8 @@ void VelocityToMomentum (const MultiFab& xvel_in, const Array4& vely = yvel_in.const_array(mfi); const Array4& velz = zvel_in.const_array(mfi); + // ******************************************************************************************** + ParallelFor(tbx, tby, tbz, [=] AMREX_GPU_DEVICE (int i, int j, int k) { momx(i,j,k) = velx(i,j,k) * 0.5 * (dens_arr(i,j,k,Rho_comp) + dens_arr(i-1,j,k,Rho_comp)); @@ -92,6 +84,8 @@ void VelocityToMomentum (const MultiFab& xvel_in, momz(i,j,k) = velz(i,j,k) * 0.5 * (dens_arr(i,j,k,Rho_comp) + dens_arr(i,j,k-1,Rho_comp)); }); + // ******************************************************************************************** + if ( (bx.smallEnd(0) == domain.smallEnd(0)) && (bc_ptr_h[BCVars::cons_bc].lo(0) == ERFBCType::ext_dir) ) { ParallelFor(makeSlab(tbx,0,domain.smallEnd(0)), [=] AMREX_GPU_DEVICE (int i, int j, int k) {