Merge Couette and Poiseuille and do some other cleanup (#1542)

CMake/BuildERFExe.cmake

@@ -49,6 +49,18 @@ function(build_erf_lib erf_lib_name)
     target_compile_definitions(${erf_lib_name} PUBLIC ERF_USE_PARTICLES)
   endif()
 
+  if(ERF_ENABLE_EB)
+    target_sources(${erf_lib_name} PRIVATE
+                   ${SRC_DIR}/EB/Init_EB.cpp
+                   ${SRC_DIR}/EB/eb_box.cpp
+                   ${SRC_DIR}/EB/eb_cylinder.cpp
+                   ${SRC_DIR}/EB/eb_regular.cpp
+                   ${SRC_DIR}/EB/initEB.cpp
+                   ${SRC_DIR}/EB/writeEBsurface.cpp)
+    target_include_directories(${erf_lib_name} PUBLIC ${SRC_DIR}/EB)
+    target_compile_definitions(${erf_lib_name} PUBLIC ERF_USE_EB)
+  endif()
+
   if(ERF_ENABLE_NETCDF)
     target_sources(${erf_lib_name} PRIVATE
                    ${SRC_DIR}/IO/NCInterface.cpp

CMake/SetAmrexOptions.cmake

@@ -37,6 +37,11 @@ if(ERF_ENABLE_PARTICLES)
   set(AMReX_PARTICLES ON)
 endif()
 
+set(AMReX_EB OFF)
+if(ERF_ENABLE_EB)
+  set(AMReX_EB ON)
+endif()
+
 if(ERF_ENABLE_CUDA)
   set(AMReX_GPU_BACKEND CUDA CACHE STRING "AMReX GPU type" FORCE)
   set(AMReX_CUDA_WARN_CAPTURE_THIS OFF)

Docs/sphinx_doc/RegressionTests.rst

@@ -408,7 +408,7 @@ Test Location: `Tests/test_files/CouetteFlow`_
 
 Problem Location: `Exec/RegTests/CouetteFlow`_
 
-.. _`Exec/RegTests/CouetteFlow`: https://github.com/erf-model/ERF/tree/development/Exec/RegTests/CouetteFlow
+.. _`Exec/RegTests/Couette_Poiseuille`: https://github.com/erf-model/ERF/tree/development/Exec/RegTests/Couette_Poiseuille
 
 Poiseuille Flow
 ---------------
@@ -424,7 +424,7 @@ Test Location: `Tests/test_files/PoiseuilleFlow`_
 
 Problem Location: `Exec/RegTests/PoiseuilleFlow`_
 
-.. _`Exec/RegTests/PoiseuilleFlow`: https://github.com/erf-model/ERF/tree/development/Exec/RegTests/PoiseuilleFlow
+.. _`Exec/RegTests/Couette_Poiseuille`: https://github.com/erf-model/ERF/tree/development/Exec/RegTests/Couette_Poiseuille
 
 Nonlinear Density Current
 ---------------------------

Exec/CMakeLists.txt

@@ -6,12 +6,11 @@ if (ERF_ENABLE_MULTIBLOCK)
   add_subdirectory(DevTests/MultiBlock)
 elseif (ERF_ENABLE_REGRESSION_TESTS_ONLY)
 #  add_subdirectory(RegTests/Bubble)
-  add_subdirectory(RegTests/CouetteFlow)
+  add_subdirectory(RegTests/Couette_Poiseuille)
   add_subdirectory(RegTests/DensityCurrent)
   add_subdirectory(RegTests/DynamicRefinement)
   add_subdirectory(RegTests/EkmanSpiral_custom)
   add_subdirectory(RegTests/IsentropicVortex)
-  add_subdirectory(RegTests/PoiseuilleFlow)
   add_subdirectory(RegTests/ScalarAdvDiff)
   add_subdirectory(RegTests/TaylorGreenVortex)
   add_subdirectory(DevTests/MovingTerrain)
@@ -21,15 +20,14 @@ else ()
   add_subdirectory(Radiation)
   add_subdirectory(SquallLine_2D)
   add_subdirectory(RegTests/Bubble)
-  add_subdirectory(RegTests/CouetteFlow)
+  add_subdirectory(RegTests/Couette_Poiseuille)
   add_subdirectory(RegTests/DensityCurrent)
   add_subdirectory(RegTests/DynamicRefinement)
   add_subdirectory(RegTests/EkmanSpiral_custom)
   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)

Exec/RegTests/CouetteFlow/CMakeLists.txt → Exec/RegTests/Couette_Poiseuille/CMakeLists.txt

@@ -1,4 +1,4 @@
-set(erf_exe_name erf_couette_flow)
+set(erf_exe_name erf_couette_poiseuille)
 
 add_executable(${erf_exe_name} "")
 target_sources(${erf_exe_name}

Exec/RegTests/CouetteFlow/GNUmakefile → Exec/RegTests/Couette_Poiseuille/GNUmakefile

@@ -29,5 +29,5 @@ Bpack := ./Make.package
 Blocs := .
 
 ERF_HOME := ../../..
-ERF_PROBLEM_DIR = $(ERF_HOME)/Exec/RegTests/CouetteFlow
+ERF_PROBLEM_DIR = $(ERF_HOME)/Exec/RegTests/Couette_Poiseuille
 include $(ERF_HOME)/Exec/Make.ERF

Exec/RegTests/Couette_Poiseuille/README

@@ -0,0 +1,13 @@
+This directory has two problem setups: Couette flow and Poiseuille flow.
+These problems primarily test the implementation of the viscous terms,
+in particular the discretization at the boundaries.
+
+If prob_type == 1 in the inputs file, then the problem setup is for classic Couette flow.
+The domain is periodic in the x- and y-directions and has no slip walls on the top and bottom.
+The velocity is specified to be (2,0,0) (for inputs_couette_x) or
+(0,2,0) (for inputs_couette_y) at the top wall.
+
+If prob_type == 10 or 11 in the inputs file, then the problem setup is for classic Poiseuille flow.
+The domain is periodic in the x- and y-directions and has no slip walls on the top and bottom.
+The flow is driven by a pressure gradient in the x-direction if prob_type == 10, 
+or in the y-direction if prob_type = 11.

Exec/RegTests/CouetteFlow/inputs_couette_x → Exec/RegTests/Couette_Poiseuille/inputs_couette_x

@@ -1,4 +1,4 @@
-# ------------------  INPUTS TO MAIN PROGRAM  -------------------
+#------------------  INPUTS TO MAIN PROGRAM  -------------------
 max_step = 50
 
 amrex.fpe_trap_invalid = 1
@@ -49,6 +49,8 @@ erf.dynamicViscosity = 0.1
 erf.init_type = "uniform"
 
 # PROBLEM PARAMETERS
+prob.prob_type = 1
+
 prob.rho_0 = 1.0
 prob.T_0 = 300.0
 # NOTE: this u_0 should match the zhi.velocity specified above

Exec/RegTests/CouetteFlow/inputs_couette_y → Exec/RegTests/Couette_Poiseuille/inputs_couette_y

@@ -49,6 +49,8 @@ erf.dynamicViscosity = 0.1
 erf.init_type = "uniform"
 
 # PROBLEM PARAMETERS
+prob.prob_type = 1
+
 prob.rho_0 = 1.0
 prob.T_0 = 300.0
 # NOTE: this v_0 should match the zhi.velocity specified above

Exec/RegTests/PoiseuilleFlow/inputs_neumann_test → Exec/RegTests/Couette_Poiseuille/inputs_neumann_test

@@ -66,5 +66,7 @@ erf.use_coriolis = false
 erf.init_type = "uniform"
 
 # PROBLEM PARAMETERS
+prob.prob_type = 10
+
 prob.rho_0 = 2.0
 prob.T_0   = 300.0

Exec/RegTests/CouetteFlow/prob.H → Exec/RegTests/Couette_Poiseuille/prob.H

@@ -9,10 +9,12 @@
 
 struct ProbParm : ProbParmDefaults {
   amrex::Real rho_0 = 1.0;
-  amrex::Real T_0   = 300.0;
   amrex::Real u_0   = 0.0;
   amrex::Real v_0   = 0.0;
   amrex::Real w_0   = 0.0;
+  amrex::Real T_0   = 300.0;
+
+  int prob_type = -1;
 }; // namespace ProbParm
 
 class Problem : public ProblemBase
@@ -43,7 +45,15 @@ public:
         const SolverChoice& sc) override;
 
 protected:
-    std::string name() override { return "Couette Flow"; }
+    std::string name() override {
+        if ( parms.prob_type == 1 ) {
+            return "Couette Flow";
+        } else if ( parms.prob_type == 10 || parms.prob_type == 11) {
+            return "Poiseuille Flow";
+        } else {
+            return "Unknown";
+        }
+    }
 
 private:
     ProbParm parms;

Exec/RegTests/Couette_Poiseuille/prob.cpp

@@ -0,0 +1,120 @@
+#include "prob.H"
+
+using namespace amrex;
+
+std::unique_ptr<ProblemBase>
+amrex_probinit(
+    const amrex_real* /*problo*/,
+    const amrex_real* /*probhi*/)
+{
+    return std::make_unique<Problem>();
+}
+
+Problem::Problem()
+{
+  // Parse params
+  amrex::ParmParse pp("prob");
+  pp.query("rho_0", parms.rho_0);
+  pp.query("u_0", parms.u_0);
+  pp.query("v_0", parms.v_0);
+  pp.query("w_0", parms.w_0);
+  pp.query("T_0", parms.T_0);
+
+  pp.get("prob_type", parms.prob_type);
+
+  init_base_parms(parms.rho_0, parms.T_0);
+}
+
+void
+Problem::init_custom_pert(
+    const Box& bx,
+    const Box& xbx,
+    const Box& ybx,
+    const Box& zbx,
+    Array4<Real const> const& /*state*/,
+    Array4<Real      > const& state_pert,
+    Array4<Real      > const& x_vel_pert,
+    Array4<Real      > const& y_vel_pert,
+    Array4<Real      > const& z_vel_pert,
+    Array4<Real      > const& /*r_hse*/,
+    Array4<Real      > const& /*p_hse*/,
+    Array4<Real const> const& /*z_nd*/,
+    Array4<Real const> const& /*z_cc*/,
+    GeometryData const& geomdata,
+    Array4<Real const> const& /*mf_m*/,
+    Array4<Real const> const& /*mf_u*/,
+    Array4<Real const> const& /*mf_v*/,
+    const SolverChoice& sc)
+{
+    const bool use_moisture = (sc.moisture_type != MoistureType::None);
+
+    ParallelFor(bx, [=, parms=parms] AMREX_GPU_DEVICE(int i, int j, int k) noexcept
+    {
+        state_pert(i, j, k, RhoScalar_comp) = 0.0;
+
+        if (use_moisture) {
+            state_pert(i, j, k, RhoQ1_comp) = 0.0;
+            state_pert(i, j, k, RhoQ2_comp) = 0.0;
+        }
+    });
+
+    AMREX_ALWAYS_ASSERT (parms.prob_type == 1 || parms.prob_type == 10 || parms.prob_type == 11);
+
+    // Couette flow
+    if (parms.prob_type == 1) {
+
+        ParallelFor(xbx, [=, parms=parms] AMREX_GPU_DEVICE(int i, int j, int k) noexcept {
+            const auto *const prob_hi  = geomdata.ProbHi();
+            const auto *const dx       = geomdata.CellSize();
+            const Real z = (k + 0.5) * dx[2];
+            x_vel_pert(i, j, k) = parms.u_0 * z / prob_hi[2];
+        });
+
+        ParallelFor(ybx, [=, parms=parms] AMREX_GPU_DEVICE(int i, int j, int k) noexcept {
+            const auto *const prob_hi  = geomdata.ProbHi();
+            const auto *const dx       = geomdata.CellSize();
+            const Real z = (k + 0.5) * dx[2];
+            y_vel_pert(i, j, k) = parms.v_0 * z / prob_hi[2];
+        });
+
+        ParallelFor(zbx, [=, parms=parms] AMREX_GPU_DEVICE(int i, int j, int k) noexcept {
+            z_vel_pert(i, j, k) = parms.w_0;
+        });
+
+    // Poiseuille flow
+    } else if (parms.prob_type == 10 || parms.prob_type == 11) {
+
+        ParallelFor(xbx, [=, parms=parms] AMREX_GPU_DEVICE(int i, int j, int k) noexcept
+        {
+            const Real* prob_lo = geomdata.ProbLo();
+            const Real* dx      = geomdata.CellSize();
+            const Real z_h = prob_lo[2] + (k + 0.5) *  dx[2];
+
+            // Set the x-velocity to be a parabolic profile with max 1 at z = 0 and 0 at z = +/-1
+            if (parms.prob_type == 10) {
+                x_vel_pert(i, j, k) = 1.0 - z_h * z_h;
+            } else {
+                x_vel_pert(i, j, k) = 0.0;
+            }
+        });
+
+        ParallelFor(ybx, [=, parms=parms] AMREX_GPU_DEVICE(int i, int j, int k) noexcept
+        {
+            const Real* prob_lo = geomdata.ProbLo();
+            const Real* dx      = geomdata.CellSize();
+            const Real z_h = prob_lo[2] + (k + 0.5) *  dx[2];
+
+            // Set the x-velocity to be a parabolic profile with max 1 at z = 0 and 0 at z = +/-1
+            if (parms.prob_type == 11) {
+               y_vel_pert(i, j, k) = 1.0 - z_h * z_h;
+            } else {
+               y_vel_pert(i, j, k) = 0.0;
+            }
+        });
+
+        ParallelFor(zbx, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept
+        {
+            z_vel_pert(i, j, k) = 0.0;
+        });
+    } // prob_type
+}