From 08eda47da1dcb4b30c8439e9cccc77e08d7bec70 Mon Sep 17 00:00:00 2001
From: Axelle <axelle.drouard@orange.fr>
Date: Thu, 30 Jan 2025 15:12:25 +0100
Subject: [PATCH] Add nodal schemes + BC to 2nd order

---
 src/language/modules/KineticSchemeModule.cpp  |   20 +-
 src/scheme/EulerKineticSolverMultiD.cpp       |  280 +++--
 src/scheme/EulerKineticSolverMultiD.hpp       |    1 +
 src/scheme/EulerKineticSolverMultiDOrder2.cpp | 1107 ++++++++++++++---
 src/scheme/EulerKineticSolverMultiDOrder2.hpp |    1 +
 5 files changed, 1134 insertions(+), 275 deletions(-)

diff --git a/src/language/modules/KineticSchemeModule.cpp b/src/language/modules/KineticSchemeModule.cpp
index a111e0522..bfe0d148e 100644
--- a/src/language/modules/KineticSchemeModule.cpp
+++ b/src/language/modules/KineticSchemeModule.cpp
@@ -564,7 +564,7 @@ KineticSchemeModule::KineticSchemeModule()
                             std::function(
 
                               [](const double& dt, const double& gamma, const std::vector<TinyVector<2>>& lambda_vector,
-                                 const double& eps, const size_t& SpaceOrder,
+                                 const double& eps, const size_t& SpaceOrder, const size_t& scheme,
                                  const std::shared_ptr<const DiscreteFunctionVariant>& F_rho,
                                  const std::shared_ptr<const DiscreteFunctionVariant>& F_rho_u,
                                  const std::shared_ptr<const DiscreteFunctionVariant>& F_E,
@@ -572,8 +572,8 @@ KineticSchemeModule::KineticSchemeModule()
                                    bc_descriptor_list) -> std::tuple<std::shared_ptr<const DiscreteFunctionVariant>,
                                                                      std::shared_ptr<const DiscreteFunctionVariant>,
                                                                      std::shared_ptr<const DiscreteFunctionVariant>> {
-                                return eulerKineticSolverMultiD(dt, gamma, lambda_vector, eps, SpaceOrder, F_rho,
-                                                                F_rho_u, F_E, bc_descriptor_list);
+                                return eulerKineticSolverMultiD(dt, gamma, lambda_vector, eps, SpaceOrder, scheme,
+                                                                F_rho, F_rho_u, F_E, bc_descriptor_list);
                               }
 
                               ));
@@ -582,7 +582,7 @@ KineticSchemeModule::KineticSchemeModule()
                             std::function(
 
                               [](const double& dt, const double& gamma, const std::vector<TinyVector<3>>& lambda_vector,
-                                 const double& eps, const size_t& SpaceOrder,
+                                 const double& eps, const size_t& SpaceOrder, const size_t& scheme,
                                  const std::shared_ptr<const DiscreteFunctionVariant>& F_rho,
                                  const std::shared_ptr<const DiscreteFunctionVariant>& F_rho_u,
                                  const std::shared_ptr<const DiscreteFunctionVariant>& F_E,
@@ -590,8 +590,8 @@ KineticSchemeModule::KineticSchemeModule()
                                    bc_descriptor_list) -> std::tuple<std::shared_ptr<const DiscreteFunctionVariant>,
                                                                      std::shared_ptr<const DiscreteFunctionVariant>,
                                                                      std::shared_ptr<const DiscreteFunctionVariant>> {
-                                return eulerKineticSolverMultiD(dt, gamma, lambda_vector, eps, SpaceOrder, F_rho,
-                                                                F_rho_u, F_E, bc_descriptor_list);
+                                return eulerKineticSolverMultiD(dt, gamma, lambda_vector, eps, SpaceOrder, scheme,
+                                                                F_rho, F_rho_u, F_E, bc_descriptor_list);
                               }
 
                               ));
@@ -650,7 +650,7 @@ KineticSchemeModule::KineticSchemeModule()
 
                               [](const double& dt, const double& gamma, const std::vector<TinyVector<2>>& lambda_vector,
                                  const double& eps, const size_t& SpaceOrder, const size_t& TimeOrder,
-                                 const std::shared_ptr<const DiscreteFunctionVariant>& F_rho,
+                                 const size_t& scheme, const std::shared_ptr<const DiscreteFunctionVariant>& F_rho,
                                  const std::shared_ptr<const DiscreteFunctionVariant>& F_rho_u,
                                  const std::shared_ptr<const DiscreteFunctionVariant>& F_E,
                                  const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>&
@@ -658,7 +658,7 @@ KineticSchemeModule::KineticSchemeModule()
                                                                      std::shared_ptr<const DiscreteFunctionVariant>,
                                                                      std::shared_ptr<const DiscreteFunctionVariant>> {
                                 return eulerKineticSolverMultiDOrder2(dt, gamma, lambda_vector, eps, SpaceOrder,
-                                                                      TimeOrder, F_rho, F_rho_u, F_E,
+                                                                      TimeOrder, scheme, F_rho, F_rho_u, F_E,
                                                                       bc_descriptor_list);
                               }
 
@@ -669,7 +669,7 @@ KineticSchemeModule::KineticSchemeModule()
 
                               [](const double& dt, const double& gamma, const std::vector<TinyVector<3>>& lambda_vector,
                                  const double& eps, const size_t& SpaceOrder, const size_t& TimeOrder,
-                                 const std::shared_ptr<const DiscreteFunctionVariant>& F_rho,
+                                 const size_t& scheme, const std::shared_ptr<const DiscreteFunctionVariant>& F_rho,
                                  const std::shared_ptr<const DiscreteFunctionVariant>& F_rho_u,
                                  const std::shared_ptr<const DiscreteFunctionVariant>& F_E,
                                  const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>&
@@ -677,7 +677,7 @@ KineticSchemeModule::KineticSchemeModule()
                                                                      std::shared_ptr<const DiscreteFunctionVariant>,
                                                                      std::shared_ptr<const DiscreteFunctionVariant>> {
                                 return eulerKineticSolverMultiDOrder2(dt, gamma, lambda_vector, eps, SpaceOrder,
-                                                                      TimeOrder, F_rho, F_rho_u, F_E,
+                                                                      TimeOrder, scheme, F_rho, F_rho_u, F_E,
                                                                       bc_descriptor_list);
                               }
 
diff --git a/src/scheme/EulerKineticSolverMultiD.cpp b/src/scheme/EulerKineticSolverMultiD.cpp
index af1c1eb63..010ddd67e 100644
--- a/src/scheme/EulerKineticSolverMultiD.cpp
+++ b/src/scheme/EulerKineticSolverMultiD.cpp
@@ -228,6 +228,7 @@ class EulerKineticSolverMultiD
   const double m_eps;
   const std::vector<TinyVector<Dimension>>& m_lambda_vector;
   const size_t m_SpaceOrder;
+  const size_t m_scheme;
   const CellArray<const double> m_Fn_rho;
   const CellArray<const TinyVector<Dimension>> m_Fn_rho_u;
   const CellArray<const double> m_Fn_rho_E;
@@ -868,14 +869,11 @@ class EulerKineticSolverMultiD
               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);
-              double sum_Cjr                            = 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];
 
                 nr[node_id] += Cjr(cell_id, i_node_cell);
-
-                sum_Cjr += sqrt(dot(Cjr(cell_id, i_node_cell), 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]);
@@ -956,14 +954,11 @@ class EulerKineticSolverMultiD
               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);
 
-              double sum_Cjr = 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];
 
                 nr[node_id] += Cjr(cell_id, i_node_cell);
-
-                sum_Cjr += sqrt(dot(Cjr(cell_id, i_node_cell), 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]);
@@ -1031,6 +1026,66 @@ class EulerKineticSolverMultiD
                 }
               }
             }
