MFEM v4.7.0
Finite element discretization library
Loading...
Searching...
No Matches
hdiv_linear_solver.hpp
Go to the documentation of this file.
1// Copyright (c) 2010-2024, 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 HDIV_LINEAR_SOLVER_HPP
13#define HDIV_LINEAR_SOLVER_HPP
14
15#include "mfem.hpp"
16#include "change_basis.hpp"
17#include <memory>
18
19namespace mfem
20{
21
22/// @brief Solve the H(div) saddle-point system using MINRES with matrix-free
23/// block-diagonal preconditioning.
24///
25/// See HdivSaddlePointSolver::HdivSaddlePointSolver for the problem
26/// description.
28{
29public:
30 /// Which type of saddle-point problem is being solved?
31 enum Mode
32 {
33 GRAD_DIV, ///< Grad-div problem.
34 DARCY ///< Darcy/mixed Poisson problem.
35 };
36private:
37 MINRESSolver minres;
38
39 static constexpr int b1 = BasisType::GaussLobatto;
40 static constexpr int b2 = BasisType::IntegratedGLL;
41 static constexpr int mt = FiniteElement::INTEGRAL;
42
43 const int order;
44
45 // L2 and RT spaces, using the interpolation-histopolation bases
46 L2_FECollection fec_l2;
48
49 RT_FECollection fec_rt;
51
52 const Array<int> &ess_rt_dofs; ///< Essential BCs (in the RT space only).
53
54 // Change of basis operators
55 ChangeOfBasis_L2 basis_l2;
56 ChangeOfBasis_RT basis_rt;
57
58 /// Whether conversion from map type VALUE to INTEGRAL is required.
59 const bool convert_map_type;
60
61 ParBilinearForm mass_l2, mass_rt;
62
63 // Components needed for the block operator
64 OperatorHandle L, R, R_e; ///< Mass matrices.
65 std::unique_ptr<HypreParMatrix> D, Dt, D_e; ///< Divergence matrices.
66 std::shared_ptr<DGMassInverse> L_inv; ///< Inverse of the DG mass matrix.
67 std::shared_ptr<Operator> A_11; ///< (1,1)-block of the matrix
68
69 /// Diagonals of the mass matrices
70 Vector L_diag, R_diag, L_diag_unweighted;
71
72 // Components needed for the preconditioner
73
74 /// Jacobi preconditioner for the RT mass matrix.
75 std::unique_ptr<OperatorJacobiSmoother> R_inv;
76 std::unique_ptr<HypreParMatrix> S; ///< Approximate Schur complement.
77 HypreBoomerAMG S_inv; ///< AMG preconditioner for #S.
78
79 Array<int> offsets, empty;
80 /// The 2x2 block operator.
81 std::unique_ptr<BlockOperator> A_block;
82 /// The block-diagonal preconditioner.
83 std::unique_ptr<BlockDiagonalPreconditioner> D_prec;
84
85 Coefficient &L_coeff, &R_coeff;
86
87 const Mode mode;
88 bool zero_l2_block = false;
90 QuadratureFunction W_coeff_qf, W_mix_coeff_qf;
91 QuadratureFunctionCoefficient W_coeff, W_mix_coeff;
92
94
95 // Work vectors
96 mutable Vector b_prime, x_prime, x_bc, w, z;
97public:
98 /// @brief Creates a solver for the H(div) saddle-point system.
99 ///
100 /// The associated matrix is given by
101 ///
102 /// [ L B ]
103 /// [ B^T -R ]
104 ///
105 /// where L is the L2 mass matrix, R is the RT mass matrix, and B is the
106 /// divergence form (VectorFEDivergenceIntegrator).
107 ///
108 /// Essential boundary conditions in the RT space are given by @a
109 /// ess_rt_dofs_. (Rows and columns are eliminated from R and columns are
110 /// eliminated from B).
111 ///
112 /// The L block has coefficient @a L_coeff_ and the R block has coefficient
113 /// @a R_coeff_.
114 ///
115 /// The parameter @a mode_ determines whether the block system corresponds to
116 /// a grad-div problem or a Darcy problem. Specifically, if @a mode_ is
117 /// Mode::GRAD_DIV, then the B and B^T blocks are also scaled by @a L_coeff_,
118 /// and if @a mode_ is Mode::DARCY, then the B and B^T blocks are unweighted.
119 ///
120 /// Mode::GRAD_DIV corresponds to the grad-div problem
121 ///
122 /// alpha u - grad ( beta div ( u )) = f,
123 ///
124 /// where alpha is @a R_coeff_ and beta is @a L_coeff_.
125 ///
126 /// Mode::DARCY corresponds to the Darcy-type problem
127 ///
128 /// alpha p - div ( beta grad ( p )) = f,
129 ///
130 /// where alpha is @a L_coeff and beta is @a R_coeff_. In this case, the
131 /// coefficient alpha is allowed to be zero.
133 ParFiniteElementSpace &fes_rt_,
134 ParFiniteElementSpace &fes_l2_,
135 Coefficient &L_coeff_,
136 Coefficient &R_coeff_,
137 const Array<int> &ess_rt_dofs_,
138 Mode mode_);
139
140 /// @brief Creates a linear solver for the case when the L2 diagonal block is
141 /// zero (for Darcy problems).
142 ///
143 /// Equivalent to passing ConstantCoefficient(0.0) as @a L_coeff_ and
144 /// Mode::DARCY as @a mode_ to the primary constructor (see above).
146 ParFiniteElementSpace &fes_rt_,
147 ParFiniteElementSpace &fes_l2_,
148 Coefficient &R_coeff_,
149 const Array<int> &ess_rt_dofs_);
150
151 /// @brief Build the linear operator and solver. Must be called when the
152 /// coefficients change.
153 void Setup();
154 /// Sets the Dirichlet boundary conditions at the RT essential DOFs.
155 void SetBC(const Vector &x_rt) { x_bc = x_rt; }
156 /// @brief Solve the linear system for L2 (scalar) and RT (flux) unknowns.
157 ///
158 /// If the problem has essential boundary conditions (i.e. if @a ess_rt_dofs
159 /// is not empty), then SetBC() must be called before Mult().
160 void Mult(const Vector &b, Vector &x) const override;
161 /// No-op.
162 void SetOperator(const Operator &op) override { }
163 /// Get the number of MINRES iterations.
164 int GetNumIterations() const { return minres.GetNumIterations(); }
165 /// Eliminates the BCs (called internally, not public interface).
166 void EliminateBC(Vector &) const;
167 /// Return the offsets of the block system.
168 const Array<int> &GetOffsets() const { return offsets; }
169 /// Returns the internal MINRES solver.
170 MINRESSolver &GetMINRES() { return minres; }
171};
172
173} // namespace mfem
174
175#endif
@ GaussLobatto
Closed type.
Definition fe_base.hpp:32
@ IntegratedGLL
Integrated GLL indicator functions.
Definition fe_base.hpp:39
Change of basis operator between L2 spaces.
Base class Coefficients that optionally depend on space and time. These are used by the BilinearFormI...
A coefficient that is constant across space and time.
Solve the H(div) saddle-point system using MINRES with matrix-free block-diagonal preconditioning.
HdivSaddlePointSolver(ParMesh &mesh_, ParFiniteElementSpace &fes_rt_, ParFiniteElementSpace &fes_l2_, Coefficient &L_coeff_, Coefficient &R_coeff_, const Array< int > &ess_rt_dofs_, Mode mode_)
Creates a solver for the H(div) saddle-point system.
void Setup()
Build the linear operator and solver. Must be called when the coefficients change.
const Array< int > & GetOffsets() const
Return the offsets of the block system.
MINRESSolver & GetMINRES()
Returns the internal MINRES solver.
void EliminateBC(Vector &) const
Eliminates the BCs (called internally, not public interface).
Mode
Which type of saddle-point problem is being solved?
@ DARCY
Darcy/mixed Poisson problem.
void Mult(const Vector &b, Vector &x) const override
Solve the linear system for L2 (scalar) and RT (flux) unknowns.
void SetBC(const Vector &x_rt)
Sets the Dirichlet boundary conditions at the RT essential DOFs.
void SetOperator(const Operator &op) override
No-op.
int GetNumIterations() const
Get the number of MINRES iterations.
The BoomerAMG solver in hypre.
Definition hypre.hpp:1691
int GetNumIterations() const
Returns the number of iterations taken during the last call to Mult()
Definition solvers.hpp:260
Arbitrary order "L2-conforming" discontinuous finite elements.
Definition fe_coll.hpp:330
MINRES method.
Definition solvers.hpp:628
Pointer to an Operator of a specified type.
Definition handle.hpp:34
Abstract operator.
Definition operator.hpp:25
Class for parallel bilinear form.
Abstract parallel finite element space.
Definition pfespace.hpp:29
Class for parallel meshes.
Definition pmesh.hpp:34
Quadrature function coefficient which requires that the quadrature rules used for this coefficient be...
Represents values or vectors of values at quadrature points on a mesh.
Definition qfunction.hpp:24
Class representing the storage layout of a QuadratureFunction.
Definition qspace.hpp:120
Arbitrary order H(div)-conforming Raviart-Thomas finite elements.
Definition fe_coll.hpp:386
Base class for solvers.
Definition operator.hpp:683
Vector data type.
Definition vector.hpp:80
real_t b
Definition lissajous.cpp:42