diff --git a/src/language/modules/LocalFSIModule.cpp b/src/language/modules/LocalFSIModule.cpp
index eeffc83f5486b0dcfe09480e3c300ff91f268532..7597da325c343a0c8e3422c8f922688908c9ca13 100644
--- a/src/language/modules/LocalFSIModule.cpp
+++ b/src/language/modules/LocalFSIModule.cpp
@@ -45,6 +45,8 @@ LocalFSIModule::LocalFSIModule()
                                  const std::shared_ptr<const DiscreteFunctionVariant>& p,
                                  const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>&
                                    bc_descriptor_list,
+                                 const std::shared_ptr<const DiscreteFunctionVariant>& gamma,
+                                 const std::shared_ptr<const DiscreteFunctionVariant>& pinf,
                                  const double& dt) -> std::tuple<std::shared_ptr<const MeshVariant>,
                                                                  std::shared_ptr<const DiscreteFunctionVariant>,
                                                                  std::shared_ptr<const DiscreteFunctionVariant>,
@@ -52,7 +54,7 @@ LocalFSIModule::LocalFSIModule()
                                 return Order2AcousticSolverHandler{getCommonMesh({rho, u, E, c, p})}
                                   .solver()
                                   .apply(Order2AcousticSolverHandler::SolverType::Eucclhyd, dt, rho, u, E, c, p,
-                                         bc_descriptor_list);
+                                         bc_descriptor_list, gamma, pinf);
                               }
 
                               ));
@@ -65,6 +67,8 @@ LocalFSIModule::LocalFSIModule()
                                  const std::shared_ptr<const DiscreteFunctionVariant>& E,
                                  const std::shared_ptr<const DiscreteFunctionVariant>& c,
                                  const std::shared_ptr<const DiscreteFunctionVariant>& p,
+                                 const std::shared_ptr<const DiscreteFunctionVariant>& gamma,
+                                 const std::shared_ptr<const DiscreteFunctionVariant>& pinf,
                                  const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>&
                                    bc_descriptor_list,
                                  const double& dt) -> std::tuple<std::shared_ptr<const MeshVariant>,
@@ -74,7 +78,7 @@ LocalFSIModule::LocalFSIModule()
                                 return Order2AcousticSolverHandler{getCommonMesh({rho, u, E, c, p})}
                                   .solver()
                                   .apply(Order2AcousticSolverHandler::SolverType::Glace, dt, rho, u, E, c, p,
-                                         bc_descriptor_list);
+                                         bc_descriptor_list, gamma, pinf);
                               }
 
                               ));
diff --git a/src/scheme/Order2AcousticSolver.cpp b/src/scheme/Order2AcousticSolver.cpp
index 0c2c166ef77412a986d594f65c25373df544a1cf..c01a423cdca606474af757c357428cd02ebca264 100644
--- a/src/scheme/Order2AcousticSolver.cpp
+++ b/src/scheme/Order2AcousticSolver.cpp
@@ -20,6 +20,7 @@
 #include <scheme/FixedBoundaryConditionDescriptor.hpp>
 #include <scheme/IBoundaryConditionDescriptor.hpp>
 #include <scheme/IDiscreteFunctionDescriptor.hpp>
+#include <scheme/Order2Limiters.hpp>
 #include <scheme/PolynomialReconstruction.hpp>
 #include <scheme/PolynomialReconstructionDescriptor.hpp>
 #include <scheme/SymmetryBoundaryConditionDescriptor.hpp>
@@ -172,6 +173,39 @@ class Order2AcousticSolverHandler::Order2AcousticSolver final : public Order2Aco
     return Ar;
   }
 
+  NodeValue<Rd>
+  _computeBr(const Mesh<Dimension>& mesh,
+             const NodeValuePerCell<const Rdxd>& Ajr,
+             const DiscreteVectorFunction& u,
+             const DiscreteScalarFunction& p) const
+  {
+    MeshDataType& mesh_data = MeshDataManager::instance().getMeshData(mesh);
+
+    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();
+
+    NodeValue<Rd> b{mesh.connectivity()};
+
+    parallel_for(
+      mesh.numberOfNodes(), PUGS_LAMBDA(NodeId r) {
+        const auto& node_to_cell                   = node_to_cell_matrix[r];
+        const auto& node_local_number_in_its_cells = node_local_numbers_in_their_cells.itemArray(r);
+
+        Rd br = zero;
+        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];
+          br += Ajr(J, R) * u[J] + p[J] * Cjr(J, R);
+        }
+
+        b[r] = br;
+      });
+
+    return b;
+  }
+
   NodeValue<Rd>
   _computeBr(const Mesh<Dimension>& mesh,
                    const NodeValuePerCell<const Rdxd>& Ajr,
@@ -308,6 +342,61 @@ class Order2AcousticSolverHandler::Order2AcousticSolver final : public Order2Aco
     this->_applyVelocityBC(bc_list, Ar, br);
   }
 
