diff --git a/src/scheme/EulerKineticSolverMultiD.cpp b/src/scheme/EulerKineticSolverMultiD.cpp
index 20077749b71094bed6ebf343795b264fdd3e42c6..0f9acb6513f58837c480dbd2e5271fc42c5de64a 100644
--- a/src/scheme/EulerKineticSolverMultiD.cpp
+++ b/src/scheme/EulerKineticSolverMultiD.cpp
@@ -535,7 +535,7 @@ class EulerKineticSolverMultiD
                 const unsigned int i_node = node_local_number_in_its_cell[i_cell];
 
                 double li_njr = dot(m_lambda_vector[i_wave], njr(cell_id, i_node));
-                if (li_njr > 0) {
+                if (li_njr > 1e-14) {
                   Fr[node_id][i_wave] += li_njr * Fn[cell_id][i_wave];
                   sum_li_njr += li_njr;
                 }
@@ -874,6 +874,7 @@ class EulerKineticSolverMultiD
               nr[node_id] = 1. / sqrt(dot(nr[node_id], nr[node_id])) * nr[node_id];
               auto nxn    = tensorProduct(nr[node_id], nr[node_id]);
               auto txt    = I - nxn;
+              std::cout << node_id << ", " << nr[node_id] << "\n";
 
               for (size_t i_wave = 0; i_wave < nb_waves; ++i_wave) {
                 double sum                = 0;
diff --git a/src/scheme/EulerKineticSolverMultiDOrder2.cpp b/src/scheme/EulerKineticSolverMultiDOrder2.cpp
index acefc2734a3475e14ca6a3269a6e9c22fe795a8f..030a0bab331a41ac32dc85a479fa3af32f4e34b6 100644
--- a/src/scheme/EulerKineticSolverMultiDOrder2.cpp
+++ b/src/scheme/EulerKineticSolverMultiDOrder2.cpp
@@ -343,111 +343,90 @@ class EulerKineticSolverMultiDOrder2
     return PFn;
   }
 
-  template <typename T>
-  NodeArray<T>
-  compute_Flux1(CellArray<T>& Fn)
+  std::tuple<FaceArray<double>, FaceArray<TinyVector<Dimension>>, FaceArray<double>>
+  compute_Flux_face_equilibrium(const DiscreteFunctionDPk<Dimension, const double>& DPk_rho,
+                                const DiscreteFunctionDPk<Dimension, const TinyVector<Dimension>>& DPk_rho_u,
+                                const DiscreteFunctionDPk<Dimension, const double>& DPk_rho_E,
+                                const CellArray<const double>& Fn_rho,
+                                const CellArray<const TinyVector<Dimension>>& Fn_rho_u,
+                                const CellArray<const double>& Fn_rho_E,
+                                FaceValue<bool>& is_wall_boundary_face)
   {
-    if constexpr (Dimension > 1) {
-      const NodeValuePerCell<const TinyVector<Dimension>> njr = MeshDataManager::instance().getMeshData(m_mesh).njr();
-      const size_t nb_waves                                   = m_lambda_vector.size();
-      NodeArray<T> Fr(m_mesh.connectivity(), nb_waves);
-      const auto& node_local_numbers_in_their_cells = p_mesh->connectivity().nodeLocalNumbersInTheirCells();
-      const auto& node_to_cell_matrix               = p_mesh->connectivity().nodeToCellMatrix();
-
-      if constexpr (is_tiny_vector_v<T>) {
-        Fr.fill(zero);
-      } else {
-        Fr.fill(0.);
-      }
-
-      parallel_for(
-        p_mesh->numberOfNodes(), PUGS_LAMBDA(NodeId node_id) {
-          const auto& node_local_number_in_its_cell = node_local_numbers_in_their_cells.itemArray(node_id);
-          const auto& node_to_cell                  = node_to_cell_matrix[node_id];
-          for (size_t i_wave = 0; i_wave < nb_waves; ++i_wave) {
-            double sum_li_njr = 0;
-
-            for (size_t i_cell = 0; i_cell < node_to_cell.size(); ++i_cell) {
-              const CellId cell_id      = node_to_cell[i_cell];
-              const unsigned int i_node = node_local_number_in_its_cell[i_cell];
+    const size_t nb_waves = m_lambda_vector.size();
+    FaceArray<double> Fl_rho(m_mesh.connectivity(), nb_waves);
+    FaceArray<TinyVector<Dimension>> Fl_rho_u(m_mesh.connectivity(), nb_waves);
+    FaceArray<double> Fl_rho_E(m_mesh.connectivity(), nb_waves);
+    const auto& face_local_numbers_in_their_cells = p_mesh->connectivity().faceLocalNumbersInTheirCells();
+    const auto& face_to_cell_matrix               = p_mesh->connectivity().faceToCellMatrix();
+    const auto& face_cell_is_reversed             = p_mesh->connectivity().cellFaceIsReversed();
+    TinyMatrix<Dimension> I                       = identity;
 
-              double li_njr = dot(m_lambda_vector[i_wave], njr(cell_id, i_node));
-              if (li_njr > 0) {
-                Fr[node_id][i_wave] += Fn[cell_id][i_wave] * li_njr;
-                sum_li_njr += li_njr;
-              }
-            }
-            if (sum_li_njr != 0) {
-              Fr[node_id][i_wave] /= sum_li_njr;
-            }
-          }
-        });
+    Fl_rho.fill(0.);
+    Fl_rho_u.fill(zero);
+    Fl_rho_E.fill(0.);
 
-      return Fr;
-    } else {
-      throw NormalError("Glace not defined in 1d");
+    TinyVector<MeshType::Dimension> inv_S = zero;
+    for (size_t d = 0; d < MeshType::Dimension; ++d) {
+      for (size_t i = 0; i < nb_waves; ++i) {
+        inv_S[d] += m_lambda_vector[i][d] * m_lambda_vector[i][d];
+      }
+      inv_S[d] = 1. / inv_S[d];
     }
-  }
 
-  template <typename T>
-  FaceArrayPerNode<T>
-  compute_Flux1_eucclhyd(CellArray<T> Fn)
-  {
-    if constexpr (Dimension > 1) {
-      const NodeValuePerFace<const TinyVector<Dimension>> nlr = MeshDataManager::instance().getMeshData(m_mesh).nlr();
-      const size_t nb_waves                                   = m_lambda_vector.size();
-      FaceArrayPerNode<T> Flr(m_mesh.connectivity(), nb_waves);
-      const auto& node_local_numbers_in_their_faces = p_mesh->connectivity().nodeLocalNumbersInTheirFaces();
-      const auto& face_local_numbers_in_their_cells = p_mesh->connectivity().faceLocalNumbersInTheirCells();
-      const auto& node_to_face_matrix               = p_mesh->connectivity().nodeToFaceMatrix();
-      const auto& face_to_cell_matrix               = p_mesh->connectivity().faceToCellMatrix();
-      const auto& face_cell_is_reversed             = p_mesh->connectivity().cellFaceIsReversed();
+    parallel_for(
+      p_mesh->numberOfFaces(), PUGS_LAMBDA(FaceId face_id) {
+        for (size_t i_wave = 0; i_wave < nb_waves; ++i_wave) {
+          const auto& face_local_number_in_its_cell = face_local_numbers_in_their_cells.itemArray(face_id);
+          const auto& face_to_cell                  = face_to_cell_matrix[face_id];
 
-      if constexpr (is_tiny_vector_v<T>) {
-        Flr.fill(zero);
-      } else {
-        Flr.fill(0.);
-      }
+          for (size_t i_cell = 0; i_cell < face_to_cell.size(); ++i_cell) {
+            const CellId cell_id     = face_to_cell[i_cell];
+            const size_t i_face_cell = face_local_number_in_its_cell[i_cell];
 
-      parallel_for(
-        p_mesh->numberOfNodes(), PUGS_LAMBDA(NodeId node_id) {
-          const auto& node_local_number_in_its_face = node_local_numbers_in_their_faces.itemArray(node_id);
-          const auto& node_to_face                  = node_to_face_matrix[node_id];
+            const double rho_jl                  = DPk_rho[cell_id](m_xl[face_id]);
+            const TinyVector<Dimension> rho_u_jl = DPk_rho_u[cell_id](m_xl[face_id]);
+            const double rho_E_jl                = DPk_rho_E[cell_id](m_xl[face_id]);
 
-          for (size_t i_wave = 0; i_wave < nb_waves; ++i_wave) {
-            for (size_t i_face = 0; i_face < node_to_face.size(); ++i_face) {
-              const FaceId face_id                      = node_to_face[i_face];
-              const auto& face_local_number_in_its_cell = face_local_numbers_in_their_cells.itemArray(face_id);
-              const auto& face_to_cell                  = face_to_cell_matrix[face_id];
-              const size_t i_node_face                  = node_local_number_in_its_face[i_face];
+            const double rho_e_jl = rho_E_jl - 0.5 * (1. / rho_jl) * dot(rho_u_jl, rho_u_jl);
+            const double p_jl     = (m_gamma - 1) * rho_e_jl;
 
-              double sum = 0;
+            const TinyVector<Dimension> A_rho_jl   = rho_u_jl;
+            const TinyMatrix<Dimension> A_rho_u_jl = 1. / rho_jl * tensorProduct(rho_u_jl, rho_u_jl) + p_jl * I;
+            const TinyVector<Dimension> A_rho_E_jl = (rho_E_jl + p_jl) / rho_jl * rho_u_jl;
 
-              for (size_t i_cell = 0; i_cell < face_to_cell.size(); ++i_cell) {
-                const CellId cell_id     = face_to_cell[i_cell];
-                const size_t i_face_cell = face_local_number_in_its_cell[i_cell];
+            double Fn_rho_jl                  = rho_jl / nb_waves;
+            TinyVector<Dimension> Fn_rho_u_jl = 1. / nb_waves * rho_u_jl;
+            double Fn_rho_E_jl                = rho_E_jl / nb_waves;
 
-                TinyVector<Dimension> n_face = nlr(face_id, i_node_face);
-                const double sign            = face_cell_is_reversed(cell_id, i_face_cell) ? -1 : 1;
+            for (size_t d1 = 0; d1 < Dimension; ++d1) {
+              Fn_rho_jl += inv_S[d1] * m_lambda_vector[i_wave][d1] * A_rho_jl[d1];
+              for (size_t d2 = 0; d2 < Dimension; ++d2) {
+                Fn_rho_u_jl[d1] += inv_S[d2] * m_lambda_vector[i_wave][d2] * A_rho_u_jl(d2, d1);
+              }
+              Fn_rho_E_jl += inv_S[d1] * m_lambda_vector[i_wave][d1] * A_rho_E_jl[d1];
+            }
 
-                double li_nlr = sign * dot(m_lambda_vector[i_wave], n_face);
+            TinyVector<Dimension> n_face = nl[face_id];
+            const double sign            = face_cell_is_reversed(cell_id, i_face_cell) ? -1 : 1;
 
-                if (li_nlr > 1e-14) {
-                  Flr[node_id][i_face][i_wave] += Fn[cell_id][i_wave] * li_nlr;
-                  sum += li_nlr;
-                }
-              }
-              if (sum != 0) {
-                Flr[node_id][i_face][i_wave] /= sum;
+            double li_nl = sign * dot(m_lambda_vector[i_wave], n_face);
+
+            if (li_nl > 1e-14) {
+              if ((not is_wall_boundary_face[face_id])) {
+                Fl_rho[face_id][i_wave] += Fn_rho_jl;
+                Fl_rho_u[face_id][i_wave] += Fn_rho_u_jl;
+                Fl_rho_E[face_id][i_wave] += Fn_rho_E_jl;
+              } else {
+                Fl_rho[face_id][i_wave] += Fn_rho[cell_id][i_wave];
+                Fl_rho_u[face_id][i_wave] += Fn_rho_u[cell_id][i_wave];
+                Fl_rho_E[face_id][i_wave] += Fn_rho_E[cell_id][i_wave];
               }
             }
           }
-        });
+        }
+      });
 
-      return Flr;
-    } else {
-      throw NormalError("Eucclhyd not defined in 1d");
-    }
+    return std::make_tuple(Fl_rho, Fl_rho_u, Fl_rho_E);
   }
 
   template <typename T>
@@ -803,6 +782,211 @@ class EulerKineticSolverMultiDOrder2
     }
   }
 
+  void
+  apply_face_bc_equilibrium(FaceArray<double> Fl_rho,
+                            FaceArray<TinyVector<Dimension>> Fl_rho_u,
+                            FaceArray<double> Fl_rho_E,
+                            DiscreteFunctionDPk<Dimension, double> rho,
+                            DiscreteFunctionDPk<Dimension, TinyVector<Dimension>> rho_u,
+                            DiscreteFunctionDPk<Dimension, double> rho_E,
+                            const CellArray<const double> Fj_rho,
+                            const CellArray<const TinyVector<Dimension>> Fj_rho_u,
+                            const CellArray<const double> Fj_rho_E,
+                            const BoundaryConditionList& bc_list)
+  {
+    const size_t nb_waves                         = m_lambda_vector.size();
+    const auto& face_local_numbers_in_their_cells = p_mesh->connectivity().faceLocalNumbersInTheirCells();
+    const auto& face_to_cell_matrix               = p_mesh->connectivity().faceToCellMatrix();
+    const auto& face_cell_is_reversed             = p_mesh->connectivity().cellFaceIsReversed();
+
+    TinyVector<MeshType::Dimension> inv_S = zero;
+    for (size_t d = 0; d < MeshType::Dimension; ++d) {
+      for (size_t i = 0; i < nb_waves; ++i) {
+        inv_S[d] += m_lambda_vector[i][d] * m_lambda_vector[i][d];
+      }
+      inv_S[d] = 1. / inv_S[d];
+    }
+
+    for (auto&& bc_v : bc_list) {
+      std::visit(
+        [=, this](auto&& bc) {
+          using BCType = std::decay_t<decltype(bc)>;
+          if constexpr (std::is_same_v<BCType, SymmetryBoundaryCondition>) {
+            auto face_list          = bc.faceList();
+            auto n                  = bc.outgoingNormal();
+            auto nxn                = tensorProduct(n, n);
+            TinyMatrix<Dimension> I = identity;
+            auto txt                = I - nxn;
+
+            for (size_t i_face = 0; i_face < face_list.size(); ++i_face) {
+              const FaceId face_id = face_list[i_face];
+
+              for (size_t i_wave = 0; i_wave < nb_waves; ++i_wave) {
+                const auto& face_local_number_in_its_cell = face_local_numbers_in_their_cells.itemArray(face_id);
+                const auto& face_to_cell                  = face_to_cell_matrix[face_id];
+
+                for (size_t i_cell = 0; i_cell < face_to_cell.size(); ++i_cell) {
+                  const CellId cell_id     = face_to_cell[i_cell];
+                  const size_t i_face_cell = face_local_number_in_its_cell[i_cell];
+
+                  double rhol                  = rho[cell_id](m_xl[face_id]);
+                  TinyVector<Dimension> rho_ul = rho_u[cell_id](m_xl[face_id]);
+                  double rho_El                = rho_E[cell_id](m_xl[face_id]);
+
+                  TinyVector<Dimension> n_face = nl[face_id];
+                  const double sign            = face_cell_is_reversed(cell_id, i_face_cell) ? -1 : 1;
+
+                  double li_nl = sign * dot(m_lambda_vector[i_wave], n_face);
+                  if (li_nl < -1e-14) {
+                    double rhol_prime                  = rhol;
+                    TinyVector<Dimension> rho_ul_prime = txt * rho_ul - nxn * rho_ul;
+                    double rho_El_prime                = rho_El;
+
+                    double rho_el = rho_El - 0.5 * (1. / rhol) * dot(rho_ul, rho_ul);
+                    double pl     = (m_gamma - 1) * rho_el;
+
+                    TinyVector<Dimension> A_rho_prime = rho_ul_prime;
+                    TinyMatrix<Dimension> A_rho_u_prime =
+                      1. / rhol_prime * tensorProduct(rho_ul_prime, rho_ul_prime) + pl * I;
+                    TinyVector<Dimension> A_rho_E_prime = (rho_El_prime + pl) / rhol_prime * rho_ul_prime;
+
+                    double Fn_rho_prime                  = rhol_prime / nb_waves;
+                    TinyVector<Dimension> Fn_rho_u_prime = 1. / nb_waves * rho_ul_prime;
+                    double Fn_rho_E_prime                = rho_El_prime / nb_waves;
+
+                    for (size_t d1 = 0; d1 < Dimension; ++d1) {
+                      Fn_rho_prime += inv_S[d1] * m_lambda_vector[i_wave][d1] * A_rho_prime[d1];
+                      for (size_t d2 = 0; d2 < Dimension; ++d2) {
+                        Fn_rho_u_prime[d1] += inv_S[d2] * m_lambda_vector[i_wave][d2] * A_rho_u_prime(d2, d1);
+                      }
+                      Fn_rho_E_prime += inv_S[d1] * m_lambda_vector[i_wave][d1] * A_rho_E_prime[d1];
+                    }
+                    Fl_rho[face_id][i_wave] += Fn_rho_prime;
+                    Fl_rho_u[face_id][i_wave] += Fn_rho_u_prime;
+                    Fl_rho_E[face_id][i_wave] += Fn_rho_E_prime;
+                  }
+                }
+              }
+            }
+          } else if constexpr (std::is_same_v<BCType, WallBoundaryCondition>) {
+            auto face_list = bc.faceList();
+
+            for (size_t i_face = 0; i_face < face_list.size(); ++i_face) {
+              const FaceId face_id = face_list[i_face];
+
+              for (size_t i_wave = 0; i_wave < nb_waves; ++i_wave) {
+                const auto& face_local_number_in_its_cell = face_local_numbers_in_their_cells.itemArray(face_id);
+                const auto& face_to_cell                  = face_to_cell_matrix[face_id];
+
+                for (size_t i_cell = 0; i_cell < face_to_cell.size(); ++i_cell) {
+                  const CellId cell_id     = face_to_cell[i_cell];
+                  const size_t i_face_cell = face_local_number_in_its_cell[i_cell];
+
+                  double rhol                  = 0;
+                  TinyVector<Dimension> rho_ul = zero;
+                  double rho_El                = 0;
+
+                  for (size_t k_wave = 0; k_wave < nb_waves; ++k_wave) {
+                    rhol += Fj_rho[cell_id][k_wave];
+                    rho_ul += Fj_rho_u[cell_id][k_wave];
+                    rho_El += Fj_rho_E[cell_id][k_wave];
+                  }
+
+                  const double sign       = face_cell_is_reversed(cell_id, i_face_cell) ? -1 : 1;
+                  auto n                  = sign * nl[face_id];
+                  auto nxn                = tensorProduct(n, n);
+                  TinyMatrix<Dimension> I = identity;
+                  auto txt                = I - nxn;
+
+                  double li_nl = dot(m_lambda_vector[i_wave], n);
+                  if (li_nl < -1e-14) {
+                    double rhol_prime                  = rhol;
+                    TinyVector<Dimension> rho_ul_prime = txt * rho_ul - nxn * rho_ul;
+                    double rho_El_prime                = rho_El;
+
+                    double rho_el = rho_El - 0.5 * (1. / rhol) * dot(rho_ul, rho_ul);
+                    double pl     = (m_gamma - 1) * rho_el;
+
+                    TinyVector<Dimension> A_rho_prime = rho_ul_prime;
+                    TinyMatrix<Dimension> A_rho_u_prime =
+                      1. / rhol_prime * tensorProduct(rho_ul_prime, rho_ul_prime) + pl * I;
+                    TinyVector<Dimension> A_rho_E_prime = (rho_El_prime + pl) / rhol_prime * rho_ul_prime;
+
+                    double Fn_rho_prime                  = rhol_prime / nb_waves;
+                    TinyVector<Dimension> Fn_rho_u_prime = 1. / nb_waves * rho_ul_prime;
+                    double Fn_rho_E_prime                = rho_El_prime / nb_waves;
+
+                    for (size_t d1 = 0; d1 < Dimension; ++d1) {
+                      Fn_rho_prime += inv_S[d1] * m_lambda_vector[i_wave][d1] * A_rho_prime[d1];
+                      for (size_t d2 = 0; d2 < Dimension; ++d2) {
+                        Fn_rho_u_prime[d1] += inv_S[d2] * m_lambda_vector[i_wave][d2] * A_rho_u_prime(d2, d1);
+                      }
+                      Fn_rho_E_prime += inv_S[d1] * m_lambda_vector[i_wave][d1] * A_rho_E_prime[d1];
+                    }
+
+                    Fl_rho[face_id][i_wave] += Fn_rho_prime;
+                    Fl_rho_u[face_id][i_wave] += Fn_rho_u_prime;
+                    Fl_rho_E[face_id][i_wave] += Fn_rho_E_prime;
+                  }
+                }
+              }
+            }
+          } else if constexpr (std::is_same_v<BCType, OutflowBoundaryCondition>) {
+            auto face_list          = bc.faceList();
+            TinyMatrix<Dimension> I = identity;
+
+            for (size_t i_face = 0; i_face < face_list.size(); ++i_face) {
+              const FaceId face_id = face_list[i_face];
+
+              for (size_t i_wave = 0; i_wave < nb_waves; ++i_wave) {
+                const auto& face_local_number_in_its_cell = face_local_numbers_in_their_cells.itemArray(face_id);
+                const auto& face_to_cell                  = face_to_cell_matrix[face_id];
+
+                for (size_t i_cell = 0; i_cell < face_to_cell.size(); ++i_cell) {
+                  const CellId cell_id     = face_to_cell[i_cell];
+                  const size_t i_face_cell = face_local_number_in_its_cell[i_cell];
+
+                  TinyVector<Dimension> n_face = nl[face_id];
+                  const double sign            = face_cell_is_reversed(cell_id, i_face_cell) ? -1 : 1;
+
+                  double li_nl = sign * dot(m_lambda_vector[i_wave], n_face);
+                  if (li_nl < -1e-14) {
+                    double rhol                  = rho[cell_id](m_xl[face_id]);
+                    TinyVector<Dimension> rho_ul = rho_u[cell_id](m_xl[face_id]);
+                    double rho_El                = rho_E[cell_id](m_xl[face_id]);
+
+                    double rho_el = rho_El - 0.5 * (1. / rhol) * dot(rho_ul, rho_ul);
+                    double pl     = (m_gamma - 1) * rho_el;
+
+                    TinyVector<Dimension> A_rho   = rho_ul;
+                    TinyMatrix<Dimension> A_rho_u = 1. / rhol * tensorProduct(rho_ul, rho_ul) + pl * I;
+                    TinyVector<Dimension> A_rho_E = (rho_El + pl) / rhol * rho_ul;
+
+                    double Fn_rho                  = rhol / nb_waves;
+                    TinyVector<Dimension> Fn_rho_u = 1. / nb_waves * rho_ul;
+                    double Fn_rho_E                = rho_El / nb_waves;
+
+                    for (size_t d1 = 0; d1 < Dimension; ++d1) {
+                      Fn_rho += inv_S[d1] * m_lambda_vector[i_wave][d1] * A_rho[d1];
+                      for (size_t d2 = 0; d2 < Dimension; ++d2) {
+                        Fn_rho_u[d1] += inv_S[d2] * m_lambda_vector[i_wave][d2] * A_rho_u(d2, d1);
+                      }
+                      Fn_rho_E += inv_S[d1] * m_lambda_vector[i_wave][d1] * A_rho_E[d1];
+                    }
+
+                    Fl_rho[face_id][i_wave] += Fn_rho;
+                    Fl_rho_u[face_id][i_wave] += Fn_rho_u;
+                    Fl_rho_E[face_id][i_wave] += Fn_rho_E;
+                  }
+                }
+              }
+            }
+          }
+        },
+        bc_v);
+    }
+  }
+
   void
   apply_glace_bc(NodeArray<double> Fr_rho,
                  NodeArray<TinyVector<Dimension>> Fr_rho_u,
@@ -1520,7 +1704,8 @@ class EulerKineticSolverMultiDOrder2
   }
 
   void
-  scalar_limiter(const CellArray<const double>& f, DiscreteFunctionDPkVector<Dimension, double>& DPk_fh) const
+  scalar_limiter_distributions(const CellArray<const double>& f,
+                               DiscreteFunctionDPkVector<Dimension, double>& DPk_fh) const
   {
     StencilManager::BoundaryDescriptorList symmetry_boundary_descriptor_list;
     StencilDescriptor stencil_descriptor{1, StencilDescriptor::ConnectionType::by_nodes};
@@ -1569,8 +1754,8 @@ class EulerKineticSolverMultiDOrder2
   }
 
   void
-  vector_limiter(const CellArray<const TinyVector<Dimension>>& f,
-                 DiscreteFunctionDPkVector<Dimension, TinyVector<Dimension>>& DPk_fh) const
+  vector_limiter_distributions(const CellArray<const TinyVector<Dimension>>& f,
+                               DiscreteFunctionDPkVector<Dimension, TinyVector<Dimension>>& DPk_fh) const
   {
     StencilManager::BoundaryDescriptorList symmetry_boundary_descriptor_list;
     StencilDescriptor stencil_descriptor{1, StencilDescriptor::ConnectionType::by_nodes};
@@ -1624,6 +1809,105 @@ class EulerKineticSolverMultiDOrder2
       });
   }
 
+  void
+  scalar_limiter_equilibrium(const CellValue<const double>& u, DiscreteFunctionDPk<Dimension, double>& DPk_uh) const
+  {
+    StencilManager::BoundaryDescriptorList symmetry_boundary_descriptor_list;
+    StencilDescriptor stencil_descriptor{1, StencilDescriptor::ConnectionType::by_nodes};
+    auto stencil = StencilManager::instance().getCellToCellStencilArray(m_mesh.connectivity(), stencil_descriptor,
+                                                                        symmetry_boundary_descriptor_list);
+
+    parallel_for(
+      m_mesh.numberOfCells(), PUGS_LAMBDA(const CellId cell_id) {
+        const double uj = u[cell_id];
+
+        double u_min = uj;
+        double u_max = uj;
+
+        const auto cell_stencil = stencil[cell_id];
+        for (size_t i_cell = 0; i_cell < cell_stencil.size(); ++i_cell) {
+          u_min = std::min(u_min, u[cell_stencil[i_cell]]);
+          u_max = std::max(u_max, u[cell_stencil[i_cell]]);
+        }
+
+        const double eps = 1E-14;
+        double coef      = 1;
+        double lambda    = 1.;
+
+        for (size_t i_cell = 0; i_cell < cell_stencil.size(); ++i_cell) {
+          const CellId cell_i_id = cell_stencil[i_cell];
+          const double u_bar_i   = DPk_uh[cell_id](m_xj[cell_i_id]);
+          const double Delta     = (u_bar_i - uj);
+
+          if (Delta > eps) {
+            coef = std::min(1., (u_max - uj) / Delta);
+          } else if (Delta < -eps) {
+            coef = std::min(1., (u_min - uj) / Delta);
+          }
+          lambda = std::min(lambda, coef);
+        }
+
+        auto coefficients = DPk_uh.coefficients(cell_id);
+        coefficients[0]   = (1 - lambda) * u[cell_id] + lambda * coefficients[0];
+        for (size_t i = 1; i < coefficients.size(); ++i) {
+          coefficients[i] *= lambda;
+        }
+      });
+  }
+
+  void
+  vector_limiter_equilibrium(const CellValue<const TinyVector<Dimension>>& u,
+                             DiscreteFunctionDPk<Dimension, TinyVector<Dimension>>& DPk_uh) const
+  {
+    StencilManager::BoundaryDescriptorList symmetry_boundary_descriptor_list;
+    StencilDescriptor stencil_descriptor{1, StencilDescriptor::ConnectionType::by_nodes};
+    auto stencil = StencilManager::instance().getCellToCellStencilArray(m_mesh.connectivity(), stencil_descriptor,
+                                                                        symmetry_boundary_descriptor_list);
+
+    parallel_for(
+      m_mesh.numberOfCells(), PUGS_LAMBDA(const CellId cell_id) {
+        const TinyVector<Dimension> uj = u[cell_id];
+
+        TinyVector<Dimension> u_min = uj;
+        TinyVector<Dimension> u_max = uj;
+
+        const auto cell_stencil = stencil[cell_id];
+        for (size_t i_cell = 0; i_cell < cell_stencil.size(); ++i_cell) {
+          for (size_t dim = 0; dim < Dimension; ++dim) {
+            u_min[dim] = std::min(u_min[dim], u[cell_stencil[i_cell]][dim]);
+            u_max[dim] = std::max(u_max[dim], u[cell_stencil[i_cell]][dim]);
+          }
+        }
+
+        const double eps = 1E-14;
+        TinyVector<Dimension> coef;
+        double lambda = 1.;
+
+        for (size_t i_cell = 0; i_cell < cell_stencil.size(); ++i_cell) {
+          const CellId cell_i_id              = cell_stencil[i_cell];
+          const TinyVector<Dimension> u_bar_i = DPk_uh[cell_id](m_xj[cell_i_id]);
+          TinyVector<Dimension> Delta;
+
+          for (size_t dim = 0; dim < Dimension; ++dim) {
+            Delta[dim] = (u_bar_i[dim] - uj[dim]);
+            coef[dim]  = 1.;
+            if (Delta[dim] > eps) {
+              coef[dim] = std::min(1., (u_max[dim] - uj[dim]) / Delta[dim]);
+            } else if (Delta[dim] < -eps) {
+              coef[dim] = std::min(1., (u_min[dim] - uj[dim]) / Delta[dim]);
+            }
+            lambda = std::min(lambda, coef[dim]);
+          }
+        }
+        auto coefficients = DPk_uh.coefficients(cell_id);
+
+        coefficients[0] = (1 - lambda) * u[cell_id] + lambda * coefficients[0];
+        for (size_t i = 1; i < coefficients.size(); ++i) {
+          coefficients[i] *= lambda;
+        }
+      });
+  }
+
   std::tuple<std::shared_ptr<const DiscreteFunctionVariant>,
              std::shared_ptr<const DiscreteFunctionVariant>,
              std::shared_ptr<const DiscreteFunctionVariant>>
@@ -1635,6 +1919,8 @@ class EulerKineticSolverMultiDOrder2
 
     CellValue<bool> is_inflow_boundary_cell = getInflowBoundaryCells(bc_list);
 
+    const bool equilibrium_reconstruction = true;
+
     const CellValue<double> rho_ex{p_mesh->connectivity()};
     const CellValue<TinyVector<Dimension>> rho_u_ex{p_mesh->connectivity()};
     const CellValue<double> rho_E_ex{p_mesh->connectivity()};
@@ -1692,8 +1978,6 @@ class EulerKineticSolverMultiDOrder2
           if constexpr (std::is_same_v<BCType, WallBoundaryCondition>) {
             auto face_list = bc.faceList();
 
-            is_wall_boundary_face.fill(0);
-
             for (size_t i_face = 0; i_face < face_list.size(); ++i_face) {
               is_wall_boundary_face[face_list[i_face]] = 1;
             }
@@ -1702,6 +1986,29 @@ class EulerKineticSolverMultiDOrder2
         bc_v);
     }
 