+          } 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) {
+                    Fr_rho[node_id][i_wave] += Fn_rho[cell_id][i_wave] * li_njr;
+                    Fr_rho_u[node_id][i_wave] += li_njr * Fn_rho_u[cell_id][i_wave];
+                    Fr_rho_E[node_id][i_wave] += Fn_rho_E[cell_id][i_wave] * 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) {
+                    Fr_rho[node_id][i_wave] += Fn_rho[cell_id][i_wave] * li_njr_prime;
+                    Fr_rho_u[node_id][i_wave] += li_njr_prime * Fn_rho_u[cell_id][i_wave];
+                    Fr_rho_E[node_id][i_wave] += Fn_rho_E[cell_id][i_wave] * 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);
@@ -1093,7 +1148,7 @@ class EulerKineticSolverMultiD
                 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] += njr(cell_id, i_node_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]);
@@ -1187,14 +1242,11 @@ class EulerKineticSolverMultiD
               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);
 
-              double sum_Cjr = 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];
 
                 nr[node_id] += Cjr(cell_id, i_node_cell);
-
-                sum_Cjr += sqrt(dot(Cjr(cell_id, i_node_cell), 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]);
@@ -1275,6 +1327,80 @@ class EulerKineticSolverMultiD
                 }
               }
             }
+          } 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) {
+                const auto& node_to_face                  = node_to_face_matrix[node_id];
+                const auto& node_local_number_in_its_face = node_local_numbers_in_their_faces.itemArray(node_id);
+
+                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_face = 0; i_face < node_to_face.size(); ++i_face) {
+                  const FaceId face_id                      = node_to_face[i_face];
+                  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 auto& face_local_number_in_its_cell = face_local_numbers_in_their_cells.itemArray(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];
+
+                    const double sign            = face_cell_is_reversed(cell_id, i_face_cell) ? -1 : 1;
+                    TinyVector<Dimension> n_face = sign * nlr(face_id, i_node_face);
+
+                    double li_nlr = dot(m_lambda_vector[i_wave], n_face);
+
+                    if (li_nlr > 0) {
+                      Fr_rho[node_id][i_wave] += Fn_rho[cell_id][i_wave] * li_nlr;
+                      Fr_rho_u[node_id][i_wave] += li_nlr * Fn_rho_u[cell_id][i_wave];
+                      Fr_rho_E[node_id][i_wave] += Fn_rho_E[cell_id][i_wave] * li_nlr;
+
+                      sum += li_nlr;
+                    }
+
+                    TinyVector<Dimension> nlr_prime = txt * n_face - nxn * n_face;
+                    double li_nlr_prime             = dot(m_lambda_vector[i_wave], nlr_prime);
+
+                    if (li_nlr_prime > 0) {
+                      Fr_rho[node_id][i_wave] += Fn_rho[cell_id][i_wave] * li_nlr_prime;
+                      Fr_rho_u[node_id][i_wave] += li_nlr_prime * Fn_rho_u[cell_id][i_wave];
+                      Fr_rho_E[node_id][i_wave] += Fn_rho_E[cell_id][i_wave] * li_nlr_prime;
+
+                      sum += li_nlr_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);
@@ -1317,66 +1443,6 @@ class EulerKineticSolverMultiD
     }
   }
 
-  template <typename T>
-  CellArray<T>
-  compute_deltaLFn_eucclhyd(FaceArrayPerNode<T> Flr)
-  {
-    if constexpr (Dimension > 1) {
-      const size_t nb_waves = m_lambda_vector.size();
-      CellArray<T> deltaLFn(p_mesh->connectivity(), nb_waves);
-
-      if constexpr (is_tiny_vector_v<T>) {
-        deltaLFn.fill(zero);
-      } else {
-        deltaLFn.fill(0.);
-      }
-
-      const NodeValuePerFace<const TinyVector<Dimension>> Nlr = MeshDataManager::instance().getMeshData(m_mesh).Nlr();
-
-      const auto& cell_to_face_matrix   = p_mesh->connectivity().cellToFaceMatrix();
-      const auto& face_to_node_matrix   = p_mesh->connectivity().faceToNodeMatrix();
-      const auto& node_to_face_matrix   = p_mesh->connectivity().nodeToFaceMatrix();
-      const auto& face_cell_is_reversed = p_mesh->connectivity().cellFaceIsReversed();
-
-      parallel_for(
-        p_mesh->numberOfCells(), PUGS_LAMBDA(CellId cell_id) {
-          const auto& cell_to_face = cell_to_face_matrix[cell_id];
-
-          for (size_t i_wave = 0; i_wave < nb_waves; ++i_wave) {
-            for (size_t i_face_cell = 0; i_face_cell < cell_to_face.size(); ++i_face_cell) {
-              const FaceId face_id     = cell_to_face[i_face_cell];
-              const auto& face_to_node = face_to_node_matrix[face_id];
-
-              for (size_t i_node_face = 0; i_node_face < face_to_node.size(); ++i_node_face) {
-                const NodeId node_id     = face_to_node[i_node_face];
-                const auto& node_to_face = node_to_face_matrix[node_id];
-
-                auto local_face_number_in_node = [&](FaceId face_number) {
-                  for (size_t i_face = 0; i_face < node_to_face.size(); ++i_face) {
-                    if (face_number == node_to_face[i_face]) {
-                      return i_face;
-                    }
-                  }
-                  return std::numeric_limits<size_t>::max();
-                };
-
-                const size_t i_face_node = local_face_number_in_node(face_id);
-
-                const double sign = face_cell_is_reversed(cell_id, i_face_cell) ? -1 : 1;
-
-                const double li_Nlr = sign * dot(m_lambda_vector[i_wave], Nlr(face_id, i_node_face));
-
-                deltaLFn[cell_id][i_wave] += li_Nlr * Flr[node_id][i_face_node][i_wave];
-              }
-            }
-          }
-        });
-      return deltaLFn;
-    } else {
-      throw NormalError("Eucclhyd not defined in 1d");
-    }
-  }
-
   template <typename T>
   CellArray<T>
   compute_deltaLFn_face(FaceArray<T> Fl)
@@ -1523,35 +1589,56 @@ class EulerKineticSolverMultiD
     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);
 
-    NodeArray<double> Fr_rho                  = compute_Flux1_glace<double>(Fn_rho, is_boundary_node);
-    NodeArray<TinyVector<Dimension>> Fr_rho_u = compute_Flux1_glace<TinyVector<Dimension>>(Fn_rho_u, is_boundary_node);
-    NodeArray<double> Fr_rho_E                = compute_Flux1_glace<double>(Fn_rho_E, is_boundary_node);
+    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_Flux1_face<double>(Fn_rho);
+      FaceArray<TinyVector<Dimension>> Fl_rho_u = compute_Flux1_face<TinyVector<Dimension>>(Fn_rho_u);
+      FaceArray<double> Fl_rho_E                = compute_Flux1_face<double>(Fn_rho_E);
 
-    // NodeArray<double> Fr_rho = compute_Flux1_eucclhyd<double>(Fn_rho, is_boundary_node);
-    // NodeArray<TinyVector<Dimension>> Fr_rho_u =
-    //   compute_Flux1_eucclhyd<TinyVector<Dimension>>(Fn_rho_u, is_boundary_node);
-    // NodeArray<double> Fr_rho_E = compute_Flux1_eucclhyd<double>(Fn_rho_E, is_boundary_node);
+      apply_face_bc(Fl_rho, Fl_rho_u, Fl_rho_E, rho, rho_u, rho_E, Fn_rho, Fn_rho_u, Fn_rho_E, bc_list);
 
-    // FaceArray<double> Fl_rho                  = compute_Flux1_face<double>(Fn_rho);
-    // FaceArray<TinyVector<Dimension>> Fl_rho_u = compute_Flux1_face<TinyVector<Dimension>>(Fn_rho_u);
-    // FaceArray<double> Fl_rho_E                = compute_Flux1_face<double>(Fn_rho_E);
+      synchronize(Fl_rho);
+      synchronize(Fl_rho_u);
+      synchronize(Fl_rho_E);
 
