diff --git a/src/language/modules/SchemeModule.cpp b/src/language/modules/SchemeModule.cpp
index 3376c89de8da3ae851850aab70ad4db0e7ba00f0..e1af0c003f55eff3855b1a16ae9a0c452d6ce8d4 100644
--- a/src/language/modules/SchemeModule.cpp
+++ b/src/language/modules/SchemeModule.cpp
@@ -55,6 +55,7 @@
 #include <scheme/IBoundaryConditionDescriptor.hpp>
 #include <scheme/IDiscreteFunctionDescriptor.hpp>
 #include <scheme/InflowBoundaryConditionDescriptor.hpp>
+#include <scheme/InflowListBoundaryConditionDescriptor.hpp>
 #include <scheme/NeumannBoundaryConditionDescriptor.hpp>
 #include <scheme/OutflowBoundaryConditionDescriptor.hpp>
 #include <scheme/Scalar2WavesKineticSolver.hpp>
@@ -1236,6 +1237,18 @@ SchemeModule::SchemeModule()
 
                                         ));
 
+  this->_addBuiltinFunction("inflow_list",
+                            std::function(
+
+                              [](std::shared_ptr<const IBoundaryDescriptor> boundary,
+                                 const std::vector<FunctionSymbolId>& function_id_list)
+                                -> std::shared_ptr<const IBoundaryConditionDescriptor> {
+                                return std::make_shared<InflowListBoundaryConditionDescriptor>(boundary,
+                                                                                               function_id_list);
+                              }
+
+                              ));
+
   this->_addBuiltinFunction("outflow", std::function(
 
                                          [](std::shared_ptr<const IBoundaryDescriptor> boundary)
diff --git a/src/scheme/EulerKineticSolver2D.cpp b/src/scheme/EulerKineticSolver2D.cpp
index 0f2513f67ce26de57d583909ba5c1696675c3336..d7df345b4462136b8ad14190570000d4d0229d3c 100644
--- a/src/scheme/EulerKineticSolver2D.cpp
+++ b/src/scheme/EulerKineticSolver2D.cpp
@@ -1,6 +1,7 @@
 #include <scheme/EulerKineticSolver2D.hpp>
 
 #include <language/utils/EvaluateAtPoints.hpp>
+#include <language/utils/InterpolateItemArray.hpp>
 #include <mesh/Connectivity.hpp>
 #include <mesh/Mesh.hpp>
 #include <mesh/MeshFlatNodeBoundary.hpp>
@@ -9,9 +10,9 @@
 #include <scheme/DiscreteFunctionDescriptorP0Vector.hpp>
 #include <scheme/DiscreteFunctionP0Vector.hpp>
 #include <scheme/DiscreteFunctionUtils.hpp>
-#include <scheme/FixedBoundaryConditionDescriptor.hpp>
 #include <scheme/IBoundaryConditionDescriptor.hpp>
 #include <scheme/IDiscreteFunctionDescriptor.hpp>
+#include <scheme/InflowListBoundaryConditionDescriptor.hpp>
 #include <scheme/SymmetryBoundaryConditionDescriptor.hpp>
 
 #include <analysis/GaussLegendreQuadratureDescriptor.hpp>
@@ -202,10 +203,10 @@ class EulerKineticSolver2D
   constexpr static size_t Dimension = MeshType::Dimension;
   using Rd                          = TinyVector<Dimension>;
 
-  class FixedBoundaryCondition;
   class SymmetryBoundaryCondition;
+  class InflowListBoundaryCondition;
 
-  using BoundaryCondition     = std::variant<FixedBoundaryCondition, SymmetryBoundaryCondition>;
+  using BoundaryCondition     = std::variant<SymmetryBoundaryCondition, InflowListBoundaryCondition>;
   using BoundaryConditionList = std::vector<BoundaryCondition>;
 
   const double m_dt;
@@ -247,8 +248,49 @@ class EulerKineticSolver2D
           SymmetryBoundaryCondition(getMeshFlatNodeBoundary(mesh, bc_descriptor->boundaryDescriptor())));
         break;
       }
