Refluxing

Docs/sphinx_doc/Inputs.rst

@@ -860,7 +860,7 @@ List of Parameters
 | **erf.use_terrain**         | use terrain-fitted |  true / false      | false      |
 |                             | coordinates?       |                    |            |
 +-----------------------------+--------------------+--------------------+------------+
-| **erf.terrain_type**        | static or moving?  |  0 / 1             | 0          |
+| **erf.terrain_type**        | static or moving?  |  Static / Moving   | Static     |
 +-----------------------------+--------------------+--------------------+------------+
 | **erf.terrain_smoothing**   | specify terrain    | 0,                 | 0          |
 |                             | following          | 1,                 |            |

Docs/sphinx_doc/RegressionTests.rst

@@ -45,7 +45,7 @@ The following problems are currently tested in the CI:
 | DensityCurrent_detJ2_MT       | 256 4 64 | Symmetry | Periodic | SlipWall   | None  | use_terrain = true    |
 |                               |          | Outflow  |          | SlipWall   |       | uses zlevels          |
 |                               |          | Outflow  |          | SlipWall   |       | detJ = 2 everywhere   |
-|                               |          |          |          |            |       | terrain_type = 1      |
+|                               |          |          |          |            |       | terrain_type = Moving |
 +-------------------------------+----------+----------+----------+------------+-------+-----------------------+
 | EkmanSpiral                   | 4 4 400  | Periodic | Periodic | NoSlipWall | Geo   | +Coriolis             |
 |                               |          |          |          | SlipWall   |       | +gravity              |

Docs/sphinx_doc/theory/WetEquations.rst

@@ -38,7 +38,7 @@ The governing equations for this model are
   \frac{\partial \rho_d}{\partial t} &= - \nabla \cdot (\rho_d \mathbf{u})
 
   \frac{\partial (\rho_d \mathbf{u})}{\partial t} &= - \nabla \cdot (\rho_d \mathbf{u} \mathbf{u}) -
-          \frac{1}{1 + q_v + q_c}  \nabla p^\prime - \nabla \cdot \tau + \mathbf{F} + \delta_{i,3}\mathbf{B}
+          \frac{1}{1 + q_v + q_c} ( \nabla p^\prime  + \delta_{i,3}\mathbf{B} ) - \nabla \cdot \tau + \mathbf{F}
 
   \frac{\partial (\rho_d \theta_d)}{\partial t} &= - \nabla \cdot (\rho_d \mathbf{u} \theta_d)
                 + \nabla \cdot ( \rho_d \alpha_{T}\ \nabla \theta_d) + \frac{\theta_d L_v}{T_d C_{pd}} f_{cond}

Exec/DevTests/MovingTerrain/inputs

@@ -19,8 +19,8 @@ zlo.type = "SlipWall"
 zhi.type = "SlipWall"
 
 # TERRRAIN GRID TYPE
-erf.use_terrain = 1         # enable terrain stencils
-erf.terrain_type = 1        # moving terrain
+erf.use_terrain  = true     # enable terrain stencils
+erf.terrain_type = Moving   # moving terrain
 erf.terrain_smoothing = 2   # Sullivan 2004 approach
 
 # FOR NO SUBSTEPPING

Source/Advection/Advection.H

@@ -28,21 +28,24 @@ void AdvectionSrcForRhoAndTheta (const amrex::Box& bx, const amrex::Box& valid_b
                                  const amrex::Array4<const amrex::Real>& mf_u,
                                  const amrex::Array4<const amrex::Real>& mf_v,
                                  const AdvType horiz_adv_type, const AdvType vert_adv_type,
-                                 const int use_terrain);
+                                 const bool use_terrain);
 
 /** Compute advection tendency for all scalars other than density and potential temperature */
-void AdvectionSrcForScalars (const amrex::Box& bx,
+void AdvectionSrcForScalars (int level, int finest_level,
+                             const amrex::MFIter& mfi,
+                             const amrex::Box& bx,
                              const int icomp, const int ncomp,
                              const amrex::Array4<const amrex::Real>& avg_xmom,
                              const amrex::Array4<const amrex::Real>& avg_ymom,
                              const amrex::Array4<const amrex::Real>& avg_zmom,
                              const amrex::Array4<const amrex::Real>& cell_prim,
                              const amrex::Array4<amrex::Real>& src,
                              const amrex::Array4<const amrex::Real>& detJ,
+                             const amrex::Real dt,
                              const amrex::GpuArray<amrex::Real, AMREX_SPACEDIM>& cellSize,
                              const amrex::Array4<const amrex::Real>& mf_m,
                              const AdvType horiz_adv_type, const AdvType vert_adv_type,
-                             const int use_terrain);
+                             const bool use_terrain, const bool is_two_way_coupling);
 
 /** Compute advection tendencies for all components of momentum */
 void AdvectionSrcForMom (const amrex::Box& bxx, const amrex::Box& bxy, const amrex::Box& bxz,
@@ -58,7 +61,7 @@ void AdvectionSrcForMom (const amrex::Box& bxx, const amrex::Box& bxy, const amr
                          const amrex::Array4<const amrex::Real>& mf_u,
                          const amrex::Array4<const amrex::Real>& mf_v,
                          const AdvType horiz_adv_type, const AdvType vert_adv_type,
-                         const int use_terrain, const int domhi_z);
+                         const bool use_terrain, const int domhi_z);
 
 AMREX_GPU_HOST_DEVICE
 AMREX_FORCE_INLINE

