diff --git a/src/scheme/HyperplasticSolver.cpp b/src/scheme/HyperplasticSolver.cpp
index a5a7b4a618c38bead4b89fdd95bfb4a181b57e5e..d0c45c2239a34ce4ba8751e4f40cec69b9a62919 100644
--- a/src/scheme/HyperplasticSolver.cpp
+++ b/src/scheme/HyperplasticSolver.cpp
@@ -170,7 +170,11 @@ class HyperplasticSolverHandler::HyperplasticSolver final : public HyperplasticS
     double chi          = 0;
     const double normS2 = frobeniusNorm(Sd) * frobeniusNorm(Sd);
     double limit2       = 2. / 3 * yield * yield;
-    chi                 = std::max(0., std::sqrt(normS2) - std::sqrt(limit2));
+    const double power  = 0.5;
+    if ((normS2 - limit2) > 0) {
+      chi = std::pow((normS2 - limit2), power);
+    }
+    //    chi                 = std::max(0., std::sqrt(normS2) - std::sqrt(limit2));
     return chi;
   }
 
@@ -617,7 +621,7 @@ class HyperplasticSolverHandler::HyperplasticSolver final : public HyperplasticS
         new_Fe[j] += dt_over_Vj * elastic_fluxes;
         const double fenorm0 = frobeniusNorm(new_Fe[j]);
         chi[j]               = _compute_chi(sigma[j], yield[j]);
-        const R3x3 M         = -dt_over_Vj * chi[j] * zeta[j] * sigmas;
+        const R3x3 M         = -dt * chi[j] * zeta[j] * sigmas;
 
         new_Fe[j]            = EigenvalueSolver{}.computeExpForTinyMatrix(M) * new_Fe[j];
         const double fenorm1 = frobeniusNorm(new_Fe[j]);
@@ -704,14 +708,14 @@ class HyperplasticSolverHandler::HyperplasticSolver final : public HyperplasticS
         new_E[j] += dt_over_Mj * energy_fluxes;
         new_Fe[j] += dt_over_Vj * elastic_fluxes;
         chi[j]     = _compute_chi(sigma[j], yield[j]);
-        int sub_it = static_cast<int>(std::ceil(dt_over_Vj * chi[j] * zeta[j] * frobeniusNorm(sigmas)));
+        int sub_it = static_cast<int>(std::ceil(dt * chi[j] * zeta[j] * frobeniusNorm(sigmas)));
         if (sub_it > 1) {
-          std::cout << sub_it << " dt_over_Vj " << dt_over_Vj << " chi[j] " << chi[j] << " zeta[j] " << zeta[j]
+          std::cout << sub_it << " dt " << dt << " chi[j] " << chi[j] << " zeta[j] " << zeta[j]
                     << " frobeniusNorm(sigmas) " << frobeniusNorm(sigmas) << "\n";
         }
         for (int i = 0; i < sub_it; i++) {
-          const R3x3 M = Vj[j] * I + dt / static_cast<double>(sub_it) * chi[j] * zeta[j] * sigmas;
-          new_Fe[j]    = Vj[j] * inverse(M) * new_Fe[j];
+          const R3x3 M = I + dt / static_cast<double>(sub_it) * chi[j] * zeta[j] * sigmas;
+          new_Fe[j]    = inverse(M) * new_Fe[j];
         }
       });
     std::cout << "sum_chi " << sum(chi) << "\n";
@@ -795,7 +799,7 @@ class HyperplasticSolverHandler::HyperplasticSolver final : public HyperplasticS
         // new_Be[j] += dt_over_Vj * elastic_fluxes;
         // const double fenorm0 = frobeniusNorm(new_Be[j]);
         chi[j]           = _compute_chi(sigma[j], yield[j]);
