From cf99210e9aa993be0bf6f9eae18dae26c75ff2f8 Mon Sep 17 00:00:00 2001 From: Quang Hoang Date: Sun, 22 Feb 2026 22:29:33 -0500 Subject: [PATCH 01/37] create BBlockOperatorclass --- .../simulator/solver/matrix_free_operators.h | 61 +++++++++++++++++++ 1 file changed, 61 insertions(+) diff --git a/include/aspect/simulator/solver/matrix_free_operators.h b/include/aspect/simulator/solver/matrix_free_operators.h index 67bff4ee07b..940f9ae9562 100644 --- a/include/aspect/simulator/solver/matrix_free_operators.h +++ b/include/aspect/simulator/solver/matrix_free_operators.h @@ -340,6 +340,67 @@ namespace aspect */ const OperatorCellData *cell_data; }; + /** + * Operator for the B^T block. + */ + template + class BBlockOperator + : public MatrixFreeOperators::Base> + { + public: + + /** + * Constructor. + */ + BBlockOperator (); + + /** + * Reset object. + */ + void clear () override; + + /** + * Pass in a reference to the problem data. + */ + void set_cell_data (const OperatorCellData &data); + + /** + * Computes the diagonal of the matrix. Since matrix-free operators have not access + * to matrix elements, we must apply the matrix-free operator to the unit vectors to + * recover the diagonal. + */ + void compute_diagonal () override; + + private: + + /** + * Performs the application of the matrix-free operator. This function is called by + * vmult() functions MatrixFreeOperators::Base. + */ + void apply_add (dealii::LinearAlgebra::distributed::BlockVector &dst, + const dealii::LinearAlgebra::distributed::BlockVector &src) const override; + + /** + * Defines the application of the cell matrix. + */ + void local_apply (const dealii::MatrixFree &data, + dealii::LinearAlgebra::distributed::BlockVector &dst, + const dealii::LinearAlgebra::distributed::BlockVector &src, + const std::pair &cell_range) const; + + /** + * This function doesn't do anything, it's created to use the matrixfree loop. + */ + void local_apply_face (const dealii::MatrixFree &data, + dealii::LinearAlgebra::distributed::BlockVector &dst, + const dealii::LinearAlgebra::distributed::BlockVector &src, + const std::pair &face_range) const; + + /** + * A pointer to the current cell data that contains viscosity and other required parameters per cell. + */ + const OperatorCellData *cell_data; + }; From d4bf50126e706a50d00c6485052cbb5c06538051 Mon Sep 17 00:00:00 2001 From: Quang Hoang Date: Sun, 22 Feb 2026 23:10:13 -0500 Subject: [PATCH 02/37] implement functions --- .../simulator/solver/matrix_free_operators.cc | 45 +++++++++++++++++++ 1 file changed, 45 insertions(+) diff --git a/source/simulator/solver/matrix_free_operators.cc b/source/simulator/solver/matrix_free_operators.cc index c1f82d75f5a..fbce431bb9b 100644 --- a/source/simulator/solver/matrix_free_operators.cc +++ b/source/simulator/solver/matrix_free_operators.cc @@ -519,6 +519,51 @@ namespace aspect data->cell_loop(&BTBlockOperator::local_apply, this, dst, src); } + + + /** + *Operator for B block. + */ + template + MatrixFreeStokesOperators::BBlockOperator::BBlockOperator (): + MatrixFreeOperators::Base>() + {} + + template + void + MatrixFreeStokesOperators::BBlockOperator::clear () + { + this->cell_data = nullptr; + MatrixFreeOperators::Base>::clear(); + } + + template + void + MatrixFreeStokesOperators::BBlockOperator::clear () + { + this->cell_data = nullptr; + MatrixFreeOperators::Base>::clear(); + } + + + template + void + MatrixFreeStokesOperators::BBlockOperator:: + set_cell_data (const OperatorCellData &data) + { + this->cell_data = &data; + } + + + template + void + MatrixFreeStokesOperators::BBlockOperator + ::compute_diagonal () + { + // There is no need in the code for this diagonal. + Assert(false, ExcNotImplemented()); + } + /** * Mass matrix operator on pressure */ From e78b376fb33d0f5b776d5a69df72303daece1ad7 Mon Sep 17 00:00:00 2001 From: Quang Hoang Date: Mon, 23 Feb 2026 00:15:19 -0500 Subject: [PATCH 03/37] implement more functions --- .../simulator/solver/matrix_free_operators.cc | 184 ++++++++++++++++++ 1 file changed, 184 insertions(+) diff --git a/source/simulator/solver/matrix_free_operators.cc b/source/simulator/solver/matrix_free_operators.cc index fbce431bb9b..459397ed349 100644 --- a/source/simulator/solver/matrix_free_operators.cc +++ b/source/simulator/solver/matrix_free_operators.cc @@ -564,6 +564,190 @@ namespace aspect Assert(false, ExcNotImplemented()); } + + + template + void + MatrixFreeStokesOperators::BBlockOperator + ::local_apply(const dealii::MatrixFree &data, + dealii::LinearAlgebra::distributed::BlockVector &dst, + const dealii::LinearAlgebra::distributed::BlockVector &src, + const std::pair &cell_range) const + { + FEEvaluation u_eval(data, 0); + FEEvaluation p_eval(data, /*dofh*/1); + + const bool use_viscosity_at_quadrature_points + = (cell_data->viscosity.size(1) == u_eval.n_q_points); + + for (unsigned int cell=cell_range.first; cell u_eval(data, 0); + FEEvaluation p_eval(data, /*dofh*/1); + + const bool use_viscosity_at_quadrature_points + = (cell_data->viscosity.size(1) == u_eval.n_q_points); + + for (unsigned int cell=cell_range.first; cell viscosity_x_2 = 2. * cell_data->viscosity(cell,0); + + u_eval.reinit(cell); + u_eval.gather_evaluate(src.block(0), EvaluationFlags::gradients); + + p_eval.reinit(cell); + p_eval.gather_evaluate(src.block(1), EvaluationFlags::values); + + // Derivative terms related to the Newton solver + VectorizedArray deta_deps_times_sym_grad_u(0.); + VectorizedArray deta_dp_times_p(0.); + VectorizedArray eps_times_sym_grad_u(0.); + VectorizedArray eps_times_sym_grad_u_JxW(0.); + VectorizedArray theta_times_div_u(0.); + VectorizedArray theta_times_div_u_JxW(0.); + + // Pre-compute the averaged values + if (cell_data->enable_newton_derivatives + && cell_data->average_newton_factors) + { + SymmetricTensor<2,dim,VectorizedArray> sym_grad_u; + VectorizedArray val_p; + for (const unsigned int q : u_eval.quadrature_point_indices()) + { + sym_grad_u = u_eval.get_symmetric_gradient(q); + val_p = p_eval.get_value(q); + deta_deps_times_sym_grad_u += cell_data->newton_factor_wrt_strain_rate_table(cell,q) + * sym_grad_u; + deta_dp_times_p += cell_data->newton_factor_wrt_pressure_table(cell,q) * val_p; + if (cell_data->symmetrize_newton_system) + { + // The symmetrization is more complicated than it looks: + // The averaged Newton term is expressed as + // 2 (\sum_q JxW_q symgrad_phi_Iq : epsilon_q) (\sum_r w_r deta_depsilon_r : symgrad_phi_Jr) + // when symmetrized, it becomes + // (\sum_q JxW_q symgrad_phi_Iq : epsilon_q) (\sum_r w_r deta_depsilon_r : symgrad_phi_Jr) + + // (\sum_r w_r symgrad_phi_Ir : deta_depsilon_r) (\sum_q JxW_q epsilon_q : symgrad_phi_Jq) + // Let us focus on the second term. In the matrix-free framework, we should submit a + // symmetric tensor for each quadrature point q to be multiplied by the test function + // symgrad_phi_Iq and the quadrature weight JxW_q. In order to do so, we first rewrite the + // second term by swapping q and r: + // \sum_q\sum_r JxW_r w_q (symgrad_phi_Iq : deta_depsilon_q) (epsilon_r : symgrad_phi_Jr) + // In the above expression, the test function is separated out, but the quadrature weight + // is not. To separate JxW_q, we further rewrite the second term by + // \sum_q\sum_r JxW_q (JxW_r / JxW_q) w_q (symgrad_phi_Iq : deta_depsilon_q) (epsilon_r : symgrad_phi_Jr) + // Thus, the term to be submitted for each quadrature point q is + // \sum_r (JxW_r / JxW_q) w_q deta_depsilon_q (epsilon_r : symgrad_phi_Jr) + // That is the reason why we need the term (epsilon_r : symgrad_phi_Jr) JxW_r. + eps_times_sym_grad_u_JxW += cell_data->strain_rate_table(cell,q) * sym_grad_u * u_eval.JxW(q); + // The same goes for the term corresponding to compressibility/dilation + if (cell_data->is_compressible || cell_data->enable_prescribed_dilation) + theta_times_div_u_JxW += trace(cell_data->strain_rate_table(cell,q)) * trace(sym_grad_u) * u_eval.JxW(q); + } + } + } + + for (const unsigned int q : u_eval.quadrature_point_indices()) + { + // Only update the viscosity if a Q1 projection is used. + if (use_viscosity_at_quadrature_points) + viscosity_x_2 = 2. * cell_data->viscosity(cell,q); + + const SymmetricTensor<2,dim,VectorizedArray> + sym_grad_u = u_eval.get_symmetric_gradient(q); + const VectorizedArray div_u = trace(sym_grad_u); + const VectorizedArray val_p = p_eval.get_value(q); + + // Terms to be tested by phi_p: + VectorizedArray pressure_terms = + -cell_data->pressure_scaling * div_u; + + if (cell_data->enable_prescribed_dilation) + pressure_terms -= cell_data->pressure_scaling * + cell_data->pressure_scaling * + cell_data->dilation_lhs_term_table(cell,q) * + val_p; + + // Terms to be tested by the symmetric gradients of phi_u: + SymmetricTensor<2,dim,VectorizedArray> + velocity_terms = viscosity_x_2 * sym_grad_u; + + for (unsigned int d=0; dpressure_scaling * val_p; + + if (cell_data->is_compressible || + cell_data->enable_prescribed_dilation) + for (unsigned int d=0; denable_newton_derivatives) + { + if (cell_data->average_newton_factors) + { + if (cell_data->symmetrize_newton_system) + eps_times_sym_grad_u = eps_times_sym_grad_u_JxW / u_eval.JxW(q); + } + else + { + deta_deps_times_sym_grad_u = cell_data->newton_factor_wrt_strain_rate_table(cell,q) + * sym_grad_u; + deta_dp_times_p = cell_data->newton_factor_wrt_pressure_table(cell,q) * val_p; + if (cell_data->symmetrize_newton_system) + eps_times_sym_grad_u = cell_data->strain_rate_table(cell,q) * sym_grad_u; + } + + velocity_terms += + ( cell_data->symmetrize_newton_system ? + ( cell_data->strain_rate_table(cell,q) * deta_deps_times_sym_grad_u + + cell_data->newton_factor_wrt_strain_rate_table(cell,q) * eps_times_sym_grad_u ) : + 2. * cell_data->strain_rate_table(cell,q) * deta_deps_times_sym_grad_u ) + + + 2. * cell_data->strain_rate_table(cell,q) * deta_dp_times_p; + + if (cell_data->is_compressible || + cell_data->enable_prescribed_dilation) + { + constexpr number one_third = 1.0 / 3.0; + if (cell_data->symmetrize_newton_system) + { + if (cell_data->average_newton_factors) + theta_times_div_u = theta_times_div_u_JxW / u_eval.JxW(q); + else + theta_times_div_u = trace(cell_data->strain_rate_table(cell,q)) * div_u; + + velocity_terms -= one_third * cell_data->newton_factor_wrt_strain_rate_table(cell,q) * theta_times_div_u; + for (unsigned int d = 0; d < dim; ++d) + velocity_terms[d][d] -= one_third * trace(cell_data->strain_rate_table(cell,q)) * deta_deps_times_sym_grad_u; + } + else + { + for (unsigned int d = 0; d < dim; ++d) + velocity_terms[d][d] -= 2.0 * one_third * trace(cell_data->strain_rate_table(cell,q)) * deta_deps_times_sym_grad_u; + } + + if (cell_data->enable_prescribed_dilation) + { + pressure_terms += ( ( cell_data->dilation_derivative_wrt_strain_rate_table(cell,q) + * sym_grad_u ) + + + ( cell_data->dilation_derivative_wrt_pressure_table(cell,q) + * cell_data->pressure_scaling * val_p ) + ) + * cell_data->pressure_scaling; + } + } + } + + u_eval.submit_symmetric_gradient(velocity_terms, q); + p_eval.submit_value(pressure_terms, q); + } + + u_eval.integrate_scatter(EvaluationFlags::gradients, dst.block(0)); + p_eval.integrate_scatter(EvaluationFlags::values, dst.block(1)); + } + } + /** * Mass matrix operator on pressure */ From a7c3a1a212d9e0dfe7a17c501f8f08d864ae546a Mon Sep 17 00:00:00 2001 From: Quang Hoang Date: Mon, 23 Feb 2026 01:24:09 -0500 Subject: [PATCH 04/37] implement local_apply --- .../simulator/solver/matrix_free_operators.cc | 154 +----------------- 1 file changed, 7 insertions(+), 147 deletions(-) diff --git a/source/simulator/solver/matrix_free_operators.cc b/source/simulator/solver/matrix_free_operators.cc index 459397ed349..ebf4eeb2c49 100644 --- a/source/simulator/solver/matrix_free_operators.cc +++ b/source/simulator/solver/matrix_free_operators.cc @@ -577,173 +577,33 @@ namespace aspect FEEvaluation u_eval(data, 0); FEEvaluation p_eval(data, /*dofh*/1); - const bool use_viscosity_at_quadrature_points - = (cell_data->viscosity.size(1) == u_eval.n_q_points); - - for (unsigned int cell=cell_range.first; cell u_eval(data, 0); - FEEvaluation p_eval(data, /*dofh*/1); - - const bool use_viscosity_at_quadrature_points - = (cell_data->viscosity.size(1) == u_eval.n_q_points); for (unsigned int cell=cell_range.first; cell viscosity_x_2 = 2. * cell_data->viscosity(cell,0); u_eval.reinit(cell); u_eval.gather_evaluate(src.block(0), EvaluationFlags::gradients); p_eval.reinit(cell); - p_eval.gather_evaluate(src.block(1), EvaluationFlags::values); - - // Derivative terms related to the Newton solver - VectorizedArray deta_deps_times_sym_grad_u(0.); - VectorizedArray deta_dp_times_p(0.); - VectorizedArray eps_times_sym_grad_u(0.); - VectorizedArray eps_times_sym_grad_u_JxW(0.); - VectorizedArray theta_times_div_u(0.); - VectorizedArray theta_times_div_u_JxW(0.); - - // Pre-compute the averaged values - if (cell_data->enable_newton_derivatives - && cell_data->average_newton_factors) - { - SymmetricTensor<2,dim,VectorizedArray> sym_grad_u; - VectorizedArray val_p; - for (const unsigned int q : u_eval.quadrature_point_indices()) - { - sym_grad_u = u_eval.get_symmetric_gradient(q); - val_p = p_eval.get_value(q); - deta_deps_times_sym_grad_u += cell_data->newton_factor_wrt_strain_rate_table(cell,q) - * sym_grad_u; - deta_dp_times_p += cell_data->newton_factor_wrt_pressure_table(cell,q) * val_p; - if (cell_data->symmetrize_newton_system) - { - // The symmetrization is more complicated than it looks: - // The averaged Newton term is expressed as - // 2 (\sum_q JxW_q symgrad_phi_Iq : epsilon_q) (\sum_r w_r deta_depsilon_r : symgrad_phi_Jr) - // when symmetrized, it becomes - // (\sum_q JxW_q symgrad_phi_Iq : epsilon_q) (\sum_r w_r deta_depsilon_r : symgrad_phi_Jr) + - // (\sum_r w_r symgrad_phi_Ir : deta_depsilon_r) (\sum_q JxW_q epsilon_q : symgrad_phi_Jq) - // Let us focus on the second term. In the matrix-free framework, we should submit a - // symmetric tensor for each quadrature point q to be multiplied by the test function - // symgrad_phi_Iq and the quadrature weight JxW_q. In order to do so, we first rewrite the - // second term by swapping q and r: - // \sum_q\sum_r JxW_r w_q (symgrad_phi_Iq : deta_depsilon_q) (epsilon_r : symgrad_phi_Jr) - // In the above expression, the test function is separated out, but the quadrature weight - // is not. To separate JxW_q, we further rewrite the second term by - // \sum_q\sum_r JxW_q (JxW_r / JxW_q) w_q (symgrad_phi_Iq : deta_depsilon_q) (epsilon_r : symgrad_phi_Jr) - // Thus, the term to be submitted for each quadrature point q is - // \sum_r (JxW_r / JxW_q) w_q deta_depsilon_q (epsilon_r : symgrad_phi_Jr) - // That is the reason why we need the term (epsilon_r : symgrad_phi_Jr) JxW_r. - eps_times_sym_grad_u_JxW += cell_data->strain_rate_table(cell,q) * sym_grad_u * u_eval.JxW(q); - // The same goes for the term corresponding to compressibility/dilation - if (cell_data->is_compressible || cell_data->enable_prescribed_dilation) - theta_times_div_u_JxW += trace(cell_data->strain_rate_table(cell,q)) * trace(sym_grad_u) * u_eval.JxW(q); - } - } - } + + SymmetricTensor<2,dim,VectorizedArray> sym_grad_u; + for (const unsigned int q : u_eval.quadrature_point_indices()) { - // Only update the viscosity if a Q1 projection is used. - if (use_viscosity_at_quadrature_points) - viscosity_x_2 = 2. * cell_data->viscosity(cell,q); - - const SymmetricTensor<2,dim,VectorizedArray> sym_grad_u = u_eval.get_symmetric_gradient(q); - const VectorizedArray div_u = trace(sym_grad_u); - const VectorizedArray val_p = p_eval.get_value(q); - // Terms to be tested by phi_p: + const VectorizedArray div_u = trace(sym_grad_u); VectorizedArray pressure_terms = -cell_data->pressure_scaling * div_u; + p_eval.submit_value(pressure_terms, q); - if (cell_data->enable_prescribed_dilation) - pressure_terms -= cell_data->pressure_scaling * - cell_data->pressure_scaling * - cell_data->dilation_lhs_term_table(cell,q) * - val_p; - - // Terms to be tested by the symmetric gradients of phi_u: - SymmetricTensor<2,dim,VectorizedArray> - velocity_terms = viscosity_x_2 * sym_grad_u; - - for (unsigned int d=0; dpressure_scaling * val_p; - - if (cell_data->is_compressible || - cell_data->enable_prescribed_dilation) - for (unsigned int d=0; denable_newton_derivatives) - { - if (cell_data->average_newton_factors) - { - if (cell_data->symmetrize_newton_system) - eps_times_sym_grad_u = eps_times_sym_grad_u_JxW / u_eval.JxW(q); - } - else - { - deta_deps_times_sym_grad_u = cell_data->newton_factor_wrt_strain_rate_table(cell,q) - * sym_grad_u; - deta_dp_times_p = cell_data->newton_factor_wrt_pressure_table(cell,q) * val_p; - if (cell_data->symmetrize_newton_system) - eps_times_sym_grad_u = cell_data->strain_rate_table(cell,q) * sym_grad_u; - } - - velocity_terms += - ( cell_data->symmetrize_newton_system ? - ( cell_data->strain_rate_table(cell,q) * deta_deps_times_sym_grad_u + - cell_data->newton_factor_wrt_strain_rate_table(cell,q) * eps_times_sym_grad_u ) : - 2. * cell_data->strain_rate_table(cell,q) * deta_deps_times_sym_grad_u ) - + - 2. * cell_data->strain_rate_table(cell,q) * deta_dp_times_p; - - if (cell_data->is_compressible || - cell_data->enable_prescribed_dilation) - { - constexpr number one_third = 1.0 / 3.0; - if (cell_data->symmetrize_newton_system) - { - if (cell_data->average_newton_factors) - theta_times_div_u = theta_times_div_u_JxW / u_eval.JxW(q); - else - theta_times_div_u = trace(cell_data->strain_rate_table(cell,q)) * div_u; - - velocity_terms -= one_third * cell_data->newton_factor_wrt_strain_rate_table(cell,q) * theta_times_div_u; - for (unsigned int d = 0; d < dim; ++d) - velocity_terms[d][d] -= one_third * trace(cell_data->strain_rate_table(cell,q)) * deta_deps_times_sym_grad_u; - } - else - { - for (unsigned int d = 0; d < dim; ++d) - velocity_terms[d][d] -= 2.0 * one_third * trace(cell_data->strain_rate_table(cell,q)) * deta_deps_times_sym_grad_u; - } - - if (cell_data->enable_prescribed_dilation) - { - pressure_terms += ( ( cell_data->dilation_derivative_wrt_strain_rate_table(cell,q) - * sym_grad_u ) - + - ( cell_data->dilation_derivative_wrt_pressure_table(cell,q) - * cell_data->pressure_scaling * val_p ) - ) - * cell_data->pressure_scaling; - } - } - } - u_eval.submit_symmetric_gradient(velocity_terms, q); - p_eval.submit_value(pressure_terms, q); + } + - u_eval.integrate_scatter(EvaluationFlags::gradients, dst.block(0)); p_eval.integrate_scatter(EvaluationFlags::values, dst.block(1)); } } From e4b31c0c45d0b889b91583fda496a4440c7c7d24 Mon Sep 17 00:00:00 2001 From: Quang Hoang Date: Mon, 2 Mar 2026 14:00:07 -0500 Subject: [PATCH 05/37] cases --- source/simulator/solver.cc | 6 ++++-- 1 file changed, 4 insertions(+), 2 deletions(-) diff --git a/source/simulator/solver.cc b/source/simulator/solver.cc index ceba95be9b1..2a47b1a9f4f 100644 --- a/source/simulator/solver.cc +++ b/source/simulator/solver.cc @@ -578,8 +578,10 @@ namespace aspect const std::string name = [&]() -> std::string { - if (parameters.stokes_solver_type == Parameters::StokesSolverType::block_gmg) + if (parameters.stokes_solver_type == Parameters::StokesSolverType::block_gmg && parameters.use_bfbt == false) return stokes_matrix_free->name(); + if (parameters.stokes_solver_type == Parameters::StokesSolverType::block_gmg && parameters.use_bfbt == true) + return "GMG-BFBT"; if (parameters.use_direct_stokes_solver) return "direct"; if (parameters.use_bfbt) @@ -772,7 +774,7 @@ namespace aspect solver_control_expensive.enable_history_data(); std::unique_ptr schur; - if (parameters.use_bfbt) + if (parameters.use_bfbt && ) { schur = std::make_unique>( system_preconditioner_matrix.block(pressure_block_index,pressure_block_index), From fd7462bfd3f9782fce658230019359ceadf540f9 Mon Sep 17 00:00:00 2001 From: Quang Hoang Date: Mon, 2 Mar 2026 14:53:32 -0500 Subject: [PATCH 06/37] fixes --- .../simulator/solver/matrix_free_operators.h | 2 +- source/simulator/solver.cc | 5 +++- .../simulator/solver/matrix_free_operators.cc | 29 +++++++------------ 3 files changed, 16 insertions(+), 20 deletions(-) diff --git a/include/aspect/simulator/solver/matrix_free_operators.h b/include/aspect/simulator/solver/matrix_free_operators.h index 940f9ae9562..a3ad3db14c2 100644 --- a/include/aspect/simulator/solver/matrix_free_operators.h +++ b/include/aspect/simulator/solver/matrix_free_operators.h @@ -343,7 +343,7 @@ namespace aspect /** * Operator for the B^T block. */ - template + template class BBlockOperator : public MatrixFreeOperators::Base> { diff --git a/source/simulator/solver.cc b/source/simulator/solver.cc index 2a47b1a9f4f..21d83ca02fa 100644 --- a/source/simulator/solver.cc +++ b/source/simulator/solver.cc @@ -180,6 +180,9 @@ namespace aspect }; + + + /** * 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}, @@ -774,7 +777,7 @@ namespace aspect solver_control_expensive.enable_history_data(); std::unique_ptr schur; - if (parameters.use_bfbt && ) + if (parameters.use_bfbt) { schur = std::make_unique>( system_preconditioner_matrix.block(pressure_block_index,pressure_block_index), diff --git a/source/simulator/solver/matrix_free_operators.cc b/source/simulator/solver/matrix_free_operators.cc index ebf4eeb2c49..eb0958ba481 100644 --- a/source/simulator/solver/matrix_free_operators.cc +++ b/source/simulator/solver/matrix_free_operators.cc @@ -524,40 +524,33 @@ namespace aspect /** *Operator for B block. */ - template - MatrixFreeStokesOperators::BBlockOperator::BBlockOperator (): + template + MatrixFreeStokesOperators::BBlockOperator::BBlockOperator (): MatrixFreeOperators::Base>() {} - template + template void - MatrixFreeStokesOperators::BBlockOperator::clear () + MatrixFreeStokesOperators::BBlockOperator::clear () { this->cell_data = nullptr; MatrixFreeOperators::Base>::clear(); } - template - void - MatrixFreeStokesOperators::BBlockOperator::clear () - { - this->cell_data = nullptr; - MatrixFreeOperators::Base>::clear(); - } - - template + + template void - MatrixFreeStokesOperators::BBlockOperator:: + MatrixFreeStokesOperators::BBlockOperator:: set_cell_data (const OperatorCellData &data) { this->cell_data = &data; } - template + template void - MatrixFreeStokesOperators::BBlockOperator + MatrixFreeStokesOperators::BBlockOperator ::compute_diagonal () { // There is no need in the code for this diagonal. @@ -566,9 +559,9 @@ namespace aspect - template + template void - MatrixFreeStokesOperators::BBlockOperator + MatrixFreeStokesOperators::BBlockOperator ::local_apply(const dealii::MatrixFree &data, dealii::LinearAlgebra::distributed::BlockVector &dst, const dealii::LinearAlgebra::distributed::BlockVector &src, From 4531bdc51c24deb7300354a22e7423354e8f824e Mon Sep 17 00:00:00 2001 From: Quang Hoang Date: Mon, 2 Mar 2026 18:49:33 -0500 Subject: [PATCH 07/37] create diagBFBT class --- source/simulator/solver.cc | 32 ++++++++++++++++++++++++++++++++ 1 file changed, 32 insertions(+) diff --git a/source/simulator/solver.cc b/source/simulator/solver.cc index 21d83ca02fa..a8d680c45ab 100644 --- a/source/simulator/solver.cc +++ b/source/simulator/solver.cc @@ -224,6 +224,38 @@ namespace aspect const LinearAlgebra::BlockSparseMatrix &system_matrix; }; + + template + class DiagBFBT: public SchurComplementOperator + { + public: + /** + * Constructor. + * @param laplace_preconditioner The preconditioner for @p pressure_laplace_matrix + * @param solver_tolerance The relative solver tolerance for the inner solve + * @param diag_A Diagonal of A used in the BFBT preconditioner. + * @param system_matrix Sparse block matrix storing the Stokes system of the form + * [A B^T + * B 0]. + */ + DiagBFBT(const PreconditionerMp &mp_preconditioner, + const double solver_tolerance, + const LinearAlgebra::Vector &diag_A, + const StokesMatrixType &system_matrix); + + void vmult(VectorType &dst, + const VectorType &src) const override; + + unsigned int n_iterations() const override; + + private: + mutable unsigned int n_iterations_; + const PreconditionerMp &laplace_preconditioner; + const double solver_tolerance; + const LinearAlgebra::Vector &diag_A; + const StokesMatrixType &system_matrix; + }; + template WeightedBFBT::WeightedBFBT( const LinearAlgebra::SparseMatrix &pressure_laplace_matrix, From b4a4e63ac09acb3b858387814843008ad1b54b3e Mon Sep 17 00:00:00 2001 From: Quang Hoang Date: Mon, 2 Mar 2026 21:14:23 -0500 Subject: [PATCH 08/37] move stuff --- .../solver/block_stokes_preconditioner.h | 31 +++++++++++++++++++ 1 file changed, 31 insertions(+) diff --git a/include/aspect/simulator/solver/block_stokes_preconditioner.h b/include/aspect/simulator/solver/block_stokes_preconditioner.h index 6f21a785a6e..a63f0b2fcad 100644 --- a/include/aspect/simulator/solver/block_stokes_preconditioner.h +++ b/include/aspect/simulator/solver/block_stokes_preconditioner.h @@ -304,6 +304,37 @@ namespace aspect n_iterations_Schur_complement_(0) {} + template + class DiagBFBT: public SchurComplementOperator + { + public: + /** + * Constructor. + * @param laplace_preconditioner The preconditioner for @p pressure_laplace_matrix + * @param solver_tolerance The relative solver tolerance for the inner solve + * @param diag_A Diagonal of A used in the BFBT preconditioner. + * @param system_matrix Sparse block matrix storing the Stokes system of the form + * [A B^T + * B 0]. + */ + DiagBFBT(const PreconditionerMp &mp_preconditioner, + const double solver_tolerance, + const LinearAlgebra::Vector &diag_A, + const StokesMatrixType &system_matrix); + + void vmult(VectorType &dst, + const VectorType &src) const override; + + unsigned int n_iterations() const override; + + private: + mutable unsigned int n_iterations_; + const PreconditionerMp &laplace_preconditioner; + const double solver_tolerance; + const LinearAlgebra::Vector &diag_A; + const StokesMatrixType &system_matrix; + }; + template void SchurApproximation::vmult(VectorType &dst, From d5dea6d50b078d1477e0807931fa1dcfdb887288 Mon Sep 17 00:00:00 2001 From: Quang Hoang Date: Wed, 4 Mar 2026 16:53:18 -0500 Subject: [PATCH 09/37] move schur operator --- .../solver/block_stokes_preconditioner.h | 80 ++++++++++++------- .../stokes_matrix_free_local_smoothing.cc | 30 +++++++ 2 files changed, 82 insertions(+), 28 deletions(-) diff --git a/include/aspect/simulator/solver/block_stokes_preconditioner.h b/include/aspect/simulator/solver/block_stokes_preconditioner.h index a63f0b2fcad..c8453a4f8db 100644 --- a/include/aspect/simulator/solver/block_stokes_preconditioner.h +++ b/include/aspect/simulator/solver/block_stokes_preconditioner.h @@ -304,36 +304,10 @@ namespace aspect n_iterations_Schur_complement_(0) {} - template - class DiagBFBT: public SchurComplementOperator - { - public: - /** - * Constructor. - * @param laplace_preconditioner The preconditioner for @p pressure_laplace_matrix - * @param solver_tolerance The relative solver tolerance for the inner solve - * @param diag_A Diagonal of A used in the BFBT preconditioner. - * @param system_matrix Sparse block matrix storing the Stokes system of the form - * [A B^T - * B 0]. - */ - DiagBFBT(const PreconditionerMp &mp_preconditioner, - const double solver_tolerance, - const LinearAlgebra::Vector &diag_A, - const StokesMatrixType &system_matrix); + - void vmult(VectorType &dst, - const VectorType &src) const override; - unsigned int n_iterations() const override; - - private: - mutable unsigned int n_iterations_; - const PreconditionerMp &laplace_preconditioner; - const double solver_tolerance; - const LinearAlgebra::Vector &diag_A; - const StokesMatrixType &system_matrix; - }; + template @@ -381,6 +355,9 @@ namespace aspect n_iterations_Schur_complement_ += 1; } } + + + template unsigned int SchurApproximation::n_iterations() const { @@ -388,6 +365,53 @@ namespace aspect } + + /** + * Base class for Schur Complement operators. + */ + template + class SchurComplementOperator + { + public: + virtual ~SchurComplementOperator() = default; + + virtual void vmult(VectorType &dst, + const VectorType &src) const=0; + virtual unsigned int n_iterations() const=0; + + }; + template + class DiagBFBT: public SchurComplementOperator + { + public: + /** + * Constructor. + * @param laplace_preconditioner The preconditioner for @p pressure_laplace_matrix + * @param solver_tolerance The relative solver tolerance for the inner solve + * @param diag_A Diagonal of A used in the BFBT preconditioner. + * @param system_matrix Sparse block matrix storing the Stokes system of the form + * [A B^T + * B 0]. + */ + DiagBFBT(const PreconditionerMp &mp_preconditioner, + const double solver_tolerance, + const LinearAlgebra::Vector &diag_A, + const StokesMatrixType &system_matrix); + + void vmult(VectorType &dst, + const VectorType &src) const override; + + unsigned int n_iterations() const override; + + private: + mutable unsigned int n_iterations_; + const PreconditionerMp &laplace_preconditioner; + const double solver_tolerance; + const LinearAlgebra::Vector &diag_A; + const StokesMatrixType &system_matrix; + }; + + } } diff --git a/source/simulator/solver/stokes_matrix_free_local_smoothing.cc b/source/simulator/solver/stokes_matrix_free_local_smoothing.cc index 9df74f2001b..abee46125ad 100644 --- a/source/simulator/solver/stokes_matrix_free_local_smoothing.cc +++ b/source/simulator/solver/stokes_matrix_free_local_smoothing.cc @@ -45,6 +45,36 @@ namespace aspect { + // template + // class DiagBFBT: public SchurComplementOperator + // { + // public: + // /** + // * Constructor. + // * @param laplace_preconditioner The preconditioner for @p pressure_laplace_matrix + // * @param solver_tolerance The relative solver tolerance for the inner solve + // * @param diag_A Diagonal of A used in the BFBT preconditioner. + // * @param system_matrix Sparse block matrix storing the Stokes system of the form + // * [A B^T + // * B 0]. + // */ + // DiagBFBT(const PreconditionerMp &mp_preconditioner, + // const double solver_tolerance, + // const LinearAlgebra::Vector &diag_A, + // const StokesMatrixType &system_matrix); + + // void vmult(VectorType &dst, + // const VectorType &src) const override; + + // unsigned int n_iterations() const override; + + // private: + // mutable unsigned int n_iterations_; + // const PreconditionerMp &laplace_preconditioner; + // const double solver_tolerance; + // const LinearAlgebra::Vector &diag_A; + // const StokesMatrixType &system_matrix; + // }; template void StokesMatrixFreeHandlerLocalSmoothingImplementation::declare_parameters(ParameterHandler &prm) From 0af2e4c135661548c6a6cfe09f018498f40d2863 Mon Sep 17 00:00:00 2001 From: Quang Hoang Date: Tue, 10 Mar 2026 14:41:54 -0400 Subject: [PATCH 10/37] refactor --- source/simulator/solver.cc | 53 +++++--------------------------------- 1 file changed, 6 insertions(+), 47 deletions(-) diff --git a/source/simulator/solver.cc b/source/simulator/solver.cc index a8d680c45ab..c5112bc3451 100644 --- a/source/simulator/solver.cc +++ b/source/simulator/solver.cc @@ -166,19 +166,7 @@ namespace aspect return dst.l2_norm(); } - /** - * Base class for Schur Complement operators. - */ - class SchurComplementOperator - { - public: - virtual ~SchurComplementOperator() = default; - virtual void vmult(LinearAlgebra::Vector &dst, - const LinearAlgebra::Vector &src) const=0; - virtual unsigned int n_iterations() const=0; - - }; @@ -191,7 +179,7 @@ namespace aspect * velocity mass matrix. */ template - class WeightedBFBT: public SchurComplementOperator + class WeightedBFBT: public SchurComplementOperator { public: /** @@ -225,36 +213,7 @@ namespace aspect }; - template - class DiagBFBT: public SchurComplementOperator - { - public: - /** - * Constructor. - * @param laplace_preconditioner The preconditioner for @p pressure_laplace_matrix - * @param solver_tolerance The relative solver tolerance for the inner solve - * @param diag_A Diagonal of A used in the BFBT preconditioner. - * @param system_matrix Sparse block matrix storing the Stokes system of the form - * [A B^T - * B 0]. - */ - DiagBFBT(const PreconditionerMp &mp_preconditioner, - const double solver_tolerance, - const LinearAlgebra::Vector &diag_A, - const StokesMatrixType &system_matrix); - - void vmult(VectorType &dst, - const VectorType &src) const override; - - unsigned int n_iterations() const override; - - private: - mutable unsigned int n_iterations_; - const PreconditionerMp &laplace_preconditioner; - const double solver_tolerance; - const LinearAlgebra::Vector &diag_A; - const StokesMatrixType &system_matrix; - }; + template WeightedBFBT::WeightedBFBT( @@ -344,7 +303,7 @@ namespace aspect * PreconditionerMp passed to the constructor. */ template - class InverseWeightedMassMatrix: public SchurComplementOperator + class InverseWeightedMassMatrix: public SchurComplementOperator { public: /** @@ -808,7 +767,7 @@ namespace aspect solver_control_cheap.enable_history_data(); solver_control_expensive.enable_history_data(); - std::unique_ptr schur; + std::unique_ptr> schur; if (parameters.use_bfbt) { schur = std::make_unique>( @@ -835,7 +794,7 @@ namespace aspect stokes_A_block_is_symmetric(), parameters.linear_solver_A_block_tolerance); const internal::BlockSchurPreconditioner, - internal::SchurComplementOperator, LinearAlgebra::SparseMatrix, LinearAlgebra::BlockVector> + internal::SchurComplementOperator, LinearAlgebra::SparseMatrix, LinearAlgebra::BlockVector> preconditioner_cheap ( inverse_velocity_block_cheap, *schur, @@ -849,7 +808,7 @@ namespace aspect stokes_A_block_is_symmetric(), parameters.linear_solver_A_block_tolerance); const internal::BlockSchurPreconditioner, - internal::SchurComplementOperator, LinearAlgebra::SparseMatrix, LinearAlgebra::BlockVector> + internal::SchurComplementOperator, LinearAlgebra::SparseMatrix, LinearAlgebra::BlockVector> preconditioner_expensive ( inverse_velocity_block_expensive, *schur, From f2c2c058e288bc0d8e4103ebe85a303fe3195378 Mon Sep 17 00:00:00 2001 From: Quang Hoang Date: Tue, 10 Mar 2026 15:22:51 -0400 Subject: [PATCH 11/37] create constructor diag BFBT --- .../solver/block_stokes_preconditioner.h | 4 +- .../stokes_matrix_free_local_smoothing.cc | 47 +++++++------------ 2 files changed, 19 insertions(+), 32 deletions(-) diff --git a/include/aspect/simulator/solver/block_stokes_preconditioner.h b/include/aspect/simulator/solver/block_stokes_preconditioner.h index c8453a4f8db..f581b17bfd2 100644 --- a/include/aspect/simulator/solver/block_stokes_preconditioner.h +++ b/include/aspect/simulator/solver/block_stokes_preconditioner.h @@ -393,9 +393,9 @@ namespace aspect * [A B^T * B 0]. */ - DiagBFBT(const PreconditionerMp &mp_preconditioner, + DiagBFBT(const PreconditionerMp &laplace_preconditioner, const double solver_tolerance, - const LinearAlgebra::Vector &diag_A, + const VectorType &diag_A, const StokesMatrixType &system_matrix); void vmult(VectorType &dst, diff --git a/source/simulator/solver/stokes_matrix_free_local_smoothing.cc b/source/simulator/solver/stokes_matrix_free_local_smoothing.cc index abee46125ad..fa5bd7d64fe 100644 --- a/source/simulator/solver/stokes_matrix_free_local_smoothing.cc +++ b/source/simulator/solver/stokes_matrix_free_local_smoothing.cc @@ -45,36 +45,23 @@ namespace aspect { - // template - // class DiagBFBT: public SchurComplementOperator - // { - // public: - // /** - // * Constructor. - // * @param laplace_preconditioner The preconditioner for @p pressure_laplace_matrix - // * @param solver_tolerance The relative solver tolerance for the inner solve - // * @param diag_A Diagonal of A used in the BFBT preconditioner. - // * @param system_matrix Sparse block matrix storing the Stokes system of the form - // * [A B^T - // * B 0]. - // */ - // DiagBFBT(const PreconditionerMp &mp_preconditioner, - // const double solver_tolerance, - // const LinearAlgebra::Vector &diag_A, - // const StokesMatrixType &system_matrix); - - // void vmult(VectorType &dst, - // const VectorType &src) const override; - - // unsigned int n_iterations() const override; - - // private: - // mutable unsigned int n_iterations_; - // const PreconditionerMp &laplace_preconditioner; - // const double solver_tolerance; - // const LinearAlgebra::Vector &diag_A; - // const StokesMatrixType &system_matrix; - // }; + namespace internal{ + template + DiagBFBT::DiagBFBT(const PreconditionerMp &mp_preconditioner, + const double solver_tolerance, + const VectorType &diag_A, + const StokesMatrixType & system_matrix): + n_iterations_(0), + laplace_preconditioner(laplace_preconditioner), + solver_tolerance(solver_tolerance), + diag_A(diag_A), + system_matrix(system_matrix) + + {} + + + } + template void StokesMatrixFreeHandlerLocalSmoothingImplementation::declare_parameters(ParameterHandler &prm) From e9cb65cef1d1ca1600d7a369e7e0dcf67fe28a29 Mon Sep 17 00:00:00 2001 From: Quang Hoang Date: Tue, 10 Mar 2026 20:38:53 -0400 Subject: [PATCH 12/37] implement n iterations function --- .../simulator/solver/stokes_matrix_free_local_smoothing.cc | 7 +++++++ 1 file changed, 7 insertions(+) diff --git a/source/simulator/solver/stokes_matrix_free_local_smoothing.cc b/source/simulator/solver/stokes_matrix_free_local_smoothing.cc index fa5bd7d64fe..de1d097d029 100644 --- a/source/simulator/solver/stokes_matrix_free_local_smoothing.cc +++ b/source/simulator/solver/stokes_matrix_free_local_smoothing.cc @@ -58,6 +58,13 @@ namespace aspect system_matrix(system_matrix) {} + + + + template + unsigned int DiagBFBT::n_iterations() const{ + return n_iterations_; + } } From b56213b6c32060e58abf80eaf55024846df19ca5 Mon Sep 17 00:00:00 2001 From: Quang Hoang Date: Tue, 10 Mar 2026 21:58:43 -0400 Subject: [PATCH 13/37] initialize --- .../stokes_matrix_free_local_smoothing.h | 2 ++ .../simulator/solver/matrix_free_operators.cc | 2 ++ .../stokes_matrix_free_local_smoothing.cc | 28 +++++++++++++++++++ 3 files changed, 32 insertions(+) diff --git a/include/aspect/simulator/solver/stokes_matrix_free_local_smoothing.h b/include/aspect/simulator/solver/stokes_matrix_free_local_smoothing.h index bd00576843c..d0c5d1ae886 100644 --- a/include/aspect/simulator/solver/stokes_matrix_free_local_smoothing.h +++ b/include/aspect/simulator/solver/stokes_matrix_free_local_smoothing.h @@ -212,12 +212,14 @@ namespace aspect using SchurComplementMatrixType = MatrixFreeStokesOperators::MassMatrixOperator; using ABlockMatrixType = MatrixFreeStokesOperators::ABlockOperator; using BTBlockOperatorType = MatrixFreeStokesOperators::BTBlockOperator; + using BBlockOperatorType = MatrixFreeStokesOperators::BBlockOperator; using GMGSchurComplementMatrixType = MatrixFreeStokesOperators::MassMatrixOperator; using GMGABlockMatrixType = MatrixFreeStokesOperators::ABlockOperator; StokesMatrixType stokes_matrix; ABlockMatrixType A_block_matrix; BTBlockOperatorType BT_block; + BBlockOperatorType B_block; SchurComplementMatrixType Schur_complement_block_matrix; AffineConstraints constraints_v; diff --git a/source/simulator/solver/matrix_free_operators.cc b/source/simulator/solver/matrix_free_operators.cc index eb0958ba481..f5cb5dc5562 100644 --- a/source/simulator/solver/matrix_free_operators.cc +++ b/source/simulator/solver/matrix_free_operators.cc @@ -1016,6 +1016,8 @@ namespace aspect template class MatrixFreeStokesOperators::StokesOperator; \ template class MatrixFreeStokesOperators::BTBlockOperator; \ template class MatrixFreeStokesOperators::BTBlockOperator; \ + template class MatrixFreeStokesOperators::BBlockOperator;\ + template class MatrixFreeStokesOperators::BBlockOperator;\ template class MatrixFreeStokesOperators::MassMatrixOperator; \ template class MatrixFreeStokesOperators::MassMatrixOperator; \ template struct MatrixFreeStokesOperators::OperatorCellData; diff --git a/source/simulator/solver/stokes_matrix_free_local_smoothing.cc b/source/simulator/solver/stokes_matrix_free_local_smoothing.cc index de1d097d029..5323e8d37e4 100644 --- a/source/simulator/solver/stokes_matrix_free_local_smoothing.cc +++ b/source/simulator/solver/stokes_matrix_free_local_smoothing.cc @@ -67,6 +67,28 @@ namespace aspect } + template + void DiagBFBT::vmult(VectorType &dst, const VectorType &src){ + SolverControl solver_control(1000, src.l2_norm() * solver_tolerance); + PrimitiveVectorMemory mem; + SolverCG solver(solver_control, mem); //These vector types might be wrong + try{ + VectorType utmp; + utmp.reinit(diag_A); + VectorType ptmp; + ptmp.reinit(src); + VectorType wtmpt; + wtmp.reinit(diag_A); + { + SolverControl solver_control(5000,1e-6*src.l2_norm(),false,true); + SolverCG solver(solver_control); + const auto Op_A=LinearOperator(system_matrix.block(0,0)); + + } + } + } + + } template @@ -1829,6 +1851,12 @@ namespace aspect BT_block.initialize(matrix_free); } + //B block matrix + { + B_block.clear(); + B_block.initialize(matrix_free); + } + // Schur complement block matrix { Schur_complement_block_matrix.clear(); From a1b9835a90ad94a00d1331c318bc13c05bbeb7ae Mon Sep 17 00:00:00 2001 From: Quang Hoang Date: Wed, 11 Mar 2026 11:47:48 -0400 Subject: [PATCH 14/37] implement apply add BOperator --- .../simulator/solver/matrix_free_operators.cc | 10 +++++ .../stokes_matrix_free_local_smoothing.cc | 38 +++++++++---------- 2 files changed, 29 insertions(+), 19 deletions(-) diff --git a/source/simulator/solver/matrix_free_operators.cc b/source/simulator/solver/matrix_free_operators.cc index f5cb5dc5562..196720d0343 100644 --- a/source/simulator/solver/matrix_free_operators.cc +++ b/source/simulator/solver/matrix_free_operators.cc @@ -601,6 +601,16 @@ namespace aspect } } + template + void + MatrixFreeStokesOperators::BBlockOperator + ::apply_add (dealii::LinearAlgebra::distributed::BlockVector &dst, + const dealii::LinearAlgebra::distributed::BlockVector &src) const + { + MatrixFreeOperators::Base>:: + data->cell_loop(&BBlockOperator::local_apply, this, dst, src); + } + /** * Mass matrix operator on pressure */ diff --git a/source/simulator/solver/stokes_matrix_free_local_smoothing.cc b/source/simulator/solver/stokes_matrix_free_local_smoothing.cc index 5323e8d37e4..12329a49afd 100644 --- a/source/simulator/solver/stokes_matrix_free_local_smoothing.cc +++ b/source/simulator/solver/stokes_matrix_free_local_smoothing.cc @@ -67,26 +67,26 @@ namespace aspect } - template - void DiagBFBT::vmult(VectorType &dst, const VectorType &src){ - SolverControl solver_control(1000, src.l2_norm() * solver_tolerance); - PrimitiveVectorMemory mem; - SolverCG solver(solver_control, mem); //These vector types might be wrong - try{ - VectorType utmp; - utmp.reinit(diag_A); - VectorType ptmp; - ptmp.reinit(src); - VectorType wtmpt; - wtmp.reinit(diag_A); - { - SolverControl solver_control(5000,1e-6*src.l2_norm(),false,true); - SolverCG solver(solver_control); - const auto Op_A=LinearOperator(system_matrix.block(0,0)); + // template + // void DiagBFBT::vmult(VectorType &dst, const VectorType &src){ + // SolverControl solver_control(1000, src.l2_norm() * solver_tolerance); + // PrimitiveVectorMemory mem; + // SolverCG solver(solver_control, mem); //These vector types might be wrong + // try{ + // VectorType utmp; + // utmp.reinit(diag_A); + // VectorType ptmp; + // ptmp.reinit(src); + // VectorType wtmpt; + // wtmp.reinit(diag_A); + // { + // SolverControl solver_control(5000,1e-6*src.l2_norm(),false,true); + // SolverCG solver(solver_control); + // const auto Op_A=LinearOperator(system_matrix.block(0,0)); - } - } - } + // } + // } + // } } From 86b26ab3464250003dc56e4bf5da40a93df11186 Mon Sep 17 00:00:00 2001 From: Quang Hoang Date: Wed, 11 Mar 2026 20:18:47 -0400 Subject: [PATCH 15/37] BC_invBT operator --- .../solver/block_stokes_preconditioner.h | 23 ++++++++++++++++ .../stokes_matrix_free_local_smoothing.cc | 27 +++++++++++++++++++ 2 files changed, 50 insertions(+) diff --git a/include/aspect/simulator/solver/block_stokes_preconditioner.h b/include/aspect/simulator/solver/block_stokes_preconditioner.h index f581b17bfd2..ff0e5720b5e 100644 --- a/include/aspect/simulator/solver/block_stokes_preconditioner.h +++ b/include/aspect/simulator/solver/block_stokes_preconditioner.h @@ -380,6 +380,29 @@ namespace aspect virtual unsigned int n_iterations() const=0; }; + + + template + class BC_invBT_Operator{ + public: + BC_invBT_Operator(const StokesMatrixType &system_matrix, + const BOperatorType &B_operator, + const BTOperatorType &BT_operator, + const dealii::LinearAlgebra::distributed::Vector &diag_A_inv): + system_matrix(system_matrix), + B_operator(B_operator), + BT_operator(BT_operator), + diag_A_inv(diag_A_inv) + {} + void vmult(dealii::LinearAlgebra::distributed::Vector &dst, + const dealii::LinearAlgebra::distributed::Vector &src) const; + private: + const StokesMatrixType &system_matrix; + const BOperatorType &B_operator; + const BTOperatorType &BT_operator; + const dealii::LinearAlgebra::distributed::Vector &diag_A_inv; + }; + template class DiagBFBT: public SchurComplementOperator { diff --git a/source/simulator/solver/stokes_matrix_free_local_smoothing.cc b/source/simulator/solver/stokes_matrix_free_local_smoothing.cc index 12329a49afd..1ddc76d7490 100644 --- a/source/simulator/solver/stokes_matrix_free_local_smoothing.cc +++ b/source/simulator/solver/stokes_matrix_free_local_smoothing.cc @@ -46,6 +46,33 @@ namespace aspect { namespace internal{ + + + template + void BC_invBT_Operator::vmult(dealii::LinearAlgebra::distributed::Vector &dst, + const dealii::LinearAlgebra::distributed::Vector &src) const + { + dealii::LinearAlgebra::distributed::BlockVector block_src; + dealii::LinearAlgebra::distributed::BlockVector block_dst; + system_matrix.initialize_dof_vector(block_src); + system_matrix.initialize_dof_vector(block_dst); + + block_src.block(1)=src; + block_src.block(0)=0; + block_dst=0; + BT_operator.vmult(block_dst,block_src); + + block_dst.block(0).scale(diag_A_inv); + + block_src.block(0)=block_dst.block(0); + block_src.block(1)=0; + block_dst=0; + B_operator.vmult(block_dst,block_src); + dst=block_dst.block(1); + } + + + template DiagBFBT::DiagBFBT(const PreconditionerMp &mp_preconditioner, const double solver_tolerance, From 7e902213e15186b2cfbfaf392540f3e66d98623e Mon Sep 17 00:00:00 2001 From: Quang Hoang Date: Wed, 11 Mar 2026 20:49:00 -0400 Subject: [PATCH 16/37] vmult for diagbfbt --- .../solver/block_stokes_preconditioner.h | 70 +++++----- .../stokes_matrix_free_local_smoothing.cc | 120 ++++++++++++------ 2 files changed, 123 insertions(+), 67 deletions(-) diff --git a/include/aspect/simulator/solver/block_stokes_preconditioner.h b/include/aspect/simulator/solver/block_stokes_preconditioner.h index ff0e5720b5e..ae98a784f84 100644 --- a/include/aspect/simulator/solver/block_stokes_preconditioner.h +++ b/include/aspect/simulator/solver/block_stokes_preconditioner.h @@ -403,36 +403,46 @@ namespace aspect const dealii::LinearAlgebra::distributed::Vector &diag_A_inv; }; - template - class DiagBFBT: public SchurComplementOperator - { - public: - /** - * Constructor. - * @param laplace_preconditioner The preconditioner for @p pressure_laplace_matrix - * @param solver_tolerance The relative solver tolerance for the inner solve - * @param diag_A Diagonal of A used in the BFBT preconditioner. - * @param system_matrix Sparse block matrix storing the Stokes system of the form - * [A B^T - * B 0]. - */ - DiagBFBT(const PreconditionerMp &laplace_preconditioner, - const double solver_tolerance, - const VectorType &diag_A, - const StokesMatrixType &system_matrix); - - void vmult(VectorType &dst, - const VectorType &src) const override; - - unsigned int n_iterations() const override; - - private: - mutable unsigned int n_iterations_; - const PreconditionerMp &laplace_preconditioner; - const double solver_tolerance; - const LinearAlgebra::Vector &diag_A; - const StokesMatrixType &system_matrix; - }; + template +class DiagBFBT: public SchurComplementOperator +{ + public: + /** + * Constructor. + * @param mp_preconditioner The preconditioner for BC^{-1}B^T operator. + * @param solver_tolerance The relative solver tolerance for the inner solve. + * @param diag_A Diagonal of A used in the BFBT preconditioner. + * @param system_matrix Sparse block matrix storing the Stokes system of the form + * [A B^T + * B 0]. + * @param A_operator The velocity block operator. + * @param B_operator The B block operator. + * @param BT_operator The B^T block operator. + */ + DiagBFBT(const PreconditionerMp &mp_preconditioner, + const double solver_tolerance, + const dealii::LinearAlgebra::distributed::Vector &diag_A, + const StokesMatrixType &system_matrix, + const AOperatorType &A_operator, + const BOperatorType &B_operator, + const BTOperatorType &BT_operator); + + void vmult(VectorType &dst, + const VectorType &src) const override; + + unsigned int n_iterations() const override; + + private: + mutable unsigned int n_iterations_; + const PreconditionerMp &mp_preconditioner; + const double solver_tolerance; + const dealii::LinearAlgebra::distributed::Vector &diag_A; + dealii::LinearAlgebra::distributed::Vector diag_A_inv; + const StokesMatrixType &system_matrix; + const AOperatorType &A_operator; + const BOperatorType &B_operator; + const BTOperatorType &BT_operator; +}; } diff --git a/source/simulator/solver/stokes_matrix_free_local_smoothing.cc b/source/simulator/solver/stokes_matrix_free_local_smoothing.cc index 1ddc76d7490..97800ebfcec 100644 --- a/source/simulator/solver/stokes_matrix_free_local_smoothing.cc +++ b/source/simulator/solver/stokes_matrix_free_local_smoothing.cc @@ -73,47 +73,93 @@ namespace aspect - template - DiagBFBT::DiagBFBT(const PreconditionerMp &mp_preconditioner, - const double solver_tolerance, - const VectorType &diag_A, - const StokesMatrixType & system_matrix): - n_iterations_(0), - laplace_preconditioner(laplace_preconditioner), - solver_tolerance(solver_tolerance), - diag_A(diag_A), - system_matrix(system_matrix) - - {} - - - - template - unsigned int DiagBFBT::n_iterations() const{ + template +DiagBFBT::DiagBFBT( + const PreconditionerMp &mp_preconditioner, + const double solver_tolerance, + const dealii::LinearAlgebra::distributed::Vector &diag_A, + const StokesMatrixType &system_matrix, + const AOperatorType &A_operator, + const BOperatorType &B_operator, + const BTOperatorType &BT_operator) + : n_iterations_(0), + mp_preconditioner(mp_preconditioner), + solver_tolerance(solver_tolerance), + diag_A(diag_A), + diag_A_inv(diag_A), + system_matrix(system_matrix), + A_operator(A_operator), + B_operator(B_operator), + BT_operator(BT_operator) +{ + for (auto &el : diag_A_inv) + el = 1.0 / el; +} + + +template +void DiagBFBT::vmult( + VectorType &dst, const VectorType &src) const +{ + try + { + BC_invBT_Operator + Op_BC_invBT(system_matrix, B_operator, BT_operator, diag_A_inv); + + VectorType ptmp, ptmp2; + ptmp.reinit(src); + ptmp2.reinit(src); + + SolverControl solver_control(5000, src.l2_norm() * solver_tolerance, false, true); + SolverCG solver(solver_control); + + ptmp = 0; + solver.solve(Op_BC_invBT, ptmp, src, mp_preconditioner); + n_iterations_ += solver_control.last_step(); + + { + dealii::LinearAlgebra::distributed::BlockVector block_src, block_dst; + system_matrix.initialize_dof_vector(block_src); + system_matrix.initialize_dof_vector(block_dst); + + block_src.block(1) = ptmp; + block_src.block(0) = 0; + block_dst = 0; + BT_operator.vmult(block_dst, block_src); + + block_dst.block(0).scale(diag_A_inv); + + A_operator.vmult(block_src.block(0), block_dst.block(0)); + + block_src.block(0).scale(diag_A_inv); + + block_src.block(1) = 0; + block_dst = 0; + B_operator.vmult(block_dst, block_src); + ptmp2 = block_dst.block(1); + } + + dst = 0; + solver.solve(Op_BC_invBT, dst, ptmp2, mp_preconditioner); + n_iterations_ += solver_control.last_step(); + } + catch (const std::exception &exc) + { + Utilities::throw_linear_solver_failure_exception("iterative (DiagBFBT) solver", + "DiagBFBT::vmult", + std::vector {}, + exc, + src.get_mpi_communicator()); + } +} + + template + unsigned int DiagBFBT::n_iterations() const{ return n_iterations_; } - // template - // void DiagBFBT::vmult(VectorType &dst, const VectorType &src){ - // SolverControl solver_control(1000, src.l2_norm() * solver_tolerance); - // PrimitiveVectorMemory mem; - // SolverCG solver(solver_control, mem); //These vector types might be wrong - // try{ - // VectorType utmp; - // utmp.reinit(diag_A); - // VectorType ptmp; - // ptmp.reinit(src); - // VectorType wtmpt; - // wtmp.reinit(diag_A); - // { - // SolverControl solver_control(5000,1e-6*src.l2_norm(),false,true); - // SolverCG solver(solver_control); - // const auto Op_A=LinearOperator(system_matrix.block(0,0)); - - // } - // } - // } + } From 7391a9fa12a2b0c8eda2e32c0577bcb381136373 Mon Sep 17 00:00:00 2001 From: Quang Hoang Date: Wed, 11 Mar 2026 23:32:42 -0400 Subject: [PATCH 17/37] implement diag bfbt --- .../solver/block_stokes_preconditioner.h | 249 +++++++------- .../simulator/solver/matrix_free_operators.h | 2 +- .../stokes_matrix_free_local_smoothing.h | 2 +- source/simulator/solver.cc | 34 +- .../simulator/solver/matrix_free_operators.cc | 22 +- .../stokes_matrix_free_local_smoothing.cc | 304 +++++++++++------- 6 files changed, 327 insertions(+), 286 deletions(-) diff --git a/include/aspect/simulator/solver/block_stokes_preconditioner.h b/include/aspect/simulator/solver/block_stokes_preconditioner.h index ae98a784f84..ea06fa8aa7a 100644 --- a/include/aspect/simulator/solver/block_stokes_preconditioner.h +++ b/include/aspect/simulator/solver/block_stokes_preconditioner.h @@ -159,91 +159,82 @@ namespace aspect { return n_iterations_; } - /** - * Implement the block Schur preconditioner - * (A B^T; 0 S)^{-1}. - */ - template + * Base class for Schur Complement operators. + */ + template + class SchurComplementOperator + { + public: + virtual ~SchurComplementOperator() = default; + + virtual void vmult(VectorType &dst, + const VectorType &src) const=0; + virtual unsigned int n_iterations() const=0; + + }; + + + + template class BlockSchurPreconditioner : public #if DEAL_II_VERSION_GTE(9,7,0) EnableObserverPointer #else Subscriptor #endif - { public: - /** - * @brief Constructor - * @param A_inverse_operator Approximation of the inverse of the velocity block. - * @param S_inverse_operator Approximation for the inverse Schur complement. - * @param BT_operator Operator for the B^T block of the Stokes system. - */ - BlockSchurPreconditioner ( - const AInvOperator &A_inverse_operator, - const SInvOperator &S_inverse_operator, - const BTOperator &BT_operator); + BlockSchurPreconditioner( + const AInvOperator &A_inverse_operator, + const SchurComplementOperator &S_inverse_operator, + const BTOperator &BT_operator); - /** - * Matrix vector product with this preconditioner object. - */ - void vmult (VectorType &dst, - const VectorType &src) const; + void vmult(VectorType &dst, + const VectorType &src) const; private: - /** - * References to the various operators this preconditioner works with. - */ - - const AInvOperator &A_inverse_operator; - const SInvOperator &S_inverse_operator; - const BTOperator &BT_operator; - mutable VectorType tmp; + const AInvOperator &A_inverse_operator; + const SchurComplementOperator &S_inverse_operator; + const BTOperator &BT_operator; + mutable VectorType tmp; }; - - template - BlockSchurPreconditioner:: - BlockSchurPreconditioner ( - const AInvOperator &A_inverse_operator, - const SInvOperator &S_inverse_operator, - const BTOperator &BT_operator) + template + BlockSchurPreconditioner:: + BlockSchurPreconditioner( + const AInvOperator &A_inverse_operator, + const SchurComplementOperator &S_inverse_operator, + const BTOperator &BT_operator) : - A_inverse_operator (A_inverse_operator), - S_inverse_operator (S_inverse_operator), - BT_operator (BT_operator) + A_inverse_operator(A_inverse_operator), + S_inverse_operator(S_inverse_operator), + BT_operator(BT_operator) {} - template + + + template void - BlockSchurPreconditioner:: - vmult (VectorType &dst, - const VectorType &src) const + BlockSchurPreconditioner:: + vmult(VectorType &dst, + const VectorType &src) const { if (tmp.size() == 0) tmp.reinit(src); + dst = 0.0; - dst=0.0; - - // first apply the Schur Complement inverse operator. { - S_inverse_operator.vmult(dst.block(1),src.block(1)); + S_inverse_operator.vmult(dst.block(1), src.block(1)); dst.block(1) *= -1.0; } - - // apply the top right block: B^T or J^{up} { - // For matrix-free usage, the BT operator needs to operate on the whole vector, - // but the matrix-based version is a SparseMatrix::vmult() that requires passing the - // individual blocks. if constexpr (std::is_same_v>) BT_operator.vmult(tmp, dst); - else BT_operator.vmult(tmp.block(0), dst.block(1)); @@ -254,27 +245,18 @@ namespace aspect A_inverse_operator.vmult(dst.block(0), tmp.block(0)); } - - template - class SchurApproximation + class SchurApproximation : public SchurComplementOperator { public: - /** - * @brief Constructor - * @param schur_preconditioner Preconditioner for the Schur Complement. - * @param stokes_matrix Stokes system. This is not necessarily a matrix. - * @param Schur_complement_block Operator for the Schur complement. - * @param do_solve_Schur_complement Flag that determines whether to do a full solve with the Schur complement or a v-cycle. - * @param Schur_complement_tolerance Tolerance in case a full solve for the Schur complement is used. - */ SchurApproximation(const OperatorType &schur_preconditioner, const StokesMatrixType &stokes_matrix, const SchurComplementMatrixType &Schur_complement_block, const bool do_solve_Schur_complement, const double Schur_complement_tolerance); - void vmult( VectorType &dst, const VectorType &src) const; - unsigned int n_iterations() const; + + void vmult(VectorType &dst, const VectorType &src) const override; + unsigned int n_iterations() const override; private: const OperatorType &schur_preconditioner; @@ -283,8 +265,6 @@ namespace aspect const bool do_solve_Schur_complement; const double Schur_complement_tolerance; mutable unsigned int n_iterations_Schur_complement_; - - }; @@ -304,10 +284,10 @@ namespace aspect n_iterations_Schur_complement_(0) {} - - + + template @@ -366,83 +346,76 @@ namespace aspect - /** - * Base class for Schur Complement operators. - */ - template - class SchurComplementOperator - { - public: - virtual ~SchurComplementOperator() = default; - - virtual void vmult(VectorType &dst, - const VectorType &src) const=0; - virtual unsigned int n_iterations() const=0; - }; template - class BC_invBT_Operator{ + class BC_invBT_Operator + { public: - BC_invBT_Operator(const StokesMatrixType &system_matrix, - const BOperatorType &B_operator, - const BTOperatorType &BT_operator, - const dealii::LinearAlgebra::distributed::Vector &diag_A_inv): - system_matrix(system_matrix), - B_operator(B_operator), - BT_operator(BT_operator), - diag_A_inv(diag_A_inv) - {} - void vmult(dealii::LinearAlgebra::distributed::Vector &dst, - const dealii::LinearAlgebra::distributed::Vector &src) const; + BC_invBT_Operator(const StokesMatrixType &system_matrix, + const BOperatorType &B_operator, + const BTOperatorType &BT_operator, + const dealii::LinearAlgebra::distributed::Vector &diag_A_inv): + system_matrix(system_matrix), + B_operator(B_operator), + BT_operator(BT_operator), + diag_A_inv(diag_A_inv) + {} + void vmult(dealii::LinearAlgebra::distributed::Vector &dst, + const dealii::LinearAlgebra::distributed::Vector &src) const; private: - const StokesMatrixType &system_matrix; - const BOperatorType &B_operator; - const BTOperatorType &BT_operator; - const dealii::LinearAlgebra::distributed::Vector &diag_A_inv; + const StokesMatrixType &system_matrix; + const BOperatorType &B_operator; + const BTOperatorType &BT_operator; + const dealii::LinearAlgebra::distributed::Vector &diag_A_inv; }; template -class DiagBFBT: public SchurComplementOperator -{ - public: - /** - * Constructor. - * @param mp_preconditioner The preconditioner for BC^{-1}B^T operator. - * @param solver_tolerance The relative solver tolerance for the inner solve. - * @param diag_A Diagonal of A used in the BFBT preconditioner. - * @param system_matrix Sparse block matrix storing the Stokes system of the form - * [A B^T - * B 0]. - * @param A_operator The velocity block operator. - * @param B_operator The B block operator. - * @param BT_operator The B^T block operator. - */ - DiagBFBT(const PreconditionerMp &mp_preconditioner, - const double solver_tolerance, - const dealii::LinearAlgebra::distributed::Vector &diag_A, - const StokesMatrixType &system_matrix, - const AOperatorType &A_operator, - const BOperatorType &B_operator, - const BTOperatorType &BT_operator); - - void vmult(VectorType &dst, - const VectorType &src) const override; - - unsigned int n_iterations() const override; - - private: - mutable unsigned int n_iterations_; - const PreconditionerMp &mp_preconditioner; - const double solver_tolerance; - const dealii::LinearAlgebra::distributed::Vector &diag_A; - dealii::LinearAlgebra::distributed::Vector diag_A_inv; - const StokesMatrixType &system_matrix; - const AOperatorType &A_operator; - const BOperatorType &B_operator; - const BTOperatorType &BT_operator; -}; + class DiagBFBT: public SchurComplementOperator + { + public: + /** + * Constructor. + * @param mp_preconditioner The preconditioner for the BC^{-1}B^T operator + * used in the inner CG solves. + * @param do_solve_schur_complement Full solve with Schur complement or just vmult. + * @param solver_tolerance The relative solver tolerance for the inner CG solves + * with BC^{-1}B^T. + * @param diag_A_inv Diagonal of A used as C and D in the diag A-BFBT preconditioner. + * @param system_matrix The Stokes operator of the form + * [A B^T + * B 0]. + * @param A_operator The velocity block operator. + * @param B_operator The B block operator. + * @param BT_operator The B^T block operator. + */ + DiagBFBT(const PreconditionerMp &mp_preconditioner, + const bool do_solve_schur_complement, + const double solver_tolerance, + const dealii::LinearAlgebra::distributed::Vector &diag_A_inv, + const StokesMatrixType &system_matrix, + const AOperatorType &A_operator, + const BOperatorType &B_operator, + const BTOperatorType &BT_operator); + + void vmult(VectorType &dst, + const VectorType &src) const override; + + unsigned int n_iterations() const override; + + private: + mutable unsigned int n_iterations_; + const PreconditionerMp &mp_preconditioner; + const bool do_solve_schur_complement; + const double solver_tolerance; + const dealii::LinearAlgebra::distributed::Vector &diag_A_inv; + const StokesMatrixType &system_matrix; + const AOperatorType &A_operator; + const BOperatorType &B_operator; + const BTOperatorType &BT_operator; + }; + } diff --git a/include/aspect/simulator/solver/matrix_free_operators.h b/include/aspect/simulator/solver/matrix_free_operators.h index a3ad3db14c2..dcafd7d83ef 100644 --- a/include/aspect/simulator/solver/matrix_free_operators.h +++ b/include/aspect/simulator/solver/matrix_free_operators.h @@ -340,7 +340,7 @@ namespace aspect */ const OperatorCellData *cell_data; }; - /** + /** * Operator for the B^T block. */ template diff --git a/include/aspect/simulator/solver/stokes_matrix_free_local_smoothing.h b/include/aspect/simulator/solver/stokes_matrix_free_local_smoothing.h index d0c5d1ae886..34a8ada3489 100644 --- a/include/aspect/simulator/solver/stokes_matrix_free_local_smoothing.h +++ b/include/aspect/simulator/solver/stokes_matrix_free_local_smoothing.h @@ -212,7 +212,7 @@ namespace aspect using SchurComplementMatrixType = MatrixFreeStokesOperators::MassMatrixOperator; using ABlockMatrixType = MatrixFreeStokesOperators::ABlockOperator; using BTBlockOperatorType = MatrixFreeStokesOperators::BTBlockOperator; - using BBlockOperatorType = MatrixFreeStokesOperators::BBlockOperator; + using BBlockOperatorType = MatrixFreeStokesOperators::BBlockOperator; using GMGSchurComplementMatrixType = MatrixFreeStokesOperators::MassMatrixOperator; using GMGABlockMatrixType = MatrixFreeStokesOperators::ABlockOperator; diff --git a/source/simulator/solver.cc b/source/simulator/solver.cc index c5112bc3451..60b95c262b5 100644 --- a/source/simulator/solver.cc +++ b/source/simulator/solver.cc @@ -169,7 +169,7 @@ namespace aspect - + /** * This class approximates the Schur Complement inverse operator @@ -213,7 +213,7 @@ namespace aspect }; - + template WeightedBFBT::WeightedBFBT( @@ -793,12 +793,15 @@ namespace aspect /* do_solve_A = */ false, stokes_A_block_is_symmetric(), parameters.linear_solver_A_block_tolerance); - const internal::BlockSchurPreconditioner, - internal::SchurComplementOperator, LinearAlgebra::SparseMatrix, LinearAlgebra::BlockVector> - preconditioner_cheap ( - inverse_velocity_block_cheap, - *schur, - system_matrix.block(0,1)); + const internal::BlockSchurPreconditioner< + internal::InverseVelocityBlock, + LinearAlgebra::SparseMatrix, + LinearAlgebra::BlockVector, + LinearAlgebra::Vector> + 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( @@ -807,12 +810,15 @@ namespace aspect /* do_solve_A = */ true, stokes_A_block_is_symmetric(), parameters.linear_solver_A_block_tolerance); - const internal::BlockSchurPreconditioner, - internal::SchurComplementOperator, LinearAlgebra::SparseMatrix, LinearAlgebra::BlockVector> - preconditioner_expensive ( - inverse_velocity_block_expensive, - *schur, - system_matrix.block(0,1)); + + const internal::BlockSchurPreconditioner< + internal::InverseVelocityBlock, + LinearAlgebra::SparseMatrix, + LinearAlgebra::BlockVector, + LinearAlgebra::Vector> preconditioner_expensive( + inverse_velocity_block_expensive, + *schur, + system_matrix.block(0,1)); // step 1a: try if the simple and fast solver // succeeds in n_cheap_stokes_solver_steps steps or less. try diff --git a/source/simulator/solver/matrix_free_operators.cc b/source/simulator/solver/matrix_free_operators.cc index 196720d0343..32953789de6 100644 --- a/source/simulator/solver/matrix_free_operators.cc +++ b/source/simulator/solver/matrix_free_operators.cc @@ -521,16 +521,16 @@ namespace aspect - /** + /** *Operator for B block. - */ + */ template MatrixFreeStokesOperators::BBlockOperator::BBlockOperator (): - MatrixFreeOperators::Base>() + MatrixFreeOperators::Base>() {} template - void + void MatrixFreeStokesOperators::BBlockOperator::clear () { this->cell_data = nullptr; @@ -538,7 +538,7 @@ namespace aspect } - + template void MatrixFreeStokesOperators::BBlockOperator:: @@ -549,7 +549,7 @@ namespace aspect template - void + void MatrixFreeStokesOperators::BBlockOperator ::compute_diagonal () { @@ -560,7 +560,7 @@ namespace aspect template - void + void MatrixFreeStokesOperators::BBlockOperator ::local_apply(const dealii::MatrixFree &data, dealii::LinearAlgebra::distributed::BlockVector &dst, @@ -580,9 +580,9 @@ namespace aspect p_eval.reinit(cell); - + SymmetricTensor<2,dim,VectorizedArray> sym_grad_u; - + for (const unsigned int q : u_eval.quadrature_point_indices()) { sym_grad_u = u_eval.get_symmetric_gradient(q); @@ -593,9 +593,9 @@ namespace aspect p_eval.submit_value(pressure_terms, q); - + } - + p_eval.integrate_scatter(EvaluationFlags::values, dst.block(1)); } diff --git a/source/simulator/solver/stokes_matrix_free_local_smoothing.cc b/source/simulator/solver/stokes_matrix_free_local_smoothing.cc index 97800ebfcec..402a03fcb08 100644 --- a/source/simulator/solver/stokes_matrix_free_local_smoothing.cc +++ b/source/simulator/solver/stokes_matrix_free_local_smoothing.cc @@ -45,125 +45,150 @@ namespace aspect { - namespace internal{ + namespace internal + { template void BC_invBT_Operator::vmult(dealii::LinearAlgebra::distributed::Vector &dst, - const dealii::LinearAlgebra::distributed::Vector &src) const - { - dealii::LinearAlgebra::distributed::BlockVector block_src; - dealii::LinearAlgebra::distributed::BlockVector block_dst; - system_matrix.initialize_dof_vector(block_src); - system_matrix.initialize_dof_vector(block_dst); - - block_src.block(1)=src; - block_src.block(0)=0; - block_dst=0; - BT_operator.vmult(block_dst,block_src); - - block_dst.block(0).scale(diag_A_inv); - - block_src.block(0)=block_dst.block(0); - block_src.block(1)=0; - block_dst=0; - B_operator.vmult(block_dst,block_src); - dst=block_dst.block(1); - } + const dealii::LinearAlgebra::distributed::Vector &src) const + { + dealii::LinearAlgebra::distributed::BlockVector block_src; + dealii::LinearAlgebra::distributed::BlockVector block_dst; + block_src.reinit(2); + block_dst.reinit(2); + + system_matrix.initialize_dof_vector(block_src); + system_matrix.initialize_dof_vector(block_dst); + + block_src.block(1)=src; + block_src.block(0)=0; + block_dst=0; + BT_operator.vmult(block_dst,block_src); + + block_dst.block(0).scale(diag_A_inv); + + block_src.block(0)=block_dst.block(0); + block_src.block(1)=0; + block_dst=0; + B_operator.vmult(block_dst,block_src); + dst=block_dst.block(1); + } template -DiagBFBT::DiagBFBT( - const PreconditionerMp &mp_preconditioner, - const double solver_tolerance, - const dealii::LinearAlgebra::distributed::Vector &diag_A, - const StokesMatrixType &system_matrix, - const AOperatorType &A_operator, - const BOperatorType &B_operator, - const BTOperatorType &BT_operator) - : n_iterations_(0), - mp_preconditioner(mp_preconditioner), - solver_tolerance(solver_tolerance), - diag_A(diag_A), - diag_A_inv(diag_A), - system_matrix(system_matrix), - A_operator(A_operator), - B_operator(B_operator), - BT_operator(BT_operator) -{ - for (auto &el : diag_A_inv) - el = 1.0 / el; -} + DiagBFBT::DiagBFBT( + const PreconditionerMp &mp_preconditioner, + const bool do_solve_schur_complement, + const double solver_tolerance, + const dealii::LinearAlgebra::distributed::Vector &diag_A_inv, + const StokesMatrixType &system_matrix, + const AOperatorType &A_operator, + const BOperatorType &B_operator, + const BTOperatorType &BT_operator) + : n_iterations_(0), + mp_preconditioner(mp_preconditioner), + do_solve_schur_complement(do_solve_schur_complement), + solver_tolerance(solver_tolerance), + diag_A_inv(diag_A_inv), + system_matrix(system_matrix), + A_operator(A_operator), + B_operator(B_operator), + BT_operator(BT_operator) + { + } -template -void DiagBFBT::vmult( - VectorType &dst, const VectorType &src) const -{ - try + template + void DiagBFBT::vmult( + VectorType &dst, const VectorType &src) const { - BC_invBT_Operator - Op_BC_invBT(system_matrix, B_operator, BT_operator, diag_A_inv); - - VectorType ptmp, ptmp2; - ptmp.reinit(src); - ptmp2.reinit(src); + try + { + BC_invBT_Operator + Op_BC_invBT(system_matrix, B_operator, BT_operator, diag_A_inv); - SolverControl solver_control(5000, src.l2_norm() * solver_tolerance, false, true); - SolverCG solver(solver_control); + VectorType ptmp; + VectorType ptmp2; + ptmp.reinit(src); + ptmp2.reinit(src); - ptmp = 0; - solver.solve(Op_BC_invBT, ptmp, src, mp_preconditioner); - n_iterations_ += solver_control.last_step(); + if (!do_solve_schur_complement) + { + Op_BC_invBT.vmult(ptmp, src); + n_iterations_ += 1; + } + else + { + SolverControl solver_control(5000, src.l2_norm() * solver_tolerance, false, true); + SolverCG solver(solver_control); + ptmp = 0; + solver.solve(Op_BC_invBT, ptmp, src, mp_preconditioner); + n_iterations_ += solver_control.last_step(); + } - { - dealii::LinearAlgebra::distributed::BlockVector block_src, block_dst; - system_matrix.initialize_dof_vector(block_src); - system_matrix.initialize_dof_vector(block_dst); + { + dealii::LinearAlgebra::distributed::BlockVector block_src; + dealii::LinearAlgebra::distributed::BlockVector block_dst; + block_src.reinit(2); + block_dst.reinit(2); + system_matrix.initialize_dof_vector(block_src); + system_matrix.initialize_dof_vector(block_dst); - block_src.block(1) = ptmp; - block_src.block(0) = 0; - block_dst = 0; - BT_operator.vmult(block_dst, block_src); + block_src.block(1) = ptmp; + block_src.block(0) = 0; + block_dst = 0; + BT_operator.vmult(block_dst, block_src); - block_dst.block(0).scale(diag_A_inv); + block_dst.block(0).scale(diag_A_inv); - A_operator.vmult(block_src.block(0), block_dst.block(0)); + A_operator.vmult(block_src.block(0), block_dst.block(0)); - block_src.block(0).scale(diag_A_inv); + block_src.block(0).scale(diag_A_inv); - block_src.block(1) = 0; - block_dst = 0; - B_operator.vmult(block_dst, block_src); - ptmp2 = block_dst.block(1); - } + block_src.block(1) = 0; + block_dst = 0; + B_operator.vmult(block_dst, block_src); + ptmp2 = block_dst.block(1); + } - dst = 0; - solver.solve(Op_BC_invBT, dst, ptmp2, mp_preconditioner); - n_iterations_ += solver_control.last_step(); - } - catch (const std::exception &exc) - { - Utilities::throw_linear_solver_failure_exception("iterative (DiagBFBT) solver", - "DiagBFBT::vmult", - std::vector {}, - exc, - src.get_mpi_communicator()); + if (!do_solve_schur_complement) + { + Op_BC_invBT.vmult(dst, ptmp2); + n_iterations_ += 1; + } + else + { + SolverControl solver_control(5000, src.l2_norm() * solver_tolerance, false, true); + SolverCG solver(solver_control); + dst = 0; + solver.solve(Op_BC_invBT, dst, ptmp2, mp_preconditioner); + n_iterations_ += solver_control.last_step(); + } + } + catch (const std::exception &exc) + { + Utilities::throw_linear_solver_failure_exception("iterative (DiagBFBT) solver", + "DiagBFBT::vmult", + std::vector {}, + exc, + src.get_mpi_communicator()); + } } -} template - unsigned int DiagBFBT::n_iterations() const{ + unsigned int DiagBFBT::n_iterations() const + { return n_iterations_; } - - - + + + } - + template void StokesMatrixFreeHandlerLocalSmoothingImplementation::declare_parameters(ParameterHandler &prm) @@ -475,6 +500,7 @@ void DiagBFBTget_parameters().n_expensive_stokes_solver_steps > 0) @@ -1335,6 +1361,12 @@ void DiagBFBT>; + using BlockSchurPreconditionerType = internal::BlockSchurPreconditioner< + internal::InverseVelocityBlock, + BTBlockOperatorType, + dealii::LinearAlgebra::distributed::BlockVector, + VectorType>; + internal::InverseVelocityBlock inverse_velocity_block_cheap( A_block_matrix, prec_A, @@ -1342,41 +1374,71 @@ void DiagBFBTget_parameters().linear_solver_A_block_tolerance); - internal::InverseVelocityBlock inverse_velocity_block_expensive( + internal::InverseVelocityBlock inverse_velocity_block_expensive( A_block_matrix, prec_A, /* do_solve_A = */ true, sim.stokes_A_block_is_symmetric(), this->get_parameters().linear_solver_A_block_tolerance); - using SchurApproximationType = internal::SchurApproximation; - internal::SchurApproximation schur_approximation_cheap( - prec_Schur, - stokes_matrix, - Schur_complement_block_matrix, - /*do_solve_Schur*/ false, - this->get_parameters().linear_solver_S_block_tolerance); + dealii::LinearAlgebra::distributed::Vector diag_A; + std::unique_ptr> schur_approximation_cheap; + std::unique_ptr> schur_approximation_expensive; - internal::SchurApproximation schur_approximation_expensive( - prec_Schur, - stokes_matrix, - Schur_complement_block_matrix, - /*do_solve_Schur*/ true, - this->get_parameters().linear_solver_S_block_tolerance); + if (this->get_parameters().use_bfbt) + { + A_block_matrix.compute_diagonal(); + const dealii::LinearAlgebra::distributed::Vector &diag_A_inv = + A_block_matrix.get_matrix_diagonal_inverse()->get_vector(); + + using DiagBFBTType = internal::DiagBFBT; + schur_approximation_cheap = std::make_unique( + prec_Schur, + /*do_solve_schur_complement*/ true, + this->get_parameters().linear_solver_S_block_tolerance, + diag_A_inv, + stokes_matrix, + A_block_matrix, + B_block, + BT_block); //hack - the full schur solves do not seem to converge + + schur_approximation_expensive = std::make_unique( + prec_Schur, + /*do_solve_schur_complement*/ true, + this->get_parameters().linear_solver_S_block_tolerance, + diag_A_inv, + stokes_matrix, + A_block_matrix, + B_block, + BT_block); + } + else + { + using SchurApproximationType = internal::SchurApproximation; + schur_approximation_cheap = std::make_unique( + prec_Schur, + stokes_matrix, + Schur_complement_block_matrix, + /*do_solve_Schur*/ false, + this->get_parameters().linear_solver_S_block_tolerance); + + schur_approximation_expensive = std::make_unique( + prec_Schur, + stokes_matrix, + Schur_complement_block_matrix, + /*do_solve_Schur*/ true, + this->get_parameters().linear_solver_S_block_tolerance); + } - const internal::BlockSchurPreconditioner, - SchurApproximationType,BTBlockOperatorType, dealii::LinearAlgebra::distributed::BlockVector> - preconditioner_cheap ( - inverse_velocity_block_cheap, - schur_approximation_cheap, - BT_block); + const BlockSchurPreconditionerType preconditioner_cheap( + inverse_velocity_block_cheap, + *schur_approximation_cheap, + BT_block); - const internal::BlockSchurPreconditioner, - SchurApproximationType, BTBlockOperatorType, dealii::LinearAlgebra::distributed::BlockVector> - preconditioner_expensive ( - inverse_velocity_block_expensive, - schur_approximation_expensive, - BT_block); + const BlockSchurPreconditionerType preconditioner_expensive( + inverse_velocity_block_expensive, + *schur_approximation_expensive, + BT_block); @@ -1630,7 +1692,7 @@ void DiagBFBTget_signals().post_stokes_solver(sim, - schur_approximation_cheap.n_iterations() + schur_approximation_expensive.n_iterations(), + schur_approximation_cheap->n_iterations() + schur_approximation_expensive->n_iterations(), inverse_velocity_block_cheap.n_iterations() + inverse_velocity_block_expensive.n_iterations(), solver_control_cheap, solver_control_expensive); @@ -1652,7 +1714,7 @@ void DiagBFBTget_signals().post_stokes_solver(sim, - schur_approximation_cheap.n_iterations() + schur_approximation_expensive.n_iterations(), + schur_approximation_cheap->n_iterations() + schur_approximation_expensive->n_iterations(), inverse_velocity_block_cheap.n_iterations() + inverse_velocity_block_expensive.n_iterations(), solver_control_cheap, solver_control_expensive); @@ -1681,9 +1743,9 @@ void DiagBFBTget_pcout() << " Schur complement preconditioner: " << schur_approximation_cheap.n_iterations() + this->get_pcout() << " Schur complement preconditioner: " << schur_approximation_cheap->n_iterations() << '+' - << schur_approximation_expensive.n_iterations() + << schur_approximation_expensive->n_iterations() << " iterations." << std::endl; this->get_pcout() << " A block preconditioner: " << inverse_velocity_block_cheap.n_iterations() << '+' From 1bc2c5bdfd2eb5cd8902e05b218adf3ac92c0b6d Mon Sep 17 00:00:00 2001 From: Quang Hoang Date: Wed, 11 Mar 2026 23:49:22 -0400 Subject: [PATCH 18/37] fix cheap solve and l2 norm --- .../simulator/solver/stokes_matrix_free_local_smoothing.cc | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/source/simulator/solver/stokes_matrix_free_local_smoothing.cc b/source/simulator/solver/stokes_matrix_free_local_smoothing.cc index 402a03fcb08..730e8019acd 100644 --- a/source/simulator/solver/stokes_matrix_free_local_smoothing.cc +++ b/source/simulator/solver/stokes_matrix_free_local_smoothing.cc @@ -116,7 +116,7 @@ namespace aspect if (!do_solve_schur_complement) { - Op_BC_invBT.vmult(ptmp, src); + mp_preconditioner.vmult(ptmp, src); n_iterations_ += 1; } else @@ -155,12 +155,12 @@ namespace aspect if (!do_solve_schur_complement) { - Op_BC_invBT.vmult(dst, ptmp2); + mp_preconditioner.vmult(dst, ptmp2); n_iterations_ += 1; } else { - SolverControl solver_control(5000, src.l2_norm() * solver_tolerance, false, true); + SolverControl solver_control(5000, ptmp.l2_norm() * solver_tolerance, false, true); SolverCG solver(solver_control); dst = 0; solver.solve(Op_BC_invBT, dst, ptmp2, mp_preconditioner); From 606074c40afce85cd879191f03c42812e5418ac3 Mon Sep 17 00:00:00 2001 From: Quang Hoang Date: Thu, 12 Mar 2026 14:15:24 -0400 Subject: [PATCH 19/37] fix l2_norm --- source/simulator/solver/stokes_matrix_free_local_smoothing.cc | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/source/simulator/solver/stokes_matrix_free_local_smoothing.cc b/source/simulator/solver/stokes_matrix_free_local_smoothing.cc index 730e8019acd..3928fbce2d6 100644 --- a/source/simulator/solver/stokes_matrix_free_local_smoothing.cc +++ b/source/simulator/solver/stokes_matrix_free_local_smoothing.cc @@ -160,7 +160,7 @@ namespace aspect } else { - SolverControl solver_control(5000, ptmp.l2_norm() * solver_tolerance, false, true); + SolverControl solver_control(5000, ptmp2.l2_norm() * solver_tolerance, false, true); SolverCG solver(solver_control); dst = 0; solver.solve(Op_BC_invBT, dst, ptmp2, mp_preconditioner); From 4150c60e296c13497acff8375e1ad06ef4f0e522 Mon Sep 17 00:00:00 2001 From: Quang Hoang Date: Thu, 12 Mar 2026 15:07:57 -0400 Subject: [PATCH 20/37] use primitive memory --- .../solver/stokes_matrix_free_local_smoothing.cc | 8 +++++--- 1 file changed, 5 insertions(+), 3 deletions(-) diff --git a/source/simulator/solver/stokes_matrix_free_local_smoothing.cc b/source/simulator/solver/stokes_matrix_free_local_smoothing.cc index 3928fbce2d6..2f80e3eaba3 100644 --- a/source/simulator/solver/stokes_matrix_free_local_smoothing.cc +++ b/source/simulator/solver/stokes_matrix_free_local_smoothing.cc @@ -113,6 +113,7 @@ namespace aspect VectorType ptmp2; ptmp.reinit(src); ptmp2.reinit(src); + PrimitiveVectorMemory mem; if (!do_solve_schur_complement) { @@ -121,8 +122,9 @@ namespace aspect } else { + SolverControl solver_control(5000, src.l2_norm() * solver_tolerance, false, true); - SolverCG solver(solver_control); + SolverCG solver(solver_control,mem); ptmp = 0; solver.solve(Op_BC_invBT, ptmp, src, mp_preconditioner); n_iterations_ += solver_control.last_step(); @@ -161,7 +163,7 @@ namespace aspect else { SolverControl solver_control(5000, ptmp2.l2_norm() * solver_tolerance, false, true); - SolverCG solver(solver_control); + SolverCG solver(solver_control,mem); dst = 0; solver.solve(Op_BC_invBT, dst, ptmp2, mp_preconditioner); n_iterations_ += solver_control.last_step(); @@ -1400,7 +1402,7 @@ namespace aspect stokes_matrix, A_block_matrix, B_block, - BT_block); //hack - the full schur solves do not seem to converge + BT_block); //hack - the vmults do not seem to vonverge. schur_approximation_expensive = std::make_unique( prec_Schur, From 57f2eb57257f4a697cd01efbe207b33c2bf2b92e Mon Sep 17 00:00:00 2001 From: Quang Hoang Date: Sat, 21 Mar 2026 20:06:55 -0400 Subject: [PATCH 21/37] add run all ffile --- benchmarks/nsinker/diag_A_run_all.sh | 40 ++++++++++++++++++++++++++++ 1 file changed, 40 insertions(+) create mode 100644 benchmarks/nsinker/diag_A_run_all.sh diff --git a/benchmarks/nsinker/diag_A_run_all.sh b/benchmarks/nsinker/diag_A_run_all.sh new file mode 100644 index 00000000000..f6a8fdcd05e --- /dev/null +++ b/benchmarks/nsinker/diag_A_run_all.sh @@ -0,0 +1,40 @@ +#!/bin/bash + +for averaging in none; do # arithmetic/geometric/harmonic average + for nsinkers in 4 8 16; do + for viscosity in 1e1 1e3 1e5 1e7 ; do + for refinement in 4 5 6 ; do + echo "subsection Material model" > current.prm + echo " set Material averaging = $averaging" >> current.prm + echo " subsection NSinker" >> current.prm + echo " set Number of sinkers = $nsinkers" >> current.prm + echo " set Dynamic viscosity ratio = $viscosity" >> current.prm + echo " end" >> current.prm + echo "end" >> current.prm + echo "subsection Solver parameters">>current.prm + echo "subsection Stokes solver parameters">>current.prm + echo "set Stokes solver type= block AMG">>current.prm + echo "set Use full A block as preconditioner=true">>current.prm + echo "set Use weighted BFBT for Schur complement=true">>current.prm + echo "set Number of cheap Stokes solver steps=500">>current.prm + echo "set Maximum number of expensive Stokes solver steps=0">>current.prm + echo "end">>current.prm + echo "subsection AMG parameters">>current.prm + echo "set AMG aggregation threshold = 0.02">>current.prm + echo "end">>current.prm + echo "end">>current.prm + echo "subsection Mesh refinement" >> current.prm + echo " set Initial global refinement = $refinement" >> current.prm + echo "end" >> current.prm + echo "subsection Discretization">>current.prm + echo "set Use locally conservative discretization=false">>current.prm + echo "end">>current.prm + + current_model="averaging${averaging}_nsinkers${nsinkers}_viscosity${viscosity}_refinement${refinement}" + echo "set Output directory = /home/qxhoang/aspect/benchmarks/nsinker/diag_A_mem_fix/output-${current_model}" >> current.prm + echo "Starting ${current_model}" + cat nsinker.prm current.prm | mpirun -np 32 ./aspect-release -- + done + done + done +done From 430924ece25cd20922fde5f8aacd2d29c5c5b9c2 Mon Sep 17 00:00:00 2001 From: Quang Hoang Date: Mon, 23 Mar 2026 16:19:51 -0400 Subject: [PATCH 22/37] fix cheap and expensive solves and change pressure scaling --- benchmarks/nsinker/nsinker.cc | 4 ++-- .../stokes_matrix_free_local_smoothing.cc | 24 +++++-------------- 2 files changed, 8 insertions(+), 20 deletions(-) diff --git a/benchmarks/nsinker/nsinker.cc b/benchmarks/nsinker/nsinker.cc index 55f42ba729d..1080f6631e8 100644 --- a/benchmarks/nsinker/nsinker.cc +++ b/benchmarks/nsinker/nsinker.cc @@ -317,11 +317,11 @@ namespace aspect // Change pressure scaling to 1.0: -double pressure_scaling_signal(const double /*pressure_scaling*/, +double pressure_scaling_signal(const double pressure_scaling, const double /*reference_viscosity*/, const double /*length_scale*/) { - return 1.0; + return pressure_scaling; } template diff --git a/source/simulator/solver/stokes_matrix_free_local_smoothing.cc b/source/simulator/solver/stokes_matrix_free_local_smoothing.cc index 2f80e3eaba3..f54bef3b842 100644 --- a/source/simulator/solver/stokes_matrix_free_local_smoothing.cc +++ b/source/simulator/solver/stokes_matrix_free_local_smoothing.cc @@ -115,20 +115,15 @@ namespace aspect ptmp2.reinit(src); PrimitiveVectorMemory mem; - if (!do_solve_schur_complement) - { - mp_preconditioner.vmult(ptmp, src); - n_iterations_ += 1; - } - else - { + + SolverControl solver_control(5000, src.l2_norm() * solver_tolerance, false, true); SolverCG solver(solver_control,mem); ptmp = 0; solver.solve(Op_BC_invBT, ptmp, src, mp_preconditioner); n_iterations_ += solver_control.last_step(); - } + { dealii::LinearAlgebra::distributed::BlockVector block_src; @@ -155,20 +150,13 @@ namespace aspect ptmp2 = block_dst.block(1); } - if (!do_solve_schur_complement) - { - mp_preconditioner.vmult(dst, ptmp2); - n_iterations_ += 1; - } - else - { - SolverControl solver_control(5000, ptmp2.l2_norm() * solver_tolerance, false, true); - SolverCG solver(solver_control,mem); + + solver_control.set_tolerance(1e-6*ptmp.l2_norm()); dst = 0; solver.solve(Op_BC_invBT, dst, ptmp2, mp_preconditioner); n_iterations_ += solver_control.last_step(); - } } + catch (const std::exception &exc) { Utilities::throw_linear_solver_failure_exception("iterative (DiagBFBT) solver", From 534a0e92354686e8b2b8d41b8533272152af9c9a Mon Sep 17 00:00:00 2001 From: Quang Hoang Date: Mon, 23 Mar 2026 17:31:33 -0400 Subject: [PATCH 23/37] set cell data A --- .../simulator/solver/stokes_matrix_free_local_smoothing.cc | 5 ++--- 1 file changed, 2 insertions(+), 3 deletions(-) diff --git a/source/simulator/solver/stokes_matrix_free_local_smoothing.cc b/source/simulator/solver/stokes_matrix_free_local_smoothing.cc index f54bef3b842..f105da76adf 100644 --- a/source/simulator/solver/stokes_matrix_free_local_smoothing.cc +++ b/source/simulator/solver/stokes_matrix_free_local_smoothing.cc @@ -493,11 +493,10 @@ namespace aspect B_block.set_cell_data(active_cell_data); BT_block.set_cell_data(active_cell_data); - if (this->get_parameters().n_expensive_stokes_solver_steps > 0) - { + A_block_matrix.set_cell_data(active_cell_data); Schur_complement_block_matrix.set_cell_data(active_cell_data); - } + const unsigned int n_levels = this->get_triangulation().n_global_levels(); level_cell_data.resize(0,n_levels-1); From c07debb6d10258b54438c765c911f10fec3f8007 Mon Sep 17 00:00:00 2001 From: Quang Hoang Date: Fri, 3 Apr 2026 02:33:25 -0400 Subject: [PATCH 24/37] use diag preconditioner --- .../solver/block_stokes_preconditioner.h | 11 +++++++++++ .../solver/stokes_matrix_free_local_smoothing.cc | 16 +++++++++++----- 2 files changed, 22 insertions(+), 5 deletions(-) diff --git a/include/aspect/simulator/solver/block_stokes_preconditioner.h b/include/aspect/simulator/solver/block_stokes_preconditioner.h index ea06fa8aa7a..efd3dec1932 100644 --- a/include/aspect/simulator/solver/block_stokes_preconditioner.h +++ b/include/aspect/simulator/solver/block_stokes_preconditioner.h @@ -403,6 +403,17 @@ namespace aspect const VectorType &src) const override; unsigned int n_iterations() const override; + class DiagonalPreconditioner{ + public: + DiagonalPreconditioner(const dealii::DiagonalMatrix &diagonal_matrix): + diagonal_matrix(diagonal_matrix) + {} + void vmult(VectorType &dst, const VectorType &src) const{ + diagonal_matrix.vmult(dst,src); + } + private: + const dealii::DiagonalMatrix &diagonal_matrix; + }; private: mutable unsigned int n_iterations_; diff --git a/source/simulator/solver/stokes_matrix_free_local_smoothing.cc b/source/simulator/solver/stokes_matrix_free_local_smoothing.cc index f105da76adf..2ad2351d2ec 100644 --- a/source/simulator/solver/stokes_matrix_free_local_smoothing.cc +++ b/source/simulator/solver/stokes_matrix_free_local_smoothing.cc @@ -151,7 +151,7 @@ namespace aspect } - solver_control.set_tolerance(1e-6*ptmp.l2_norm()); + solver_control.set_tolerance(1e-6*ptmp2.l2_norm()); dst = 0; solver.solve(Op_BC_invBT, dst, ptmp2, mp_preconditioner); n_iterations_ += solver_control.last_step(); @@ -1377,12 +1377,18 @@ namespace aspect if (this->get_parameters().use_bfbt) { A_block_matrix.compute_diagonal(); + Schur_complement_block_matrix.compute_diagonal(); + const dealii::LinearAlgebra::distributed::Vector &diag_A_inv = A_block_matrix.get_matrix_diagonal_inverse()->get_vector(); - - using DiagBFBTType = internal::DiagBFBT; + const dealii::DiagonalMatrix &diag_mp=*Schur_complement_block_matrix.get_matrix_diagonal_inverse(); + using DiagPreconditionerType=typename internal::DiagBFBT< StokesMatrixType,ABlockMatrixType, + BBlockOperatorType,BTBlockOperatorType,VectorType,GMGPreconditioner>::DiagonalPreconditioner; + using DiagBFBTType = internal::DiagBFBT; + DiagPreconditionerType diag_prec(diag_mp); + schur_approximation_cheap = std::make_unique( - prec_Schur, + diag_prec, /*do_solve_schur_complement*/ true, this->get_parameters().linear_solver_S_block_tolerance, diag_A_inv, @@ -1392,7 +1398,7 @@ namespace aspect BT_block); //hack - the vmults do not seem to vonverge. schur_approximation_expensive = std::make_unique( - prec_Schur, + diag_prec, /*do_solve_schur_complement*/ true, this->get_parameters().linear_solver_S_block_tolerance, diag_A_inv, From 7fead789571c07d0d89af350c57c8fe1f15fe1e5 Mon Sep 17 00:00:00 2001 From: Quang Hoang Date: Fri, 3 Apr 2026 20:52:41 -0400 Subject: [PATCH 25/37] set up tests --- benchmarks/nsinker/diag_A_run_all.sh | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) mode change 100644 => 100755 benchmarks/nsinker/diag_A_run_all.sh diff --git a/benchmarks/nsinker/diag_A_run_all.sh b/benchmarks/nsinker/diag_A_run_all.sh old mode 100644 new mode 100755 index f6a8fdcd05e..4f6d0fea7df --- a/benchmarks/nsinker/diag_A_run_all.sh +++ b/benchmarks/nsinker/diag_A_run_all.sh @@ -13,7 +13,7 @@ for averaging in none; do # arithmetic/geometric/harmonic average echo "end" >> current.prm echo "subsection Solver parameters">>current.prm echo "subsection Stokes solver parameters">>current.prm - echo "set Stokes solver type= block AMG">>current.prm + echo "set Stokes solver type= block GMG">>current.prm echo "set Use full A block as preconditioner=true">>current.prm echo "set Use weighted BFBT for Schur complement=true">>current.prm echo "set Number of cheap Stokes solver steps=500">>current.prm @@ -31,7 +31,7 @@ for averaging in none; do # arithmetic/geometric/harmonic average echo "end">>current.prm current_model="averaging${averaging}_nsinkers${nsinkers}_viscosity${viscosity}_refinement${refinement}" - echo "set Output directory = /home/qxhoang/aspect/benchmarks/nsinker/diag_A_mem_fix/output-${current_model}" >> current.prm + echo "set Output directory = home/quang-hoang/Documents/aspect/benchmarks/nsinker/diag_A_diag_mp/output-${current_model}" >> current.prm echo "Starting ${current_model}" cat nsinker.prm current.prm | mpirun -np 32 ./aspect-release -- done From cbec3ac77c5f760ba0d627142a2c3231a1a9637d Mon Sep 17 00:00:00 2001 From: Quang Hoang Date: Mon, 6 Apr 2026 14:20:19 -0400 Subject: [PATCH 26/37] diagonal wrapper class --- .../simulator/solver/block_stokes_preconditioner.h | 13 +------------ .../solver/stokes_matrix_free_local_smoothing.cc | 9 +++------ 2 files changed, 4 insertions(+), 18 deletions(-) diff --git a/include/aspect/simulator/solver/block_stokes_preconditioner.h b/include/aspect/simulator/solver/block_stokes_preconditioner.h index efd3dec1932..a6b9965a0bf 100644 --- a/include/aspect/simulator/solver/block_stokes_preconditioner.h +++ b/include/aspect/simulator/solver/block_stokes_preconditioner.h @@ -403,18 +403,7 @@ namespace aspect const VectorType &src) const override; unsigned int n_iterations() const override; - class DiagonalPreconditioner{ - public: - DiagonalPreconditioner(const dealii::DiagonalMatrix &diagonal_matrix): - diagonal_matrix(diagonal_matrix) - {} - void vmult(VectorType &dst, const VectorType &src) const{ - diagonal_matrix.vmult(dst,src); - } - private: - const dealii::DiagonalMatrix &diagonal_matrix; - }; - + private: mutable unsigned int n_iterations_; const PreconditionerMp &mp_preconditioner; diff --git a/source/simulator/solver/stokes_matrix_free_local_smoothing.cc b/source/simulator/solver/stokes_matrix_free_local_smoothing.cc index 2ad2351d2ec..286ddca0c23 100644 --- a/source/simulator/solver/stokes_matrix_free_local_smoothing.cc +++ b/source/simulator/solver/stokes_matrix_free_local_smoothing.cc @@ -1382,13 +1382,10 @@ namespace aspect const dealii::LinearAlgebra::distributed::Vector &diag_A_inv = A_block_matrix.get_matrix_diagonal_inverse()->get_vector(); const dealii::DiagonalMatrix &diag_mp=*Schur_complement_block_matrix.get_matrix_diagonal_inverse(); - using DiagPreconditionerType=typename internal::DiagBFBT< StokesMatrixType,ABlockMatrixType, - BBlockOperatorType,BTBlockOperatorType,VectorType,GMGPreconditioner>::DiagonalPreconditioner; - using DiagBFBTType = internal::DiagBFBT; - DiagPreconditionerType diag_prec(diag_mp); + using DiagBFBTType = internal::DiagBFBT>; schur_approximation_cheap = std::make_unique( - diag_prec, + diag_mp, /*do_solve_schur_complement*/ true, this->get_parameters().linear_solver_S_block_tolerance, diag_A_inv, @@ -1398,7 +1395,7 @@ namespace aspect BT_block); //hack - the vmults do not seem to vonverge. schur_approximation_expensive = std::make_unique( - diag_prec, + diag_mp, /*do_solve_schur_complement*/ true, this->get_parameters().linear_solver_S_block_tolerance, diag_A_inv, From 810a9be41d9e49d83f007e3bae98a6e7ebccad42 Mon Sep 17 00:00:00 2001 From: Quang Hoang Date: Fri, 10 Apr 2026 23:08:00 -0400 Subject: [PATCH 27/37] remove mean value --- .../solver/stokes_matrix_free_local_smoothing.cc | 12 +++++++----- 1 file changed, 7 insertions(+), 5 deletions(-) diff --git a/source/simulator/solver/stokes_matrix_free_local_smoothing.cc b/source/simulator/solver/stokes_matrix_free_local_smoothing.cc index 286ddca0c23..282c7725f74 100644 --- a/source/simulator/solver/stokes_matrix_free_local_smoothing.cc +++ b/source/simulator/solver/stokes_matrix_free_local_smoothing.cc @@ -117,8 +117,9 @@ namespace aspect - - SolverControl solver_control(5000, src.l2_norm() * solver_tolerance, false, true); + VectorType rhs1=src; //nullspace removal + rhs1.add(-rhs1.mean_value()); + SolverControl solver_control(5000, rhs1.l2_norm() * solver_tolerance, false, true); SolverCG solver(solver_control,mem); ptmp = 0; solver.solve(Op_BC_invBT, ptmp, src, mp_preconditioner); @@ -150,10 +151,11 @@ namespace aspect ptmp2 = block_dst.block(1); } - - solver_control.set_tolerance(1e-6*ptmp2.l2_norm()); + VectorType rhs2=ptmp2; + rhs2.add(-rhs2.mean_value()); + solver_control.set_tolerance(1e-6*rhs2.l2_norm()); dst = 0; - solver.solve(Op_BC_invBT, dst, ptmp2, mp_preconditioner); + solver.solve(Op_BC_invBT, dst, rhs2, mp_preconditioner); n_iterations_ += solver_control.last_step(); } From 06999aa98ebd00e48b865bc0c2e725b24b2fea59 Mon Sep 17 00:00:00 2001 From: Quang Hoang Date: Fri, 10 Apr 2026 23:39:38 -0400 Subject: [PATCH 28/37] operator ns removal --- .../solver/block_stokes_preconditioner.h | 2 +- .../stokes_matrix_free_local_smoothing.cc | 91 +++++++++++++++---- 2 files changed, 74 insertions(+), 19 deletions(-) diff --git a/include/aspect/simulator/solver/block_stokes_preconditioner.h b/include/aspect/simulator/solver/block_stokes_preconditioner.h index a6b9965a0bf..ea06fa8aa7a 100644 --- a/include/aspect/simulator/solver/block_stokes_preconditioner.h +++ b/include/aspect/simulator/solver/block_stokes_preconditioner.h @@ -403,7 +403,7 @@ namespace aspect const VectorType &src) const override; unsigned int n_iterations() const override; - + private: mutable unsigned int n_iterations_; const PreconditionerMp &mp_preconditioner; diff --git a/source/simulator/solver/stokes_matrix_free_local_smoothing.cc b/source/simulator/solver/stokes_matrix_free_local_smoothing.cc index 282c7725f74..cb84b4148de 100644 --- a/source/simulator/solver/stokes_matrix_free_local_smoothing.cc +++ b/source/simulator/solver/stokes_matrix_free_local_smoothing.cc @@ -47,6 +47,33 @@ namespace aspect namespace internal { + template + struct Nullspace + { + std::vector basis; + }; + + + + + template + LinearOperator remove_mean_value(LinearOperator &exemplar) + { + 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.add(-dest.mean_value()); + }; + return return_op; + } template @@ -108,6 +135,34 @@ namespace aspect { BC_invBT_Operator Op_BC_invBT(system_matrix, B_operator, BT_operator, diag_A_inv); + dealii::LinearOperator op_BC_invBT; + op_BC_invBT.reinit_range_vector=[&](VectorType &v, bool) + { + v.reinit(src); + }; + op_BC_invBT.reinit_domain_vector=[&](VectorType &v, bool) + { + v.reinit(src); + }; + op_BC_invBT.vmult=[&](VectorType &dst, const VectorType &src) + { + Op_BC_invBT.vmult(dst,src); + }; + + dealii::LinearOperator op_preconditioner; + op_preconditioner.reinit_range_vector = [&](VectorType &v, bool) + { + v.reinit(src); + }; + op_preconditioner.reinit_domain_vector = [&](VectorType &v, bool) + { + v.reinit(src); + }; + op_preconditioner.vmult = [&](VectorType &dst, const VectorType &src) + { + mp_preconditioner.vmult(dst, src); + }; + auto rmv=remove_mean_value<>(op_BC_invBT); VectorType ptmp; VectorType ptmp2; @@ -115,16 +170,16 @@ namespace aspect ptmp2.reinit(src); PrimitiveVectorMemory mem; - - + + VectorType rhs1=src; //nullspace removal rhs1.add(-rhs1.mean_value()); - SolverControl solver_control(5000, rhs1.l2_norm() * solver_tolerance, false, true); - SolverCG solver(solver_control,mem); - ptmp = 0; - solver.solve(Op_BC_invBT, ptmp, src, mp_preconditioner); - n_iterations_ += solver_control.last_step(); - + SolverControl solver_control(5000, rhs1.l2_norm() * solver_tolerance, false, true); + SolverCG solver(solver_control,mem); + ptmp = 0; + solver.solve(rmv*op_BC_invBT, ptmp, rhs1, rmv*op_preconditioner); + n_iterations_ += solver_control.last_step(); + { dealii::LinearAlgebra::distributed::BlockVector block_src; @@ -153,12 +208,12 @@ namespace aspect VectorType rhs2=ptmp2; rhs2.add(-rhs2.mean_value()); - solver_control.set_tolerance(1e-6*rhs2.l2_norm()); - dst = 0; - solver.solve(Op_BC_invBT, dst, rhs2, mp_preconditioner); - n_iterations_ += solver_control.last_step(); + solver_control.set_tolerance(1e-6*rhs2.l2_norm()); + dst = 0; + solver.solve(rmv*op_BC_invBT, dst, rhs2, rmv*op_preconditioner); + n_iterations_ += solver_control.last_step(); } - + catch (const std::exception &exc) { Utilities::throw_linear_solver_failure_exception("iterative (DiagBFBT) solver", @@ -495,10 +550,10 @@ namespace aspect B_block.set_cell_data(active_cell_data); BT_block.set_cell_data(active_cell_data); - - A_block_matrix.set_cell_data(active_cell_data); - Schur_complement_block_matrix.set_cell_data(active_cell_data); - + + A_block_matrix.set_cell_data(active_cell_data); + Schur_complement_block_matrix.set_cell_data(active_cell_data); + const unsigned int n_levels = this->get_triangulation().n_global_levels(); level_cell_data.resize(0,n_levels-1); @@ -1385,7 +1440,7 @@ namespace aspect A_block_matrix.get_matrix_diagonal_inverse()->get_vector(); const dealii::DiagonalMatrix &diag_mp=*Schur_complement_block_matrix.get_matrix_diagonal_inverse(); using DiagBFBTType = internal::DiagBFBT>; - + schur_approximation_cheap = std::make_unique( diag_mp, /*do_solve_schur_complement*/ true, From a072e8e5e199cbfc9f2748aaacb4611cef7cec87 Mon Sep 17 00:00:00 2001 From: Quang Hoang Date: Wed, 15 Apr 2026 17:11:20 -0400 Subject: [PATCH 29/37] debug code --- .../solver/stokes_matrix_free_local_smoothing.cc | 10 ++++++++++ 1 file changed, 10 insertions(+) diff --git a/source/simulator/solver/stokes_matrix_free_local_smoothing.cc b/source/simulator/solver/stokes_matrix_free_local_smoothing.cc index cb84b4148de..88988f06d2a 100644 --- a/source/simulator/solver/stokes_matrix_free_local_smoothing.cc +++ b/source/simulator/solver/stokes_matrix_free_local_smoothing.cc @@ -124,6 +124,16 @@ namespace aspect B_operator(B_operator), BT_operator(BT_operator) { + std::cout << "diag_A_inv: "; + for (auto i : diag_A_inv.locally_owned_elements()) + std::cout << diag_A_inv[i] << " "; + std::cout << std::endl; + + const auto &diag_mp_vec = mp_preconditioner.get_vector(); + std::cout << "diag_mp_inv: "; + for (auto i : diag_mp_vec.locally_owned_elements()) + std::cout << diag_mp_vec[i] << " "; + std::cout << std::endl; } From d3c8c779db691df41d16f59d043e2c96186626eb Mon Sep 17 00:00:00 2001 From: Quang Hoang Date: Fri, 24 Apr 2026 22:16:00 -0400 Subject: [PATCH 30/37] debugging --- benchmarks/nsinker/nsinker.cc | 4 +-- benchmarks/nsinker/nsinker.prm | 9 +++--- .../stokes_matrix_free_local_smoothing.cc | 32 +++++++++++++++++++ 3 files changed, 39 insertions(+), 6 deletions(-) diff --git a/benchmarks/nsinker/nsinker.cc b/benchmarks/nsinker/nsinker.cc index 1080f6631e8..55f42ba729d 100644 --- a/benchmarks/nsinker/nsinker.cc +++ b/benchmarks/nsinker/nsinker.cc @@ -317,11 +317,11 @@ namespace aspect // Change pressure scaling to 1.0: -double pressure_scaling_signal(const double pressure_scaling, +double pressure_scaling_signal(const double /*pressure_scaling*/, const double /*reference_viscosity*/, const double /*length_scale*/) { - return pressure_scaling; + return 1.0; } template diff --git a/benchmarks/nsinker/nsinker.prm b/benchmarks/nsinker/nsinker.prm index a2e194c0c5c..a96755fe547 100644 --- a/benchmarks/nsinker/nsinker.prm +++ b/benchmarks/nsinker/nsinker.prm @@ -16,14 +16,15 @@ set Use years instead of seconds = false # Follow as closely as possible the parameters from Rudi et al. (2017) subsection Solver parameters subsection Stokes solver parameters - set Stokes solver type = block AMG + set Stokes solver type = block GMG set Use full A block as preconditioner = true set Number of cheap Stokes solver steps = 500 - set Maximum number of expensive Stokes solver steps = 1000 + set Maximum number of expensive Stokes solver steps = 0 set Linear solver tolerance = 1e-6 set GMRES solver restart length = 100 - set Use weighted BFBT for Schur complement = false + set Use weighted BFBT for Schur complement = true set Krylov method for cheap solver steps = GMRES + set Linear solver S block tolerance = 1e-2 end subsection AMG parameters @@ -72,7 +73,7 @@ subsection Initial temperature model end subsection Mesh refinement - set Initial global refinement = 4 + set Initial global refinement = 1 set Initial adaptive refinement = 0 end diff --git a/source/simulator/solver/stokes_matrix_free_local_smoothing.cc b/source/simulator/solver/stokes_matrix_free_local_smoothing.cc index 88988f06d2a..aaeb6ab82cb 100644 --- a/source/simulator/solver/stokes_matrix_free_local_smoothing.cc +++ b/source/simulator/solver/stokes_matrix_free_local_smoothing.cc @@ -184,12 +184,27 @@ namespace aspect VectorType rhs1=src; //nullspace removal rhs1.add(-rhs1.mean_value()); + + //DEBUG CODE + std::cout<<"rhs1 before solve:"; + for(auto i : rhs1.locally_owned_elements()){ + std::cout< solver(solver_control,mem); ptmp = 0; solver.solve(rmv*op_BC_invBT, ptmp, rhs1, rmv*op_preconditioner); n_iterations_ += solver_control.last_step(); + //DEBUG CODE + std::cout<<"ptmp after first solve: "; + for(auto i: ptmp.locally_owned_elements()){ + std::cout< block_src; @@ -218,10 +233,27 @@ namespace aspect VectorType rhs2=ptmp2; rhs2.add(-rhs2.mean_value()); + + //DEBUG CODE + + std::cout<<"rhs 2 before solve: "; + for(auto i: rhs2.locally_owned_elements()){ + std::cout< Date: Fri, 24 Apr 2026 22:22:46 -0400 Subject: [PATCH 31/37] precondition with identity --- .../stokes_matrix_free_local_smoothing.cc | 56 ++++++++++--------- 1 file changed, 29 insertions(+), 27 deletions(-) diff --git a/source/simulator/solver/stokes_matrix_free_local_smoothing.cc b/source/simulator/solver/stokes_matrix_free_local_smoothing.cc index aaeb6ab82cb..bc2394d6059 100644 --- a/source/simulator/solver/stokes_matrix_free_local_smoothing.cc +++ b/source/simulator/solver/stokes_matrix_free_local_smoothing.cc @@ -185,25 +185,27 @@ namespace aspect VectorType rhs1=src; //nullspace removal rhs1.add(-rhs1.mean_value()); - //DEBUG CODE - std::cout<<"rhs1 before solve:"; - for(auto i : rhs1.locally_owned_elements()){ - std::cout< solver(solver_control,mem); ptmp = 0; - solver.solve(rmv*op_BC_invBT, ptmp, rhs1, rmv*op_preconditioner); + solver.solve(rmv*op_BC_invBT, ptmp, rhs1, identity); n_iterations_ += solver_control.last_step(); - //DEBUG CODE - std::cout<<"ptmp after first solve: "; - for(auto i: ptmp.locally_owned_elements()){ - std::cout< Date: Fri, 24 Apr 2026 22:30:12 -0400 Subject: [PATCH 32/37] comment out debug code --- .../stokes_matrix_free_local_smoothing.cc | 20 +++++++++---------- 1 file changed, 10 insertions(+), 10 deletions(-) diff --git a/source/simulator/solver/stokes_matrix_free_local_smoothing.cc b/source/simulator/solver/stokes_matrix_free_local_smoothing.cc index bc2394d6059..43803ac259c 100644 --- a/source/simulator/solver/stokes_matrix_free_local_smoothing.cc +++ b/source/simulator/solver/stokes_matrix_free_local_smoothing.cc @@ -124,16 +124,16 @@ namespace aspect B_operator(B_operator), BT_operator(BT_operator) { - std::cout << "diag_A_inv: "; - for (auto i : diag_A_inv.locally_owned_elements()) - std::cout << diag_A_inv[i] << " "; - std::cout << std::endl; - - const auto &diag_mp_vec = mp_preconditioner.get_vector(); - std::cout << "diag_mp_inv: "; - for (auto i : diag_mp_vec.locally_owned_elements()) - std::cout << diag_mp_vec[i] << " "; - std::cout << std::endl; + // std::cout << "diag_A_inv: "; + // for (auto i : diag_A_inv.locally_owned_elements()) + // std::cout << diag_A_inv[i] << " "; + // std::cout << std::endl; + + // const auto &diag_mp_vec = mp_preconditioner.get_vector(); + // std::cout << "diag_mp_inv: "; + // for (auto i : diag_mp_vec.locally_owned_elements()) + // std::cout << diag_mp_vec[i] << " "; + // std::cout << std::endl; } From aff3b1957904720d5e9ced149ac2595b33643291 Mon Sep 17 00:00:00 2001 From: Quang Hoang Date: Tue, 28 Apr 2026 16:32:27 -0400 Subject: [PATCH 33/37] full solve mp --- .../solver/block_stokes_preconditioner.h | 7 ++- .../stokes_matrix_free_local_smoothing.cc | 50 +++++++++++-------- 2 files changed, 33 insertions(+), 24 deletions(-) diff --git a/include/aspect/simulator/solver/block_stokes_preconditioner.h b/include/aspect/simulator/solver/block_stokes_preconditioner.h index ea06fa8aa7a..d6a76e5976d 100644 --- a/include/aspect/simulator/solver/block_stokes_preconditioner.h +++ b/include/aspect/simulator/solver/block_stokes_preconditioner.h @@ -371,7 +371,7 @@ namespace aspect const dealii::LinearAlgebra::distributed::Vector &diag_A_inv; }; - template + template class DiagBFBT: public SchurComplementOperator { public: @@ -389,6 +389,7 @@ namespace aspect * @param A_operator The velocity block operator. * @param B_operator The B block operator. * @param BT_operator The B^T block operator. + * @param mp_matrix Pressure mass matrix used as preconditioner for BC^{-1}B^T. */ DiagBFBT(const PreconditionerMp &mp_preconditioner, const bool do_solve_schur_complement, @@ -397,7 +398,8 @@ namespace aspect const StokesMatrixType &system_matrix, const AOperatorType &A_operator, const BOperatorType &B_operator, - const BTOperatorType &BT_operator); + const BTOperatorType &BT_operator, + const SchurComplementMatrixType &mp_matrix); void vmult(VectorType &dst, const VectorType &src) const override; @@ -414,6 +416,7 @@ namespace aspect const AOperatorType &A_operator; const BOperatorType &B_operator; const BTOperatorType &BT_operator; + const SchurComplementMatrixType &mp_matrix; }; diff --git a/source/simulator/solver/stokes_matrix_free_local_smoothing.cc b/source/simulator/solver/stokes_matrix_free_local_smoothing.cc index 43803ac259c..26daee12c19 100644 --- a/source/simulator/solver/stokes_matrix_free_local_smoothing.cc +++ b/source/simulator/solver/stokes_matrix_free_local_smoothing.cc @@ -104,8 +104,8 @@ namespace aspect - template - DiagBFBT::DiagBFBT( + template + DiagBFBT::DiagBFBT( const PreconditionerMp &mp_preconditioner, const bool do_solve_schur_complement, const double solver_tolerance, @@ -113,7 +113,8 @@ namespace aspect const StokesMatrixType &system_matrix, const AOperatorType &A_operator, const BOperatorType &B_operator, - const BTOperatorType &BT_operator) + const BTOperatorType &BT_operator, + const SchurComplementMatrixType &mp_matrix) : n_iterations_(0), mp_preconditioner(mp_preconditioner), do_solve_schur_complement(do_solve_schur_complement), @@ -122,7 +123,8 @@ namespace aspect system_matrix(system_matrix), A_operator(A_operator), B_operator(B_operator), - BT_operator(BT_operator) + BT_operator(BT_operator), + mp_matrix(mp_matrix) { // std::cout << "diag_A_inv: "; // for (auto i : diag_A_inv.locally_owned_elements()) @@ -137,8 +139,8 @@ namespace aspect } - template - void DiagBFBT::vmult( + template + void DiagBFBT::vmult( VectorType &dst, const VectorType &src) const { try @@ -159,18 +161,21 @@ namespace aspect Op_BC_invBT.vmult(dst,src); }; - dealii::LinearOperator op_preconditioner; - op_preconditioner.reinit_range_vector = [&](VectorType &v, bool) - { + dealii::LinearOperator op_mp_preconditioner; + op_mp_preconditioner.reinit_range_vector=[&](VectorType &v, bool){ v.reinit(src); }; - op_preconditioner.reinit_domain_vector = [&](VectorType &v, bool) - { + op_mp_preconditioner.reinit_domain_vector=[&](VectorType &v, bool){ v.reinit(src); }; - op_preconditioner.vmult = [&](VectorType &dst, const VectorType &src) - { - mp_preconditioner.vmult(dst, src); + + //precondition with solve + op_mp_preconditioner.vmult=[&](VectorType &dst, const VectorType &src){ + PrimitiveVectorMemory mp_mem; + SolverControl solver_control(1000,src.l2_norm()*1e-6); + SolverCG solver(solver_control,mp_mem); + dst=0.0; + solver.solve(mp_matrix,dst,src,mp_preconditioner); }; auto rmv=remove_mean_value<>(op_BC_invBT); @@ -193,11 +198,10 @@ namespace aspect // std::cout< solver(solver_control,mem); ptmp = 0; - solver.solve(rmv*op_BC_invBT, ptmp, rhs1, identity); + solver.solve(rmv*op_BC_invBT, ptmp, rhs1, op_mp_preconditioner); n_iterations_ += solver_control.last_step(); // //DEBUG CODE @@ -246,7 +250,7 @@ namespace aspect solver_control.set_tolerance(1e-6*rhs2.l2_norm()); dst = 0; - solver.solve(rmv*op_BC_invBT, dst, rhs2, identity); + solver.solve(rmv*op_BC_invBT, dst, rhs2, op_mp_preconditioner); n_iterations_ += solver_control.last_step(); // //DEBUG CODE @@ -268,8 +272,8 @@ namespace aspect } } - template - unsigned int DiagBFBT::n_iterations() const + template + unsigned int DiagBFBT::n_iterations() const { return n_iterations_; } @@ -1483,7 +1487,7 @@ namespace aspect const dealii::LinearAlgebra::distributed::Vector &diag_A_inv = A_block_matrix.get_matrix_diagonal_inverse()->get_vector(); const dealii::DiagonalMatrix &diag_mp=*Schur_complement_block_matrix.get_matrix_diagonal_inverse(); - using DiagBFBTType = internal::DiagBFBT>; + using DiagBFBTType = internal::DiagBFBT>; schur_approximation_cheap = std::make_unique( diag_mp, @@ -1493,7 +1497,8 @@ namespace aspect stokes_matrix, A_block_matrix, B_block, - BT_block); //hack - the vmults do not seem to vonverge. + BT_block, + Schur_complement_block_matrix); //hack - the vmults do not seem to vonverge. schur_approximation_expensive = std::make_unique( diag_mp, @@ -1503,7 +1508,8 @@ namespace aspect stokes_matrix, A_block_matrix, B_block, - BT_block); + BT_block, + Schur_complement_block_matrix); } else { From 4983f64065e8cffa1a28dd397e33906b5e303546 Mon Sep 17 00:00:00 2001 From: Quang Hoang Date: Wed, 29 Apr 2026 15:13:04 -0400 Subject: [PATCH 34/37] vmults --- .../simulator/solver/stokes_matrix_free_local_smoothing.cc | 6 ++++-- 1 file changed, 4 insertions(+), 2 deletions(-) diff --git a/source/simulator/solver/stokes_matrix_free_local_smoothing.cc b/source/simulator/solver/stokes_matrix_free_local_smoothing.cc index 26daee12c19..6408b4cd422 100644 --- a/source/simulator/solver/stokes_matrix_free_local_smoothing.cc +++ b/source/simulator/solver/stokes_matrix_free_local_smoothing.cc @@ -201,7 +201,8 @@ namespace aspect SolverControl solver_control(5000, rhs1.l2_norm() * solver_tolerance, false, true); SolverCG solver(solver_control,mem); ptmp = 0; - solver.solve(rmv*op_BC_invBT, ptmp, rhs1, op_mp_preconditioner); + op_mp_preconditioner.vmult(ptmp,rhs1); + // solver.solve(rmv*op_BC_invBT, ptmp, rhs1, op_mp_preconditioner); n_iterations_ += solver_control.last_step(); // //DEBUG CODE @@ -250,7 +251,8 @@ namespace aspect solver_control.set_tolerance(1e-6*rhs2.l2_norm()); dst = 0; - solver.solve(rmv*op_BC_invBT, dst, rhs2, op_mp_preconditioner); + // solver.solve(rmv*op_BC_invBT, dst, rhs2, op_mp_preconditioner); + op_mp_preconditioner.vmult(dst,rhs2); n_iterations_ += solver_control.last_step(); // //DEBUG CODE From e9c47ca925ca906c113f088f4bfb63a1a0f3ad5a Mon Sep 17 00:00:00 2001 From: Quang Hoang Date: Wed, 29 Apr 2026 15:34:10 -0400 Subject: [PATCH 35/37] use prec_schur --- .../solver/stokes_matrix_free_local_smoothing.cc | 12 +++++------- 1 file changed, 5 insertions(+), 7 deletions(-) diff --git a/source/simulator/solver/stokes_matrix_free_local_smoothing.cc b/source/simulator/solver/stokes_matrix_free_local_smoothing.cc index 6408b4cd422..6ca8aa10a23 100644 --- a/source/simulator/solver/stokes_matrix_free_local_smoothing.cc +++ b/source/simulator/solver/stokes_matrix_free_local_smoothing.cc @@ -201,8 +201,7 @@ namespace aspect SolverControl solver_control(5000, rhs1.l2_norm() * solver_tolerance, false, true); SolverCG solver(solver_control,mem); ptmp = 0; - op_mp_preconditioner.vmult(ptmp,rhs1); - // solver.solve(rmv*op_BC_invBT, ptmp, rhs1, op_mp_preconditioner); + solver.solve(rmv*op_BC_invBT, ptmp, rhs1, op_mp_preconditioner); n_iterations_ += solver_control.last_step(); // //DEBUG CODE @@ -251,8 +250,7 @@ namespace aspect solver_control.set_tolerance(1e-6*rhs2.l2_norm()); dst = 0; - // solver.solve(rmv*op_BC_invBT, dst, rhs2, op_mp_preconditioner); - op_mp_preconditioner.vmult(dst,rhs2); + solver.solve(rmv*op_BC_invBT, dst, rhs2, op_mp_preconditioner); n_iterations_ += solver_control.last_step(); // //DEBUG CODE @@ -1489,10 +1487,10 @@ namespace aspect const dealii::LinearAlgebra::distributed::Vector &diag_A_inv = A_block_matrix.get_matrix_diagonal_inverse()->get_vector(); const dealii::DiagonalMatrix &diag_mp=*Schur_complement_block_matrix.get_matrix_diagonal_inverse(); - using DiagBFBTType = internal::DiagBFBT>; + using DiagBFBTType = internal::DiagBFBT; schur_approximation_cheap = std::make_unique( - diag_mp, + prec_Schur, /*do_solve_schur_complement*/ true, this->get_parameters().linear_solver_S_block_tolerance, diag_A_inv, @@ -1503,7 +1501,7 @@ namespace aspect Schur_complement_block_matrix); //hack - the vmults do not seem to vonverge. schur_approximation_expensive = std::make_unique( - diag_mp, + prec_Schur, /*do_solve_schur_complement*/ true, this->get_parameters().linear_solver_S_block_tolerance, diag_A_inv, From 9585277c360bc7b8acd6740eeb5980d6b72720e9 Mon Sep 17 00:00:00 2001 From: Quang Hoang Date: Tue, 12 May 2026 12:02:53 -0400 Subject: [PATCH 36/37] remove ns removal --- source/simulator/solver.cc | 11 +++++++++++ .../solver/stokes_matrix_free_local_smoothing.cc | 16 ++++++++-------- 2 files changed, 19 insertions(+), 8 deletions(-) diff --git a/source/simulator/solver.cc b/source/simulator/solver.cc index 60b95c262b5..ccc392058e9 100644 --- a/source/simulator/solver.cc +++ b/source/simulator/solver.cc @@ -786,6 +786,17 @@ namespace aspect } + //DEBUG CODE print diag(A) + LinearAlgebra::Vector diag_A_amg(system_matrix.block(0,0).locally_owned_range_indices(),this->mpi_communicator); + std::cout<<"diag_A (matrix based)"; + for(auto i:diag_A_amg.locally_owned_elements()){ + diag_A_amg[i]=system_matrix.block(0,0).diag_element(i); + } + for(auto i:diag_A_amg.locally_owned_elements()){ + std::cout< inverse_velocity_block_cheap( system_matrix.block(velocity_block_index,velocity_block_index), diff --git a/source/simulator/solver/stokes_matrix_free_local_smoothing.cc b/source/simulator/solver/stokes_matrix_free_local_smoothing.cc index 6ca8aa10a23..af923947add 100644 --- a/source/simulator/solver/stokes_matrix_free_local_smoothing.cc +++ b/source/simulator/solver/stokes_matrix_free_local_smoothing.cc @@ -126,10 +126,10 @@ namespace aspect BT_operator(BT_operator), mp_matrix(mp_matrix) { - // std::cout << "diag_A_inv: "; - // for (auto i : diag_A_inv.locally_owned_elements()) - // std::cout << diag_A_inv[i] << " "; - // std::cout << std::endl; + std::cout << "diag_A_inv (matrix free): "; + for (auto i : diag_A_inv.locally_owned_elements()) + std::cout << 1.0/diag_A_inv[i] << " "; + std::cout << std::endl; // const auto &diag_mp_vec = mp_preconditioner.get_vector(); // std::cout << "diag_mp_inv: "; @@ -188,7 +188,7 @@ namespace aspect VectorType rhs1=src; //nullspace removal - rhs1.add(-rhs1.mean_value()); + // rhs1.add(-rhs1.mean_value()); // //DEBUG CODE // std::cout<<"rhs1 before solve:"; @@ -201,7 +201,7 @@ namespace aspect SolverControl solver_control(5000, rhs1.l2_norm() * solver_tolerance, false, true); SolverCG solver(solver_control,mem); ptmp = 0; - solver.solve(rmv*op_BC_invBT, ptmp, rhs1, op_mp_preconditioner); + solver.solve(/*rmv**/op_BC_invBT, ptmp, rhs1, op_mp_preconditioner); n_iterations_ += solver_control.last_step(); // //DEBUG CODE @@ -238,7 +238,7 @@ namespace aspect } VectorType rhs2=ptmp2; - rhs2.add(-rhs2.mean_value()); + // rhs2.add(-rhs2.mean_value()); // //DEBUG CODE @@ -250,7 +250,7 @@ namespace aspect solver_control.set_tolerance(1e-6*rhs2.l2_norm()); dst = 0; - solver.solve(rmv*op_BC_invBT, dst, rhs2, op_mp_preconditioner); + solver.solve(/*rmv**/op_BC_invBT, dst, rhs2, op_mp_preconditioner); n_iterations_ += solver_control.last_step(); // //DEBUG CODE From e64cc70cc81a7dfa9add4176c31ad8fcab2e0f50 Mon Sep 17 00:00:00 2001 From: Quang Hoang Date: Tue, 12 May 2026 12:17:42 -0400 Subject: [PATCH 37/37] remove debug comments --- source/simulator/solver.cc | 16 ++++++++-------- .../solver/stokes_matrix_free_local_smoothing.cc | 8 ++++---- 2 files changed, 12 insertions(+), 12 deletions(-) diff --git a/source/simulator/solver.cc b/source/simulator/solver.cc index ccc392058e9..b69b1eba464 100644 --- a/source/simulator/solver.cc +++ b/source/simulator/solver.cc @@ -788,14 +788,14 @@ namespace aspect //DEBUG CODE print diag(A) LinearAlgebra::Vector diag_A_amg(system_matrix.block(0,0).locally_owned_range_indices(),this->mpi_communicator); - std::cout<<"diag_A (matrix based)"; - for(auto i:diag_A_amg.locally_owned_elements()){ - diag_A_amg[i]=system_matrix.block(0,0).diag_element(i); - } - for(auto i:diag_A_amg.locally_owned_elements()){ - std::cout< inverse_velocity_block_cheap( diff --git a/source/simulator/solver/stokes_matrix_free_local_smoothing.cc b/source/simulator/solver/stokes_matrix_free_local_smoothing.cc index af923947add..6300be8f99d 100644 --- a/source/simulator/solver/stokes_matrix_free_local_smoothing.cc +++ b/source/simulator/solver/stokes_matrix_free_local_smoothing.cc @@ -126,10 +126,10 @@ namespace aspect BT_operator(BT_operator), mp_matrix(mp_matrix) { - std::cout << "diag_A_inv (matrix free): "; - for (auto i : diag_A_inv.locally_owned_elements()) - std::cout << 1.0/diag_A_inv[i] << " "; - std::cout << std::endl; + // std::cout << "diag_A_inv (matrix free): "; + // for (auto i : diag_A_inv.locally_owned_elements()) + // std::cout << 1.0/diag_A_inv[i] << " "; + // std::cout << std::endl; // const auto &diag_mp_vec = mp_preconditioner.get_vector(); // std::cout << "diag_mp_inv: ";