diff --git a/src/scheme/EulerKineticSolverMultiD.cpp b/src/scheme/EulerKineticSolverMultiD.cpp
index 0f9acb6513f58837c480dbd2e5271fc42c5de64a..1db8117321328f070f2b085fa309c6be4e21514e 100644
--- a/src/scheme/EulerKineticSolverMultiD.cpp
+++ b/src/scheme/EulerKineticSolverMultiD.cpp
@@ -242,6 +242,7 @@ class EulerKineticSolverMultiD
   CellValue<const double> m_inv_Vj;
   const CellValue<const TinyVector<Dimension>> m_xj = MeshDataManager::instance().getMeshData(m_mesh).xj();
   const NodeValue<const TinyVector<Dimension>> m_xr = m_mesh.xr();
+  const FaceValue<const TinyVector<Dimension>> m_xl = MeshDataManager::instance().getMeshData(m_mesh).xl();
 
   BoundaryConditionList
   _getBCList(const MeshType& mesh,
@@ -649,6 +650,16 @@ class EulerKineticSolverMultiD
 
             if (li_nl > 1e-14) {
               Fl[face_id][i_wave] += Fn[cell_id][i_wave];
+
+              if ((face_id == 51) and (i_wave == 1 or i_wave == 2)) {
+                T test;
+                if constexpr (not is_tiny_vector_v<T>) {
+                  test = 0;
+                }
+                for (size_t k_wave = 0; k_wave < nb_waves; ++k_wave) {
+                  test += Fn[cell_id][k_wave];
+                }
+              }
             }
           }
         }
@@ -704,32 +715,33 @@ class EulerKineticSolverMultiD
 
                   double li_nl = dot(m_lambda_vector[i_wave], n);
 
-                  if (li_nl < -1e-14) {
-                    double rhoj_prime                  = rho[cell_id];
-                    TinyVector<Dimension> rho_uj       = rho_u[cell_id];
-                    TinyVector<Dimension> rho_uj_prime = txt * rho_uj - nxn * rho_uj;
-                    double rho_Ej_prime                = rho_E[cell_id];
+                  //                  if (li_nl < -1e-14) {
+                  double rhoj_prime                  = rho[cell_id];
+                  TinyVector<Dimension> rho_uj       = rho_u[cell_id];
+                  TinyVector<Dimension> rho_uj_prime = txt * rho_uj - nxn * rho_uj;
 
-                    double rho_e = rho_Ej_prime - 0.5 * (1. / rhoj_prime) * dot(rho_uj, rho_uj);
-                    double p     = (m_gamma - 1) * rho_e;
+                  double rho_Ej_prime = rho_E[cell_id];
 
-                    TinyVector<Dimension> A_rho_prime = rho_uj_prime;
-                    TinyMatrix<Dimension> A_rho_u_prime =
-                      1. / rhoj_prime * tensorProduct(rho_uj_prime, rho_uj_prime) + p * I;
-                    TinyVector<Dimension> A_rho_E_prime = (rho_Ej_prime + p) / rhoj_prime * rho_uj_prime;
+                  double rho_e = rho_Ej_prime - 0.5 * (1. / rhoj_prime) * dot(rho_uj, rho_uj);
+                  double p     = (m_gamma - 1) * rho_e;
 
-                    double Fn_rho_prime                  = rhoj_prime / nb_waves;
-                    TinyVector<Dimension> Fn_rho_u_prime = 1. / nb_waves * rho_uj_prime;
-                    double Fn_rho_E_prime                = rho_Ej_prime / nb_waves;
+                  TinyVector<Dimension> A_rho_prime = rho_uj_prime;
+                  TinyMatrix<Dimension> A_rho_u_prime =
+                    1. / rhoj_prime * tensorProduct(rho_uj_prime, rho_uj_prime) + p * I;
+                  TinyVector<Dimension> A_rho_E_prime = (rho_Ej_prime + p) / rhoj_prime * rho_uj_prime;
 
-                    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];
-                    }
+                  double Fn_rho_prime                  = rhoj_prime / nb_waves;
+                  TinyVector<Dimension> Fn_rho_u_prime = 1. / nb_waves * rho_uj_prime;
+                  double Fn_rho_E_prime                = rho_Ej_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];
+                  }
+                  if (li_nl < -1e-14) {
                     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;
@@ -874,7 +886,6 @@ 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;
@@ -1474,6 +1485,7 @@ class EulerKineticSolverMultiD
             }
           }
         });
+
       return deltaLFn;
     } else {
       throw NormalError("Nl in meaningless in 1d");
@@ -1515,6 +1527,7 @@ class EulerKineticSolverMultiD
     const size_t nb_waves = m_lambda_vector.size();
 
     CellValue<bool> is_inflow_boundary_cell = getInflowBoundaryCells(bc_list);
+    // const auto& cell_to_face_matrix         = p_mesh->connectivity().cellToFaceMatrix();
 
     const CellValue<double> rho_ex{p_mesh->connectivity()};
     const CellValue<TinyVector<Dimension>> rho_u_ex{p_mesh->connectivity()};
diff --git a/src/scheme/EulerKineticSolverMultiDOrder2.cpp b/src/scheme/EulerKineticSolverMultiDOrder2.cpp
index 030a0bab331a41ac32dc85a479fa3af32f4e34b6..e6827ea1a20c59946c772509db1a580a9c9e4c1d 100644
--- a/src/scheme/EulerKineticSolverMultiDOrder2.cpp
+++ b/src/scheme/EulerKineticSolverMultiDOrder2.cpp
@@ -383,35 +383,35 @@ class EulerKineticSolverMultiDOrder2
             const CellId cell_id     = face_to_cell[i_cell];
             const size_t i_face_cell = face_local_number_in_its_cell[i_cell];
 
-            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]);
+            const TinyVector<Dimension> n_face = nl[face_id];
+            const double sign                  = face_cell_is_reversed(cell_id, i_face_cell) ? -1 : 1;
 
-            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;
+            const double li_nl = sign * dot(m_lambda_vector[i_wave], n_face);
 
-            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;
+            if (li_nl > 1e-14) {
+              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]);
 
-            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;
+              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;
 
-            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];
-            }
+              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;
 
-            TinyVector<Dimension> n_face = nl[face_id];
-            const double sign            = face_cell_is_reversed(cell_id, i_face_cell) ? -1 : 1;
+              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;
 
-            double li_nl = sign * dot(m_lambda_vector[i_wave], n_face);
+              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];
+              }
 
-            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;
@@ -524,6 +524,95 @@ class EulerKineticSolverMultiDOrder2
     }
   }
 
