From d15023eef5e2bbc163cbaeddb24b6b80c4ca3b9a Mon Sep 17 00:00:00 2001
From: labourasse <labourassee@gmail.com>
Date: Wed, 20 Nov 2024 17:03:01 +0100
Subject: [PATCH] fix a bug in the plasticity management
---
src/algebra/EigenvalueSolver.cpp | 28 ++++++++++++----
src/algebra/EigenvalueSolver.hpp | 2 +-
src/scheme/HyperplasticSolver.cpp | 53 +++++++++----------------------
tests/test_EigenvalueSolver.cpp | 2 +-
4 files changed, 38 insertions(+), 47 deletions(-)
diff --git a/src/algebra/EigenvalueSolver.cpp b/src/algebra/EigenvalueSolver.cpp
index 28ec7131d..dc653c188 100644
--- a/src/algebra/EigenvalueSolver.cpp
+++ b/src/algebra/EigenvalueSolver.cpp
@@ -368,6 +368,10 @@ EigenvalueSolver::findEigenvectors(const TinyMatrix<3, 3>& A, TinyVector<3> eige
TinyMatrix<3> Try = findEigenvector(A, eigenvalues[count], space_size[count], epsilon);
// std::cout << eigenvalues[count] << " Frobenius norm try: " << frobeniusNorm(Try) << "\n";
Assert(frobeniusNorm(Try) > 1e-3, " Error : no eigenvector corresponds to eigenvalue ");
+ // if (frobeniusNorm(Try) < 1e-3) {
+ // std::cout << " Error : no eigenvector corresponds to eigenvalue \n";
+ // exit(0);
+ // }
TinyVector<3> eigenvector = zero;
int count2 = 0;
@@ -385,6 +389,10 @@ EigenvalueSolver::findEigenvectors(const TinyMatrix<3, 3>& A, TinyVector<3> eige
}
}
Assert(count2 == space_size[count - count2], " eigenvector space size is not what was expected ");
+ // if (count2 != space_size[count - count2]) {
+ // std::cout << " eigenvector space size is not what was expected \n";
+ // exit(0);
+ // }
it++;
}
// std ::cout << " eigenMatrix " << eigenMatrix << "\n";
@@ -477,10 +485,10 @@ EigenvalueSolver::_findFirstRoot(Polynomial<3> P) const
double old_guess = -P.coefficients()[2] / 3.;
// std::cout << " coefs P " << P.coefficients() << "\n";
// double delta = pow(2.*P.coefficients()[2],2) - 4*3.*P.coefficients()[3]*P.coefficients()[1];
- Polynomial<2> Q = derivative(P);
- Polynomial<1> R = derivative(Q);
- size_t iteration = 0;
- double residu = 0.;
+ Polynomial<2> Q = derivative(P);
+ Polynomial<1> R = derivative(Q);
+ // size_t iteration = 0;
+ double residu = 0.;
do {
// guess -= P(guess) / Q(guess);
old_guess = guess;
@@ -489,7 +497,7 @@ EigenvalueSolver::_findFirstRoot(Polynomial<3> P) const
// std::cout << "g(guess) = " << Q(guess) << "\n";
residu = P(guess);
// std::cout << "residu = " << residu << "\n";
- iteration += 1;
+ // iteration += 1;
} while ((fabs(residu) > 1.e-16) or (fabs(old_guess - guess) > 1.e-10));
// std::cout << " nb Newton iterations " << iteration << "\n";
return guess;
@@ -500,6 +508,10 @@ EigenvalueSolver::_findLastRoots(Polynomial<2> P, const double epsilon) const
{
TinyVector<2> roots(0., 0.);
double delta = pow(P.coefficients()[1], 2) - 4 * P.coefficients()[2] * P.coefficients()[0];
+ // if (fabs(delta) > epsilon) {
+ // std::cout << "Find roots is only for symmetric matrix \n";
+ // exit(0);
+ // }
Assert(delta > -epsilon, "Find roots is only for symmetric matrix");
if (fabs(delta) < epsilon) {
roots[0] = -P.coefficients()[1] / (2. * P.coefficients()[2]);
@@ -512,15 +524,17 @@ EigenvalueSolver::_findLastRoots(Polynomial<2> P, const double epsilon) const
return roots;
}
-void
-EigenvalueSolver::computeExpForTinyMatrix(const TinyMatrix<3>& A, TinyMatrix<3>& expA)
+TinyMatrix<3>
+EigenvalueSolver::computeExpForTinyMatrix(const TinyMatrix<3>& A)
{
+ TinyMatrix<3> expA;
auto [eigenvalues, eigenvectors] = findEigen(A);
TinyMatrix<3> D = zero;
for (size_t i = 0; i < 3; ++i) {
D(i, i) = std::exp(eigenvalues[i]);
}
expA = eigenvectors * D * transpose(eigenvectors);
+ return expA;
}
EigenvalueSolver::EigenvalueSolver() {}
diff --git a/src/algebra/EigenvalueSolver.hpp b/src/algebra/EigenvalueSolver.hpp
index 728cacdfa..e01134ac8 100644
--- a/src/algebra/EigenvalueSolver.hpp
+++ b/src/algebra/EigenvalueSolver.hpp
@@ -100,7 +100,7 @@ class EigenvalueSolver
TinyVector<2> _findLastRoots(Polynomial<2> P, const double epsilon) const;
- void computeExpForTinyMatrix(const TinyMatrix<3>& A, TinyMatrix<3>& expA);
+ TinyMatrix<3> computeExpForTinyMatrix(const TinyMatrix<3>& A);
EigenvalueSolver();
~EigenvalueSolver() = default;
diff --git a/src/scheme/HyperplasticSolver.cpp b/src/scheme/HyperplasticSolver.cpp
index 9e8a27278..78e742d80 100644
--- a/src/scheme/HyperplasticSolver.cpp
+++ b/src/scheme/HyperplasticSolver.cpp
@@ -170,46 +170,13 @@ class HyperplasticSolverHandler::HyperplasticSolver final : public HyperplasticS
return vector;
}
- R3x3
- matrix_exp(const R3x3 A) const
- {
- double normA = frobeniusNorm(A);
- R3x3 expA = identity;
- if (normA <= 1.e-12) {
- return expA;
- }
- SmallMatrix<double> expB{3};
- SmallMatrix<double> B(3, 3);
- // B.fill(0.);
- for (size_t i = 0; i < 3; i++) {
- for (size_t j = 0; j < 3; j++) {
- B(i, i) = A(i, j);
- }
- }
-
- std::cout << B << "\n";
- std::cout << std::flush;
-
- EigenvalueSolver{}.computeExpForSymmetricMatrix(B, expB);
- for (size_t i = 0; i < 3; i++) {
- for (size_t j = 0; j < 3; j++) {
- expA(i, j) = expB(i, j);
- }
- }
- std::cout << expA << "\n";
- return expA;
- }
-
double
_compute_chi(const R3x3 S, const double yield) const
{
const R3x3 Id = identity;
const R3x3 Sd = S - 1. / 3 * trace(S) * Id;
// auto [eigenvalues2, eigenmatrix2] = EigenvalueSolver{}.findEigen(TestB);
- R3x3 B;
- EigenvalueSolver{}.computeExpForTinyMatrix(TestA2, B);
- std::cout << B << "\n";
- exit(0);
+ // exit(0);
double chi = 0;
const double normS = frobeniusNorm(Sd);
if ((normS - std::sqrt(2. / 3) * yield) > 0.) {
@@ -612,7 +579,9 @@ class HyperplasticSolverHandler::HyperplasticSolver final : public HyperplasticS
CellValue<double> new_yield = copy(yield.cellValues());
CellValue<R3x3> sigma_s = copy(sigma.cellValues());
CellValue<double> chi{mesh.connectivity()};
+ CellValue<double> rationorm{mesh.connectivity()};
chi.fill(0.);
+ rationorm.fill(0.);
parallel_for(
mesh.numberOfCells(), PUGS_LAMBDA(CellId j) {
const auto& cell_nodes = cell_to_node_matrix[j];
@@ -626,18 +595,26 @@ class HyperplasticSolverHandler::HyperplasticSolver final : public HyperplasticS
momentum_fluxes += Fjr(j, R);
energy_fluxes += dot(Fjr(j, R), ur[r]);
}
+ R3x3 I = identity;
R3x3 gradv3d = to3D(gradv);
+ R3x3 sigmas = sigma[j] - 1. / 3 * trace(sigma[j]) * I;
const R3x3 elastic_fluxes = gradv3d * Fe[j];
const double dt_over_Mj = dt / (rho[j] * Vj[j]);
const double dt_over_Vj = dt / Vj[j];
new_u[j] += dt_over_Mj * momentum_fluxes;
new_E[j] += dt_over_Mj * energy_fluxes;
new_Fe[j] += dt_over_Vj * elastic_fluxes;
- chi[j] = _compute_chi(sigma[j], yield[j]);
- const R3x3 M = dt * chi[j] * zeta[j] * new_Fe[j];
- new_Fe[j] = new_Fe[j] * matrix_exp(M);
+ const double fenorm0 = frobeniusNorm(new_Fe[j]);
+ chi[j] = _compute_chi(sigma[j], yield[j]);
+ 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]);
+ rationorm[j] = fenorm1 / fenorm0;
});
- std::cout << "sum_chi " << sum(chi) << "\n";
+ std::cout << "sum_chi " << sum(chi) << " min norm ratio " << min(rationorm) << " max norm ratio " << max(rationorm)
+ << "\n";
+
CellValue<const double> new_Vj = MeshDataManager::instance().getMeshData(*new_mesh).Vj();
parallel_for(
diff --git a/tests/test_EigenvalueSolver.cpp b/tests/test_EigenvalueSolver.cpp
index d2316d884..f9946d72e 100644
--- a/tests/test_EigenvalueSolver.cpp
+++ b/tests/test_EigenvalueSolver.cpp
@@ -190,7 +190,7 @@ TEST_CASE("EigenvalueSolver", "[algebra]")
REQUIRE((B(i, j) - TestA2(i, j)) / frobeniusNorm(TestA2) == Catch::Approx(0).margin(1E-8));
}
}
- EigenvalueSolver{}.computeExpForTinyMatrix(TestA2, expA2);
+ expA2 = EigenvalueSolver{}.computeExpForTinyMatrix(TestA2);
for (size_t i = 0; i < 3; ++i) {
Diag(i, i) = std::exp(eigenvalues2[i]);
}
--
GitLab