Source/Advection/AdvectionSrcForMom.cpp

@@ -51,7 +51,7 @@ AdvectionSrcForMom (const Box& bxx, const Box& bxy, const Box& bxz,
                     const Array4<const Real>& mf_v,
                     const AdvType horiz_adv_type,
                     const AdvType vert_adv_type,
-                    const int use_terrain,
+                    const bool use_terrain,
                     const int domhi_z)
 {
     BL_PROFILE_VAR("AdvectionSrcForMom", AdvectionSrcForMom);

Source/Advection/AdvectionSrcForMom_N.H

@@ -42,22 +42,27 @@ AdvectionSrcForXMom_N (int i, int j, int k,
     amrex::Real interp_hi(0.), interp_lo(0.);
 
     rho_u_avg_hi = 0.5 * (rho_u(i+1, j, k) * mf_u_inv(i+1,j,0) + rho_u(i, j, k) * mf_u_inv(i,j,0));
-    rho_u_avg_lo = 0.5 * (rho_u(i-1, j, k) * mf_u_inv(i-1,j,0) + rho_u(i, j, k) * mf_u_inv(i,j,0));
-    interp_u_h.InterpolateInX(i,j,k,0,interp_hi,interp_lo,rho_u_avg_hi,rho_u_avg_lo);
+    interp_u_h.InterpolateInX(i+1,j,k,0,interp_hi,rho_u_avg_hi);
     xflux_hi = rho_u_avg_hi * interp_hi;
+
+    rho_u_avg_lo = 0.5 * (rho_u(i-1, j, k) * mf_u_inv(i-1,j,0) + rho_u(i, j, k) * mf_u_inv(i,j,0));
+    interp_u_h.InterpolateInX(i,j,k,0,interp_lo,rho_u_avg_lo);
     xflux_lo = rho_u_avg_lo * interp_lo;
 
     rho_v_avg_hi = 0.5 * (rho_v(i, j+1, k) * mf_v_inv(i,j+1,0) + rho_v(i-1, j+1, k) * mf_v_inv(i-1,j+1,0));
-    rho_v_avg_lo = 0.5 * (rho_v(i, j  , k) * mf_v_inv(i,j  ,0) + rho_v(i-1, j  , k) * mf_v_inv(i-1,j  ,0));
-    interp_u_h.InterpolateInY(i,j,k,0,interp_hi,interp_lo,rho_v_avg_hi,rho_v_avg_lo);
+    interp_u_h.InterpolateInY(i,j+1,k,0,interp_hi,rho_v_avg_hi);
     yflux_hi = rho_v_avg_hi * interp_hi;
+
+    rho_v_avg_lo = 0.5 * (rho_v(i, j  , k) * mf_v_inv(i,j  ,0) + rho_v(i-1, j  , k) * mf_v_inv(i-1,j  ,0));
+    interp_u_h.InterpolateInY(i,j,k,0,interp_lo,rho_v_avg_lo);
     yflux_lo = rho_v_avg_lo * interp_lo;
 
     rho_w_avg_hi = 0.5 * (rho_w(i, j, k+1) + rho_w(i-1, j, k+1));
-    rho_w_avg_lo = 0.5 * (rho_w(i, j, k  ) + rho_w(i-1, j, k  ));
-    interp_u_v.InterpolateInZ_hi(i,j,k,0,interp_hi,rho_w_avg_hi);
-    interp_u_v.InterpolateInZ_lo(i,j,k,0,interp_lo,rho_w_avg_lo);
+    interp_u_v.InterpolateInZ(i,j,k+1,0,interp_hi,rho_w_avg_hi);
     zflux_hi = rho_w_avg_hi * interp_hi;
+
+    rho_w_avg_lo = 0.5 * (rho_w(i, j, k  ) + rho_w(i-1, j, k  ));
+    interp_u_v.InterpolateInZ(i,j,k,0,interp_lo,rho_w_avg_lo);
     zflux_lo = rho_w_avg_lo * interp_lo;
 
     amrex::Real mfsq = 1 / (mf_u_inv(i,j,0) * mf_u_inv(i,j,0));
@@ -110,22 +115,27 @@ AdvectionSrcForYMom_N (int i, int j, int k,
     amrex::Real interp_hi(0.), interp_lo(0.);
 
     rho_u_avg_hi = 0.5 * (rho_u(i+1, j, k) * mf_u_inv(i+1,j  ,0) + rho_u(i+1, j-1, k) * mf_u_inv(i+1,j-1,0));
-    rho_u_avg_lo = 0.5 * (rho_u(i  , j, k) * mf_u_inv(i  ,j  ,0) + rho_u(i  , j-1, k) * mf_u_inv(i  ,j-1,0));
-    interp_v_h.InterpolateInX(i,j,k,0,interp_hi,interp_lo,rho_u_avg_hi,rho_u_avg_lo);
+    interp_v_h.InterpolateInX(i+1,j,k,0,interp_hi,rho_u_avg_hi);
     xflux_hi = rho_u_avg_hi * interp_hi;
+
+    rho_u_avg_lo = 0.5 * (rho_u(i  , j, k) * mf_u_inv(i  ,j  ,0) + rho_u(i  , j-1, k) * mf_u_inv(i  ,j-1,0));
+    interp_v_h.InterpolateInX(i,j,k,0,interp_lo,rho_u_avg_lo);
     xflux_lo = rho_u_avg_lo * interp_lo;
 
     rho_v_avg_hi = 0.5 * (rho_v(i, j, k) * mf_v_inv(i  ,j  ,0) + rho_v(i, j+1, k) * mf_v_inv(i  ,j+1,0));
-    rho_v_avg_lo = 0.5 * (rho_v(i, j, k) * mf_v_inv(i  ,j  ,0) + rho_v(i, j-1, k) * mf_v_inv(i  ,j-1,0));
-    interp_v_h.InterpolateInY(i,j,k,0,interp_hi,interp_lo,rho_v_avg_hi,rho_v_avg_lo);
+    interp_v_h.InterpolateInY(i,j+1,k,0,interp_hi,rho_v_avg_hi);
     yflux_hi = rho_v_avg_hi * interp_hi;
+
+    rho_v_avg_lo = 0.5 * (rho_v(i, j, k) * mf_v_inv(i  ,j  ,0) + rho_v(i, j-1, k) * mf_v_inv(i  ,j-1,0));
+    interp_v_h.InterpolateInY(i,j,k,0,interp_lo,rho_v_avg_lo);
     yflux_lo = rho_v_avg_lo * interp_lo;
 
     rho_w_avg_hi = 0.5 * (rho_w(i, j, k+1) + rho_w(i, j-1, k+1));
-    rho_w_avg_lo = 0.5 * (rho_w(i, j, k  ) + rho_w(i, j-1, k  ));
-    interp_v_v.InterpolateInZ_hi(i,j,k,0,interp_hi,rho_w_avg_hi);
-    interp_v_v.InterpolateInZ_lo(i,j,k,0,interp_lo,rho_w_avg_lo);
+    interp_v_v.InterpolateInZ(i,j,k+1,0,interp_hi,rho_w_avg_hi);
     zflux_hi = rho_w_avg_hi * interp_hi;
+
+    rho_w_avg_lo = 0.5 * (rho_w(i, j, k  ) + rho_w(i, j-1, k  ));
+    interp_v_v.InterpolateInZ(i,j,k  ,0,interp_lo,rho_w_avg_lo);
     zflux_lo = rho_w_avg_lo * interp_lo;
 
     amrex::Real mfsq = 1 / (mf_v_inv(i,j,0) * mf_v_inv(i,j,0));
@@ -186,15 +196,19 @@ AdvectionSrcForZMom_N (int i, int j, int k,
     amrex::Real interp_hi(0.), interp_lo(0.);
 
     rho_u_avg_hi = 0.5 * (rho_u(i+1, j, k) + rho_u(i+1, j, k-1)) * mf_u_inv(i+1,j  ,0);
-    rho_u_avg_lo = 0.5 * (rho_u(i  , j, k) + rho_u(i  , j, k-1)) * mf_u_inv(i  ,j  ,0);
-    interp_w_h.InterpolateInX(i,j,k,0,interp_hi,interp_lo,rho_u_avg_hi,rho_u_avg_lo);
+    interp_w_h.InterpolateInX(i+1,j,k,0,interp_hi,rho_u_avg_hi);
     xflux_hi = rho_u_avg_hi * interp_hi;
+
+    rho_u_avg_lo = 0.5 * (rho_u(i  , j, k) + rho_u(i  , j, k-1)) * mf_u_inv(i  ,j  ,0);
+    interp_w_h.InterpolateInX(i,j,k,0,interp_lo,rho_u_avg_lo);
     xflux_lo = rho_u_avg_lo * interp_lo;
 
     rho_v_avg_hi = 0.5 * (rho_v(i, j+1, k) + rho_v(i, j+1, k-1)) * mf_v_inv(i  ,j+1,0);
-    rho_v_avg_lo = 0.5 * (rho_v(i, j  , k) + rho_v(i, j  , k-1)) * mf_v_inv(i  ,j  ,0);
-    interp_w_h.InterpolateInY(i,j,k,0,interp_hi,interp_lo,rho_v_avg_hi,rho_v_avg_lo);
+    interp_w_h.InterpolateInY(i,j+1,k,0,interp_hi,rho_v_avg_hi);
     yflux_hi = rho_v_avg_hi * interp_hi;
+
+    rho_v_avg_lo = 0.5 * (rho_v(i, j  , k) + rho_v(i, j  , k-1)) * mf_v_inv(i  ,j  ,0);
+    interp_w_h.InterpolateInY(i,j,k,0,interp_lo,rho_v_avg_lo);
     yflux_lo = rho_v_avg_lo * interp_lo;
 
     // int l_spatial_order_hi = std::min(std::min(vert_spatial_order, 2*(domhi_z+1-k)), 2*(k+1));
@@ -208,15 +222,15 @@ AdvectionSrcForZMom_N (int i, int j, int k,
         rho_w_avg_hi = 0.5 * (rho_w(i,j,k) + rho_w(i,j,k+1));
         if (k == domhi_z)
         {
-            interp_w_wall.InterpolateInZ_hi(i,j,k,0,interp_hi,rho_w_avg_hi,AdvType::Centered_2nd);
+            interp_w_wall.InterpolateInZ(i,j,k+1,0,interp_hi,rho_w_avg_hi,AdvType::Centered_2nd);
         } else if (k == domhi_z-1 || k == 1) {
             if (vert_adv_type != AdvType::Centered_2nd && vert_adv_type != AdvType::Upwind_3rd) {
-               interp_w_wall.InterpolateInZ_hi(i,j,k,0,interp_hi,rho_w_avg_hi,AdvType::Centered_4th);
+               interp_w_wall.InterpolateInZ(i,j,k+1,0,interp_hi,rho_w_avg_hi,AdvType::Centered_4th);
             } else {
-               interp_w_wall.InterpolateInZ_hi(i,j,k,0,interp_hi,rho_w_avg_hi,vert_adv_type);
+               interp_w_wall.InterpolateInZ(i,j,k+1,0,interp_hi,rho_w_avg_hi,vert_adv_type);
             }
         } else {
-            interp_w_v.InterpolateInZ_hi(i,j,k,0,interp_hi,rho_w_avg_hi);
+            interp_w_v.InterpolateInZ(i,j,k+1,0,interp_hi,rho_w_avg_hi);
         }
         zflux_hi = rho_w_avg_hi * interp_hi;
     }
@@ -231,15 +245,15 @@ AdvectionSrcForZMom_N (int i, int j, int k,
     } else {
         rho_w_avg_lo = 0.5 * (rho_w(i,j,k) + rho_w(i,j,k-1));
         if (k == 1) {
-            interp_w_wall.InterpolateInZ_lo(i,j,k,0,interp_lo,rho_w_avg_lo,AdvType::Centered_2nd);
+            interp_w_wall.InterpolateInZ(i,j,k,0,interp_lo,rho_w_avg_lo,AdvType::Centered_2nd);
         } else if (k == 2 || k == domhi_z) {
             if (vert_adv_type != AdvType::Centered_2nd && vert_adv_type != AdvType::Upwind_3rd) {
-                interp_w_wall.InterpolateInZ_lo(i,j,k,0,interp_lo,rho_w_avg_lo,AdvType::Centered_4th);
+                interp_w_wall.InterpolateInZ(i,j,k,0,interp_lo,rho_w_avg_lo,AdvType::Centered_4th);
             } else {
-                interp_w_wall.InterpolateInZ_lo(i,j,k,0,interp_lo,rho_w_avg_lo,vert_adv_type);
+                interp_w_wall.InterpolateInZ(i,j,k,0,interp_lo,rho_w_avg_lo,vert_adv_type);
             }
         } else {
-            interp_w_v.InterpolateInZ_lo(i,j,k,0,interp_lo,rho_w_avg_lo);
+            interp_w_v.InterpolateInZ(i,j,k,0,interp_lo,rho_w_avg_lo);
         }
         zflux_lo = rho_w_avg_lo * interp_lo;
     }