Inflow-outflow boundary condition (Exawind#1066)

---------

Co-authored-by: Marc T. Henry de Frahan <[email protected]>

amr-wind/boundary_conditions/BCInterface.H

@@ -49,7 +49,7 @@ public:
     virtual void operator()(const amrex::Real value = 0.0);
 
     //! User-defined functions for Dirichlet-type boundaries
-    std::pair<const std::string, const std::string> get_dirichlet_udfs();
+    amrex::Array<const std::string, 3> get_dirichlet_udfs();
 
 protected:
     //! Setup AMReX mathematical BC types

amr-wind/boundary_conditions/BCInterface.cpp

@@ -1,6 +1,7 @@
 #include "amr-wind/boundary_conditions/BCInterface.H"
 #include "amr-wind/core/FieldRepo.H"
 #include "amr-wind/boundary_conditions/FixedGradientBC.H"
+#include "amr-wind/boundary_conditions/MassInflowOutflowBC.H"
 #include "AMReX_ParmParse.H"
 
 namespace amr_wind {
@@ -64,6 +65,8 @@ void BCIface::read_bctype()
             ibctype[ori] = BC::pressure_outflow;
         } else if ((bcstr == "mass_inflow") || (bcstr == "mi")) {
             ibctype[ori] = BC::mass_inflow;
+        } else if ((bcstr == "mass_inflow_outflow") || (bcstr == "mio")) {
+            ibctype[ori] = BC::mass_inflow_outflow;
         } else if ((bcstr == "no_slip_wall") || (bcstr == "nsw")) {
             ibctype[ori] = BC::no_slip_wall;
         } else if ((bcstr == "slip_wall") || (bcstr == "sw")) {
@@ -100,18 +103,27 @@ void BCIface::set_bcfuncs()
         if (bct == BC::fixed_gradient) {
             m_field.register_custom_bc<FixedGradientBC>(ori);
         }
+
+        if ((m_field.name() == "velocity") // only velocity for now
+            && (bct == BC::mass_inflow_outflow)) {
+
+            m_field.register_custom_bc<MassInflowOutflowBC>(ori);
+        }
     }
 }
 
-std::pair<const std::string, const std::string> BCIface::get_dirichlet_udfs()
+amrex::Array<const std::string, 3> BCIface::get_dirichlet_udfs()
 {
     const auto& fname = m_field.name();
     const auto& bctype = m_field.bc_type();
     const std::string inflow_key = fname + ".inflow_type";
+    const std::string inflow_outflow_key = fname + ".inflow_outflow_type";
     const std::string wall_key = fname + ".wall_type";
     std::string inflow_udf{"ConstDirichlet"};
+    std::string inflow_outflow_udf{"ConstDirichlet"};
     std::string wall_udf{"ConstDirichlet"};
     bool has_inflow_udf = false;
+    bool has_inflow_outflow_udf = false;
     bool has_wall_udf = false;
 
     for (amrex::OrientationIter oit; oit != nullptr; ++oit) {
@@ -135,6 +147,22 @@ std::pair<const std::string, const std::string> BCIface::get_dirichlet_udfs()
             }
         }
 
+        if (bct == BC::mass_inflow_outflow) {
+            if (pp.contains(inflow_outflow_key.c_str())) {
+                std::string val;
+                pp.get(inflow_outflow_key.c_str(), val);
+
+                if (has_inflow_outflow_udf && (inflow_outflow_udf != val)) {
+                    amrex::Abort(
+                        "BC: Inflow-outflow UDF must be same for all "
+                        "inflow-outflow faces");
+                } else {
+                    inflow_outflow_udf = val;
+                    has_inflow_outflow_udf = true;
+                }
+            }
+        }
+
         if (bct == BC::slip_wall) {
             if (pp.contains(wall_key.c_str())) {
                 std::string val;
@@ -151,7 +179,7 @@ std::pair<const std::string, const std::string> BCIface::get_dirichlet_udfs()
         }
     }
 
-    return {inflow_udf, wall_udf};
+    return {inflow_udf, inflow_outflow_udf, wall_udf};
 }
 
 void BCVelocity::set_bcrec()
