Add inlet outlet bc

src/Control/Inciter/InputDeck/InputDeck.hpp

@@ -47,7 +47,6 @@ using ncomp_t = std::size_t;
 using bclist = tk::TaggedTuple< brigand::list<
   tag::dirichlet,   std::vector< std::size_t >,
   tag::symmetry,    std::vector< std::size_t >,
-  tag::inlet,       std::vector< std::size_t >,
   tag::outlet,      std::vector< std::size_t >,
   tag::farfield,    std::vector< std::size_t >,
   tag::extrapolate, std::vector< std::size_t >,
@@ -128,23 +127,31 @@ using materialList = tk::TaggedTuple< brigand::list<
 
 // Boundary conditions block
 using bcList = tk::TaggedTuple< brigand::list<
-  tag::mesh,           std::vector< std::size_t >,
-  tag::dirichlet,      std::vector< std::size_t >,
-  tag::symmetry,       std::vector< std::size_t >,
-  tag::inlet,          std::vector< std::size_t >,
-  tag::outlet,         std::vector< std::size_t >,
-  tag::farfield,       std::vector< std::size_t >,
-  tag::extrapolate,    std::vector< std::size_t >,
-  tag::noslipwall,     std::vector< std::size_t >,
-  tag::stag_point,     std::vector< tk::real >,
-  tag::radius,         tk::real,
-  tag::velocity,       std::vector< tk::real >,
-  tag::pressure,       tk::real,
-  tag::density,        tk::real,
-  tag::temperature,    tk::real,
+  tag::mesh,        std::vector< std::size_t >,
+  tag::dirichlet,   std::vector< std::size_t >,
+  tag::symmetry,    std::vector< std::size_t >,
+  tag::outlet,      std::vector< std::size_t >,
+  tag::farfield,    std::vector< std::size_t >,
+  tag::extrapolate, std::vector< std::size_t >,
+  tag::noslipwall,  std::vector< std::size_t >,
+  tag::stag_point,  std::vector< tk::real >,
+  tag::radius,      tk::real,
+  tag::velocity,    std::vector< tk::real >,
+  tag::pressure,    tk::real,
+  tag::density,     tk::real,
+  tag::temperature, tk::real,
   tag::mass_fractions, std::vector< tk::real >,
-  tag::materialid,     std::size_t,
-  tag::timedep,        std::vector<
+  tag::materialid,  std::size_t,
+  tag::inlet,       std::vector<
+    tk::TaggedTuple< brigand::list<
+      tag::sideset,      std::vector< uint64_t >,
+      tag::velocity,     std::vector< tk::real >,
+      tag::pressure,     tk::real,
+      tag::temperature,     tk::real,
+      tag::materialid,   std::size_t
+    > >
+  >,
+  tag::timedep,     std::vector<
     tk::TaggedTuple< brigand::list<
       tag::sideset,    std::vector< uint64_t >,
       tag::fn,         std::vector< tk::real >

src/Control/Inciter/InputDeck/LuaParser.cpp