-    //    apply_face_bc(Fl_rho, Fl_rho_u, Fl_rho_E, rho, rho_u, rho_E, Fn_rho, Fn_rho_u, Fn_rho_E, bc_list);
-    apply_glace_bc(Fr_rho, Fr_rho_u, Fr_rho_E, rho, rho_u, rho_E, Fn_rho, Fn_rho_u, Fn_rho_E, bc_list);
-    // apply_eucclhyd_bc(Fr_rho, Fr_rho_u, Fr_rho_E, rho, rho_u, rho_E, Fn_rho, Fn_rho_u, Fn_rho_E, bc_list);
+      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);
 
-    synchronize(Fr_rho);
-    synchronize(Fr_rho_u);
-    synchronize(Fr_rho_E);
+    } else if (m_scheme == 2) {
+      NodeArray<double> Fr_rho = compute_Flux1_glace<double>(Fn_rho, is_boundary_node);
+      NodeArray<TinyVector<Dimension>> Fr_rho_u =
+        compute_Flux1_glace<TinyVector<Dimension>>(Fn_rho_u, is_boundary_node);
+      NodeArray<double> Fr_rho_E = compute_Flux1_glace<double>(Fn_rho_E, is_boundary_node);
 
-    const CellArray<const double> deltaLFn_rho                  = compute_deltaLFn_node(Fr_rho);
-    const CellArray<const TinyVector<Dimension>> deltaLFn_rho_u = compute_deltaLFn_node(Fr_rho_u);
-    const CellArray<const double> deltaLFn_rho_E                = compute_deltaLFn_node(Fr_rho_E);
+      apply_glace_bc(Fr_rho, Fr_rho_u, Fr_rho_E, rho, rho_u, rho_E, Fn_rho, Fn_rho_u, Fn_rho_E, bc_list);
+      synchronize(Fr_rho);
+      synchronize(Fr_rho_u);
+      synchronize(Fr_rho_E);
 
-    // const CellArray<const double> deltaLFn_rho = compute_deltaLFn_face<double>(Fl_rho);
-    // const CellArray<const TinyVector<Dimension>> deltaLFn_rho_u =
-    //   compute_deltaLFn_face<TinyVector<Dimension>>(Fl_rho_u);
-    // const CellArray<const double> 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 == 3) {
+      NodeArray<double> Fr_rho = compute_Flux1_eucclhyd<double>(Fn_rho, is_boundary_node);
+      NodeArray<TinyVector<Dimension>> Fr_rho_u =
+        compute_Flux1_eucclhyd<TinyVector<Dimension>>(Fn_rho_u, is_boundary_node);
+      NodeArray<double> Fr_rho_E = compute_Flux1_eucclhyd<double>(Fn_rho_E, is_boundary_node);
+
+      apply_eucclhyd_bc(Fr_rho, Fr_rho_u, Fr_rho_E, rho, rho_u, rho_E, Fn_rho, Fn_rho_u, Fn_rho_E, 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);
+    }
 
     const CellValue<const TinyVector<Dimension>> A_rho   = getA_rho(rho_u);
     const CellValue<const TinyMatrix<Dimension>> A_rho_u = getA_rho_u(rho, rho_u, rho_E);
@@ -1609,6 +1696,7 @@ class EulerKineticSolverMultiD
                            const double& gamma,
                            const std::vector<TinyVector<MeshType::Dimension>>& lambda_vector,
                            const size_t& SpaceOrder,
+                           const size_t& scheme,
                            const CellArray<const double>& Fn_rho,
                            const CellArray<const TinyVector<MeshType::Dimension>>& Fn_rho_u,
                            const CellArray<const double>& Fn_rho_E)
@@ -1616,6 +1704,7 @@ class EulerKineticSolverMultiD
       m_eps{eps},
       m_lambda_vector{lambda_vector},
       m_SpaceOrder{SpaceOrder},
+      m_scheme{scheme},
       m_Fn_rho{Fn_rho},
       m_Fn_rho_u{Fn_rho_u},
       m_Fn_rho_E{Fn_rho_E},
@@ -1834,6 +1923,7 @@ eulerKineticSolverMultiD(const double& dt,
                          const std::vector<TinyVector<Dimension>>& lambda_vector,
                          const double& eps,
                          const size_t& SpaceOrder,
+                         const size_t& scheme,
                          const std::shared_ptr<const DiscreteFunctionVariant>& F_rho,
                          const std::shared_ptr<const DiscreteFunctionVariant>& F_rho_u,
                          const std::shared_ptr<const DiscreteFunctionVariant>& F_rho_E,
@@ -1857,7 +1947,7 @@ eulerKineticSolverMultiD(const double& dt,
         auto Fn_rho_u = F_rho_u->get<DiscreteFunctionP0Vector<const TinyVector<MeshType::Dimension>>>();
         auto Fn_rho_E = F_rho_E->get<DiscreteFunctionP0Vector<const double>>();
 
-        EulerKineticSolverMultiD solver(p_mesh, dt, eps, gamma, lambda_vector, SpaceOrder, Fn_rho.cellArrays(),
+        EulerKineticSolverMultiD solver(p_mesh, dt, eps, gamma, lambda_vector, SpaceOrder, scheme, Fn_rho.cellArrays(),
                                         Fn_rho_u.cellArrays(), Fn_rho_E.cellArrays());
 
         return solver.apply(bc_descriptor_list);
@@ -1876,6 +1966,7 @@ eulerKineticSolverMultiD(const double& dt,
                          const std::vector<TinyVector<2>>& lambda_vector,
                          const double& eps,
                          const size_t& SpaceOrder,
+                         const size_t& scheme,
                          const std::shared_ptr<const DiscreteFunctionVariant>& F_rho,
                          const std::shared_ptr<const DiscreteFunctionVariant>& F_rho_u,
                          const std::shared_ptr<const DiscreteFunctionVariant>& F_rho_E,
@@ -1889,6 +1980,7 @@ eulerKineticSolverMultiD(const double& dt,
                          const std::vector<TinyVector<3>>& lambda_vector,
                          const double& eps,
                          const size_t& SpaceOrder,
+                         const size_t& scheme,
                          const std::shared_ptr<const DiscreteFunctionVariant>& F_rho,
                          const std::shared_ptr<const DiscreteFunctionVariant>& F_rho_u,
                          const std::shared_ptr<const DiscreteFunctionVariant>& F_rho_E,
diff --git a/src/scheme/EulerKineticSolverMultiD.hpp b/src/scheme/EulerKineticSolverMultiD.hpp
index 987f10f8a..4053facd0 100644
--- a/src/scheme/EulerKineticSolverMultiD.hpp
+++ b/src/scheme/EulerKineticSolverMultiD.hpp
@@ -32,6 +32,7 @@ eulerKineticSolverMultiD(const double& dt,
                          const std::vector<TinyVector<Dimension>>& lambda_vector,
                          const double& eps,
                          const size_t& SpaceOrder,
+                         const size_t& scheme,
                          const std::shared_ptr<const DiscreteFunctionVariant>& F_rho,
                          const std::shared_ptr<const DiscreteFunctionVariant>& F_rho_u,
                          const std::shared_ptr<const DiscreteFunctionVariant>& F_rho_E,
diff --git a/src/scheme/EulerKineticSolverMultiDOrder2.cpp b/src/scheme/EulerKineticSolverMultiDOrder2.cpp
index 82776378c..02c995fdd 100644
--- a/src/scheme/EulerKineticSolverMultiDOrder2.cpp
+++ b/src/scheme/EulerKineticSolverMultiDOrder2.cpp
@@ -50,6 +50,7 @@ class EulerKineticSolverMultiDOrder2
   const std::vector<TinyVector<Dimension>>& m_lambda_vector;
   const size_t m_SpaceOrder;
   const size_t m_TimeOrder;
+  const size_t m_scheme;
   const CellArray<const double> m_Fn_rho;
   const CellArray<const TinyVector<Dimension>> m_Fn_rho_u;
   const CellArray<const double> m_Fn_rho_E;
@@ -451,7 +452,7 @@ class EulerKineticSolverMultiDOrder2
 
   template <typename T>
   FaceArray<T>
-  compute_Flux1_face(const DiscreteFunctionDPkVector<Dimension, const T> DPk_Fn)
+  compute_Flux_face(const DiscreteFunctionDPkVector<Dimension, const T> DPk_Fn)
   {
     const size_t nb_waves = m_lambda_vector.size();
     FaceArray<T> Fl(m_mesh.connectivity(), nb_waves);
@@ -490,17 +491,125 @@ class EulerKineticSolverMultiDOrder2
     return Fl;
   }
 
+  template <typename T>
+  NodeArray<T>
+  compute_Flux_glace(const DiscreteFunctionDPkVector<Dimension, const T> DPk_Fn, 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<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) {
+          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) {
+                  Fr[node_id][i_wave] += li_njr * DPk_Fn(cell_id, i_wave)(m_xr[node_id]);
+                  sum_li_njr += li_njr;
+                }
+              }
+              if (sum_li_njr != 0) {
+                Fr[node_id][i_wave] = 1. / sum_li_njr * Fr[node_id][i_wave];
+              }
+            }
+          }
+        });
+
+      return Fr;
+    } 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)
+  {
+    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();
+      NodeArray<T> Fr(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();
+
+      if constexpr (is_tiny_vector_v<T>) {
+        Fr.fill(zero);
+      } else {
+        Fr.fill(0.);
+      }
+
+      parallel_for(
+        p_mesh->numberOfNodes(), PUGS_LAMBDA(NodeId node_id) {
+          if (not is_boundary_node[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];
+
+            for (size_t i_wave = 0; i_wave < nb_waves; ++i_wave) {
+              double sum = 0;
+
+              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];
+
+                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 = nlr(face_id, i_node_face);
+                  const double sign            = face_cell_is_reversed(cell_id, i_face_cell) ? -1 : 1;
+
+                  double li_nlr = sign * dot(m_lambda_vector[i_wave], n_face);
+
+                  if (li_nlr > 0) {
+                    Fr[node_id][i_wave] += li_nlr * DPk_Fn(cell_id, i_wave)(m_xl[face_id]);
+                    sum += li_nlr;
+                  }
+                }
+              }
+              if (sum != 0) {
+                Fr[node_id][i_wave] = 1. / sum * Fr[node_id][i_wave];
+              }
+            }
+          }
+        });
+
+      return Fr;
+    } else {
+      throw NormalError("Eucclhyd not defined in 1d");
+    }
+  }
+
   void
