MFEM  v4.3.0
Finite element discretization library
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Pages
div_free_solver.hpp
Go to the documentation of this file.
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
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
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;
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)
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 
97  const DFSData& GetDFSData() const { return data_; }
98  ParFiniteElementSpace* GetHdivFES() const { return hdiv_fes_.get(); }
99  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:
107 public:
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:
122  virtual void Mult(const Vector &x, Vector &y) const { BBT_solver_.Mult(x, y); }
123  virtual void SetOperator(const Operator &op) { }
124 };
125 
126 /// Block diagonal solver for symmetric A, each block is inverted by direct solver
128 {
129 public:
131  : DirectSubBlockSolver(A, block_dof) { }
132  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 ]
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. */
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;
167  virtual void MultTranspose(const Vector &x, Vector &y) const { Mult(x, y); }
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;
180  virtual void SetOperator(const Operator &op) { }
181 };
182 
183 /// Wrapper for the block-diagonal-preconditioned MINRES defined in ex5p.cpp
185 {
186  BlockOperator op_;
188  OperatorPtr BT_;
189  OperatorPtr S_; // S_ = B diag(M)^{-1} B^T
190  MINRESSolver solver_;
191  Array<int> ess_zero_dofs_;
192 public:
194  IterSolveParameters param);
195  virtual void Mult(const Vector & x, Vector & y) const;
196  virtual void SetOperator(const Operator &op) { }
197  void SetEssZeroDofs(const Array<int>& dofs) { dofs.Copy(ess_zero_dofs_); }
198  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. */
214 {
215  const DFSData& data_;
216  DFSParameters param_;
217  OperatorPtr BT_;
218  BBTSolver BBT_solver_;
219  std::vector<Array<int>> ops_offsets_;
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;
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
Conjugate gradient method.
Definition: solvers.hpp:316
int GetNumIterations() const
Definition: solvers.hpp:103
virtual void Mult(const Vector &b, Vector &x) const
Operator application: y=A(x).
Definition: solvers.cpp:616
DivFreeSolver(const HypreParMatrix &M, const HypreParMatrix &B, const DFSData &data)
virtual void Mult(const Vector &x, Vector &y) const
Operator application: y=A(x).
Definition: solvers.cpp:3173
std::vector< OperatorPtr > Q_l2
Pointer to an Operator of a specified type.
Definition: handle.hpp:33
Abstract solver class for Darcy&#39;s flow.
void Copy(Array &copy) const
Create a copy of the internal array to the provided copy.
Definition: array.hpp:851
Data type dense matrix using column-major storage.
Definition: densemat.hpp:23
virtual void Mult(const Vector &x, Vector &y) const
Operator application: y=A(x).
Abstract parallel finite element space.
Definition: pfespace.hpp:28
SymDirectSubBlockSolver(const SparseMatrix &A, const SparseMatrix &block_dof)
SaddleSchwarzSmoother(const HypreParMatrix &M, const HypreParMatrix &B, const SparseMatrix &agg_hdivdof, const SparseMatrix &agg_l2dof, const HypreParMatrix &P_l2, const HypreParMatrix &Q_l2)
MINRES method.
Definition: solvers.hpp:426
std::vector< OperatorPtr > P_hcurl
virtual void MultTranspose(const Vector &x, Vector &y) const
Action of the transpose operator: y=A^t(x). The default behavior in class Operator is to generate an ...
IterSolveParameters coarse_solve_param
virtual void Mult(const Vector &x, Vector &y) const
Operator application: y=A(x).
LocalSolver(const DenseMatrix &M, const DenseMatrix &B)
virtual void Mult(const Vector &x, Vector &y) const
Operator application: y=A(x).
Data for the divergence free solver.
const DFSData & GetDFSData() const
A class to handle Block diagonal preconditioners in a matrix-free implementation. ...
Data type sparse matrix.
Definition: sparsemat.hpp:41
Wrapper for the block-diagonal-preconditioned MINRES defined in ex5p.cpp.
ParFiniteElementSpace * GetHdivFES() const
virtual void SetOperator(const Operator &op)
Set/update the solver for the given operator.
std::vector< OperatorPtr > P_hdiv
virtual void SetOperator(const Operator &op)
Set/update the solver for the given operator.
DarcySolver(int size0, int size1)
virtual void SetOperator(const Operator &op)
Set/update the solver for the given operator.
ParFiniteElementSpace * GetL2FES() const
Arbitrary order H(div)-conforming Raviart-Thomas finite elements.
Definition: fe_coll.hpp:341
virtual void MultTranspose(const Vector &x, Vector &y) const
Action of the transpose operator: y=A^t(x). The default behavior in class Operator is to generate an ...
BDPMinresSolver(HypreParMatrix &M, HypreParMatrix &B, IterSolveParameters param)
std::vector< OperatorPtr > agg_hdivdof
std::vector< OperatorPtr > agg_l2dof
Parameters for the divergence free solver.
void SetEssZeroDofs(const Array< int > &dofs)
int height
Dimension of the output / number of rows in the matrix.
Definition: operator.hpp:27
std::vector< OperatorPtr > C
virtual int GetNumIterations() const =0
Block diagonal solver for symmetric A, each block is inverted by direct solver.
virtual void SetOperator(const Operator &op)
Set/update the solver for the given operator.
virtual void Mult(const Vector &x, Vector &y) const
Operator application: y=A(x).
virtual void Mult(const Vector &x, Vector &y) const
Operator application: y=A(x).
Vector data type.
Definition: vector.hpp:60
Parameters for iterative solver.
virtual void SetOperator(const Operator &op)
Set/update the solver for the given operator.
Base class for solvers.
Definition: operator.hpp:648
DFSSpaces(int order, int num_refine, ParMesh *mesh, const Array< int > &ess_attr, const DFSParameters &param)
Solver for local problems in SaddleSchwarzSmoother.
Abstract operator.
Definition: operator.hpp:24
Wrapper for hypre&#39;s ParCSR matrix class.
Definition: hypre.hpp:277
Block diagonal solver for A, each block is inverted by direct solver.
Definition: solvers.hpp:911
A class to handle Block systems in a matrix-free implementation.
BBTSolver(const HypreParMatrix &B, IterSolveParameters param)
Class for parallel meshes.
Definition: pmesh.hpp:32
std::vector< OperatorPtr > P_l2
Arbitrary order &quot;L2-conforming&quot; discontinuous finite elements.
Definition: fe_coll.hpp:285