diff --git a/src/algebra/EigenvalueSolver.cpp b/src/algebra/EigenvalueSolver.cpp index 28ec7131d09db6b9f08e2616479bd912d4dc5799..dc653c1889a2c49fe970aad946a9fb1833c81db9 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 728cacdfa59843cffc6d3d45e2ed0f716e03722a..e01134ac8f53e522f8f45e9cb24a656176e90167 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 9e8a272782d708c10a5a33d7398fba5ef5599f1e..78e742d80373a9726fc55b7995d5e17d953f5c75 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 d2316d8845db03772ce0c0d7d56ec8b60bc3b96a..f9946d72ef5108658be809e531a30d0beb77c57f 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]); }