-  apply_bc(FaceArray<double> Fl_rho,
-           FaceArray<TinyVector<Dimension>> Fl_rho_u,
-           FaceArray<double> Fl_rho_E,
-           const CellValue<const double> rho,
-           const CellValue<const TinyVector<Dimension>> rho_u,
-           const CellValue<const double> rho_E,
-           const CellArray<const double> Fn_rho,
-           const CellArray<const TinyVector<Dimension>> Fn_rho_u,
-           const CellArray<const double> Fn_rho_E,
-           const BoundaryConditionList& bc_list)
+  apply_face_bc(FaceArray<double> Fl_rho,
+                FaceArray<TinyVector<Dimension>> Fl_rho_u,
+                FaceArray<double> Fl_rho_E,
+                DiscreteFunctionDPkVector<Dimension, double> Fn_rho,
+                DiscreteFunctionDPkVector<Dimension, TinyVector<Dimension>> Fn_rho_u,
+                DiscreteFunctionDPkVector<Dimension, double> Fn_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();
@@ -537,27 +646,36 @@ 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];
 
+                  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 += Fn_rho(cell_id, k_wave)(m_xl[face_id]);
+                    rho_ul += Fn_rho_u(cell_id, k_wave)(m_xl[face_id]);
+                    rho_El += Fn_rho_E(cell_id, k_wave)(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 < 0) {
-                    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];
+                    double rhol_prime                  = rhol;
+                    TinyVector<Dimension> rho_ul_prime = txt * rho_ul - nxn * rho_ul;
+                    double rho_El_prime                = rho_El;
 
-                    double rho_e = rho_E[cell_id] - 0.5 * (1. / rho[cell_id]) * dot(rho_u[cell_id], rho_u[cell_id]);
-                    double p     = (m_gamma - 1) * rho_e;
+                    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_uj_prime;
+                    TinyVector<Dimension> A_rho_prime = rho_ul_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;
+                      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                  = 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;
+                    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];
@@ -587,8 +705,18 @@ class EulerKineticSolverMultiDOrder2
                 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];
-                  const double sign        = face_cell_is_reversed(cell_id, i_face_cell) ? -1 : 1;
 
+                  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 += Fn_rho(cell_id, k_wave)(m_xl[face_id]);
+                    rho_ul += Fn_rho_u(cell_id, k_wave)(m_xl[face_id]);
+                    rho_El += Fn_rho_E(cell_id, k_wave)(m_xl[face_id]);
+                  }
+
+                  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;
@@ -596,22 +724,21 @@ class EulerKineticSolverMultiDOrder2
 
                   double li_nl = dot(m_lambda_vector[i_wave], n);
                   if (li_nl < 0) {
-                    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];
+                    double rhol_prime                  = rhol;
+                    TinyVector<Dimension> rho_ul_prime = txt * rho_ul - nxn * rho_ul;
+                    double rho_El_prime                = rho_El;
 
-                    double rho_e = rho_E[cell_id] - 0.5 * (1. / rho[cell_id]) * dot(rho_u[cell_id], rho_u[cell_id]);
-                    double p     = (m_gamma - 1) * rho_e;
+                    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_uj_prime;
+                    TinyVector<Dimension> A_rho_prime = rho_ul_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;
+                      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                  = 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;
+                    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];
@@ -647,9 +774,9 @@ class EulerKineticSolverMultiDOrder2
 
                   double li_nl = dot(m_lambda_vector[i_wave], n);
                   if (li_nl < 0) {
-                    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];
+                    Fl_rho[face_id][i_wave] += Fn_rho(cell_id, i_wave)(m_xl[face_id]);
+                    Fl_rho_u[face_id][i_wave] += Fn_rho_u(cell_id, i_wave)(m_xl[face_id]);
+                    Fl_rho_E[face_id][i_wave] += Fn_rho_E(cell_id, i_wave)(m_xl[face_id]);
                   }
                 }
               }
@@ -660,88 +787,654 @@ class EulerKineticSolverMultiDOrder2
     }
   }
 
