Skip to content

[WIP] Implement weighted BFBT with different discretization for (BC^{-1}B^T). #6440

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

Open
wants to merge 4 commits into
base: main
Choose a base branch
from
Open

[WIP] Implement weighted BFBT with different discretization for (BC^{-1}B^T). #6440

Changes from 3 commits
Commits
Show all changes
4 commits
Select commit Hold shift + click to select a range
a4dd9fe
initial commit
quangx Jun 14, 2025
61c79f3
change type to const for diag_operator
quangx Jun 14, 2025
abc7524
fixup
quangx Jun 17, 2025
c0eaef7
fixup
quangx Jun 17, 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
45 changes: 43 additions & 2 deletions source/simulator/solver.cc
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -27,15 +27,19 @@
#include <aspect/mesh_deformation/interface.h>

#include <deal.II/base/signaling_nan.h>
#include <deal.II/lac/diagonal_matrix.h>
#include <deal.II/lac/linear_operator.h>
#include <deal.II/lac/solver_gmres.h>
#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
{
namespace internal
{

/**
* Implement multiplication with Stokes part of system matrix. In essence, this
* object represents a 2x2 block matrix that corresponds to the top left
Expand Down Expand Up @@ -165,6 +169,8 @@ namespace aspect
}




/**
Copy link
Contributor

@bangerth bangerth Jun 17, 2025

Choose a reason for hiding this comment

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

We typically have 3 empty lines between top-level constructs. So take one off again here.

* Implement the block Schur preconditioner for the Stokes system.
*/
Expand Down Expand Up @@ -381,6 +387,33 @@ namespace aspect

};


/**
* Given a diagonal matrix stored as a vector,
* create an operator that represents its action.
*/

template<typename Range,
Copy link
Contributor

@bangerth bangerth Jun 17, 2025

Choose a reason for hiding this comment

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

Suggested change
/**
* Given a diagonal matrix stored as a vector,
* create an operator that represents its action.
*/
template<typename Range,
/**
* Given a diagonal matrix stored as a vector,
* create an operator that represents its action.
*/
template<typename Range,

typename Domain,
typename Payload>
LinearOperator<Range, Domain, Payload> diag_operator(const LinearOperator<Range,Domain,Payload> &exemplar, const TrilinosWrappers::MPI::Vector &diagonal)
{
LinearOperator<Range, Domain, Payload> return_op;


return_op.reinit_range_vector = exemplar.reinit_range_vector;
return_op.reinit_domain_vector = exemplar.reinit_domain_vector;

return_op.vmult = [&](Range &dest, const Domain &src)
{
dest = src;
dest.scale(diagonal);
};
return return_op;
}



/**
* This class approximates the Schur Complement inverse operator
* by S^{-1} = (BC^{-1}B^T)^{-1}(BC^{-1}AD^{-1}B^T)(BD^{-1}B^T)^{-1},
Expand Down Expand Up @@ -458,7 +491,15 @@ namespace aspect
SolverControl solver_control(5000, 1e-6 * src.l2_norm(), false, true);
SolverCG<TrilinosWrappers::MPI::Vector> solver(solver_control);
//Solve with Schur Complement approximation
solver.solve(mp_matrix,
const auto Op_A = LinearOperator<TrilinosWrappers::MPI::Vector>(system_matrix.block(0,0));
const auto Op_BT = LinearOperator<TrilinosWrappers::MPI::Vector>(system_matrix.block(0,1));
const auto Op_B = LinearOperator<TrilinosWrappers::MPI::Vector>(system_matrix.block(1,0));
const auto Op_C_inv = diag_operator(Op_A,inverse_lumped_mass_matrix);
const auto BC_invBT = Op_B*Op_C_inv*Op_BT;



solver.solve(BC_invBT,
ptmp,
src,
mp_preconditioner);
Expand All @@ -471,7 +512,7 @@ namespace aspect
system_matrix.block(1,0).vmult(ptmp,wtmp);

dst=0;
solver.solve(mp_matrix,
solver.solve(BC_invBT,
dst,
ptmp,
mp_preconditioner);
Expand Down
Loading