Skip to content

Gmg preconditioner #15

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Merged
tjhei merged 4 commits into from
Oct 24, 2025
+139 −132
Merged

Gmg preconditioner #15

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
4 commits
Select commit Hold shift + click to select a range
0a35ab6
template types ainverseoperator
quangx Sep 12, 2025
dd2cc38
fixup
quangx Oct 17, 2025
c0cd4a6
Merge branch 'main' into gmg_preconditioner
quangx Oct 20, 2025
c49b67a
fix typo
quangx Oct 21, 2025
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
128 changes: 128 additions & 0 deletions include/aspect/simulator/solver/block_stokes_preconditioner.h
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -20,11 +20,139 @@
#ifndef aspect_block_stokes_preconditioner_h
#define aspect_block_stokes_preconditioner_h

#include <deal.II/lac/solver_bicgstab.h>

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 PreconditionerA, class VectorType, class ABlockType>
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 <class PreconditionerA,class VectorType, class ABlockType>
InverseVelocityBlock<PreconditionerA,VectorType,ABlockType>::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 <class PreconditionerA, class VectorType, class ABlockType>
void InverseVelocityBlock<PreconditionerA,VectorType,ABlockType>::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<VectorType> mem;

try
{
dst = 0.0;

if (A_block_is_symmetric)
{
SolverCG<VectorType> 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<VectorType>
solver(solver_control,
mem,
typename SolverBicgstab<VectorType>::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<SolverControl> {solver_control},
exc,
src.get_mpi_communicator());
}
}
else
{
preconditioner.vmult (dst, src);
n_iterations_ += 1;
}
}



template <class PreconditionerA, class VectorType, class ABlockType>
unsigned int InverseVelocityBlock<PreconditionerA, VectorType, ABlockType>::n_iterations() const
{
return n_iterations_;
}

/**
* Implement the block Schur preconditioner
* (A B^T; 0 S)^{-1}.
Expand Down
138 changes: 7 additions & 131 deletions source/simulator/solver.cc
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -32,6 +32,7 @@
#include <deal.II/lac/solver_bicgstab.h>
#include <deal.II/lac/solver_cg.h>
#include <deal.II/fe/fe_values.h>
#include <deal.II/lac/trilinos_vector.h>

namespace aspect
{
Expand Down Expand Up @@ -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 PreconditionerA>
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 <class PreconditionerA>
InverseVelocityBlock<PreconditionerA>::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 <class PreconditionerA>
void InverseVelocityBlock<PreconditionerA>::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<LinearAlgebra::Vector> mem;

try
{
dst = 0.0;

if (A_block_is_symmetric)
{
SolverCG<LinearAlgebra::Vector> 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<LinearAlgebra::Vector>
solver(solver_control,
mem,
SolverBicgstab<LinearAlgebra::Vector>::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<SolverControl> {solver_control},
exc,
src.get_mpi_communicator());
}
}
else
{
preconditioner.vmult (dst, src);
n_iterations_ += 1;
}
}



template <class PreconditionerA>
unsigned int InverseVelocityBlock<PreconditionerA>::n_iterations() const
{
return n_iterations_;
}

/**
* Base class for Schur Complement operators.
*/
Expand Down Expand Up @@ -914,28 +790,28 @@ namespace aspect
}

// create a cheap preconditioner that consists of only a single V-cycle
internal::InverseVelocityBlock<LinearAlgebra::PreconditionAMG> inverse_velocity_block_cheap(
internal::InverseVelocityBlock<LinearAlgebra::PreconditionAMG,TrilinosWrappers::MPI::Vector,TrilinosWrappers::SparseMatrix> 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::InverseVelocityBlock<LinearAlgebra::PreconditionAMG>,
internal::SchurComplementOperator, LinearAlgebra::SparseMatrix, aspect::LinearAlgebra::BlockVector>
const internal::BlockSchurPreconditioner<internal::InverseVelocityBlock<LinearAlgebra::PreconditionAMG,TrilinosWrappers::MPI::Vector,TrilinosWrappers::SparseMatrix>,
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<LinearAlgebra::PreconditionAMG> inverse_velocity_block_expensive(
internal::InverseVelocityBlock<LinearAlgebra::PreconditionAMG,TrilinosWrappers::MPI::Vector,dealii::TrilinosWrappers::SparseMatrix> 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::InverseVelocityBlock<LinearAlgebra::PreconditionAMG>,
internal::SchurComplementOperator, LinearAlgebra::SparseMatrix, aspect::LinearAlgebra::BlockVector>
const internal::BlockSchurPreconditioner<internal::InverseVelocityBlock<LinearAlgebra::PreconditionAMG,TrilinosWrappers::MPI::Vector,TrilinosWrappers::SparseMatrix>,
internal::SchurComplementOperator, LinearAlgebra::SparseMatrix, dealii::TrilinosWrappers::MPI::BlockVector>
preconditioner_expensive (
inverse_velocity_block_expensive,
*schur,
Expand Down
1 change: 1 addition & 0 deletions source/simulator/solver/stokes_matrix_free_local_smoothing.cc
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -21,6 +21,7 @@
#include <aspect/simulator/solver/matrix_free_operators.h>
#include <aspect/simulator/solver/stokes_matrix_free.h>
#include <aspect/simulator/solver/stokes_matrix_free_local_smoothing.h>
#include <aspect/simulator/solver/block_stokes_preconditioner.h>
#include <aspect/mesh_deformation/interface.h>

#include <aspect/mesh_deformation/interface.h>
Expand Down
4 changes: 3 additions & 1 deletion tests/stokes_solver_fail/screen-output
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -30,7 +30,9 @@ Additional information:
void dealii::SolverFGMRES<VectorType>::solve(const MatrixType&,
VectorType&, const VectorType&, const PreconditionerType&) [with
MatrixType = aspect::internal::StokesBlock; PreconditionerType =
aspect::internal::BlockSchurPreconditioner<aspect::internal::InverseVelocityBlock<dealii::TrilinosWrappers::PreconditionAMG>,
aspect::internal::BlockSchurPreconditioner<aspect::internal::InverseVelocityBlock<dealii::TrilinosWrappers::PreconditionAMG,
dealii::TrilinosWrappers::MPI::Vector,
dealii::TrilinosWrappers::SparseMatrix>,
aspect::internal::SchurComplementOperator,
dealii::TrilinosWrappers::SparseMatrix,
dealii::TrilinosWrappers::MPI::BlockVector>; VectorType =
Expand Down