-  DiscreteFunctionP0Vector<double>
-  compute_deltaLFn(NodeArray<double> Fr)
+  void
+  apply_glace_bc(NodeArray<double> Fr_rho,
+                 NodeArray<TinyVector<Dimension>> Fr_rho_u,
+                 NodeArray<double> Fr_rho_E,
+                 DiscreteFunctionDPkVector<Dimension, double> Fn_rho,
+                 DiscreteFunctionDPkVector<Dimension, TinyVector<Dimension>> Fn_rho_u,
+                 DiscreteFunctionDPkVector<Dimension, double> Fn_rho_E,
+                 const BoundaryConditionList& bc_list)
   {
-    if constexpr (Dimension > 1) {
-      const NodeValuePerCell<const TinyVector<Dimension>> Cjr = MeshDataManager::instance().getMeshData(m_mesh).Cjr();
-      const size_t nb_waves                                   = m_lambda_vector.size();
-      DiscreteFunctionP0Vector<double> deltaLFn(p_mesh, nb_waves);
+    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();
 
-      const auto& cell_to_node_matrix = p_mesh->connectivity().cellToNodeMatrix();
+    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->numberOfCells(), PUGS_LAMBDA(CellId cell_id) {
-          const auto& cell_to_node = cell_to_node_matrix[cell_id];
+    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();
 
-          for (size_t i_wave = 0; i_wave < nb_waves; ++i_wave) {
-            deltaLFn[cell_id][i_wave] = 0;
-            for (size_t i_node = 0; i_node < cell_to_node.size(); ++i_node) {
-              const NodeId node_id = cell_to_node[i_node];
+            TinyMatrix<Dimension> I = identity;
 
-              double li_Cjr = dot(m_lambda_vector[i_wave], Cjr(cell_id, i_node));
-              deltaLFn[cell_id][i_wave] += Fr[node_id][i_wave] * li_Cjr;
+            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;
+
+                  for (size_t k_wave = 0; k_wave < nb_waves; ++k_wave) {
+                    rhor += Fn_rho(cell_id, k_wave)(m_xr[node_id]);
+                    rho_ur += Fn_rho_u(cell_id, k_wave)(m_xr[node_id]);
+                    rho_Er += Fn_rho_E(cell_id, k_wave)(m_xr[node_id]);
+                  }
+
+                  double li_njr = dot(m_lambda_vector[i_wave], njr(cell_id, i_node_cell));
+
+                  if (li_njr > 0) {
+                    Fr_rho[node_id][i_wave] += Fn_rho(cell_id, i_wave)(m_xr[node_id]) * li_njr;
+                    Fr_rho_u[node_id][i_wave] += li_njr * Fn_rho_u(cell_id, i_wave)(m_xr[node_id]);
+                    Fr_rho_E[node_id][i_wave] += Fn_rho_E(cell_id, i_wave)(m_xr[node_id]) * 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) {
+                    double rhor_prime                  = rhor;
+                    TinyVector<Dimension> rho_ur_prime = txt * rho_ur - nxn * rho_ur;
+                    double rho_Er_prime                = rho_Er;
+
+                    double rho_er = rho_Er - 0.5 * (1. / rhor) * dot(rho_ur, rho_ur);
+                    double pr     = (m_gamma - 1) * rho_er;
+
+                    TinyVector<Dimension> A_rho_prime = rho_ur_prime;
+                    TinyMatrix<Dimension> A_rho_u_prime =
+                      1. / rhor_prime * tensorProduct(rho_ur_prime, rho_ur_prime) + pr * I;
+                    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];
+
+                  double rhor                  = 0;
+                  TinyVector<Dimension> rho_ur = zero;
+                  double rho_Er                = 0;
+
+                  for (size_t k_wave = 0; k_wave < nb_waves; ++k_wave) {
+                    rhor += Fn_rho(cell_id, k_wave)(m_xr[node_id]);
+                    rho_ur += Fn_rho_u(cell_id, k_wave)(m_xr[node_id]);
+                    rho_Er += Fn_rho_E(cell_id, k_wave)(m_xr[node_id]);
+                  }
+
+                  double li_njr = dot(m_lambda_vector[i_wave], njr(cell_id, i_node_cell));
+
+                  if (li_njr > 0) {
+                    Fr_rho[node_id][i_wave] += Fn_rho(cell_id, i_wave)(m_xr[node_id]) * li_njr;
+                    Fr_rho_u[node_id][i_wave] += li_njr * Fn_rho_u(cell_id, i_wave)(m_xr[node_id]);
+                    Fr_rho_E[node_id][i_wave] += Fn_rho_E(cell_id, i_wave)(m_xr[node_id]) * 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) {
+                    double rhor_prime                  = rhor;
+                    TinyVector<Dimension> rho_ur_prime = txt * rho_ur - nxn * rho_ur;
+                    double rho_Er_prime                = rho_Er;
+
+                    double rho_er = rho_Er - 0.5 * (1. / rhor) * dot(rho_ur, rho_ur);
+                    double pr     = (m_gamma - 1) * rho_er;
+
+                    TinyVector<Dimension> A_rho_prime = rho_ur_prime;
+                    TinyMatrix<Dimension> A_rho_u_prime =
+                      1. / rhor_prime * tensorProduct(rho_ur_prime, rho_ur_prime) + pr * I;
+                    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) {
+                    Fr_rho[node_id][i_wave] += Fn_rho(cell_id, i_wave)(m_xr[node_id]) * li_njr;
+                    Fr_rho_u[node_id][i_wave] += li_njr * Fn_rho_u(cell_id, i_wave)(m_xr[node_id]);
+                    Fr_rho_E[node_id][i_wave] += Fn_rho_E(cell_id, i_wave)(m_xr[node_id]) * 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) {
+                    Fr_rho[node_id][i_wave] += Fn_rho(cell_id, i_wave)(m_xr[node_id]) * li_njr_prime;
+                    Fr_rho_u[node_id][i_wave] += li_njr_prime * Fn_rho_u(cell_id, i_wave)(m_xr[node_id]);
+                    Fr_rho_E[node_id][i_wave] += Fn_rho_E(cell_id, i_wave)(m_xr[node_id]) * 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;
+                }
+              }
             }
           }
-        });
-      return deltaLFn;
-    } else {
-      throw NormalError("Glace not defined in 1d");
+        },
+        bc_v);
     }
   }
 
-  DiscreteFunctionP0Vector<double>
-  compute_deltaLFn_eucclhyd(FaceArrayPerNode<double> Fr)
+  void
+  apply_eucclhyd_bc(NodeArray<double> Fr_rho,
+                    NodeArray<TinyVector<Dimension>> Fr_rho_u,
+                    NodeArray<double> Fr_rho_E,
+                    DiscreteFunctionDPkVector<Dimension, double> Fn_rho,
+                    DiscreteFunctionDPkVector<Dimension, TinyVector<Dimension>> Fn_rho_u,
+                    DiscreteFunctionDPkVector<Dimension, double> Fn_rho_E,
+                    const BoundaryConditionList& bc_list)
   {
-    if constexpr (Dimension > 1) {
-      const size_t nb_waves = m_lambda_vector.size();
-      DiscreteFunctionP0Vector<double> deltaLFn(p_mesh, nb_waves);
+    const size_t nb_waves = m_lambda_vector.size();
 
-      const NodeValuePerFace<const TinyVector<Dimension>> Nlr = MeshDataManager::instance().getMeshData(m_mesh).Nlr();
+    const NodeValuePerFace<const TinyVector<Dimension>> nlr = MeshDataManager::instance().getMeshData(m_mesh).nlr();
+    const auto& face_cell_is_reversed                       = p_mesh->connectivity().cellFaceIsReversed();
+    const auto& node_local_numbers_in_their_cells           = p_mesh->connectivity().nodeLocalNumbersInTheirCells();
+    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& face_to_cell_matrix                         = p_mesh->connectivity().faceToCellMatrix();
+    const auto& node_to_cell_matrix                         = p_mesh->connectivity().nodeToCellMatrix();
+    const auto& node_to_face_matrix                         = p_mesh->connectivity().nodeToFaceMatrix();
+    const NodeValuePerCell<const TinyVector<Dimension>> njr = MeshDataManager::instance().getMeshData(m_mesh).njr();
+    const NodeValuePerCell<const TinyVector<Dimension>> Cjr = MeshDataManager::instance().getMeshData(m_mesh).Cjr();
 
-      const auto& cell_to_face_matrix   = p_mesh->connectivity().cellToFaceMatrix();
-      const auto& face_to_node_matrix   = p_mesh->connectivity().faceToNodeMatrix();
-      const auto& node_to_face_matrix   = p_mesh->connectivity().nodeToFaceMatrix();
-      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];
+    }
+    FaceArrayPerNode<double> sum_li_nlr(m_mesh.connectivity(), nb_waves);
+    sum_li_nlr.fill(0);
 
