diff --git a/Exec/CMakeLists.txt b/Exec/CMakeLists.txt index ea5fc102d..5246ae6e1 100644 --- a/Exec/CMakeLists.txt +++ b/Exec/CMakeLists.txt @@ -28,13 +28,13 @@ else () add_subdirectory(RegTests/EkmanSpiral_ideal) add_subdirectory(RegTests/EkmanSpiral_input_sounding) add_subdirectory(RegTests/IsentropicVortex) + add_subdirectory(RegTests/ParticlesOverWoA) add_subdirectory(RegTests/PoiseuilleFlow) add_subdirectory(RegTests/ScalarAdvDiff) add_subdirectory(RegTests/TaylorGreenVortex) add_subdirectory(RegTests/WitchOfAgnesi) add_subdirectory(RegTests/WPS_Test) add_subdirectory(DevTests/MovingTerrain) - add_subdirectory(DevTests/ParticlesOverWoA) add_subdirectory(DevTests/MiguelDev) add_subdirectory(DevTests/MetGrid) add_subdirectory(DevTests/LandSurfaceModel) diff --git a/Exec/DevTests/ParticlesOverWoA/CMakeLists.txt b/Exec/RegTests/ParticlesOverWoA/CMakeLists.txt similarity index 100% rename from Exec/DevTests/ParticlesOverWoA/CMakeLists.txt rename to Exec/RegTests/ParticlesOverWoA/CMakeLists.txt diff --git a/Exec/DevTests/ParticlesOverWoA/GNUmakefile b/Exec/RegTests/ParticlesOverWoA/GNUmakefile similarity index 88% rename from Exec/DevTests/ParticlesOverWoA/GNUmakefile rename to Exec/RegTests/ParticlesOverWoA/GNUmakefile index d00e9f848..8bc175ef1 100644 --- a/Exec/DevTests/ParticlesOverWoA/GNUmakefile +++ b/Exec/RegTests/ParticlesOverWoA/GNUmakefile @@ -31,5 +31,5 @@ Bpack := ./Make.package Blocs := . ERF_HOME := ../../.. -ERF_PROBLEM_DIR = $(ERF_HOME)/Exec/DevTests/ParticlesOverWoA +ERF_PROBLEM_DIR = $(ERF_HOME)/Exec/RegTests/ParticlesOverWoA include $(ERF_HOME)/Exec/Make.ERF diff --git a/Exec/DevTests/ParticlesOverWoA/Make.package b/Exec/RegTests/ParticlesOverWoA/Make.package similarity index 100% rename from Exec/DevTests/ParticlesOverWoA/Make.package rename to Exec/RegTests/ParticlesOverWoA/Make.package diff --git a/Exec/DevTests/ParticlesOverWoA/README b/Exec/RegTests/ParticlesOverWoA/README similarity index 100% rename from Exec/DevTests/ParticlesOverWoA/README rename to Exec/RegTests/ParticlesOverWoA/README diff --git a/Exec/DevTests/ParticlesOverWoA/inputs b/Exec/RegTests/ParticlesOverWoA/inputs similarity index 97% rename from Exec/DevTests/ParticlesOverWoA/inputs rename to Exec/RegTests/ParticlesOverWoA/inputs index 69ff18cb1..725789c3f 100644 --- a/Exec/DevTests/ParticlesOverWoA/inputs +++ b/Exec/RegTests/ParticlesOverWoA/inputs @@ -1,5 +1,5 @@ # ------------------ INPUTS TO MAIN PROGRAM ------------------- -max_step = 200 +max_step = 40 amrex.fpe_trap_invalid = 1 @@ -27,7 +27,7 @@ zhi.type = "SlipWall" erf.use_tracer_particles = 1 # TIME STEP CONTROL -erf.fixed_dt = 2E-4 +erf.fixed_dt = 1E-3 # DIAGNOSTICS & VERBOSITY erf.sum_interval = 1 # timesteps between computing mass @@ -81,5 +81,5 @@ erf.alpha_T = 0.0 # [m^2/s] # PROBLEM PARAMETERS (optional) prob.T_0 = 300.0 -prob.U_0 = 0.0 +prob.U_0 = 10.0 prob.rho_0 = 1.16 diff --git a/Exec/DevTests/ParticlesOverWoA/prob.H b/Exec/RegTests/ParticlesOverWoA/prob.H similarity index 100% rename from Exec/DevTests/ParticlesOverWoA/prob.H rename to Exec/RegTests/ParticlesOverWoA/prob.H diff --git a/Exec/DevTests/ParticlesOverWoA/prob.cpp b/Exec/RegTests/ParticlesOverWoA/prob.cpp similarity index 96% rename from Exec/DevTests/ParticlesOverWoA/prob.cpp rename to Exec/RegTests/ParticlesOverWoA/prob.cpp index 19b5ce3bd..422823dde 100644 --- a/Exec/DevTests/ParticlesOverWoA/prob.cpp +++ b/Exec/RegTests/ParticlesOverWoA/prob.cpp @@ -48,6 +48,7 @@ Problem::init_custom_pert( Array4 const& /*mf_v*/, const SolverChoice& sc) { + const int klo = geomdata.Domain().smallEnd()[2]; const int khi = geomdata.Domain().bigEnd()[2]; const bool use_moisture = (sc.moisture_type != MoistureType::None); @@ -94,7 +95,9 @@ Problem::init_custom_pert( // Set the x-velocity amrex::ParallelFor(xbx, [=, parms=parms] AMREX_GPU_DEVICE(int i, int j, int k) noexcept { - x_vel_pert(i, j, k) = parms.U_0; + Real ztop = z_nd(i,j,khi+1); + Real zht = z_nd(i,j,klo); + x_vel_pert(i, j, k) = parms.U_0 * ztop / (ztop - zht); }); // Set the y-velocity diff --git a/Source/BoundaryConditions/ERF_FillPatch.cpp b/Source/BoundaryConditions/ERF_FillPatch.cpp index 63fbaec77..8b92438b4 100644 --- a/Source/BoundaryConditions/ERF_FillPatch.cpp +++ b/Source/BoundaryConditions/ERF_FillPatch.cpp @@ -50,6 +50,8 @@ ERF::FillPatch (int lev, Real time, *mfs_mom[IntVars::xmom], *mfs_mom[IntVars::ymom], *mfs_mom[IntVars::zmom], + Geom(lev).Domain(), + domain_bcs_type, solverChoice.use_NumDiff); FPr_u[lev-1].FillSet(*mfs_mom[IntVars::xmom], time, null_bc, domain_bcs_type); FPr_v[lev-1].FillSet(*mfs_mom[IntVars::ymom], time, null_bc, domain_bcs_type); @@ -59,7 +61,9 @@ ERF::FillPatch (int lev, Real time, *mfs_vel[Vars::cons], *mfs_mom[IntVars::xmom], *mfs_mom[IntVars::ymom], - *mfs_mom[IntVars::zmom]); + *mfs_mom[IntVars::zmom], + Geom(lev).Domain(), + domain_bcs_type); } } @@ -257,7 +261,8 @@ ERF::FillIntermediatePatch (int lev, Real time, // This only fills VALID region of velocity MomentumToVelocity(*mfs_vel[Vars::xvel], *mfs_vel[Vars::yvel], *mfs_vel[Vars::zvel], *mfs_vel[Vars::cons], - *mfs_mom[IntVars::xmom], *mfs_mom[IntVars::ymom], *mfs_mom[IntVars::zmom]); + *mfs_mom[IntVars::xmom], *mfs_mom[IntVars::ymom], *mfs_mom[IntVars::zmom], + Geom(lev).Domain(), domain_bcs_type); } // We now start working on VELOCITY @@ -366,6 +371,8 @@ ERF::FillIntermediatePatch (int lev, Real time, *mfs_vel[Vars::zvel], mfs_vel[Vars::zvel]->nGrowVect(), *mfs_vel[Vars::cons], *mfs_mom[IntVars::xmom], *mfs_mom[IntVars::ymom], *mfs_mom[IntVars::zmom], + Geom(lev).Domain(), + domain_bcs_type, solverChoice.use_NumDiff); } } @@ -417,6 +424,8 @@ ERF::FillCoarsePatch (int lev, Real time) rU_new[which_lev], rV_new[which_lev], rW_new[which_lev], + Geom(lev).Domain(), + domain_bcs_type, true); } @@ -446,7 +455,9 @@ ERF::FillCoarsePatch (int lev, Real time) vars_new[which_lev][Vars::cons], rU_new[which_lev], rV_new[which_lev], - rW_new[which_lev]); + rW_new[which_lev], + Geom(lev).Domain(), + domain_bcs_type); } vars_new[lev][Vars::cons].FillBoundary(geom[lev].periodicity()); diff --git a/Source/ERF.cpp b/Source/ERF.cpp index 72e64034a..5fefe3c5d 100644 --- a/Source/ERF.cpp +++ b/Source/ERF.cpp @@ -1576,6 +1576,8 @@ ERF::AverageDownTo (int crse_lev, int scomp, int ncomp) // NOLINT rU_new[lev], rV_new[lev], rW_new[lev], + Geom(lev).Domain(), + domain_bcs_type, true); } @@ -1590,7 +1592,9 @@ ERF::AverageDownTo (int crse_lev, int scomp, int ncomp) // NOLINT vars_new[lev][Vars::cons], rU_new[lev], rV_new[lev], - rW_new[lev]); + rW_new[lev], + Geom(lev).Domain(), + domain_bcs_type); } } diff --git a/Source/TimeIntegration/ERF_advance_dycore.cpp b/Source/TimeIntegration/ERF_advance_dycore.cpp index 78e30b8c0..0a6dc286f 100644 --- a/Source/TimeIntegration/ERF_advance_dycore.cpp +++ b/Source/TimeIntegration/ERF_advance_dycore.cpp @@ -43,7 +43,6 @@ void ERF::advance_dycore(int level, const Real dt_advance, const Real old_time) { BL_PROFILE_VAR("erf_advance_dycore()",erf_advance_dycore); - if (verbose) Print() << "Starting advance_dycore at level " << level << std::endl; DiffChoice dc = solverChoice.diffChoice; TurbChoice tc = solverChoice.turbChoice[level]; @@ -222,14 +221,15 @@ void ERF::advance_dycore(int level, state_old[IntVars::xmom], state_old[IntVars::ymom], state_old[IntVars::zmom], + Geom(level).Domain(), domain_bcs_type, solverChoice.use_NumDiff); - MultiFab::Copy(xvel_new,xvel_old,0,0,1,xvel_old.nGrowVect()); MultiFab::Copy(yvel_new,yvel_old,0,0,1,yvel_old.nGrowVect()); MultiFab::Copy(zvel_new,zvel_old,0,0,1,zvel_old.nGrowVect()); bool fast_only = false; bool vel_and_mom_synced = true; + apply_bcs(state_old, old_time, state_old[IntVars::cons].nGrow(), state_old[IntVars::xmom].nGrow(), fast_only, vel_and_mom_synced); diff --git a/Source/Utils/MomentumToVelocity.cpp b/Source/Utils/MomentumToVelocity.cpp index ca8d4ca5b..281c9567b 100644 --- a/Source/Utils/MomentumToVelocity.cpp +++ b/Source/Utils/MomentumToVelocity.cpp @@ -17,15 +17,21 @@ using namespace amrex; * @param[in] xmom_in x-component of momentum * @param[in] ymom_in y-component of momentum * @param[in] zmom_in z-component of momentum + * @param[in] domain Domain at this level + * @param[in] domain_bcs_type_h host vector for domain boundary conditions */ void MomentumToVelocity (MultiFab& xvel, MultiFab& yvel, MultiFab& zvel, const MultiFab& density, - const MultiFab& xmom_in, const MultiFab& ymom_in, const MultiFab& zmom_in) + const MultiFab& xmom_in, const MultiFab& ymom_in, const MultiFab& zmom_in, + const Box& domain, + const Vector& domain_bcs_type_h) { BL_PROFILE_VAR("MomentumToVelocity()",MomentumToVelocity); + const BCRec* bc_ptr_h = domain_bcs_type_h.data(); + #ifdef _OPENMP #pragma omp parallel if (amrex::Gpu::notInLaunchRegion()) #endif @@ -53,16 +59,38 @@ MomentumToVelocity (MultiFab& xvel, MultiFab& yvel, MultiFab& zvel, ParallelFor(tbx, tby, tbz, [=] AMREX_GPU_DEVICE (int i, int j, int k) { - Real rho_x_inv = 2.0 / (dens_arr(i,j,k,Rho_comp) + dens_arr(i-1,j,k,Rho_comp)); - velx(i,j,k) = momx(i,j,k) * rho_x_inv; + velx(i,j,k) = momx(i,j,k) * 2.0 / (dens_arr(i,j,k,Rho_comp) + dens_arr(i-1,j,k,Rho_comp)); }, [=] AMREX_GPU_DEVICE (int i, int j, int k) { - Real rho_y_inv = 2.0 / (dens_arr(i,j,k,Rho_comp) + dens_arr(i,j-1,k,Rho_comp)); - vely(i,j,k) = momy(i,j,k) * rho_y_inv; + vely(i,j,k) = momy(i,j,k) * 2.0 / (dens_arr(i,j,k,Rho_comp) + dens_arr(i,j-1,k,Rho_comp)); }, [=] AMREX_GPU_DEVICE (int i, int j, int k) { - Real rho_z_inv = 2.0 / (dens_arr(i,j,k,Rho_comp) + dens_arr(i,j,k-1,Rho_comp)); - velz(i,j,k) = momz(i,j,k) * rho_z_inv; + velz(i,j,k) = momz(i,j,k) * 2.0 / (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) { + velx(i,j,k) = momx(i,j,k) / dens_arr(i-1,j,k,Rho_comp); + }); + } + if ( (bx.bigEnd(0) == domain.bigEnd(0)) && + (bc_ptr_h[BCVars::cons_bc].hi(0) == ERFBCType::ext_dir) ) { + ParallelFor(makeSlab(tbx,0,domain.bigEnd(0)+1), [=] AMREX_GPU_DEVICE (int i, int j, int k) { + velx(i,j,k) = momx(i,j,k) / dens_arr(i,j,k,Rho_comp); + }); + } + if ( (bx.smallEnd(1) == domain.smallEnd(1)) && + (bc_ptr_h[BCVars::cons_bc].lo(1) == ERFBCType::ext_dir) ) { + ParallelFor(makeSlab(tby,1,domain.smallEnd(1)), [=] AMREX_GPU_DEVICE (int i, int j, int k) { + vely(i,j,k) = momy(i,j,k) / dens_arr(i,j-1,k,Rho_comp); + }); + } + if ( (bx.bigEnd(0) == domain.bigEnd(0)) && + (bc_ptr_h[BCVars::cons_bc].hi(0) == ERFBCType::ext_dir) ) { + ParallelFor(makeSlab(tby,1,domain.bigEnd(1)+1), [=] AMREX_GPU_DEVICE (int i, int j, int k) { + vely(i,j,k) = momy(i,j,k) / dens_arr(i,j,k,Rho_comp); + }); + } } // end MFIter } diff --git a/Source/Utils/Utils.H b/Source/Utils/Utils.H index 0058f085f..e3afec3e6 100644 --- a/Source/Utils/Utils.H +++ b/Source/Utils/Utils.H @@ -32,7 +32,9 @@ void MomentumToVelocity (amrex::MultiFab& xvel_out, const amrex::MultiFab& cons_in, const amrex::MultiFab& xmom_in, const amrex::MultiFab& ymom_in, - const amrex::MultiFab& zmom_in); + const amrex::MultiFab& zmom_in, + const amrex::Box& domain, + const amrex::Vector& domain_bcs_type_h); /* * Convert velocity to momentum by multiplying by density averaged onto faces @@ -47,6 +49,8 @@ void VelocityToMomentum (const amrex::MultiFab& xvel_in, amrex::MultiFab& xmom_out, amrex::MultiFab& ymom_out, amrex::MultiFab& zmom_out, + const amrex::Box& domain, + const amrex::Vector& domain_bcs_type_h, bool l_use_ndiff); /* diff --git a/Source/Utils/VelocityToMomentum.cpp b/Source/Utils/VelocityToMomentum.cpp index 81642c4c5..cc0a12b23 100644 --- a/Source/Utils/VelocityToMomentum.cpp +++ b/Source/Utils/VelocityToMomentum.cpp @@ -4,6 +4,7 @@ #include #include #include +#include using namespace amrex; @@ -19,6 +20,8 @@ using namespace amrex; * @param[out] xmom x-component of momentum * @param[out] ymom y-component of momentum * @param[out] zmom z-component of momentum + * @param[in] domain Domain at this level + * @param[in] domain_bcs_type_h host vector for domain boundary conditions * @param[in] l_use_ndiff flag describing whether we will later add explicit numerical diffusion */ @@ -30,10 +33,14 @@ void VelocityToMomentum (const MultiFab& xvel_in, const IntVect& zvel_ngrow, const MultiFab& density, MultiFab& xmom, MultiFab& ymom, MultiFab& zmom, + const Box& domain, + const Vector& domain_bcs_type_h, bool l_use_ndiff) { BL_PROFILE_VAR("VelocityToMomentum()",VelocityToMomentum); + const BCRec* bc_ptr_h = domain_bcs_type_h.data(); + // Loop over boxes = valid boxes grown by ngrow #ifdef _OPENMP #pragma omp parallel if (amrex::Gpu::notInLaunchRegion()) @@ -41,6 +48,7 @@ void VelocityToMomentum (const MultiFab& xvel_in, for ( MFIter mfi(density,TilingIfNotGPU()); mfi.isValid(); ++mfi) { // We need momentum in the interior ghost cells (init == real) + const Box& bx = mfi.tilebox(); Box tbx, tby, tbz; if (l_use_ndiff) { tbx = mfi.tilebox(IntVect(1,0,0),xvel_ngrow); @@ -67,16 +75,39 @@ void VelocityToMomentum (const MultiFab& xvel_in, ParallelFor(tbx, tby, tbz, [=] AMREX_GPU_DEVICE (int i, int j, int k) { - Real rho_x_face = 0.5 * (dens_arr(i,j,k,Rho_comp) + dens_arr(i-1,j,k,Rho_comp)); - momx(i,j,k) = velx(i,j,k) * rho_x_face; + Real rho_x_face = + 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)); }, [=] AMREX_GPU_DEVICE (int i, int j, int k) { - Real rho_y_face = 0.5 * (dens_arr(i,j,k,Rho_comp) + dens_arr(i,j-1,k,Rho_comp)); - momy(i,j,k) = vely(i,j,k) * rho_y_face; + momy(i,j,k) = vely(i,j,k) * 0.5 * (dens_arr(i,j,k,Rho_comp) + dens_arr(i,j-1,k,Rho_comp)); }, [=] AMREX_GPU_DEVICE (int i, int j, int k) { - Real rho_z_face = 0.5 * (dens_arr(i,j,k,Rho_comp) + dens_arr(i,j,k-1,Rho_comp)); - momz(i,j,k) = velz(i,j,k) * rho_z_face; + 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) { + momx(i,j,k) = velx(i,j,k) * dens_arr(i-1,j,k,Rho_comp); + }); + } + if ( (bx.bigEnd(0) == domain.bigEnd(0)) && + (bc_ptr_h[BCVars::cons_bc].hi(0) == ERFBCType::ext_dir) ) { + ParallelFor(makeSlab(tbx,0,domain.bigEnd(0)+1), [=] AMREX_GPU_DEVICE (int i, int j, int k) { + momx(i,j,k) = velx(i,j,k) * dens_arr(i,j,k,Rho_comp); + }); + } + if ( (bx.smallEnd(1) == domain.smallEnd(1)) && + (bc_ptr_h[BCVars::cons_bc].lo(1) == ERFBCType::ext_dir) ) { + ParallelFor(makeSlab(tby,1,domain.smallEnd(1)), [=] AMREX_GPU_DEVICE (int i, int j, int k) { + momy(i,j,k) = vely(i,j,k) * dens_arr(i,j-1,k,Rho_comp); + }); + } + if ( (bx.bigEnd(0) == domain.bigEnd(0)) && + (bc_ptr_h[BCVars::cons_bc].hi(0) == ERFBCType::ext_dir) ) { + ParallelFor(makeSlab(tby,1,domain.bigEnd(1)+1), [=] AMREX_GPU_DEVICE (int i, int j, int k) { + momy(i,j,k) = vely(i,j,k) * dens_arr(i,j,k,Rho_comp); + }); + } } // end MFIter }