+  std::tuple<NodeArray<double>, NodeArray<TinyVector<Dimension>>, NodeArray<double>>
+  compute_Flux_glace_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,
+                                 NodeValue<bool> is_boundary_node)
+  {
+    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<double> Fr_rho(m_mesh.connectivity(), nb_waves);
+      NodeArray<TinyVector<Dimension>> Fr_rho_u(m_mesh.connectivity(), nb_waves);
+      NodeArray<double> Fr_rho_E(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();
+      TinyMatrix<Dimension> I                       = identity;
+
+      Fr_rho.fill(0.);
+      Fr_rho_u.fill(zero);
+      Fr_rho_E.fill(0.);
+
+      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];
+      }
+
+      parallel_for(
+        p_mesh->numberOfNodes(), PUGS_LAMBDA(NodeId node_id) {
+          if (not is_boundary_node[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];
+
+                double li_njr = dot(m_lambda_vector[i_wave], njr(cell_id, i_node));
+
+                if (li_njr > 0) {
+                  const double rho_jr                  = DPk_rho[cell_id](m_xr[node_id]);
+                  const TinyVector<Dimension> rho_u_jr = DPk_rho_u[cell_id](m_xr[node_id]);
+                  const double rho_E_jr                = DPk_rho_E[cell_id](m_xr[node_id]);
+
+                  const double rho_e_jr = rho_E_jr - 0.5 * (1. / rho_jr) * dot(rho_u_jr, rho_u_jr);
+                  const double p_jr     = (m_gamma - 1) * rho_e_jr;
+
+                  const TinyVector<Dimension> A_rho_jr   = rho_u_jr;
+                  const TinyMatrix<Dimension> A_rho_u_jr = 1. / rho_jr * tensorProduct(rho_u_jr, rho_u_jr) + p_jr * I;
+                  const TinyVector<Dimension> A_rho_E_jr = (rho_E_jr + p_jr) / rho_jr * rho_u_jr;
+
+                  double Fn_rho_jr                  = rho_jr / nb_waves;
+                  TinyVector<Dimension> Fn_rho_u_jr = 1. / nb_waves * rho_u_jr;
+                  double Fn_rho_E_jr                = rho_E_jr / nb_waves;
+
+                  for (size_t d1 = 0; d1 < Dimension; ++d1) {
+                    Fn_rho_jr += inv_S[d1] * m_lambda_vector[i_wave][d1] * A_rho_jr[d1];
+                    for (size_t d2 = 0; d2 < Dimension; ++d2) {
+                      Fn_rho_u_jr[d1] += inv_S[d2] * m_lambda_vector[i_wave][d2] * A_rho_u_jr(d2, d1);
+                    }
+                    Fn_rho_E_jr += inv_S[d1] * m_lambda_vector[i_wave][d1] * A_rho_E_jr[d1];
+                  }
+
+                  Fr_rho[node_id][i_wave] += li_njr * Fn_rho_jr;
+                  Fr_rho_u[node_id][i_wave] += li_njr * Fn_rho_u_jr;
+                  Fr_rho_E[node_id][i_wave] += li_njr * Fn_rho_E_jr;
+
+                  sum_li_njr += li_njr;
+                }
+              }
+              if (sum_li_njr != 0) {
+                Fr_rho[node_id][i_wave]   = 1. / sum_li_njr * Fr_rho[node_id][i_wave];
+                Fr_rho_u[node_id][i_wave] = 1. / sum_li_njr * Fr_rho_u[node_id][i_wave];
+                Fr_rho_E[node_id][i_wave] = 1. / sum_li_njr * Fr_rho_E[node_id][i_wave];
+              }
+            }
+          }
+        });
+
+      return std::make_tuple(Fr_rho, Fr_rho_u, Fr_rho_E);
+    } else {
+      throw NormalError("Glace not defined in 1d");
+    }
+  }
+
   template <typename T>
   NodeArray<T>
   compute_Flux_eucclhyd(const DiscreteFunctionDPkVector<Dimension, const T> DPk_Fn, NodeValue<bool> is_boundary_node)
@@ -648,45 +737,38 @@ class EulerKineticSolverMultiDOrder2
                   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;
+                  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;
+                    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;
+                    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;
+                    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);
+                    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];
                     }
-                    Fn_rho_E_prime += inv_S[d1] * m_lambda_vector[i_wave][d1] * A_rho_E_prime[d1];
-                  }
-                  if (li_nl < -1e-14) {
                     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;
                   }
-                  // if (cell_id >= 500 and cell_id < 502) {
-                  //   std::cout << "rho'(" << cell_id << ") = " << rhol_prime << "\n";
-                  //   std::cout << "rho_u'(" << cell_id << ") = " << rho_ul_prime << "\n";
-                  //   std::cout << "rho_E'(" << cell_id << ") = " << rho_El_prime << "\n";
-                  //   std::cout << "F_rho_u'(" << i_wave << ") = " << Fn_rho_u_prime << "\n";
-                  // }
                 }
               }
             }
-            //            std::cout << "\n";
+
           } else if constexpr (std::is_same_v<BCType, WallBoundaryCondition>) {
             auto face_list = bc.faceList();
 
@@ -1270,6 +1352,365 @@ class EulerKineticSolverMultiDOrder2
     }
   }
 