-      parallel_for(
-        p_mesh->numberOfCells(), PUGS_LAMBDA(CellId cell_id) {
-          const auto& cell_to_face = cell_to_face_matrix[cell_id];
+    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, WallBoundaryCondition>) {
+            auto node_list = bc.nodeList();
 
-          for (size_t i_wave = 0; i_wave < nb_waves; ++i_wave) {
-            deltaLFn[cell_id][i_wave] = 0;
+            TinyMatrix<Dimension> I = identity;
 
-            for (size_t i_face_cell = 0; i_face_cell < cell_to_face.size(); ++i_face_cell) {
-              const FaceId face_id     = cell_to_face[i_face_cell];
-              const auto& face_to_node = face_to_node_matrix[face_id];
+            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) {
+                const auto& node_to_face                  = node_to_face_matrix[node_id];
+                const auto& node_local_number_in_its_face = node_local_numbers_in_their_faces.itemArray(node_id);
+
+                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_face = 0; i_face < node_to_face.size(); ++i_face) {
+                  const FaceId face_id                      = node_to_face[i_face];
+                  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 auto& face_local_number_in_its_cell = face_local_numbers_in_their_cells.itemArray(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 += Fn_rho(cell_id, k_wave)(m_xl[face_id]);
+                      rho_ul += Fn_rho_u(cell_id, k_wave)(m_xl[face_id]);
+                      rho_El += Fn_rho_E(cell_id, k_wave)(m_xl[face_id]);
+                    }
+
+                    const double sign            = face_cell_is_reversed(cell_id, i_face_cell) ? -1 : 1;
+                    TinyVector<Dimension> n_face = sign * nlr(face_id, i_node_face);
+
+                    double li_nlr = dot(m_lambda_vector[i_wave], n_face);
+
+                    if (li_nlr > 0) {
+                      Fr_rho[node_id][i_wave] += Fn_rho(cell_id, i_wave)(m_xl[face_id]) * li_nlr;
+                      Fr_rho_u[node_id][i_wave] += li_nlr * Fn_rho_u(cell_id, i_wave)(m_xl[face_id]);
+                      Fr_rho_E[node_id][i_wave] += Fn_rho_E(cell_id, i_wave)(m_xl[face_id]) * li_nlr;
+
+                      sum += li_nlr;
+                    }
 
-              for (size_t i_node_face = 0; i_node_face < face_to_node.size(); ++i_node_face) {
-                const NodeId node_id     = face_to_node[i_node_face];
-                const auto& node_to_face = node_to_face_matrix[node_id];
+                    TinyVector<Dimension> nlr_prime = txt * n_face - nxn * n_face;
+                    double li_nlr_prime             = dot(m_lambda_vector[i_wave], nlr_prime);
 
-                auto local_face_number_in_node = [&](FaceId face_number) {
-                  for (size_t i_face = 0; i_face < node_to_face.size(); ++i_face) {
-                    if (face_number == node_to_face[i_face]) {
-                      return i_face;
+                    if (li_nlr_prime > 0) {
+                      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];
+                      }
+
+                      Fr_rho[node_id][i_wave] += Fn_rho_prime * li_nlr_prime;
+                      Fr_rho_u[node_id][i_wave] += li_nlr_prime * Fn_rho_u_prime;
+                      Fr_rho_E[node_id][i_wave] += Fn_rho_E_prime * li_nlr_prime;
+
+                      sum += li_nlr_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, 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) {
+                const auto& node_to_face                  = node_to_face_matrix[node_id];
+                const auto& node_local_number_in_its_face = node_local_numbers_in_their_faces.itemArray(node_id);
+
+                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_face = 0; i_face < node_to_face.size(); ++i_face) {
+                  const FaceId face_id                      = node_to_face[i_face];
+                  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 auto& face_local_number_in_its_cell = face_local_numbers_in_their_cells.itemArray(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 += Fn_rho(cell_id, k_wave)(m_xl[face_id]);
+                      rho_ul += Fn_rho_u(cell_id, k_wave)(m_xl[face_id]);
+                      rho_El += Fn_rho_E(cell_id, k_wave)(m_xl[face_id]);
+                    }
+
+                    const double sign            = face_cell_is_reversed(cell_id, i_face_cell) ? -1 : 1;
+                    TinyVector<Dimension> n_face = sign * nlr(face_id, i_node_face);
+
+                    double li_nlr = dot(m_lambda_vector[i_wave], n_face);
+
+                    if (li_nlr > 0) {
+                      Fr_rho[node_id][i_wave] += Fn_rho(cell_id, i_wave)(m_xl[face_id]) * li_nlr;
+                      Fr_rho_u[node_id][i_wave] += li_nlr * Fn_rho_u(cell_id, i_wave)(m_xl[face_id]);
+                      Fr_rho_E[node_id][i_wave] += Fn_rho_E(cell_id, i_wave)(m_xl[face_id]) * li_nlr;
+
+                      sum += li_nlr;
+                    }
+
+                    TinyVector<Dimension> nlr_prime = txt * n_face - nxn * n_face;
+                    double li_nlr_prime             = dot(m_lambda_vector[i_wave], nlr_prime);
+
+                    if (li_nlr_prime > 0) {
+                      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];
+                      }
+
+                      Fr_rho[node_id][i_wave] += Fn_rho_prime * li_nlr_prime;
+                      Fr_rho_u[node_id][i_wave] += li_nlr_prime * Fn_rho_u_prime;
+                      Fr_rho_E[node_id][i_wave] += Fn_rho_E_prime * li_nlr_prime;
+
+                      sum += li_nlr_prime;
                     }
                   }
-                  return std::numeric_limits<size_t>::max();
-                };
+                }
+                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) {
+                const auto& node_to_face                  = node_to_face_matrix[node_id];
+                const auto& node_local_number_in_its_face = node_local_numbers_in_their_faces.itemArray(node_id);
+
+                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_face = 0; i_face < node_to_face.size(); ++i_face) {
+                  const FaceId face_id                      = node_to_face[i_face];
+                  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 auto& face_local_number_in_its_cell = face_local_numbers_in_their_cells.itemArray(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];
 
-                const size_t i_face_node = local_face_number_in_node(face_id);
+                    const double sign            = face_cell_is_reversed(cell_id, i_face_cell) ? -1 : 1;
+                    TinyVector<Dimension> n_face = sign * nlr(face_id, i_node_face);
 
-                const double sign = face_cell_is_reversed(cell_id, i_face_cell) ? -1 : 1;
+                    double li_nlr = dot(m_lambda_vector[i_wave], n_face);
 
-                const double li_Nlr = sign * dot(m_lambda_vector[i_wave], Nlr(face_id, i_node_face));
+                    if (li_nlr > 0) {
+                      Fr_rho[node_id][i_wave] += Fn_rho(cell_id, i_wave)(m_xl[face_id]) * li_nlr;
+                      Fr_rho_u[node_id][i_wave] += li_nlr * Fn_rho_u(cell_id, i_wave)(m_xl[face_id]);
+                      Fr_rho_E[node_id][i_wave] += Fn_rho_E(cell_id, i_wave)(m_xl[face_id]) * li_nlr;
 
-                deltaLFn[cell_id][i_wave] += Fr[node_id][i_face_node][i_wave] * li_Nlr;
+                      sum += li_nlr;
+                    }
+
+                    TinyVector<Dimension> nlr_prime = txt * n_face - nxn * n_face;
+                    double li_nlr_prime             = dot(m_lambda_vector[i_wave], nlr_prime);
+
+                    if (li_nlr_prime > 0) {
+                      Fr_rho[node_id][i_wave] += Fn_rho(cell_id, i_wave)(m_xl[face_id]) * li_nlr_prime;
+                      Fr_rho_u[node_id][i_wave] += li_nlr_prime * Fn_rho_u(cell_id, i_wave)(m_xl[face_id]);
+                      Fr_rho_E[node_id][i_wave] += Fn_rho_E(cell_id, i_wave)(m_xl[face_id]) * li_nlr_prime;
+
+                      sum += li_nlr_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);
+    }
+  }
+
+  template <typename T>
+  CellArray<T>
+  compute_deltaLFn_node(NodeArray<T> Fr)
+  {
+    if constexpr (Dimension > 1) {
+      const NodeValuePerCell<const TinyVector<Dimension>> Cjr = MeshDataManager::instance().getMeshData(m_mesh).Cjr();
+      const size_t nb_waves                                   = m_lambda_vector.size();
+      CellArray<T> deltaLFn(p_mesh->connectivity(), nb_waves);
+
+      const auto& cell_to_node_matrix = p_mesh->connectivity().cellToNodeMatrix();
+
+      if constexpr (is_tiny_vector_v<T>) {
+        deltaLFn.fill(zero);
+      } else {
+        deltaLFn.fill(0.);
+      }
+
+      parallel_for(
+        p_mesh->numberOfCells(), PUGS_LAMBDA(CellId cell_id) {
+          const auto& cell_to_node = cell_to_node_matrix[cell_id];
+
+          for (size_t i_wave = 0; i_wave < nb_waves; ++i_wave) {
+            for (size_t i_node = 0; i_node < cell_to_node.size(); ++i_node) {
+              const NodeId node_id = cell_to_node[i_node];
+
+              double li_Cjr = dot(m_lambda_vector[i_wave], Cjr(cell_id, i_node));
+              deltaLFn[cell_id][i_wave] += li_Cjr * Fr[node_id][i_wave];
+            }
+          }
         });
       return deltaLFn;
     } else {
-      throw NormalError("Eucclhyd not defined in 1d");
+      throw NormalError("Nodal solver not defined in 1d");
     }
   }
 