@@ -1087,8 +1087,30 @@ LuaParser::storeInputDeck(
       storeVecIfSpecd< uint64_t >(sol_bc[i+1], "symmetry",
         bc_deck[i].get< tag::symmetry >(), {});
 
-      storeVecIfSpecd< uint64_t >(sol_bc[i+1], "inlet",
-        bc_deck[i].get< tag::inlet >(), {});
+      if (sol_bc[i+1]["inlet"].valid()) {
+        const sol::table& sol_inbc = sol_bc[i+1]["inlet"];
+        auto& inbc_deck = bc_deck[i].get< tag::inlet >();
+        inbc_deck.resize(sol_inbc.size());
+
+        for (std::size_t j=0; j<inbc_deck.size(); ++j) {
+          storeVecIfSpecd< uint64_t >(sol_inbc[j+1], "sideset",
+            inbc_deck[j].get< tag::sideset >(), {});
+
+          storeVecIfSpecd< tk::real >(sol_inbc[j+1], "velocity",
+            inbc_deck[j].get< tag::velocity >(), {0.0, 0.0, 0.0});
+          if (inbc_deck[j].get< tag::velocity >().size() != 3)
+            Throw("Inlet velocity requires 3 components.");
+
+          storeIfSpecd< tk::real >(sol_inbc[j+1], "pressure",
+            inbc_deck[j].get< tag::pressure >(), 0.0);
+
+          storeIfSpecd< tk::real >(sol_inbc[j+1], "temperature",
+            inbc_deck[j].get< tag::temperature >(), 0.0);
+
+          storeIfSpecd< std::size_t >(sol_inbc[j+1], "materialid",
+            inbc_deck[j].get< tag::materialid >(), 1);
+        }
+      }
 
       storeVecIfSpecd< uint64_t >(sol_bc[i+1], "outlet",
         bc_deck[i].get< tag::outlet >(), {});

src/Inciter/Transporter.cpp

@@ -475,14 +475,19 @@ Transporter::matchBCs( std::map< int, std::vector< std::size_t > >& bnd )
 // *****************************************************************************
 {
   // Query side set ids at which BCs assigned for all BC types for all PDEs
-  using bclist = ctr::bclist::Keys;
   std::unordered_set< int > usedsets;
+
+  using bclist = ctr::bclist::Keys;
+  const auto& bcs = g_inputdeck.get< tag::bc >();
   brigand::for_each< bclist >( UserBC( g_inputdeck, usedsets ) );
 
-  // Query side sets of time dependent BCs (since tag::bctimedep is not a part
+  // Query side sets of time dependent and inlet BCs (since these are not a part
   // of tag::bc)
-  const auto& bcs = g_inputdeck.get< tag::bc >();
   for (const auto& bci : bcs) {
+    for (const auto& b : bci.get< tag::inlet >()) {
+      for (auto i : b.get< tag::sideset >())
+        usedsets.insert(static_cast<int>(i));
+    }
     for (const auto& b : bci.get< tag::timedep >()) {
       for (auto i : b.get< tag::sideset >())
         usedsets.insert(static_cast<int>(i));

src/PDE/CompFlow/DGCompFlow.hpp

@@ -72,23 +72,21 @@ class CompFlow {
       // associate boundary condition configurations with state functions, the
       // order in which the state functions listed matters, see ctr::bc::Keys
       brigand::for_each< ctr::bclist::Keys >( ConfigBC( m_bc,
-        // BC State functions
-        { dirichlet
-        , symmetry
-        , invalidBC         // Inlet BC not implemented
-        , invalidBC         // Outlet BC not implemented
-        , farfield
-        , extrapolate
-        , invalidBC },      // No slip wall BC not implemented
-        // BC Gradient functions
-        { noOpGrad
-        , noOpGrad
-        , noOpGrad
-        , noOpGrad
-        , noOpGrad
-        , noOpGrad
-        , noOpGrad }
-        ) );
+      // BC State functions
+      { dirichlet
+      , symmetry
+      , invalidBC         // Outlet BC not implemented
+      , farfield
+      , extrapolate
+      , invalidBC },      // No slip wall BC not implemented
+      // BC Gradient functions
+      { noOpGrad
+      , noOpGrad
+      , noOpGrad
+      , noOpGrad
+      , noOpGrad
+      , noOpGrad }
+      ) );
 
       // EoS initialization
       const auto& matprop =

src/PDE/ConfigureTransport.cpp

@@ -101,8 +101,10 @@ infoTransport( std::map< ctr::PDEType, tk::ncomp_t >& cnt )
     const auto& bcinlet =
       ib.get< tag::inlet >();
     if (!bcinlet.empty())
-      nfo.emplace_back( "Inlet boundary [" + std::to_string( ncomp ) + "]",
-        parameters( bcinlet ) );
+      for (const auto& bndry : bcinlet) {
+        nfo.emplace_back( "Inlet boundary [" + std::to_string( ncomp ) + "]",
+          parameters(bndry.get< tag::sideset >()) );
+    }
 
     const auto& bcoutlet =
       ib.get< tag::outlet >();

src/PDE/DGPDE.cpp

@@ -11,9 +11,37 @@
 // *****************************************************************************
 
 #include "DGPDE.hpp"
+namespace inciter {
+
+void ConfigInletBC( BCStateFn& s,
+                          const tk::StateFn& f,
+                          const tk::StateFn& gf )
+// *****************************************************************************
+// Extract information from input deck on inlet boundary conditions, and append
+// to BC configuration vector
+//! \param[inout] s BC configuration vector
+//! \param[in] f Function to evaluate the left and right solution state at
+//!   boundaries
+//! \param[in] gf Function to evaluate the left and right solution gradients
+//!   at boundaries
+// *****************************************************************************
+{
+  const auto& bc = g_inputdeck.get< tag::bc >();
+  std::vector< std::size_t > v;
+  for (const auto& ib : bc) {
+    const auto& in = ib.get< tag::inlet >();
+    if (!in.empty()) {
+      for (const auto& bndry : in) {
+        const auto& sideset = bndry.get< tag::sideset >();
+        v.insert(v.end(), sideset.begin(), sideset.end());
+        s.push_back( { v, f, gf } );
+      }
+    }
+  }
+}
 
 [[noreturn]] tk::StateFn::result_type
-inciter::invalidBC( ncomp_t, const std::vector< EOS >&,
+invalidBC( ncomp_t, const std::vector< EOS >&,
            const std::vector< tk::real >&, tk::real, tk::real, tk::real,
            tk::real, const std::array< tk::real, 3>& )
 // *****************************************************************************
@@ -24,3 +52,5 @@ inciter::invalidBC( ncomp_t, const std::vector< EOS >&,
   Throw( "Invalid boundary condition set up in input file or the PDE does not "
           "support this BC type" );
 }
+
+}  // inciter::

src/PDE/DGPDE.hpp

@@ -71,6 +71,12 @@ struct ConfigBC {
   }
 };
 
+//! \brief Extract information on inlet BCs, which have a different structure
+//! than other BCs
+void ConfigInletBC( BCStateFn&,
+                    const tk::StateFn&,
+                    const tk::StateFn& );
+
 //! State function for invalid/un-configured boundary conditions
 [[noreturn]] tk::StateFn::result_type
 invalidBC( ncomp_t, const std::vector< EOS >&,

src/PDE/MultiMat/BCFunctions.hpp

@@ -116,6 +116,106 @@ namespace inciter {
     return {{ std::move(ul), std::move(ur) }};
   }
 
+  //! \brief Boundary state function providing the left and right state of a
+  //!   face at inlet boundaries
+  //! \param[in] ncomp Number of scalar components in this PDE system
+  //! \param[in] ul Left (domain-internal) state
+  //! \param[in] fn Unit face normal
+  //! \return Left and right states for all scalar components in this PDE
+  //!   system
+  //! \details The inlet boundary condition specifies a velocity at a
+  //!   sideset and assumes a zero bulk pressure and density gradient
+  //! \note The function signature must follow tk::StateFn
+  static tk::StateFn::result_type
+  inlet( ncomp_t ncomp,
+            const std::vector< EOS >& mat_blk,
+            const std::vector< tk::real >& ul,
+            tk::real, tk::real, tk::real, tk::real,
+            const std::array< tk::real, 3 >& fn )
+  {
+    auto nmat = g_inputdeck.get< tag::multimat, tag::nmat >();
+    const auto& solidx = g_inputdeck.get< tag::matidxmap, tag::solidx >();
+    auto& inbc = g_inputdeck.get< tag::bc >()[0].get< tag::inlet >();
+
+    // inlet velocity and material
+    auto u_in = inbc[0].get< tag::velocity >();
+    auto mat_in = inbc[0].get< tag::materialid >() - 1;
+    auto p_in = inbc[0].get< tag::pressure >();
+    auto t_in = inbc[0].get< tag::temperature >();
+
+    [[maybe_unused]] auto nsld = numSolids(nmat, solidx);
+
+    Assert( ul.size() == ncomp+nmat+3+nsld*6, "Incorrect size for appended "
+            "internal state vector" );
+
+    auto ur = ul;
+
+    // Internal cell velocity components
+    auto v1l = ul[ncomp+velocityIdx(nmat, 0)];
+    auto v2l = ul[ncomp+velocityIdx(nmat, 1)];
+    auto v3l = ul[ncomp+velocityIdx(nmat, 2)];
+
+    // External cell velocity, such that velocity = v_in at face
+    auto v1r = 2.0*u_in[0] - v1l;
+    auto v2r = 2.0*u_in[1] - v2l;
+    auto v3r = 2.0*u_in[2] - v3l;
+
+    // Normal inlet velocity
+    auto vn = u_in[0]*fn[0] + u_in[1]*fn[1] + u_in[2]*fn[2];
+
+    // Acoustic speed
+    tk::real a(0.0);
+    for (std::size_t k=0; k<nmat; ++k)
+      if (ul[volfracIdx(nmat, k)] > 1.0e-04)
+        a = std::max( a, mat_blk[k].compute< EOS::soundspeed >(
+          ul[densityIdx(nmat, k)], ul[ncomp+pressureIdx(nmat, k)],
+          ul[volfracIdx(nmat, k)], k ) );
+
+    // Mach number
+    auto Ma = vn / a;
+
+    tk::real alphamin = 1e-12;
+    tk::real pk(0.0);
+    tk::real rho(0.0);
+    for (std::size_t k=0; k<nmat; ++k) {
+      if (k == mat_in)
+        ur[volfracIdx(nmat,k)] = 1.0 -
+          (static_cast< tk::real >(nmat-1))*alphamin;
+      else
+        ur[volfracIdx(nmat,k)] = alphamin;
+
+      // Material pressure, which, for supersonic inflow, is the exterior
+      // pressure and the interior pressure for subsonic
+      if(Ma <= -1)
+        pk = p_in;
+      else
+        pk = ul[ncomp+pressureIdx(nmat,k)]/ul[volfracIdx(nmat,k)];
+      auto rhok = mat_blk[k].compute< EOS::density >(pk, t_in);
+
+      ur[ncomp+pressureIdx(nmat, k)] = ur[volfracIdx(nmat,k)] * pk;
+      ur[densityIdx(nmat,k)] = ur[volfracIdx(nmat,k)] * rhok;
+      ur[energyIdx(nmat,k)] = ur[volfracIdx(nmat,k)] *
+        mat_blk[k].compute< EOS::totalenergy >(rhok, v1r, v2r, v3r, pk);
+
+      // bulk density
+      rho += ur[densityIdx(nmat,k)];
+    }
+
+    ur[momentumIdx(nmat, 0)] = rho * v1r;
+    ur[momentumIdx(nmat, 1)] = rho * v2r;
+    ur[momentumIdx(nmat, 2)] = rho * v3r;
+
+    // velocity
+    ur[ncomp+velocityIdx(nmat, 0)] = v1r;
+    ur[ncomp+velocityIdx(nmat, 1)] = v2r;
+    ur[ncomp+velocityIdx(nmat, 2)] = v3r;
+
+    Assert( ur.size() == ncomp+nmat+3+nsld*6, "Incorrect size for appended "
+            "boundary state vector" );
+
+    return {{ std::move(ul), std::move(ur) }};
+  }
+
   //! \brief Boundary state function providing the left and right state of a
   //!   face at farfield boundaries
   //! \param[in] ncomp Number of scalar components in this PDE system

src/PDE/MultiMat/DGMultiMat.hpp

@@ -76,23 +76,25 @@ class MultiMat {
     {
       // associate boundary condition configurations with state functions
       brigand::for_each< ctr::bclist::Keys >( ConfigBC( m_bc,
-        // BC State functions
-        { dirichlet
-        , symmetry
-        , invalidBC         // Inlet BC not implemented
-        , invalidBC         // Outlet BC not implemented
-        , farfield
-        , extrapolate
-        , noslipwall },
-        // BC Gradient functions
-        { noOpGrad
-        , symmetryGrad
-        , noOpGrad
-        , noOpGrad
-        , noOpGrad
-        , noOpGrad
-        , noOpGrad }
-        ) );
+      // BC State functions
+      { dirichlet
+      , symmetry
+      , invalidBC         // Outlet BC not implemented
+      , farfield
+      , extrapolate
+      , noslipwall },
+      // BC Gradient functions
+      { noOpGrad
+      , symmetryGrad
+      , noOpGrad
+      , noOpGrad
+      , noOpGrad
+      , noOpGrad }
+      ) );
+
+      // Inlet BC has a different structure than above BCs, so it must be 
+      // handled differently than with ConfigBC
+      ConfigInletBC(m_bc, inlet, noOpGrad);
 
       // EoS initialization
       initializeMaterialEoS( m_mat_blk );

src/PDE/MultiMat/FVMultiMat.hpp

@@ -75,7 +75,6 @@ class MultiMat {
         // BC State functions
         { dirichlet
         , symmetry
-        , invalidBC         // Inlet BC not implemented
         , invalidBC         // Outlet BC not implemented
         , farfield
         , extrapolate
@@ -86,10 +85,13 @@ class MultiMat {
         , noOpGrad
         , noOpGrad
         , noOpGrad
-        , noOpGrad
         , noOpGrad }
         ) );
 
+      // Inlet BC has a different structure than above BCs, so it must be 
+      // handled differently than with ConfigBC
+      ConfigInletBC(m_bc, inlet, noOpGrad);
+
       // EoS initialization
       initializeMaterialEoS( m_mat_blk );
     }