Boundary condition treatment for wallmodel and symmetry wall; fixes f…

…vm problems (#844)



Co-authored-by: prakash <prakash.mohan@nrel.gov>

amr-wind/boundary_conditions/BCInterface.cpp

@@ -76,6 +76,8 @@ void BCIface::read_bctype()
             ibctype[ori] = BC::fixed_gradient;
         } else if ((bcstr == "wave_generation") || (bcstr == "wg")) {
             ibctype[ori] = BC::wave_generation;
+        } else if ((bcstr == "symmetric_wall") || (bcstr == "symw")) {
+            ibctype[ori] = BC::symmetric_wall;
         } else {
             ibctype[ori] = BC::undefined;
         }
@@ -206,6 +208,19 @@ void BCVelocity::set_bcrec()
                 bcrec[dir].setHi(dir, amrex::BCType::ext_dir);
             }
             break;
+        case BC::symmetric_wall:
+            if (side == amrex::Orientation::low) {
+                // Tangential directions use first-order extrapolation
+                set_bcrec_lo(dir, amrex::BCType::foextrap);
+                // Normal direction uses dirichlet (ext-dir)
+                bcrec[dir].setLo(dir, amrex::BCType::ext_dir);
+            } else {
+                // Tangential directions use first-order extrapolation
+                set_bcrec_hi(dir, amrex::BCType::foextrap);
+                // Normal direction uses dirichlet (ext-dir)
+                bcrec[dir].setHi(dir, amrex::BCType::ext_dir);
+            }
+            break;
 
         default:
             amrex::Abort("Invalid incflo BC type encountered");
@@ -254,6 +269,7 @@ void BCScalar::set_bcrec()
         case BC::pressure_inflow:
         case BC::pressure_outflow:
         case BC::zero_gradient:
