diff --git a/amr-wind/boundary_conditions/BCInterface.H b/amr-wind/boundary_conditions/BCInterface.H index b6d815e82d..45c6e4e0ef 100644 --- a/amr-wind/boundary_conditions/BCInterface.H +++ b/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 get_dirichlet_udfs(); + amrex::Array get_dirichlet_udfs(); protected: //! Setup AMReX mathematical BC types diff --git a/amr-wind/boundary_conditions/BCInterface.cpp b/amr-wind/boundary_conditions/BCInterface.cpp index b6ca50ab21..b363f38d0f 100644 --- a/amr-wind/boundary_conditions/BCInterface.cpp +++ b/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(ori); } + + if ((m_field.name() == "velocity") // only velocity for now + && (bct == BC::mass_inflow_outflow)) { + + m_field.register_custom_bc(ori); + } } } -std::pair BCIface::get_dirichlet_udfs() +amrex::Array 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 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 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,7 +356,8 @@ 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(); @@ -318,7 +365,9 @@ void BCScalar::read_values() 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); diff --git a/amr-wind/boundary_conditions/CMakeLists.txt b/amr-wind/boundary_conditions/CMakeLists.txt index bdc98bd3f6..26af39ca57 100644 --- a/amr-wind/boundary_conditions/CMakeLists.txt +++ b/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) diff --git a/amr-wind/boundary_conditions/MassInflowOutflowBC.H b/amr-wind/boundary_conditions/MassInflowOutflowBC.H new file mode 100644 index 0000000000..ba8d2c53bf --- /dev/null +++ b/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 */ diff --git a/amr-wind/boundary_conditions/MassInflowOutflowBC.cpp b/amr-wind/boundary_conditions/MassInflowOutflowBC.cpp new file mode 100644 index 0000000000..1262c096b5 --- /dev/null +++ b/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(idim == 0), static_cast(idim == 1), + static_cast(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(idim != 0), static_cast(idim != 1), + static_cast(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 \ No newline at end of file diff --git a/amr-wind/boundary_conditions/scalar_bcs.H b/amr-wind/boundary_conditions/scalar_bcs.H index d27adaac6c..1a4ebd2e10 100644 --- a/amr-wind/boundary_conditions/scalar_bcs.H +++ b/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& udfs); + const amrex::Array& udfs); } // namespace amr_wind::scalar_bc diff --git a/amr-wind/boundary_conditions/scalar_bcs.cpp b/amr-wind/boundary_conditions/scalar_bcs.cpp index e2e9adbf81..3158e4fa8e 100644 --- a/amr-wind/boundary_conditions/scalar_bcs.cpp +++ b/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& udfs) + const amrex::Array& 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( - field, inflow_udf, mesh, time); + if (inflow_udf != "ConstDirichlet") { + register_inflow_scalar_dirichlet( + field, inflow_udf, mesh, time); + } + + if (inflow_outflow_udf != "ConstDirichlet") { + register_inflow_scalar_dirichlet( + field, inflow_outflow_udf, mesh, time); + } } } // namespace amr_wind::scalar_bc diff --git a/amr-wind/boundary_conditions/velocity_bcs.H b/amr-wind/boundary_conditions/velocity_bcs.H index 6d36fa1e83..a2a901e8c8 100644 --- a/amr-wind/boundary_conditions/velocity_bcs.H +++ b/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; field.register_fill_patch_op>( mesh, time, InflowOp(field)); + } else if (inflow_udf == "TwoLayer") { + using InflowOp = BCOpCreator; + field.register_fill_patch_op>( + 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& udfs) + const amrex::Array& 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( - field, inflow_udf, mesh, time); + if (inflow_udf != "ConstDirichlet") { + register_inflow_vel_dirichlet( + field, inflow_udf, mesh, time); + } + + if (inflow_outflow_udf != "ConstDirichlet") { + register_inflow_vel_dirichlet( + field, inflow_outflow_udf, mesh, time); + } } } // namespace amr_wind::vel_bc diff --git a/amr-wind/core/Field.H b/amr-wind/core/Field.H index beaf535844..d00261ded6 100644 --- a/amr-wind/core/Field.H +++ b/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 diff --git a/amr-wind/core/FieldBCOps.H b/amr-wind/core/FieldBCOps.H index bdd2f43121..4f0d2c3dc8 100644 --- a/amr-wind/core/FieldBCOps.H +++ b/amr-wind/core/FieldBCOps.H @@ -166,9 +166,10 @@ struct DirichletOp auto ori = oit(); const int idir = ori.coordDir(); - // Check if this boundary is a dirichlet type + // Check if this boundary is a dirichlet or mixed type const auto bctyp = (ori.isLow() ? bc.lo(idir) : bc.hi(idir)); - if (bctyp != amrex::BCType::ext_dir) { + if ((bctyp != amrex::BCType::ext_dir) && + (bctyp != amrex::BCType::direction_dependent)) { continue; } @@ -182,7 +183,8 @@ struct DirichletOp continue; } - if (m_bc_type[ori] == BC::mass_inflow) { + if ((m_bc_type[ori] == BC::mass_inflow) || + (m_bc_type[ori] == BC::mass_inflow_outflow)) { m_inflow_op( iv, field, geom, time, ori, n, dcomp, orig_comp); } else { diff --git a/amr-wind/core/FieldFillPatchOps.H b/amr-wind/core/FieldFillPatchOps.H index 56bde83b2d..90561875fa 100644 --- a/amr-wind/core/FieldFillPatchOps.H +++ b/amr-wind/core/FieldFillPatchOps.H @@ -370,7 +370,8 @@ public: for (amrex::OrientationIter oit; oit != nullptr; ++oit) { const auto ori = oit(); - if (bctype[ori] != BC::mass_inflow) { + if ((bctype[ori] != BC::mass_inflow) && + (bctype[ori] != BC::mass_inflow_outflow)) { continue; } @@ -414,7 +415,8 @@ public: for (amrex::OrientationIter oit; oit != nullptr; ++oit) { const auto ori = oit(); - if (bctype[ori] != BC::mass_inflow) { + if ((bctype[ori] != BC::mass_inflow) && + (bctype[ori] != BC::mass_inflow_outflow)) { continue; } @@ -441,12 +443,28 @@ public: } const auto& marr = mfab[mfi].array(); + amrex::FArrayBox tmp_fab( + bx, mfab.nComp(), amrex::The_Async_Arena()); + amrex::Array4 const& tmp_marr = + tmp_fab.array(); + amrex::ParallelFor( bx, [=] AMREX_GPU_DEVICE( const int i, const int j, const int k) noexcept { bcfunc.set_inflow( - {i, j, k}, marr, gdata, time, ori, 0, 0, fdir); + {i, j, k}, tmp_marr, gdata, time, ori, 0, 0, + fdir); + + const bool is_outflow = + (ori.isLow() && tmp_marr(i, j, k) < 0) || + (ori.isHigh() && tmp_marr(i, j, k) > 0); + if ((bctype[ori] == BC::mass_inflow_outflow) && + is_outflow) { + marr(i, j, k) = marr(i, j, k); + } else { + marr(i, j, k) = tmp_marr(i, j, k); + } }); } } diff --git a/amr-wind/diffusion/incflo_diffusion.cpp b/amr-wind/diffusion/incflo_diffusion.cpp index a26afaac53..da951dfa4c 100644 --- a/amr-wind/diffusion/incflo_diffusion.cpp +++ b/amr-wind/diffusion/incflo_diffusion.cpp @@ -29,6 +29,7 @@ Vector> get_diffuse_tensor_bc( } case BC::wave_generation: case BC::mass_inflow: + case BC::mass_inflow_outflow: case BC::no_slip_wall: { // All three components are Dirichlet r[0][dir] = LinOpBCType::Dirichlet; @@ -90,6 +91,7 @@ get_diffuse_scalar_bc(amr_wind::Field& scalar, Orientation::Side side) noexcept } case BC::wave_generation: case BC::mass_inflow: + case BC::mass_inflow_outflow: case BC::no_slip_wall: { r[dir] = LinOpBCType::Dirichlet; break; diff --git a/amr-wind/equation_systems/icns/icns_advection.H b/amr-wind/equation_systems/icns/icns_advection.H index a72d6d24d0..934ae092fd 100644 --- a/amr-wind/equation_systems/icns/icns_advection.H +++ b/amr-wind/equation_systems/icns/icns_advection.H @@ -49,6 +49,10 @@ private: void init_projector(const FaceFabPtrVec& /*beta*/) noexcept; void init_projector(const amrex::Real /*beta*/) noexcept; + void enforce_inout_solvability( + const amrex::Vector>& + a_umac) noexcept; + FieldRepo& m_repo; PhysicsMgr& m_phy_mgr; std::unique_ptr m_mac_proj; diff --git a/amr-wind/equation_systems/icns/icns_advection.cpp b/amr-wind/equation_systems/icns/icns_advection.cpp index e381e54c36..e191801e3c 100644 --- a/amr-wind/equation_systems/icns/icns_advection.cpp +++ b/amr-wind/equation_systems/icns/icns_advection.cpp @@ -64,6 +64,17 @@ MacProjOp::MacProjOp( m_has_overset = m_has_overset && !disable_ovst_mac; } +void MacProjOp::enforce_inout_solvability( + const amrex::Vector>& + a_umac) noexcept +{ + auto& velocity = m_repo.get_field("velocity"); + amrex::BCRec const* bc_type = velocity.bcrec_device().data(); + const amrex::Vector& geom = m_repo.mesh().Geom(); + + HydroUtils::enforceInOutSolvability(a_umac, bc_type, geom); +} + void MacProjOp::init_projector(const MacProjOp::FaceFabPtrVec& beta) noexcept { m_mac_proj = std::make_unique( @@ -266,6 +277,12 @@ void MacProjOp::operator()(const FieldState fstate, const amrex::Real dt) } } + const bool has_inout_bndry = + (m_repo.get_field("velocity")).has_inout_bndry(); + if (has_inout_bndry) { + enforce_inout_solvability(mac_vec); + } + m_mac_proj->setUMAC(mac_vec); if (m_has_overset) { diff --git a/amr-wind/incflo_enums.H b/amr-wind/incflo_enums.H index b73ff06b11..888ff0ae70 100644 --- a/amr-wind/incflo_enums.H +++ b/amr-wind/incflo_enums.H @@ -7,6 +7,7 @@ enum struct BC { pressure_inflow, pressure_outflow, mass_inflow, + mass_inflow_outflow, zero_gradient, no_slip_wall, slip_wall, diff --git a/amr-wind/physics/udfs/CMakeLists.txt b/amr-wind/physics/udfs/CMakeLists.txt index f29413ba3b..1ab72fd445 100644 --- a/amr-wind/physics/udfs/CMakeLists.txt +++ b/amr-wind/physics/udfs/CMakeLists.txt @@ -8,4 +8,5 @@ target_sources(${amr_wind_lib_name} Rankine.cpp CustomVelocity.cpp CustomScalar.cpp + TwoLayer.cpp ) diff --git a/amr-wind/physics/udfs/TwoLayer.H b/amr-wind/physics/udfs/TwoLayer.H new file mode 100644 index 0000000000..5de0418d5a --- /dev/null +++ b/amr-wind/physics/udfs/TwoLayer.H @@ -0,0 +1,67 @@ +#ifndef TWO_LAYER_H +#define TWO_LAYER_H + +#include "AMReX_Geometry.H" +#include "AMReX_Gpu.H" + +namespace amr_wind { + +class Field; + +namespace udf { + +struct TwoLayer +{ + struct DeviceOp + { + // velocities of the top and bottom layer respectively + amrex::GpuArray top_vel = {0.0}; + amrex::GpuArray bottom_vel = {0.0}; + amrex::Real z_part{0.5}; // z-coordinate that divides top and bottom + amrex::Real init_perturb{1.0}; // initial perturbation for testing + + AMREX_GPU_DEVICE + inline void operator()( + const amrex::IntVect& iv, + amrex::Array4 const& field, + amrex::GeometryData const& geom, + const amrex::Real time, + amrex::Orientation /* ori */, + const int comp, + const int dcomp, + const int orig_comp) const + { + const auto* problo = geom.ProbLo(); + const auto* dx = geom.CellSize(); + const auto z = problo[2] + (iv[2] + 0.5) * dx[2]; + + amrex::GpuArray vel; + if (z >= z_part) { + vel = {top_vel[0], top_vel[1], 0.0}; + } else { + vel = {bottom_vel[0], bottom_vel[1], 0.0}; + } + + if (time == 0.0) { + vel[0] *= init_perturb; + vel[1] *= init_perturb; + } + + field(iv[0], iv[1], iv[2], dcomp + comp) = vel[orig_comp + comp]; + } + }; + using DeviceType = DeviceOp; + + static std::string identifier() { return "TwoLayer"; } + + explicit TwoLayer(const Field& fld); + + DeviceType device_instance() const { return m_op; } + + DeviceOp m_op; +}; + +} // namespace udf +} // namespace amr_wind + +#endif /* TWO_LAYER_H */ diff --git a/amr-wind/physics/udfs/TwoLayer.cpp b/amr-wind/physics/udfs/TwoLayer.cpp new file mode 100644 index 0000000000..050b7b590a --- /dev/null +++ b/amr-wind/physics/udfs/TwoLayer.cpp @@ -0,0 +1,38 @@ +#include "amr-wind/physics/udfs/TwoLayer.H" +#include "amr-wind/core/Field.H" +#include "amr-wind/core/FieldRepo.H" +#include "amr-wind/core/vs/vector.H" +#include "amr-wind/equation_systems/icns/icns.H" + +#include "AMReX_ParmParse.H" + +namespace amr_wind::udf { + +TwoLayer::TwoLayer(const Field& fld) +{ + // For a 2-layer flow with a top and a bottom layer + // divided at a specific z-coordinate (default: 0.5) + // and an optional initial perturbation (default: 1.0) + + { + const int ncomp = fld.num_comp(); + amrex::ParmParse pp("TwoLayer"); + + amrex::Vector top_vel(0.0, ncomp); + amrex::Vector bottom_vel(0.0, ncomp); + pp.getarr("top_vel", top_vel); + pp.getarr("bottom_vel", bottom_vel); + + AMREX_ALWAYS_ASSERT(top_vel.size() == ncomp); + AMREX_ALWAYS_ASSERT(bottom_vel.size() == ncomp); + for (int i = 0; i < ncomp; ++i) { + m_op.top_vel[i] = top_vel[i]; + m_op.bottom_vel[i] = bottom_vel[i]; + } + + pp.query("init_perturb", m_op.init_perturb); + pp.query("z_partition", m_op.z_part); + } +} + +} // namespace amr_wind::udf diff --git a/amr-wind/physics/udfs/UDF.cpp b/amr-wind/physics/udfs/UDF.cpp index 5372b27976..7a3427de26 100644 --- a/amr-wind/physics/udfs/UDF.cpp +++ b/amr-wind/physics/udfs/UDF.cpp @@ -7,6 +7,7 @@ #include "amr-wind/physics/udfs/Rankine.H" #include "amr-wind/physics/udfs/CustomVelocity.H" #include "amr-wind/physics/udfs/CustomScalar.H" +#include "amr-wind/physics/udfs/TwoLayer.H" #include "AMReX_ParmParse.H" @@ -63,5 +64,6 @@ template class UDFImpl; template class UDFImpl; template class UDFImpl; template class UDFImpl; +template class UDFImpl; } // namespace amr_wind::udf diff --git a/amr-wind/projection/incflo_apply_nodal_projection.cpp b/amr-wind/projection/incflo_apply_nodal_projection.cpp index e8bb944559..ad12059c8c 100644 --- a/amr-wind/projection/incflo_apply_nodal_projection.cpp +++ b/amr-wind/projection/incflo_apply_nodal_projection.cpp @@ -5,6 +5,7 @@ #include "amr-wind/utilities/console_io.H" #include "amr-wind/core/field_ops.H" #include "amr-wind/projection/nodal_projection_ops.H" +#include "hydro_utils.H" using namespace amrex; @@ -46,6 +47,7 @@ amr_wind::nodal_projection::get_projection_bc( r[dir] = LinOpBCType::Dirichlet; break; } + case BC::mass_inflow_outflow: case BC::mass_inflow: { r[dir] = LinOpBCType::inflow; break; @@ -76,6 +78,23 @@ void amr_wind::nodal_projection::apply_dirichlet_vel( }); } +void amr_wind::nodal_projection::enforce_inout_solvability( + amr_wind::Field& velocity, + const Vector& geom, + const int num_levels) +{ + BCRec const* bc_type = velocity.bcrec_device().data(); + Vector> vel_vec(num_levels); + + for (int lev = 0; lev < num_levels; ++lev) { + vel_vec[lev][0] = new MultiFab(velocity(lev), amrex::make_alias, 0, 1); + vel_vec[lev][1] = new MultiFab(velocity(lev), amrex::make_alias, 1, 1); + vel_vec[lev][2] = new MultiFab(velocity(lev), amrex::make_alias, 2, 1); + } + + HydroUtils::enforceInOutSolvability(vel_vec, bc_type, geom, true); +} + /** Perform nodal projection * * Computes the following decomposition: @@ -329,9 +348,22 @@ void incflo::ApplyProjection( if (!proj_for_small_dt and !incremental) { amr_wind::nodal_projection::set_inflow_velocity( m_sim.physics_manager(), velocity, lev, time, *vel[lev], 1); + + // fill periodic boundaries to avoid corner cell issues + vel[lev]->FillBoundary(geom[lev].periodicity()); } } + // Need to apply custom Neumann funcs for inflow-outflow BC + // after setting the inflow vels above + // and then enforce solvability by matching outflow to inflow. + if (!proj_for_small_dt and !incremental and velocity.has_inout_bndry()) { + velocity.apply_bc_funcs(amr_wind::FieldState::New); + + amr_wind::nodal_projection::enforce_inout_solvability( + velocity, m_repo.mesh().Geom(), m_repo.num_active_levels()); + } + if (is_anelastic) { for (int lev = 0; lev <= finest_level; ++lev) { for (int idim = 0; idim < velocity.num_comp(); ++idim) { diff --git a/amr-wind/projection/nodal_projection_ops.H b/amr-wind/projection/nodal_projection_ops.H index 5158fd17e7..c437146463 100644 --- a/amr-wind/projection/nodal_projection_ops.H +++ b/amr-wind/projection/nodal_projection_ops.H @@ -26,6 +26,9 @@ Array get_projection_bc( void apply_dirichlet_vel( amrex::MultiFab& mf_velocity, amrex::iMultiFab& mf_iblank); +void enforce_inout_solvability( + amr_wind::Field& velocity, const Vector& geom, int num_levels); + } // namespace amr_wind::nodal_projection #endif \ No newline at end of file diff --git a/docs/sphinx/user/inputs_Boundary_conditions.rst b/docs/sphinx/user/inputs_Boundary_conditions.rst index eaae52bb51..326714c2e5 100644 --- a/docs/sphinx/user/inputs_Boundary_conditions.rst +++ b/docs/sphinx/user/inputs_Boundary_conditions.rst @@ -1,21 +1,21 @@ Section: Boundary conditions ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ - + This section controls the boundary conditions. Only non-periodic BC's need to be defined here. .. input_param:: xlo.type (or ylo.type, zlo.type, xhi.type, yhi.type, zhi.type) **type:** String, mandatory if non-periodic - - 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, symmetric_wall and wall_model. + + Boundary condition type on the lo (or hi) side of the domain. + Current options are: periodic, pressure_inflow, pressure_outflow, mass_inflow, + mass_inflow_outflow, no_slip_wall, slip_wall, symmetric_wall and wall_model. .. input_param:: xlo.temperature (or ylo.temperature, etc) **type:** Real, optional, default = 0.0 - - Specifies a temperature gradient at the wall, only activated for slip_wall and wall_model BC types. + + Specifies a temperature gradient at the wall, only activated for slip_wall and wall_model BC types. Mass inflow boundary conditions ``````````````````````````````` @@ -42,3 +42,46 @@ If the user wants to define their own boundary conditions, this is done by editi CustomVelocity.foo = 1.0 They do not both need to be defined at the same time. It is the user's responsibility to ensure that the source files are appropriately edited for their use case. Examples of how these files can be edited are found through comparison of the other mass_inflow functions in the `udfs` folder. + +Mass inflow-outflow boundary conditions +``````````````````````````````````````` + +The mass_inflow_outflow boundary condition is designed to handle both inflow and outflow at the same boundary. +For the advection schemes, it implements a Neumann type behavior at the outflow cells and a Dirichlet behavior at the inflow cells. +It uses Neumann conditions for the MAC and nodal projections and +enforces solvability before the projections +by correcting the outflow to match with the inflow within the specified mass_inflow_outflow boundaries. +It uses a Dirichlet condition for the diffusion solver. + +Both the approaches mentioned above for the mass inflow condition, +constant values and UDFs, can be used to specify the boundary values. +The outflow values will be automatically replaced by a value from the interior cell +to enforce the Neumann type behavior. +See the ``freestream_godunov_inout`` test for an example that uses the TwoLayer UDF. +This test involves two z-layers of the flow along opposite x-directions. +The input file options are copied here:: + + geometry.is_periodic = 0 1 0 # Periodic in y + + # Boundary conditions + TwoLayer.bottom_vel = -1.0 0.0 0.0 + TwoLayer.top_vel = 1.0 0.0 0.0 + TwoLayer.init_perturb = 0.9 + TwoLayer.z_partition = 0.5 + + xlo.type = "mass_inflow_outflow" + xlo.density = 1.0 + xlo.velocity.inflow_outflow_type = TwoLayer + + xhi.type = "mass_inflow_outflow" + xhi.density = 1.0 + xhi.velocity.inflow_outflow_type = TwoLayer + + zlo.type = "slip_wall" + zhi.type = "slip_wall" + + +The most applicable use case for this boundary condition is with the +:ref:`amrwind-abl-bndry-io` for flows that change directions +across the vertical coordinate or with time. +The work to integrate this condition with the ABL class is under progress. diff --git a/test/CMakeLists.txt b/test/CMakeLists.txt index 8f29309f07..47cac81ffd 100644 --- a/test/CMakeLists.txt +++ b/test/CMakeLists.txt @@ -248,6 +248,7 @@ add_test_re(rankine) add_test_re(rankine-sym) add_test_re(box_refinement) add_test_re(cylinder_refinement) +add_test_re(freestream_godunov_inout) if(NOT AMR_WIND_ENABLE_CUDA) add_test_re(ctv_godunov_plm) diff --git a/test/test_files/freestream_godunov_inout/freestream_godunov_inout.inp b/test/test_files/freestream_godunov_inout/freestream_godunov_inout.inp new file mode 100644 index 0000000000..44fad0a075 --- /dev/null +++ b/test/test_files/freestream_godunov_inout/freestream_godunov_inout.inp @@ -0,0 +1,70 @@ +# This is a 2D poiseuille flow when run to steady state + +#¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨# +# SIMULATION STOP # +#.......................................# +time.stop_time = 10.0 # Max (simulated) time to evolve +time.max_step = 10 #-1000 # Max number of time steps + +#¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨# +# TIME STEP COMPUTATION # +#.......................................# +time.fixed_dt = 0.02 # Use this constant dt if > 0 +time.cfl = 0.95 # CFL factor + +#¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨# +# INPUT AND OUTPUT # +#.......................................# +time.plot_interval = 1 # 100 # 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 = 0.0 +turbulence.model = Laminar + +ICNS.source_terms = BodyForce +BodyForce.magnitude = 0 0 0 + +incflo.physics = FreeStream +FreeStream.velocity_type = TwoLayer + +io.output_default_variables = 1 + +#¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨# +# ADAPTIVE MESH REFINEMENT # +#.......................................# +amr.n_cell = 16 8 16 # Grid cells at coarsest AMRlevel +amr.blocking_factor = 4 +amr.max_level = 0 # Max AMR level in hierarchy +fabarray.mfiter_tile_size = 1024 1024 1024 + +#¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨# +# GEOMETRY # +#.......................................# +geometry.prob_lo = 0.0 0.0 0.0 # Lo corner coordinates +geometry.prob_hi = 1.0 0.5 1.0 # Hi corner coordinates +geometry.is_periodic = 0 1 0 # Periodicity x y z (0/1) + +# Boundary conditions +TwoLayer.bottom_vel = -1.0 0.0 0.0 +TwoLayer.top_vel = 1.0 0.0 0.0 +TwoLayer.init_perturb = 0.9 +TwoLayer.z_partition = 0.5 + +xlo.type = "mass_inflow_outflow" +xlo.density = 1.0 +xlo.velocity.inflow_outflow_type = TwoLayer + +xhi.type = "mass_inflow_outflow" +xhi.density = 1.0 +xhi.velocity.inflow_outflow_type = TwoLayer + +zlo.type = "slip_wall" +zhi.type = "slip_wall"