enforce inflow value of density when converting between velocity and … (

#1509)

* enforce inflow value of density when converting between velocity and momentum

* move ParticlesOverWoA and minor cleanup

* whitespace

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)

Exec/DevTests/ParticlesOverWoA/GNUmakefile → 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

Exec/DevTests/ParticlesOverWoA/inputs → 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

Exec/DevTests/ParticlesOverWoA/prob.cpp → Exec/RegTests/ParticlesOverWoA/prob.cpp

@@ -48,6 +48,7 @@ Problem::init_custom_pert(
     Array4<Real const> 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

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());

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);
     }
 }
 

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);

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<BCRec>& 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
 }

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<amrex::BCRec>& 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<amrex::BCRec>& domain_bcs_type_h,
                          bool l_use_ndiff);
 
 /*

Source/Utils/VelocityToMomentum.cpp

@@ -4,6 +4,7 @@
 #include <AMReX.H>
 #include <AMReX_MultiFab.H>
 #include <Utils.H>
+#include <AMReX_MultiFabUtil.H>
 
 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,17 +33,22 @@ void VelocityToMomentum (const MultiFab& xvel_in,
                          const  IntVect& zvel_ngrow,
                          const MultiFab& density,
                          MultiFab& xmom, MultiFab& ymom, MultiFab& zmom,
+                         const Box& domain,
+                         const Vector<BCRec>& 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())
 #endif
     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
 }