1 // Copyright (c) 2010-2021, Lawrence Livermore National Security, LLC. Produced 2 // at the Lawrence Livermore National Laboratory. All Rights reserved. See files 3 // LICENSE and NOTICE for details. LLNL-CODE-806117. 4 // 5 // This file is part of the MFEM library. For more information and source code 6 // availability visit https://mfem.org. 7 // 8 // MFEM is free software; you can redistribute it and/or modify it under the 9 // terms of the BSD-3 license. We welcome feedback and contributions, see file 10 // CONTRIBUTING.md for details. 11 12 #ifndef MFEM_DIVFREE_SOLVER_HPP 13 #define MFEM_DIVFREE_SOLVER_HPP 14 15 #include "mfem.hpp" 16 #include <memory> 17 #include <vector> 18 19 namespace mfem 20 { 21 namespace blocksolvers 22 { 23 24 /// Parameters for iterative solver 25 struct IterSolveParameters 26 { 27 int print_level = 0; 28 int max_iter = 500; 29 double abs_tol = 1e-12; 30 double rel_tol = 1e-9; 31 }; 32 33 /// Parameters for the divergence free solver 34 struct DFSParameters : IterSolveParameters 35 { 36 /** There are three components in the solver: a particular solution 37 satisfying the divergence constraint, the remaining div-free component of 38 the flux, and the pressure. When coupled_solve == false, the three 39 components will be solved one by one in the aforementioned order. 40 Otherwise, they will be solved at the same time. */ 41 bool coupled_solve = false; 42 bool verbose = false; 43 IterSolveParameters coarse_solve_param; 44 IterSolveParameters BBT_solve_param; 45 }; 46 47 /// Data for the divergence free solver 48 struct DFSData 49 { 50 std::vector<OperatorPtr> agg_hdivdof; // agglomerates to H(div) dofs table 51 std::vector<OperatorPtr> agg_l2dof; // agglomerates to L2 dofs table 52 std::vector<OperatorPtr> P_hdiv; // Interpolation matrix for H(div) space 53 std::vector<OperatorPtr> P_l2; // Interpolation matrix for L2 space 54 std::vector<OperatorPtr> P_hcurl; // Interpolation for kernel space of div 55 std::vector<OperatorPtr> Q_l2; // Q_l2[l] = (W_{l+1})^{-1} P_l2[l]^T W_l 56 Array<int> coarsest_ess_hdivdofs; // coarsest level essential H(div) dofs 57 std::vector<OperatorPtr> C; // discrete curl: ND -> RT, map to Null(B) 58 DFSParameters param; 59 }; 60 61 /// Finite element spaces concerning divergence free solver. 62 /// The main usage of this class is to collect data needed for the solver. 63 class DFSSpaces 64 { 65 RT_FECollection hdiv_fec_; 66 L2_FECollection l2_fec_; 67 std::unique_ptr<FiniteElementCollection> hcurl_fec_; 68 L2_FECollection l2_0_fec_; 69 70 std::unique_ptr<ParFiniteElementSpace> coarse_hdiv_fes_; 71 std::unique_ptr<ParFiniteElementSpace> coarse_l2_fes_; 72 std::unique_ptr<ParFiniteElementSpace> coarse_hcurl_fes_; 73 std::unique_ptr<ParFiniteElementSpace> l2_0_fes_; 74 75 std::unique_ptr<ParFiniteElementSpace> hdiv_fes_; 76 std::unique_ptr<ParFiniteElementSpace> l2_fes_; 77 std::unique_ptr<ParFiniteElementSpace> hcurl_fes_; 78 79 std::vector<SparseMatrix> el_l2dof_; 80 const Array<int>& ess_bdr_attr_; 81 Array<int> all_bdr_attr_; 82 83 int level_; 84 DFSData data_; 85 86 void MakeDofRelationTables(int level); 87 void DataFinalize(); 88 public: 89 DFSSpaces(int order, int num_refine, ParMesh *mesh, 90 const Array<int>& ess_attr, const DFSParameters& param); 91 92 /** This should be called each time when the mesh (where the FE spaces are 93 defined) is refined. The spaces will be updated, and the prolongation for 94 the spaces and other data needed for the div-free solver are stored. */ 95 void CollectDFSData(); 96 GetDFSData() const97 const DFSData& GetDFSData() const { return data_; } GetHdivFES() const98 ParFiniteElementSpace* GetHdivFES() const { return hdiv_fes_.get(); } GetL2FES() const99 ParFiniteElementSpace* GetL2FES() const { return l2_fes_.get(); } 100 }; 101 102 /// Abstract solver class for Darcy's flow 103 class DarcySolver : public Solver 104 { 105 protected: 106 Array<int> offsets_; 107 public: DarcySolver(int size0,int size1)108 DarcySolver(int size0, int size1) : Solver(size0 + size1), offsets_(3) 109 { offsets_[0] = 0; offsets_[1] = size0; offsets_[2] = height; } 110 virtual int GetNumIterations() const = 0; 111 }; 112 113 /// Solver for B * B^T 114 /// Compute the product B * B^T and solve it with CG preconditioned by BoomerAMG 115 class BBTSolver : public Solver 116 { 117 OperatorPtr BBT_; 118 OperatorPtr BBT_prec_; 119 CGSolver BBT_solver_; 120 public: 121 BBTSolver(const HypreParMatrix &B, IterSolveParameters param); Mult(const Vector & x,Vector & y) const122 virtual void Mult(const Vector &x, Vector &y) const { BBT_solver_.Mult(x, y); } SetOperator(const Operator & op)123 virtual void SetOperator(const Operator &op) { } 124 }; 125 126 /// Block diagonal solver for symmetric A, each block is inverted by direct solver 127 class SymDirectSubBlockSolver : public DirectSubBlockSolver 128 { 129 public: SymDirectSubBlockSolver(const SparseMatrix & A,const SparseMatrix & block_dof)130 SymDirectSubBlockSolver(const SparseMatrix& A, const SparseMatrix& block_dof) 131 : DirectSubBlockSolver(A, block_dof) { } MultTranspose(const Vector & x,Vector & y) const132 virtual void MultTranspose(const Vector &x, Vector &y) const { Mult(x, y); } 133 }; 134 135 /// non-overlapping additive Schwarz smoother for saddle point systems 136 /// [ M B^T ] 137 /// [ B 0 ] 138 class SaddleSchwarzSmoother : public Solver 139 { 140 const SparseMatrix& agg_hdivdof_; 141 const SparseMatrix& agg_l2dof_; 142 OperatorPtr coarse_l2_projector_; 143 144 Array<int> offsets_; 145 mutable Array<int> offsets_loc_; 146 mutable Array<int> hdivdofs_loc_; 147 mutable Array<int> l2dofs_loc_; 148 std::vector<OperatorPtr> solvers_loc_; 149 public: 150 /** SaddleSchwarzSmoother solves local saddle point problems defined on a 151 list of non-overlapping aggregates (of elements). 152 @param M the [1,1] block of the saddle point system 153 @param B the [2,1] block of the saddle point system 154 @param agg_hdivdof aggregate to H(div) dofs relation table (boolean matrix) 155 @param agg_l2dof aggregate to L2 dofs relation table (boolean matrix) 156 @param P_l2 prolongation matrix of the discrete L2 space 157 @param Q_l2 projection matrix of the discrete L2 space: 158 Q_l2 := (P_l2 W P_l2)^{-1} * P_l2 * W, 159 where W is the mass matrix of the discrete L2 space. */ 160 SaddleSchwarzSmoother(const HypreParMatrix& M, 161 const HypreParMatrix& B, 162 const SparseMatrix& agg_hdivdof, 163 const SparseMatrix& agg_l2dof, 164 const HypreParMatrix& P_l2, 165 const HypreParMatrix& Q_l2); 166 virtual void Mult(const Vector &x, Vector &y) const; MultTranspose(const Vector & x,Vector & y) const167 virtual void MultTranspose(const Vector &x, Vector &y) const { Mult(x, y); } SetOperator(const Operator & op)168 virtual void SetOperator(const Operator &op) { } 169 }; 170 171 /// Solver for local problems in SaddleSchwarzSmoother 172 class LocalSolver : public Solver 173 { 174 DenseMatrix local_system_; 175 DenseMatrixInverse local_solver_; 176 const int offset_; 177 public: 178 LocalSolver(const DenseMatrix &M, const DenseMatrix &B); 179 virtual void Mult(const Vector &x, Vector &y) const; SetOperator(const Operator & op)180 virtual void SetOperator(const Operator &op) { } 181 }; 182 183 /// Wrapper for the block-diagonal-preconditioned MINRES defined in ex5p.cpp 184 class BDPMinresSolver : public DarcySolver 185 { 186 BlockOperator op_; 187 BlockDiagonalPreconditioner prec_; 188 OperatorPtr BT_; 189 OperatorPtr S_; // S_ = B diag(M)^{-1} B^T 190 MINRESSolver solver_; 191 Array<int> ess_zero_dofs_; 192 public: 193 BDPMinresSolver(HypreParMatrix& M, HypreParMatrix& B, 194 IterSolveParameters param); 195 virtual void Mult(const Vector & x, Vector & y) const; SetOperator(const Operator & op)196 virtual void SetOperator(const Operator &op) { } SetEssZeroDofs(const Array<int> & dofs)197 void SetEssZeroDofs(const Array<int>& dofs) { dofs.Copy(ess_zero_dofs_); } GetNumIterations() const198 virtual int GetNumIterations() const { return solver_.GetNumIterations(); } 199 }; 200 201 /** Divergence free solver. 202 The basic idea of the solver is to exploit a multilevel decomposition of 203 Raviart-Thomas space to find a particular solution satisfying the divergence 204 constraint, and then solve the remaining (divergence-free) component in the 205 kernel space of the discrete divergence operator. 206 207 For more details see 208 1. Vassilevski, Multilevel Block Factorization Preconditioners (Appendix F.3), 209 Springer, 2008. 210 2. Voronin, Lee, Neumuller, Sepulveda, Vassilevski, Space-time discretizations 211 using constrained first-order system least squares (CFOSLS). 212 J. Comput. Phys. 373: 863-876, 2018. */ 213 class DivFreeSolver : public DarcySolver 214 { 215 const DFSData& data_; 216 DFSParameters param_; 217 OperatorPtr BT_; 218 BBTSolver BBT_solver_; 219 std::vector<Array<int>> ops_offsets_; 220 Array<BlockOperator*> ops_; 221 Array<BlockOperator*> blk_Ps_; 222 Array<Solver*> smoothers_; 223 OperatorPtr prec_; 224 OperatorPtr solver_; 225 226 void SolveParticular(const Vector& rhs, Vector& sol) const; 227 void SolveDivFree(const Vector& rhs, Vector& sol) const; 228 void SolvePotential(const Vector &rhs, Vector& sol) const; 229 public: 230 DivFreeSolver(const HypreParMatrix& M, const HypreParMatrix &B, 231 const DFSData& data); 232 ~DivFreeSolver(); 233 virtual void Mult(const Vector &x, Vector &y) const; SetOperator(const Operator & op)234 virtual void SetOperator(const Operator &op) { } 235 virtual int GetNumIterations() const; 236 }; 237 238 } // namespace blocksolvers 239 240 } // namespace mfem 241 242 #endif // MFEM_DIVFREE_SOLVER_HPP 243