diff --git a/src/scheme/ImplicitAcousticSolver.cpp b/src/scheme/ImplicitAcousticSolver.cpp
index 411c4112c0f1af52061ffb791094732cdd03f538..df37a7591ee3112debc34200fcac7c02729231da 100644
--- a/src/scheme/ImplicitAcousticSolver.cpp
+++ b/src/scheme/ImplicitAcousticSolver.cpp
@@ -1133,25 +1133,8 @@ class ImplicitAcousticSolverHandler::ImplicitAcousticSolver final
 
     CRSMatrix Hess_J_crs = Hess_J.getCRSMatrix();
 
-    // std::cout << "Hess J = " << Hess_J_crs << '\n';
-
     CRSMatrix gradient_f = Hess_J_crs - A;
 
-    // Vector<double> x(gradient_f.numberOfRows());
-    // Vector<double> y(gradient_f.numberOfRows());
-    // double max_err = 0;
-    // for (size_t i = 0; i < gradient_f.numberOfRows(); ++i) {
-    //   x.fill(0);
-    //   x[i] = 1;
-    //   for (size_t j = 0; j < gradient_f.numberOfRows(); ++j) {
-    //     y.fill(0);
-    //     y[j] = 1;
-
-    //     max_err = std::max(max_err, std::abs(dot(x, gradient_f * y) - dot(y, gradient_f * x)));
-    //   }
-    // }
-    // std::cout << "MAX NON SYM=" << rang::fgB::cyan << max_err << rang::fg::reset << '\n';
-
     HJ_A_t.pause();
 
     return gradient_f;
@@ -1628,6 +1611,47 @@ class ImplicitAcousticSolverHandler::ImplicitAcousticSolver final
     parallel_for(
       m_mesh.numberOfNodes(), PUGS_LAMBDA(NodeId r) { computed_ur[r] = m_inv_Ar[r] * b[r]; });
 
+    for (const auto& boundary_condition : m_boundary_condition_list) {
+      std::visit(
+        [&](auto&& bc) {
+          using T = std::decay_t<decltype(bc)>;
+
+          if constexpr ((Dimension > 1) and (std::is_same_v<AxisBoundaryCondition, T>)) {
+            const Rd& t    = bc.direction();
+            const Rdxd txt = tensorProduct(t, t);
+
+            const Array<const NodeId>& node_list = bc.nodeList();
+
+            for (size_t i_node = 0; i_node < node_list.size(); ++i_node) {
+              const NodeId node_id = node_list[i_node];   // on fixe le sommet r
+              computed_ur[node_id] = txt * computed_ur[node_id];
+            }
+          } else 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();
+
+            for (size_t i_node = 0; i_node < node_list.size(); ++i_node) {
+              const NodeId node_id = node_list[i_node];   // on fixe le sommet r
+              computed_ur[node_id] = P * computed_ur[node_id];
+            }
+          } else if constexpr (std::is_same_v<VelocityBoundaryCondition, T>) {
+            const Array<const NodeId>& node_list = bc.nodeList();
+            const Array<const Rd>& value_list    = bc.valueList();
+
+            for (size_t i_node = 0; i_node < node_list.size(); ++i_node) {
+              const NodeId node_id = node_list[i_node];
+              computed_ur[node_id] = value_list[i_node];
+            }
+          }
+        },
+        boundary_condition);
+    }
+
     return computed_ur;
   }