diff --git a/src/language/modules/SchemeModule.cpp b/src/language/modules/SchemeModule.cpp
index ea6a01fbf0752ba0a26d63711c796057d072d1fb..d4326d14b5b5fc94d459131a2942f7a5df57ebc4 100644
--- a/src/language/modules/SchemeModule.cpp
+++ b/src/language/modules/SchemeModule.cpp
@@ -495,6 +495,39 @@ SchemeModule::SchemeModule()
 
                               ));
 
+    this->_addBuiltinFunction("local_dt_eucclhyd_solver_order2",
+                            std::function(
+
+                              [](const std::shared_ptr<const DiscreteFunctionVariant>& rho1,
+                                 const std::shared_ptr<const DiscreteFunctionVariant>& u1,
+                                 const std::shared_ptr<const DiscreteFunctionVariant>& E1,
+                                 const std::shared_ptr<const DiscreteFunctionVariant>& c1,
+                                 const std::shared_ptr<const DiscreteFunctionVariant>& p1,
+                                 const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>&
+                                   bc_descriptor_list1,
+                                 const double& dt1, const std::shared_ptr<const DiscreteFunctionVariant>& rho2,
+                                 const std::shared_ptr<const DiscreteFunctionVariant>& u2,
+                                 const std::shared_ptr<const DiscreteFunctionVariant>& E2,
+                                 const std::shared_ptr<const DiscreteFunctionVariant>& c2,
+                                 const std::shared_ptr<const DiscreteFunctionVariant>& p2,
+                                 const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>&
+                                   bc_descriptor_list2)
+                                -> std::tuple<
+                                  std::shared_ptr<const IMesh>, std::shared_ptr<const DiscreteFunctionVariant>,
+                                  std::shared_ptr<const DiscreteFunctionVariant>,
+                                  std::shared_ptr<const DiscreteFunctionVariant>, std::shared_ptr<const IMesh>,
+                                  std::shared_ptr<const DiscreteFunctionVariant>,
+                                  std::shared_ptr<const DiscreteFunctionVariant>,
+                                  std::shared_ptr<const DiscreteFunctionVariant>> {
+                                return LocalDtAcousticSolverHandler{getCommonMesh({rho1, u1, E1, c1, p1}),
+                                                                    getCommonMesh({rho2, u2, E2, c2, p2})}
+                                  .solver()
+                                  .applyOrder2(LocalDtAcousticSolverHandler::SolverType::Eucclhyd, dt1, rho1, rho2, u1, u2,
+                                         E1, E2, c1, c2, p1, p2, bc_descriptor_list1, bc_descriptor_list2);
+                              }
+
+                              ));
+
   this->_addBuiltinFunction("local_dt_glace_solver",
                             std::function(
 
@@ -529,6 +562,39 @@ SchemeModule::SchemeModule()
 
                               ));
 
+  this->_addBuiltinFunction("local_dt_glace_solver",
+                            std::function(
+
+                              [](const std::shared_ptr<const DiscreteFunctionVariant>& rho1,
+                                 const std::shared_ptr<const DiscreteFunctionVariant>& u1,
+                                 const std::shared_ptr<const DiscreteFunctionVariant>& E1,
+                                 const std::shared_ptr<const DiscreteFunctionVariant>& c1,
+                                 const std::shared_ptr<const DiscreteFunctionVariant>& p1,
+                                 const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>&
+                                   bc_descriptor_list1,
+                                 const double& dt1, const std::shared_ptr<const DiscreteFunctionVariant>& rho2,
+                                 const std::shared_ptr<const DiscreteFunctionVariant>& u2,
+                                 const std::shared_ptr<const DiscreteFunctionVariant>& E2,
+                                 const std::shared_ptr<const DiscreteFunctionVariant>& c2,
+                                 const std::shared_ptr<const DiscreteFunctionVariant>& p2,
+                                 const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>&
+                                   bc_descriptor_list2)
+                                -> std::tuple<
+                                  std::shared_ptr<const IMesh>, std::shared_ptr<const DiscreteFunctionVariant>,
+                                  std::shared_ptr<const DiscreteFunctionVariant>,
+                                  std::shared_ptr<const DiscreteFunctionVariant>, std::shared_ptr<const IMesh>,
+                                  std::shared_ptr<const DiscreteFunctionVariant>,
+                                  std::shared_ptr<const DiscreteFunctionVariant>,
+                                  std::shared_ptr<const DiscreteFunctionVariant>> {
+                                return LocalDtAcousticSolverHandler{getCommonMesh({rho1, u1, E1, c1, p1}),
+                                                                    getCommonMesh({rho2, u2, E2, c2, p2})}
+                                  .solver()
+                                  .apply(LocalDtAcousticSolverHandler::SolverType::Glace, dt1,rho1, rho2, u1, u2,
+                                         E1, E2, c1, c2, p1, p2, bc_descriptor_list1, bc_descriptor_list2);
+                              }
+
+                              ));
+
   this->_addBuiltinFunction("local_dt_eucclhyd_solver",
                             std::function(
 
@@ -563,6 +629,39 @@ SchemeModule::SchemeModule()
 
                               ));
 
+    this->_addBuiltinFunction("local_dt_eucclhyd_solver",
+                            std::function(
+
+                              [](const std::shared_ptr<const DiscreteFunctionVariant>& rho1,
+                                 const std::shared_ptr<const DiscreteFunctionVariant>& u1,
+                                 const std::shared_ptr<const DiscreteFunctionVariant>& E1,
+                                 const std::shared_ptr<const DiscreteFunctionVariant>& c1,
+                                 const std::shared_ptr<const DiscreteFunctionVariant>& p1,
+                                 const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>&
+                                   bc_descriptor_list1,
+                                 const double& dt1, const std::shared_ptr<const DiscreteFunctionVariant>& rho2,
+                                 const std::shared_ptr<const DiscreteFunctionVariant>& u2,
+                                 const std::shared_ptr<const DiscreteFunctionVariant>& E2,
+                                 const std::shared_ptr<const DiscreteFunctionVariant>& c2,
+                                 const std::shared_ptr<const DiscreteFunctionVariant>& p2,
+                                 const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>&
+                                   bc_descriptor_list2)
+                                -> std::tuple<
+                                  std::shared_ptr<const IMesh>, std::shared_ptr<const DiscreteFunctionVariant>,
+                                  std::shared_ptr<const DiscreteFunctionVariant>,
+                                  std::shared_ptr<const DiscreteFunctionVariant>, std::shared_ptr<const IMesh>,
+                                  std::shared_ptr<const DiscreteFunctionVariant>,
+                                  std::shared_ptr<const DiscreteFunctionVariant>,
+                                  std::shared_ptr<const DiscreteFunctionVariant>> {
+                                return LocalDtAcousticSolverHandler{getCommonMesh({rho1, u1, E1, c1, p1}),
+                                                                    getCommonMesh({rho2, u2, E2, c2, p2})}
+                                  .solver()
+                                  .apply(LocalDtAcousticSolverHandler::SolverType::Eucclhyd, dt1,rho1, rho2, u1, u2,
+                                         E1, E2, c1, c2, p1, p2, bc_descriptor_list1, bc_descriptor_list2);
+                              }
+
+                              ));
+
   this->_addBuiltinFunction("apply_acoustic_fluxes",
                             std::function(
 
diff --git a/src/scheme/LocalDtAcousticSolver.cpp b/src/scheme/LocalDtAcousticSolver.cpp
index 978593bdfef59017ada43da915cca6cc276c395d..62091415cda41c2aac061bb041320a236d35197b 100644
--- a/src/scheme/LocalDtAcousticSolver.cpp
+++ b/src/scheme/LocalDtAcousticSolver.cpp
@@ -16,6 +16,7 @@
 #include <scheme/IDiscreteFunctionDescriptor.hpp>
 #include <scheme/SymmetryBoundaryConditionDescriptor.hpp>
 #include <scheme/CouplingBoundaryConditionDescriptor.hpp>
+#include <scheme/AcousticSolver.hpp>
 #include <utils/Socket.hpp>
 
 #include <variant>
@@ -311,6 +312,12 @@ class LocalDtAcousticSolverHandler::LocalDtAcousticSolver final
 			NodeValuePerCell<Rd>& Fjr,
 			NodeValuePerCell<Rd>& CR_Fjr,
 			const std::vector<NodeId>& map2) const;
+  void _applyCouplingBC2(NodeValue<Rdxd>& Ar1,
+			NodeValue<Rdxd>& Ar2,
+			NodeValue<Rd>& ur1,
+			NodeValue<Rd>& ur2,
+			const std::vector<NodeId>& map1,
+			const std::vector<NodeId>& map2) const;
 
   void
   _applyBoundaryConditions(const BoundaryConditionList& bc_list,
@@ -590,6 +597,91 @@ class LocalDtAcousticSolverHandler::LocalDtAcousticSolver final
                            std::make_shared<const SubItemValuePerItemVariant>(Fjr));
   }
 