@@ -833,30 +1526,23 @@ class EulerKineticSolverMultiDOrder2
             f_max = std::max(f_max, f[cell_stencil[i_cell]][i_wave]);
           }
 
-          double f_bar_min = fj;
-          double f_bar_max = fj;
-
-          for (size_t i_cell = 0; i_cell < cell_stencil.size(); ++i_cell) {
-            const CellId cell_k_id = cell_stencil[i_cell];
-            const double f_xk      = DPk_fh(cell_id, i_wave)(m_xj[cell_k_id]);
-
-            f_bar_min = std::min(f_bar_min, f_xk);
-            f_bar_max = std::max(f_bar_max, f_xk);
-          }
-
           const double eps = 1E-14;
-          double coef1     = 1;
-          if (std::abs(f_bar_max - fj) > eps) {
-            coef1 = (f_max - fj) / ((f_bar_max - fj));
-          }
+          double coef      = 1;
+          double lambda    = 1.;
 
-          double coef2 = 1.;
-          if (std::abs(f_bar_min - fj) > eps) {
-            coef2 = (f_min - fj) / ((f_bar_min - fj));
+          for (size_t i_cell = 0; i_cell < cell_stencil.size(); ++i_cell) {
+            const CellId cell_i_id = cell_stencil[i_cell];
+            const double f_bar_i   = DPk_fh(cell_id, i_wave)(m_xj[cell_i_id]);
+            const double Delta     = (f_bar_i - fj);
+
+            if (Delta > eps) {
+              coef = std::min(1., (f_max - fj) / Delta);
+            } else if (Delta < -eps) {
+              coef = std::min(1., (f_min - fj) / Delta);
+            }
+            lambda = std::min(lambda, coef);
           }
 
-          const double lambda = std::max(0., std::min(1., std::min(coef1, coef2)));
-
           auto coefficients = DPk_fh.componentCoefficients(cell_id, i_wave);
           coefficients[0]   = (1 - lambda) * f[cell_id][i_wave] + lambda * coefficients[0];
           for (size_t i = 1; i < coefficients.size(); ++i) {
@@ -892,45 +1578,26 @@ class EulerKineticSolverMultiDOrder2
             }
           }
 
-          TinyVector<Dimension> f_bar_min = fj;
-          TinyVector<Dimension> f_bar_max = fj;
+          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_k_id           = cell_stencil[i_cell];
-            const TinyVector<Dimension> f_xk = DPk_fh(cell_id, i_wave)(m_xj[cell_k_id]);
+            const CellId cell_i_id              = cell_stencil[i_cell];
+            const TinyVector<Dimension> f_bar_i = DPk_fh(cell_id, i_wave)(m_xj[cell_i_id]);
+            TinyVector<Dimension> Delta;
 
             for (size_t dim = 0; dim < Dimension; ++dim) {
-              f_bar_min[dim] = std::min(f_bar_min[dim], f_xk[dim]);
-              f_bar_max[dim] = std::max(f_bar_max[dim], f_xk[dim]);
-            }
-          }
-
-          const double eps = 1E-14;
-          TinyVector<Dimension> coef1;
-          for (size_t dim = 0; dim < Dimension; ++dim) {
-            coef1[dim] = 1;
-            if (std::abs(f_bar_max[dim] - fj[dim]) > eps) {
-              coef1[dim] = (f_max[dim] - fj[dim]) / ((f_bar_max[dim] - fj[dim]));
-            }
-          }
-
-          TinyVector<Dimension> coef2;
-          for (size_t dim = 0; dim < Dimension; ++dim) {
-            coef2[dim] = 1;
-            if (std::abs(f_bar_min[dim] - fj[dim]) > eps) {
-              coef2[dim] = (f_min[dim] - fj[dim]) / ((f_bar_min[dim] - fj[dim]));
+              Delta[dim] = (f_bar_i[dim] - fj[dim]);
+              coef[dim]  = 1.;
+              if (Delta[dim] > eps) {
+                coef[dim] = std::min(1., (f_max[dim] - fj[dim]) / Delta[dim]);
+              } else if (Delta[dim] < -eps) {
+                coef[dim] = std::min(1., (f_min[dim] - fj[dim]) / Delta[dim]);
+              }
+              lambda = std::min(lambda, coef[dim]);
             }
           }
-
-          double min_coef1 = coef1[0];
-          double min_coef2 = coef2[0];
-          for (size_t dim = 1; dim < Dimension; ++dim) {
-            min_coef1 = std::min(min_coef1, coef1[dim]);
-            min_coef2 = std::min(min_coef2, coef2[dim]);
-          }
-
-          const double lambda = std::max(0., std::min(1., std::min(min_coef1, min_coef2)));
-
           auto coefficients = DPk_fh.componentCoefficients(cell_id, i_wave);
 
           coefficients[0] = (1 - lambda) * f[cell_id][i_wave] + lambda * coefficients[0];
@@ -979,6 +1646,27 @@ class EulerKineticSolverMultiDOrder2
         bc_v);
     }
 