@@ -194,6 +222,15 @@ void BCVelocity::set_bcrec()
             }
             break;
 
+        case BC::mass_inflow_outflow:
+            m_field.set_inout_bndry();
+            if (side == amrex::Orientation::low) {
+                set_bcrec_lo(dir, amrex::BCType::direction_dependent);
+            } else {
+                set_bcrec_hi(dir, amrex::BCType::direction_dependent);
+            }
+            break;
+
         case BC::slip_wall:
         case BC::wall_model:
             if (side == amrex::Orientation::low) {
@@ -287,6 +324,15 @@ void BCScalar::set_bcrec()
             }
             break;
 
+        case BC::mass_inflow_outflow:
+            m_field.set_inout_bndry();
+            if (side == amrex::Orientation::low) {
+                set_bcrec_lo(dir, amrex::BCType::direction_dependent);
+            } else {
+                set_bcrec_hi(dir, amrex::BCType::direction_dependent);
+            }
+            break;
+
         case BC::slip_wall:
         case BC::wall_model:
         case BC::fixed_gradient:
@@ -310,15 +356,18 @@ void BCScalar::read_values()
     auto& bcval = m_field.bc_values();
     const int ndim = m_field.num_comp();
     const auto udfs = get_dirichlet_udfs();
-    const auto const_dirichlet_inflow = udfs.first == "ConstDirichlet";
+    const auto const_dirichlet_inflow = udfs[0] == "ConstDirichlet";
+    const auto const_dirichlet_inflow_outflow = udfs[1] == "ConstDirichlet";
 
     for (amrex::OrientationIter oit; oit != nullptr; ++oit) {
         auto ori = oit();
         const auto& bcid = bcnames[ori];
         const auto bct = bctype[ori];
 
         amrex::ParmParse pp(bcid);
-        if ((bct == BC::mass_inflow) && (const_dirichlet_inflow)) {
+        if (((bct == BC::mass_inflow) && (const_dirichlet_inflow)) ||
+            ((bct == BC::mass_inflow_outflow) &&
+             (const_dirichlet_inflow_outflow))) {
             pp.getarr(fname.c_str(), bcval[ori], 0, ndim);
         } else {
             pp.queryarr(fname.c_str(), bcval[ori], 0, ndim);

amr-wind/boundary_conditions/CMakeLists.txt

@@ -4,6 +4,7 @@ target_sources(${amr_wind_lib_name}
       BCInterface.cpp
       FixedGradientBC.cpp
       scalar_bcs.cpp
+      MassInflowOutflowBC.cpp
    )
 
 add_subdirectory(wall_models)

amr-wind/boundary_conditions/MassInflowOutflowBC.H

@@ -0,0 +1,31 @@
+#ifndef MASSINFLOWOUTFLOWBC_H
+#define MASSINFLOWOUTFLOWBC_H
+
+#include "amr-wind/core/FieldBCOps.H"
+#include "amr-wind/core/FieldRepo.H"
+
+#include "AMReX_Orientation.H"
+
+namespace amr_wind {
+
+/** Custom Neumann fills for the inflow-outflow BC
+ *  \ingroup field_bc
+ *
+ *  Used to fill the outflow boundary cells for mass-inflow-outflow BC
+ */
+class MassInflowOutflowBC : public FieldBCIface
+{
+public:
+    MassInflowOutflowBC(Field& field, amrex::Orientation ori);
+
+    void operator()(Field& field, const FieldState /*rho_state*/) override;
+
+private:
+    Field& m_field;
+
+    amrex::Orientation m_ori;
+};
+
+} // namespace amr_wind
+
+#endif /* MASSINFLOWOUTFLOWBC_H */

amr-wind/boundary_conditions/MassInflowOutflowBC.cpp

@@ -0,0 +1,75 @@
+
+#include "amr-wind/boundary_conditions/MassInflowOutflowBC.H"
+
+namespace amr_wind {
+
+MassInflowOutflowBC::MassInflowOutflowBC(Field& field, amrex::Orientation ori)
+    : m_field(field), m_ori(ori)
+{}
+
+void MassInflowOutflowBC::operator()(
+    Field& /*field*/, const FieldState /*rho_state*/)
+{
+    const auto& repo = m_field.repo();
+    const auto& velocity = repo.get_field("velocity");
+    const int ncomp = m_field.num_comp();
+    const int idim = m_ori.coordDir();
+    const auto islow = m_ori.isLow();
+    const auto ishigh = m_ori.isHigh();
+    const int nlevels = m_field.repo().num_active_levels();
+    const amrex::IntVect iv_dir = {
+        static_cast<int>(idim == 0), static_cast<int>(idim == 1),
+        static_cast<int>(idim == 2)};
+
+    for (int lev = 0; lev < nlevels; ++lev) {
+        const auto& domain = repo.mesh().Geom(lev).Domain();
+
+        amrex::MFItInfo mfi_info{};
+        if (amrex::Gpu::notInLaunchRegion()) {
+            mfi_info.SetDynamic(true);
+        }
+#ifdef AMREX_USE_OMP
+#pragma omp parallel if (amrex::Gpu::notInLaunchRegion())
+#endif
+        for (amrex::MFIter mfi(m_field(lev), mfi_info); mfi.isValid(); ++mfi) {
+            auto bx = mfi.validbox();
+            bx.grow(
+                {static_cast<int>(idim != 0), static_cast<int>(idim != 1),
+                 static_cast<int>(idim != 2)});
+            const auto& bc_a = m_field(lev).array(mfi);
+            const auto& vel = velocity(lev).array(mfi);
+
+            if (islow && (bx.smallEnd(idim) == domain.smallEnd(idim))) {
+                amrex::ParallelFor(
+                    amrex::bdryLo(bx, idim),
+                    [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept {
+                        const amrex::IntVect iv = {i, j, k};
+                        const amrex::IntVect ivm = iv - iv_dir;
+                        if (vel(ivm[0], ivm[1], ivm[2], idim) < 0) {
+                            for (int n = 0; n < ncomp; n++) {
+                                bc_a(ivm[0], ivm[1], ivm[2], n) =
+                                    bc_a(i, j, k, n);
+                            }
+                        }
+                    });
+            }
+
+            if (ishigh && (bx.bigEnd(idim) == domain.bigEnd(idim))) {
+                amrex::ParallelFor(
+                    amrex::bdryHi(bx, idim),
+                    [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept {
+                        const amrex::IntVect iv = {i, j, k};
+                        const amrex::IntVect ivm = iv - iv_dir;
+                        if (vel(i, j, k, idim) > 0) {
+                            for (int n = 0; n < ncomp; n++) {
+                                bc_a(i, j, k, n) =
+                                    bc_a(ivm[0], ivm[1], ivm[2], n);
+                            }
+                        }
+                    });
+            }
+        }
+    }
+}
+
+} // namespace amr_wind

amr-wind/boundary_conditions/scalar_bcs.H

@@ -29,7 +29,7 @@ void register_scalar_dirichlet(
     Field& field,
     const amrex::AmrCore& mesh,
     const SimTime& time,
-    const std::pair<const std::string, const std::string>& udfs);
+    const amrex::Array<const std::string, 3>& udfs);
 
 } // namespace amr_wind::scalar_bc
 

amr-wind/boundary_conditions/scalar_bcs.cpp

@@ -5,12 +5,15 @@ void register_scalar_dirichlet(
     Field& field,
     const amrex::AmrCore& mesh,
     const SimTime& time,
-    const std::pair<const std::string, const std::string>& udfs)
+    const amrex::Array<const std::string, 3>& udfs)
 {
-    const std::string inflow_udf = udfs.first;
-    const std::string wall_udf = udfs.second;
+    const std::string& inflow_udf = udfs[0];
+    const std::string& inflow_outflow_udf = udfs[1];
+    const std::string& wall_udf = udfs[2];
 
-    if ((inflow_udf == "ConstDirichlet") && (wall_udf == "ConstDirichlet")) {
+    if ((inflow_udf == "ConstDirichlet") &&
+        (inflow_outflow_udf == "ConstDirichlet") &&
+        (wall_udf == "ConstDirichlet")) {
         return;
     }
 
@@ -19,7 +22,14 @@ void register_scalar_dirichlet(
             "Scalar BC: Only constant dirichlet supported for Wall BC");
     }
 
-    register_inflow_scalar_dirichlet<ConstDirichlet>(
-        field, inflow_udf, mesh, time);
+    if (inflow_udf != "ConstDirichlet") {
+        register_inflow_scalar_dirichlet<ConstDirichlet>(
+            field, inflow_udf, mesh, time);
+    }
+
+    if (inflow_outflow_udf != "ConstDirichlet") {
+        register_inflow_scalar_dirichlet<ConstDirichlet>(
+            field, inflow_outflow_udf, mesh, time);
+    }
 }
 } // namespace amr_wind::scalar_bc

amr-wind/boundary_conditions/velocity_bcs.H

@@ -10,6 +10,7 @@
 #include "amr-wind/physics/udfs/BurggrafLid.H"
 #include "amr-wind/physics/udfs/Rankine.H"
 #include "amr-wind/physics/udfs/CustomVelocity.H"
+#include "amr-wind/physics/udfs/TwoLayer.H"
 
 namespace amr_wind::vel_bc {
 
@@ -40,6 +41,10 @@ void register_inflow_vel_dirichlet(
         using InflowOp = BCOpCreator<udf::CustomVelocity, WallOp>;
         field.register_fill_patch_op<FieldFillPatchOps<InflowOp>>(
             mesh, time, InflowOp(field));
+    } else if (inflow_udf == "TwoLayer") {
+        using InflowOp = BCOpCreator<udf::TwoLayer, WallOp>;
+        field.register_fill_patch_op<FieldFillPatchOps<InflowOp>>(
+            mesh, time, InflowOp(field));
     } else {
         amrex::Abort("Velocity BC: Invalid dirichlet BC type = " + inflow_udf);
     }
@@ -49,12 +54,15 @@ void register_velocity_dirichlet(
     Field& field,
     const amrex::AmrCore& mesh,
     const SimTime& time,
-    const std::pair<const std::string, const std::string>& udfs)
+    const amrex::Array<const std::string, 3>& udfs)
 {
-    const std::string inflow_udf = udfs.first;
-    const std::string wall_udf = udfs.second;
+    const std::string& inflow_udf = udfs[0];
+    const std::string& inflow_outflow_udf = udfs[1];
+    const std::string& wall_udf = udfs[2];
 
-    if ((inflow_udf == "ConstDirichlet") && (wall_udf == "ConstDirichlet")) {
+    if ((inflow_udf == "ConstDirichlet") &&
+        (inflow_outflow_udf == "ConstDirichlet") &&
+        (wall_udf == "ConstDirichlet")) {
         return;
     }
 
@@ -63,8 +71,15 @@ void register_velocity_dirichlet(
             "Velocity BC: Only constant dirichlet supported for Wall BC");
     }
 
-    register_inflow_vel_dirichlet<ConstDirichlet>(
-        field, inflow_udf, mesh, time);
+    if (inflow_udf != "ConstDirichlet") {
+        register_inflow_vel_dirichlet<ConstDirichlet>(
+            field, inflow_udf, mesh, time);
+    }
+
+    if (inflow_outflow_udf != "ConstDirichlet") {
+        register_inflow_vel_dirichlet<ConstDirichlet>(
+            field, inflow_outflow_udf, mesh, time);
+    }
 }
 
 } // namespace amr_wind::vel_bc

amr-wind/core/Field.H

@@ -396,6 +396,12 @@ public:
     inline bool& in_uniform_space() { return m_mesh_mapped; }
     inline bool in_uniform_space() const { return m_mesh_mapped; }
 
+    //! Check if any of the boundaries is a mass-inflow-outflow
+    inline bool has_inout_bndry() const { return m_inout_bndry; }
+
+    //! Set the inout_bndry flag
+    void set_inout_bndry() { m_inout_bndry = true; }
+
 protected:
     Field(
         FieldRepo& repo,
@@ -425,6 +431,9 @@ protected:
 
     //! Flag to track mesh mapping (to uniform space) of field
     bool m_mesh_mapped{false};
+
+    //! Flag to indicate whether any of the boundaries is mass-inflow-outflow
+    bool m_inout_bndry{false};
 };
 
 } // namespace amr_wind