+  void
+  apply_glace_bc_equilibrium(NodeArray<double> Fr_rho,
+                             NodeArray<TinyVector<Dimension>> Fr_rho_u,
+                             NodeArray<double> Fr_rho_E,
+                             DiscreteFunctionDPk<Dimension, double> rho,
+                             DiscreteFunctionDPk<Dimension, TinyVector<Dimension>> rho_u,
+                             DiscreteFunctionDPk<Dimension, double> rho_E,
+                             const BoundaryConditionList& bc_list)
+  {
+    const size_t nb_waves                                   = m_lambda_vector.size();
+    const auto& node_local_numbers_in_their_cells           = p_mesh->connectivity().nodeLocalNumbersInTheirCells();
+    const auto& node_to_cell_matrix                         = p_mesh->connectivity().nodeToCellMatrix();
+    const NodeValuePerCell<const TinyVector<Dimension>> njr = MeshDataManager::instance().getMeshData(m_mesh).njr();
+    const NodeValuePerCell<const TinyVector<Dimension>> Cjr = MeshDataManager::instance().getMeshData(m_mesh).Cjr();
+
+    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 node_list = bc.nodeList();
+
+            TinyMatrix<Dimension> I = identity;
+
+            NodeValue<TinyVector<Dimension>> nr{p_mesh->connectivity()};
+            nr.fill(zero);
+
+            for (size_t i_node = 0; i_node < node_list.size(); ++i_node) {
+              const NodeId node_id                      = node_list[i_node];
+              const auto& node_to_cell                  = node_to_cell_matrix[node_id];
+              const auto& node_local_number_in_its_cell = node_local_numbers_in_their_cells.itemArray(node_id);
+
+              for (size_t i_cell = 0; i_cell < node_to_cell.size(); ++i_cell) {
+                const CellId cell_id     = node_to_cell[i_cell];
+                const size_t i_node_cell = node_local_number_in_its_cell[i_cell];
+
+                nr[node_id] += Cjr(cell_id, i_node_cell);
+              }
+              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;
+
+              for (size_t i_wave = 0; i_wave < nb_waves; ++i_wave) {
+                double sum                = 0;
+                Fr_rho[node_id][i_wave]   = 0;
+                Fr_rho_u[node_id][i_wave] = zero;
+                Fr_rho_E[node_id][i_wave] = 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 size_t i_node_cell = node_local_number_in_its_cell[i_cell];
+
+                  double rhor                  = 0;
+                  TinyVector<Dimension> rho_ur = zero;
+                  double rho_Er                = 0;
+
+                  rhor   = rho[cell_id](m_xr[node_id]);
+                  rho_ur = rho_u[cell_id](m_xr[node_id]);
+                  rho_Er = rho_E[cell_id](m_xr[node_id]);
+
+                  double li_njr = dot(m_lambda_vector[i_wave], njr(cell_id, i_node_cell));
+
+                  if (li_njr > 0) {
+                    const double rho_er = rho_Er - 0.5 * (1. / rhor) * dot(rho_ur, rho_ur);
+                    const double pr     = (m_gamma - 1) * rho_er;
+
+                    const TinyVector<Dimension> A_rho   = rho_ur;
+                    const TinyMatrix<Dimension> A_rho_u = 1. / rhor * tensorProduct(rho_ur, rho_ur) + pr * I;
+                    const TinyVector<Dimension> A_rho_E = (rho_Er + pr) / rhor * rho_ur;
+
+                    double Fn_rho                  = rhor / nb_waves;
+                    TinyVector<Dimension> Fn_rho_u = 1. / nb_waves * rho_ur;
+                    double Fn_rho_E                = rho_Er / 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];
+                    }
+
+                    Fr_rho[node_id][i_wave] += Fn_rho * li_njr;
+                    Fr_rho_u[node_id][i_wave] += li_njr * Fn_rho_u;
+                    Fr_rho_E[node_id][i_wave] += Fn_rho_E * li_njr;
+
+                    sum += li_njr;
+                  }
+
+                  TinyVector<Dimension> njr_prime = txt * njr(cell_id, i_node_cell) - nxn * njr(cell_id, i_node_cell);
+                  double li_njr_prime             = dot(m_lambda_vector[i_wave], njr_prime);
+
+                  if (li_njr_prime > 0) {
+                    const double rhor_prime                  = rhor;
+                    const TinyVector<Dimension> rho_ur_prime = txt * rho_ur - nxn * rho_ur;
+                    const double rho_Er_prime                = rho_Er;
+
+                    const double rho_er = rho_Er - 0.5 * (1. / rhor) * dot(rho_ur, rho_ur);
+                    const double pr     = (m_gamma - 1) * rho_er;
+
+                    const TinyVector<Dimension> A_rho_prime = rho_ur_prime;
+                    const TinyMatrix<Dimension> A_rho_u_prime =
+                      1. / rhor_prime * tensorProduct(rho_ur_prime, rho_ur_prime) + pr * I;
+                    const TinyVector<Dimension> A_rho_E_prime = (rho_Er_prime + pr) / rhor_prime * rho_ur_prime;
+
+                    double Fn_rho_prime                  = rhor_prime / nb_waves;
+                    TinyVector<Dimension> Fn_rho_u_prime = 1. / nb_waves * rho_ur_prime;
+                    double Fn_rho_E_prime                = rho_Er_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];
+                    }
+
+                    Fr_rho[node_id][i_wave] += Fn_rho_prime * li_njr_prime;
+                    Fr_rho_u[node_id][i_wave] += li_njr_prime * Fn_rho_u_prime;
+                    Fr_rho_E[node_id][i_wave] += Fn_rho_E_prime * li_njr_prime;
+
+                    sum += li_njr_prime;
+                  }
+                }
+                if (sum != 0) {
+                  Fr_rho[node_id][i_wave] /= sum;
+                  Fr_rho_u[node_id][i_wave] = 1. / sum * Fr_rho_u[node_id][i_wave];
+                  Fr_rho_E[node_id][i_wave] /= sum;
+                }
+              }
+            }
+          } else if constexpr (std::is_same_v<BCType, WallBoundaryCondition>) {
+            auto node_list = bc.nodeList();
+
+            TinyMatrix<Dimension> I = identity;
+
+            NodeValue<TinyVector<Dimension>> nr{p_mesh->connectivity()};
+            nr.fill(zero);
+
+            for (size_t i_node = 0; i_node < node_list.size(); ++i_node) {
+              const NodeId node_id                      = node_list[i_node];
+              const auto& node_to_cell                  = node_to_cell_matrix[node_id];
+              const auto& node_local_number_in_its_cell = node_local_numbers_in_their_cells.itemArray(node_id);
+
+              for (size_t i_cell = 0; i_cell < node_to_cell.size(); ++i_cell) {
+                const CellId cell_id     = node_to_cell[i_cell];
+                const size_t i_node_cell = node_local_number_in_its_cell[i_cell];
+
+                nr[node_id] += Cjr(cell_id, i_node_cell);
+              }
+              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;
+
+              for (size_t i_wave = 0; i_wave < nb_waves; ++i_wave) {
+                double sum                = 0;
+                Fr_rho[node_id][i_wave]   = 0;
+                Fr_rho_u[node_id][i_wave] = zero;
+                Fr_rho_E[node_id][i_wave] = 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 size_t i_node_cell = node_local_number_in_its_cell[i_cell];
+
+                  const double rhor                  = rho[cell_id](m_xr[node_id]);
+                  const TinyVector<Dimension> rho_ur = rho_u[cell_id](m_xr[node_id]);
+                  const double rho_Er                = rho_E[cell_id](m_xr[node_id]);
+
+                  double li_njr = dot(m_lambda_vector[i_wave], njr(cell_id, i_node_cell));
+
+                  if (li_njr > 0) {
+                    const double rho_er = rho_Er - 0.5 * (1. / rhor) * dot(rho_ur, rho_ur);
+                    const double pr     = (m_gamma - 1) * rho_er;
+
+                    const TinyVector<Dimension> A_rho   = rho_ur;
+                    const TinyMatrix<Dimension> A_rho_u = 1. / rhor * tensorProduct(rho_ur, rho_ur) + pr * I;
+                    const TinyVector<Dimension> A_rho_E = (rho_Er + pr) / rhor * rho_ur;
+
+                    double Fn_rho                  = rhor / nb_waves;
+                    TinyVector<Dimension> Fn_rho_u = 1. / nb_waves * rho_ur;
+                    double Fn_rho_E                = rho_Er / 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];
+                    }
+
+                    Fr_rho[node_id][i_wave] += Fn_rho * li_njr;
+                    Fr_rho_u[node_id][i_wave] += li_njr * Fn_rho_u;
+                    Fr_rho_E[node_id][i_wave] += Fn_rho_E * li_njr;
+
+                    sum += li_njr;
+                  }
+
+                  TinyVector<Dimension> njr_prime = txt * njr(cell_id, i_node_cell) - nxn * njr(cell_id, i_node_cell);
+                  double li_njr_prime             = dot(m_lambda_vector[i_wave], njr_prime);
+
+                  if (li_njr_prime > 0) {
+                    const double rhor_prime                  = rhor;
+                    const TinyVector<Dimension> rho_ur_prime = txt * rho_ur - nxn * rho_ur;
+                    const double rho_Er_prime                = rho_Er;
+
+                    const double rho_er = rho_Er - 0.5 * (1. / rhor) * dot(rho_ur, rho_ur);
+                    const double pr     = (m_gamma - 1) * rho_er;
+
+                    const TinyVector<Dimension> A_rho_prime = rho_ur_prime;
+                    const TinyMatrix<Dimension> A_rho_u_prime =
+                      1. / rhor_prime * tensorProduct(rho_ur_prime, rho_ur_prime) + pr * I;
+                    const TinyVector<Dimension> A_rho_E_prime = (rho_Er_prime + pr) / rhor_prime * rho_ur_prime;
+
+                    double Fn_rho_prime                  = rhor_prime / nb_waves;
+                    TinyVector<Dimension> Fn_rho_u_prime = 1. / nb_waves * rho_ur_prime;
+                    double Fn_rho_E_prime                = rho_Er_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];
+                    }
+
+                    Fr_rho[node_id][i_wave] += Fn_rho_prime * li_njr_prime;
+                    Fr_rho_u[node_id][i_wave] += li_njr_prime * Fn_rho_u_prime;
+                    Fr_rho_E[node_id][i_wave] += Fn_rho_E_prime * li_njr_prime;
+
+                    sum += li_njr_prime;
+                  }
+                }
+                if (sum != 0) {
+                  Fr_rho[node_id][i_wave] /= sum;
+                  Fr_rho_u[node_id][i_wave] = 1. / sum * Fr_rho_u[node_id][i_wave];
+                  Fr_rho_E[node_id][i_wave] /= sum;
+                }
+              }
+            }
+          } else if constexpr (std::is_same_v<BCType, OutflowBoundaryCondition>) {
+            auto node_list          = bc.nodeList();
+            TinyMatrix<Dimension> I = identity;
+
+            NodeValue<TinyVector<Dimension>> nr{p_mesh->connectivity()};
+            nr.fill(zero);
+
+            for (size_t i_node = 0; i_node < node_list.size(); ++i_node) {
+              const NodeId node_id                      = node_list[i_node];
+              const auto& node_to_cell                  = node_to_cell_matrix[node_id];
+              const auto& node_local_number_in_its_cell = node_local_numbers_in_their_cells.itemArray(node_id);
+
+              for (size_t i_cell = 0; i_cell < node_to_cell.size(); ++i_cell) {
+                const CellId cell_id     = node_to_cell[i_cell];
+                const size_t i_node_cell = node_local_number_in_its_cell[i_cell];
+
+                nr[node_id] += Cjr(cell_id, i_node_cell);
+              }
+              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;
+
+              for (size_t i_wave = 0; i_wave < nb_waves; ++i_wave) {
+                double sum                = 0;
+                Fr_rho[node_id][i_wave]   = 0;
+                Fr_rho_u[node_id][i_wave] = zero;
+                Fr_rho_E[node_id][i_wave] = 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 size_t i_node_cell = node_local_number_in_its_cell[i_cell];
+
+                  double li_njr = dot(m_lambda_vector[i_wave], njr(cell_id, i_node_cell));
+
+                  if (li_njr > 0) {
+                    const double rhor                  = rho[cell_id](m_xr[node_id]);
+                    const TinyVector<Dimension> rho_ur = rho_u[cell_id](m_xr[node_id]);
+                    const double rho_Er                = rho_E[cell_id](m_xr[node_id]);
+
+                    const double rho_er = rho_Er - 0.5 * (1. / rhor) * dot(rho_ur, rho_ur);
+                    const double pr     = (m_gamma - 1) * rho_er;
+
+                    const TinyVector<Dimension> A_rho   = rho_ur;
+                    const TinyMatrix<Dimension> A_rho_u = 1. / rhor * tensorProduct(rho_ur, rho_ur) + pr * I;
+                    const TinyVector<Dimension> A_rho_E = (rho_Er + pr) / rhor * rho_ur;
+
+                    double Fn_rho                  = rhor / nb_waves;
+                    TinyVector<Dimension> Fn_rho_u = 1. / nb_waves * rho_ur;
+                    double Fn_rho_E                = rho_Er / 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];
+                    }
+
+                    Fr_rho[node_id][i_wave] += Fn_rho * li_njr;
+                    Fr_rho_u[node_id][i_wave] += li_njr * Fn_rho_u;
+                    Fr_rho_E[node_id][i_wave] += Fn_rho_E * li_njr;
+
+                    sum += li_njr;
+                  }
+
+                  TinyVector<Dimension> njr_prime = txt * njr(cell_id, i_node_cell) - nxn * njr(cell_id, i_node_cell);
+                  double li_njr_prime             = dot(m_lambda_vector[i_wave], njr_prime);
+
+                  if (li_njr_prime > 0) {
+                    const double rhor                  = rho[cell_id](m_xr[node_id]);
+                    const TinyVector<Dimension> rho_ur = rho_u[cell_id](m_xr[node_id]);
+                    const double rho_Er                = rho_E[cell_id](m_xr[node_id]);
+
+                    const double rho_er = rho_Er - 0.5 * (1. / rhor) * dot(rho_ur, rho_ur);
+                    const double pr     = (m_gamma - 1) * rho_er;
+
+                    const TinyVector<Dimension> A_rho   = rho_ur;
+                    const TinyMatrix<Dimension> A_rho_u = 1. / rhor * tensorProduct(rho_ur, rho_ur) + pr * I;
+                    const TinyVector<Dimension> A_rho_E = (rho_Er + pr) / rhor * rho_ur;
+
+                    double Fn_rho                  = rhor / nb_waves;
+                    TinyVector<Dimension> Fn_rho_u = 1. / nb_waves * rho_ur;
+                    double Fn_rho_E                = rho_Er / 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];
+                    }
+
+                    Fr_rho[node_id][i_wave] += Fn_rho * li_njr_prime;
+                    Fr_rho_u[node_id][i_wave] += li_njr_prime * Fn_rho_u;
+                    Fr_rho_E[node_id][i_wave] += Fn_rho_E * li_njr_prime;
+
+                    sum += li_njr_prime;
+                  }
+                }
+                if (sum != 0) {
+                  Fr_rho[node_id][i_wave] /= sum;
+                  Fr_rho_u[node_id][i_wave] = 1. / sum * Fr_rho_u[node_id][i_wave];
+                  Fr_rho_E[node_id][i_wave] /= sum;
+                }
+              }
+            }
+          }
+        },
+        bc_v);
+    }
+  }
+
   void
   apply_eucclhyd_bc(NodeArray<double> Fr_rho,
                     NodeArray<TinyVector<Dimension>> Fr_rho_u,
@@ -1908,6 +2349,144 @@ class EulerKineticSolverMultiDOrder2
       });
   }
 