+    std::tuple<const std::shared_ptr<const ItemValueVariant>,
+	     const std::shared_ptr<const SubItemValuePerItemVariant>,
+	     const std::shared_ptr<const ItemValueVariant>,
+	     const std::shared_ptr<const SubItemValuePerItemVariant>,
+	     NodeValue<Rd>,
+	     NodeValuePerCell<Rd>>
+  compute_fluxes2(const SolverType& solver_type,
+                 const std::shared_ptr<const DiscreteFunctionVariant>& rho1_v,
+                 const std::shared_ptr<const DiscreteFunctionVariant>& c1_v,
+                 const std::shared_ptr<const DiscreteFunctionVariant>& u1_v,
+                 const std::shared_ptr<const DiscreteFunctionVariant>& p1_v,
+                 const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>& bc_descriptor_list1,
+                 const std::shared_ptr<const DiscreteFunctionVariant>& rho2_v,
+                 const std::shared_ptr<const DiscreteFunctionVariant>& c2_v,
+                 const std::shared_ptr<const DiscreteFunctionVariant>& u2_v,
+                 const std::shared_ptr<const DiscreteFunctionVariant>& p2_v,
+                 const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>& bc_descriptor_list2,
+		 const std::vector<NodeId>& map1,
+		 const std::vector<NodeId>& map2) const
+  {
+    std::shared_ptr i_mesh1 = getCommonMesh({rho1_v, c1_v, u1_v, p1_v});
+    std::shared_ptr i_mesh2 = getCommonMesh({rho2_v, c2_v, u2_v, p2_v});
+    if (not i_mesh1) {
+      throw NormalError("discrete functions are not defined on the same mesh");
+    }
+
+    if (not checkDiscretizationType({rho1_v, c1_v, u1_v, p1_v}, DiscreteFunctionType::P0)) {
+      throw NormalError("acoustic solver expects P0 functions");
+    }
+
+    if (not i_mesh2) {
+      throw NormalError("discrete functions are not defined on the same mesh");
+    }
+
+    if (not checkDiscretizationType({rho2_v, c2_v, u2_v, p2_v}, DiscreteFunctionType::P0)) {
+      throw NormalError("acoustic solver expects P0 functions");
+    }
+
+
+    const MeshType& mesh1              = dynamic_cast<const MeshType&>(*i_mesh1);
+    const DiscreteScalarFunction& rho1 = rho1_v->get<DiscreteScalarFunction>();
+    const DiscreteScalarFunction& c1   = c1_v->get<DiscreteScalarFunction>();
+    const DiscreteVectorFunction& u1   = u1_v->get<DiscreteVectorFunction>();
+    const DiscreteScalarFunction& p1   = p1_v->get<DiscreteScalarFunction>();
+
+    const MeshType& mesh2              = dynamic_cast<const MeshType&>(*i_mesh2);
+    const DiscreteScalarFunction& rho2 = rho2_v->get<DiscreteScalarFunction>();
+    const DiscreteScalarFunction& c2   = c2_v->get<DiscreteScalarFunction>();
+    const DiscreteVectorFunction& u2   = u2_v->get<DiscreteVectorFunction>();
+    const DiscreteScalarFunction& p2   = p2_v->get<DiscreteScalarFunction>();
+
+    NodeValuePerCell<const Rdxd> Ajr1 = this->_computeAjr(solver_type, mesh1, rho1 * c1);
+    NodeValuePerCell<const Rdxd> Ajr2 = this->_computeAjr(solver_type, mesh2, rho2 * c2);
+
+    NodeValue<Rdxd> Ar1 = this->_computeAr(mesh1, Ajr1);
+    NodeValue<Rd> br1   = this->_computeBr(mesh1, Ajr1, u1, p1);
+    NodeValue<Rdxd> Ar2 = this->_computeAr(mesh2, Ajr2);
+    NodeValue<Rd> br2   = this->_computeBr(mesh2, Ajr2, u2, p2);
+
+    const BoundaryConditionList bc_list1 = this->_getBCList(mesh1, bc_descriptor_list1);
+    const BoundaryConditionList bc_list2 = this->_getBCList(mesh2, bc_descriptor_list2);
+    this->_applyBoundaryConditions(bc_list1, mesh1, Ar1, br1);
+    this->_applyBoundaryConditions(bc_list2, mesh2, Ar2, br2);
+    this->_applyExternalVelocityBC(bc_list1, p1, Ar1, br1);
+    this->_applyExternalVelocityBC(bc_list2, p2, Ar2, br2);
+
+    synchronize(Ar1);
+    synchronize(br1);
+    synchronize(Ar2);
+    synchronize(br2);
+
+    NodeValue<Rd> ur1         = this->_computeUr(mesh1, Ar1, br1);
+    NodeValue<Rd> ur2         = this->_computeUr(mesh2, Ar2, br2);
+    this->_applyCouplingBC2(Ar1,Ar2,ur1,ur2,map1,map2);
+    NodeValuePerCell<Rd> Fjr1 = this->_computeFjr(mesh1, Ajr1, ur1, u1, p1);
+    NodeValuePerCell<Rd> Fjr2 = this->_computeFjr(mesh2, Ajr2, ur2, u2, p2);
+
+    return std::make_tuple(std::make_shared<const ItemValueVariant>(ur1),
+                           std::make_shared<const SubItemValuePerItemVariant>(Fjr1),
+			   std::make_shared<const ItemValueVariant>(ur2),
+                           std::make_shared<const SubItemValuePerItemVariant>(Fjr2),
+			   ur2,
+                           Fjr2);
+  }
+
   std::tuple<std::shared_ptr<const IMesh>,
              std::shared_ptr<const DiscreteFunctionVariant>,
              std::shared_ptr<const DiscreteFunctionVariant>,
@@ -824,136 +916,780 @@ class LocalDtAcousticSolverHandler::LocalDtAcousticSolver final
     return {new_m1,new_rho1,new_u1,new_E1,new_m2,new_rho2,new_u2,new_E2};
   }
 