+        case BC::symmetric_wall:
             if (side == amrex::Orientation::low) {
                 set_bcrec_lo(dir, amrex::BCType::foextrap);
             } else {

amr-wind/boundary_conditions/README.org

@@ -42,14 +42,17 @@ The boundary conditions that can be specified in the input file are:
   if the user does not specify any inputs in the input file.
 
 - =zero_gradient= :: Sets =BCRec= to =foextrap= and =LinOpBCType= to Neumann.
-  Like =no_slip_wall=, the normal component of velocity is automatically reset
-  to zero. The value on the boundary face is the same as the value at the cell
+  Note, unlike the =no_slip_wall=, the normal component of velocity is also zero
+  Neumann. The value on the boundary face is the same as the value at the cell
   center. No additional inputs are necessary when this BC is specified.
 
 - =slip_wall= :: =BCRec = hoextrap= and =LinOpBCType = Neumann=. The value on
   the boundary face is /extrapolated/ from the two interior cells using a
   one-sided stencil. No additional inputs are necessary when this BC is specified.
 
+- =symmetric_wall= ::  Sets =BCRec= to =foextrap= and =LinOpBCType= to Neumann.
+  Similar to =zero_gradient= but the normal component of velocity is set to 0.
+
 - =fixed_gradient= :: =BCRec = hoextrap= and =LinOpBCType = inhomogNeumann=.
   Like =slip_wall= the value on the face is extrapolated (during fillpatch)
   using a one-sided stencil. However, for linear solvers the value on the face
@@ -86,7 +89,8 @@ possible for scalars and tangential component of velocity.
 | pressure_outflow | po    | foextrap    | foextrap | foextrap | foextrap |
 | mass_inflow      | mi    | ext_dir     | ext_dir  | ext_dir  | foextrap |
 | no_slip_wall     | nsw   | ext_dir (0) | ext_dir  | ext_dir  | foextrap |
-| zero_gradient    | zg    | ext_dir (0) | foextrap | foextrap | foextrap |
+| symmetric_wall   | symw  | ext_dir (0) | foextrap | foextrap | foextrap |
+| zero_gradient    | zg    | foextrap    | foextrap | foextrap | foextrap |
 | slip_wall        | sw    | ext_dir (0) | hoextrap | hoextrap | foextrap |
 | wall_model       | wm    | ext_dir (0) | hoextrap | hoextrap | foextrap |
 | fixed_gradient   | fg    | hoextrap    | hoextrap | hoextrap | foextrap |
@@ -103,6 +107,7 @@ projection (pressure) and MAC projections.
 | pressure_outflow | po    | Dirichlet        |
 | mass_inflow      | mi    | Neumann          |
 | zero_gradient    | zg    | Neumann          |
+| symmetric_wall   | symw  | Neumann          |
 | no_slip_wall     | nsw   | Neumann          |
 | slip_wall        | sw    | Neumann          |
 | wall_model       | wm    | Neumann          |
@@ -118,6 +123,7 @@ projection (pressure) and MAC projections.
 | mass_inflow      | mi    | Dirichlet        | Dirichlet            |
 | zero_gradient    | zg    | Dirichlet (0)    | Neumann              |
 | no_slip_wall     | nsw   | Dirichlet (0)    | Dirichlet            |
+| symmetric_wall   | symw  | Dirichlet (0)    | Neumann              |
 | slip_wall        | sw    | Dirichlet (0)    | Neumann              |
 | wall_model       | wm    | Dirichlet (0)    | Inhomog Neumann      |
 | fixed_gradient   | fg    | Inhomog Neumann  | Inhomog Neumann      |
@@ -132,6 +138,7 @@ projection (pressure) and MAC projections.
 | pressure_outflow | po    | Neumann          |
 | mass_inflow      | mi    | Dirichlet        |
 | zero_gradient    | zg    | Neumann          |
+| symmetric_wall   | symw  | Neumann          |
 | no_slip_wall     | nsw   | Dirichlet        |
 | slip_wall        | sw    | Neumann          |
 | wall_model       | wm    | Inhomog Neumann  |

amr-wind/diffusion/incflo_diffusion.cpp

@@ -36,6 +36,7 @@ Vector<Array<LinOpBCType, AMREX_SPACEDIM>> get_diffuse_tensor_bc(
                 r[2][dir] = LinOpBCType::Dirichlet;
                 break;
             }
+            case BC::symmetric_wall:
             case BC::slip_wall: {
                 // Tangential components are Neumann
                 // Normal     component  is  Dirichlet
@@ -77,6 +78,7 @@ get_diffuse_scalar_bc(amr_wind::Field& scalar, Orientation::Side side) noexcept
             case BC::pressure_inflow:
             case BC::pressure_outflow:
             case BC::zero_gradient:
+            case BC::symmetric_wall:
             case BC::slip_wall: {
                 r[dir] = LinOpBCType::Neumann;
                 break;

amr-wind/equation_systems/vof/split_advection.H

@@ -124,15 +124,15 @@ AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE void fluxes_bc_save(
     if (dir == 0) {
         // For wall BCs, do not allow flow into or out of domain
         if (bclo == BC::no_slip_wall || bclo == BC::slip_wall ||
-            bclo == BC::wall_model) {
+            bclo == BC::wall_model || bclo == BC::symmetric_wall) {
             if (i == domlo) {
                 vofR(domlo - 1, j, k) = 0.0;
                 vofL(domlo, j, k) = 0.0;
             }
         }
 
         if (bchi == BC::no_slip_wall || bchi == BC::slip_wall ||
-            bchi == BC::wall_model) {
+            bchi == BC::wall_model || bchi == BC::symmetric_wall) {
             if (i == domhi + 1) {
                 vofL(domhi + 1, j, k) = 0.0;
                 vofR(domhi, j, k) = 0.0;
@@ -146,14 +146,14 @@ AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE void fluxes_bc_save(
         // eulerian_implicit(), so this sum is not for the sake of accumulation
     } else if (dir == 1) {
         if (bclo == BC::no_slip_wall || bclo == BC::slip_wall ||
-            bclo == BC::wall_model) {
+            bclo == BC::wall_model || bclo == BC::symmetric_wall) {
             if (j == domlo) {
                 vofR(i, domlo - 1, k) = 0.0;
                 vofL(i, domlo, k) = 0.0;
             }
         }
         if (bchi == BC::no_slip_wall || bchi == BC::slip_wall ||
-            bchi == BC::wall_model) {
+            bchi == BC::wall_model || bchi == BC::symmetric_wall) {
             if (j == domhi + 1) {
                 vofL(i, domhi + 1, k) = 0.0;
                 vofR(i, domhi, k) = 0.0;
@@ -163,14 +163,14 @@ AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE void fluxes_bc_save(
         f_f(i, j, k) = advalpha_f(i, j, k) * disp;
     } else if (dir == 2) {
         if (bclo == BC::no_slip_wall || bclo == BC::slip_wall ||
-            bclo == BC::wall_model) {
+            bclo == BC::wall_model || bclo == BC::symmetric_wall) {
             if (k == domlo) {
                 vofR(i, j, domlo - 1) = 0.0;
                 vofL(i, j, domlo) = 0.0;
             }
         }
         if (bchi == BC::no_slip_wall || bchi == BC::slip_wall ||
-            bchi == BC::wall_model) {
+            bchi == BC::wall_model || bchi == BC::symmetric_wall) {
             if (k == domhi + 1) {
                 vofL(i, j, domhi + 1) = 0.0;
                 vofR(i, j, domhi) = 0.0;

amr-wind/equation_systems/vof/vof_bciface.H

@@ -41,6 +41,7 @@ protected:
             case BC::pressure_outflow:
             case BC::zero_gradient:
             case BC::no_slip_wall:
+            case BC::symmetric_wall:
                 if (side == amrex::Orientation::low) {
                     set_bcrec_lo(dir, amrex::BCType::foextrap);
                 } else {

amr-wind/incflo_enums.H

@@ -14,6 +14,7 @@ enum struct BC {
     fixed_gradient,
     periodic,
     wave_generation,
+    symmetric_wall,
     undefined
 };
 

amr-wind/physics/ChannelFlow.H

@@ -81,6 +81,8 @@ private:
 
     amrex::Real m_mean_vel;
 
+    bool m_half{false};
+
     //! direction of mean velocity
     int m_mean_vel_dir{0};
 

amr-wind/physics/ChannelFlow.cpp

@@ -31,6 +31,7 @@ ChannelFlow::ChannelFlow(CFDSim& sim)
         if (m_laminar) {
             pp.query("density", m_rho);
             pp.query("Mean_Velocity", m_mean_vel);
+            pp.query("half_channel", m_half);
         } else {
             pp.query("density", m_rho);
             pp.query("re_tau", m_re_tau);
@@ -159,6 +160,10 @@ amrex::Real ChannelFlow::compute_error(const IndexSelector& idxOp)
             }
         }
         ht = probhi_physical[norm_dir] - problo[norm_dir];
+        // For half channel, channel height is double the domain height
+        if (m_half) {
+            ht *= 2.0;
+        }
     }
 
     const auto& velocity = m_repo.get_field("velocity");

docs/sphinx/user/inputs_Boundary_conditions.rst

@@ -9,7 +9,7 @@ This section controls the boundary conditions. Only non-periodic BC's need to be
 
    Boundary condition type on the lo (or hi) side of the domain. 
    Current options are: periodic, pressure_inflow, pressure_outflow, mass_inflow, 
-   no_slip_wall, slip_wall, and wall_model. 
+   no_slip_wall, slip_wall, symmetric_wall and wall_model.
 
 .. input_param:: xlo.temperature (or ylo.temperature, etc)
 

test/CMakeLists.txt

@@ -191,6 +191,8 @@ add_test_re(joukowsky_disk)
 add_test_re(channel_kwsst)
 add_test_re(channel_kwsstiddes)
 add_test_re(channel_godunov_laminar)
+add_test_re(halfchannel_zerogradient)
+add_test_re(halfchannel_symmetricwall)
 add_test_re(ekman_spiral)
 add_test_re(rayleigh_taylor_godunov)
 add_test_re(rayleigh_taylor_mol)

test/test_files/halfchannel_symmetricwall/halfchannel_symmetricwall.inp

@@ -0,0 +1,57 @@
+# This is a 2D poiseuille flow when run to steady state
+
+#¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨#
+#            SIMULATION STOP            #
+#.......................................#
+time.stop_time               = 5.00     # Max (simulated) time to evolve
+time.max_step                = -1000       # Max number of time steps
+
+#¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨#
+#         TIME STEP COMPUTATION         #
+#.......................................#
+time.fixed_dt         =   0.005       # Use this constant dt if > 0
+time.cfl              =   0.95        # CFL factor
+
+#¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨#
+#            INPUT AND OUTPUT           #
+#.......................................#
+time.plot_interval            =  10       # Steps between plot files
+time.checkpoint_interval      =  -100       # Steps between checkpoint files
+
+#¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨#
+#               PHYSICS                 #
+#.......................................#
+incflo.density        = 1.0             # Reference density
+incflo.use_godunov    = 1
+incflo.diffusion_type = 1
+incflo.godunov_type = "weno_z"
+
+transport.viscosity = 1.0
+turbulence.model = Laminar
+
+ICNS.source_terms = BodyForce
+BodyForce.magnitude = 12.0 0 0
+incflo.physics = ChannelFlow
+ChannelFlow.density = 1.0
+ChannelFlow.Mean_Velocity = 1.0
+ChannelFlow.half_channel = true
+
+io.output_default_variables = 1
+
+#¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨#
+#        ADAPTIVE MESH REFINEMENT       #
+#.......................................#
+amr.n_cell              = 16 8 4     # Grid cells at coarsest AMRlevel
+amr.blocking_factor     = 4
+amr.max_level           = 0           # Max AMR level in hierarchy
+
+#¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨#
+#              GEOMETRY                 #
+#.......................................#
+geometry.prob_lo        =   0.0 0.0 0.0  # Lo corner coordinates
+geometry.prob_hi        =   1.0 0.5 0.25  # Hi corner coordinates
+geometry.is_periodic    =   1   0   1   # Periodicity x y z (0/1)
+
+# Boundary conditions
+ylo.type =   "no_slip_wall"
+yhi.type =   "symmetric_wall"

test/test_files/halfchannel_zerogradient/halfchannel_zerogradient.inp

@@ -0,0 +1,57 @@
+# This is a 2D poiseuille flow when run to steady state
+
+#¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨#
+#            SIMULATION STOP            #
+#.......................................#
+time.stop_time               = 5.00     # Max (simulated) time to evolve
+time.max_step                = -1000       # Max number of time steps
+
+#¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨#
+#         TIME STEP COMPUTATION         #
+#.......................................#
+time.fixed_dt         =   0.005       # Use this constant dt if > 0
+time.cfl              =   0.95        # CFL factor
+
+#¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨#
+#            INPUT AND OUTPUT           #
+#.......................................#
+time.plot_interval            =  10       # Steps between plot files
+time.checkpoint_interval      =  -100       # Steps between checkpoint files
+
+#¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨#
+#               PHYSICS                 #
+#.......................................#
+incflo.density        = 1.0             # Reference density
+incflo.use_godunov    = 1
+incflo.diffusion_type = 1
+incflo.godunov_type = "weno_z"
+
+transport.viscosity = 1.0
+turbulence.model = Laminar
+
+ICNS.source_terms = BodyForce
+BodyForce.magnitude = 12.0 0 0
+incflo.physics = ChannelFlow
+ChannelFlow.density = 1.0
+ChannelFlow.Mean_Velocity = 1.0
+ChannelFlow.half_channel = true
+
+io.output_default_variables = 1
+
+#¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨#
+#        ADAPTIVE MESH REFINEMENT       #
+#.......................................#
+amr.n_cell              = 16 8 4     # Grid cells at coarsest AMRlevel
+amr.blocking_factor     = 4
+amr.max_level           = 0           # Max AMR level in hierarchy
+
+#¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨#
+#              GEOMETRY                 #
+#.......................................#
+geometry.prob_lo        =   0.0 0.0 0.0  # Lo corner coordinates
+geometry.prob_hi        =   1.0 0.5 0.25  # Hi corner coordinates
+geometry.is_periodic    =   1   0   1   # Periodicity x y z (0/1)
+
+# Boundary conditions
+ylo.type =   "no_slip_wall"
+yhi.type =   "zero_gradient"

unit_tests/turbulence/CMakeLists.txt

@@ -3,4 +3,5 @@ target_sources(${amr_wind_unit_test_exe_name}
 
   test_turbulence_init.cpp
   test_turbulence_LES.cpp
+  test_turbulence_LES_bc.cpp
   )