+  void
+  velocity_limiter(const CellValue<const TinyVector<Dimension>>& u,
+                   DiscreteFunctionDPk<Dimension, const double>& DPk_p,
+                   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) {
+        auto coefficients_p = DPk_p.coefficients(cell_id);
+        TinyVector<Dimension> n;
+        for (size_t dim = 1; dim < coefficients_p.size(); ++dim) {
+          n[dim - 1] = coefficients_p[dim];
+        }
+        if (dot(n, n) != 0) {
+          n = 1. / sqrt(dot(n, n)) * n;
+
+          TinyVector<Dimension> v;
+          bool is_colinear = false;
+          v[0]             = 1;
+          for (size_t dim = 1; dim < Dimension; ++dim) {
+            v[dim] = 0;
+            if (n[0] != 0 and n[dim] == 0) {
+              is_colinear = true;
+            } else {
+              is_colinear = false;
+            }
+          }
+          if (is_colinear) {
+            for (size_t dim = 1; dim < Dimension; ++dim) {
+              if (dim != 1) {
+                v[dim] = 0;
+              } else {
+                v[dim] = 1;
+              }
+            }
+          }
+
+          TinyVector<Dimension> t = v - dot(v, n) * n;
+          t                       = 1. / sqrt(dot(t, t)) * t;
+          TinyVector<Dimension> l = zero;
+
+          if constexpr (Dimension == 3) {
+            l[0] = n[1] * t[2] - n[2] * t[1];
+            l[1] = n[2] * t[0] - n[0] * t[2];
+            l[0] = n[0] * t[1] - n[1] * t[0];
+            l    = 1. / sqrt(dot(l, l)) * l;
+          }
+
+          const double un = dot(u[cell_id], n);
+          const double ut = dot(u[cell_id], t);
+          const double ul = dot(u[cell_id], l);
+
+          double un_min = un;
+          double un_max = un;
+          double ut_min = ut;
+          double ut_max = ut;
+          double ul_min = ul;
+          double ul_max = ul;
+
+          const auto cell_stencil = stencil[cell_id];
+          for (size_t i_cell = 0; i_cell < cell_stencil.size(); ++i_cell) {
+            un_min = std::min(un_min, dot(u[cell_stencil[i_cell]], n));
+            un_max = std::max(un_max, dot(u[cell_stencil[i_cell]], n));
+            ut_min = std::min(ut_min, dot(u[cell_stencil[i_cell]], t));
+            ut_max = std::max(ut_max, dot(u[cell_stencil[i_cell]], t));
+            ul_min = std::min(ul_min, dot(u[cell_stencil[i_cell]], l));
+            ul_max = std::max(ul_max, dot(u[cell_stencil[i_cell]], l));
+          }
+
+          const double eps = 1E-14;
+          double coef_n    = 1;
+          double lambda_n  = 1.;
+          double coef_t    = 1;
+          double lambda_t  = 1.;
+          double coef_l    = 1;
+          double lambda_l  = 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 un_bar_i  = dot(DPk_uh[cell_id](m_xj[cell_i_id]), n);
+            const double Delta_n   = (un_bar_i - un);
+            const double ut_bar_i  = dot(DPk_uh[cell_id](m_xj[cell_i_id]), t);
+            const double Delta_t   = (ut_bar_i - ut);
+            const double ul_bar_i  = dot(DPk_uh[cell_id](m_xj[cell_i_id]), l);
+            const double Delta_l   = (ul_bar_i - ul);
+
+            if (Delta_n > eps) {
+              coef_n = std::min(1., (un_max - un) / Delta_n);
+            } else if (Delta_n < -eps) {
+              coef_n = std::min(1., (un_min - un) / Delta_n);
+            }
+            lambda_n = std::min(lambda_n, coef_n);
+
+            if (Delta_t > eps) {
+              coef_t = std::min(1., (ut_max - ut) / Delta_t);
+            } else if (Delta_t < -eps) {
+              coef_t = std::min(1., (ut_min - ut) / Delta_t);
+            }
+            lambda_t = std::min(lambda_t, coef_t);
+
+            if (Delta_l > eps) {
+              coef_l = std::min(1., (ul_max - ul) / Delta_l);
+            } else if (Delta_l < -eps) {
+              coef_l = std::min(1., (ul_min - ul) / Delta_l);
+            }
+            lambda_l = std::min(lambda_l, coef_l);
+          }
+
+          auto coefficients = DPk_uh.coefficients(cell_id);
+          coefficients[0]   = (1 - lambda_n) * dot(u[cell_id], n) * n + lambda_n * dot(coefficients[0], n) * n +
+                            (1 - lambda_t) * dot(u[cell_id], t) * t + lambda_t * dot(coefficients[0], t) * t +
+                            (1 - lambda_l) * dot(u[cell_id], l) * l + lambda_l * dot(coefficients[0], l) * l;
+          for (size_t i = 1; i < coefficients.size(); ++i) {
+            coefficients[i] = lambda_n * dot(coefficients[i], n) * n + lambda_t * dot(coefficients[i], t) * t +
+                              lambda_l * dot(coefficients[i], l) * l;
+          }
+        }
+      });
+  }
+
+  const CellValue<const double>
+  build_pressure(const CellValue<const double> rho,
+                 const CellValue<const TinyVector<Dimension>> rho_u,
+                 const CellValue<const double> rho_E)
+  {
+    CellValue<double> p{p_mesh->connectivity()};
+    parallel_for(
+      m_mesh.numberOfCells(), PUGS_LAMBDA(const CellId cell_id) {
+        double rho_e = rho_E[cell_id] - 0.5 * (1. / rho[cell_id]) * dot(rho_u[cell_id], rho_u[cell_id]);
+        p[cell_id]   = (m_gamma - 1) * rho_e;
+      });
+    return p;
+  }
+
   std::tuple<std::shared_ptr<const DiscreteFunctionVariant>,
              std::shared_ptr<const DiscreteFunctionVariant>,
              std::shared_ptr<const DiscreteFunctionVariant>>