+    CellValue<bool> is_wall_boundary_cell{p_mesh->connectivity()};
+    is_wall_boundary_cell.fill(false);
+    for (auto&& bc_v : bc_list) {
+      std::visit(
+        [=](auto&& bc) {
+          using BCType = std::decay_t<decltype(bc)>;
+          if constexpr (std::is_same_v<BCType, WallBoundaryCondition>) {
+            auto face_list                  = bc.faceList();
+            const auto& face_to_cell_matrix = m_mesh.connectivity().faceToCellMatrix();
+
+            for (size_t i_face = 0; i_face < face_list.size(); ++i_face) {
+              FaceId face_id          = face_list[i_face];
+              const auto face_to_cell = face_to_cell_matrix[face_id];
+              for (size_t i_cell = 0; i_cell < face_to_cell.size(); ++i_cell) {
+                CellId cell_id                 = face_to_cell[i_cell];
+                is_wall_boundary_cell[cell_id] = 1;
+              }
+            }
+          }
+        },
+        bc_v);
+    }
+
     const CellValue<const TinyVector<Dimension>> A_rho_ex   = getA_rho(rho_u_ex);
     const CellValue<const TinyMatrix<Dimension>> A_rho_u_ex = getA_rho_u(rho_ex, rho_u_ex, rho_E_ex);
     const CellValue<const TinyVector<Dimension>> A_rho_E_ex = getA_rho_E(rho_ex, rho_u_ex, rho_E_ex);