+  CellValue<Rdxd>
+  _computeGradU(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 mesh_v = getCommonMesh({rho_v, c_v, u_v, p_v});
+    if (not mesh_v) {
+      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("hyperelastic solver expects P0 functions");
+    }
+
+    const MeshType& mesh              = *mesh_v->get<MeshType>();
+    const DiscreteScalarFunction& rho = rho_v->get<DiscreteScalarFunction>();
+    const DiscreteVectorFunction& u   = u_v->get<DiscreteVectorFunction>();
+    const DiscreteScalarFunction& c   = c_v->get<DiscreteScalarFunction>();
+    const DiscreteScalarFunction& p   = p_v->get<DiscreteScalarFunction>();
+
+    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);
+
+    synchronize(Ar);
+    synchronize(br);
+
+    const auto& cell_to_node_matrix = mesh.connectivity().cellToNodeMatrix();
+
+    const NodeValuePerCell<const Rd> Cjr = MeshDataManager::instance().getMeshData(mesh).Cjr();
+    const NodeValue<const Rd> ur         = this->_computeUr(mesh, Ar, br);
+
+    CellValue<Rdxd> grad_u{mesh.connectivity()};
+    parallel_for(
+      mesh.numberOfCells(), PUGS_LAMBDA(CellId j) {
+        auto cell_nodes = cell_to_node_matrix[j];
+
+        Rdxd gradv           = zero;
+        for (size_t R = 0; R < cell_nodes.size(); ++R) {
+          NodeId r = cell_nodes[R];
+          gradv += tensorProduct(ur[r], Cjr(j, R));
+        }
+        grad_u[j] = gradv;
+      });
+
+    return grad_u;
+  }
+
   NodeValue<const Rd>
   _computeUr(const MeshType& mesh, const NodeValue<Rdxd>& Ar, const NodeValue<Rd>& br) const
   {
@@ -349,228 +438,6 @@ class Order2AcousticSolverHandler::Order2AcousticSolver final : public Order2Aco
     return F;
   }
 