@@ -1918,8 +2497,10 @@ class EulerKineticSolverMultiDOrder2
     const size_t nb_waves = m_lambda_vector.size();
 
     CellValue<bool> is_inflow_boundary_cell = getInflowBoundaryCells(bc_list);
+    // const auto& cell_to_face_matrix = p_mesh->connectivity().cellToFaceMatrix();
 
     const bool equilibrium_reconstruction = true;
+    const bool limitation                 = true;
 
     const CellValue<double> rho_ex{p_mesh->connectivity()};
     const CellValue<TinyVector<Dimension>> rho_u_ex{p_mesh->connectivity()};
@@ -2045,58 +2626,43 @@ class EulerKineticSolverMultiDOrder2
 
     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));
-
-    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_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_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{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 (equilibrium_reconstruction) {
+      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, TinyVector<Dimension>> rho_u_lim = copy(DPk_rho_u);
+      DiscreteFunctionDPk<Dimension, double> rho_E_lim                = copy(DPk_rho_E);
+
+      if (limitation) {
+        const CellValue<const double> p = build_pressure(rho, rho_u, rho_E);
+        auto reconstruction_pressure =
+          PolynomialReconstruction{reconstruction_descriptor}.build(DiscreteFunctionP0(p_mesh, p));
+        auto DPk_p = reconstruction_pressure[0]->template get<DiscreteFunctionDPk<Dimension, const double>>();
+
+        scalar_limiter_equilibrium(rho, rho_lim);
+        scalar_limiter_equilibrium(rho_E, rho_E_lim);
+        vector_limiter_equilibrium(rho_u, rho_u_lim);
+        // velocity_limiter(rho_u, DPk_p, rho_u_lim);
+      }
+
+      if (m_scheme == 1) {
+        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);
 
@@ -2109,53 +2675,97 @@ class EulerKineticSolverMultiDOrder2
 
         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);
+
+        deltaLFn_rho   = compute_deltaLFn_face<double>(Fl_rho);
+        deltaLFn_rho_u = compute_deltaLFn_face<TinyVector<Dimension>>(Fl_rho_u);
+        deltaLFn_rho_E = compute_deltaLFn_face<double>(Fl_rho_E);
+
+      } else if (m_scheme == 2) {
+        NodeArray<double> Fr_rho{p_mesh->connectivity(), nb_waves};
+        NodeArray<TinyVector<Dimension>> Fr_rho_u{p_mesh->connectivity(), nb_waves};
+        NodeArray<double> Fr_rho_E{p_mesh->connectivity(), nb_waves};
+
+        std::tuple<NodeArray<double>, NodeArray<TinyVector<Dimension>>, NodeArray<double>> vec_Fr(Fr_rho, Fr_rho_u,
+                                                                                                  Fr_rho_E);
+
+        vec_Fr = compute_Flux_glace_equilibrium(rho_lim, rho_u_lim, rho_E_lim, is_boundary_node);
+
+        Fr_rho   = get<0>(vec_Fr);
+        Fr_rho_u = get<1>(vec_Fr);
+        Fr_rho_E = get<2>(vec_Fr);
+
+        apply_glace_bc_equilibrium(Fr_rho, Fr_rho_u, Fr_rho_E, rho_lim, rho_u_lim, rho_E_lim, bc_list);
+
+        synchronize(Fr_rho);
+        synchronize(Fr_rho_u);
+        synchronize(Fr_rho_E);
+
+        deltaLFn_rho   = compute_deltaLFn_node(Fr_rho);
+        deltaLFn_rho_u = compute_deltaLFn_node(Fr_rho_u);
+        deltaLFn_rho_E = compute_deltaLFn_node(Fr_rho_E);
+      }
+    } else {
+      auto reconstruction_distributions =
+        PolynomialReconstruction{reconstruction_descriptor}.build(DiscreteFunctionP0Vector(p_mesh, Fn_rho),
+                                                                  DiscreteFunctionP0Vector(p_mesh, Fn_rho_u),
+                                                                  DiscreteFunctionP0Vector(p_mesh, Fn_rho_E));
+
+      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, TinyVector<Dimension>> Fn_rho_u_lim = copy(DPk_Fn_rho_u);
+      DiscreteFunctionDPkVector<Dimension, double> Fn_rho_E_lim                = copy(DPk_Fn_rho_E);
+
+      if (limitation) {
+        scalar_limiter_distributions(Fn_rho, Fn_rho_lim);
+        scalar_limiter_distributions(Fn_rho_E, Fn_rho_E_lim);
+        vector_limiter_distributions(Fn_rho_u, Fn_rho_u_lim);
       }
-      synchronize(Fl_rho);
-      synchronize(Fl_rho_u);
-      synchronize(Fl_rho_E);
 
-      deltaLFn_rho   = compute_deltaLFn_face<double>(Fl_rho);
-      deltaLFn_rho_u = compute_deltaLFn_face<TinyVector<Dimension>>(Fl_rho_u);
-      deltaLFn_rho_E = compute_deltaLFn_face<double>(Fl_rho_E);
+      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);
 