@@ -1718,8 +2025,6 @@ class EulerKineticSolverMultiDOrder2
     CellArray<TinyVector<Dimension>> Fnp1_rho_u = copy(Fn_rho_u);
     CellArray<double> Fnp1_rho_E                = copy(Fn_rho_E);
 
-    // Compute first order F
-
     const CellValue<const double> rho                  = compute_PFn<double>(Fn_rho);
     const CellValue<const TinyVector<Dimension>> rho_u = compute_PFn<TinyVector<Dimension>>(Fn_rho_u);
     const CellValue<const double> rho_E                = compute_PFn<double>(Fn_rho_E);
@@ -1738,91 +2043,80 @@ class EulerKineticSolverMultiDOrder2
 
     std::vector<std::shared_ptr<const IBoundaryDescriptor>> boundary_descriptor_list;
 
-    PolynomialReconstructionDescriptor reconstruction_descriptor_scalar(IntegrationMethodType::cell_center, 1,
-                                                                        symmetry_boundary_descriptor_list);
+    PolynomialReconstructionDescriptor reconstruction_descriptor(IntegrationMethodType::cell_center, 1,
+                                                                 symmetry_boundary_descriptor_list);
+    auto reconstruction_distributions =
+      PolynomialReconstruction{reconstruction_descriptor}.build(DiscreteFunctionP0Vector(p_mesh, Fn_rho),
+                                                                DiscreteFunctionP0Vector(p_mesh, Fn_rho_u),
+                                                                DiscreteFunctionP0Vector(p_mesh, Fn_rho_E));
 