-        const R3x3 M     = dt_over_Vj * (gradv3d - chi[j] * zeta[j] * sigmas);
+        const R3x3 M     = dt_over_Vj * gradv3d - dt * chi[j] * zeta[j] * sigmas;
         const R3x3 expM  = EigenvalueSolver{}.computeExpForTinyMatrix(M);
         const R3x3 expMT = EigenvalueSolver{}.computeExpForTinyMatrix(transpose(M));
         new_Be[j]        = expM * Be[j] * expMT;
@@ -906,32 +910,28 @@ class HyperplasticSolverHandler::HyperplasticSolver final : public HyperplasticS
     case RelaxationType::Exponential: {
       parallel_for(
         mesh->numberOfCells(), PUGS_LAMBDA(CellId j) {
-          const double dt_over_Vj = dt / Vj[j];
-          R3x3 sigmasn            = sigman[j] - 1. / 3 * trace(sigman[j]) * I;
-          R3x3 sigmasnp1          = sigmanp1[j] - 1. / 3 * trace(sigmanp1[j]) * I;
-          chi[j]                  = _compute_chi(0.5 * (sigmasn + sigmasnp1), yield[j]);
-          const R3x3 M            = -dt_over_Vj * chi[j] * zeta[j] * 0.5 * (sigmasn + sigmasnp1);
-          new_Fe[j]               = EigenvalueSolver{}.computeExpForTinyMatrix(M) * new_Fe[j];
+          R3x3 sigmasn   = sigman[j] - 1. / 3 * trace(sigman[j]) * I;
+          R3x3 sigmasnp1 = sigmanp1[j] - 1. / 3 * trace(sigmanp1[j]) * I;
+          chi[j]         = _compute_chi(0.5 * (sigmasn + sigmasnp1), yield[j]);
+          const R3x3 M   = -dt * chi[j] * zeta[j] * 0.5 * (sigmasn + sigmasnp1);
+          new_Fe[j]      = EigenvalueSolver{}.computeExpForTinyMatrix(M) * new_Fe[j];
         });
       break;
     }
     case RelaxationType::Implicit: {
       parallel_for(
         mesh->numberOfCells(), PUGS_LAMBDA(CellId j) {
-          const double dt_over_Vj = dt / Vj[j];
-          R3x3 sigmasn            = sigman[j] - 1. / 3 * trace(sigman[j]) * I;
-          R3x3 sigmasnp1          = sigmanp1[j] - 1. / 3 * trace(sigmanp1[j]) * I;
-          chi[j]                  = _compute_chi(0.5 * (sigmasn + sigmasnp1), yield[j]);
-          int sub_it =
-            static_cast<int>(std::ceil(dt_over_Vj * chi[j] * zeta[j] * frobeniusNorm(0.5 * (sigmasn + sigmasnp1))));
+          R3x3 sigmasn   = sigman[j] - 1. / 3 * trace(sigman[j]) * I;
+          R3x3 sigmasnp1 = sigmanp1[j] - 1. / 3 * trace(sigmanp1[j]) * I;
+          chi[j]         = _compute_chi(0.5 * (sigmasn + sigmasnp1), yield[j]);
+          int sub_it = static_cast<int>(std::ceil(dt * chi[j] * zeta[j] * frobeniusNorm(0.5 * (sigmasn + sigmasnp1))));
           if (sub_it > 1) {
-            std::cout << sub_it << " dt_over_Vj " << dt_over_Vj << " chi[j] " << chi[j] << " zeta[j] " << zeta[j]
+            std::cout << sub_it << " dt " << dt << " chi[j] " << chi[j] << " zeta[j] " << zeta[j]
                       << " frobeniusNorm(sigmas) " << frobeniusNorm(0.5 * (sigmasn + sigmasnp1)) << "\n";
           }
           for (int i = 0; i < sub_it; i++) {
-            const R3x3 M =
-              Vj[j] * I + 0.5 * dt / static_cast<double>(sub_it) * chi[j] * zeta[j] * (sigmasn + sigmasnp1);
-            new_Fe[j] = Vj[j] * inverse(M) * new_Fe[j];
+            const R3x3 M = I + 0.5 * dt / static_cast<double>(sub_it) * chi[j] * zeta[j] * (sigmasn + sigmasnp1);
+            new_Fe[j]    = inverse(M) * new_Fe[j];
           }
         });
       break;
diff --git a/src/scheme/HyperplasticSolverO2.cpp b/src/scheme/HyperplasticSolverO2.cpp
index 1d6174dd139213a263755dd1fb311484a74b4c86..41543fca0c6fe4f683aef5850af8b2d4802debfa 100644
--- a/src/scheme/HyperplasticSolverO2.cpp
+++ b/src/scheme/HyperplasticSolverO2.cpp
@@ -793,7 +793,7 @@ class HyperplasticSolverO2Handler::HyperplasticSolverO2 final
         new_Fe[j] += dt_over_Vj * elastic_fluxes;
         const double fenorm0 = frobeniusNorm(new_Fe[j]);
         chi[j]               = _compute_chi(sigma[j], yield[j]);
-        const R3x3 M         = -dt_over_Vj * chi[j] * zeta[j] * sigmas;
+        const R3x3 M         = -dt * chi[j] * zeta[j] * sigmas;
 
         new_Fe[j]            = EigenvalueSolver{}.computeExpForTinyMatrix(M) * new_Fe[j];
         const double fenorm1 = frobeniusNorm(new_Fe[j]);
@@ -880,14 +880,14 @@ class HyperplasticSolverO2Handler::HyperplasticSolverO2 final
         new_E[j] += dt_over_Mj * energy_fluxes;
         new_Fe[j] += dt_over_Vj * elastic_fluxes;
         chi[j]     = _compute_chi(sigma[j], yield[j]);
-        int sub_it = static_cast<int>(std::ceil(dt_over_Vj * chi[j] * zeta[j] * frobeniusNorm(sigmas)));
+        int sub_it = static_cast<int>(std::ceil(dt * chi[j] * zeta[j] * frobeniusNorm(sigmas)));
         if (sub_it > 1) {
-          std::cout << sub_it << " dt_over_Vj " << dt_over_Vj << " chi[j] " << chi[j] << " zeta[j] " << zeta[j]
+          std::cout << sub_it << " dt " << dt << " chi[j] " << chi[j] << " zeta[j] " << zeta[j]
                     << " frobeniusNorm(sigmas) " << frobeniusNorm(sigmas) << "\n";
         }
         for (int i = 0; i < sub_it; i++) {
-          const R3x3 M = Vj[j] * I + dt / static_cast<double>(sub_it) * chi[j] * zeta[j] * sigmas;
-          new_Fe[j]    = Vj[j] * inverse(M) * new_Fe[j];
+          const R3x3 M = I + dt / static_cast<double>(sub_it) * chi[j] * zeta[j] * sigmas;
+          new_Fe[j]    = inverse(M) * new_Fe[j];
         }
       });
     std::cout << "sum_chi " << sum(chi) << "\n";