-    } else if (m_scheme == 2) {
-      NodeArray<double> Fr_rho = compute_Flux_glace<double>(Fn_rho_lim, is_boundary_node);
-      NodeArray<TinyVector<Dimension>> Fr_rho_u =
-        compute_Flux_glace<TinyVector<Dimension>>(Fn_rho_u_lim, is_boundary_node);
-      NodeArray<double> Fr_rho_E = compute_Flux_glace<double>(Fn_rho_E_lim, is_boundary_node);
+        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);
 
-      apply_glace_bc(Fr_rho, Fr_rho_u, Fr_rho_E, Fn_rho_lim, Fn_rho_u_lim, Fn_rho_E_lim, bc_list);
+        synchronize(Fl_rho);
+        synchronize(Fl_rho_u);
+        synchronize(Fl_rho_E);
 
-      synchronize(Fr_rho);
-      synchronize(Fr_rho_u);
-      synchronize(Fr_rho_E);
+        deltaLFn_rho   = compute_deltaLFn_face<double>(Fl_rho);
+        deltaLFn_rho_u = compute_deltaLFn_face<TinyVector<Dimension>>(Fl_rho_u);
+        deltaLFn_rho_E = compute_deltaLFn_face<double>(Fl_rho_E);
 
-      deltaLFn_rho   = compute_deltaLFn_node(Fr_rho);
-      deltaLFn_rho_u = compute_deltaLFn_node(Fr_rho_u);
-      deltaLFn_rho_E = compute_deltaLFn_node(Fr_rho_E);
+      } else if (m_scheme == 2) {
+        NodeArray<double> Fr_rho = compute_Flux_glace<double>(Fn_rho_lim, is_boundary_node);
+        NodeArray<TinyVector<Dimension>> Fr_rho_u =
+          compute_Flux_glace<TinyVector<Dimension>>(Fn_rho_u_lim, is_boundary_node);
+        NodeArray<double>
 
-    } else if (m_scheme == 3) {
-      NodeArray<double> Fr_rho = compute_Flux_eucclhyd<double>(Fn_rho_lim, is_boundary_node);
-      NodeArray<TinyVector<Dimension>> Fr_rho_u =
-        compute_Flux_eucclhyd<TinyVector<Dimension>>(Fn_rho_u_lim, is_boundary_node);
-      NodeArray<double> Fr_rho_E = compute_Flux_eucclhyd<double>(Fn_rho_E_lim, is_boundary_node);
+          Fr_rho_E = compute_Flux_glace<double>(Fn_rho_E_lim, is_boundary_node);
 
-      apply_eucclhyd_bc(Fr_rho, Fr_rho_u, Fr_rho_E, Fn_rho_lim, Fn_rho_u_lim, Fn_rho_E_lim, bc_list);
+        apply_glace_bc(Fr_rho, Fr_rho_u, Fr_rho_E, Fn_rho_lim, Fn_rho_u_lim, Fn_rho_E_lim, bc_list);
 
-      synchronize(Fr_rho);
-      synchronize(Fr_rho_u);
-      synchronize(Fr_rho_E);
+        synchronize(Fr_rho);
+        synchronize(Fr_rho_u);
+        synchronize(Fr_rho_E);
 
-      deltaLFn_rho   = compute_deltaLFn_node(Fr_rho);
-      deltaLFn_rho_u = compute_deltaLFn_node(Fr_rho_u);
-      deltaLFn_rho_E = compute_deltaLFn_node(Fr_rho_E);
+        deltaLFn_rho   = compute_deltaLFn_node(Fr_rho);
+        deltaLFn_rho_u = compute_deltaLFn_node(Fr_rho_u);
+        deltaLFn_rho_E = compute_deltaLFn_node(Fr_rho_E);
+      }
     }
 
     const CellValue<const TinyVector<Dimension>> A_rho   = getA_rho(rho_u);
