From a4dd9fe4aa62ba3afbe09831fe1d6077c010d753 Mon Sep 17 00:00:00 2001 From: Quang Hoang Date: Fri, 13 Jun 2025 18:12:08 -0600 Subject: [PATCH 1/4] initial commit --- source/simulator/solver.cc | 37 +++++++++++++++++++++++++++++++++++-- 1 file changed, 35 insertions(+), 2 deletions(-) diff --git a/source/simulator/solver.cc b/source/simulator/solver.cc index 03ca896b4aa..f0a58b391b4 100644 --- a/source/simulator/solver.cc +++ b/source/simulator/solver.cc @@ -27,15 +27,19 @@ #include #include +#include +#include #include #include #include #include +#include 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 @@ -163,6 +167,8 @@ namespace aspect return dst.l2_norm(); } + + /** @@ -380,6 +386,25 @@ namespace aspect virtual unsigned int n_iterations() const=0; }; + template + LinearOperator diag_operator(LinearOperator &exemplar,TrilinosWrappers::MPI::Vector &diagonal) + { + LinearOperator 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); + }; + // std::cout << "ok" << std::endl; + return return_op; + } /** * This class approximates the Schur Complement inverse operator @@ -458,7 +483,15 @@ namespace aspect SolverControl solver_control(5000, 1e-6 * src.l2_norm(), false, true); SolverCG solver(solver_control); //Solve with Schur Complement approximation - solver.solve(mp_matrix, + auto Op_A=LinearOperator(system_matrix.block(0,0)); + auto Op_BT = LinearOperator(system_matrix.block(0,1)); + auto Op_B = LinearOperator(system_matrix.block(1,0)); + auto Op_C = diag_operator(Op_A,inverse_lumped_mass_matrix); + auto matrix = Op_B*Op_C*Op_BT; + + + + solver.solve(matrix, ptmp, src, mp_preconditioner); @@ -471,7 +504,7 @@ namespace aspect system_matrix.block(1,0).vmult(ptmp,wtmp); dst=0; - solver.solve(mp_matrix, + solver.solve(matrix, dst, ptmp, mp_preconditioner); From 61c79f393125aaf0921c0ef85d040ca226f0908b Mon Sep 17 00:00:00 2001 From: Quang Hoang Date: Fri, 13 Jun 2025 20:29:05 -0600 Subject: [PATCH 2/4] change type to const for diag_operator --- source/simulator/solver.cc | 28 ++++++++++++++++++---------- 1 file changed, 18 insertions(+), 10 deletions(-) diff --git a/source/simulator/solver.cc b/source/simulator/solver.cc index f0a58b391b4..2d5a8a37b22 100644 --- a/source/simulator/solver.cc +++ b/source/simulator/solver.cc @@ -167,7 +167,7 @@ namespace aspect return dst.l2_norm(); } - + @@ -386,10 +386,17 @@ namespace aspect virtual unsigned int n_iterations() const=0; }; - template - LinearOperator diag_operator(LinearOperator &exemplar,TrilinosWrappers::MPI::Vector &diagonal) + LinearOperator diag_operator(LinearOperator &exemplar, const TrilinosWrappers::MPI::Vector &diagonal) { LinearOperator return_op; @@ -402,10 +409,11 @@ namespace aspect dest = src; dest.scale(diagonal); }; - // std::cout << "ok" << std::endl; 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}, @@ -483,11 +491,11 @@ namespace aspect SolverControl solver_control(5000, 1e-6 * src.l2_norm(), false, true); SolverCG solver(solver_control); //Solve with Schur Complement approximation - auto Op_A=LinearOperator(system_matrix.block(0,0)); - auto Op_BT = LinearOperator(system_matrix.block(0,1)); - auto Op_B = LinearOperator(system_matrix.block(1,0)); - auto Op_C = diag_operator(Op_A,inverse_lumped_mass_matrix); - auto matrix = Op_B*Op_C*Op_BT; + const auto Op_A = LinearOperator(system_matrix.block(0,0)); + const auto Op_BT = LinearOperator(system_matrix.block(0,1)); + const auto Op_B = LinearOperator(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; @@ -504,7 +512,7 @@ namespace aspect system_matrix.block(1,0).vmult(ptmp,wtmp); dst=0; - solver.solve(matrix, + solver.solve(BC_invBT, dst, ptmp, mp_preconditioner); From abc7524123c19e2092e89504de5b59964eb7c713 Mon Sep 17 00:00:00 2001 From: Quang Hoang Date: Tue, 17 Jun 2025 11:12:45 -0600 Subject: [PATCH 3/4] fixup --- source/simulator/solver.cc | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/source/simulator/solver.cc b/source/simulator/solver.cc index 2d5a8a37b22..b0668d4adf9 100644 --- a/source/simulator/solver.cc +++ b/source/simulator/solver.cc @@ -396,7 +396,7 @@ namespace aspect template - LinearOperator diag_operator(LinearOperator &exemplar, const TrilinosWrappers::MPI::Vector &diagonal) + LinearOperator diag_operator(const LinearOperator &exemplar, const TrilinosWrappers::MPI::Vector &diagonal) { LinearOperator return_op; @@ -499,7 +499,7 @@ namespace aspect - solver.solve(matrix, + solver.solve(BC_invBT, ptmp, src, mp_preconditioner); From c0eaef7b682187772f4340e45dd659c69fe7705f Mon Sep 17 00:00:00 2001 From: Quang Hoang Date: Tue, 17 Jun 2025 11:30:53 -0600 Subject: [PATCH 4/4] fixup --- source/simulator/solver.cc | 10 +++------- 1 file changed, 3 insertions(+), 7 deletions(-) diff --git a/source/simulator/solver.cc b/source/simulator/solver.cc index b0668d4adf9..7190a888302 100644 --- a/source/simulator/solver.cc +++ b/source/simulator/solver.cc @@ -170,7 +170,6 @@ namespace aspect - /** * Implement the block Schur preconditioner for the Stokes system. */ @@ -391,16 +390,15 @@ namespace aspect /** * Given a diagonal matrix stored as a vector, * create an operator that represents its action. - */ + */ template - LinearOperator diag_operator(const LinearOperator &exemplar, const TrilinosWrappers::MPI::Vector &diagonal) + LinearOperator diag_operator(const LinearOperator &exemplar, + const TrilinosWrappers::MPI::Vector &diagonal) { LinearOperator return_op; - - return_op.reinit_range_vector = exemplar.reinit_range_vector; return_op.reinit_domain_vector = exemplar.reinit_domain_vector; @@ -497,8 +495,6 @@ namespace aspect 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,