-    PolynomialReconstructionDescriptor reconstruction_descriptor_vector(IntegrationMethodType::cell_center, 1,
-                                                                        symmetry_boundary_descriptor_list);
-
-    auto reconstruction_scalar =
-      PolynomialReconstruction{reconstruction_descriptor_scalar}.build(DiscreteFunctionP0Vector(p_mesh, Fn_rho),
-                                                                       // DiscreteFunctionP0Vector(p_mesh, Fn_rho_u),
-                                                                       DiscreteFunctionP0Vector(p_mesh, Fn_rho_E));
-
-    auto reconstruction_vector =
-      PolynomialReconstruction{reconstruction_descriptor_vector}.build(DiscreteFunctionP0Vector(p_mesh, Fn_rho_u));
-
-    auto DPk_Fn_rho   = reconstruction_scalar[0]->template get<DiscreteFunctionDPkVector<Dimension, const double>>();
-    auto DPk_Fn_rho_E = reconstruction_scalar[1]->template get<DiscreteFunctionDPkVector<Dimension, const double>>();
-
-    auto DPk_Fn_rho_u =
-      reconstruction_vector[0]->template get<DiscreteFunctionDPkVector<Dimension, const TinyVector<Dimension>>>();
-
-    // const auto& cell_to_face_matrix = m_mesh.connectivity().cellToFaceMatrix();
-    // // const auto& face_cell_is_reversed = p_mesh->connectivity().cellFaceIsReversed();
-
-    // for (CellId cell_id = 0; cell_id < m_mesh.numberOfCells(); ++cell_id) {
-    //   const auto cell_to_face = cell_to_face_matrix[cell_id];
-    //   if ((cell_id >= 500) and (cell_id < 510)) {
-    //     size_t i_face = 2;   // for (size_t i_face = 0; i_face < cell_to_face.size(); ++i_face) {
-
-    //     FaceId face_id = cell_to_face[i_face];
-
-    //     // const double sign = face_cell_is_reversed(cell_id, i_face) ? -1 : 1;
-    //     // const TinyVector<Dimension> n = sign * nl[face_id];
-
-    //     // for (size_t i_wave = 0; i_wave < m_lambda_vector.size(); ++i_wave) {
-    //     size_t i_wave = 0;
-    //     // std::cout << "n = " << n << ", x = " << m_xj[cell_id] << ", lambda_1 = " << m_lambda_vector[i_wave] <<
-    //     "\n"; std::cout << "F^rho_" << i_wave << "(" << cell_id << "," << face_id
-    //               << ") = " << DPk_Fn_rho(cell_id, i_wave)(m_xl[face_id]) << "; F^rho_u_" << i_wave << "(" << cell_id
-    //               << "," << face_id << ") = " << DPk_Fn_rho_u(cell_id, i_wave)(m_xl[face_id]) << "; F^rho_E_" <<
-    //               i_wave
-    //               << "(" << cell_id << "," << face_id << ") = " << DPk_Fn_rho_E(cell_id, i_wave)(m_xl[face_id]) <<
-    //               "\n";
-    //     // }
-    //     //}
-    //   }
-    // }
-    // std::cout << "\n";
-    // //    exit(0);
+    auto DPk_Fn_rho =
+      reconstruction_distributions[0]->template get<DiscreteFunctionDPkVector<Dimension, const double>>();
+    auto DPk_Fn_rho_u = reconstruction_distributions[1]
+                          ->template get<DiscreteFunctionDPkVector<Dimension, const TinyVector<Dimension>>>();
+    auto DPk_Fn_rho_E =
+      reconstruction_distributions[2]->template get<DiscreteFunctionDPkVector<Dimension, const double>>();
 
     DiscreteFunctionDPkVector<Dimension, double> Fn_rho_lim   = copy(DPk_Fn_rho);
     DiscreteFunctionDPkVector<Dimension, double> Fn_rho_E_lim = copy(DPk_Fn_rho_E);
 