@@ -2167,64 +2777,48 @@ 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_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_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);
-
-      scalar_limiter_distributions(Fnp1_rho, Fnp1_rho_lim);
-      scalar_limiter_distributions(Fnp1_rho_E, Fnp1_rho_E_lim);
-
-      DiscreteFunctionDPkVector<Dimension, TinyVector<Dimension>> Fnp1_rho_u_lim = copy(DPk_Fnp1_rho_u);
-
-      vector_limiter_distributions(Fnp1_rho_u, Fnp1_rho_u_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);
 
-      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);
-
-      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{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 (equilibrium_reconstruction) {
+        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, TinyVector<Dimension>> rho_u_np1_lim = copy(DPk_rho_u_np1);
+        DiscreteFunctionDPk<Dimension, double> rho_E_np1_lim                = copy(DPk_rho_E_np1);
+
+        if (limitation) {
+          const CellValue<const double> p_np1 = build_pressure(rho_np1, rho_u_np1, rho_E_np1);
+          auto reconstruction_pressure_np1 =
+            PolynomialReconstruction{reconstruction_descriptor}.build(DiscreteFunctionP0(p_mesh, p_np1));
+          auto DPk_p_np1 = reconstruction_pressure_np1[0]->template get<DiscreteFunctionDPk<Dimension, const double>>();
+
+          scalar_limiter_equilibrium(rho_np1, rho_np1_lim);
+          scalar_limiter_equilibrium(rho_E_np1, rho_E_np1_lim);
+          vector_limiter_equilibrium(rho_u_np1, rho_u_np1_lim);
+          // velocity_limiter(rho_u_np1, DPk_p_np1, rho_u_np1_lim);
+        }
+
+        if (m_scheme == 1) {
+          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);
@@ -2236,54 +2830,97 @@ class EulerKineticSolverMultiDOrder2
           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_equilibrium(Fl_rho_np1, Fl_rho_u_np1, Fl_rho_E_np1, rho_np1_lim, rho_u_np1_lim, rho_E_np1_lim,
+                                    Fnp1_rho, Fnp1_rho_u, Fnp1_rho_E, bc_list);
+
+          synchronize(Fl_rho_np1);
+          synchronize(Fl_rho_u_np1);
+          synchronize(Fl_rho_E_np1);
+
+          deltaLFnp1_rho   = compute_deltaLFn_face<double>(Fl_rho_np1);
+          deltaLFnp1_rho_u = compute_deltaLFn_face<TinyVector<Dimension>>(Fl_rho_u_np1);
+          deltaLFnp1_rho_E = compute_deltaLFn_face<double>(Fl_rho_E_np1);
+        } else if (m_scheme == 2) {
+          NodeArray<double> Fr_rho_np1{p_mesh->connectivity(), nb_waves};
+          NodeArray<TinyVector<Dimension>> Fr_rho_u_np1{p_mesh->connectivity(), nb_waves};
+          NodeArray<double> Fr_rho_E_np1{p_mesh->connectivity(), nb_waves};
+
+          std::tuple<NodeArray<double>, NodeArray<TinyVector<Dimension>>, NodeArray<double>> vec_Fr_np1(Fr_rho_np1,
+                                                                                                        Fr_rho_u_np1,
+                                                                                                        Fr_rho_E_np1);
+
+          vec_Fr_np1 = compute_Flux_glace_equilibrium(rho_np1_lim, rho_u_np1_lim, rho_E_np1_lim, is_boundary_node);
+
+          Fr_rho_np1   = get<0>(vec_Fr_np1);
+          Fr_rho_u_np1 = get<1>(vec_Fr_np1);
+          Fr_rho_E_np1 = get<2>(vec_Fr_np1);
 
-          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);
+          apply_glace_bc_equilibrium(Fr_rho_np1, Fr_rho_u_np1, Fr_rho_E_np1, rho_np1_lim, rho_u_np1_lim, rho_E_np1_lim,
+                                     bc_list);
+
+          synchronize(Fr_rho_np1);
+          synchronize(Fr_rho_u_np1);
+          synchronize(Fr_rho_E_np1);
+
+          deltaLFnp1_rho   = compute_deltaLFn_node<double>(Fr_rho_np1);
+          deltaLFnp1_rho_u = compute_deltaLFn_node<TinyVector<Dimension>>(Fr_rho_u_np1);
+          deltaLFnp1_rho_E = compute_deltaLFn_node<double>(Fr_rho_E_np1);
+        }
+      } else {
+        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_distributions_np1[2]->template get<DiscreteFunctionDPkVector<Dimension, const double>>();
+
+        DiscreteFunctionDPkVector<Dimension, double> Fnp1_rho_lim                  = copy(DPk_Fnp1_rho);
+        DiscreteFunctionDPkVector<Dimension, TinyVector<Dimension>> Fnp1_rho_u_lim = copy(DPk_Fnp1_rho_u);
+        DiscreteFunctionDPkVector<Dimension, double> Fnp1_rho_E_lim                = copy(DPk_Fnp1_rho_E);
+
+        if (limitation) {
+          scalar_limiter_distributions(Fnp1_rho, Fnp1_rho_lim);
+          scalar_limiter_distributions(Fnp1_rho_E, Fnp1_rho_E_lim);
+          vector_limiter_distributions(Fnp1_rho_u, Fnp1_rho_u_lim);
         }
-        synchronize(Fl_rho_np1);
-        synchronize(Fl_rho_u_np1);
-        synchronize(Fl_rho_E_np1);
 
-        deltaLFnp1_rho   = compute_deltaLFn_face<double>(Fl_rho_np1);
-        deltaLFnp1_rho_u = compute_deltaLFn_face<TinyVector<Dimension>>(Fl_rho_u_np1);
-        deltaLFnp1_rho_E = compute_deltaLFn_face<double>(Fl_rho_E_np1);
-      } else if (m_scheme == 2) {
-        NodeArray<double> Fr_rho_np1 = compute_Flux_glace<double>(Fnp1_rho_lim, is_boundary_node);
-        NodeArray<TinyVector<Dimension>> Fr_rho_u_np1 =
-          compute_Flux_glace<TinyVector<Dimension>>(Fnp1_rho_u_lim, is_boundary_node);
-        NodeArray<double> Fr_rho_E_np1 = compute_Flux_glace<double>(Fnp1_rho_E_lim, is_boundary_node);
-
-        apply_glace_bc(Fr_rho_np1, Fr_rho_u_np1, Fr_rho_E_np1, Fnp1_rho_lim, Fnp1_rho_u_lim, Fnp1_rho_E_lim, bc_list);
-
-        synchronize(Fr_rho_np1);
-        synchronize(Fr_rho_u_np1);
-        synchronize(Fr_rho_E_np1);
-
-        deltaLFnp1_rho   = compute_deltaLFn_node<double>(Fr_rho_np1);
-        deltaLFnp1_rho_u = compute_deltaLFn_node<TinyVector<Dimension>>(Fr_rho_u_np1);
-        deltaLFnp1_rho_E = compute_deltaLFn_node<double>(Fr_rho_E_np1);
-      } else if (m_scheme == 3) {
-        NodeArray<double> Fr_rho_np1 = compute_Flux_eucclhyd<double>(Fnp1_rho_lim, is_boundary_node);
-        NodeArray<TinyVector<Dimension>> Fr_rho_u_np1 =
-          compute_Flux_eucclhyd<TinyVector<Dimension>>(Fnp1_rho_u_lim, is_boundary_node);
-        NodeArray<double> Fr_rho_E_np1 = compute_Flux_eucclhyd<double>(Fnp1_rho_E_lim, is_boundary_node);
-
-        apply_eucclhyd_bc(Fr_rho_np1, Fr_rho_u_np1, Fr_rho_E_np1, Fnp1_rho_lim, Fnp1_rho_u_lim, Fnp1_rho_E_lim,
-                          bc_list);
-
-        synchronize(Fr_rho_np1);
-        synchronize(Fr_rho_u_np1);
-        synchronize(Fr_rho_E_np1);
-
-        deltaLFnp1_rho   = compute_deltaLFn_node<double>(Fr_rho_np1);
-        deltaLFnp1_rho_u = compute_deltaLFn_node<TinyVector<Dimension>>(Fr_rho_u_np1);
-        deltaLFnp1_rho_E = compute_deltaLFn_node<double>(Fr_rho_E_np1);
+        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, Fnp1_rho,
+                        Fnp1_rho_u, Fnp1_rho_E, bc_list);
+
+          synchronize(Fl_rho_np1);
+          synchronize(Fl_rho_u_np1);
+          synchronize(Fl_rho_E_np1);
+
+          deltaLFnp1_rho   = compute_deltaLFn_face<double>(Fl_rho_np1);
+          deltaLFnp1_rho_u = compute_deltaLFn_face<TinyVector<Dimension>>(Fl_rho_u_np1);
+          deltaLFnp1_rho_E = compute_deltaLFn_face<double>(Fl_rho_E_np1);
+        } else if (m_scheme == 2) {
+          NodeArray<double> Fr_rho_np1 = compute_Flux_glace<double>(Fnp1_rho_lim, is_boundary_node);
+          NodeArray<TinyVector<Dimension>> Fr_rho_u_np1 =
+            compute_Flux_glace<TinyVector<Dimension>>(Fnp1_rho_u_lim, is_boundary_node);
+          NodeArray<double> Fr_rho_E_np1 = compute_Flux_glace<double>(Fnp1_rho_E_lim, is_boundary_node);
+
+          apply_glace_bc(Fr_rho_np1, Fr_rho_u_np1, Fr_rho_E_np1, Fnp1_rho_lim, Fnp1_rho_u_lim, Fnp1_rho_E_lim, bc_list);
+
+          synchronize(Fr_rho_np1);
+          synchronize(Fr_rho_u_np1);
+          synchronize(Fr_rho_E_np1);
+
+          deltaLFnp1_rho   = compute_deltaLFn_node<double>(Fr_rho_np1);
+          deltaLFnp1_rho_u = compute_deltaLFn_node<TinyVector<Dimension>>(Fr_rho_u_np1);
+          deltaLFnp1_rho_E = compute_deltaLFn_node<double>(Fr_rho_E_np1);
+        }
       }
 
       rho_np1   = compute_PFnp1<double>(Fn_rho, deltaLFn_rho, deltaLFnp1_rho);
@@ -2617,6 +3254,9 @@ eulerKineticSolverMultiDOrder2(
   if (SpaceOrder > 2 or TimeOrder > 2) {
     throw NotImplementedError("Euler kinetic solver in Multi D not implemented at order > 2");
   }
+  if (scheme > 2) {
+    throw NotImplementedError("Eucclhyd scheme not implemented");
+  }
 
   return std::visit(
     [&](auto&& p_mesh)