diff --git a/benchmarks/nsinker/diag_A_run_all.sh b/benchmarks/nsinker/diag_A_run_all.sh new file mode 100755 index 00000000000..4f6d0fea7df --- /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 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 + 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/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 + done + done +done 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/include/aspect/simulator/solver/block_stokes_preconditioner.h b/include/aspect/simulator/solver/block_stokes_preconditioner.h index 6f21a785a6e..d6a76e5976d 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_; - - }; @@ -305,6 +285,11 @@ namespace aspect {} + + + + + template void SchurApproximation::vmult(VectorType &dst, const VectorType &src) const @@ -350,6 +335,9 @@ namespace aspect n_iterations_Schur_complement_ += 1; } } + + + template unsigned int SchurApproximation::n_iterations() const { @@ -357,6 +345,82 @@ namespace aspect } + + + + + 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 + { + 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. + * @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, + 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, + const SchurComplementMatrixType &mp_matrix); + + 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; + const SchurComplementMatrixType &mp_matrix; + }; + + + } } diff --git a/include/aspect/simulator/solver/matrix_free_operators.h b/include/aspect/simulator/solver/matrix_free_operators.h index 67bff4ee07b..dcafd7d83ef 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; + }; 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..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,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.cc b/source/simulator/solver.cc index ceba95be9b1..b69b1eba464 100644 --- a/source/simulator/solver.cc +++ b/source/simulator/solver.cc @@ -166,19 +166,10 @@ 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; - }; + + /** * This class approximates the Schur Complement inverse operator @@ -188,7 +179,7 @@ namespace aspect * velocity mass matrix. */ template - class WeightedBFBT: public SchurComplementOperator + class WeightedBFBT: public SchurComplementOperator { public: /** @@ -221,6 +212,9 @@ namespace aspect const LinearAlgebra::BlockSparseMatrix &system_matrix; }; + + + template WeightedBFBT::WeightedBFBT( const LinearAlgebra::SparseMatrix &pressure_laplace_matrix, @@ -309,7 +303,7 @@ namespace aspect * PreconditionerMp passed to the constructor. */ template - class InverseWeightedMassMatrix: public SchurComplementOperator + class InverseWeightedMassMatrix: public SchurComplementOperator { public: /** @@ -578,8 +572,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) @@ -771,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>( @@ -790,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), @@ -797,12 +804,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( @@ -811,12 +821,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 c1f82d75f5a..32953789de6 100644 --- a/source/simulator/solver/matrix_free_operators.cc +++ b/source/simulator/solver/matrix_free_operators.cc @@ -519,6 +519,98 @@ 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:: + 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()); + } + + + + 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); + + + + for (unsigned int cell=cell_range.first; cell> sym_grad_u; + + for (const unsigned int q : u_eval.quadrature_point_indices()) + { + sym_grad_u = u_eval.get_symmetric_gradient(q); + + const VectorizedArray div_u = trace(sym_grad_u); + VectorizedArray pressure_terms = + -cell_data->pressure_scaling * div_u; + p_eval.submit_value(pressure_terms, q); + + + + } + + + p_eval.integrate_scatter(EvaluationFlags::values, dst.block(1)); + } + } + + 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 */ @@ -934,6 +1026,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 9df74f2001b..6300be8f99d 100644 --- a/source/simulator/solver/stokes_matrix_free_local_smoothing.cc +++ b/source/simulator/solver/stokes_matrix_free_local_smoothing.cc @@ -45,6 +45,245 @@ 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 + 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; + 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 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, + const SchurComplementMatrixType &mp_matrix) + : 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), + 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; + + // 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; + + } + + 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); + 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_mp_preconditioner; + op_mp_preconditioner.reinit_range_vector=[&](VectorType &v, bool){ + v.reinit(src); + }; + op_mp_preconditioner.reinit_domain_vector=[&](VectorType &v, bool){ + v.reinit(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); + + VectorType ptmp; + VectorType ptmp2; + ptmp.reinit(src); + ptmp2.reinit(src); + PrimitiveVectorMemory mem; + + + + 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, op_mp_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; + 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_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); + } + + 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< {}, + exc, + src.get_mpi_communicator()); + } + } + + template + unsigned int DiagBFBT::n_iterations() const + { + return n_iterations_; + } + + + + + + } + template void StokesMatrixFreeHandlerLocalSmoothingImplementation::declare_parameters(ParameterHandler &prm) @@ -356,13 +595,13 @@ namespace aspect // Store viscosity tables and other data into the active level matrix-free objects. stokes_matrix.set_cell_data(active_cell_data); + 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); - } + + 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); @@ -1216,6 +1455,12 @@ namespace aspect solver_control_expensive.enable_history_data(); using GMGPreconditioner = PreconditionMG>; + using BlockSchurPreconditionerType = internal::BlockSchurPreconditioner< + internal::InverseVelocityBlock, + BTBlockOperatorType, + dealii::LinearAlgebra::distributed::BlockVector, + VectorType>; + internal::InverseVelocityBlock inverse_velocity_block_cheap( A_block_matrix, prec_A, @@ -1223,41 +1468,76 @@ namespace aspect sim.stokes_A_block_is_symmetric(), this->get_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(); + Schur_complement_block_matrix.compute_diagonal(); + + 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; + + 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, + Schur_complement_block_matrix); //hack - the vmults do not seem to vonverge. + + 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, + Schur_complement_block_matrix); + } + 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); @@ -1511,7 +1791,7 @@ namespace aspect catch (const std::exception &exc) { this->get_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); @@ -1533,7 +1813,7 @@ namespace aspect //signal successful solver this->get_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); @@ -1562,9 +1842,9 @@ namespace aspect if (print_details) { - this->get_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() << '+' @@ -1805,6 +2085,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();