MFEM  v4.5.1
Finite element discretization library
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Pages
auxiliary.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_AUXILIARY
13 #define MFEM_AUXILIARY
14 
15 #include "../config/config.hpp"
16 
17 #ifdef MFEM_USE_MPI
18 
19 #include "../general/tic_toc.hpp"
20 #include "solvers.hpp"
21 
22 namespace mfem
23 {
24 
25 // forward declarations
26 class Coefficient;
27 class MatrixCoefficient;
28 class ParMesh;
29 class ParBilinearForm;
30 class ParDiscreteLinearOperator;
31 
32 /** @brief Auxiliary space solvers for MatrixFreeAMS preconditioner
33 
34  Given an operator A and a transfer G, create a solver that approximates
35  (G^T A G)^{-1}. Used for two different auxiliary spaces in the AMS cycle.
36 
37  The produced solver is based on a low-order refined discretization for the
38  high-order H1 problem. */
40 {
41 public:
42  /** @brief Pi space constructor
43 
44  In the AMS framework this auxiliary space has two coefficients.
45 
46  @param mesh_lor Low-order refined auxiliary mesh
47  @param alpha_coeff Coefficient on curl-curl term (1 if null)
48  @param beta_coeff Coefficient on mass term (1 if null)
49  @param beta_mcoeff Matrix coefficient on mass term
50  @param ess_bdr Attributes for essential boundaries
51  @param curlcurl_oper High-order operator for the system
52  @param pi Identity interpolation operator
53  @param useAmgX_ Use AmgX instead of hypre for auxiliary solves
54  @param cg_iterations Number of CG iterations used to invert auxiliary
55  system, choosing 0 means to use a single V-cycle
56  */
58  ParMesh& mesh_lor, Coefficient* alpha_coeff,
59  Coefficient* beta_coeff, MatrixCoefficient* beta_mcoeff,
60  Array<int>& ess_bdr, Operator& curlcurl_oper, Operator& pi,
61 #ifdef MFEM_USE_AMGX
62  bool useAmgX_,
63 #endif
64  int cg_iterations = 0);
65 
66  /** @brief G space constructor
67 
68  This has one coefficient in the AMS framework.
69 
70  @param mesh_lor Low-order refined auxiliary mesh
71  @param beta_coeff Coefficient on mass term (1 if null)
72  @param beta_mcoeff Matrix coefficient on mass term
73  @param ess_bdr Attributes for essential boundaries
74  @param curlcurl_oper High-order operator for the system
75  @param g Gradient interpolation operator
76  @param useAmgX_ Use AmgX instead of hypre for auxiliary solves
77  @param cg_iterations Number of CG iterations used to invert auxiliary
78  system, choosing 0 means to use a single V-cycle
79  */
81  ParMesh& mesh_lor, Coefficient* beta_coeff,
82  MatrixCoefficient* beta_mcoeff, Array<int>& ess_bdr,
83  Operator& curlcurl_oper, Operator& g,
84 #ifdef MFEM_USE_AMGX
85  bool useAmgX_,
86 #endif
87  int cg_iterations = 1);
88 
90 
91  void Mult(const Vector& x, Vector& y) const;
92 
93  void SetOperator(const Operator& op) {}
94 
95 private:
96  /** @brief Helper routine for constructors.
97 
98  @param system_dimension is passed to HypreBoomerAMG::SetSystemsOptions
99  */
100  void SetupAMG(int system_dimension);
101  void SetupVCycle();
102 
103  /// inner_cg_iterations > 99 applies an exact solve here
104  void SetupCG(Operator& curlcurl_oper, Operator& conn,
105  int inner_cg_iterations);
106 
107  MPI_Comm comm;
108  Array<int> ess_tdof_list;
109  HypreParMatrix * lor_matrix;
110  Solver * lor_pc;
111  Operator* matfree;
112  CGSolver* cg;
113  Operator* aspacewrapper;
114 #ifdef MFEM_USE_AMGX
115  const bool useAmgX;
116 #endif
117  mutable int inner_aux_iterations;
118 };
119 
120 
121 /** @brief Perform AMS cycle with generic Operator objects.
122 
123  Most users should use MatrixFreeAMS, which wraps this. */
124 class GeneralAMS : public Solver
125 {
126 public:
127  /** @brief Constructor.
128 
129  Most of these arguments just need a Mult() operation, but pi and g also
130  require MultTranspose() */
131  GeneralAMS(const Operator& curlcurl_op_,
132  const Operator& pi_,
133  const Operator& gradient_,
134  const Operator& pispacesolver_,
135  const Operator& gspacesolver_,
136  const Operator& smoother_,
137  const Array<int>& ess_tdof_list_);
138  virtual ~GeneralAMS();
139 
140  /// in principle this should set A_ = op;
141  void SetOperator(const Operator &op) {}
142 
143  virtual void Mult(const Vector& x, Vector& y) const;
144 
145 private:
146  const Operator& curlcurl_op;
147  const Operator& pi;
148  const Operator& gradient;
149  const Operator& pispacesolver;
150  const Operator& gspacesolver;
151  const Operator& smoother;
152  const Array<int> ess_tdof_list;
153 
154  void FormResidual(const Vector& rhs, const Vector& x,
155  Vector& residual) const;
156 };
157 
158 
159 /** @brief An auxiliary Maxwell solver for a high-order curl-curl system without
160  high-order assembly.
161 
162  The auxiliary space solves are done using a low-order refined approach, but
163  all the interpolation operators, residuals, etc. are done in a matrix-free
164  manner.
165 
166  See Barker and Kolev, Matrix-free preconditioning for high-order H(curl)
167  discretizations (https://doi.org/10.1002/nla.2348) */
168 class MatrixFreeAMS : public Solver
169 {
170 public:
171  /** @brief Construct matrix-free AMS preconditioner
172 
173  @param aform BilinearForm for curl-curl problem, generally will
174  have a CurlCurlIntegrator and possibly a
175  VectorFEMassIntegrator.
176  @param oper Operator to precondition.
177  @param nd_fespace Underlying Nedelec finite element space.
178  @param alpha_coeff Coefficient on curl-curl term in Maxwell problem
179  (can be null, in which case constant 1 is assumed)
180  @param beta_coeff (scalar) coefficient on mass term in Maxwell problem
181  @param beta_mcoeff (matrix) coefficient on mass term
182  @param ess_bdr Boundary *attributes* that are marked essential. In
183  contrast to other MFEM cases, these are *attributes*
184  not dofs, because we need to apply these boundary
185  conditions to different bilinear forms.
186  @param useAmgX Use AmgX (instead of hypre) for LOR problems
187  @param inner_pi_its Number of CG iterations on auxiliary pi space,
188  may need more for difficult coefficients
189  @param inner_g_its Number of CG iterations on auxiliary g space,
190  may need more for difficult coefficients
191  @param nd_smoother Optional user-provided smoother for Nedelec space,
192  this object takes ownership and will delete.
193  */
194  MatrixFreeAMS(ParBilinearForm& aform, Operator& oper,
195  ParFiniteElementSpace& nd_fespace, Coefficient* alpha_coeff,
196  Coefficient* beta_coeff, MatrixCoefficient* beta_mcoeff,
197  Array<int>& ess_bdr,
198 #ifdef MFEM_USE_AMGX
199  bool useAmgX = false,
200 #endif
201  int inner_pi_its = 0, int inner_g_its = 1,
202  Solver* nd_smoother = NULL);
203 
204  ~MatrixFreeAMS();
205 
206  void SetOperator(const Operator &op) {}
207 
208  void Mult(const Vector& x, Vector& y) const { general_ams->Mult(x, y); }
209 
210 private:
211  GeneralAMS * general_ams;
212 
213  Solver * smoother;
214  ParDiscreteLinearOperator * pa_grad;
215  OperatorPtr Gradient;
216  ParDiscreteLinearOperator * pa_interp;
217  OperatorPtr Pi;
218 
219  Solver * Gspacesolver;
220  Solver * Pispacesolver;
221 
222  ParFiniteElementSpace * h1_fespace;
223  ParFiniteElementSpace * h1_fespace_d;
224 };
225 
226 } // namespace mfem
227 
228 #endif // MFEM_USE_MPI
229 
230 #endif
Conjugate gradient method.
Definition: solvers.hpp:465
Pointer to an Operator of a specified type.
Definition: handle.hpp:33
void SetOperator(const Operator &op)
in principle this should set A_ = op;
Definition: auxiliary.hpp:141
Abstract parallel finite element space.
Definition: pfespace.hpp:28
MatrixFreeAuxiliarySpace(ParMesh &mesh_lor, Coefficient *alpha_coeff, Coefficient *beta_coeff, MatrixCoefficient *beta_mcoeff, Array< int > &ess_bdr, Operator &curlcurl_oper, Operator &pi,#ifdef MFEM_USE_AMGX bool useAmgX_,#endif int cg_iterations=0)
Pi space constructor.
Definition: auxiliary.cpp:120
An auxiliary Maxwell solver for a high-order curl-curl system without high-order assembly.
Definition: auxiliary.hpp:168
Perform AMS cycle with generic Operator objects.
Definition: auxiliary.hpp:124
void Mult(const Vector &x, Vector &y) const
Operator application: y=A(x).
Definition: auxiliary.hpp:208
void Mult(const Vector &x, Vector &y) const
Operator application: y=A(x).
Definition: auxiliary.cpp:399
void SetOperator(const Operator &op)
Set/update the solver for the given operator.
Definition: auxiliary.hpp:206
void SetOperator(const Operator &op)
Set/update the solver for the given operator.
Definition: auxiliary.hpp:93
Base class Coefficients that optionally depend on space and time. These are used by the BilinearFormI...
Definition: coefficient.hpp:41
Base class for Matrix Coefficients that optionally depend on time and space.
Auxiliary space solvers for MatrixFreeAMS preconditioner.
Definition: auxiliary.hpp:39
GeneralAMS(const Operator &curlcurl_op_, const Operator &pi_, const Operator &gradient_, const Operator &pispacesolver_, const Operator &gspacesolver_, const Operator &smoother_, const Array< int > &ess_tdof_list_)
Constructor.
Definition: auxiliary.cpp:23
Class for parallel bilinear form.
Vector data type.
Definition: vector.hpp:60
MatrixFreeAMS(ParBilinearForm &aform, Operator &oper, ParFiniteElementSpace &nd_fespace, Coefficient *alpha_coeff, Coefficient *beta_coeff, MatrixCoefficient *beta_mcoeff, Array< int > &ess_bdr,#ifdef MFEM_USE_AMGX bool useAmgX=false,#endif int inner_pi_its=0, int inner_g_its=1, Solver *nd_smoother=NULL)
Construct matrix-free AMS preconditioner.
Definition: auxiliary.cpp:426
Base class for solvers.
Definition: operator.hpp:655
Abstract operator.
Definition: operator.hpp:24
Wrapper for hypre&#39;s ParCSR matrix class.
Definition: hypre.hpp:343
Class for parallel meshes.
Definition: pmesh.hpp:32
virtual void Mult(const Vector &x, Vector &y) const
Operator application: y=A(x).
Definition: auxiliary.cpp:70
virtual ~GeneralAMS()
Definition: auxiliary.cpp:42