diff --git a/src/language/modules/SchemeModule.cpp b/src/language/modules/SchemeModule.cpp
index 9a367ef65a4770ac51b4b701b5fe39fe70848f59..b1d6068829fb5df8bfa681976a75605ce8a6f469 100644
--- a/src/language/modules/SchemeModule.cpp
+++ b/src/language/modules/SchemeModule.cpp
@@ -180,7 +180,7 @@ struct GlaceScheme
                 << " time=" << rang::fg::green << t << rang::style::reset << " dt=" << rang::fgB::blue << dt
                 << rang::style::reset << '\n';
 
-      acoustic_solver.computeNextStep(t, dt, unknowns);
+      m_mesh = acoustic_solver.computeNextStep(dt, unknowns);
 
       block_eos.updatePandCFromRhoE();
 
diff --git a/src/main.cpp b/src/main.cpp
index aa9e2e3ff6378c1f1e20fe27388ea6019e765b08..daf0d4edb97c9939c5b25217649452207a54816f 100644
--- a/src/main.cpp
+++ b/src/main.cpp
@@ -155,7 +155,7 @@ main(int argc, char* argv[])
         if (t + dt > tmax) {
           dt = tmax - t;
         }
-        acoustic_solver.computeNextStep(t, dt, unknowns);
+        mesh = acoustic_solver.computeNextStep(dt, unknowns);
 
         block_eos.updatePandCFromRhoE();
 
@@ -272,7 +272,7 @@ main(int argc, char* argv[])
         if (t + dt > tmax) {
           dt = tmax - t;
         }
-        acoustic_solver.computeNextStep(t, dt, unknowns);
+        mesh = acoustic_solver.computeNextStep(dt, unknowns);
 
         block_eos.updatePandCFromRhoE();
 
@@ -378,7 +378,7 @@ main(int argc, char* argv[])
         if (t + dt > tmax) {
           dt = tmax - t;
         }
-        acoustic_solver.computeNextStep(t, dt, unknowns);
+        mesh = acoustic_solver.computeNextStep(dt, unknowns);
         block_eos.updatePandCFromRhoE();
 
         t += dt;
diff --git a/src/mesh/Mesh.hpp b/src/mesh/Mesh.hpp
index 0fe378e7e4c51e8e3a4a20b077d0f9025a0e818a..d04807650cd08c0b96a7259cbe5756596d54aad1 100644
--- a/src/mesh/Mesh.hpp
+++ b/src/mesh/Mesh.hpp
@@ -20,7 +20,6 @@ class Mesh final : public IMesh
  private:
   const std::shared_ptr<const Connectivity> m_connectivity;
   NodeValue<const Rd> m_xr;
-  NodeValue<Rd> m_mutable_xr;
 
  public:
   PUGS_INLINE
@@ -72,16 +71,9 @@ class Mesh final : public IMesh
     return m_xr;
   }
 
-  PUGS_INLINE
-  NodeValue<Rd>
-  mutableXr() const
-  {
-    return m_mutable_xr;
-  }
-
   PUGS_INLINE
   Mesh(const std::shared_ptr<const Connectivity>& connectivity, NodeValue<Rd>& xr)
-    : m_connectivity{connectivity}, m_xr{xr}, m_mutable_xr{xr}
+    : m_connectivity{connectivity}, m_xr{xr}
   {
     ;
   }
diff --git a/src/scheme/AcousticSolver.hpp b/src/scheme/AcousticSolver.hpp
index 0a5ef970d4a65242def2fe0fbad937e5fbab02ad..c531154bd5db52a3f102dfc8e0dd455157e01351 100644
--- a/src/scheme/AcousticSolver.hpp
+++ b/src/scheme/AcousticSolver.hpp
@@ -34,7 +34,7 @@ class AcousticSolver
   using MeshDataType = MeshData<Dimension>;
   using UnknownsType = FiniteVolumesEulerUnknowns<MeshDataType>;
 
-  const std::shared_ptr<MeshType> m_mesh;
+  std::shared_ptr<const MeshType> m_mesh;
   const typename MeshType::Connectivity& m_connectivity;
   const std::vector<BoundaryConditionHandler>& m_boundary_condition_list;
 
@@ -229,7 +229,7 @@ class AcousticSolver
   CellValue<double> m_Vj_over_cj;
 
  public:
-  AcousticSolver(std::shared_ptr<MeshType>& p_mesh, const std::vector<BoundaryConditionHandler>& bc_list)
+  AcousticSolver(std::shared_ptr<const MeshType> p_mesh, const std::vector<BoundaryConditionHandler>& bc_list)
     : m_mesh(p_mesh),
       m_connectivity(m_mesh->connectivity()),
       m_boundary_condition_list(bc_list),
@@ -268,8 +268,8 @@ class AcousticSolver
     return min(m_Vj_over_cj);
   }
 
-  void
-  computeNextStep(double, double dt, UnknownsType& unknowns)
+  [[nodiscard]] std::shared_ptr<const MeshType>
+  computeNextStep(double dt, UnknownsType& unknowns)
   {
     CellValue<double>& rhoj = unknowns.rhoj();
     CellValue<Rd>& uj       = unknowns.uj();
@@ -311,15 +311,20 @@ class AcousticSolver
     parallel_for(
       m_mesh->numberOfCells(), PUGS_LAMBDA(CellId j) { ej[j] = Ej[j] - 0.5 * (uj[j], uj[j]); });
 
-    NodeValue<Rd> new_xr = m_mesh->mutableXr();
+    NodeValue<Rd> new_xr = copy(m_mesh->xr());
     parallel_for(
       m_mesh->numberOfNodes(), PUGS_LAMBDA(NodeId r) { new_xr[r] += dt * ur[r]; });
 
     MeshDataManager::instance().getMeshData(*m_mesh).updateAllData();
 
+    m_mesh                         = std::make_shared<MeshType>(m_mesh->shared_connectivity(), new_xr);
+    CellValue<const double> new_Vj = MeshDataManager::instance().getMeshData(*m_mesh).Vj();
+
     const CellValue<const double> mj = unknowns.mj();
     parallel_for(
-      m_mesh->numberOfCells(), PUGS_LAMBDA(CellId j) { rhoj[j] = mj[j] / Vj[j]; });
+      m_mesh->numberOfCells(), PUGS_LAMBDA(CellId j) { rhoj[j] = mj[j] / new_Vj[j]; });
+
+    return m_mesh;
   }
 };