From d4676b10634b28051ac638f56fbdb4ee04ba5c7a Mon Sep 17 00:00:00 2001
From: labourasse <labourassee@gmail.com>
Date: Wed, 8 Jan 2025 14:24:59 +0100
Subject: [PATCH] fix a bug in the relaxation
---
src/scheme/HyperplasticSolver.cpp | 46 ++++++++++++++---------------
src/scheme/HyperplasticSolverO2.cpp | 40 +++++++++++--------------
2 files changed, 41 insertions(+), 45 deletions(-)
diff --git a/src/scheme/HyperplasticSolver.cpp b/src/scheme/HyperplasticSolver.cpp
index a5a7b4a61..d0c45c223 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 1d6174dd1..41543fca0 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;
--
GitLab