@@ -971,7 +971,7 @@ class HyperplasticSolverO2Handler::HyperplasticSolverO2 final
         // new_Be[j] += dt_over_Vj * elastic_fluxes;
         // const double fenorm0 = frobeniusNorm(new_Be[j]);
         chi[j]           = _compute_chi(sigma[j], yield[j]);
-        const R3x3 M     = dt_over_Vj * (gradv3d - chi[j] * zeta[j] * sigmas);
+        const R3x3 M     = dt_over_Vj * gradv3d - dt * chi[j] * zeta[j] * sigmas;
         const R3x3 expM  = EigenvalueSolver{}.computeExpForTinyMatrix(M);
         const R3x3 expMT = EigenvalueSolver{}.computeExpForTinyMatrix(transpose(M));
         new_Be[j]        = expM * Be[j] * expMT;
@@ -1082,32 +1082,28 @@ class HyperplasticSolverO2Handler::HyperplasticSolverO2 final
     case RelaxationType::Exponential: {
       parallel_for(
         mesh->numberOfCells(), PUGS_LAMBDA(CellId j) {
-          const double dt_over_Vj = dt / Vj[j];
-          R3x3 sigmasn            = sigman[j] - 1. / 3 * trace(sigman[j]) * I;
-          R3x3 sigmasnp1          = sigmanp1[j] - 1. / 3 * trace(sigmanp1[j]) * I;
-          chi[j]                  = _compute_chi(0.5 * (sigmasn + sigmasnp1), yield[j]);
-          const R3x3 M            = -dt_over_Vj * chi[j] * zeta[j] * 0.5 * (sigmasn + sigmasnp1);
-          new_Fe[j]               = EigenvalueSolver{}.computeExpForTinyMatrix(M) * new_Fe[j];
+          R3x3 sigmasn   = sigman[j] - 1. / 3 * trace(sigman[j]) * I;
+          R3x3 sigmasnp1 = sigmanp1[j] - 1. / 3 * trace(sigmanp1[j]) * I;
+          chi[j]         = _compute_chi(0.5 * (sigmasn + sigmasnp1), yield[j]);
+          const R3x3 M   = -dt * chi[j] * zeta[j] * 0.5 * (sigmasn + sigmasnp1);
+          new_Fe[j]      = EigenvalueSolver{}.computeExpForTinyMatrix(M) * new_Fe[j];
         });
       break;
     }
     case RelaxationType::Implicit: {
       parallel_for(
         mesh->numberOfCells(), PUGS_LAMBDA(CellId j) {
-          const double dt_over_Vj = dt / Vj[j];
-          R3x3 sigmasn            = sigman[j] - 1. / 3 * trace(sigman[j]) * I;
-          R3x3 sigmasnp1          = sigmanp1[j] - 1. / 3 * trace(sigmanp1[j]) * I;
-          chi[j]                  = _compute_chi(0.5 * (sigmasn + sigmasnp1), yield[j]);
-          int sub_it =
-            static_cast<int>(std::ceil(dt_over_Vj * chi[j] * zeta[j] * frobeniusNorm(0.5 * (sigmasn + sigmasnp1))));
+          R3x3 sigmasn   = sigman[j] - 1. / 3 * trace(sigman[j]) * I;
+          R3x3 sigmasnp1 = sigmanp1[j] - 1. / 3 * trace(sigmanp1[j]) * I;
+          chi[j]         = _compute_chi(0.5 * (sigmasn + sigmasnp1), yield[j]);
+          int sub_it = static_cast<int>(std::ceil(dt * chi[j] * zeta[j] * frobeniusNorm(0.5 * (sigmasn + sigmasnp1))));
           if (sub_it > 1) {
-            std::cout << sub_it << " dt_over_Vj " << dt_over_Vj << " chi[j] " << chi[j] << " zeta[j] " << zeta[j]
+            std::cout << sub_it << " dt " << dt << " chi[j] " << chi[j] << " zeta[j] " << zeta[j]
                       << " frobeniusNorm(sigmas) " << frobeniusNorm(0.5 * (sigmasn + sigmasnp1)) << "\n";
           }
           for (int i = 0; i < sub_it; i++) {
-            const R3x3 M =
-              Vj[j] * I + 0.5 * dt / static_cast<double>(sub_it) * chi[j] * zeta[j] * (sigmasn + sigmasnp1);
-            new_Fe[j] = Vj[j] * inverse(M) * new_Fe[j];
+            const R3x3 M = I + 0.5 * dt / static_cast<double>(sub_it) * chi[j] * zeta[j] * (sigmasn + sigmasnp1);
+            new_Fe[j]    = inverse(M) * new_Fe[j];
           }
         });
       break;