-  LocalDtAcousticSolver()                        = default;
-  LocalDtAcousticSolver(LocalDtAcousticSolver&&) = default;
-  ~LocalDtAcousticSolver()                       = default;
-};
+  std::tuple<std::shared_ptr<const IMesh>,
+             std::shared_ptr<const DiscreteFunctionVariant>,
+             std::shared_ptr<const DiscreteFunctionVariant>,
+             std::shared_ptr<const DiscreteFunctionVariant>,
+	     std::shared_ptr<const IMesh>,
+	     std::shared_ptr<const DiscreteFunctionVariant>,
+	     std::shared_ptr<const DiscreteFunctionVariant>,
+	     std::shared_ptr<const DiscreteFunctionVariant>>
+  apply(const SolverType& solver_type,
+        const double& dt1,                                     
+        const std::shared_ptr<const DiscreteFunctionVariant>& rho1,
+	const std::shared_ptr<const DiscreteFunctionVariant>& rho2,
+        const std::shared_ptr<const DiscreteFunctionVariant>& u1,
+	const std::shared_ptr<const DiscreteFunctionVariant>& u2,
+        const std::shared_ptr<const DiscreteFunctionVariant>& E1,
+	const std::shared_ptr<const DiscreteFunctionVariant>& E2,
+        const std::shared_ptr<const DiscreteFunctionVariant>& c1,
+	const std::shared_ptr<const DiscreteFunctionVariant>& c2,
+        const std::shared_ptr<const DiscreteFunctionVariant>& p1,
+	const std::shared_ptr<const DiscreteFunctionVariant>& p2,
+        const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>& bc_descriptor_list1,
+	const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>& bc_descriptor_list2) const
+  {
+    
 
-template <size_t Dimension>
-void
-LocalDtAcousticSolverHandler::LocalDtAcousticSolver<Dimension>::_applyPressureBC(const BoundaryConditionList& bc_list,
-                                                                                 const MeshType& mesh,
-                                                                                 NodeValue<Rd>& br) const
-{
-  for (const auto& boundary_condition : bc_list) {
-    std::visit(
-      [&](auto&& bc) {
-        using T = std::decay_t<decltype(bc)>;
-        if constexpr (std::is_same_v<PressureBoundaryCondition, T>) {
-          MeshData<Dimension>& mesh_data = MeshDataManager::instance().getMeshData(mesh);
-          if constexpr (Dimension == 1) {
-            const NodeValuePerCell<const Rd> Cjr = mesh_data.Cjr();
+    const double& gamma = 1.4;
 
-            const auto& node_to_cell_matrix               = mesh.connectivity().nodeToCellMatrix();
-            const auto& node_local_numbers_in_their_cells = mesh.connectivity().nodeLocalNumbersInTheirCells();
+    std::shared_ptr<const DiscreteFunctionVariant> new_rho2 = rho2;
+    std::shared_ptr<const DiscreteFunctionVariant> new_u2   = u2;
+    std::shared_ptr<const DiscreteFunctionVariant> new_E2   = E2;
+    std::shared_ptr<const DiscreteFunctionVariant> new_c2   = c2;
+    std::shared_ptr<const DiscreteFunctionVariant> new_p2   = p2;
+    std::shared_ptr<const IMesh> new_m2 = getCommonMesh({new_rho2, new_c2, new_u2, new_p2});
+    std::shared_ptr<const IMesh> m1 = getCommonMesh({rho1, c1, u1, p1});
 
-            const auto& node_list  = bc.nodeList();
-            const auto& value_list = bc.valueList();
-            parallel_for(
-              node_list.size(), PUGS_LAMBDA(size_t i_node) {
-                const NodeId node_id       = node_list[i_node];
-                const auto& node_cell_list = node_to_cell_matrix[node_id];
-                Assert(node_cell_list.size() == 1);
+    double dt2 = 0.4 * acoustic_dt(new_c2);
 
-                CellId node_cell_id              = node_cell_list[0];
-                size_t node_local_number_in_cell = node_local_numbers_in_their_cells(node_id, 0);
+    auto [map1, map2] = _computeMapping(m1,new_m2,bc_descriptor_list1,bc_descriptor_list2);
 
-                br[node_id] -= value_list[i_node] * Cjr(node_cell_id, node_local_number_in_cell);
-              });
-          } else {
-            const NodeValuePerFace<const Rd> Nlr = mesh_data.Nlr();
+    auto [ur1, Fjr1, ur2, Fjr2,CR_ur,CR_Fjr] = compute_fluxes2(solver_type, rho1, c1, u1, p1, bc_descriptor_list1, new_rho2, new_c2, new_u2, new_p2, bc_descriptor_list2,map1,map2);
 
-            const auto& face_to_cell_matrix               = mesh.connectivity().faceToCellMatrix();
-            const auto& face_to_node_matrix               = mesh.connectivity().faceToNodeMatrix();
-            const auto& face_local_numbers_in_their_cells = mesh.connectivity().faceLocalNumbersInTheirCells();
-            const auto& face_cell_is_reversed             = mesh.connectivity().cellFaceIsReversed();
+    const auto [new_m1,new_rho1,new_u1,new_E1] = apply_fluxes(dt1, rho1, u1, E1, ur1, Fjr1);
+    std::tie(new_m2,new_rho2,new_u2,new_E2) = apply_fluxes(dt2,new_rho2,new_u2,new_E2,ur2,Fjr2);
+    
+    double sum_dt = dt2;
 
-            const auto& face_list  = bc.faceList();
-            const auto& value_list = bc.valueList();
-            for (size_t i_face = 0; i_face < face_list.size(); ++i_face) {
-              const FaceId face_id       = face_list[i_face];
-              const auto& face_cell_list = face_to_cell_matrix[face_id];
-              Assert(face_cell_list.size() == 1);
+    while(sum_dt < dt1){
 
-              CellId face_cell_id              = face_cell_list[0];
-              size_t face_local_number_in_cell = face_local_numbers_in_their_cells(face_id, 0);
+      const DiscreteScalarFunction& rho_d = new_rho2->get<DiscreteScalarFunction>();
+      const DiscreteVectorFunction& u_d   = new_u2->get<DiscreteVectorFunction>();
+      const DiscreteScalarFunction& E_d   = new_E2->get<DiscreteScalarFunction>();
+      const DiscreteScalarFunction& eps   = E_d - 0.5 * dot(u_d,u_d);
+      const DiscreteScalarFunction& p_d   = (gamma-1) * rho_d * eps;
+      const DiscreteScalarFunction& c_d   = sqrt(gamma * p_d / rho_d);
 
-              const double sign = face_cell_is_reversed(face_cell_id, face_local_number_in_cell) ? -1 : 1;
+      new_p2 = std::make_shared<DiscreteFunctionVariant>(p_d);
+      new_c2 = std::make_shared<DiscreteFunctionVariant>(c_d);
 
-              const auto& face_nodes = face_to_node_matrix[face_id];
+      double dt2 = 0.4 * acoustic_dt(new_c2);
+      //std::cout << "dt2 = " << dt2 << std::endl;
+      if(sum_dt + dt2 > dt1){
+	    dt2 = dt1 - sum_dt;
+      }
 
-              for (size_t i_node = 0; i_node < face_nodes.size(); ++i_node) {
-                NodeId node_id = face_nodes[i_node];
-                br[node_id] -= sign * value_list[i_face] * Nlr(face_id, i_node);
-              }
-            }
-          }
-        }
-      },
-      boundary_condition);
+      auto [ur2,Fjr2] = compute_fluxes(solver_type,new_rho2,new_c2,new_u2,new_p2,bc_descriptor_list2,CR_ur,CR_Fjr,map2); 
+      std::tie(new_m2,new_rho2,new_u2,new_E2) = apply_fluxes(dt2,new_rho2,new_u2,new_E2,ur2,Fjr2);
+      sum_dt += dt2;
+    }
+    std::tie(new_m2,new_rho2,new_u2,new_E2) = mesh_correction(new_m1,new_m2,new_rho2,new_u2,new_E2,map1,map2);
+    return {new_m1,new_rho1,new_u1,new_E1,new_m2,new_rho2,new_u2,new_E2};
   }
-}
 
-template <size_t Dimension>
-void
-LocalDtAcousticSolverHandler::LocalDtAcousticSolver<Dimension>::_applySymmetryBC(const BoundaryConditionList& bc_list,
-                                                                                 NodeValue<Rdxd>& Ar,
-                                                                                 NodeValue<Rd>& br) const
-{
-  for (const auto& boundary_condition : bc_list) {
-    std::visit(
-      [&](auto&& bc) {
-        using T = std::decay_t<decltype(bc)>;
-        if constexpr (std::is_same_v<SymmetryBoundaryCondition, T>) {
-          const Rd& n = bc.outgoingNormal();
+  //************************************************
+  //**** Nouvelles fcts pour l'ordre 2 en temps ****
+  //************************************************
 
-          const Rdxd I   = identity;
-          const Rdxd nxn = tensorProduct(n, n);
-          const Rdxd P   = I - nxn;
+  std::tuple<std::shared_ptr<const IMesh>,
+             std::shared_ptr<const DiscreteFunctionVariant>,
+             std::shared_ptr<const DiscreteFunctionVariant>,
+	     std::shared_ptr<const DiscreteFunctionVariant>>
+  apply_fluxes(const double& dt,
+               const double& tk, 
+               const MeshType& mesh,
+               const DiscreteScalarFunction& rho,
+               const DiscreteVectorFunction& u,
+               const DiscreteScalarFunction& E,
+               const NodeValue<const Rd>& ur_n,
+               const NodeValuePerCell<const Rd>& Fjr_n,
+               const NodeValue<const Rd>& ur_app,
+               const NodeValuePerCell<const Rd>& Fjr_app) const
+  {
+    const auto& cell_to_node_matrix = mesh.connectivity().cellToNodeMatrix();
 
-          const Array<const NodeId>& node_list = bc.nodeList();
-          parallel_for(
-            bc.numberOfNodes(), PUGS_LAMBDA(int r_number) {
-              const NodeId r = node_list[r_number];
+    if ((mesh.shared_connectivity() != ur_n.connectivity_ptr()) or
+        (mesh.shared_connectivity() != Fjr_n.connectivity_ptr())) {
+    //  throw NormalError("fluxes are not defined on the same connectivity than the mesh");
+    }
 
-              Ar[r] = P * Ar[r] * P + nxn;
-              br[r] = P * br[r];
-            });
+    NodeValue<Rd> ur{mesh.connectivity()};
+    parallel_for(
+      mesh.numberOfNodes(), PUGS_LAMBDA(NodeId r){
+
+        for(size_t i = 0; i<Dimension; ++i){
+          ur[r][i] = (1-tk)*ur_n[r][i] + tk*ur_app[r][i];
         }
-      },
-      boundary_condition);
-  }
-}
+      }
+    );
 
-template <size_t Dimension>
-void
-LocalDtAcousticSolverHandler::LocalDtAcousticSolver<Dimension>::_applyVelocityBC(const BoundaryConditionList& bc_list,
-                                                                                 NodeValue<Rdxd>& Ar,
-                                                                                 NodeValue<Rd>& br) const
-{
-  for (const auto& boundary_condition : bc_list) {
-    std::visit(
-      [&](auto&& bc) {
-        using T = std::decay_t<decltype(bc)>;
-        if constexpr (std::is_same_v<VelocityBoundaryCondition, T>) {
-          const auto& node_list  = bc.nodeList();
-          const auto& value_list = bc.valueList();
+    NodeValuePerCell<Rd> Fjr{mesh.connectivity()};
+    parallel_for(
+      mesh.numberOfCells(), PUGS_LAMBDA(CellId j){
+        const auto& cell_nodes = cell_to_node_matrix[j];
 
-          parallel_for(
-            node_list.size(), PUGS_LAMBDA(size_t i_node) {
-              NodeId node_id    = node_list[i_node];
-              const auto& value = value_list[i_node];
+        for (size_t r = 0; r < cell_nodes.size(); ++r){
+          for(size_t i = 0; i<Dimension; i++){
+            Fjr(j,r)[i] = (1-tk)*Fjr_n(j,r)[i] + tk*Fjr_app(j,r)[i];
+          }
+        }
+      }
+    );
 
-              Ar[node_id] = identity;
-              br[node_id] = value;
-            });
-        } else if constexpr (std::is_same_v<FixedBoundaryCondition, T>) {
-          const auto& node_list = bc.nodeList();
-          parallel_for(
-            node_list.size(), PUGS_LAMBDA(size_t i_node) {
-              NodeId node_id = node_list[i_node];
+    NodeValue<Rd> new_xr = copy(mesh.xr());
+    parallel_for(
+      mesh.numberOfNodes(), PUGS_LAMBDA(NodeId r) { new_xr[r] += dt * ur[r]; });
 
-              Ar[node_id] = identity;
+    std::shared_ptr<const MeshType> new_mesh = std::make_shared<MeshType>(mesh.shared_connectivity(), new_xr);
+
+    CellValue<const double> Vj = MeshDataManager::instance().getMeshData(mesh).Vj();
+
+    CellValue<double> new_rho = copy(rho.cellValues());
+    CellValue<Rd> new_u       = copy(u.cellValues());
+    CellValue<double> new_E   = copy(E.cellValues());
+
+    parallel_for(
+      mesh.numberOfCells(), PUGS_LAMBDA(CellId j) {
+        const auto& cell_nodes = cell_to_node_matrix[j];
+
+        Rd momentum_fluxes   = zero;
+        double energy_fluxes = 0 ;
+        //double energy_fluxes_n = 0;
+        //double energy_fluxes_app = 0;
+        for (size_t R = 0; R < cell_nodes.size(); ++R) {
+          const NodeId r = cell_nodes[R];
+          momentum_fluxes += Fjr(j, R);
+          energy_fluxes   += ((1-tk)*dot(Fjr_n(j, R), ur_n[r]) + tk*dot(Fjr_app(j,R),ur_app[r]));
+          //energy_fluxes_n += dot(Fjr_n(j, R), ur_n[r]);
+          //energy_fluxes_app += dot(Fjr_app(j, R), ur_app[r]);
+        }
+        //double energy_fluxes = (1-tk)*energy_fluxes_n + tk*energy_fluxes_app;
+        const double dt_over_Mj =  dt / (rho[j] * Vj[j]);
+        new_u[j] -= dt_over_Mj * momentum_fluxes;
+        new_E[j] -= dt_over_Mj * energy_fluxes;
+      });
+
+    CellValue<const double> new_Vj = MeshDataManager::instance().getMeshData(*new_mesh).Vj();
+
+    parallel_for(
+      mesh.numberOfCells(), PUGS_LAMBDA(CellId j) { new_rho[j] *= Vj[j] / new_Vj[j]; });
+
+    return {new_mesh, std::make_shared<DiscreteFunctionVariant>(DiscreteScalarFunction(new_mesh, new_rho)),
+            std::make_shared<DiscreteFunctionVariant>(DiscreteVectorFunction(new_mesh, new_u)),
+            std::make_shared<DiscreteFunctionVariant>(DiscreteScalarFunction(new_mesh, new_E))};
+  }
+
+  std::tuple<std::shared_ptr<const IMesh>,
+             std::shared_ptr<const DiscreteFunctionVariant>,
+             std::shared_ptr<const DiscreteFunctionVariant>,
+	     std::shared_ptr<const DiscreteFunctionVariant>>
+  apply_fluxes(const double& dt,
+               const MeshType& mesh,
+               const DiscreteScalarFunction& rho,
+               const DiscreteVectorFunction& u,
+               const DiscreteScalarFunction& E,
+               const NodeValue<const Rd>& ur_n,
+               const NodeValuePerCell<const Rd>& Fjr_n,
+               const NodeValue<const Rd>& ur_app,
+               const NodeValuePerCell<const Rd>& Fjr_app) const
+  {
+    const auto& cell_to_node_matrix = mesh.connectivity().cellToNodeMatrix();
+
+    if ((mesh.shared_connectivity() != ur_n.connectivity_ptr()) or
+        (mesh.shared_connectivity() != Fjr_n.connectivity_ptr())) {
+    //  throw NormalError("fluxes are not defined on the same connectivity than the mesh");
+    }
+
+    NodeValue<Rd> ur{mesh.connectivity()};
+    parallel_for(
+      mesh.numberOfNodes(), PUGS_LAMBDA(NodeId r){
+        ur[r] = 0.5*(ur_n[r] + ur_app[r]);
+      }
+    );
+
+    NodeValuePerCell<Rd> Fjr{mesh.connectivity()};
+    parallel_for(
+      mesh.numberOfCells(), PUGS_LAMBDA(CellId j){
+        const auto& cell_nodes = cell_to_node_matrix[j];
+
+        for (size_t r = 0; r < cell_nodes.size(); ++r){
+          Fjr(j,r) = 0.5*(Fjr_n(j,r) + Fjr_app(j,r));
+        }
+      }
+    );
+
+    NodeValue<Rd> new_xr = copy(mesh.xr());
+    parallel_for(
+      mesh.numberOfNodes(), PUGS_LAMBDA(NodeId r) { new_xr[r] += dt * ur[r]; });
+
+    std::shared_ptr<const MeshType> new_mesh = std::make_shared<MeshType>(mesh.shared_connectivity(), new_xr);
+
+    CellValue<const double> Vj = MeshDataManager::instance().getMeshData(mesh).Vj();
+
+    CellValue<double> new_rho = copy(rho.cellValues());
+    CellValue<Rd> new_u       = copy(u.cellValues());
+    CellValue<double> new_E   = copy(E.cellValues());
+
+    parallel_for(
+      mesh.numberOfCells(), PUGS_LAMBDA(CellId j) {
+        const auto& cell_nodes = cell_to_node_matrix[j];
+
+        Rd momentum_fluxes   = zero;
+        double energy_fluxes = 0 ;
+        for (size_t R = 0; R < cell_nodes.size(); ++R) {
+          const NodeId r = cell_nodes[R];
+          momentum_fluxes += Fjr(j, R);
+          energy_fluxes += 0.5*dot(Fjr_n(j,R),ur_n[r]) + 0.5*dot(Fjr_app(j,R),ur_app[r]);
+        }
+        const double dt_over_Mj =  dt / (rho[j] * Vj[j]);
+        new_u[j] -= dt_over_Mj * momentum_fluxes;
+        new_E[j] -= dt_over_Mj * energy_fluxes;
+      });
+
+    CellValue<const double> new_Vj = MeshDataManager::instance().getMeshData(*new_mesh).Vj();
+
+    parallel_for(
+      mesh.numberOfCells(), PUGS_LAMBDA(CellId j) { new_rho[j] *= Vj[j] / new_Vj[j]; });
+
+    return {new_mesh, std::make_shared<DiscreteFunctionVariant>(DiscreteScalarFunction(new_mesh, new_rho)),
+            std::make_shared<DiscreteFunctionVariant>(DiscreteVectorFunction(new_mesh, new_u)),
+            std::make_shared<DiscreteFunctionVariant>(DiscreteScalarFunction(new_mesh, new_E))};
+  }
+
+  std::tuple<std::shared_ptr<const IMesh>,
+             std::shared_ptr<const DiscreteFunctionVariant>,
+             std::shared_ptr<const DiscreteFunctionVariant>,
+             std::shared_ptr<const DiscreteFunctionVariant>>
+  apply_fluxes(const double& dt,
+               const double& tk,
+               const std::shared_ptr<const DiscreteFunctionVariant>& rho_v,
+               const std::shared_ptr<const DiscreteFunctionVariant>& u_v,
+               const std::shared_ptr<const DiscreteFunctionVariant>& E_v,
+               const std::shared_ptr<const ItemValueVariant>& ur_n,
+               const std::shared_ptr<const SubItemValuePerItemVariant>& Fjr_n,
+               const std::shared_ptr<const ItemValueVariant>& ur_app,
+               const std::shared_ptr<const SubItemValuePerItemVariant>& Fjr_app) const
+  {
+    std::shared_ptr i_mesh = getCommonMesh({rho_v, u_v, E_v});
+    if (not i_mesh) {
+      throw NormalError("discrete functions are not defined on the same mesh");
+    }
+
+    if (not checkDiscretizationType({rho_v, u_v, E_v}, DiscreteFunctionType::P0)) {
+      throw NormalError("acoustic solver expects P0 functions");
+    }
+
+
+    return this->apply_fluxes(dt,
+                              tk,                                       //
+                              dynamic_cast<const MeshType&>(*i_mesh),   //
+                              rho_v->get<DiscreteScalarFunction>(),     //
+                              u_v->get<DiscreteVectorFunction>(),       //
+                              E_v->get<DiscreteScalarFunction>(),       //
+                              ur_n->get<NodeValue<const Rd>>(),           //
+                              Fjr_n->get<NodeValuePerCell<const Rd>>(),
+                              ur_app->get<NodeValue<const Rd>>(),           //
+                              Fjr_app->get<NodeValuePerCell<const Rd>>());
+  }
+
+  std::tuple<std::shared_ptr<const IMesh>,
+             std::shared_ptr<const DiscreteFunctionVariant>,
+             std::shared_ptr<const DiscreteFunctionVariant>,
+             std::shared_ptr<const DiscreteFunctionVariant>>
+  apply_fluxes(const double& dt,
+               const std::shared_ptr<const DiscreteFunctionVariant>& rho_v,
+               const std::shared_ptr<const DiscreteFunctionVariant>& u_v,
+               const std::shared_ptr<const DiscreteFunctionVariant>& E_v,
+               const std::shared_ptr<const ItemValueVariant>& ur_n,
+               const std::shared_ptr<const SubItemValuePerItemVariant>& Fjr_n,
+               const std::shared_ptr<const ItemValueVariant>& ur_app,
+               const std::shared_ptr<const SubItemValuePerItemVariant>& Fjr_app) const
+  {
+    std::shared_ptr i_mesh = getCommonMesh({rho_v, u_v, E_v});
+    if (not i_mesh) {
+      throw NormalError("discrete functions are not defined on the same mesh");
+    }
+
+    if (not checkDiscretizationType({rho_v, u_v, E_v}, DiscreteFunctionType::P0)) {
+      throw NormalError("acoustic solver expects P0 functions");
+    }
+
+
+    return this->apply_fluxes(dt,                                       //
+                              dynamic_cast<const MeshType&>(*i_mesh),   //
+                              rho_v->get<DiscreteScalarFunction>(),     //
+                              u_v->get<DiscreteVectorFunction>(),       //
+                              E_v->get<DiscreteScalarFunction>(),       //
+                              ur_n->get<NodeValue<const Rd>>(),           //
+                              Fjr_n->get<NodeValuePerCell<const Rd>>(),
+                              ur_app->get<NodeValue<const Rd>>(),           //
+                              Fjr_app->get<NodeValuePerCell<const Rd>>());
+  }
+
+  std::tuple<const std::shared_ptr<const ItemValueVariant>,
+	     const std::shared_ptr<const SubItemValuePerItemVariant>,
+	     const std::shared_ptr<const ItemValueVariant>,
+	     const std::shared_ptr<const SubItemValuePerItemVariant>>
+  compute_fluxes2(const SolverType& solver_type,
+                 const std::shared_ptr<const DiscreteFunctionVariant>& rho1_v,
+                 const std::shared_ptr<const DiscreteFunctionVariant>& c1_v,
+                 const std::shared_ptr<const DiscreteFunctionVariant>& u1_v,
+                 const std::shared_ptr<const DiscreteFunctionVariant>& p1_v,
+                 const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>& bc_descriptor_list1,
+                 const std::shared_ptr<const DiscreteFunctionVariant>& rho2_v,
+                 const std::shared_ptr<const DiscreteFunctionVariant>& c2_v,
+                 const std::shared_ptr<const DiscreteFunctionVariant>& u2_v,
+                 const std::shared_ptr<const DiscreteFunctionVariant>& p2_v,
+                 const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>& bc_descriptor_list2) const
+  {
+    std::shared_ptr i_mesh1 = getCommonMesh({rho1_v, c1_v, u1_v, p1_v});
+    std::shared_ptr i_mesh2 = getCommonMesh({rho2_v, c2_v, u2_v, p2_v});
+    if (not i_mesh1) {
+      throw NormalError("discrete functions are not defined on the same mesh");
+    }
+
+    if (not checkDiscretizationType({rho1_v, c1_v, u1_v, p1_v}, DiscreteFunctionType::P0)) {
+      throw NormalError("acoustic solver expects P0 functions");
+    }
+
+    if (not i_mesh2) {
+      throw NormalError("discrete functions are not defined on the same mesh");
+    }
+
+    if (not checkDiscretizationType({rho2_v, c2_v, u2_v, p2_v}, DiscreteFunctionType::P0)) {
+      throw NormalError("acoustic solver expects P0 functions");
+    }
+
+    auto [map1,map2] = _computeMapping(i_mesh1,i_mesh2,bc_descriptor_list1,bc_descriptor_list2);
+
+    const MeshType& mesh1              = dynamic_cast<const MeshType&>(*i_mesh1);
+    const DiscreteScalarFunction& rho1 = rho1_v->get<DiscreteScalarFunction>();
+    const DiscreteScalarFunction& c1   = c1_v->get<DiscreteScalarFunction>();
+    const DiscreteVectorFunction& u1   = u1_v->get<DiscreteVectorFunction>();
+    const DiscreteScalarFunction& p1   = p1_v->get<DiscreteScalarFunction>();
+
+    const MeshType& mesh2              = dynamic_cast<const MeshType&>(*i_mesh2);
+    const DiscreteScalarFunction& rho2 = rho2_v->get<DiscreteScalarFunction>();
+    const DiscreteScalarFunction& c2   = c2_v->get<DiscreteScalarFunction>();
+    const DiscreteVectorFunction& u2   = u2_v->get<DiscreteVectorFunction>();
+    const DiscreteScalarFunction& p2   = p2_v->get<DiscreteScalarFunction>();
+
+    NodeValuePerCell<const Rdxd> Ajr1 = this->_computeAjr(solver_type, mesh1, rho1 * c1);
+    NodeValuePerCell<const Rdxd> Ajr2 = this->_computeAjr(solver_type, mesh2, rho2 * c2);
+
+    NodeValue<Rdxd> Ar1 = this->_computeAr(mesh1, Ajr1);
+    NodeValue<Rd> br1   = this->_computeBr(mesh1, Ajr1, u1, p1);
+    NodeValue<Rdxd> Ar2 = this->_computeAr(mesh2, Ajr2);
+    NodeValue<Rd> br2   = this->_computeBr(mesh2, Ajr2, u2, p2);
+
+    const BoundaryConditionList bc_list1 = this->_getBCList(mesh1, bc_descriptor_list1);
+    const BoundaryConditionList bc_list2 = this->_getBCList(mesh2, bc_descriptor_list2);
+    this->_applyBoundaryConditions(bc_list1, mesh1, Ar1, br1);
+    this->_applyBoundaryConditions(bc_list2, mesh2, Ar2, br2);
+    this->_applyExternalVelocityBC(bc_list1, p1, Ar1, br1);
+    this->_applyExternalVelocityBC(bc_list2, p2, Ar2, br2);
+
+    synchronize(Ar1);
+    synchronize(br1);
+    synchronize(Ar2);
+    synchronize(br2);
+
+    NodeValue<Rd> ur1         = this->_computeUr(mesh1, Ar1, br1);
+    NodeValue<Rd> ur2         = this->_computeUr(mesh2, Ar2, br2);
+    this->_applyCouplingBC2(Ar1,Ar2,ur1,ur2,map1,map2);
+    NodeValuePerCell<Rd> Fjr1 = this->_computeFjr(mesh1, Ajr1, ur1, u1, p1);
+    NodeValuePerCell<Rd> Fjr2 = this->_computeFjr(mesh2, Ajr2, ur2, u2, p2);
+
+    return std::make_tuple(std::make_shared<const ItemValueVariant>(ur1),
+                           std::make_shared<const SubItemValuePerItemVariant>(Fjr1),
+			   std::make_shared<const ItemValueVariant>(ur2),
+                           std::make_shared<const SubItemValuePerItemVariant>(Fjr2));
+  }
+
+  std::tuple<const std::shared_ptr<const ItemValueVariant>, const std::shared_ptr<const SubItemValuePerItemVariant>>
+  compute_fluxes(const SolverType& solver_type,
+                 const std::shared_ptr<const DiscreteFunctionVariant>& rho_v,
+                 const std::shared_ptr<const DiscreteFunctionVariant>& c_v,
+                 const std::shared_ptr<const DiscreteFunctionVariant>& u_v,
+                 const std::shared_ptr<const DiscreteFunctionVariant>& p_v,
+                 const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>& bc_descriptor_list,
+		             const std::shared_ptr<const ItemValueVariant>& CR_ur,
+                 const std::shared_ptr<const SubItemValuePerItemVariant>& CR_Fjr) const 
+  {
+    std::shared_ptr i_mesh = getCommonMesh({rho_v, c_v, u_v, p_v});
+    if (not i_mesh) {
+      throw NormalError("discrete functions are not defined on the same mesh");
+    }
+
+    if (not checkDiscretizationType({rho_v, c_v, u_v, p_v}, DiscreteFunctionType::P0)) {
+      throw NormalError("acoustic solver expects P0 functions");
+    }
+
+    const MeshType& mesh              = dynamic_cast<const MeshType&>(*i_mesh);
+    const DiscreteScalarFunction& rho = rho_v->get<DiscreteScalarFunction>();
+    const DiscreteScalarFunction& c   = c_v->get<DiscreteScalarFunction>();
+    const DiscreteVectorFunction& u   = u_v->get<DiscreteVectorFunction>();
+    const DiscreteScalarFunction& p   = p_v->get<DiscreteScalarFunction>();
+
+    std::vector<NodeId> map2;
+    const BoundaryConditionList bc_list2 = this->_getBCList(mesh,bc_descriptor_list);
+
+    for(const auto& boundary_condition2 : bc_list2){
+      std::visit([&](auto&& bc2){
+        using T2 = std::decay_t<decltype(bc2)>;
+        if constexpr (std::is_same_v<CouplingBoundaryCondition,T2>){
+          const auto& node_list2 = bc2.nodeList();
+          for(size_t i = 0; i<node_list2.size(); i++){
+            map2.push_back(node_list2[i]);
+          }
+        }
+      },boundary_condition2);
+    }
+
+    NodeValuePerCell<const Rdxd> Ajr = this->_computeAjr(solver_type, mesh, rho * c);
+
+    NodeValue<Rdxd> Ar = this->_computeAr(mesh, Ajr);
+    NodeValue<Rd> br   = this->_computeBr(mesh, Ajr, u, p);
+    
+    const BoundaryConditionList bc_list = this->_getBCList(mesh, bc_descriptor_list);
+    this->_applyBoundaryConditions(bc_list, mesh, Ar, br);
+    this->_applyExternalVelocityBC(bc_list, p, Ar, br);
+
+    synchronize(Ar);
+    synchronize(br);
+
+    NodeValue<Rd> ur         = this->_computeUr(mesh, Ar, br);
+    NodeValuePerCell<Rd> Fjr = this->_computeFjr(mesh, Ajr, ur, u, p);
+
+    NodeValue<const Rd> ur_c          = CR_ur->get<NodeValue<const Rd>>();
+    NodeValuePerCell<const Rd>  Fjr_c = CR_Fjr->get<NodeValuePerCell<const Rd>>();
+
+    const size_t& n = map2.size();
+
+  const auto& node_to_cell_matrix               = mesh.connectivity().nodeToCellMatrix();
+  const auto& node_local_numbers_in_their_cells = mesh.connectivity().nodeLocalNumbersInTheirCells();
+
+  for(size_t i = 0; i<n; i++){
+    
+    const NodeId node_id                       = map2[i];
+    const auto& node_to_cell                   = node_to_cell_matrix[node_id];
+    const auto& node_local_number_in_its_cells = node_local_numbers_in_their_cells.itemArray(node_id);
+
+    for(size_t j = 0; j < node_to_cell.size(); j++){
+      const CellId J       = node_to_cell[j];
+      const unsigned int R = node_local_number_in_its_cells[j];
+      Fjr(J,R)             = Fjr_c(J,R);
+    }
+    ur[node_id] = ur_c[node_id];
+   }
+    
+    //this->_applyCouplingBC(mesh,ur,ur_c,Fjr,Fjr_c,map2);
+    //this->_applyCouplingBC(mesh,ur,CR_ur,Fjr,CR_Fjr,map2);
+
+    return std::make_tuple(std::make_shared<const ItemValueVariant>(ur),
+                           std::make_shared<const SubItemValuePerItemVariant>(Fjr));
+  }
+
+  std::tuple<std::shared_ptr<const IMesh>,
+             std::shared_ptr<const DiscreteFunctionVariant>,
+             std::shared_ptr<const DiscreteFunctionVariant>,
+             std::shared_ptr<const DiscreteFunctionVariant>,
+	     std::shared_ptr<const IMesh>,
+	     std::shared_ptr<const DiscreteFunctionVariant>,
+	     std::shared_ptr<const DiscreteFunctionVariant>,
+	     std::shared_ptr<const DiscreteFunctionVariant>>
+  applyOrder2(const SolverType& solver_type,
+        const double& dt1,                                     
+        const std::shared_ptr<const DiscreteFunctionVariant>& rho1,
+	const std::shared_ptr<const DiscreteFunctionVariant>& rho2,
+        const std::shared_ptr<const DiscreteFunctionVariant>& u1,
+	const std::shared_ptr<const DiscreteFunctionVariant>& u2,
+        const std::shared_ptr<const DiscreteFunctionVariant>& E1,
+	const std::shared_ptr<const DiscreteFunctionVariant>& E2,
+        const std::shared_ptr<const DiscreteFunctionVariant>& c1,
+	const std::shared_ptr<const DiscreteFunctionVariant>& c2,
+        const std::shared_ptr<const DiscreteFunctionVariant>& p1,
+	const std::shared_ptr<const DiscreteFunctionVariant>& p2,
+        const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>& bc_descriptor_list1,
+	const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>& bc_descriptor_list2) const
+  {
+    
+
+    const double& gamma = 1.4;
+
+    std::shared_ptr<const DiscreteFunctionVariant> new_rho2 = rho2;
+    std::shared_ptr<const DiscreteFunctionVariant> new_u2   = u2;
+    std::shared_ptr<const DiscreteFunctionVariant> new_E2   = E2;
+    std::shared_ptr<const DiscreteFunctionVariant> new_c2   = c2;
+    std::shared_ptr<const DiscreteFunctionVariant> new_p2   = p2;
+    std::shared_ptr<const IMesh> new_m2 = getCommonMesh({new_rho2, new_c2, new_u2, new_p2});
+    std::shared_ptr<const IMesh> m1 = getCommonMesh({rho1, c1, u1, p1});
+
+    const MeshType& mesh1 = dynamic_cast<const MeshType&>(*m1);
+    const MeshType& mesh2  = dynamic_cast<const MeshType&>(*new_m2);
+
+    //double dt2 = 0.4 * acoustic_dt(new_c2);
+    //double dt2_k = dt1/2;
+    double dt2_k = dt1/2;
+
+    auto [map1, map2] = _computeMapping(m1,new_m2,bc_descriptor_list1,bc_descriptor_list2);
+
+    const auto [ur1_n, Fjr1_n, ur2_n, Fjr2_n,CR_ur,CR_Fjr] = compute_fluxes2(solver_type, rho1, c1, u1, p1, bc_descriptor_list1, new_rho2, new_c2, new_u2, new_p2, bc_descriptor_list2,map1,map2);
+
+
+    std::shared_ptr<const DiscreteFunctionVariant> rho1_app;
+    std::shared_ptr<const DiscreteFunctionVariant> u1_app;
+    std::shared_ptr<const DiscreteFunctionVariant> E1_app;
+    std::shared_ptr<const DiscreteFunctionVariant> rho2_app;
+    std::shared_ptr<const DiscreteFunctionVariant> u2_app;
+    std::shared_ptr<const DiscreteFunctionVariant> E2_app;
+    std::shared_ptr<const IMesh> m1_app;
+    std::shared_ptr<const IMesh> m2_app;
+
+    std::tie(m1_app,rho1_app,u1_app,E1_app) = apply_fluxes(dt1, rho1, u1, E1, ur1_n, Fjr1_n);
+    std::tie(m2_app,rho2_app,u2_app,E2_app) = apply_fluxes(dt2_k, rho2, u2, E2, ur2_n, Fjr2_n);
+
+    double sum_dt = dt2_k;
+
+    while(sum_dt < dt1){
+
+      const DiscreteScalarFunction& rho_d = rho2_app->get<DiscreteScalarFunction>();
+      const DiscreteVectorFunction& u_d   = u2_app->get<DiscreteVectorFunction>();
+      const DiscreteScalarFunction& E_d   = E2_app->get<DiscreteScalarFunction>();
+      const DiscreteScalarFunction& eps   = E_d - 0.5 * dot(u_d,u_d);
+      const DiscreteScalarFunction& p_d   = (gamma-1) * rho_d * eps;
+      const DiscreteScalarFunction& c_d   = sqrt(gamma * p_d / rho_d);
+
+      const std::shared_ptr<const DiscreteFunctionVariant> p2_k = std::make_shared<DiscreteFunctionVariant>(p_d);
+      const std::shared_ptr<const DiscreteFunctionVariant> c2_k = std::make_shared<DiscreteFunctionVariant>(c_d);
+
+      //double dt2 = 0.4 * acoustic_dt(new_c2);
+      dt2_k = dt1/2;
+      //std::cout << "dt2 = " << dt2 << std::endl;
+      if(sum_dt + dt2_k > dt1){
+	    dt2_k = dt1 - sum_dt;
+      }
+
+      auto [ur_k,Fjr_k] = compute_fluxes(solver_type,rho2_app,c2_k,u2_app,p2_k,bc_descriptor_list2,CR_ur,CR_Fjr,map2); 
+      std::tie(m2_app,rho2_app,u2_app,E2_app) = apply_fluxes(dt2_k,rho2_app,u2_app,E2_app,ur_k,Fjr_k);
+      sum_dt += dt2_k;
+    }
+
+    const DiscreteScalarFunction& rho1_d = rho1_app->get<DiscreteScalarFunction>();
+    const DiscreteVectorFunction& u1_d   = u1_app->get<DiscreteVectorFunction>();
+    const DiscreteScalarFunction& E1_d   = E1_app->get<DiscreteScalarFunction>();
+    const DiscreteScalarFunction& eps1   = E1_d - 0.5 * dot(u1_d,u1_d);
+    const DiscreteScalarFunction& p1_d   = (gamma-1) * rho1_d * eps1;
+    const DiscreteScalarFunction& c1_d   = sqrt(gamma * p1_d / rho1_d);
+    const DiscreteScalarFunction& rho2_d = rho2_app->get<DiscreteScalarFunction>();
+    const DiscreteVectorFunction& u2_d   = u2_app->get<DiscreteVectorFunction>();
+    const DiscreteScalarFunction& E2_d   = E2_app->get<DiscreteScalarFunction>();
+    const DiscreteScalarFunction& eps2   = E2_d - 0.5 * dot(u2_d,u2_d);
+    const DiscreteScalarFunction& p2_d   = (gamma-1) * rho2_d * eps2;
+    const DiscreteScalarFunction& c2_d   = sqrt(gamma * p2_d / rho2_d);
+
+    const std::shared_ptr<const DiscreteFunctionVariant> c1_app = std::make_shared<const DiscreteFunctionVariant>(c1_d);
+    const std::shared_ptr<const DiscreteFunctionVariant> p1_app = std::make_shared<const DiscreteFunctionVariant>(p1_d);
+    const std::shared_ptr<const DiscreteFunctionVariant> c2_app = std::make_shared<const DiscreteFunctionVariant>(c2_d);
+    const std::shared_ptr<const DiscreteFunctionVariant> p2_app = std::make_shared<const DiscreteFunctionVariant>(p2_d);
+
+
+    const auto [ur1_app, Fjr1_app, ur2_app, Fjr2_app, CRapp1, CRapp2] = compute_fluxes2(solver_type, rho1_app, c1_app, u1_app, p1_app, bc_descriptor_list1, rho2_app, c2_app, u2_app, p2_app, bc_descriptor_list2,map1,map2);
+    //double dt2 = dt1/2;
+    double dt2 = dt1/2;
+    double tk = (dt2*0.5)/dt1;
+    //double tk = 1;
+
+    //std::shared_ptr<const ItemValueVariant> ur_inter =_InterpolateUr(mesh2,ur2_n,ur2_app,tk);
+    //std::shared_ptr<const SubItemValuePerItemVariant> Fjr_inter= _InterpolateFjr(mesh2,Fjr2_n,Fjr2_app,tk);
+    //std::shared_ptr<const ItemValueVariant> ur2 = _computeOrder2Ur(mesh2,ur_inter,ur2_n);
+    //std::shared_ptr<const SubItemValuePerItemVariant> Fjr2 = _computeOrder2Fjr(mesh2,Fjr_inter,Fjr2_n);
+    //std::shared_ptr<const ItemValueVariant> ur1 = _computeOrder2Ur(mesh1,ur1_app,ur1_n);
+    //std::shared_ptr<const SubItemValuePerItemVariant> Fjr1 = _computeOrder2Fjr(mesh1,Fjr1_app,Fjr1_n);
+
+    const auto [new_m1,new_rho1,new_u1,new_E1] = apply_fluxes(dt1,rho1,u1,E1,ur1_n,Fjr1_n,ur1_app,Fjr1_app);
+    //const auto [new_m1,new_rho1,new_u1,new_E1] = apply_fluxes(dt1,rho1,u1,E1,ur1,Fjr1);
+    //std::tie(new_m2,new_rho2,new_u2,new_E2) = apply_fluxes(dt2,rho2,u2,E2,ur2,Fjr2);
+    std::tie(new_m2,new_rho2,new_u2,new_E2) = apply_fluxes(dt2,tk,rho2,u2,E2,ur2_n,Fjr2_n,ur2_app,Fjr2_app);
+    //std::tie(new_m2,new_rho2,new_u2,new_E2) = apply_fluxes(dt2,tk,rho2,u2,E2,ur2_n,Fjr2_n,ur2_app,Fjr2_app,ur2_n,Fjr2_n);
+  
+    sum_dt = dt2;
+
+    while(sum_dt < dt1){
+
+      const DiscreteScalarFunction& rho_d = new_rho2->get<DiscreteScalarFunction>();
+      const DiscreteVectorFunction& u_d   = new_u2->get<DiscreteVectorFunction>();
+      const DiscreteScalarFunction& E_d   = new_E2->get<DiscreteScalarFunction>();
+      const DiscreteScalarFunction& eps   = E_d - 0.5 * dot(u_d,u_d);
+      const DiscreteScalarFunction& p_d   = (gamma-1) * rho_d * eps;
+      const DiscreteScalarFunction& c_d   = sqrt(gamma * p_d / rho_d);
+
+      new_p2 = std::make_shared<DiscreteFunctionVariant>(p_d);
+      new_c2 = std::make_shared<DiscreteFunctionVariant>(c_d);
+
+      //double dt2 = 0.4 * acoustic_dt(new_c2);
+      dt2 = dt1/2;
+      //dt2 = dt1;
+      //std::cout << "boucle" << std::endl;
+      if(sum_dt + dt2 > dt1){
+	    dt2 = dt1 - sum_dt;
+      }
+
+      tk = (sum_dt + dt2*0.5)/dt1; 
+      //tk = 0.5*3;   
+
+      //const MeshType& new_mesh2 = dynamic_cast<const MeshType&>(*new_m2);
+//
+//
+      //auto ur_inter            = _InterpolateUr(new_mesh2,ur2_n,ur2_app,tk);
+      //auto Fjr_inter      = _InterpolateFjr(new_mesh2,Fjr2_n,Fjr2_app,tk);
+      //auto [ur2_k,Fjr2_k] = compute_fluxes(solver_type,new_rho2,new_c2,new_u2,new_p2,bc_descriptor_list2,CR_ur,CR_Fjr,map2);
+      //auto new_ur2                 = _computeOrder2Ur(new_mesh2,ur_inter,ur2_k);
+      //auto new_Fjr2                = _computeOrder2Fjr(new_mesh2,Fjr_inter,Fjr2_k);
+
+      //std::tie(new_m2,new_rho2,new_u2,new_E2) = apply_fluxes(dt2,new_rho2,new_u2,new_E2,new_ur2,new_Fjr2);
+      std::tie(new_m2,new_rho2,new_u2,new_E2) = apply_fluxes(dt2,tk,new_rho2,new_u2,new_E2,ur2_n,Fjr2_n,ur2_app,Fjr2_app);
+      sum_dt += dt2;
+    }
+    std::tie(new_m2,new_rho2,new_u2,new_E2) = mesh_correction(new_m1,new_m2,new_rho2,new_u2,new_E2,map1,map2);
+    return {new_m1,new_rho1,new_u1,new_E1,new_m2,new_rho2,new_u2,new_E2};
+  }
+
+  //************************************************
+  //************************************************ Fin
+  //************************************************
+
+  LocalDtAcousticSolver()                        = default;
+  LocalDtAcousticSolver(LocalDtAcousticSolver&&) = default;
+  ~LocalDtAcousticSolver()                       = default;
+};
+
+template <size_t Dimension>
+void
+LocalDtAcousticSolverHandler::LocalDtAcousticSolver<Dimension>::_applyPressureBC(const BoundaryConditionList& bc_list,
+                                                                                 const MeshType& mesh,
+                                                                                 NodeValue<Rd>& br) const
+{
+  for (const auto& boundary_condition : bc_list) {
+    std::visit(
+      [&](auto&& bc) {
+        using T = std::decay_t<decltype(bc)>;
+        if constexpr (std::is_same_v<PressureBoundaryCondition, T>) {
+          MeshData<Dimension>& mesh_data = MeshDataManager::instance().getMeshData(mesh);
+          if constexpr (Dimension == 1) {
+            const NodeValuePerCell<const Rd> Cjr = mesh_data.Cjr();
+
+            const auto& node_to_cell_matrix               = mesh.connectivity().nodeToCellMatrix();
+            const auto& node_local_numbers_in_their_cells = mesh.connectivity().nodeLocalNumbersInTheirCells();
+
+            const auto& node_list  = bc.nodeList();
+            const auto& value_list = bc.valueList();
+            parallel_for(
+              node_list.size(), PUGS_LAMBDA(size_t i_node) {
+                const NodeId node_id       = node_list[i_node];
+                const auto& node_cell_list = node_to_cell_matrix[node_id];
+                Assert(node_cell_list.size() == 1);
+
+                CellId node_cell_id              = node_cell_list[0];
+                size_t node_local_number_in_cell = node_local_numbers_in_their_cells(node_id, 0);
+
+                br[node_id] -= value_list[i_node] * Cjr(node_cell_id, node_local_number_in_cell);
+              });
+          } else {
+            const NodeValuePerFace<const Rd> Nlr = mesh_data.Nlr();
+
+            const auto& face_to_cell_matrix               = mesh.connectivity().faceToCellMatrix();
+            const auto& face_to_node_matrix               = mesh.connectivity().faceToNodeMatrix();
+            const auto& face_local_numbers_in_their_cells = mesh.connectivity().faceLocalNumbersInTheirCells();
+            const auto& face_cell_is_reversed             = mesh.connectivity().cellFaceIsReversed();
+
+            const auto& face_list  = bc.faceList();
+            const auto& value_list = bc.valueList();
+            for (size_t i_face = 0; i_face < face_list.size(); ++i_face) {
+              const FaceId face_id       = face_list[i_face];
+              const auto& face_cell_list = face_to_cell_matrix[face_id];
+              Assert(face_cell_list.size() == 1);
+
+              CellId face_cell_id              = face_cell_list[0];
+              size_t face_local_number_in_cell = face_local_numbers_in_their_cells(face_id, 0);
+
+              const double sign = face_cell_is_reversed(face_cell_id, face_local_number_in_cell) ? -1 : 1;
+
+              const auto& face_nodes = face_to_node_matrix[face_id];
+
+              for (size_t i_node = 0; i_node < face_nodes.size(); ++i_node) {
+                NodeId node_id = face_nodes[i_node];
+                br[node_id] -= sign * value_list[i_face] * Nlr(face_id, i_node);
+              }
+            }
+          }
+        }
+      },
+      boundary_condition);
+  }
+}
+
+template <size_t Dimension>
+void
+LocalDtAcousticSolverHandler::LocalDtAcousticSolver<Dimension>::_applySymmetryBC(const BoundaryConditionList& bc_list,
+                                                                                 NodeValue<Rdxd>& Ar,
+                                                                                 NodeValue<Rd>& br) const
+{
+  for (const auto& boundary_condition : bc_list) {
+    std::visit(
+      [&](auto&& bc) {
+        using T = std::decay_t<decltype(bc)>;
+        if constexpr (std::is_same_v<SymmetryBoundaryCondition, T>) {
+          const Rd& n = bc.outgoingNormal();
+
+          const Rdxd I   = identity;
+          const Rdxd nxn = tensorProduct(n, n);
+          const Rdxd P   = I - nxn;
+
+          const Array<const NodeId>& node_list = bc.nodeList();
+          parallel_for(
+            bc.numberOfNodes(), PUGS_LAMBDA(int r_number) {
+              const NodeId r = node_list[r_number];
+
+              Ar[r] = P * Ar[r] * P + nxn;
+              br[r] = P * br[r];
+            });
+        }
+      },
+      boundary_condition);
+  }
+}
+
+template <size_t Dimension>
+void
+LocalDtAcousticSolverHandler::LocalDtAcousticSolver<Dimension>::_applyVelocityBC(const BoundaryConditionList& bc_list,
+                                                                                 NodeValue<Rdxd>& Ar,
+                                                                                 NodeValue<Rd>& br) const
+{
+  for (const auto& boundary_condition : bc_list) {
+    std::visit(
+      [&](auto&& bc) {
+        using T = std::decay_t<decltype(bc)>;
+        if constexpr (std::is_same_v<VelocityBoundaryCondition, T>) {
+          const auto& node_list  = bc.nodeList();
+          const auto& value_list = bc.valueList();
+
+          parallel_for(
+            node_list.size(), PUGS_LAMBDA(size_t i_node) {
+              NodeId node_id    = node_list[i_node];
+              const auto& value = value_list[i_node];
+
+              Ar[node_id] = identity;
+              br[node_id] = value;
+            });
+        } else if constexpr (std::is_same_v<FixedBoundaryCondition, T>) {
+          const auto& node_list = bc.nodeList();
+          parallel_for(
+            node_list.size(), PUGS_LAMBDA(size_t i_node) {
+              NodeId node_id = node_list[i_node];
+
+              Ar[node_id] = identity;
               br[node_id] = zero;
             });
         }
@@ -1016,6 +1752,29 @@ void LocalDtAcousticSolverHandler::LocalDtAcousticSolver<Dimension>::_applyCoupl
   }
 }
 
+template <size_t Dimension>
+void LocalDtAcousticSolverHandler::LocalDtAcousticSolver<Dimension>::_applyCouplingBC2(NodeValue<Rdxd>& Ar1,
+										      NodeValue<Rdxd>& Ar2,
+										      NodeValue<Rd>& ur1,
+										      NodeValue<Rd>& ur2,
+										      const std::vector<NodeId>& map1,
+										      const std::vector<NodeId>& map2) const
+{
+  size_t n = map1.size();
+  Rd lambda;
+
+  for(size_t i = 0; i<n; i++){
+
+    NodeId node_id1 = map1[i];
+    NodeId node_id2 = map2[i];
+    
+    lambda = inverse(inverse(Ar2[node_id2]) + inverse(Ar1[node_id1]))*(ur2[node_id2]-ur1[node_id1]); 
+
+    ur1[node_id1] = ur1[node_id1] + inverse(Ar1[node_id1])*lambda;
+    ur2[node_id2] = ur2[node_id2] - inverse(Ar2[node_id2])*lambda; 
+  }
+}
+
 template <size_t Dimension>
 void LocalDtAcousticSolverHandler::LocalDtAcousticSolver<Dimension>::_applyCouplingBC(const MeshType& mesh,
  										      NodeValue<Rd>& ur,
diff --git a/src/scheme/LocalDtAcousticSolver.hpp b/src/scheme/LocalDtAcousticSolver.hpp
index 3f84ae440822e22c8be14ddcb6dc60d1858e3723..c8d23daedde06a12c26a4a1086d26060c4cd15bb 100644
--- a/src/scheme/LocalDtAcousticSolver.hpp
+++ b/src/scheme/LocalDtAcousticSolver.hpp
@@ -33,6 +33,29 @@ class LocalDtAcousticSolverHandler
       const std::shared_ptr<const DiscreteFunctionVariant>& p,
       const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>& bc_descriptor_list) const = 0;
 
+    virtual std::tuple<std::shared_ptr<const IMesh>,
+                       std::shared_ptr<const DiscreteFunctionVariant>,
+                       std::shared_ptr<const DiscreteFunctionVariant>,
+                       std::shared_ptr<const DiscreteFunctionVariant>,
+	                     std::shared_ptr<const IMesh>,
+	                     std::shared_ptr<const DiscreteFunctionVariant>,
+	                     std::shared_ptr<const DiscreteFunctionVariant>,
+	                     std::shared_ptr<const DiscreteFunctionVariant>>
+    applyOrder2(const SolverType& solver_type,
+        const double& dt1,                                     
+        const std::shared_ptr<const DiscreteFunctionVariant>& rho1,
+	      const std::shared_ptr<const DiscreteFunctionVariant>& rho2,
+        const std::shared_ptr<const DiscreteFunctionVariant>& u1,
+	      const std::shared_ptr<const DiscreteFunctionVariant>& u2,
+        const std::shared_ptr<const DiscreteFunctionVariant>& E1,
+	      const std::shared_ptr<const DiscreteFunctionVariant>& E2,
+        const std::shared_ptr<const DiscreteFunctionVariant>& c1,
+	      const std::shared_ptr<const DiscreteFunctionVariant>& c2,
+        const std::shared_ptr<const DiscreteFunctionVariant>& p1,
+	      const std::shared_ptr<const DiscreteFunctionVariant>& p2,
+        const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>& bc_descriptor_list1,
+	      const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>& bc_descriptor_list2) const = 0;
+
     virtual std::tuple<std::shared_ptr<const IMesh>,
                        std::shared_ptr<const DiscreteFunctionVariant>,
                        std::shared_ptr<const DiscreteFunctionVariant>,
@@ -68,6 +91,29 @@ class LocalDtAcousticSolverHandler
 	  const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>& bc_descriptor_list1,
 	  const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>& bc_descriptor_list2) const = 0;
 
+    virtual std::tuple<std::shared_ptr<const IMesh>,
+             std::shared_ptr<const DiscreteFunctionVariant>,
+             std::shared_ptr<const DiscreteFunctionVariant>,
+             std::shared_ptr<const DiscreteFunctionVariant>,
+	           std::shared_ptr<const IMesh>,
+	           std::shared_ptr<const DiscreteFunctionVariant>,
+	           std::shared_ptr<const DiscreteFunctionVariant>,
+	           std::shared_ptr<const DiscreteFunctionVariant>>
+  apply(const SolverType& solver_type,
+        const double& dt1,                                     
+        const std::shared_ptr<const DiscreteFunctionVariant>& rho1,
+	      const std::shared_ptr<const DiscreteFunctionVariant>& rho2,
+        const std::shared_ptr<const DiscreteFunctionVariant>& u1,
+	      const std::shared_ptr<const DiscreteFunctionVariant>& u2,
+        const std::shared_ptr<const DiscreteFunctionVariant>& E1,
+	      const std::shared_ptr<const DiscreteFunctionVariant>& E2,
+        const std::shared_ptr<const DiscreteFunctionVariant>& c1,
+	      const std::shared_ptr<const DiscreteFunctionVariant>& c2,
+        const std::shared_ptr<const DiscreteFunctionVariant>& p1,
+	      const std::shared_ptr<const DiscreteFunctionVariant>& p2,
+        const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>& bc_descriptor_list1,
+	      const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>& bc_descriptor_list2) const = 0;
+
     ILocalDtAcousticSolver()                         = default;
     ILocalDtAcousticSolver(ILocalDtAcousticSolver&&) = default;
     ILocalDtAcousticSolver& operator=(ILocalDtAcousticSolver&&) = default;