Skip to content

Create SchurApproximation class so BlockSchurPreconditioner can be used for both AMG and GMG.#6806

Open
quangx wants to merge 2 commits intogeodynamics:mainfrom
quangx:gmg_preconditioner_2
Open

Create SchurApproximation class so BlockSchurPreconditioner can be used for both AMG and GMG.#6806
quangx wants to merge 2 commits intogeodynamics:mainfrom
quangx:gmg_preconditioner_2

Conversation

@quangx
Copy link
Contributor

@quangx quangx commented Dec 9, 2025

Currently ASPECT uses separate classes, BlockSchurPreconditioner and BlockSchurGMGPreconditioner, for AMG and GMG respectively. This pull request refactors the code so only BlockSchurPreconditioner is used in both cases.

@tjhei
Copy link
Member

tjhei commented Dec 9, 2025

Sadly, this requires c++17:
../include/aspect/simulator/solver/block_stokes_preconditioner.h:242:12: error: 'if constexpr' only available with '-std=c++17' or '-std=gnu++17' [-Werror]

@quangx quangx force-pushed the gmg_preconditioner_2 branch 2 times, most recently from cf29d33 to ca6b8b9 Compare December 10, 2025 23:19
const SchurComplementMatrixType &Schur_complement_block,
const bool do_solve_Schur_complement,
const double Schur_complement_tolerance
):
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

formatting

Copy link
Member

@gassmoeller gassmoeller left a comment

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I like the intention, but I dont quite understand. If this PR unifies the use of BlockSchurPreconditioner and BlockSchurGMGPreconditioner, shouldnt the second be deleted? But I dont see this as a change in this PR. Or is it still used in the GMG-GC solver?

Some comments otherwise, but just about code style.

Comment on lines 87 to 88


Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

This is not a strict rule, but we usually use 3 empty lines between function implementations, even if they are in header files. So please keep these empty lines.

Comment on lines +1215 to +1221
internal::InverseVelocityBlock<GMGPreconditioner,VectorType,ABlockMatrixType> inverse_velocity_block_cheap(
A_block_matrix,
prec_A,
/* do_solve_A = */ false,
sim.stokes_A_block_is_symmetric(),
this->get_parameters().linear_solver_A_block_tolerance);

Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I think these lines need indentation.

Comment on lines 1230 to 1234
internal::SchurApproximation<GMGPreconditioner, StokesMatrixType, SchurComplementMatrixType, VectorType> schur_approximation_cheap(prec_Schur,stokes_matrix,
Schur_complement_block_matrix, /*do_solve_Schur*/ false, this->get_parameters().linear_solver_S_block_tolerance);

internal::SchurApproximation<GMGPreconditioner, StokesMatrixType, SchurComplementMatrixType, VectorType> schur_approximation_expensive(prec_Schur,stokes_matrix,
Schur_complement_block_matrix, /*do_solve_Schur*/ true, this->get_parameters().linear_solver_S_block_tolerance);
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

break these lines similar to the ones above to make them easier to read

@quangx
Copy link
Contributor Author

quangx commented Jan 11, 2026

I like the intention, but I dont quite understand. If this PR unifies the use of BlockSchurPreconditioner and BlockSchurGMGPreconditioner, shouldnt the second be deleted? But I dont see this as a change in this PR. Or is it still used in the GMG-GC solver?

Some comments otherwise, but just about code style.

I do intend to delete the BlockSchrGMGPreconditioner class. I am debugging this PR currently because tests are failing. I will address the other feedback you gave, thank you for looking at this!

// but the matrix-based version is a SparseMatrix::vmult() that requires passing the
// individual blocks.
if constexpr (std::is_same_v<VectorType, dealii::LinearAlgebra::distributed::BlockVector<double>>)
BT_operator.vmult(tmp, dst);
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

try to set tmp=0; before this vmult.

Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

If you need tmp=0 before the application of vmult, then you have a bug in the implementation of vmult: It represents an assignment to tmp, and whatever is in that vector before the call to vmult should not matter.

}
else
{
schur_preconditioner.vmult(dst,src);
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

or try dst=0 here

if (tmp.size() == 0)
{
tmp.reinit(src);
}
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

In BlockSchurGMGPreconditioner we do dst = 0.0; here

@tjhei
Copy link
Member

tjhei commented Jan 24, 2026

@quangx Wolfgang is of course right: we should try to avoid zeroing out memory as much as possible for performance reasons. We now know that zeroing out data is enough to make all tests pass in this PR. Now you can try to remove as many of them as possible.
Conceptually, we only need to zero out source vectors passed to any solver.solve() call. Everything else should work as is. Ideally, we remove the calls to zero out temporary vectors from the block preconditioner and only add them where necessary furthest down the call stack (for example in InverseSchurComplement and only if we actually do a solve).

@quangx quangx force-pushed the gmg_preconditioner_2 branch from edccb4a to cc117b3 Compare January 30, 2026 01:43
@quangx quangx force-pushed the gmg_preconditioner_2 branch from 4847120 to 9340e45 Compare February 12, 2026 21:42
@tjhei
Copy link
Member

tjhei commented Feb 13, 2026

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment

Projects

None yet

Development

Successfully merging this pull request may close these issues.

4 participants