+    NodeValue<bool> is_boundary_node{p_mesh->connectivity()};
+    is_boundary_node.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, SymmetryBoundaryCondition>) {
+            auto node_list = bc.nodeList();
+            for (size_t i_node = 0; i_node < node_list.size(); ++i_node) {
+              is_boundary_node[node_list[i_node]] = 1;
+            }
+          } else if constexpr (std::is_same_v<BCType, WallBoundaryCondition>) {
+            auto node_list = bc.nodeList();
+            for (size_t i_node = 0; i_node < node_list.size(); ++i_node) {
+              is_boundary_node[node_list[i_node]] = 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);
@@ -1035,20 +1723,57 @@ class EulerKineticSolverMultiDOrder2
 
     vector_limiter(Fn_rho_u, Fn_rho_u_lim);
 
-    FaceArray<double> Fl_rho                  = compute_Flux1_face<double>(Fn_rho_lim);
-    FaceArray<TinyVector<Dimension>> Fl_rho_u = compute_Flux1_face<TinyVector<Dimension>>(Fn_rho_u_lim);
-    FaceArray<double> Fl_rho_E                = compute_Flux1_face<double>(Fn_rho_E_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);
+      FaceArray<TinyVector<Dimension>> Fl_rho_u = compute_Flux_face<TinyVector<Dimension>>(Fn_rho_u_lim);
+      FaceArray<double> Fl_rho_E                = compute_Flux_face<double>(Fn_rho_E_lim);
+
+      apply_face_bc(Fl_rho, Fl_rho_u, Fl_rho_E, Fn_rho_lim, Fn_rho_u_lim, Fn_rho_E_lim, bc_list);
 
-    apply_bc(Fl_rho, Fl_rho_u, Fl_rho_E, rho, rho_u, rho_E, Fn_rho, Fn_rho_u, Fn_rho_E, bc_list);
+      synchronize(Fl_rho);
+      synchronize(Fl_rho_u);
+      synchronize(Fl_rho_E);
 
-    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);
 
-    const CellArray<const double> deltaLFn_rho = compute_deltaLFn_face<double>(Fl_rho);
-    const CellArray<const TinyVector<Dimension>> deltaLFn_rho_u =
-      compute_deltaLFn_face<TinyVector<Dimension>>(Fl_rho_u);
-    const CellArray<const double> deltaLFn_rho_E = compute_deltaLFn_face<double>(Fl_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> Fr_rho_E = compute_Flux_glace<double>(Fn_rho_E_lim, is_boundary_node);
+
+      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);
+
+      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 == 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);
+
+      apply_eucclhyd_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);
+
+      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);
     const CellValue<const TinyMatrix<Dimension>> A_rho_u = getA_rho_u(rho, rho_u, rho_E);
@@ -1078,21 +1803,56 @@ class EulerKineticSolverMultiDOrder2
 
       vector_limiter(Fnp1_rho_u, Fnp1_rho_u_lim);
 
-      FaceArray<double> Fl_rho_np1                  = compute_Flux1_face<double>(Fnp1_rho_lim);
-      FaceArray<TinyVector<Dimension>> Fl_rho_u_np1 = compute_Flux1_face<TinyVector<Dimension>>(Fnp1_rho_u_lim);
-      FaceArray<double> Fl_rho_E_np1                = compute_Flux1_face<double>(Fnp1_rho_E_lim);
-
-      apply_bc(Fl_rho_np1, Fl_rho_u_np1, Fl_rho_E_np1, rho_np1, rho_u_np1, rho_E_np1, Fnp1_rho, Fnp1_rho_u, Fnp1_rho_E,
-               bc_list);
-
-      synchronize(Fl_rho_np1);
-      synchronize(Fl_rho_u_np1);
-      synchronize(Fl_rho_E_np1);
-
-      const CellArray<const double> deltaLFnp1_rho = compute_deltaLFn_face<double>(Fl_rho_np1);
-      const CellArray<const TinyVector<Dimension>> deltaLFnp1_rho_u =
-        compute_deltaLFn_face<TinyVector<Dimension>>(Fl_rho_u_np1);
-      const CellArray<const double> deltaLFnp1_rho_E = compute_deltaLFn_face<double>(Fl_rho_E_np1);
+      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);
+        FaceArray<TinyVector<Dimension>> Fl_rho_u_np1 = compute_Flux_face<TinyVector<Dimension>>(Fnp1_rho_u_lim);
+        FaceArray<double> Fl_rho_E_np1                = compute_Flux_face<double>(Fnp1_rho_E_lim);
+
+        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, 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);
+      } 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);
+      }
 
       rho_np1   = compute_PFnp1<double>(Fn_rho, deltaLFn_rho, deltaLFnp1_rho);
       rho_u_np1 = compute_PFnp1<TinyVector<Dimension>>(Fn_rho_u, deltaLFn_rho_u, deltaLFnp1_rho_u);
@@ -1182,6 +1942,7 @@ class EulerKineticSolverMultiDOrder2
                                  const std::vector<TinyVector<MeshType::Dimension>>& lambda_vector,
                                  const size_t& SpaceOrder,
                                  const size_t& TimeOrder,
+                                 const size_t& scheme,
                                  const CellArray<const double>& Fn_rho,
                                  const CellArray<const TinyVector<MeshType::Dimension>>& Fn_rho_u,
                                  const CellArray<const double>& Fn_rho_E)
@@ -1190,6 +1951,7 @@ class EulerKineticSolverMultiDOrder2
       m_lambda_vector{lambda_vector},
       m_SpaceOrder{SpaceOrder},
       m_TimeOrder{TimeOrder},
+      m_scheme{scheme},
       m_Fn_rho{Fn_rho},
       m_Fn_rho_u{Fn_rho_u},
       m_Fn_rho_E{Fn_rho_E},
@@ -1410,6 +2172,7 @@ eulerKineticSolverMultiDOrder2(
   const double& eps,
   const size_t& SpaceOrder,
   const size_t& TimeOrder,
+  const size_t& scheme,
   const std::shared_ptr<const DiscreteFunctionVariant>& F_rho,
   const std::shared_ptr<const DiscreteFunctionVariant>& F_rho_u,
   const std::shared_ptr<const DiscreteFunctionVariant>& F_rho_E,
@@ -1433,7 +2196,7 @@ eulerKineticSolverMultiDOrder2(
         auto Fn_rho_u = F_rho_u->get<DiscreteFunctionP0Vector<const TinyVector<MeshType::Dimension>>>();
         auto Fn_rho_E = F_rho_E->get<DiscreteFunctionP0Vector<const double>>();
 
-        EulerKineticSolverMultiDOrder2 solver(p_mesh, dt, eps, gamma, lambda_vector, SpaceOrder, TimeOrder,
+        EulerKineticSolverMultiDOrder2 solver(p_mesh, dt, eps, gamma, lambda_vector, SpaceOrder, TimeOrder, scheme,
                                               Fn_rho.cellArrays(), Fn_rho_u.cellArrays(), Fn_rho_E.cellArrays());
 
         return solver.apply(bc_descriptor_list);
@@ -1454,6 +2217,7 @@ eulerKineticSolverMultiDOrder2(
   const double& eps,
   const size_t& SpaceOrder,
   const size_t& TimeOrder,
+  const size_t& scheme,
   const std::shared_ptr<const DiscreteFunctionVariant>& F_rho,
   const std::shared_ptr<const DiscreteFunctionVariant>& F_rho_u,
   const std::shared_ptr<const DiscreteFunctionVariant>& F_rho_E,
@@ -1469,6 +2233,7 @@ eulerKineticSolverMultiDOrder2(
   const double& eps,
   const size_t& SpaceOrder,
   const size_t& TimeOrder,
+  const size_t& scheme,
   const std::shared_ptr<const DiscreteFunctionVariant>& F_rho,
   const std::shared_ptr<const DiscreteFunctionVariant>& F_rho_u,
   const std::shared_ptr<const DiscreteFunctionVariant>& F_rho_E,
diff --git a/src/scheme/EulerKineticSolverMultiDOrder2.hpp b/src/scheme/EulerKineticSolverMultiDOrder2.hpp
index c533e9d12..5eebcaf5a 100644
--- a/src/scheme/EulerKineticSolverMultiDOrder2.hpp
+++ b/src/scheme/EulerKineticSolverMultiDOrder2.hpp
@@ -17,6 +17,7 @@ eulerKineticSolverMultiDOrder2(
   const double& eps,
   const size_t& SpaceOrder,
   const size_t& TimeOrder,
+  const size_t& scheme,
   const std::shared_ptr<const DiscreteFunctionVariant>& F_rho,
   const std::shared_ptr<const DiscreteFunctionVariant>& F_rho_u,
   const std::shared_ptr<const DiscreteFunctionVariant>& F_rho_E,
-- 
GitLab