-  void
-  _scalar_limiter(const MeshType& mesh,
-                  const DiscreteFunctionP0<const double>& f,
-                  DiscreteFunctionDPk<Dimension, double>& DPk_fh) const
-  {
-    MeshData<MeshType>& mesh_data = MeshDataManager::instance().getMeshData(mesh);
-    StencilManager::BoundaryDescriptorList symmetry_boundary_descriptor_list;
-    StencilDescriptor stencil_descriptor{1, StencilDescriptor::ConnectionType::by_nodes};
-    auto stencil = StencilManager::instance().getCellToCellStencilArray(mesh.connectivity(), stencil_descriptor, symmetry_boundary_descriptor_list);
-
-    const auto xj = mesh_data.xj();
-
-    parallel_for(mesh.numberOfCells(), PUGS_LAMBDA(const CellId cell_id){
-      const double fj = f[cell_id];
-
-      double f_min = fj;
-      double f_max = fj;
-
-      const auto cell_stencil = stencil[cell_id];
-      for(size_t i_cell = 0; i_cell < cell_stencil.size(); ++i_cell){
-        f_min = std::min(f_min, f[cell_stencil[i_cell]]);
-        f_max = std::max(f_max, f[cell_stencil[i_cell]]);
-      }
-
-      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](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 coef2 = 1.;
-      if (std::abs(f_bar_min - fj) > eps) {
-        coef2 = (f_min - fj) / ((f_bar_min - fj));
-      }
-
-      const double lambda = std::max(0., std::min(1., std::min(coef1, coef2)));
-
-      auto coefficients = DPk_fh.coefficients(cell_id);
-
-      coefficients[0] = (1 - lambda) * f[cell_id] + lambda * coefficients[0];
-      for (size_t i = 1; i < coefficients.size(); ++i) {
-        coefficients[i] *= lambda;
-      }
-    });
-  }
-
-  void 
-  _vector_limiter(const MeshType& mesh,
-                  const DiscreteFunctionP0<const Rd>& f,
-                  DiscreteFunctionDPk<Dimension, Rd>& DPK_fh,
-                  const DiscreteFunctionDPk<Dimension, const double>& DPk_ph) const 
-  {
-    MeshData<MeshType>& mesh_data = MeshDataManager::instance().getMeshData(mesh);
-    StencilManager::BoundaryDescriptorList symmetry_boundary_descriptor_list;
-    StencilDescriptor stencil_descriptor{1, StencilDescriptor::ConnectionType::by_nodes};
-    auto stencil = StencilManager::instance().getCellToCellStencilArray(mesh.connectivity(), stencil_descriptor,
-                                                                        symmetry_boundary_descriptor_list);
-    
-    const auto xj = mesh_data.xj();
-
-    CellValue<Rd> n{mesh.connectivity()};
-    parallel_for(
-      mesh.numberOfCells(), PUGS_LAMBDA(const CellId cell_id) {
-        auto coefficients_p = DPk_ph.coefficients(cell_id);
-        Rd grad_p = zero;
-        for(size_t i = 1; i < coefficients_p.size(); ++i){
-          grad_p[i-1] = coefficients_p[i];
-        }
-        const double norm_grad_p = l2Norm(grad_p);
-        if(norm_grad_p == 0){
-          n[cell_id] = zero;
-        } 
-        else {
-          n[cell_id] = (1/norm_grad_p) * grad_p;
-        }
-
-        
-      });
-
-    const CellValue<Rd> t{mesh.connectivity()};
-    const CellValue<Rd> l{mesh.connectivity()};
-    parallel_for(
-      mesh.numberOfCells(), PUGS_LAMBDA(const CellId cell_id) {
-        const Rd nj = n[cell_id];
-        Rd a = zero;
-        if(l2Norm(nj) != 0){
-          if((nj[0] / l2Norm(nj) != 1) and (nj[0] / l2Norm(nj) != -1)){
-            a[0] = 1.;
-          }
-          else { 
-            a[1] = 1.;
-          }
-        }
-
-        Rd tj = a - dot(a,nj) * nj;
-        const double& norm_tj = l2Norm(tj);
-        if(norm_tj == 0){
-          tj = zero;
-        }
-        else {
-          tj = (1/norm_tj) * tj;
-        }
-        t[cell_id] = tj;
-
-        Rd lj = zero;
-        if(Dimension == 3){
-          lj[0] = nj[1]*tj[2] - nj[2]*tj[1];
-          lj[1] = nj[2]*tj[0] - nj[0]*tj[2];
-          lj[2] = nj[0]*tj[1] - nj[1]*tj[0];
-          const double& norm_lj = l2Norm(lj);
-          lj = (1/norm_lj) * lj;
-        }
-        l[cell_id] = lj;
-      });
-
-    parallel_for(
-      mesh.numberOfCells(), PUGS_LAMBDA(const CellId cell_id){
-        const double fn = dot(f[cell_id],n[cell_id]);
-        const double ft = dot(f[cell_id],t[cell_id]);
-        const double fl = dot(f[cell_id],l[cell_id]);
-
-        double fn_min = fn;
-        double fn_max = fn;
-        double ft_min = ft;
-        double ft_max = ft;
-        double fl_min = fl;
-        double fl_max = fl;
-
-        const auto cell_stencil = stencil[cell_id];
-        for(size_t i_cell = 0; i_cell < cell_stencil.size(); ++i_cell){
-          const double fn_k = dot(f[cell_stencil[i_cell]],n[cell_stencil[i_cell]]);
-          fn_min            = std::min(fn_min,fn_k);
-          fn_max            = std::max(fn_max,fn_k);
-
-          const double ft_k = dot(f[cell_stencil[i_cell]],t[cell_stencil[i_cell]]);
-          ft_min            = std::min(ft_min,ft_k);
-          ft_max            = std::max(ft_max,ft_k);
-
-          const double fl_k = dot(f[cell_stencil[i_cell]],l[cell_stencil[i_cell]]);
-          fl_min            = std::min(fl_min,fl_k);
-          fl_max            = std::max(fl_max,fl_k);
-        }
-
-        double fn_bar_min = fn;
-        double fn_bar_max = fn;
-        double ft_bar_min = ft;
-        double ft_bar_max = ft;
-        double fl_bar_min = fl;
-        double fl_bar_max = fl;
-
-        for(size_t i_cell = 0; i_cell < cell_stencil.size(); ++i_cell){
-          const CellId cell_k_id = cell_stencil[i_cell];
-          const Rd f_xk          = DPK_fh[cell_id](xj[cell_k_id]);
-
-          const double fn_xk = dot(f_xk,n[cell_k_id]);
-          fn_bar_min         = std::min(fn_bar_min,fn_xk);
-          fn_bar_max         = std::max(fn_bar_max,fn_xk);
-
-          const double ft_xk = dot(f_xk,t[cell_k_id]);
-          ft_bar_min         = std::min(ft_bar_min,ft_xk);
-          ft_bar_max         = std::max(ft_bar_max,ft_xk);
-
-          const double fl_xk = dot(f_xk,l[cell_k_id]);
-          fl_bar_min         = std::min(fl_bar_min,fl_xk);
-          fl_bar_max         = std::max(fl_bar_max,fl_xk);
-        }
-
-        const double eps = 1E-14;
-        double coef1_n   = 1;
-        if(std::abs(fn_bar_max - fn) > eps){
-          coef1_n = (fn_max - fn) / (fn_bar_max - fn);
-        } 
-        double coef2_n = 1;
-        if(std::abs(fn_bar_min - fn) > eps){
-          coef2_n = (fn_min - fn) / (fn_bar_min - fn);
-        }
-        
-        double coef1_t = 1;
-        if(std::abs(ft_bar_max - ft) > eps){
-          coef1_t = (ft_max - ft) / (ft_bar_max - ft);
-        }
-        double coef2_t = 1;
-        if(std::abs(ft_bar_min - ft) > eps){
-          coef2_t = (ft_min - ft) / (ft_bar_min - ft);
-        }
-
-        double coef1_l = 1;
-        if(std::abs(fl_bar_max - fl) > eps){
-          coef1_l = (fl_max - fl) / (fl_bar_max - fl);
-        }
-        double coef2_l = 1;
-        if(std::abs(fl_bar_min - fl) > eps){
-          coef2_l = (fl_min - fl) / (fl_bar_min - fl);
-        }
-
-        const double lambda_n = std::max(0., std::min(1., std::min(coef1_n,coef2_n)));
-        const double lambda_t = std::max(0., std::min(1., std::min(coef1_t,coef2_t)));
-        const double lambda_l = std::max(0., std::min(1., std::min(coef1_l,coef2_l)));
-
-        auto coefficients = DPK_fh.coefficients(cell_id);
-        coefficients[0] = (1 - lambda_n)*fn*n[cell_id] + lambda_n*dot(coefficients[0],n[cell_id])*n[cell_id] 
-                        + (1 - lambda_t)*ft*t[cell_id] + lambda_t*dot(coefficients[0],t[cell_id])*t[cell_id] 
-                        + (1 - lambda_l)*fl*l[cell_id] + lambda_l*dot(coefficients[0],l[cell_id])*l[cell_id];
-        for(size_t i = 1; i < coefficients.size(); ++i){
-          coefficients[i] = lambda_n*dot(coefficients[i],n[cell_id])*n[cell_id] 
-                          + lambda_t*dot(coefficients[i],t[cell_id])*t[cell_id] 
-                          + lambda_l*dot(coefficients[i],l[cell_id])*l[cell_id];
-        }
-    });
-  }
-
  public:
   std::tuple<const std::shared_ptr<const ItemValueVariant>, const std::shared_ptr<const SubItemValuePerItemVariant>>
   compute_fluxes(const SolverType& solver_type,
@@ -617,8 +484,9 @@ class Order2AcousticSolverHandler::Order2AcousticSolver final : public Order2Aco
     DiscreteFunctionDPk<Dimension, Rd> u_lim     = copy(DPk_uh);
     DiscreteFunctionDPk<Dimension, double> p_lim = copy(DPk_ph);
 
-    this->_scalar_limiter(mesh,p,p_lim);
-    this->_vector_limiter(mesh,u,u_lim,DPk_ph);
+    const CellValue<const Rdxd> grad_u = this->_computeGradU(solver_type, rho_v, c_v, u_v, p_v, bc_descriptor_list);
+    scalar_limiter(mesh,p,p_lim);
+    vector_limiter(mesh,u,u_lim,grad_u);
 
     NodeValuePerCell<const Rdxd> Ajr = this->_computeAjr(solver_type, mesh, rho * c);
 
@@ -737,7 +605,9 @@ class Order2AcousticSolverHandler::Order2AcousticSolver final : public Order2Aco
         const std::shared_ptr<const DiscreteFunctionVariant>& E,
         const std::shared_ptr<const DiscreteFunctionVariant>& c,
         const std::shared_ptr<const DiscreteFunctionVariant>& p,
-        const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>& bc_descriptor_list) const
+        const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>& bc_descriptor_list,
+        const std::shared_ptr<const DiscreteFunctionVariant>& gamma,
+        const std::shared_ptr<const DiscreteFunctionVariant>& pinf) const
   {
     std::shared_ptr m_app = getCommonMesh({rho, c, u, p});
     std::shared_ptr<const DiscreteFunctionVariant> rho_app = rho;
@@ -747,14 +617,17 @@ class Order2AcousticSolverHandler::Order2AcousticSolver final : public Order2Aco
     auto [ur, Fjr] = compute_fluxes(solver_type, rho, c, u, p, bc_descriptor_list);
     std::tie(m_app,rho_app,u_app,E_app) = apply_fluxes(dt/2, rho, u, E, ur, Fjr);
 
-    const double& gamma = 1.4;
-
-    const DiscreteScalarFunction& rho_d = rho_app->get<DiscreteScalarFunction>();
-    const DiscreteVectorFunction& u_d   = u_app->get<DiscreteVectorFunction>();
-    const DiscreteScalarFunction& E_d   = E_app->get<DiscreteScalarFunction>();
-    const DiscreteScalarFunction& eps_d = E_d - 0.5 * dot(u_d, u_d);
-    const DiscreteScalarFunction& p_d   = (gamma - 1) * rho_d * eps_d;
-    const DiscreteScalarFunction& c_d   = sqrt(gamma * p_d / rho_d);
+    const std::shared_ptr<const DiscreteFunctionVariant>& gamma_app = shallowCopy(m_app,gamma);
+    const std::shared_ptr<const DiscreteFunctionVariant>& pinf_app  = shallowCopy(m_app,pinf);
+
+    const DiscreteScalarFunction& rho_d   = rho_app->get<DiscreteScalarFunction>();
+    const DiscreteVectorFunction& u_d     = u_app->get<DiscreteVectorFunction>();
+    const DiscreteScalarFunction& E_d     = E_app->get<DiscreteScalarFunction>();
+    const DiscreteScalarFunction& gamma_d = gamma_app->get<DiscreteScalarFunction>();
+    const DiscreteScalarFunction& pinf_d  = pinf_app->get<DiscreteScalarFunction>();
+    const DiscreteScalarFunction& eps_d   = E_d - 0.5 * dot(u_d, u_d);
+    const DiscreteScalarFunction& p_d     = (gamma_d - 1) * rho_d * eps_d - gamma_d * pinf_d;
+    const DiscreteScalarFunction& c_d     = sqrt(gamma_d * (p_d + pinf_d)/ rho_d);
 
     const std::shared_ptr<const DiscreteFunctionVariant>& p_app =
       std::make_shared<const DiscreteFunctionVariant>(p_d);
diff --git a/src/scheme/Order2AcousticSolver.hpp b/src/scheme/Order2AcousticSolver.hpp
index fec4185dac437a53747484a5e5fb7cb31f719c88..73924b0854124c3448c95aa43826ac98c083a1e6 100644
--- a/src/scheme/Order2AcousticSolver.hpp
+++ b/src/scheme/Order2AcousticSolver.hpp
@@ -57,7 +57,9 @@ class Order2AcousticSolverHandler
           const std::shared_ptr<const DiscreteFunctionVariant>& E,
           const std::shared_ptr<const DiscreteFunctionVariant>& c,
           const std::shared_ptr<const DiscreteFunctionVariant>& p,
-          const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>& bc_descriptor_list) const = 0;
+          const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>& bc_descriptor_list,
+          const std::shared_ptr<const DiscreteFunctionVariant>& gamma,
+          const std::shared_ptr<const DiscreteFunctionVariant>& pinf) const = 0;
 
     IOrder2AcousticSolver()                  = default;
     IOrder2AcousticSolver(IOrder2AcousticSolver&&) = default;