From 909c6d5ff8e78f1b2239030d89a5dc515a255283 Mon Sep 17 00:00:00 2001 From: labourasse <labourassee@gmail.com> Date: Mon, 24 Feb 2025 11:07:24 +0100 Subject: [PATCH] projection on the symmetry boundary condition --- src/algebra/EigenvalueSolver.cpp | 51 +++------------ src/algebra/EigenvalueSolver.hpp | 10 +-- src/language/modules/SchemeModule.cpp | 41 ++++-------- src/language/modules/SchemeModule.hpp | 32 ---------- src/mesh/MeshFlatNodeBoundary.cpp | 12 ++-- src/scheme/HyperplasticSolverO2.cpp | 56 ++++++++++++++--- src/scheme/HyperplasticSolverO2.hpp | 3 +- src/scheme/PolynomialReconstruction.cpp | 9 --- tests/test_EigenvalueSolver.cpp | 83 +------------------------ 9 files changed, 79 insertions(+), 218 deletions(-) diff --git a/src/algebra/EigenvalueSolver.cpp b/src/algebra/EigenvalueSolver.cpp index edd217399..fbe0e344b 100644 --- a/src/algebra/EigenvalueSolver.cpp +++ b/src/algebra/EigenvalueSolver.cpp @@ -234,43 +234,6 @@ EigenvalueSolver::_slepscComputeForSymmetricMatrix(const CRSMatrix<double>& A, S { Internals::slepscComputeForSymmetricMatrix(A, eigenvalues); } - -void -EigenvalueSolver::computeExpForSymmetricMatrix(const PETScAijMatrixEmbedder& A, SmallMatrix<double>& expA) -{ - EPS eps; - - EPSCreate(PETSC_COMM_SELF, &eps); - EPSSetOperators(eps, A, nullptr); - EPSSetProblemType(eps, EPS_HEP); - - Internals::computeAllEigenvaluesOfSymmetricMatrix(eps); - - PetscInt nb_eigenvalues; - EPSGetDimensions(eps, &nb_eigenvalues, nullptr, nullptr); - - SmallArray<double> eigenvalues(nb_eigenvalues); - SmallMatrix<double> P(nb_eigenvalues, nb_eigenvalues); - SmallMatrix<double> D(nb_eigenvalues, nb_eigenvalues); - D = zero; - expA = SmallMatrix<double>(nb_eigenvalues, nb_eigenvalues); - Array<double> eigenvector(nb_eigenvalues); - for (PetscInt i = 0; i < nb_eigenvalues; ++i) { - Vec Vr; - VecCreateSeqWithArray(PETSC_COMM_SELF, 1, A.numberOfRows(), &(eigenvector[0]), &Vr); - EPSGetEigenpair(eps, i, &(eigenvalues[i]), nullptr, Vr, nullptr); - VecDestroy(&Vr); - for (size_t j = 0; j < eigenvector.size(); ++j) { - P(j, i) = eigenvector[j]; - } - } - - for (size_t i = 0; i < eigenvalues.size(); ++i) { - D(i, i) = std::exp(eigenvalues[i]); - } - expA = P * D * transpose(P); - EPSDestroy(&eps); -} #endif // PUGS_HAS_SLEPC template <typename T> @@ -329,7 +292,7 @@ EigenvalueSolver::findEigenvector(const TinyMatrix<3, 3>& A, if (space_size == 3) { // std::cout << "space_size = 3 " // << "\n"; - return eigenvectors0; + return eigenvectors; } TinyMatrix<3, 3> B = A; @@ -485,15 +448,17 @@ EigenvalueSolver::findEigenvectors(const TinyMatrix<3, 3>& A, TinyVector<3> eige std::tuple<TinyVector<3>, TinyMatrix<3>> EigenvalueSolver::findEigen(const TinyMatrix<3> A) { - const double epsilon = 1.e-6; // * l2Norm(eigenvalues); + constexpr TinyVector<3> eigenvalues0 = zero; + constexpr TinyMatrix<3> eigenvectors0 = identity; + const double epsilon = 1.e-6; // * l2Norm(eigenvalues); if (frobeniusNorm(A) < 1.e-15) { // std::cout << " Frobenius norm 0 " << frobeniusNorm(A) << "\n"; return std::make_tuple(eigenvalues0, eigenvectors0); } TinyMatrix<3> C = 1. / frobeniusNorm(A) * A; if (isDiagonal(C, epsilon)) { - // std::cout << "Matrix C " << C << " is diagonal " - // << "\n"; + std::cout << "Matrix C " << C << " is diagonal " + << "\n"; const TinyVector<3> eigenvalues(A(0, 0), A(1, 1), A(2, 2)); TinyMatrix<3> Diag = zero; for (size_t i = 0; i < 3; ++i) { @@ -630,7 +595,9 @@ EigenvalueSolver::computeExpForTinyMatrix(const TinyMatrix<3>& A) return expA; } -EigenvalueSolver::EigenvalueSolver() {} +#ifdef PUGS_HAS_SLEPC + +void EigenvalueSolver::_slepscComputeForSymmetricMatrix(const DenseMatrixWrapper<double>& A, SmallArray<double>& eigenvalues) { Internals::slepscComputeForSymmetricMatrix(A, eigenvalues); diff --git a/src/algebra/EigenvalueSolver.hpp b/src/algebra/EigenvalueSolver.hpp index 1c0024211..986961522 100644 --- a/src/algebra/EigenvalueSolver.hpp +++ b/src/algebra/EigenvalueSolver.hpp @@ -16,18 +16,11 @@ class EigenvalueSolver { private: struct Internals; - TinyMatrix<3> eigenvectors0{1, 0, 0, 0, 1, 0, 0, 0, 1}; - TinyVector<3> eigenvalues0{0, 0, 0}; const EigenvalueSolverOptions m_options; void _slepscComputeForSymmetricMatrix(const CRSMatrix<double>& A, SmallArray<double>& eigenvalues); - void computeForSymmetricMatrix(const PETScAijMatrixEmbedder& A, - SmallArray<double>& eigenvalues, - SmallMatrix<double>& P); - void computeExpForSymmetricMatrix(const PETScAijMatrixEmbedder& A, SmallMatrix<double>& expA); -#endif // PUGS_HAS_SLEPC void _slepscComputeForSymmetricMatrix(const DenseMatrixWrapper<double>& A, SmallArray<double>& eigenvalues); void _slepscComputeForSymmetricMatrix(const CRSMatrix<double>& A, @@ -142,7 +135,7 @@ class EigenvalueSolver #ifdef PUGS_HAS_SLEPC this->computeExpForSymmetricMatrix(PETScAijMatrixEmbedder{A}, expA); #else // PUGS_HAS_SLEPC - throw NotImplementedError("SLEPc is required to solve eigenvalue problems"); + throw NotImplementedError("SLEPc is required to solve eigenvalue problems"); #endif // PUGS_HAS_SLEPC } @@ -170,7 +163,6 @@ class EigenvalueSolver TinyMatrix<3> computeExpForTinyMatrix(const TinyMatrix<3>& A); - EigenvalueSolver(); EigenvalueSolver(const EigenvalueSolverOptions& options = EigenvalueSolverOptions::default_options); ~EigenvalueSolver() = default; }; diff --git a/src/language/modules/SchemeModule.cpp b/src/language/modules/SchemeModule.cpp index 1d5346666..23456fc1d 100644 --- a/src/language/modules/SchemeModule.cpp +++ b/src/language/modules/SchemeModule.cpp @@ -77,7 +77,7 @@ SchemeModule::SchemeModule() this->_addTypeDescriptor(ast_node_data_type_from<std::shared_ptr<const IBoundaryConditionDescriptor>>); this->_addTypeDescriptor(ast_node_data_type_from<std::shared_ptr<const VariableBCDescriptor>>); - this->_addTypeDescriptor(ast_node_data_type_from<std::shared_ptr<const BasisType>>); + // this->_addTypeDescriptor(ast_node_data_type_from<std::shared_ptr<const BasisType>>); this->_addBuiltinFunction("P0", std::function( @@ -527,30 +527,6 @@ SchemeModule::SchemeModule() )); - this->_addBuiltinFunction("canonicalBasis", std::function( - - []() -> std::shared_ptr<const BasisType> { - return std::make_shared<BasisType>(BasisType::canonical); - } - - )); - - this->_addBuiltinFunction("taylorBasis", std::function( - - []() -> std::shared_ptr<const BasisType> { - return std::make_shared<BasisType>(BasisType::taylor); - } - - )); - - this->_addBuiltinFunction("lagrangeBasis", std::function( - - []() -> std::shared_ptr<const BasisType> { - return std::make_shared<BasisType>(BasisType::lagrange); - } - - )); - this->_addBuiltinFunction("apply_acoustic_fluxes", std::function( @@ -1356,6 +1332,8 @@ SchemeModule::SchemeModule() const std::shared_ptr<const DiscreteFunctionVariant>& sigma, // const std::shared_ptr<const ItemValueVariant>& ur, // const std::shared_ptr<const SubItemValuePerItemVariant>& Fjr, // + const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>& + bc_descriptor_list, // const double& dt) -> std::tuple<std::shared_ptr<const MeshVariant>, std::shared_ptr<const DiscreteFunctionVariant>, std::shared_ptr<const DiscreteFunctionVariant>, @@ -1364,7 +1342,8 @@ SchemeModule::SchemeModule() return HyperplasticSolverO2Handler{getCommonMesh({rho, u, E, Fe})} // .solver() .apply_fluxes(dt, rho, u, E, Fe, zeta, yield, sigma, ur, Fjr, - HyperplasticSolverO2Handler::RelaxationType::Implicit); + HyperplasticSolverO2Handler::RelaxationType::Implicit, + bc_descriptor_list); } )); @@ -1381,6 +1360,8 @@ SchemeModule::SchemeModule() const std::shared_ptr<const DiscreteFunctionVariant>& sigma, // const std::shared_ptr<const ItemValueVariant>& ur, // const std::shared_ptr<const SubItemValuePerItemVariant>& Fjr, // + const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>& + bc_descriptor_list, // const double& dt) -> std::tuple<std::shared_ptr<const MeshVariant>, std::shared_ptr<const DiscreteFunctionVariant>, std::shared_ptr<const DiscreteFunctionVariant>, @@ -1389,7 +1370,8 @@ SchemeModule::SchemeModule() return HyperplasticSolverO2Handler{getCommonMesh({rho, u, E, Fe})} // .solver() .apply_fluxes(dt, rho, u, E, Fe, zeta, yield, sigma, ur, Fjr, - HyperplasticSolverO2Handler::RelaxationType::CauchyGreen); + HyperplasticSolverO2Handler::RelaxationType::CauchyGreen, + bc_descriptor_list); } )); @@ -1406,6 +1388,8 @@ SchemeModule::SchemeModule() const std::shared_ptr<const DiscreteFunctionVariant>& sigma, // const std::shared_ptr<const ItemValueVariant>& ur, // const std::shared_ptr<const SubItemValuePerItemVariant>& Fjr, // + const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>& + bc_descriptor_list, // const double& dt) -> std::tuple<std::shared_ptr<const MeshVariant>, std::shared_ptr<const DiscreteFunctionVariant>, std::shared_ptr<const DiscreteFunctionVariant>, @@ -1414,7 +1398,8 @@ SchemeModule::SchemeModule() return HyperplasticSolverO2Handler{getCommonMesh({rho, u, E, Fe})} // .solver() .apply_fluxes(dt, rho, u, E, Fe, zeta, yield, sigma, ur, Fjr, - HyperplasticSolverO2Handler::RelaxationType::Exponential); + HyperplasticSolverO2Handler::RelaxationType::Exponential, + bc_descriptor_list); } )); diff --git a/src/language/modules/SchemeModule.hpp b/src/language/modules/SchemeModule.hpp index defa9a18f..6cd0b0342 100644 --- a/src/language/modules/SchemeModule.hpp +++ b/src/language/modules/SchemeModule.hpp @@ -3,38 +3,6 @@ #include <analysis/PolynomialBasis.hpp> #include <language/modules/BuiltinModule.hpp> -#include <language/utils/ASTNodeDataTypeTraits.hpp> -#include <utils/PugsMacros.hpp> - -class IBoundaryConditionDescriptor; -template <> -inline ASTNodeDataType ast_node_data_type_from<std::shared_ptr<const IBoundaryConditionDescriptor>> = - ASTNodeDataType::build<ASTNodeDataType::type_id_t>("boundary_condition"); - -template <> -inline ASTNodeDataType ast_node_data_type_from<std::shared_ptr<const BasisType>> = - ASTNodeDataType::build<ASTNodeDataType::type_id_t>("basis_type"); - -class VariableBCDescriptor; - -template <> -inline ASTNodeDataType ast_node_data_type_from<std::shared_ptr<const VariableBCDescriptor>> = - ASTNodeDataType::build<ASTNodeDataType::type_id_t>("variable_boundary_condition"); - -class DiscreteFunctionVariant; -template <> -inline ASTNodeDataType ast_node_data_type_from<std::shared_ptr<const DiscreteFunctionVariant>> = - ASTNodeDataType::build<ASTNodeDataType::type_id_t>("Vh"); - -class IDiscreteFunctionDescriptor; -template <> -inline ASTNodeDataType ast_node_data_type_from<std::shared_ptr<const IDiscreteFunctionDescriptor>> = - ASTNodeDataType::build<ASTNodeDataType::type_id_t>("Vh_type"); - -class IQuadratureDescriptor; -template <> -inline ASTNodeDataType ast_node_data_type_from<std::shared_ptr<const IQuadratureDescriptor>> = - ASTNodeDataType::build<ASTNodeDataType::type_id_t>("quadrature"); class SchemeModule : public BuiltinModule { diff --git a/src/mesh/MeshFlatNodeBoundary.cpp b/src/mesh/MeshFlatNodeBoundary.cpp index 5621a6230..647298e97 100644 --- a/src/mesh/MeshFlatNodeBoundary.cpp +++ b/src/mesh/MeshFlatNodeBoundary.cpp @@ -24,12 +24,12 @@ MeshFlatNodeBoundary<MeshType>::_checkBoundaryIsFlat(const TinyVector<MeshType:: } }); - // if (parallel::allReduceOr(is_bad)) { - // std::ostringstream ost; - // ost << "invalid boundary \"" << rang::fgB::yellow << this->m_ref_node_list.refId() << rang::style::reset - // << "\": boundary is not flat!"; - // throw NormalError(ost.str()); - // } + if (parallel::allReduceOr(is_bad)) { + std::ostringstream ost; + ost << "invalid boundary \"" << rang::fgB::yellow << this->m_ref_node_list.refId() << rang::style::reset + << "\": boundary is not flat!"; + throw NormalError(ost.str()); + } } template <> diff --git a/src/scheme/HyperplasticSolverO2.cpp b/src/scheme/HyperplasticSolverO2.cpp index c4efdd128..9a9eff713 100644 --- a/src/scheme/HyperplasticSolverO2.cpp +++ b/src/scheme/HyperplasticSolverO2.cpp @@ -614,6 +614,29 @@ class HyperplasticSolverO2Handler::HyperplasticSolverO2 final this->_applyVelocityBC(bc_list, Ar, br); } + void + _projectOnBoundary(const BoundaryConditionList& bc_list, NodeValue<Rd>& xr) const + { + for (const auto& boundary_condition : bc_list) { + std::visit( + [&](auto&& bc) { + using T = std::decay_t<decltype(bc)>; + if constexpr (std::is_same_v<SymmetryBoundaryCondition, T>) { + const Rd& n = bc.outgoingNormal(); + const Rd& o = bc.origin(); + + const Array<const NodeId>& node_list = bc.nodeList(); + parallel_for( + bc.numberOfNodes(), PUGS_LAMBDA(int r_number) { + const NodeId r = node_list[r_number]; + Rd& x = xr[node_list[r]]; + x -= dot(x - o, n) * n; + }); + } + }, + boundary_condition); + } + } NodeValue<const Rd> _computeUr(const MeshType& mesh, const NodeValue<Rdxd>& Ar, const NodeValue<Rd>& br) const { @@ -740,7 +763,8 @@ class HyperplasticSolverO2Handler::HyperplasticSolverO2 final const DiscreteScalarFunction& yield, const DiscreteTensorFunction3d& sigma, const NodeValue<const Rd>& ur, - const NodeValuePerCell<const Rd>& Fjr) const + const NodeValuePerCell<const Rd>& Fjr, + const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>& bc_descriptor_list) const { const auto& cell_to_node_matrix = mesh.connectivity().cellToNodeMatrix(); @@ -752,7 +776,8 @@ class HyperplasticSolverO2Handler::HyperplasticSolverO2 final NodeValue<Rd> new_xr = copy(mesh.xr()); parallel_for( mesh.numberOfNodes(), PUGS_LAMBDA(NodeId r) { new_xr[r] += dt * ur[r]; }); - + const BoundaryConditionList bc_list = this->_getBCList(mesh, bc_descriptor_list); + this->_projectOnBoundary(bc_list, new_xr); std::shared_ptr<const MeshType> new_mesh = std::make_shared<MeshType>(mesh.shared_connectivity(), new_xr); CellValue<const double> Vj = MeshDataManager::instance().getMeshData(mesh).Vj(); @@ -858,7 +883,8 @@ class HyperplasticSolverO2Handler::HyperplasticSolverO2 final const DiscreteScalarFunction& yield, const DiscreteTensorFunction3d& sigma, const NodeValue<const Rd>& ur, - const NodeValuePerCell<const Rd>& Fjr) const + const NodeValuePerCell<const Rd>& Fjr, + const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>& bc_descriptor_list) const { const auto& cell_to_node_matrix = mesh.connectivity().cellToNodeMatrix(); @@ -870,6 +896,8 @@ class HyperplasticSolverO2Handler::HyperplasticSolverO2 final NodeValue<Rd> new_xr = copy(mesh.xr()); parallel_for( mesh.numberOfNodes(), PUGS_LAMBDA(NodeId r) { new_xr[r] += dt * ur[r]; }); + const BoundaryConditionList bc_list = this->_getBCList(mesh, bc_descriptor_list); + this->_projectOnBoundary(bc_list, new_xr); std::shared_ptr<const MeshType> new_mesh = std::make_shared<MeshType>(mesh.shared_connectivity(), new_xr); @@ -948,7 +976,8 @@ class HyperplasticSolverO2Handler::HyperplasticSolverO2 final const DiscreteScalarFunction& yield, const DiscreteTensorFunction3d& sigma, const NodeValue<const Rd>& ur, - const NodeValuePerCell<const Rd>& Fjr) const + const NodeValuePerCell<const Rd>& Fjr, + const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>& bc_descriptor_list) const { const auto& cell_to_node_matrix = mesh.connectivity().cellToNodeMatrix(); @@ -960,6 +989,8 @@ class HyperplasticSolverO2Handler::HyperplasticSolverO2 final NodeValue<Rd> new_xr = copy(mesh.xr()); parallel_for( mesh.numberOfNodes(), PUGS_LAMBDA(NodeId r) { new_xr[r] += dt * ur[r]; }); + const BoundaryConditionList bc_list = this->_getBCList(mesh, bc_descriptor_list); + this->_projectOnBoundary(bc_list, new_xr); std::shared_ptr<const MeshType> new_mesh = std::make_shared<MeshType>(mesh.shared_connectivity(), new_xr); @@ -1036,7 +1067,8 @@ class HyperplasticSolverO2Handler::HyperplasticSolverO2 final const std::shared_ptr<const DiscreteFunctionVariant>& sigma, const std::shared_ptr<const ItemValueVariant>& ur, const std::shared_ptr<const SubItemValuePerItemVariant>& Fjr, - const RelaxationType& relaxation_type) const + const RelaxationType& relaxation_type, + const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>& bc_descriptor_list) const { std::shared_ptr mesh_v = getCommonMesh({rho, u, E}); if (not mesh_v) { @@ -1058,7 +1090,7 @@ class HyperplasticSolverO2Handler::HyperplasticSolverO2 final yield->get<DiscreteScalarFunction>(), // sigma->get<DiscreteTensorFunction3d>(), // ur->get<NodeValue<const Rd>>(), // - Fjr->get<NodeValuePerCell<const Rd>>()); + Fjr->get<NodeValuePerCell<const Rd>>(), bc_descriptor_list); } case RelaxationType::Implicit: { return this->apply_fluxes2(dt, // @@ -1071,7 +1103,7 @@ class HyperplasticSolverO2Handler::HyperplasticSolverO2 final yield->get<DiscreteScalarFunction>(), // sigma->get<DiscreteTensorFunction3d>(), // ur->get<NodeValue<const Rd>>(), // - Fjr->get<NodeValuePerCell<const Rd>>()); + Fjr->get<NodeValuePerCell<const Rd>>(), bc_descriptor_list); } case RelaxationType::CauchyGreen: { return this->apply_fluxes_B(dt, // @@ -1084,7 +1116,7 @@ class HyperplasticSolverO2Handler::HyperplasticSolverO2 final yield->get<DiscreteScalarFunction>(), // sigma->get<DiscreteTensorFunction3d>(), // ur->get<NodeValue<const Rd>>(), // - Fjr->get<NodeValuePerCell<const Rd>>()); + Fjr->get<NodeValuePerCell<const Rd>>(), bc_descriptor_list); } default: { throw UnexpectedError("invalid relaxation type"); @@ -1298,7 +1330,7 @@ class HyperplasticSolverO2Handler::HyperplasticSolverO2 final const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>& bc_descriptor_list) const { auto [ur, Fjr] = compute_fluxes(solver_type, rho, aL, aT, u, sigma, bc_descriptor_list); - return apply_fluxes(dt, rho, u, E, Fe, zeta, yield, sigma, ur, Fjr, relaxation_type); + return apply_fluxes(dt, rho, u, E, Fe, zeta, yield, sigma, ur, Fjr, relaxation_type, bc_descriptor_list); } std::tuple<std::shared_ptr<const MeshVariant>, @@ -1695,6 +1727,12 @@ class HyperplasticSolverO2Handler::HyperplasticSolverO2<MeshType>::SymmetryBound return m_mesh_flat_node_boundary.outgoingNormal(); } + const Rd& + origin() const + { + return m_mesh_flat_node_boundary.origin(); + } + size_t numberOfNodes() const { diff --git a/src/scheme/HyperplasticSolverO2.hpp b/src/scheme/HyperplasticSolverO2.hpp index ff89780d1..063dcc7d6 100644 --- a/src/scheme/HyperplasticSolverO2.hpp +++ b/src/scheme/HyperplasticSolverO2.hpp @@ -60,7 +60,8 @@ class HyperplasticSolverO2Handler const std::shared_ptr<const DiscreteFunctionVariant>& sigma, const std::shared_ptr<const ItemValueVariant>& ur, const std::shared_ptr<const SubItemValuePerItemVariant>& Fjr, - const RelaxationType& relaxation_type) const = 0; + const RelaxationType& relaxation_type, + const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>& bc_descriptor_list) const = 0; virtual std::tuple<std::shared_ptr<const MeshVariant>, std::shared_ptr<const DiscreteFunctionVariant>, diff --git a/src/scheme/PolynomialReconstruction.cpp b/src/scheme/PolynomialReconstruction.cpp index 110fa9c07..1a2b1e072 100644 --- a/src/scheme/PolynomialReconstruction.cpp +++ b/src/scheme/PolynomialReconstruction.cpp @@ -25,15 +25,6 @@ #include <scheme/DiscreteFunctionUtils.hpp> #include <scheme/DiscreteFunctionVariant.hpp> -template <size_t Dimension> -PUGS_INLINE auto -symmetrize_matrix(const TinyVector<Dimension>& normal, const TinyMatrix<Dimension>& M) -{ - const TinyMatrix<Dimension> Id = identity; - const TinyMatrix<Dimension> S = Id - 2. * tensorProduct(normal, normal); - return S * M * S; -} - template <size_t Dimension> PUGS_INLINE auto symmetrize_vector(const TinyVector<Dimension>& normal, const TinyVector<Dimension>& u) diff --git a/tests/test_EigenvalueSolver.cpp b/tests/test_EigenvalueSolver.cpp index 26acafc6f..e63f5d2e1 100644 --- a/tests/test_EigenvalueSolver.cpp +++ b/tests/test_EigenvalueSolver.cpp @@ -146,33 +146,6 @@ TEST_CASE("EigenvalueSolver", "[algebra]") SmallArray<double> eigenvalues; #ifdef PUGS_HAS_SLEPC - EigenvalueSolver{}.computeForSymmetricMatrix(A, eigenvalue_list, eigenvector_list); - - REQUIRE(eigenvalue_list[0] == Catch::Approx(-1)); - REQUIRE(eigenvalue_list[1] == Catch::Approx(-1)); - REQUIRE(eigenvalue_list[2] == Catch::Approx(8)); - - SmallMatrix<double> P{3}; - SmallMatrix<double> PT{3}; - SmallMatrix<double> D{3}; - D = zero; - - for (size_t i = 0; i < 3; ++i) { - for (size_t j = 0; j < 3; ++j) { - P(i, j) = eigenvector_list[j][i]; - PT(i, j) = eigenvector_list[i][j]; - } - D(i, i) = eigenvalue_list[i]; - } - SmallMatrix<double> eigen{3}; - eigen = zero; - - SmallMatrix PDPT = P * D * PT; - for (size_t i = 0; i < 3; ++i) { - for (size_t j = 0; j < 3; ++j) { - REQUIRE(PDPT(i, j) - S(i, j) == Catch::Approx(0).margin(1E-13)); - } - } EigenvalueSolver{options}.computeForSymmetricMatrix(A, eigenvalues); REQUIRE(eigenvalues[0] == Catch::Approx(-1)); REQUIRE(eigenvalues[1] == Catch::Approx(-1)); @@ -273,32 +246,6 @@ TEST_CASE("EigenvalueSolver", "[algebra]") SmallArray<double> eigenvalues; #ifdef PUGS_HAS_SLEPC - EigenvalueSolver{}.computeForSymmetricMatrix(A, eigenvalue_list, P); - - REQUIRE(eigenvalue_list[0] == Catch::Approx(-1)); - REQUIRE(eigenvalue_list[1] == Catch::Approx(-1)); - REQUIRE(eigenvalue_list[2] == Catch::Approx(8)); - - SmallMatrix<double> D{3}; - D = zero; - for (size_t i = 0; i < 3; ++i) { - D(i, i) = eigenvalue_list[i]; - } - SmallMatrix PT = transpose(P); - - TinyMatrix<3, 3, double> TinyA; - for (size_t i = 0; i < 3; ++i) { - for (size_t j = 0; j < 3; ++j) { - TinyA(i, j) = S(i, j); - } - } - - SmallMatrix PDPT = P * D * PT; - for (size_t i = 0; i < 3; ++i) { - for (size_t j = 0; j < 3; ++j) { - REQUIRE(PDPT(i, j) - S(i, j) == Catch::Approx(0).margin(1E-13)); - } - } EigenvalueSolver{options}.computeForSymmetricMatrix(A, eigenvalues); REQUIRE(eigenvalues[0] == Catch::Approx(-1)); REQUIRE(eigenvalues[1] == Catch::Approx(-1)); @@ -687,36 +634,8 @@ TEST_CASE("EigenvalueSolver", "[algebra]") #endif // PUGS_HAS_EIGEN3 } } - SECTION("exponential of a matrix") - { - SmallMatrix<double> expA{}; - SmallMatrix<double> expS{3}; - - expS(0, 0) = 1325.074593930307812; - expS(0, 1) = 662.353357244568285; - expS(0, 2) = 1324.70671448913637; - expS(1, 0) = expS(0, 1); - expS(1, 1) = 331.544558063455535; - expS(1, 2) = 662.353357244568185; - expS(2, 0) = expS(0, 2); - expS(2, 1) = expS(1, 2); - expS(2, 2) = expS(0, 0); - -#ifdef PUGS_HAS_SLEPC - EigenvalueSolver{}.computeExpForSymmetricMatrix(A, expA); - - for (size_t i = 0; i < 3; ++i) { - for (size_t j = 0; j < 3; ++j) { - REQUIRE(expA(i, j) - expS(i, j) == Catch::Approx(0).margin(1E-12)); - } - } - -#else // PUGS_HAS_SLEPC - REQUIRE_THROWS_WITH(EigenvalueSolver{}.computeExpForSymmetricMatrix(A, expA), - "not implemented yet: SLEPc is required to solve eigenvalue problems"); -#endif // PUGS_HAS_SLEPC - } } + SECTION("symmetric tiny matrix") { TinyMatrix<3> TestA{3e10, 2e10, 4e10, 2e10, 0, 2e10, 4e10, 2e10, 3e10}; -- GitLab