diff --git a/include/aspect/simulator/solver/block_stokes_preconditioner.h b/include/aspect/simulator/solver/block_stokes_preconditioner.h index a8900a80902..755c622e5e9 100644 --- a/include/aspect/simulator/solver/block_stokes_preconditioner.h +++ b/include/aspect/simulator/solver/block_stokes_preconditioner.h @@ -20,11 +20,139 @@ #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. This applies the action of A^{-1} by either + * performing a solve with A or using a preconditioner sweep. + */ + 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 a7bb7d8a460..dd41e24c22a 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 LinearAlgebra::SparseMatrix &matrix, - const PreconditionerA &preconditioner, - const bool do_solve_A, - const bool A_block_is_symmetric, - const double solver_tolerance); - - void vmult(LinearAlgebra::Vector &dst, - const LinearAlgebra::Vector &src) const; - - unsigned int n_iterations() const; - - private: - mutable unsigned int n_iterations_; - const LinearAlgebra::SparseMatrix &matrix; - const PreconditionerA &preconditioner; - const bool do_solve_A; - const bool A_block_is_symmetric; - const double solver_tolerance; - }; - - - - template - InverseVelocityBlock::InverseVelocityBlock( - const LinearAlgebra::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(LinearAlgebra::Vector &dst, - const LinearAlgebra::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,28 +790,28 @@ 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, - internal::SchurComplementOperator, LinearAlgebra::SparseMatrix, aspect::LinearAlgebra::BlockVector> + const internal::BlockSchurPreconditioner, + internal::SchurComplementOperator, LinearAlgebra::SparseMatrix, dealii::TrilinosWrappers::MPI::BlockVector> preconditioner_cheap ( inverse_velocity_block_cheap, *schur, 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, - internal::SchurComplementOperator, LinearAlgebra::SparseMatrix, aspect::LinearAlgebra::BlockVector> + const internal::BlockSchurPreconditioner, + internal::SchurComplementOperator, LinearAlgebra::SparseMatrix, dealii::TrilinosWrappers::MPI::BlockVector> preconditioner_expensive ( inverse_velocity_block_expensive, *schur, diff --git a/source/simulator/solver/stokes_matrix_free_local_smoothing.cc b/source/simulator/solver/stokes_matrix_free_local_smoothing.cc index cd4dc1fe9b2..0e63829d70c 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 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 =