diff --git a/src/scheme/LocalDtHyperelasticSolver.cpp b/src/scheme/LocalDtHyperelasticSolver.cpp
index 62ff4426c3d142e462067d204f3722eddc32414f..69f3ed208be443d8763475af58284187ad6d048c 100644
--- a/src/scheme/LocalDtHyperelasticSolver.cpp
+++ b/src/scheme/LocalDtHyperelasticSolver.cpp
@@ -836,26 +836,21 @@ class LocalDtHyperelasticSolverHandler::LocalDtHyperelasticSolver final : public
 	dt2 = dt1 - sum_dt;
       }
 
-      if(i == 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 DiscreteTensorFunction& sigma_d = -(gamma - 1) * rho_d * eps * I;
-	// const DiscreteScalarFunction& aL_d    = sqrt(gamma*(gamma-1)*eps);
-	// const DiscreteScalarFunction& aT_d    = 0 * aL_d;
-	const DiscreteTensorFunction& CG_d    = new_CG2->get<DiscreteTensorFunction>();
-	const DiscreteScalarFunction& p = fluid*(gamma - 1)*rho_d*eps;
-	const DiscreteTensorFunction& sigma_d = mu/sqrt(det(CG_d))* (CG_d - I) + lambda/sqrt(det(CG_d)) * log(sqrt(det(CG_d))) * I - p * I;
-	const DiscreteScalarFunction& aL_d = sqrt(lambda+2*mu/rho_d + gamma * p / rho_d);
-	const DiscreteScalarFunction& aT_d = sqrt(mu/rho_d);
-
-	const std::shared_ptr<const DiscreteFunctionVariant>& new_sigma2 = std::make_shared<const DiscreteFunctionVariant>(sigma_d);
-	const std::shared_ptr<const DiscreteFunctionVariant>& new_aL2 = std::make_shared<const DiscreteFunctionVariant>(aL_d);
-	const std::shared_ptr<const DiscreteFunctionVariant>& new_aT2 = std::make_shared<const DiscreteFunctionVariant>(aT_d);
-
-	auto [ur2,Fjr2] = compute_fluxes(solver_type,new_rho2,new_aL2,new_aT2,new_u2,new_sigma2,bc_descriptor_list2,CR_ur,CR_Fjr,map2);
-      }
+      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 DiscreteTensorFunction& CG_d    = new_CG2->get<DiscreteTensorFunction>();
+      const DiscreteScalarFunction& p = fluid*(gamma - 1)*rho_d*eps;
+      const DiscreteTensorFunction& sigma_d = mu/sqrt(det(CG_d))* (CG_d - I) + lambda/sqrt(det(CG_d)) * log(sqrt(det(CG_d))) * I - p * I;
+      const DiscreteScalarFunction& aL_d = sqrt(lambda+2*mu/rho_d + gamma * p / rho_d);
+      const DiscreteScalarFunction& aT_d = sqrt(mu/rho_d);
+
+      const std::shared_ptr<const DiscreteFunctionVariant>& new_sigma2 = std::make_shared<const DiscreteFunctionVariant>(sigma_d);
+      const std::shared_ptr<const DiscreteFunctionVariant>& new_aL2 = std::make_shared<const DiscreteFunctionVariant>(aL_d);
+      const std::shared_ptr<const DiscreteFunctionVariant>& new_aT2 = std::make_shared<const DiscreteFunctionVariant>(aT_d);
+
+      auto [ur2,Fjr2] = compute_fluxes(solver_type,new_rho2,new_aL2,new_aT2,new_u2,new_sigma2,bc_descriptor_list2,CR_ur,CR_Fjr,map2);
 	
       std::tie(new_m2,new_rho2,new_u2,new_E2,new_CG2) = apply_fluxes(dt2,new_rho2,new_u2,new_E2,new_CG2,ur2,Fjr2);