From 0a35ab6f8cbfac5613769ad58bf6d9024adc70be Mon Sep 17 00:00:00 2001 From: Quang Hoang Date: Fri, 12 Sep 2025 15:06:45 -0400 Subject: [PATCH 1/3] template types ainverseoperator --- .../solver/block_stokes_preconditioner.h | 131 +++++++++++++++++ source/simulator/solver.cc | 134 +----------------- .../stokes_matrix_free_local_smoothing.cc | 3 +- tests/stokes_solver_fail/screen-output | 4 +- 4 files changed, 141 insertions(+), 131 deletions(-) diff --git a/include/aspect/simulator/solver/block_stokes_preconditioner.h b/include/aspect/simulator/solver/block_stokes_preconditioner.h index a8900a80902..68b33b5c1f9 100644 --- a/include/aspect/simulator/solver/block_stokes_preconditioner.h +++ b/include/aspect/simulator/solver/block_stokes_preconditioner.h @@ -20,11 +20,142 @@ #ifndef aspect_block_stokes_preconditioner_h #define aspect_block_stokes_preconditioner_h +#include + namespace aspect { namespace internal { + /** + * This class is used in the implementation of the right preconditioner + * as an approximation for the inverse of the velocity (A) block. + * This operator can either just apply the preconditioner (AMG) + * or perform an inner CG / BiCGStab solve with the same preconditioner. + */ + template + class InverseVelocityBlock + { + public: + /** + * Constructor. + * @param matrix The matrix that contains A (from the system matrix) + * @param preconditioner The preconditioner to be used + * @param do_solve_A A flag indicating whether we should actually solve with + * the matrix $A$, or only apply one preconditioner step with it. + * @param A_block_is_symmetric A flag indicating whether the matrix $A$ is symmetric. + * @param solver_tolerance The tolerance for the CG solver which computes + * the inverse of the A block. + */ + InverseVelocityBlock(const ABlockType &matrix, + const PreconditionerA &preconditioner, + const bool do_solve_A, + const bool A_block_is_symmetric, + const double solver_tolerance); + + void vmult(VectorType &dst, + const VectorType &src) const; + + unsigned int n_iterations() const; + + private: + mutable unsigned int n_iterations_; + const ABlockType &matrix; + const PreconditionerA &preconditioner; + const bool do_solve_A; + const bool A_block_is_symmetric; + const double solver_tolerance; + }; + + + + template + InverseVelocityBlock::InverseVelocityBlock( + const ABlockType &matrix, + const PreconditionerA &preconditioner, + const bool do_solve_A, + const bool A_block_is_symmetric, + const double solver_tolerance) + : n_iterations_ (0), + matrix (matrix), + preconditioner (preconditioner), + do_solve_A (do_solve_A), + A_block_is_symmetric (A_block_is_symmetric), + solver_tolerance (solver_tolerance) + {} + + + + /** + *Implements the vmult for InverseVelocityBlock. + */ + template + void InverseVelocityBlock::vmult(VectorType &dst, + const VectorType &src) const + { + + // Either solve with the top left block + // or just apply one preconditioner sweep (for the first few + // iterations of our two-stage outer GMRES iteration) + if (do_solve_A == true) + { + SolverControl solver_control(10000, src.l2_norm() * solver_tolerance); + PrimitiveVectorMemory mem; + + try + { + dst = 0.0; + + if (A_block_is_symmetric) + { + SolverCG solver(solver_control, mem); + solver.solve(matrix, dst, src, preconditioner); + } + else + { + // Use BiCGStab for non-symmetric matrices. + // BiCGStab can also solve indefinite systems if necessary. + // Do not compute the exact residual, as this + // is more expensive, and we only need an approximate solution. + SolverBicgstab + solver(solver_control, + mem, + typename SolverBicgstab::AdditionalData(/*exact_residual=*/ false)); + solver.solve(matrix, dst, src, preconditioner); + } + n_iterations_ += solver_control.last_step(); + } + catch (const std::exception &exc) + { + // if the solver fails, report the error from processor 0 with some additional + // information about its location, and throw a quiet exception on all other + // processors + Utilities::throw_linear_solver_failure_exception("iterative (top left) solver", + "BlockSchurPreconditioner::vmult", + std::vector {solver_control}, + exc, + src.get_mpi_communicator()); + } + } + else + { + preconditioner.vmult (dst, src); + n_iterations_ += 1; + } + } + + + + + + + + template + unsigned int InverseVelocityBlock::n_iterations() const + { + return n_iterations_; + } + /** * Implement the block Schur preconditioner * (A B^T; 0 S)^{-1}. diff --git a/source/simulator/solver.cc b/source/simulator/solver.cc index bb1a9629753..abc086265bc 100644 --- a/source/simulator/solver.cc +++ b/source/simulator/solver.cc @@ -32,6 +32,7 @@ #include #include #include +#include namespace aspect { @@ -165,131 +166,6 @@ namespace aspect return dst.l2_norm(); } - - - - - /** - * This class is used in the implementation of the right preconditioner - * as an approximation for the inverse of the velocity (A) block. - * This operator can either just apply the preconditioner (AMG) - * or perform an inner CG solve with the same preconditioner. - */ - template - class InverseVelocityBlock - { - public: - /** - * Constructor. - * @param matrix The matrix that contains A (from the system matrix) - * @param preconditioner The preconditioner to be used - * @param do_solve_A A flag indicating whether we should actually solve with - * the matrix $A$, or only apply one preconditioner step with it. - * @param A_block_is_symmetric A flag indicating whether the matrix $A$ is symmetric. - * @param A_block_tolerance The tolerance for the CG solver which computes - * the inverse of the A block. - */ - InverseVelocityBlock(const TrilinosWrappers::SparseMatrix &matrix, - const PreconditionerA &preconditioner, - const bool do_solve_A, - const bool A_block_is_symmetric, - const double solver_tolerance); - - void vmult(TrilinosWrappers::MPI::Vector &dst, - const TrilinosWrappers::MPI::Vector &src) const; - - unsigned int n_iterations() const; - - private: - mutable unsigned int n_iterations_; - const TrilinosWrappers::SparseMatrix &matrix; - const PreconditionerA &preconditioner; - const bool do_solve_A; - const bool A_block_is_symmetric; - const double solver_tolerance; - }; - - - - template - InverseVelocityBlock::InverseVelocityBlock( - const TrilinosWrappers::SparseMatrix &matrix, - const PreconditionerA &preconditioner, - const bool do_solve_A, - const bool A_block_is_symmetric, - const double solver_tolerance) - : n_iterations_ (0), - matrix (matrix), - preconditioner (preconditioner), - do_solve_A (do_solve_A), - A_block_is_symmetric (A_block_is_symmetric), - solver_tolerance (solver_tolerance) - {} - - - - template - void InverseVelocityBlock::vmult(TrilinosWrappers::MPI::Vector &dst, - const TrilinosWrappers::MPI::Vector &src) const - { - // Either solve with the top left block - // or just apply one preconditioner sweep (for the first few - // iterations of our two-stage outer GMRES iteration) - if (do_solve_A == true) - { - SolverControl solver_control(10000, src.l2_norm() * solver_tolerance); - PrimitiveVectorMemory mem; - - try - { - dst = 0.0; - - if (A_block_is_symmetric) - { - SolverCG solver(solver_control, mem); - solver.solve(matrix, dst, src, preconditioner); - } - else - { - // Use BiCGStab for non-symmetric matrices. - // BiCGStab can also solve indefinite systems if necessary. - // Do not compute the exact residual, as this - // is more expensive, and we only need an approximate solution. - SolverBicgstab - solver(solver_control, - mem, - SolverBicgstab::AdditionalData(/*exact_residual=*/ false)); - solver.solve(matrix, dst, src, preconditioner); - } - n_iterations_ += solver_control.last_step(); - } - catch (const std::exception &exc) - { - // if the solver fails, report the error from processor 0 with some additional - // information about its location, and throw a quiet exception on all other - // processors - Utilities::throw_linear_solver_failure_exception("iterative (top left) solver", - "BlockSchurPreconditioner::vmult", - std::vector {solver_control}, - exc, - src.get_mpi_communicator()); - } - } - else - { - preconditioner.vmult (dst, src); - n_iterations_ += 1; - } - } - - - - template - unsigned int InverseVelocityBlock::n_iterations() const - { - return n_iterations_; - } - /** * Base class for Schur Complement operators. */ @@ -914,13 +790,13 @@ namespace aspect } // create a cheap preconditioner that consists of only a single V-cycle - internal::InverseVelocityBlock inverse_velocity_block_cheap( + internal::InverseVelocityBlock inverse_velocity_block_cheap( system_matrix.block(velocity_block_index,velocity_block_index), *Amg_preconditioner, /* do_solve_A = */ false, stokes_A_block_is_symmetric(), parameters.linear_solver_A_block_tolerance); - const internal::BlockSchurPreconditioner, + const internal::BlockSchurPreconditioner, internal::SchurComplementOperator, LinearAlgebra::SparseMatrix, dealii::TrilinosWrappers::MPI::BlockVector> preconditioner_cheap ( inverse_velocity_block_cheap, @@ -928,13 +804,13 @@ namespace aspect system_matrix.block(0,1)); // create an expensive preconditioner that solves for the A block with CG - internal::InverseVelocityBlock inverse_velocity_block_expensive( + internal::InverseVelocityBlock inverse_velocity_block_expensive( system_matrix.block(velocity_block_index,velocity_block_index), *Amg_preconditioner, /* do_solve_A = */ true, stokes_A_block_is_symmetric(), parameters.linear_solver_A_block_tolerance); - const internal::BlockSchurPreconditioner, + const internal::BlockSchurPreconditioner, internal::SchurComplementOperator, LinearAlgebra::SparseMatrix, dealii::TrilinosWrappers::MPI::BlockVector> preconditioner_expensive ( inverse_velocity_block_expensive, diff --git a/source/simulator/solver/stokes_matrix_free_local_smoothing.cc b/source/simulator/solver/stokes_matrix_free_local_smoothing.cc index dc2b2829237..f06f3b0c5c9 100644 --- a/source/simulator/solver/stokes_matrix_free_local_smoothing.cc +++ b/source/simulator/solver/stokes_matrix_free_local_smoothing.cc @@ -21,6 +21,7 @@ #include #include #include +#include #include #include @@ -43,6 +44,7 @@ namespace aspect { + template void StokesMatrixFreeHandlerLocalSmoothingImplementation::declare_parameters(ParameterHandler &prm) @@ -1204,7 +1206,6 @@ namespace aspect solver_control_cheap.enable_history_data(); solver_control_expensive.enable_history_data(); - // create a cheap preconditioner that consists of only a single V-cycle const internal::BlockSchurGMGPreconditioner preconditioner_cheap (stokes_matrix, A_block_matrix, BT_block, Schur_complement_block_matrix, prec_A, prec_Schur, diff --git a/tests/stokes_solver_fail/screen-output b/tests/stokes_solver_fail/screen-output index 8a89c330b64..d672755068a 100644 --- a/tests/stokes_solver_fail/screen-output +++ b/tests/stokes_solver_fail/screen-output @@ -30,7 +30,9 @@ Additional information: void dealii::SolverFGMRES::solve(const MatrixType&, VectorType&, const VectorType&, const PreconditionerType&) [with MatrixType = aspect::internal::StokesBlock; PreconditionerType = - aspect::internal::BlockSchurPreconditioner, + aspect::internal::BlockSchurPreconditioner, aspect::internal::SchurComplementOperator, dealii::TrilinosWrappers::SparseMatrix, dealii::TrilinosWrappers::MPI::BlockVector>; VectorType = From dd2cc387f06cbc8da0a9b367acea90c874fae5e0 Mon Sep 17 00:00:00 2001 From: Quang Hoang Date: Fri, 17 Oct 2025 16:14:00 -0400 Subject: [PATCH 2/3] fixup --- .../aspect/simulator/solver/block_stokes_preconditioner.h | 7 ++----- .../simulator/solver/stokes_matrix_free_local_smoothing.cc | 2 +- 2 files changed, 3 insertions(+), 6 deletions(-) diff --git a/include/aspect/simulator/solver/block_stokes_preconditioner.h b/include/aspect/simulator/solver/block_stokes_preconditioner.h index 68b33b5c1f9..50fd6b4e6b4 100644 --- a/include/aspect/simulator/solver/block_stokes_preconditioner.h +++ b/include/aspect/simulator/solver/block_stokes_preconditioner.h @@ -87,7 +87,8 @@ namespace aspect /** - *Implements the vmult for InverseVelocityBlock. + * Implements the vmult for InverseVelocityBlock. This applies the action of A^{-1} by either + * performing a solve with A or using a preconditioner sweet. */ template void InverseVelocityBlock::vmult(VectorType &dst, @@ -146,10 +147,6 @@ namespace aspect - - - - template unsigned int InverseVelocityBlock::n_iterations() const { diff --git a/source/simulator/solver/stokes_matrix_free_local_smoothing.cc b/source/simulator/solver/stokes_matrix_free_local_smoothing.cc index f06f3b0c5c9..d66aa77643f 100644 --- a/source/simulator/solver/stokes_matrix_free_local_smoothing.cc +++ b/source/simulator/solver/stokes_matrix_free_local_smoothing.cc @@ -44,7 +44,6 @@ namespace aspect { - template void StokesMatrixFreeHandlerLocalSmoothingImplementation::declare_parameters(ParameterHandler &prm) @@ -1206,6 +1205,7 @@ namespace aspect solver_control_cheap.enable_history_data(); solver_control_expensive.enable_history_data(); + // create a cheap preconditioner that consists of only a single V-cycle const internal::BlockSchurGMGPreconditioner preconditioner_cheap (stokes_matrix, A_block_matrix, BT_block, Schur_complement_block_matrix, prec_A, prec_Schur, From c49b67a59c84a25df1bc73a2ec611f317c5a3b18 Mon Sep 17 00:00:00 2001 From: Quang Hoang Date: Mon, 20 Oct 2025 21:10:30 -0400 Subject: [PATCH 3/3] fix typo --- include/aspect/simulator/solver/block_stokes_preconditioner.h | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/include/aspect/simulator/solver/block_stokes_preconditioner.h b/include/aspect/simulator/solver/block_stokes_preconditioner.h index 50fd6b4e6b4..755c622e5e9 100644 --- a/include/aspect/simulator/solver/block_stokes_preconditioner.h +++ b/include/aspect/simulator/solver/block_stokes_preconditioner.h @@ -88,7 +88,7 @@ namespace aspect /** * Implements the vmult for InverseVelocityBlock. This applies the action of A^{-1} by either - * performing a solve with A or using a preconditioner sweet. + * performing a solve with A or using a preconditioner sweep. */ template void InverseVelocityBlock::vmult(VectorType &dst,