-    scalar_limiter(Fn_rho, Fn_rho_lim);
-    scalar_limiter(Fn_rho_E, Fn_rho_E_lim);
+    scalar_limiter_distributions(Fn_rho, Fn_rho_lim);
+    scalar_limiter_distributions(Fn_rho_E, Fn_rho_E_lim);
 
     DiscreteFunctionDPkVector<Dimension, TinyVector<Dimension>> Fn_rho_u_lim = copy(DPk_Fn_rho_u);
 
-    vector_limiter(Fn_rho_u, Fn_rho_u_lim);
-
-    // const auto& cell_to_face_matrix = m_mesh.connectivity().cellToFaceMatrix();
-    // for (CellId cell_id = 500; cell_id < 502; ++cell_id) {
-    //   std::cout << "rho(" << cell_id << ") = " << rho[cell_id] << "\n";
-    //   std::cout << "rho_u(" << cell_id << ") = " << rho_u[cell_id] << "\n";
-    //   std::cout << "rho_E(" << cell_id << ") = " << rho_E[cell_id] << "\n";
-    //   const auto cell_to_face = cell_to_face_matrix[cell_id];
-    //   FaceId face_id          = cell_to_face[0];
-    //   for (size_t i_wave = 0; i_wave < nb_waves; ++i_wave) {
-    //     std::cout << "F_rho_u(" << i_wave << ") = " << Fn_rho_u_lim(cell_id, i_wave)(m_xl[face_id]) << "\n";
-    //   }
-    // }
-    // std::cout << "\n";
+    vector_limiter_distributions(Fn_rho_u, Fn_rho_u_lim);
+
+    auto reconstruction_equilibrium =
+      PolynomialReconstruction{reconstruction_descriptor}.build(DiscreteFunctionP0(p_mesh, rho),
+                                                                DiscreteFunctionP0(p_mesh, rho_u),
+                                                                DiscreteFunctionP0(p_mesh, rho_E));
+
+    auto DPk_rho = reconstruction_equilibrium[0]->template get<DiscreteFunctionDPk<Dimension, const double>>();
+    auto DPk_rho_u =
+      reconstruction_equilibrium[1]->template get<DiscreteFunctionDPk<Dimension, const TinyVector<Dimension>>>();
+    auto DPk_rho_E = reconstruction_equilibrium[2]->template get<DiscreteFunctionDPk<Dimension, const double>>();
+
+    DiscreteFunctionDPk<Dimension, double> rho_lim   = copy(DPk_rho);
+    DiscreteFunctionDPk<Dimension, double> rho_E_lim = copy(DPk_rho_E);
+
+    scalar_limiter_equilibrium(rho, rho_lim);
+    scalar_limiter_equilibrium(rho_E, rho_E_lim);
+
+    DiscreteFunctionDPk<Dimension, TinyVector<Dimension>> rho_u_lim = copy(DPk_rho_u);
+
+    vector_limiter_equilibrium(rho_u, rho_u_lim);
 
     CellArray<double> deltaLFn_rho{p_mesh->connectivity(), nb_waves};
     CellArray<TinyVector<Dimension>> deltaLFn_rho_u{p_mesh->connectivity(), nb_waves};
     CellArray<double> deltaLFn_rho_E{p_mesh->connectivity(), nb_waves};
 
     if (m_scheme == 1) {
-      FaceArray<double> Fl_rho = compute_Flux_face<double>(Fn_rho_lim, Fn_rho, is_wall_boundary_face);
-      FaceArray<TinyVector<Dimension>> Fl_rho_u =
-        compute_Flux_face<TinyVector<Dimension>>(Fn_rho_u_lim, Fn_rho_u, is_wall_boundary_face);
-      FaceArray<double> Fl_rho_E = compute_Flux_face<double>(Fn_rho_E_lim, Fn_rho_E, is_wall_boundary_face);
+      FaceArray<double> Fl_rho{p_mesh->connectivity(), nb_waves};
+      FaceArray<TinyVector<Dimension>> Fl_rho_u{p_mesh->connectivity(), nb_waves};
+      FaceArray<double> Fl_rho_E{p_mesh->connectivity(), nb_waves};
+
+      if (not equilibrium_reconstruction) {
+        std::tuple<FaceArray<double>, FaceArray<TinyVector<Dimension>>, FaceArray<double>> vec_Fl(Fl_rho, Fl_rho_u,
+                                                                                                  Fl_rho_E);
 
-      apply_face_bc(Fl_rho, Fl_rho_u, Fl_rho_E, Fn_rho_lim, Fn_rho_u_lim, Fn_rho_E_lim, Fn_rho, Fn_rho_u, Fn_rho_E,
-                    bc_list);
+        vec_Fl = compute_Flux_face_equilibrium(rho_lim, rho_u_lim, rho_E_lim, Fn_rho, Fn_rho_u, Fn_rho_E,
+                                               is_wall_boundary_face);
 
+        Fl_rho   = get<0>(vec_Fl);
+        Fl_rho_u = get<1>(vec_Fl);
+        Fl_rho_E = get<2>(vec_Fl);
+
+        apply_face_bc_equilibrium(Fl_rho, Fl_rho_u, Fl_rho_E, rho_lim, rho_u_lim, rho_E_lim, Fn_rho, Fn_rho_u, Fn_rho_E,
+                                  bc_list);
+      } else {
+        Fl_rho   = compute_Flux_face<double>(Fn_rho_lim, Fn_rho, is_wall_boundary_face);
+        Fl_rho_u = compute_Flux_face<TinyVector<Dimension>>(Fn_rho_u_lim, Fn_rho_u, is_wall_boundary_face);
+        Fl_rho_E = compute_Flux_face<double>(Fn_rho_E_lim, Fn_rho_E, is_wall_boundary_face);
+
+        apply_face_bc(Fl_rho, Fl_rho_u, Fl_rho_E, Fn_rho_lim, Fn_rho_u_lim, Fn_rho_E_lim, Fn_rho, Fn_rho_u, Fn_rho_E,
+                      bc_list);
+      }
       synchronize(Fl_rho);
       synchronize(Fl_rho_u);
       synchronize(Fl_rho_E);
@@ -1873,59 +2167,85 @@ class EulerKineticSolverMultiDOrder2
     const CellArray<const double> M_rho_E                = compute_scalar_M(rho_E, A_rho_E);
 
     for (size_t i_subtimestep = 0; i_subtimestep < m_TimeOrder; i_subtimestep++) {
-      auto reconstruction_np1_scalar =
-        PolynomialReconstruction{reconstruction_descriptor_scalar}.build(DiscreteFunctionP0Vector(p_mesh, Fnp1_rho),
-                                                                         DiscreteFunctionP0Vector(p_mesh, Fnp1_rho_E));
-      auto reconstruction_np1_vector =
-        PolynomialReconstruction{reconstruction_descriptor_vector}.build(DiscreteFunctionP0Vector(p_mesh, Fnp1_rho_u));
-
-      auto DPk_Fnp1_rho = reconstruction_scalar[0]->template get<DiscreteFunctionDPkVector<Dimension, const double>>();
+      auto reconstruction_distributions_np1 =
+        PolynomialReconstruction{reconstruction_descriptor}.build(DiscreteFunctionP0Vector(p_mesh, Fnp1_rho),
+                                                                  DiscreteFunctionP0Vector(p_mesh, Fnp1_rho_u),
+                                                                  DiscreteFunctionP0Vector(p_mesh, Fnp1_rho_E));
+
+      auto DPk_Fnp1_rho =
+        reconstruction_distributions_np1[0]->template get<DiscreteFunctionDPkVector<Dimension, const double>>();
+      auto DPk_Fnp1_rho_u = reconstruction_distributions_np1[1]
+                              ->template get<DiscreteFunctionDPkVector<Dimension, const TinyVector<Dimension>>>();
       auto DPk_Fnp1_rho_E =
-        reconstruction_scalar[1]->template get<DiscreteFunctionDPkVector<Dimension, const double>>();
-
-      auto DPk_Fnp1_rho_u =
-        reconstruction_vector[0]->template get<DiscreteFunctionDPkVector<Dimension, const TinyVector<Dimension>>>();
+        reconstruction_distributions_np1[2]->template get<DiscreteFunctionDPkVector<Dimension, const double>>();
 
       DiscreteFunctionDPkVector<Dimension, double> Fnp1_rho_lim   = copy(DPk_Fnp1_rho);
       DiscreteFunctionDPkVector<Dimension, double> Fnp1_rho_E_lim = copy(DPk_Fnp1_rho_E);
 
-      // for (CellId cell_id = 0; cell_id < m_mesh.numberOfCells(); ++cell_id) {
-      //   const auto cell_to_face = cell_to_face_matrix[cell_id];
-      //   if ((cell_id >= 500) and (cell_id < 510)) {
-      //     for (size_t i_face = 0; i_face < cell_to_face.size(); ++i_face) {
-      //       FaceId face_id = cell_to_face[i_face];
+      scalar_limiter_distributions(Fnp1_rho, Fnp1_rho_lim);
+      scalar_limiter_distributions(Fnp1_rho_E, Fnp1_rho_E_lim);
 
-      //       std::cout << "F^rho_1(" << cell_id << "," << face_id << ") = " << DPk_Fnp1_rho(cell_id,
-      //       0)(m_xl[face_id]); std::cout << "F^rho_u_1(" << cell_id << "," << face_id
-      //                 << ") = " << DPk_Fnp1_rho_u(cell_id, 0)(m_xl[face_id]);
+      DiscreteFunctionDPkVector<Dimension, TinyVector<Dimension>> Fnp1_rho_u_lim = copy(DPk_Fnp1_rho_u);
 
-      //       std::cout << "F^rho_E_1(" << cell_id << "," << face_id
-      //                 << ") = " << DPk_Fnp1_rho_E(cell_id, 0)(m_xl[face_id]);
-      //     }
-      //   }
-      // }
-      // exit(0);
+      vector_limiter_distributions(Fnp1_rho_u, Fnp1_rho_u_lim);
 
-      scalar_limiter(Fnp1_rho, Fnp1_rho_lim);
-      scalar_limiter(Fnp1_rho_E, Fnp1_rho_E_lim);
+      rho_np1   = compute_PFn<double>(Fnp1_rho);
+      rho_u_np1 = compute_PFn<TinyVector<Dimension>>(Fnp1_rho_u);
+      rho_E_np1 = compute_PFn<double>(Fnp1_rho_E);
 
-      DiscreteFunctionDPkVector<Dimension, TinyVector<Dimension>> Fnp1_rho_u_lim = copy(DPk_Fnp1_rho_u);
+      auto reconstruction_equilibrium_np1 =
+        PolynomialReconstruction{reconstruction_descriptor}.build(DiscreteFunctionP0(p_mesh, rho_np1),
+                                                                  DiscreteFunctionP0(p_mesh, rho_u_np1),
+                                                                  DiscreteFunctionP0(p_mesh, rho_E_np1));
+
+      auto DPk_rho_np1 =
+        reconstruction_equilibrium_np1[0]->template get<DiscreteFunctionDPk<Dimension, const double>>();
+      auto DPk_rho_u_np1 =
+        reconstruction_equilibrium_np1[1]->template get<DiscreteFunctionDPk<Dimension, const TinyVector<Dimension>>>();
+      auto DPk_rho_E_np1 =
+        reconstruction_equilibrium_np1[2]->template get<DiscreteFunctionDPk<Dimension, const double>>();
+
+      DiscreteFunctionDPk<Dimension, double> rho_np1_lim   = copy(DPk_rho_np1);
+      DiscreteFunctionDPk<Dimension, double> rho_E_np1_lim = copy(DPk_rho_E_np1);
 
-      vector_limiter(Fnp1_rho_u, Fnp1_rho_u_lim);
+      scalar_limiter_equilibrium(rho_np1, rho_np1_lim);
+      scalar_limiter_equilibrium(rho_E_np1, rho_E_np1_lim);
+
+      DiscreteFunctionDPk<Dimension, TinyVector<Dimension>> rho_u_np1_lim = copy(DPk_rho_u_np1);
+
+      vector_limiter_equilibrium(rho_u_np1, rho_u_np1_lim);
 
       CellArray<double> deltaLFnp1_rho{p_mesh->connectivity(), nb_waves};
       CellArray<TinyVector<Dimension>> deltaLFnp1_rho_u{p_mesh->connectivity(), nb_waves};
       CellArray<double> deltaLFnp1_rho_E{p_mesh->connectivity(), nb_waves};
 
       if (m_scheme == 1) {
-        FaceArray<double> Fl_rho_np1 = compute_Flux_face<double>(Fnp1_rho_lim, Fn_rho, is_wall_boundary_face);
-        FaceArray<TinyVector<Dimension>> Fl_rho_u_np1 =
-          compute_Flux_face<TinyVector<Dimension>>(Fnp1_rho_u_lim, Fn_rho_u, is_wall_boundary_face);
-        FaceArray<double> Fl_rho_E_np1 = compute_Flux_face<double>(Fnp1_rho_E_lim, Fn_rho_E, is_wall_boundary_face);
-
-        apply_face_bc(Fl_rho_np1, Fl_rho_u_np1, Fl_rho_E_np1, Fnp1_rho_lim, Fnp1_rho_u_lim, Fnp1_rho_E_lim, Fn_rho,
-                      Fn_rho_u, Fn_rho_E, bc_list);
-
+        FaceArray<double> Fl_rho_np1{p_mesh->connectivity(), nb_waves};
+        FaceArray<TinyVector<Dimension>> Fl_rho_u_np1{p_mesh->connectivity(), nb_waves};
+        FaceArray<double> Fl_rho_E_np1{p_mesh->connectivity(), nb_waves};
+
+        if (not equilibrium_reconstruction) {
+          std::tuple<FaceArray<double>, FaceArray<TinyVector<Dimension>>, FaceArray<double>> vec_Fl_np1(Fl_rho_np1,
+                                                                                                        Fl_rho_u_np1,
+                                                                                                        Fl_rho_E_np1);
+
+          vec_Fl_np1 = compute_Flux_face_equilibrium(rho_np1_lim, rho_u_np1_lim, rho_E_np1_lim, Fnp1_rho, Fnp1_rho_u,
+                                                     Fnp1_rho_E, is_wall_boundary_face);
+
+          Fl_rho_np1   = get<0>(vec_Fl_np1);
+          Fl_rho_u_np1 = get<1>(vec_Fl_np1);
+          Fl_rho_E_np1 = get<2>(vec_Fl_np1);
+
+          apply_face_bc_equilibrium(Fl_rho_np1, Fl_rho_u_np1, Fl_rho_E_np1, rho_lim, rho_u_lim, rho_E_lim, Fn_rho,
+                                    Fn_rho_u, Fn_rho_E, bc_list);
+        } else {
+          Fl_rho_np1   = compute_Flux_face<double>(Fnp1_rho_lim, Fn_rho, is_wall_boundary_face);
+          Fl_rho_u_np1 = compute_Flux_face<TinyVector<Dimension>>(Fnp1_rho_u_lim, Fn_rho_u, is_wall_boundary_face);
+          Fl_rho_E_np1 = compute_Flux_face<double>(Fnp1_rho_E_lim, Fn_rho_E, is_wall_boundary_face);
+
+          apply_face_bc(Fl_rho_np1, Fl_rho_u_np1, Fl_rho_E_np1, Fnp1_rho_lim, Fnp1_rho_u_lim, Fnp1_rho_E_lim, Fn_rho,
+                        Fn_rho_u, Fn_rho_E, bc_list);
+        }
         synchronize(Fl_rho_np1);
         synchronize(Fl_rho_u_np1);
         synchronize(Fl_rho_E_np1);