-      case IBoundaryConditionDescriptor::Type::fixed: {
-        bc_list.emplace_back(FixedBoundaryCondition(getMeshNodeBoundary(mesh, bc_descriptor->boundaryDescriptor())));
+      case IBoundaryConditionDescriptor::Type::inflow_list: {
+        const InflowListBoundaryConditionDescriptor& inflow_list_bc_descriptor =
+          dynamic_cast<const InflowListBoundaryConditionDescriptor&>(*bc_descriptor);
+        if (inflow_list_bc_descriptor.functionSymbolIdList().size() != 2 + Dimension) {
+          std::ostringstream error_msg;
+          error_msg << "invalid number of functions for inflow boundary "
+                    << inflow_list_bc_descriptor.boundaryDescriptor() << ", found "
+                    << inflow_list_bc_descriptor.functionSymbolIdList().size() << ", expecting " << 2 + Dimension;
+          throw NormalError(error_msg.str());
+        }
+        auto node_boundary = getMeshNodeBoundary(mesh, bc_descriptor->boundaryDescriptor());
+
+        Array<int8_t> is_boundary_cell{p_mesh->numberOfCells()};
+        auto is_boundary_node = p_mesh->connectivity().isBoundaryNode();
+
+        is_boundary_cell.fill(0);
+
+        auto node_to_cell_matrix = m_mesh.connectivity().nodeToCellMatrix();
+
+        auto node_list = node_boundary.nodeList();
+        for (size_t i_node = 0; i_node < node_list.size(); ++i_node) {
+          auto cell_list = node_to_cell_matrix[node_list[i_node]];
+          for (size_t i_cell = 0; i_cell < cell_list.size(); ++i_cell) {
+            const CellId cell_id      = cell_list[i_cell];
+            is_boundary_cell[cell_id] = 1;
+          }
+        }
+        size_t nb_boundary_cells = sum(is_boundary_cell);
+        Array<CellId> cell_list{nb_boundary_cells};
+
+        size_t index = 0;
+        for (CellId cell_id = 0; cell_id < p_mesh->numberOfCells(); ++cell_id) {
+          if (is_boundary_cell[cell_id]) {
+            cell_list[index++] = cell_id;
+          }
+        }
+        auto xj = MeshDataManager::instance().getMeshData(m_mesh).xj();
+
+        Table<const double> cell_array_list =
+          InterpolateItemArray<double(Rd)>::template interpolate<ItemType::cell>(inflow_list_bc_descriptor
+                                                                                   .functionSymbolIdList(),
+                                                                                 xj, cell_list);
+        bc_list.emplace_back(InflowListBoundaryCondition(cell_array_list, cell_list));
         break;
       }
       default: {
@@ -530,14 +572,11 @@ class EulerKineticSolver2D
       std::visit(
         [=](auto&& bc) {
           using BCType = std::decay_t<decltype(bc)>;
-          if constexpr (std::is_same_v<BCType, FixedBoundaryCondition>) {
-            auto node_list = bc.nodeList();
-            for (size_t i_node = 0; i_node < node_list.size(); ++i_node) {
-              auto cell_list = node_to_cell_matrix[node_list[i_node]];
-              for (size_t i_cell = 0; i_cell < cell_list.size(); ++i_cell) {
-                const CellId cell_id      = cell_list[i_cell];
-                is_boundary_cell[cell_id] = true;
-              }
+          if constexpr (std::is_same_v<BCType, InflowListBoundaryCondition>) {
+            auto cell_list = bc.cellList();
+            for (size_t i_cell = 0; i_cell < cell_list.size(); ++i_cell) {
+              const CellId cell_id      = cell_list[i_cell];
+              is_boundary_cell[cell_id] = true;
             }
           }
         },
@@ -559,16 +598,36 @@ class EulerKineticSolver2D
 
     CellValue<bool> is_boundary_cell = getBoundaryCells(bc_list);
 
+    DiscreteFunctionP0<double> rho_ex{p_mesh};
     DiscreteFunctionP0<double> rho_u1_ex{p_mesh};
     DiscreteFunctionP0<double> rho_u2_ex{p_mesh};
+    DiscreteFunctionP0<double> rho_E_ex{p_mesh};
 
-    parallel_for(
-      m_mesh.numberOfCells(), PUGS_LAMBDA(const CellId cell_id) {
-        rho_u1_ex[cell_id] = m_rho_u_ex[cell_id][0];
-        rho_u2_ex[cell_id] = m_rho_u_ex[cell_id][1];
-      });
+    rho_ex.fill(1);
+    rho_u1_ex.fill(0);
+    rho_u2_ex.fill(0);
+    rho_E_ex.fill(0);
+
+    for (auto&& bc_v : bc_list) {
+      std::visit(
+        [=](auto&& bc) {
+          using BCType = std::decay_t<decltype(bc)>;
+          if constexpr (std::is_same_v<BCType, InflowListBoundaryCondition>) {
+            auto cell_list       = bc.cellList();
+            auto cell_array_list = bc.cellArrayList();
+            for (size_t i_cell = 0; i_cell < cell_list.size(); ++i_cell) {
+              const CellId cell_id = cell_list[i_cell];
+              rho_ex[cell_id]      = cell_array_list[i_cell][0];
+              rho_u1_ex[cell_id]   = cell_array_list[i_cell][1];
+              rho_u2_ex[cell_id]   = cell_array_list[i_cell][2];
+              rho_E_ex[cell_id]    = cell_array_list[i_cell][3];
+            }
+          }
+        },
+        bc_v);
+    }
 
-    const CellArray<const TinyVector<Dimension>> APFn_ex = getA(m_rho_ex, rho_u1_ex, rho_u2_ex, m_rho_E_ex);
+    const CellArray<const TinyVector<Dimension>> APFn_ex = getA(rho_ex, rho_u1_ex, rho_u2_ex, rho_E_ex);
     const DiscreteFunctionP0<TinyVector<Dimension>> A_rho_ex{p_mesh};
     const DiscreteFunctionP0<TinyVector<Dimension>> A_rho_u1_ex{p_mesh};
     const DiscreteFunctionP0<TinyVector<Dimension>> A_rho_u2_ex{p_mesh};
@@ -582,10 +641,10 @@ class EulerKineticSolver2D
         A_rho_E_ex[cell_id]  = APFn_ex[cell_id][3];
       });
 
-    DiscreteFunctionP0Vector<double> Fn_rho_ex    = compute_M(m_rho_ex, A_rho_ex);
+    DiscreteFunctionP0Vector<double> Fn_rho_ex    = compute_M(rho_ex, A_rho_ex);
     DiscreteFunctionP0Vector<double> Fn_rho_u1_ex = compute_M(rho_u1_ex, A_rho_u1_ex);
     DiscreteFunctionP0Vector<double> Fn_rho_u2_ex = compute_M(rho_u2_ex, A_rho_u2_ex);
-    DiscreteFunctionP0Vector<double> Fn_rho_E_ex  = compute_M(m_rho_E_ex, A_rho_E_ex);
+    DiscreteFunctionP0Vector<double> Fn_rho_E_ex  = compute_M(rho_E_ex, A_rho_E_ex);
 
     DiscreteFunctionP0Vector<double> Fn_rho    = copy(m_Fn_rho);
     DiscreteFunctionP0Vector<double> Fn_rho_u1 = copy(m_Fn_rho_u1);
@@ -797,24 +856,6 @@ eulerKineticSolver2D(const double& dt,
     mesh_v->variant());
 }
 
-template <MeshConcept MeshType>
-class EulerKineticSolver2D<MeshType>::FixedBoundaryCondition
-{
- private:
-  const MeshNodeBoundary m_mesh_node_boundary;
-
- public:
-  const Array<const NodeId>&
-  nodeList() const
-  {
-    return m_mesh_node_boundary.nodeList();
-  }
-
-  FixedBoundaryCondition(const MeshNodeBoundary& mesh_node_boundary) : m_mesh_node_boundary{mesh_node_boundary} {}
-
-  ~FixedBoundaryCondition() = default;
-};
-
 template <MeshConcept MeshType>
 class EulerKineticSolver2D<MeshType>::SymmetryBoundaryCondition
 {
@@ -853,3 +894,43 @@ class EulerKineticSolver2D<MeshType>::SymmetryBoundaryCondition
 
   ~SymmetryBoundaryCondition() = default;
 };
+
+template <MeshConcept MeshType>
+class EulerKineticSolver2D<MeshType>::InflowListBoundaryCondition
+{
+ public:
+  static constexpr const size_t Dimension = MeshType::Dimension;
+
+  using Rd = TinyVector<Dimension, double>;
+
+ private:
+  const Table<const double> m_cell_array_list;
+  const Array<const CellId> m_cell_list;
+
+ public:
+  size_t
+  numberOfCells() const
+  {
+    return m_cell_list.size();
+  }
+
+  const Array<const CellId>&
+  cellList() const
+  {
+    return m_cell_list;
+  }
+
+  const Table<const double>&
+  cellArrayList() const
+  {
+    return m_cell_array_list;
+  }
+
+  InflowListBoundaryCondition(const Table<const double>& cell_array_list, const Array<const CellId> cell_list)
+    : m_cell_array_list(cell_array_list), m_cell_list(cell_list)
+  {
+    ;
+  }
+
+  ~InflowListBoundaryCondition() = default;
+};
diff --git a/src/scheme/IBoundaryConditionDescriptor.hpp b/src/scheme/IBoundaryConditionDescriptor.hpp
index 58e8b6076ff2c363049de6f107e6a949966d6c72..14967ae3f26778823565422cb722112d6476e23e 100644
--- a/src/scheme/IBoundaryConditionDescriptor.hpp
+++ b/src/scheme/IBoundaryConditionDescriptor.hpp
@@ -18,6 +18,7 @@ class IBoundaryConditionDescriptor
     fixed,
     free,
     inflow,
+    inflow_list,
     neumann,
     outflow,
     symmetry
diff --git a/src/scheme/InflowListBoundaryConditionDescriptor.hpp b/src/scheme/InflowListBoundaryConditionDescriptor.hpp
new file mode 100644
index 0000000000000000000000000000000000000000..ad1427cf72532b770a378c5e86ecbc1cc39ac063
--- /dev/null
+++ b/src/scheme/InflowListBoundaryConditionDescriptor.hpp
@@ -0,0 +1,48 @@
+#ifndef INFLOW_LIST_BOUNDARY_CONDITION_DESCRIPTOR_HPP
+#define INFLOW_LIST_BOUNDARY_CONDITION_DESCRIPTOR_HPP
+
+#include <language/utils/FunctionSymbolId.hpp>
+#include <mesh/IBoundaryDescriptor.hpp>
+#include <scheme/BoundaryConditionDescriptorBase.hpp>
+
+#include <memory>
+
+class InflowListBoundaryConditionDescriptor : public BoundaryConditionDescriptorBase
+{
+ private:
+  std::ostream&
+  _write(std::ostream& os) const final
+  {
+    os << "inflow_list(" << ',' << *m_boundary_descriptor << ")";
+    return os;
+  }
+
+  const std::vector<FunctionSymbolId> m_function_symbol_id_list;
+
+ public:
+  const std::vector<FunctionSymbolId>&
+  functionSymbolIdList() const
+  {
+    return m_function_symbol_id_list;
+  }
+
+  Type
+  type() const final
+  {
+    return Type::inflow_list;
+  }
+
+  InflowListBoundaryConditionDescriptor(std::shared_ptr<const IBoundaryDescriptor> boundary_descriptor,
+                                        const std::vector<FunctionSymbolId>& function_symbol_id_list)
+    : BoundaryConditionDescriptorBase{boundary_descriptor}, m_function_symbol_id_list{function_symbol_id_list}
+  {
+    ;
+  }
+
+  InflowListBoundaryConditionDescriptor(const InflowListBoundaryConditionDescriptor&) = delete;
+  InflowListBoundaryConditionDescriptor(InflowListBoundaryConditionDescriptor&&)      = delete;
+
+  ~InflowListBoundaryConditionDescriptor() = default;
+};
+
+#endif   // INFLOW_LIST_BOUNDARY_CONDITION_DESCRIPTOR_HPP