From 81b2e988c0d80209fa4586af08b0d6525b9f4de6 Mon Sep 17 00:00:00 2001 From: Max Firmbach Date: Sat, 22 Nov 2025 12:44:49 +0100 Subject: [PATCH 1/4] Work on thyra converters and iterative solver --- .../src/sparse/4C_linalg_sparsematrix.hpp | 1 + .../src/sparse/4C_linalg_sparseoperator.hpp | 1 + ...4C_linear_solver_amgnxn_preconditioner.hpp | 6 + .../4C_linear_solver_method_iterative.cpp | 44 ++-- ...4C_linear_solver_preconditioner_ifpack.cpp | 4 +- ...4C_linear_solver_preconditioner_ifpack.hpp | 7 + .../4C_linear_solver_preconditioner_muelu.cpp | 9 +- .../4C_linear_solver_preconditioner_muelu.hpp | 7 + ...inear_solver_preconditioner_projection.hpp | 7 + .../4C_linear_solver_preconditioner_teko.cpp | 4 +- .../4C_linear_solver_preconditioner_teko.hpp | 7 + .../4C_linear_solver_preconditioner_type.hpp | 3 + .../utils/4C_linear_solver_thyra_utils.cpp | 201 ++++++++++++++++++ .../utils/4C_linear_solver_thyra_utils.hpp | 24 ++- 14 files changed, 295 insertions(+), 30 deletions(-) diff --git a/src/core/linalg/src/sparse/4C_linalg_sparsematrix.hpp b/src/core/linalg/src/sparse/4C_linalg_sparsematrix.hpp index 7029be930f9..84171a69b8d 100644 --- a/src/core/linalg/src/sparse/4C_linalg_sparsematrix.hpp +++ b/src/core/linalg/src/sparse/4C_linalg_sparsematrix.hpp @@ -18,6 +18,7 @@ #include #include +#include FOUR_C_NAMESPACE_OPEN diff --git a/src/core/linalg/src/sparse/4C_linalg_sparseoperator.hpp b/src/core/linalg/src/sparse/4C_linalg_sparseoperator.hpp index 7bce011bfcb..065aaccd1d0 100644 --- a/src/core/linalg/src/sparse/4C_linalg_sparseoperator.hpp +++ b/src/core/linalg/src/sparse/4C_linalg_sparseoperator.hpp @@ -16,6 +16,7 @@ #include "4C_utils_shared_ptr_from_ref.hpp" #include +#include #include #include diff --git a/src/core/linear_solver/src/amgnxn/4C_linear_solver_amgnxn_preconditioner.hpp b/src/core/linear_solver/src/amgnxn/4C_linear_solver_amgnxn_preconditioner.hpp index ef2620ac119..1069ff58ccb 100644 --- a/src/core/linear_solver/src/amgnxn/4C_linear_solver_amgnxn_preconditioner.hpp +++ b/src/core/linear_solver/src/amgnxn/4C_linear_solver_amgnxn_preconditioner.hpp @@ -41,6 +41,12 @@ namespace Core::LinearSolver /// linear operator used for preconditioning std::shared_ptr prec_operator() const override; + Teuchos::RCP> thyra_operator() const override + { + FOUR_C_THROW("One should not end up here."); + } + + private: // Private variables std::shared_ptr p_; // The underlying preconditioner object diff --git a/src/core/linear_solver/src/method/4C_linear_solver_method_iterative.cpp b/src/core/linear_solver/src/method/4C_linear_solver_method_iterative.cpp index 7581c6732cf..641e652debc 100644 --- a/src/core/linear_solver/src/method/4C_linear_solver_method_iterative.cpp +++ b/src/core/linear_solver/src/method/4C_linear_solver_method_iterative.cpp @@ -14,15 +14,16 @@ #include "4C_linear_solver_preconditioner_muelu.hpp" #include "4C_linear_solver_preconditioner_projection.hpp" #include "4C_linear_solver_preconditioner_teko.hpp" +#include "4C_linear_solver_thyra_utils.hpp" #include "4C_utils_exceptions.hpp" #include #include #include #include -#include #include #include +#include #include #include #include @@ -30,9 +31,6 @@ FOUR_C_NAMESPACE_OPEN -using BelosVectorType = Epetra_MultiVector; -using BelosMatrixType = Epetra_Operator; - //---------------------------------------------------------------------------------- //---------------------------------------------------------------------------------- Core::LinearSolver::IterativeSolver::IterativeSolver(MPI_Comm comm, Teuchos::ParameterList& params) @@ -89,14 +87,20 @@ int Core::LinearSolver::IterativeSolver::solve(Core::LinAlg::MultiVector { Teuchos::ParameterList& belist = params().sublist("Belos Parameters"); - auto problem = Teuchos::make_rcp>( - Teuchos::rcpFromRef(a_->epetra_operator()), Teuchos::rcpFromRef(x.get_epetra_multi_vector()), - Teuchos::rcpFromRef(b_->get_epetra_multi_vector())); + using MV = Thyra::MultiVectorBase; + using OP = Thyra::LinearOpBase; + using belos_problem = Belos::LinearProblem; + + auto a = Utils::create_thyra_linear_op(*a_, LinAlg::DataAccess::Share); + auto solution = Utils::create_thyra_multi_vector(x, a->range()); + auto b = Utils::create_thyra_multi_vector(*b_, a->domain()); + + auto problem = Teuchos::make_rcp(a, solution, b); if (preconditioner_ != nullptr) { - auto belosPrec = - Teuchos::make_rcp(Teuchos::rcp(preconditioner_->prec_operator())); + // TODO: Preconditioner seems to work not that properly (Blocked MueLu and Krylov Projection) + auto belosPrec = preconditioner_->thyra_operator(); problem->setRightPrec(belosPrec); } @@ -104,7 +108,7 @@ int Core::LinearSolver::IterativeSolver::solve(Core::LinAlg::MultiVector if (set == false) FOUR_C_THROW("Core::LinearSolver::BelosSolver: Iterative solver failed to set up correctly."); - std::shared_ptr> newSolver; + std::shared_ptr> newSolver; if (belist.isParameter("SOLVER_XML_FILE")) { @@ -123,8 +127,7 @@ int Core::LinearSolver::IterativeSolver::solve(Core::LinAlg::MultiVector } newSolver = - std::make_shared>( - problem, belosSolverList); + std::make_shared>(problem, belosSolverList); } else if (belosParams.isSublist("CG")) { @@ -135,8 +138,7 @@ int Core::LinearSolver::IterativeSolver::solve(Core::LinAlg::MultiVector } newSolver = - std::make_shared>( - problem, belosSolverList); + std::make_shared>(problem, belosSolverList); } else if (belosParams.isSublist("BiCGSTAB")) { @@ -146,8 +148,7 @@ int Core::LinearSolver::IterativeSolver::solve(Core::LinAlg::MultiVector belosSolverList->set("Convergence Tolerance", belist.get("Convergence Tolerance")); } - newSolver = std::make_shared>( - problem, belosSolverList); + newSolver = std::make_shared>(problem, belosSolverList); } else FOUR_C_THROW("Core::LinearSolver::BelosSolver: Unknown iterative solver solver type chosen."); @@ -161,14 +162,13 @@ int Core::LinearSolver::IterativeSolver::solve(Core::LinAlg::MultiVector std::string solverType = belist.get("Solver Type"); if (solverType == "GMRES") - newSolver = - std::make_shared>( - problem, Teuchos::rcpFromRef(belist)); + newSolver = std::make_shared>( + problem, Teuchos::rcpFromRef(belist)); else if (solverType == "CG") - newSolver = std::make_shared>( + newSolver = std::make_shared>( problem, Teuchos::rcpFromRef(belist)); else if (solverType == "BiCGSTAB") - newSolver = std::make_shared>( + newSolver = std::make_shared>( problem, Teuchos::rcpFromRef(belist)); else FOUR_C_THROW("Core::LinearSolver::BelosSolver: Unknown iterative solver solver type chosen."); @@ -204,6 +204,8 @@ int Core::LinearSolver::IterativeSolver::solve(Core::LinAlg::MultiVector } } + x = *Utils::get_linalg_multi_vector_from_thyra(a_->domain_map(), solution); + return 0; } diff --git a/src/core/linear_solver/src/preconditioner/4C_linear_solver_preconditioner_ifpack.cpp b/src/core/linear_solver/src/preconditioner/4C_linear_solver_preconditioner_ifpack.cpp index 5b68f598685..bc79e7ebfa0 100644 --- a/src/core/linear_solver/src/preconditioner/4C_linear_solver_preconditioner_ifpack.cpp +++ b/src/core/linear_solver/src/preconditioner/4C_linear_solver_preconditioner_ifpack.cpp @@ -71,9 +71,9 @@ void Core::LinearSolver::IFPACKPreconditioner::setup( builder.createPreconditioningStrategy("Ifpack"); Teuchos::RCP> prec = Thyra::prec(*precFactory, pmatrix_); - auto inverseOp = prec->getUnspecifiedPrecOp(); + p_thyra_ = prec->getUnspecifiedPrecOp(); - p_ = std::make_shared(inverseOp); + p_ = std::make_shared(p_thyra_); } FOUR_C_NAMESPACE_CLOSE diff --git a/src/core/linear_solver/src/preconditioner/4C_linear_solver_preconditioner_ifpack.hpp b/src/core/linear_solver/src/preconditioner/4C_linear_solver_preconditioner_ifpack.hpp index 4b816e2fe59..ac7a11b481f 100644 --- a/src/core/linear_solver/src/preconditioner/4C_linear_solver_preconditioner_ifpack.hpp +++ b/src/core/linear_solver/src/preconditioner/4C_linear_solver_preconditioner_ifpack.hpp @@ -35,6 +35,11 @@ namespace Core::LinearSolver /// linear operator used for preconditioning std::shared_ptr prec_operator() const override { return p_; } + Teuchos::RCP> thyra_operator() const override + { + return p_thyra_; + } + private: //! IFPACK parameter list Teuchos::ParameterList& ifpacklist_; @@ -42,6 +47,8 @@ namespace Core::LinearSolver //! system of equations used for preconditioning used by P_ only Teuchos::RCP> pmatrix_; + Teuchos::RCP> p_thyra_; + //! preconditioner std::shared_ptr p_; }; diff --git a/src/core/linear_solver/src/preconditioner/4C_linear_solver_preconditioner_muelu.cpp b/src/core/linear_solver/src/preconditioner/4C_linear_solver_preconditioner_muelu.cpp index 5599f59fd4e..a213925202e 100644 --- a/src/core/linear_solver/src/preconditioner/4C_linear_solver_preconditioner_muelu.cpp +++ b/src/core/linear_solver/src/preconditioner/4C_linear_solver_preconditioner_muelu.cpp @@ -22,6 +22,9 @@ #include #include #include +#include +#include +#include #include #include #include @@ -116,9 +119,9 @@ void Core::LinearSolver::MueLuPreconditioner::setup( builder.createPreconditioningStrategy("MueLu"); Teuchos::RCP> prec = Thyra::prec(*precFactory, pmatrix_); - auto inverseOp = prec->getUnspecifiedPrecOp(); - - p_ = std::make_shared(inverseOp); + p_thyra_ = prec->getUnspecifiedPrecOp(); + + p_ = std::make_shared(p_thyra_); } else { diff --git a/src/core/linear_solver/src/preconditioner/4C_linear_solver_preconditioner_muelu.hpp b/src/core/linear_solver/src/preconditioner/4C_linear_solver_preconditioner_muelu.hpp index f872475c6f2..3d965b90f8a 100644 --- a/src/core/linear_solver/src/preconditioner/4C_linear_solver_preconditioner_muelu.hpp +++ b/src/core/linear_solver/src/preconditioner/4C_linear_solver_preconditioner_muelu.hpp @@ -61,6 +61,11 @@ namespace Core::LinearSolver //! linear operator used for preconditioning std::shared_ptr prec_operator() const final { return p_; } + Teuchos::RCP> thyra_operator() const override + { + return p_thyra_; + } + private: //! system of equations used for preconditioning used by P_ only Teuchos::RCP> pmatrix_; @@ -72,6 +77,8 @@ namespace Core::LinearSolver //! preconditioner std::shared_ptr p_; + Teuchos::RCP> p_thyra_; + //! MueLu hierarchy Teuchos::RCP> H_; }; diff --git a/src/core/linear_solver/src/preconditioner/4C_linear_solver_preconditioner_projection.hpp b/src/core/linear_solver/src/preconditioner/4C_linear_solver_preconditioner_projection.hpp index f313b4c950c..7dadb190eea 100644 --- a/src/core/linear_solver/src/preconditioner/4C_linear_solver_preconditioner_projection.hpp +++ b/src/core/linear_solver/src/preconditioner/4C_linear_solver_preconditioner_projection.hpp @@ -34,6 +34,11 @@ namespace Core::LinearSolver /// linear operator used for preconditioning std::shared_ptr prec_operator() const override { return p_; } + Teuchos::RCP> thyra_operator() const override + { + FOUR_C_THROW("One should not end up here."); + } + private: std::shared_ptr preconditioner_; @@ -41,6 +46,8 @@ namespace Core::LinearSolver std::shared_ptr projector_; std::shared_ptr p_; + + Teuchos::RCP> p_thyra_; }; } // namespace Core::LinearSolver diff --git a/src/core/linear_solver/src/preconditioner/4C_linear_solver_preconditioner_teko.cpp b/src/core/linear_solver/src/preconditioner/4C_linear_solver_preconditioner_teko.cpp index b68cfeecf18..cc6d95ab920 100644 --- a/src/core/linear_solver/src/preconditioner/4C_linear_solver_preconditioner_teko.cpp +++ b/src/core/linear_solver/src/preconditioner/4C_linear_solver_preconditioner_teko.cpp @@ -147,9 +147,9 @@ void Core::LinearSolver::TekoPreconditioner::setup( builder.createPreconditioningStrategy("Teko"); Teuchos::RCP> prec = Thyra::prec(*precFactory, pmatrix_); - Teko::LinearOp inverseOp = prec->getUnspecifiedPrecOp(); + p_thyra_ = prec->getUnspecifiedPrecOp(); - p_ = std::make_shared(inverseOp); + p_ = std::make_shared(p_thyra_); } //---------------------------------------------------------------------------------- diff --git a/src/core/linear_solver/src/preconditioner/4C_linear_solver_preconditioner_teko.hpp b/src/core/linear_solver/src/preconditioner/4C_linear_solver_preconditioner_teko.hpp index 7c9cfc7b215..3dcf68c32b4 100644 --- a/src/core/linear_solver/src/preconditioner/4C_linear_solver_preconditioner_teko.hpp +++ b/src/core/linear_solver/src/preconditioner/4C_linear_solver_preconditioner_teko.hpp @@ -38,6 +38,11 @@ namespace Core::LinearSolver /// linear operator used for preconditioning std::shared_ptr prec_operator() const override { return p_; } + Teuchos::RCP> thyra_operator() const override + { + return p_thyra_; + } + private: Teuchos::ParameterList& tekolist_; @@ -46,6 +51,8 @@ namespace Core::LinearSolver //! preconditioner std::shared_ptr p_; + + Teuchos::RCP> p_thyra_; }; diff --git a/src/core/linear_solver/src/preconditioner/4C_linear_solver_preconditioner_type.hpp b/src/core/linear_solver/src/preconditioner/4C_linear_solver_preconditioner_type.hpp index 63d9b255ed1..55ccbf4a51c 100644 --- a/src/core/linear_solver/src/preconditioner/4C_linear_solver_preconditioner_type.hpp +++ b/src/core/linear_solver/src/preconditioner/4C_linear_solver_preconditioner_type.hpp @@ -14,6 +14,7 @@ #include "4C_linalg_vector.hpp" #include +#include #include @@ -43,6 +44,8 @@ namespace Core::LinearSolver virtual void setup( Core::LinAlg::SparseOperator& matrix, Core::LinAlg::MultiVector& b) = 0; + virtual Teuchos::RCP> thyra_operator() const = 0; + /// linear operator used for preconditioning virtual std::shared_ptr prec_operator() const = 0; }; diff --git a/src/core/linear_solver/src/utils/4C_linear_solver_thyra_utils.cpp b/src/core/linear_solver/src/utils/4C_linear_solver_thyra_utils.cpp index dff66e153ee..6659f20d5fd 100644 --- a/src/core/linear_solver/src/utils/4C_linear_solver_thyra_utils.cpp +++ b/src/core/linear_solver/src/utils/4C_linear_solver_thyra_utils.cpp @@ -7,9 +7,13 @@ #include "4C_linear_solver_thyra_utils.hpp" +#include +#include #include #include #include +#include +#include FOUR_C_NAMESPACE_OPEN @@ -34,6 +38,121 @@ Teuchos::RCP> Core::LinearSolver::Utils::create_t } +//------------------------------------------------------------------------------ +//------------------------------------------------------------------------------ +void blockEpetraToThyra(int numVectors, double* epetraData, int leadingDim, + const Teuchos::Ptr>& mv, int& localDim) +{ + localDim = 0; + + // check the base case + const Teuchos::Ptr> prodMV = + ptr_dynamic_cast>(mv); + if (prodMV == Teuchos::null) + { + // VS object must be a SpmdMultiVector object + const Teuchos::Ptr> spmdX = + ptr_dynamic_cast>(mv, true); + const Teuchos::RCP> spmdVS = spmdX->spmdSpace(); + + int localSubDim = spmdVS->localSubDim(); + + Thyra::Ordinal thyraLeadingDim = 0; + + Teuchos::ArrayRCP thyraData_arcp; + Teuchos::ArrayView thyraData; + spmdX->getNonconstLocalData(Teuchos::outArg(thyraData_arcp), Teuchos::outArg(thyraLeadingDim)); + thyraData = thyraData_arcp(); // build array view + + for (int i = 0; i < localSubDim; i++) + { + // copy each vector + for (int v = 0; v < numVectors; v++) + { + thyraData[i + thyraLeadingDim * v] = epetraData[i + leadingDim * v]; + } + } + + localDim = localSubDim; + + return; + } + + // this keeps track of current location in the epetraData vector + double* localData = epetraData; + + // loop over all the blocks in the vector space + for (int blkIndex = 0; blkIndex < prodMV->productSpace()->numBlocks(); blkIndex++) + { + int subDim = 0; + const Teuchos::RCP> blockVec = + prodMV->getNonconstMultiVectorBlock(blkIndex); + + // perorm the recursive copy + blockEpetraToThyra(numVectors, localData, leadingDim, blockVec.ptr(), subDim); + + // shift to the next block + localData += subDim; + + // account for the size of this subblock + localDim += subDim; + } +} + +Teuchos::RCP> Core::LinearSolver::Utils::create_thyra_multi_vector( + const Core::LinAlg::MultiVector& multi_vector, + Teuchos::RCP> map) +{ + // TODO: We might need to cast if we have a product vector space and handle things differently ... + // If we have a product space, we will also need to get a product vector. + auto product_map = Teuchos::rcp_dynamic_cast>(map); + + if (product_map.is_null()) + { + auto const_thyra_vector = + Thyra::create_MultiVector(Teuchos::rcpFromRef(multi_vector.get_epetra_multi_vector()), map); + + return Teuchos::rcp_const_cast>(const_thyra_vector); + } + else + { + const int num_vectors = multi_vector.num_vectors(); + + auto thyra_multi_vector = Thyra::defaultProductMultiVector(product_map, num_vectors); + + // extract local information from the Epetra_MultiVector + int leadingDim = 0, localDim = 0; + double* epetraData = 0; + multi_vector.get_epetra_multi_vector().ExtractView(&epetraData, &leadingDim); + + blockEpetraToThyra(num_vectors, epetraData, leadingDim, thyra_multi_vector.ptr(), localDim); + + return thyra_multi_vector; + } +} + + +//------------------------------------------------------------------------------ +//------------------------------------------------------------------------------ +Teuchos::RCP> Core::LinearSolver::Utils::create_thyra_linear_op( + Core::LinAlg::SparseOperator& matrix, Core::LinAlg::DataAccess access) +{ + auto block_matrix = std::dynamic_pointer_cast( + Core::Utils::shared_ptr_from_ref(matrix)); + + if (block_matrix == nullptr) + { + auto sparse_matrix = std::dynamic_pointer_cast( + Core::Utils::shared_ptr_from_ref(matrix)); + return create_thyra_linear_op(*sparse_matrix, access); + } + else + { + return create_thyra_linear_op(*block_matrix, access); + } +} + + //------------------------------------------------------------------------------ //------------------------------------------------------------------------------ Teuchos::RCP> Core::LinearSolver::Utils::create_thyra_linear_op( @@ -104,4 +223,86 @@ Teuchos::RCP Core::LinearSolver::Utils::get_epetra_vector_f return ::Thyra::get_Epetra_Vector(map.get_epetra_map(), thyra_vector); } + +//------------------------------------------------------------------------------ +//------------------------------------------------------------------------------ +void blockThyraToEpetra(int numVectors, double* epetraData, int leadingDim, + const Teuchos::RCP>& tX, int& localDim) +{ + localDim = 0; + + // check the base case + const Teuchos::RCP> prodX = + rcp_dynamic_cast>(tX); + if (prodX == Teuchos::null) + { + // the base case + + // VS object must be a SpmdMultiVector object + Teuchos::RCP> spmdX = + rcp_dynamic_cast>(tX, true); + Teuchos::RCP> spmdVS = spmdX->spmdSpace(); + + int localSubDim = spmdVS->localSubDim(); + + Thyra::Ordinal thyraLeadingDim = 0; + + Teuchos::ArrayView thyraData; + Teuchos::ArrayRCP thyraData_arcp; + spmdX->getLocalData(Teuchos::outArg(thyraData_arcp), Teuchos::outArg(thyraLeadingDim)); + thyraData = thyraData_arcp(); // grab the array view + + for (int i = 0; i < localSubDim; i++) + { + // copy each vector + for (int v = 0; v < numVectors; v++) + epetraData[i + leadingDim * v] = thyraData[i + thyraLeadingDim * v]; + } + + localDim = localSubDim; + + return; + } + + const Teuchos::RCP> prodVS = prodX->productSpace(); + + // this keeps track of current location in the epetraData vector + double* localData = epetraData; + + // loop over all the blocks in the vector space + for (int blkIndex = 0; blkIndex < prodVS->numBlocks(); blkIndex++) + { + int subDim = 0; + + // construct the block vector + blockThyraToEpetra( + numVectors, localData, leadingDim, prodX->getMultiVectorBlock(blkIndex), subDim); + + // shift to the next block + localData += subDim; + + // account for the size of this subblock + localDim += subDim; + } +} + +Teuchos::RCP> Core::LinearSolver::Utils::get_linalg_multi_vector_from_thyra( + const Core::LinAlg::Map& map, const Teuchos::RCP>& thyra_vector) +{ + // build an Epetra_MultiVector object + int numVectors = thyra_vector->domain()->dim(); + + auto epetra_multi_vector = Epetra_MultiVector(map.get_epetra_map(), numVectors); + + // extract local information from the Epetra_MultiVector + int leadingDim = 0, localDim = 0; + double* epetraData = 0; + epetra_multi_vector.ExtractView(&epetraData, &leadingDim); + + // perform recursive copy + blockThyraToEpetra(numVectors, epetraData, leadingDim, thyra_vector, localDim); + + return Teuchos::make_rcp>(epetra_multi_vector); +} + FOUR_C_NAMESPACE_CLOSE diff --git a/src/core/linear_solver/src/utils/4C_linear_solver_thyra_utils.hpp b/src/core/linear_solver/src/utils/4C_linear_solver_thyra_utils.hpp index 315cfce74e4..3f487620c69 100644 --- a/src/core/linear_solver/src/utils/4C_linear_solver_thyra_utils.hpp +++ b/src/core/linear_solver/src/utils/4C_linear_solver_thyra_utils.hpp @@ -21,10 +21,13 @@ FOUR_C_NAMESPACE_OPEN namespace Core::LinearSolver::Utils { + // Create Thyra vector space from Map Teuchos::RCP> create_thyra_map(const LinAlg::Map& map); - Teuchos::RCP> create_thyra_multi_vector( - const LinAlg::MultiVector& multi_vector, const LinAlg::Map& map); + + // Create Thyra linear operator from SparseOperator, SparseMatrix or BlockSparseMatrixBase + Teuchos::RCP> create_thyra_linear_op( + LinAlg::SparseOperator& matrix, Core::LinAlg::DataAccess access); Teuchos::RCP> create_thyra_linear_op( const LinAlg::SparseMatrix& matrix, LinAlg::DataAccess access); @@ -32,11 +35,28 @@ namespace Core::LinearSolver::Utils Teuchos::RCP> create_thyra_linear_op( const LinAlg::BlockSparseMatrixBase& matrix, LinAlg::DataAccess access); + + // Create Thyra multi vector from MultiVector and Map or Thyra vector space + Teuchos::RCP> create_thyra_multi_vector( + const LinAlg::MultiVector& multi_vector, const LinAlg::Map& map); + + Teuchos::RCP> create_thyra_multi_vector( + const LinAlg::MultiVector& multi_vector, + Teuchos::RCP> map); + + +// Get Epetra vector from Thyra Teuchos::RCP get_epetra_vector_from_thyra( const Core::LinAlg::Map& map, const Teuchos::RCP<::Thyra::VectorBase>& thyra_vector); Teuchos::RCP get_epetra_vector_from_thyra(const Core::LinAlg::Map& map, const Teuchos::RCP>& thyra_vector); + + + // Get MultiVector from Thyra + Teuchos::RCP> get_linalg_multi_vector_from_thyra(const Core::LinAlg::Map& map, + const Teuchos::RCP>& thyra_vector); + } // namespace Core::LinearSolver::Utils FOUR_C_NAMESPACE_CLOSE From 02ebd0bf0105a340e2fc8764883dd132b6cc9d9f Mon Sep 17 00:00:00 2001 From: Max Firmbach Date: Thu, 4 Dec 2025 20:24:34 +0100 Subject: [PATCH 2/4] Add projection preconditioner based on thyra --- .../sparse/4C_linalg_projected_precond.cpp | 58 ----- .../sparse/4C_linalg_projected_precond.hpp | 201 ------------------ ...inear_solver_preconditioner_projection.cpp | 72 ++++++- ...inear_solver_preconditioner_projection.hpp | 140 +++++++++++- 4 files changed, 207 insertions(+), 264 deletions(-) delete mode 100644 src/core/linalg/src/sparse/4C_linalg_projected_precond.cpp delete mode 100644 src/core/linalg/src/sparse/4C_linalg_projected_precond.hpp diff --git a/src/core/linalg/src/sparse/4C_linalg_projected_precond.cpp b/src/core/linalg/src/sparse/4C_linalg_projected_precond.cpp deleted file mode 100644 index 57d64ee939f..00000000000 --- a/src/core/linalg/src/sparse/4C_linalg_projected_precond.cpp +++ /dev/null @@ -1,58 +0,0 @@ -// This file is part of 4C multiphysics licensed under the -// GNU Lesser General Public License v3.0 or later. -// -// See the LICENSE.md file in the top-level for license information. -// -// SPDX-License-Identifier: LGPL-3.0-or-later - -#include "4C_linalg_projected_precond.hpp" - -#include "4C_linalg_krylov_projector.hpp" -#include "4C_linalg_vector.hpp" -#include "4C_linear_solver_method_projector.hpp" -#include "4C_utils_exceptions.hpp" - - -FOUR_C_NAMESPACE_OPEN - -/* -------------------------------------------------------------------- - Constructor - -------------------------------------------------------------------- */ -Core::LinAlg::LinalgPrecondOperator::LinalgPrecondOperator(std::shared_ptr precond, - bool project, std::shared_ptr projector) - : project_(project), precond_(precond), projector_(projector) -{ - if (project_ && (projector == nullptr)) - FOUR_C_THROW("Kernel projection enabled but got no projector object"); - - return; -} // Core::LinAlg::LinalgPrecondOperator::LinalgPrecondOperator - - -/* -------------------------------------------------------------------- - (Modified) ApplyInverse call - -------------------------------------------------------------------- */ -int Core::LinAlg::LinalgPrecondOperator::ApplyInverse( - const Epetra_MultiVector& X, Epetra_MultiVector& Y) const -{ - int ierr = 0; - // Apply the inverse preconditioner to get new basis vector for the - // Krylov space - ierr = precond_->ApplyInverse(X, Y); - - // if necessary, project out matrix kernel to maintain well-posedness - // of problem - if (project_) - { - View Y_view(Y); - - FOUR_C_ASSERT_ALWAYS(Y.NumVectors() == 1, - "Expecting only one solution vector during projector call! Got {} vectors.", - Y.NumVectors()); - Y_view.underlying().get_vector(0) = projector_->to_full(Y_view.underlying().get_vector(0)); - } - - return (ierr); -} // Core::LinAlg::LinalgPrecondOperator::ApplyInverse - -FOUR_C_NAMESPACE_CLOSE diff --git a/src/core/linalg/src/sparse/4C_linalg_projected_precond.hpp b/src/core/linalg/src/sparse/4C_linalg_projected_precond.hpp deleted file mode 100644 index 12cc917406a..00000000000 --- a/src/core/linalg/src/sparse/4C_linalg_projected_precond.hpp +++ /dev/null @@ -1,201 +0,0 @@ -// This file is part of 4C multiphysics licensed under the -// GNU Lesser General Public License v3.0 or later. -// -// See the LICENSE.md file in the top-level for license information. -// -// SPDX-License-Identifier: LGPL-3.0-or-later - -#ifndef FOUR_C_LINALG_PROJECTED_PRECOND_HPP -#define FOUR_C_LINALG_PROJECTED_PRECOND_HPP - - -#include "4C_config.hpp" - -#include "4C_comm_mpi_utils.hpp" -#include "4C_linear_solver_method_projector.hpp" - -#include - -#include - -// forward declarations - -FOUR_C_NAMESPACE_OPEN - -namespace Core::LinAlg -{ - class KrylovProjector; - - - /*! - - A common interface for ifpack, ml and simpler preconditioners. - This interface allows a modification of the vector returned - by the ApplyInverse call, which is necessary to do a solution on - a Krylov space krylovized to certain (for example rigid body) - modes. - - The linalg preconditioner interface class holds a pointer (rcp) - to the actual preconditioner. All methods implemented to support - the Epetra_Operator interface just call the corresponding functions - of the actual preconditioner. - - Only the ApplyInverse method is modified and performs the - projection if desired. - - See linalg_projected_operator.H for related docu and code. - - */ - class LinalgPrecondOperator : public Epetra_Operator - { - public: - /*! - \brief Standard Constructor - - */ - LinalgPrecondOperator(std::shared_ptr precond, bool project, - std::shared_ptr projector); - - - - //! @name Attribute set methods required to support the Epetra_Operator interface - - int SetUseTranspose(bool UseTranspose) override - { - return (precond_->SetUseTranspose(UseTranspose)); - } - //! @} - - //! @name Mathematical functions required to support the Epetra_Operator interface (pass - //! through) - int Apply(const Epetra_MultiVector& X, Epetra_MultiVector& Y) const override - { - return (precond_->Apply(X, Y)); - } - - double NormInf() const override { return (precond_->NormInf()); } - //! @} - - //! @name Attribute access functions required to support the Epetra_Operator interface - const char* Label() const override { return (precond_->Label()); } - - bool UseTranspose() const override { return (precond_->UseTranspose()); } - - bool HasNormInf() const override { return (precond_->HasNormInf()); } - - const Epetra_Comm& Comm() const override { return (precond_->Comm()); } - - const Epetra_Map& OperatorDomainMap() const override { return (precond_->OperatorDomainMap()); } - - const Epetra_Map& OperatorRangeMap() const override { return (precond_->OperatorRangeMap()); } - - //! @} - - //! @name Mathematical functions required to support the Epetra_Operator interface (modified) - /* - (Modified) ApplyInverse call - - This method calls ApplyInverse on the actual preconditioner and, the - solution is krylovized against a set of weight vectors provided in a - multivector. - - This is done using a projector P defined by - - T - x * w - P x = x - ------ c - T - w * c - - w is the vector of weights, c a vector of ones (in the dofs under - consideration) corresponding to the matrix kernel. - - The system we are solving with this procedure is not Au=b for u (since A - might be singular), but we are solving - - / T \ T - | P A | P u = P b , - \ / - - for the projection of the solution Pu, i.e. in the preconditioned case - - - -+ - / T \ / -1 \ T | - | P * A | * | P * M | * xi = P * b | - \ / \ / | - -1 | - x = P * M * xi | - -+ - - - Hence, P is always associated with the apply inverse call of the - preconditioner (the right bracket) and always called after the call - to ApplyInverse. - - - Properties of P are: - - 1) c defines the kernel of P, i.e. P projects out the matrix kernel - - T - c * w - P c = c - ------- c = c - c = 0 - T - w * c - - 2) The space spanned by P x is krylov to the weight vector - - / T \ T - T / \ T | x * w | T x * w T T T - w * | P x | = w * | x - ------- c | = w * x - ------- * w * c = w * x - w * x = 0 - \ / | T | T - \ w * c / w * c - - - This modified Apply call is for singular matrices A when c is - a vector defining A's nullspace. The preceding projection - operator ensures - | | - -+- -+-T - A u = A u where u * c =0, - - even if A*c != 0 (for numerical inaccuracies during the computation - of A) - - See the following article for further reading: - - @article{1055401, - author = {Bochev,, Pavel and Lehoucq,, R. B.}, - title = {On the Finite Element Solution of the Pure Neumann Problem}, - journal = {SIAM Rev.}, - volume = {47}, - number = {1}, - year = {2005}, - issn = {0036-1445}, - pages = {50--66}, - doi = {http://dx.doi.org/10.1137/S0036144503426074}, - publisher = {Society for Industrial and Applied Mathematics}, - address = {Philadelphia, PA, USA}, - } - - */ - int ApplyInverse(const Epetra_MultiVector& X, Epetra_MultiVector& Y) const override; - - - private: - //! flag whether to do a projection or just pass through - bool project_; - - //! the actual preconditioner - std::shared_ptr precond_; - - //! projector - std::shared_ptr projector_; - }; - -} // namespace Core::LinAlg - -FOUR_C_NAMESPACE_CLOSE - -#endif diff --git a/src/core/linear_solver/src/preconditioner/4C_linear_solver_preconditioner_projection.cpp b/src/core/linear_solver/src/preconditioner/4C_linear_solver_preconditioner_projection.cpp index 7bc6425bd57..a6894789b4f 100644 --- a/src/core/linear_solver/src/preconditioner/4C_linear_solver_preconditioner_projection.cpp +++ b/src/core/linear_solver/src/preconditioner/4C_linear_solver_preconditioner_projection.cpp @@ -7,14 +7,16 @@ #include "4C_linear_solver_preconditioner_projection.hpp" -#include "4C_linalg_projected_precond.hpp" #include "4C_linear_solver_method_projector.hpp" +#include + #include FOUR_C_NAMESPACE_OPEN - +//---------------------------------------------------------------------------------- +//---------------------------------------------------------------------------------- Core::LinearSolver::ProjectionPreconditioner::ProjectionPreconditioner( std::shared_ptr preconditioner, std::shared_ptr projector) @@ -38,8 +40,70 @@ void Core::LinearSolver::ProjectionPreconditioner::setup( // actual preconditioner is called first and the projection is done // afterwards. - p_ = std::make_shared( - preconditioner_->prec_operator(), true, projector_); + p_thyra_ = Teuchos::make_rcp( + preconditioner_->thyra_operator(), true, projector_); +} + + +//---------------------------------------------------------------------------------- +//---------------------------------------------------------------------------------- +Core::LinearSolver::LinalgPrecondOperator::LinalgPrecondOperator( + Teuchos::RCP> precond, bool project, + std::shared_ptr projector) + : project_(project), precond_(precond), projector_(projector) +{ + if (project_ && (projector == nullptr)) + FOUR_C_THROW("Kernel projection enabled but got no projector object"); +} + +//---------------------------------------------------------------------------------- +//---------------------------------------------------------------------------------- +Teuchos::RCP> +Core::LinearSolver::LinalgPrecondOperator::range() const +{ + return precond_->range(); +} + +//---------------------------------------------------------------------------------- +//---------------------------------------------------------------------------------- +Teuchos::RCP> +Core::LinearSolver::LinalgPrecondOperator::domain() const +{ + return precond_->domain(); +} + +//---------------------------------------------------------------------------------- +//---------------------------------------------------------------------------------- +bool Core::LinearSolver::LinalgPrecondOperator::opSupportedImpl( + const Thyra::EOpTransp M_trans) const +{ + return (M_trans == Thyra::NOTRANS); +} + +//---------------------------------------------------------------------------------- +//---------------------------------------------------------------------------------- +void Core::LinearSolver::LinalgPrecondOperator::applyImpl(const Thyra::EOpTransp M_trans, + const Thyra::MultiVectorBase& x, const Teuchos::Ptr>& y, + const double alpha, const double beta) const +{ + // Apply the inverse preconditioner to get new basis vector for the + // Krylov space + precond_->apply(::Thyra::NOTRANS, x, y, alpha, beta); + + // if necessary, project out matrix kernel to maintain well-posedness + // of problem + if (project_) + { + auto map = Thyra::get_Epetra_Map(range()); + auto Y = Thyra::get_Epetra_MultiVector(Teuchos::rcpFromPtr(y), map); + Core::LinAlg::View Y_view(*Y); + + FOUR_C_ASSERT_ALWAYS(Y->NumVectors() == 1, + "Expecting only one solution vector during projector call! Got {} vectors.", + Y->NumVectors()); + + Y_view.underlying().get_vector(0) = projector_->to_full(Y_view.underlying().get_vector(0)); + } } FOUR_C_NAMESPACE_CLOSE diff --git a/src/core/linear_solver/src/preconditioner/4C_linear_solver_preconditioner_projection.hpp b/src/core/linear_solver/src/preconditioner/4C_linear_solver_preconditioner_projection.hpp index 7dadb190eea..23ef2c85f3c 100644 --- a/src/core/linear_solver/src/preconditioner/4C_linear_solver_preconditioner_projection.hpp +++ b/src/core/linear_solver/src/preconditioner/4C_linear_solver_preconditioner_projection.hpp @@ -36,7 +36,7 @@ namespace Core::LinearSolver Teuchos::RCP> thyra_operator() const override { - FOUR_C_THROW("One should not end up here."); + return p_thyra_; } private: @@ -49,6 +49,144 @@ namespace Core::LinearSolver Teuchos::RCP> p_thyra_; }; + + /*! + A common interface for ifpack, ml and simpler preconditioners. + This interface allows a modification of the vector returned + by the ApplyInverse call, which is necessary to do a solution on + a Krylov space krylovized to certain (for example rigid body) + modes. + + The linalg preconditioner interface class holds a pointer (rcp) + to the actual preconditioner. All methods implemented to support + the Epetra_Operator interface just call the corresponding functions + of the actual preconditioner. + + Only the ApplyInverse method is modified and performs the + projection if desired. + + See linalg_projected_operator.H for related docu and code. + + */ + class LinalgPrecondOperator : public Thyra::LinearOpBase + { + public: + LinalgPrecondOperator(Teuchos::RCP> precond, bool project, + std::shared_ptr projector); + + /** @brief Range space of this operator */ + virtual Teuchos::RCP> range() const; + + /** @brief Domain space of this operator */ + virtual Teuchos::RCP> domain() const; + + virtual bool opSupportedImpl(const Thyra::EOpTransp M_trans) const; + + //! @} + + //! @name Mathematical functions required to support the Epetra_Operator interface (modified) + /* + (Modified) ApplyInverse call + + This method calls ApplyInverse on the actual preconditioner and, the + solution is krylovized against a set of weight vectors provided in a + multivector. + + This is done using a projector P defined by + + T + x * w + P x = x - ------ c + T + w * c + + w is the vector of weights, c a vector of ones (in the dofs under + consideration) corresponding to the matrix kernel. + + The system we are solving with this procedure is not Au=b for u (since A + might be singular), but we are solving + + / T \ T + | P A | P u = P b , + \ / + + for the projection of the solution Pu, i.e. in the preconditioned case + + + -+ + / T \ / -1 \ T | + | P * A | * | P * M | * xi = P * b | + \ / \ / | + -1 | + x = P * M * xi | + -+ + + + Hence, P is always associated with the apply inverse call of the + preconditioner (the right bracket) and always called after the call + to ApplyInverse. + + + Properties of P are: + + 1) c defines the kernel of P, i.e. P projects out the matrix kernel + + T + c * w + P c = c - ------- c = c - c = 0 + T + w * c + + 2) The space spanned by P x is krylov to the weight vector + + / T \ T + T / \ T | x * w | T x * w T T T + w * | P x | = w * | x - ------- c | = w * x - ------- * w * c = w * x - w * x = 0 + \ / | T | T + \ w * c / w * c + + + This modified Apply call is for singular matrices A when c is + a vector defining A's nullspace. The preceding projection + operator ensures + | | + -+- -+-T + A u = A u where u * c =0, + + even if A*c != 0 (for numerical inaccuracies during the computation + of A) + + See the following article for further reading: + + @article{1055401, + author = {Bochev,, Pavel and Lehoucq,, R. B.}, + title = {On the Finite Element Solution of the Pure Neumann Problem}, + journal = {SIAM Rev.}, + volume = {47}, + number = {1}, + year = {2005}, + issn = {0036-1445}, + pages = {50--66}, + doi = {http://dx.doi.org/10.1137/S0036144503426074}, + publisher = {Society for Industrial and Applied Mathematics}, + address = {Philadelphia, PA, USA}, + } + + */ + virtual void applyImpl(const Thyra::EOpTransp M_trans, const Thyra::MultiVectorBase& x, + const Teuchos::Ptr>& y, const double alpha, + const double beta) const; + + private: + //! flag whether to do a projection or just pass through + bool project_; + + //! the actual preconditioner + Teuchos::RCP> precond_; + + //! projector + std::shared_ptr projector_; + }; } // namespace Core::LinearSolver FOUR_C_NAMESPACE_CLOSE From 7a0d0a3eb10b87165e3917b4574a234f17d1c619 Mon Sep 17 00:00:00 2001 From: Max Firmbach Date: Thu, 4 Dec 2025 20:35:17 +0100 Subject: [PATCH 3/4] Post cleanup preconditioner interface --- .../amgnxn/4C_linear_solver_amgnxn_preconditioner.cpp | 4 ++-- .../amgnxn/4C_linear_solver_amgnxn_preconditioner.hpp | 4 +--- .../src/method/4C_linear_solver_method_iterative.cpp | 4 ++-- .../4C_linear_solver_preconditioner_ifpack.cpp | 4 +--- .../4C_linear_solver_preconditioner_ifpack.hpp | 11 ++--------- .../4C_linear_solver_preconditioner_muelu.cpp | 6 ++---- .../4C_linear_solver_preconditioner_muelu.hpp | 11 ++--------- .../4C_linear_solver_preconditioner_projection.cpp | 4 ++-- .../4C_linear_solver_preconditioner_projection.hpp | 11 ++--------- .../4C_linear_solver_preconditioner_teko.cpp | 4 +--- .../4C_linear_solver_preconditioner_teko.hpp | 11 ++--------- .../4C_linear_solver_preconditioner_type.hpp | 5 +---- 12 files changed, 20 insertions(+), 59 deletions(-) diff --git a/src/core/linear_solver/src/amgnxn/4C_linear_solver_amgnxn_preconditioner.cpp b/src/core/linear_solver/src/amgnxn/4C_linear_solver_amgnxn_preconditioner.cpp index b0662737315..6872e6b669c 100644 --- a/src/core/linear_solver/src/amgnxn/4C_linear_solver_amgnxn_preconditioner.cpp +++ b/src/core/linear_solver/src/amgnxn/4C_linear_solver_amgnxn_preconditioner.cpp @@ -34,12 +34,12 @@ Core::LinearSolver::AmGnxnPreconditioner::AmGnxnPreconditioner(Teuchos::Paramete /*------------------------------------------------------------------------------*/ /*------------------------------------------------------------------------------*/ - +/* std::shared_ptr Core::LinearSolver::AmGnxnPreconditioner::prec_operator() const { return p_; } - +*/ /*------------------------------------------------------------------------------*/ /*------------------------------------------------------------------------------*/ diff --git a/src/core/linear_solver/src/amgnxn/4C_linear_solver_amgnxn_preconditioner.hpp b/src/core/linear_solver/src/amgnxn/4C_linear_solver_amgnxn_preconditioner.hpp index 1069ff58ccb..e23e89cebb8 100644 --- a/src/core/linear_solver/src/amgnxn/4C_linear_solver_amgnxn_preconditioner.hpp +++ b/src/core/linear_solver/src/amgnxn/4C_linear_solver_amgnxn_preconditioner.hpp @@ -39,9 +39,7 @@ namespace Core::LinearSolver virtual void setup(std::shared_ptr A); /// linear operator used for preconditioning - std::shared_ptr prec_operator() const override; - - Teuchos::RCP> thyra_operator() const override + Teuchos::RCP> prec_operator() const override { FOUR_C_THROW("One should not end up here."); } diff --git a/src/core/linear_solver/src/method/4C_linear_solver_method_iterative.cpp b/src/core/linear_solver/src/method/4C_linear_solver_method_iterative.cpp index 641e652debc..9f6074a3c87 100644 --- a/src/core/linear_solver/src/method/4C_linear_solver_method_iterative.cpp +++ b/src/core/linear_solver/src/method/4C_linear_solver_method_iterative.cpp @@ -99,8 +99,8 @@ int Core::LinearSolver::IterativeSolver::solve(Core::LinAlg::MultiVector if (preconditioner_ != nullptr) { - // TODO: Preconditioner seems to work not that properly (Blocked MueLu and Krylov Projection) - auto belosPrec = preconditioner_->thyra_operator(); + // TODO: Preconditioner seems to work not that properly (Blocked MueLu) + auto belosPrec = preconditioner_->prec_operator(); problem->setRightPrec(belosPrec); } diff --git a/src/core/linear_solver/src/preconditioner/4C_linear_solver_preconditioner_ifpack.cpp b/src/core/linear_solver/src/preconditioner/4C_linear_solver_preconditioner_ifpack.cpp index bc79e7ebfa0..6545cfc54cd 100644 --- a/src/core/linear_solver/src/preconditioner/4C_linear_solver_preconditioner_ifpack.cpp +++ b/src/core/linear_solver/src/preconditioner/4C_linear_solver_preconditioner_ifpack.cpp @@ -71,9 +71,7 @@ void Core::LinearSolver::IFPACKPreconditioner::setup( builder.createPreconditioningStrategy("Ifpack"); Teuchos::RCP> prec = Thyra::prec(*precFactory, pmatrix_); - p_thyra_ = prec->getUnspecifiedPrecOp(); - - p_ = std::make_shared(p_thyra_); + p_ = prec->getUnspecifiedPrecOp(); } FOUR_C_NAMESPACE_CLOSE diff --git a/src/core/linear_solver/src/preconditioner/4C_linear_solver_preconditioner_ifpack.hpp b/src/core/linear_solver/src/preconditioner/4C_linear_solver_preconditioner_ifpack.hpp index ac7a11b481f..45ce88394ae 100644 --- a/src/core/linear_solver/src/preconditioner/4C_linear_solver_preconditioner_ifpack.hpp +++ b/src/core/linear_solver/src/preconditioner/4C_linear_solver_preconditioner_ifpack.hpp @@ -33,12 +33,7 @@ namespace Core::LinearSolver void setup(Core::LinAlg::SparseOperator& matrix, Core::LinAlg::MultiVector& b) override; /// linear operator used for preconditioning - std::shared_ptr prec_operator() const override { return p_; } - - Teuchos::RCP> thyra_operator() const override - { - return p_thyra_; - } + Teuchos::RCP> prec_operator() const override { return p_; } private: //! IFPACK parameter list @@ -47,10 +42,8 @@ namespace Core::LinearSolver //! system of equations used for preconditioning used by P_ only Teuchos::RCP> pmatrix_; - Teuchos::RCP> p_thyra_; - //! preconditioner - std::shared_ptr p_; + Teuchos::RCP> p_; }; } // namespace Core::LinearSolver diff --git a/src/core/linear_solver/src/preconditioner/4C_linear_solver_preconditioner_muelu.cpp b/src/core/linear_solver/src/preconditioner/4C_linear_solver_preconditioner_muelu.cpp index a213925202e..c56ea4b2c65 100644 --- a/src/core/linear_solver/src/preconditioner/4C_linear_solver_preconditioner_muelu.cpp +++ b/src/core/linear_solver/src/preconditioner/4C_linear_solver_preconditioner_muelu.cpp @@ -119,9 +119,7 @@ void Core::LinearSolver::MueLuPreconditioner::setup( builder.createPreconditioningStrategy("MueLu"); Teuchos::RCP> prec = Thyra::prec(*precFactory, pmatrix_); - p_thyra_ = prec->getUnspecifiedPrecOp(); - - p_ = std::make_shared(p_thyra_); + p_ = prec->getUnspecifiedPrecOp(); } else { @@ -222,7 +220,7 @@ void Core::LinearSolver::MueLuPreconditioner::setup( } mueLuFactory.SetupHierarchy(*H_); - p_ = std::make_shared(H_); + // p_ = std::make_shared(H_); } } diff --git a/src/core/linear_solver/src/preconditioner/4C_linear_solver_preconditioner_muelu.hpp b/src/core/linear_solver/src/preconditioner/4C_linear_solver_preconditioner_muelu.hpp index 3d965b90f8a..347378b6c12 100644 --- a/src/core/linear_solver/src/preconditioner/4C_linear_solver_preconditioner_muelu.hpp +++ b/src/core/linear_solver/src/preconditioner/4C_linear_solver_preconditioner_muelu.hpp @@ -59,12 +59,7 @@ namespace Core::LinearSolver void setup(Core::LinAlg::SparseOperator& matrix, Core::LinAlg::MultiVector& b) override; //! linear operator used for preconditioning - std::shared_ptr prec_operator() const final { return p_; } - - Teuchos::RCP> thyra_operator() const override - { - return p_thyra_; - } + Teuchos::RCP> prec_operator() const override { return p_; } private: //! system of equations used for preconditioning used by P_ only @@ -75,9 +70,7 @@ namespace Core::LinearSolver Teuchos::ParameterList& muelulist_; //! preconditioner - std::shared_ptr p_; - - Teuchos::RCP> p_thyra_; + Teuchos::RCP> p_; //! MueLu hierarchy Teuchos::RCP> H_; diff --git a/src/core/linear_solver/src/preconditioner/4C_linear_solver_preconditioner_projection.cpp b/src/core/linear_solver/src/preconditioner/4C_linear_solver_preconditioner_projection.cpp index a6894789b4f..e6e362f9614 100644 --- a/src/core/linear_solver/src/preconditioner/4C_linear_solver_preconditioner_projection.cpp +++ b/src/core/linear_solver/src/preconditioner/4C_linear_solver_preconditioner_projection.cpp @@ -40,8 +40,8 @@ void Core::LinearSolver::ProjectionPreconditioner::setup( // actual preconditioner is called first and the projection is done // afterwards. - p_thyra_ = Teuchos::make_rcp( - preconditioner_->thyra_operator(), true, projector_); + p_ = Teuchos::make_rcp( + preconditioner_->prec_operator(), true, projector_); } diff --git a/src/core/linear_solver/src/preconditioner/4C_linear_solver_preconditioner_projection.hpp b/src/core/linear_solver/src/preconditioner/4C_linear_solver_preconditioner_projection.hpp index 23ef2c85f3c..baaeea29d4d 100644 --- a/src/core/linear_solver/src/preconditioner/4C_linear_solver_preconditioner_projection.hpp +++ b/src/core/linear_solver/src/preconditioner/4C_linear_solver_preconditioner_projection.hpp @@ -32,12 +32,7 @@ namespace Core::LinearSolver void setup(Core::LinAlg::SparseOperator& matrix, Core::LinAlg::MultiVector& b) override; /// linear operator used for preconditioning - std::shared_ptr prec_operator() const override { return p_; } - - Teuchos::RCP> thyra_operator() const override - { - return p_thyra_; - } + Teuchos::RCP> prec_operator() const override { return p_; } private: std::shared_ptr preconditioner_; @@ -45,9 +40,7 @@ namespace Core::LinearSolver /// projector object that does the actual work std::shared_ptr projector_; - std::shared_ptr p_; - - Teuchos::RCP> p_thyra_; + Teuchos::RCP> p_; }; /*! diff --git a/src/core/linear_solver/src/preconditioner/4C_linear_solver_preconditioner_teko.cpp b/src/core/linear_solver/src/preconditioner/4C_linear_solver_preconditioner_teko.cpp index cc6d95ab920..77d0b06c84b 100644 --- a/src/core/linear_solver/src/preconditioner/4C_linear_solver_preconditioner_teko.cpp +++ b/src/core/linear_solver/src/preconditioner/4C_linear_solver_preconditioner_teko.cpp @@ -147,9 +147,7 @@ void Core::LinearSolver::TekoPreconditioner::setup( builder.createPreconditioningStrategy("Teko"); Teuchos::RCP> prec = Thyra::prec(*precFactory, pmatrix_); - p_thyra_ = prec->getUnspecifiedPrecOp(); - - p_ = std::make_shared(p_thyra_); + p_ = prec->getUnspecifiedPrecOp(); } //---------------------------------------------------------------------------------- diff --git a/src/core/linear_solver/src/preconditioner/4C_linear_solver_preconditioner_teko.hpp b/src/core/linear_solver/src/preconditioner/4C_linear_solver_preconditioner_teko.hpp index 3dcf68c32b4..baec99297c1 100644 --- a/src/core/linear_solver/src/preconditioner/4C_linear_solver_preconditioner_teko.hpp +++ b/src/core/linear_solver/src/preconditioner/4C_linear_solver_preconditioner_teko.hpp @@ -36,12 +36,7 @@ namespace Core::LinearSolver void setup(Core::LinAlg::SparseOperator& matrix, Core::LinAlg::MultiVector& b) override; /// linear operator used for preconditioning - std::shared_ptr prec_operator() const override { return p_; } - - Teuchos::RCP> thyra_operator() const override - { - return p_thyra_; - } + Teuchos::RCP> prec_operator() const override { return p_; } private: Teuchos::ParameterList& tekolist_; @@ -50,9 +45,7 @@ namespace Core::LinearSolver Teuchos::RCP> pmatrix_; //! preconditioner - std::shared_ptr p_; - - Teuchos::RCP> p_thyra_; + Teuchos::RCP> p_; }; diff --git a/src/core/linear_solver/src/preconditioner/4C_linear_solver_preconditioner_type.hpp b/src/core/linear_solver/src/preconditioner/4C_linear_solver_preconditioner_type.hpp index 55ccbf4a51c..53a1cd9fe4d 100644 --- a/src/core/linear_solver/src/preconditioner/4C_linear_solver_preconditioner_type.hpp +++ b/src/core/linear_solver/src/preconditioner/4C_linear_solver_preconditioner_type.hpp @@ -44,10 +44,7 @@ namespace Core::LinearSolver virtual void setup( Core::LinAlg::SparseOperator& matrix, Core::LinAlg::MultiVector& b) = 0; - virtual Teuchos::RCP> thyra_operator() const = 0; - - /// linear operator used for preconditioning - virtual std::shared_ptr prec_operator() const = 0; + virtual Teuchos::RCP> prec_operator() const = 0; }; } // namespace Core::LinearSolver From dbc0bd1e6069a8dad4e360e641e31d2c87ab2c6a Mon Sep 17 00:00:00 2001 From: Max Firmbach Date: Sat, 6 Dec 2025 12:38:22 +0100 Subject: [PATCH 4/4] Multi status test with Belos --- .../input_files/fsi_pw_mono_fs_ga_ga.4C.yaml | 2 +- .../iterative_gmres_multiphysics.xml | 28 +++++++++++++++++++ 2 files changed, 29 insertions(+), 1 deletion(-) create mode 100644 tests/input_files/xml/linear_solver/iterative_gmres_multiphysics.xml diff --git a/tests/input_files/fsi_pw_mono_fs_ga_ga.4C.yaml b/tests/input_files/fsi_pw_mono_fs_ga_ga.4C.yaml index a91ba08c29c..38b88623fc0 100644 --- a/tests/input_files/fsi_pw_mono_fs_ga_ga.4C.yaml +++ b/tests/input_files/fsi_pw_mono_fs_ga_ga.4C.yaml @@ -67,7 +67,7 @@ SOLVER 2: SOLVER: "Belos" AZPREC: "Teko" AZREUSE: 10 - SOLVER_XML_FILE: "xml/linear_solver/iterative_gmres_template.xml" + SOLVER_XML_FILE: "xml/linear_solver/iterative_gmres_multiphysics.xml" TEKO_XML_FILE: "xml/block_preconditioner/fluid_solid_ale.xml" NAME: "Fsi_Solver" DESIGN SURF DIRICH CONDITIONS: diff --git a/tests/input_files/xml/linear_solver/iterative_gmres_multiphysics.xml b/tests/input_files/xml/linear_solver/iterative_gmres_multiphysics.xml new file mode 100644 index 00000000000..20d47991439 --- /dev/null +++ b/tests/input_files/xml/linear_solver/iterative_gmres_multiphysics.xml @@ -0,0 +1,28 @@ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + \ No newline at end of file