MFEM  v4.5.1
Finite element discretization library
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Pages
dgmassinv.hpp
Go to the documentation of this file.
1 // Copyright (c) 2010-2022, 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_DGMASSINV_HPP
13 #define MFEM_DGMASSINV_HPP
14 
15 #include "../linalg/operator.hpp"
16 #include "fespace.hpp"
17 
18 namespace mfem
19 {
20 
21 /// @brief Solver for the discontinuous Galerkin mass matrix.
22 ///
23 /// This class performs a @a local (diagonally preconditioned) conjugate
24 /// gradient iteration for each element. Optionally, a change of basis is
25 /// performed to iterate on a better-conditioned system. This class fully
26 /// supports execution on device (GPU).
27 class DGMassInverse : public Solver
28 {
29 protected:
30  DG_FECollection fec; ///< FE collection in requested basis.
31  FiniteElementSpace fes; ///< FE space in requested basis.
32  const DofToQuad *d2q; ///< Change of basis. Not owned.
33  Array<double> B_; ///< Inverse of change of basis.
34  Array<double> Bt_; ///< Inverse of change of basis, transposed.
35  class BilinearForm *M; ///< Mass bilinear form, owned.
36  class MassIntegrator *m; ///< Mass integrator, owned by the form @ref M.
37  Vector diag_inv; ///< Jacobi preconditioner.
38  double rel_tol = 1e-12; ///< Relative CG tolerance.
39  double abs_tol = 1e-12; ///< Absolute CG tolerance.
40  int max_iter = 100; ///< Maximum number of CG iterations;
41 
42  /// @name Intermediate vectors needed for CG three-term recurrence.
43  ///@{
44  mutable Vector r_, d_, z_, b2_;
45  ///@}
46 
47  /// @brief Protected constructor, used internally.
48  ///
49  /// Custom coefficient and integration rule are used if @a coeff and @a ir
50  /// are non-NULL.
52  const IntegrationRule *ir, int btype);
53 public:
54  /// @brief Construct the DG inverse mass operator for @a fes_.
55  ///
56  /// The basis type @a btype determines which basis should be used internally
57  /// in the solver. This <b>does not</b> have to be the same basis as @a fes_.
58  /// The best choice is typically BasisType::GaussLegendre because it is
59  /// well-preconditioned by its diagonal.
60  ///
61  /// The solution and right-hand side used for the solver are not affected by
62  /// this basis (they correspond to the basis of @a fes_). @a btype is only
63  /// used internally, and only has an effect on the convergence rate.
65  /// @brief Construct the DG inverse mass operator for @a fes_ with
66  /// Coefficient @a coeff.
67  ///
68  /// @sa DGMassInverse(FiniteElementSpace&, int) for information about @a
69  /// btype.
71  int btype=BasisType::GaussLegendre);
72  /// @brief Construct the DG inverse mass operator for @a fes_ with
73  /// Coefficient @a coeff and IntegrationRule @a ir.
74  ///
75  /// @sa DGMassInverse(FiniteElementSpace&, int) for information about @a
76  /// btype.
78  const IntegrationRule &ir, int btype=BasisType::GaussLegendre);
79  /// @brief Construct the DG inverse mass operator for @a fes_ with
80  /// IntegrationRule @a ir.
81  ///
82  /// @sa DGMassInverse(FiniteElementSpace&, int) for information about @a
83  /// btype.
85  int btype=BasisType::GaussLegendre);
86  /// @brief Solve the system M b = u.
87  ///
88  /// If @ref iterative_mode is @a true, @a u is used as an initial guess.
89  void Mult(const Vector &b, Vector &u) const;
90  /// Not implemented. Aborts.
91  void SetOperator(const Operator &op);
92  /// Set the relative tolerance.
93  void SetRelTol(const double rel_tol_);
94  /// Set the absolute tolerance.
95  void SetAbsTol(const double abs_tol_);
96  /// Set the maximum number of iterations.
97  void SetMaxIter(const double max_iter_);
98  /// Recompute operator and preconditioner (when coefficient or mesh changes).
99  void Update();
100 
101  ~DGMassInverse();
102 
103  /// @brief Solve the system M b = u. <b>Not part of the public interface.</b>
104  /// @note This member function must be public because it contains an
105  /// MFEM_FORALL kernel (nvcc limitation)
106  template<int DIM, int D1D = 0, int Q1D = 0>
107  void DGMassCGIteration(const Vector &b_, Vector &u_) const;
108 };
109 
110 } // namespace mfem
111 
112 #endif
Class for an integration rule - an Array of IntegrationPoint.
Definition: intrules.hpp:90
void SetOperator(const Operator &op)
Not implemented. Aborts.
Definition: dgmassinv.cpp:95
DGMassInverse(FiniteElementSpace &fes_, Coefficient *coeff, const IntegrationRule *ir, int btype)
Protected constructor, used internally.
Definition: dgmassinv.cpp:20
void SetAbsTol(const double abs_tol_)
Set the absolute tolerance.
Definition: dgmassinv.cpp:102
Array< double > Bt_
Inverse of change of basis, transposed.
Definition: dgmassinv.hpp:34
void SetMaxIter(const double max_iter_)
Set the maximum number of iterations.
Definition: dgmassinv.cpp:104
double abs_tol
Absolute CG tolerance.
Definition: dgmassinv.hpp:39
double b
Definition: lissajous.cpp:42
double rel_tol
Relative CG tolerance.
Definition: dgmassinv.hpp:38
Array< double > B_
Inverse of change of basis.
Definition: dgmassinv.hpp:33
Class FiniteElementSpace - responsible for providing FEM view of the mesh, mainly managing the set of...
Definition: fespace.hpp:96
Base class Coefficients that optionally depend on space and time. These are used by the BilinearFormI...
Definition: coefficient.hpp:41
DG_FECollection fec
FE collection in requested basis.
Definition: dgmassinv.hpp:30
Structure representing the matrices/tensors needed to evaluate (in reference space) the values...
Definition: fe_base.hpp:136
void Mult(const Vector &b, Vector &u) const
Solve the system M b = u.
Definition: dgmassinv.cpp:265
class BilinearForm * M
Mass bilinear form, owned.
Definition: dgmassinv.hpp:35
int max_iter
Maximum number of CG iterations;.
Definition: dgmassinv.hpp:40
void DGMassCGIteration(const Vector &b_, Vector &u_) const
Solve the system M b = u. Not part of the public interface.
Definition: dgmassinv.cpp:119
A &quot;square matrix&quot; operator for the associated FE space and BLFIntegrators The sum of all the BLFInteg...
class MassIntegrator * m
Mass integrator, owned by the form M.
Definition: dgmassinv.hpp:36
const DofToQuad * d2q
Change of basis. Not owned.
Definition: dgmassinv.hpp:32
Vector diag_inv
Jacobi preconditioner.
Definition: dgmassinv.hpp:37
Vector data type.
Definition: vector.hpp:60
double u(const Vector &xvec)
Definition: lor_mms.hpp:24
Base class for solvers.
Definition: operator.hpp:655
Abstract operator.
Definition: operator.hpp:24
FiniteElementSpace fes
FE space in requested basis.
Definition: dgmassinv.hpp:31
void SetRelTol(const double rel_tol_)
Set the relative tolerance.
Definition: dgmassinv.cpp:100
Solver for the discontinuous Galerkin mass matrix.
Definition: dgmassinv.hpp:27
void Update()
Recompute operator and preconditioner (when coefficient or mesh changes).
Definition: dgmassinv.cpp:106
Arbitrary order &quot;L2-conforming&quot; discontinuous finite elements.
Definition: fe_coll.hpp:288