MFEM  v4.4.0
Finite element discretization library
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Pages
auxiliary.cpp
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 #include "../config/config.hpp"
13 
14 #ifdef MFEM_USE_MPI
15 
16 #include "linalg.hpp"
17 #include "../fem/pfespace.hpp"
18 #include "../fem/pbilinearform.hpp"
19 
20 namespace mfem
21 {
22 
23 GeneralAMS::GeneralAMS(const Operator& curlcurl_op_,
24  const Operator& pi_,
25  const Operator& gradient_,
26  const Operator& pispacesolver_,
27  const Operator& gspacesolver_,
28  const Operator& smoother_,
29  const Array<int>& ess_tdof_list_)
30  :
31  Solver(curlcurl_op_.Height()),
32  curlcurl_op(curlcurl_op_),
33  pi(pi_),
34  gradient(gradient_),
35  pispacesolver(pispacesolver_),
36  gspacesolver(gspacesolver_),
37  smoother(smoother_),
38  ess_tdof_list(ess_tdof_list_)
39 {
40 }
41 
43 {
44 }
45 
46 void GeneralAMS::FormResidual(const Vector& rhs, const Vector& x,
47  Vector& residual) const
48 {
49  curlcurl_op.Mult(x, residual);
50  residual *= -1.0;
51  residual += rhs;
52 }
53 
54 /*
55  This implementation follows that in hypre, see hypre_ParCSRSubspacePrec()
56  in hypre/src/parcsr_ls/ams.c and also hypre_AMSSolve() in the same file.
57 
58  hypre's default cyle (cycle 1) is "01210", ie, smooth, correct in space
59  1, correct in space 2, correct in space 1, smooth. Their space 1 is G and
60  space 2 is Pi by default.
61 
62  The MFEM interface in mfem::HypreAMS though picks cycle 13, or 034515430,
63  which separates the Pi-space solve into three separate (scalar) AMG solves
64  instead of a single vector solve.
65 
66  We choose below the hypre default, but we have experimented with some other
67  cycles; the short version is that they often work but the differences are
68  generally not large.
69 */
70 void GeneralAMS::Mult(const Vector& x, Vector& y) const
71 {
72  MFEM_ASSERT(x.Size() == y.Size(), "Sizes don't match!");
73  MFEM_ASSERT(curlcurl_op.Height() == x.Size(), "Sizes don't match!");
74 
75  Vector residual(x.Size());
76  residual = 0.0;
77  y = 0.0;
78 
79  // smooth
80  smoother.Mult(x, y);
81 
82  // g-space correction
83  FormResidual(x, y, residual);
84  Vector gspacetemp(gradient.Width());
85  gradient.MultTranspose(residual, gspacetemp);
86  Vector gspacecorrection(gradient.Width());
87  gspacecorrection = 0.0;
88  gspacesolver.Mult(gspacetemp, gspacecorrection);
89  gradient.Mult(gspacecorrection, residual);
90  y += residual;
91 
92  // pi-space correction
93  FormResidual(x, y, residual);
94  Vector pispacetemp(pi.Width());
95 
96  pi.MultTranspose(residual, pispacetemp);
97 
98  Vector pispacecorrection(pi.Width());
99  pispacecorrection = 0.0;
100  pispacesolver.Mult(pispacetemp, pispacecorrection);
101  pi.Mult(pispacecorrection, residual);
102  y += residual;
103 
104  // g-space correction
105  FormResidual(x, y, residual);
106  gradient.MultTranspose(residual, gspacetemp);
107  gspacecorrection = 0.0;
108  gspacesolver.Mult(gspacetemp, gspacecorrection);
109  gradient.Mult(gspacecorrection, residual);
110  y += residual;
111 
112  // smooth
113  FormResidual(x, y, residual);
114  Vector temp(x.Size());
115  smoother.Mult(residual, temp);
116  y += temp;
117 }
118 
119 // Pi-space constructor
121  ParMesh& mesh_lor, Coefficient* alpha_coeff,
122  Coefficient* beta_coeff, MatrixCoefficient* beta_mcoeff, Array<int>& ess_bdr,
123  Operator& curlcurl_oper, Operator& pi,
124 #ifdef MFEM_USE_AMGX
125  bool useAmgX_,
126 #endif
127  int cg_iterations) :
128  Solver(pi.Width()),
129  comm(mesh_lor.GetComm()),
130  matfree(NULL),
131  cg(NULL),
132 #ifdef MFEM_USE_AMGX
133  useAmgX(useAmgX_),
134 #endif
135  inner_aux_iterations(0)
136 {
137  H1_FECollection * fec_lor = new H1_FECollection(1, mesh_lor.Dimension());
138  ParFiniteElementSpace fespace_lor_d(&mesh_lor, fec_lor, mesh_lor.Dimension(),
140 
141  // build LOR AMG v-cycle
142  if (ess_bdr.Size())
143  {
144  fespace_lor_d.GetEssentialTrueDofs(ess_bdr, ess_tdof_list);
145  }
146  ParBilinearForm a_lor(&fespace_lor_d);
147 
148  // this choice of policy is important for the G-space solver, but
149  // also can make some difference here
151  a_lor.SetDiagonalPolicy(policy);
152  if (alpha_coeff == NULL)
153  {
155  }
156  else
157  {
158  a_lor.AddDomainIntegrator(new VectorDiffusionIntegrator(*alpha_coeff));
159  }
160 
161  if (beta_mcoeff != NULL)
162  {
163  MFEM_VERIFY(beta_coeff == NULL, "Only one beta coefficient should be defined.");
164  a_lor.AddDomainIntegrator(new VectorMassIntegrator(*beta_mcoeff));
165  }
166  else if (beta_coeff != NULL)
167  {
168  a_lor.AddDomainIntegrator(new VectorMassIntegrator(*beta_coeff));
169  }
170  else
171  {
173  }
174  a_lor.UsePrecomputedSparsity();
175  a_lor.Assemble();
176  a_lor.EliminateEssentialBC(ess_bdr, policy);
177  a_lor.Finalize();
178  lor_matrix = a_lor.ParallelAssemble();
179  lor_matrix->CopyRowStarts();
180  lor_matrix->CopyColStarts();
181 
182  SetupAMG(fespace_lor_d.GetMesh()->Dimension());
183 
184  if (cg_iterations > 0)
185  {
186  SetupCG(curlcurl_oper, pi, cg_iterations);
187  }
188  else
189  {
190  SetupVCycle();
191  }
192  delete fec_lor;
193 }
194 
195 /* G-space constructor
196 
197  The auxiliary space solves in general, and this one in particular,
198  seem to be quite sensitive to handling of boundary conditions. Note
199  some careful choices for Matrix::DiagonalPolicy and the ZeroWrapAMG
200  object, as well as the use of a single CG iteration (instead of just
201  an AMG V-cycle). Just a V-cycle may be more efficient in some cases,
202  but we recommend the CG wrapper for robustness here. */
204  ParMesh& mesh_lor, Coefficient* beta_coeff,
205  MatrixCoefficient* beta_mcoeff, Array<int>& ess_bdr, Operator& curlcurl_oper,
206  Operator& g,
207 #ifdef MFEM_USE_AMGX
208  bool useAmgX_,
209 #endif
210  int cg_iterations)
211  :
212  Solver(curlcurl_oper.Height()),
213  comm(mesh_lor.GetComm()),
214  matfree(NULL),
215  cg(NULL),
216 #ifdef MFEM_USE_AMGX
217  useAmgX(useAmgX_),
218 #endif
219  inner_aux_iterations(0)
220 {
221  H1_FECollection * fec_lor = new H1_FECollection(1, mesh_lor.Dimension());
222  ParFiniteElementSpace fespace_lor(&mesh_lor, fec_lor);
223 
224  // build LOR AMG v-cycle
225  ParBilinearForm a_lor(&fespace_lor);
226 
227  // we need something like DIAG_ZERO in the solver, but explicitly doing
228  // that makes BoomerAMG setup complain, so instead we constrain the boundary
229  // in the CG solver
231 
232  a_lor.SetDiagonalPolicy(policy);
233 
234  if (beta_mcoeff != NULL)
235  {
236  MFEM_VERIFY(beta_coeff == NULL, "Only one beta coefficient should be defined.");
237  a_lor.AddDomainIntegrator(new DiffusionIntegrator(*beta_mcoeff));
238  }
239  else if (beta_coeff != NULL)
240  {
241  a_lor.AddDomainIntegrator(new DiffusionIntegrator(*beta_coeff));
242  }
243  else
244  {
246  }
247  a_lor.UsePrecomputedSparsity();
248  a_lor.Assemble();
249  if (ess_bdr.Size())
250  {
251  fespace_lor.GetEssentialTrueDofs(ess_bdr, ess_tdof_list);
252  }
253 
254  // you have to use (serial) BilinearForm eliminate routines to get
255  // diag policy DIAG_ZERO all the ParallelEliminateTDofs etc. routines
256  // implicitly have a Matrix::DIAG_KEEP policy
257  a_lor.EliminateEssentialBC(ess_bdr, policy);
258  a_lor.Finalize();
259  lor_matrix = a_lor.ParallelAssemble();
260 
261  lor_matrix->CopyRowStarts();
262  lor_matrix->CopyColStarts();
263 
264  SetupAMG(0);
265 
266  if (cg_iterations > 0)
267  {
268  SetupCG(curlcurl_oper, g, cg_iterations);
269  }
270  else
271  {
272  SetupVCycle();
273  }
274 
275  delete fec_lor;
276 }
277 
278 void MatrixFreeAuxiliarySpace::SetupCG(
279  Operator& curlcurl_oper, Operator& conn,
280  int inner_cg_iterations)
281 {
282  MFEM_ASSERT(conn.Height() == curlcurl_oper.Width(),
283  "Operators don't match!");
284  matfree = new RAPOperator(conn, curlcurl_oper, conn);
285  MFEM_ASSERT(matfree->Height() == lor_pc->Height(),
286  "Operators don't match!");
287 
288  cg = new CGSolver(comm);
289  cg->SetOperator(*matfree);
290  cg->SetPreconditioner(*lor_pc);
291  if (inner_cg_iterations > 99)
292  {
293  cg->SetRelTol(1.e-14);
294  cg->SetMaxIter(100);
295  }
296  else
297  {
298  cg->SetRelTol(0.0);
299  cg->SetMaxIter(inner_cg_iterations);
300  }
301  cg->SetPrintLevel(-1);
302 
303  aspacewrapper = cg;
304 }
305 
306 void MatrixFreeAuxiliarySpace::SetupVCycle()
307 {
308  aspacewrapper = lor_pc;
309 }
310 
311 class ZeroWrapAMG : public Solver
312 {
313 public:
314 #ifdef MFEM_USE_AMGX
315  ZeroWrapAMG(HypreParMatrix& mat, Array<int>& ess_tdof_list_,
316  const bool useAmgX) :
317 #else
318  ZeroWrapAMG(HypreParMatrix& mat, Array<int>& ess_tdof_list_) :
319 #endif
320  Solver(mat.Height()), ess_tdof_list(ess_tdof_list_)
321  {
322 #ifdef MFEM_USE_AMGX
323  if (useAmgX)
324  {
325  const bool amgx_verbose = false;
326  AmgXSolver *amgx = new AmgXSolver(mat.GetComm(),
328  amgx_verbose);
329  amgx->SetOperator(mat);
330  amg_ = amgx;
331  }
332  else
333 #endif
334  {
335  HypreBoomerAMG *amg = new HypreBoomerAMG(mat);
336  amg->SetPrintLevel(0);
337  amg_ = amg;
338  }
339  }
340 
341  void Mult(const Vector& x, Vector& y) const
342  {
343  amg_->Mult(x, y);
344 
345  auto Y = y.HostReadWrite();
346  for (int k : ess_tdof_list)
347  {
348  Y[k] = 0.0;
349  }
350  }
351 
352  void SetOperator(const Operator&) { }
353 
354  ~ZeroWrapAMG()
355  {
356  delete amg_;
357  }
358 
359 private:
360  Solver *amg_ = NULL;
361  Array<int>& ess_tdof_list;
362 };
363 
364 void MatrixFreeAuxiliarySpace::SetupAMG(int system_dimension)
365 {
366  if (system_dimension == 0)
367  {
368  // boundary condition tweak for G-space solver
369 #ifdef MFEM_USE_AMGX
370  lor_pc = new ZeroWrapAMG(*lor_matrix, ess_tdof_list, useAmgX);
371 #else
372  lor_pc = new ZeroWrapAMG(*lor_matrix, ess_tdof_list);
373 #endif
374  }
375  else
376  {
377  // systems options for Pi-space solver
378 #ifdef MFEM_USE_AMGX
379  if (useAmgX)
380  {
381  const bool amgx_verbose = false;
382  AmgXSolver *amgx = new AmgXSolver(lor_matrix->GetComm(),
384  amgx_verbose);
385  amgx->SetOperator(*lor_matrix);
386  lor_pc = amgx;
387  }
388  else
389 #endif
390  {
391  HypreBoomerAMG* hpc = new HypreBoomerAMG(*lor_matrix);
392  hpc->SetSystemsOptions(system_dimension);
393  hpc->SetPrintLevel(0);
394  lor_pc = hpc;
395  }
396  }
397 }
398 
400 {
401  int rank;
402  MPI_Comm_rank(comm, &rank);
403 
404  y = 0.0;
405  aspacewrapper->Mult(x, y);
406  if (cg && rank == 0)
407  {
408  int q = cg->GetNumIterations();
409  inner_aux_iterations += q;
410  }
411 }
412 
414 {
415  delete lor_matrix;
416  delete lor_pc;
417  delete matfree;
418  if (lor_pc != aspacewrapper) { delete aspacewrapper; }
419  if (cg != aspacewrapper) { delete cg; }
420 }
421 
422 /* As an implementation note, a lot depends on the quality of the auxiliary
423  space solves. For high-contrast coefficients, and other difficult problems,
424  inner iteration counts may need to be increased. Boundary conditions can
425  matter as well (see DIAG_ZERO policy). */
427  ParBilinearForm& aform, Operator& oper, ParFiniteElementSpace& nd_fespace,
428  Coefficient* alpha_coeff, Coefficient* beta_coeff,
429  MatrixCoefficient* beta_mcoeff, Array<int>& ess_bdr,
430 #ifdef MFEM_USE_AMGX
431  bool useAmgX,
432 #endif
433  int inner_pi_iterations, int inner_g_iterations, Solver * nd_smoother) :
434  Solver(oper.Height())
435 {
436  int order = nd_fespace.GetFE(0)->GetOrder();
437  ParMesh *mesh = nd_fespace.GetParMesh();
438  int dim = mesh->Dimension();
439 
440  // smoother
441  Array<int> ess_tdof_list;
442  nd_fespace.GetEssentialTrueDofs(ess_bdr, ess_tdof_list);
443  if (nd_smoother)
444  {
445  smoother = nd_smoother;
446  }
447  else
448  {
449  const double scale = 0.25;
450  smoother = new OperatorJacobiSmoother(aform, ess_tdof_list, scale);
451  }
452 
453  // get H1 space
454  FiniteElementCollection *h1_fec = new H1_FECollection(order, dim);
455  h1_fespace = new ParFiniteElementSpace(mesh, h1_fec);
456  h1_fespace_d = new ParFiniteElementSpace(mesh, h1_fec, dim, Ordering::byVDIM);
457 
458  // build G operator
459  pa_grad = new ParDiscreteLinearOperator(h1_fespace, &nd_fespace);
462  pa_grad->Assemble();
463  pa_grad->FormRectangularSystemMatrix(Gradient);
464 
465  // build Pi operator
466  pa_interp = new ParDiscreteLinearOperator(h1_fespace_d, &nd_fespace);
469  pa_interp->Assemble();
470  pa_interp->FormRectangularSystemMatrix(Pi);
471 
472  // build LOR space
473  ParMesh mesh_lor = ParMesh::MakeRefined(*mesh, order, BasisType::GaussLobatto);
474 
475  // build G space solver
476  Gspacesolver = new MatrixFreeAuxiliarySpace(mesh_lor, beta_coeff,
477  beta_mcoeff, ess_bdr, oper,
478  *Gradient,
479 #ifdef MFEM_USE_AMGX
480  useAmgX,
481 #endif
482  inner_g_iterations);
483 
484  // build Pi space solver
485  Pispacesolver = new MatrixFreeAuxiliarySpace(mesh_lor, alpha_coeff,
486  beta_coeff, beta_mcoeff,
487  ess_bdr, oper, *Pi,
488 #ifdef MFEM_USE_AMGX
489  useAmgX,
490 #endif
491  inner_pi_iterations);
492 
493  general_ams = new GeneralAMS(oper, *Pi, *Gradient, *Pispacesolver,
494  *Gspacesolver, *smoother, ess_tdof_list);
495 
496  delete h1_fec;
497 }
498 
500 {
501  delete smoother;
502  delete pa_grad;
503  delete pa_interp;
504  delete Gspacesolver;
505  delete Pispacesolver;
506  delete general_ams;
507  delete h1_fespace;
508  delete h1_fespace_d;
509 }
510 
511 } // namespace mfem
512 
513 #endif
virtual void FormRectangularSystemMatrix(OperatorHandle &A)
Return in A a parallel (on truedofs) version of this operator.
int Size() const
Return the logical size of the array.
Definition: array.hpp:138
virtual void GetEssentialTrueDofs(const Array< int > &bdr_attr_is_ess, Array< int > &ess_tdof_list, int component=-1)
Definition: pfespace.cpp:1032
Conjugate gradient method.
Definition: solvers.hpp:465
MPI_Comm GetComm() const
MPI communicator.
Definition: hypre.hpp:528
Solver(int s=0, bool iter_mode=false)
Initialize a square Solver with size s.
Definition: operator.hpp:661
ParMesh * GetParMesh() const
Definition: pfespace.hpp:277
static ParMesh MakeRefined(ParMesh &orig_mesh, int ref_factor, int ref_type)
Create a uniformly refined (by any factor) version of orig_mesh.
Definition: pmesh.cpp:1366
void Mult(const Table &A, const Table &B, Table &C)
C = A * B (as boolean matrices)
Definition: table.cpp:472
int Width() const
Get the width (size of input) of the Operator. Synonym with NumCols().
Definition: operator.hpp:72
int GetOrder() const
Returns the order of the finite element. In the case of anisotropic orders, returns the maximum order...
Definition: fe_base.hpp:327
int Size() const
Returns the size of the vector.
Definition: vector.hpp:199
virtual void Mult(const Vector &x, Vector &y) const =0
Operator application: y=A(x).
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
void SetAssemblyLevel(AssemblyLevel assembly_level)
Set the desired assembly level. The default is AssemblyLevel::FULL.
void AddDomainInterpolator(DiscreteInterpolator *di)
Adds a domain interpolator. Assumes ownership of di.
HypreParMatrix * ParallelAssemble()
Returns the matrix assembled on the true dofs, i.e. P^t A P.
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 ...
Definition: operator.hpp:93
void SetDiagonalPolicy(DiagonalPolicy policy)
Sets diagonal policy used upon construction of the linear system.
virtual const FiniteElement * GetFE(int i) const
Definition: pfespace.cpp:535
int Height() const
Get the height (size of output) of the Operator. Synonym with NumRows().
Definition: operator.hpp:66
Jacobi smoothing for a given bilinear form (no matrix necessary).
Definition: solvers.hpp:274
void EliminateEssentialBC(const Array< int > &bdr_attr_is_ess, const Vector &sol, Vector &rhs, DiagonalPolicy dpolicy=DIAG_ONE)
Eliminate essential boundary DOFs from the system.
Perform AMS cycle with generic Operator objects.
Definition: auxiliary.hpp:124
void Assemble(int skip_zeros=1)
Assemble the local matrix.
The operator x -&gt; R*A*P*x constructed through the actions of R^T, A and P.
Definition: operator.hpp:763
void Mult(const Vector &x, Vector &y) const
Operator application: y=A(x).
Definition: auxiliary.cpp:399
Set the diagonal value to one.
Definition: operator.hpp:49
int Dimension() const
Definition: mesh.hpp:999
void UsePrecomputedSparsity(int ps=1)
For scalar FE spaces, precompute the sparsity pattern of the matrix (assuming dense element matrices)...
Base class Coefficients that optionally depend on space and time. These are used by the BilinearFormI...
Definition: coefficient.hpp:39
Collection of finite elements from the same family in multiple dimensions. This class is used to matc...
Definition: fe_coll.hpp:26
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
virtual void Finalize(int skip_zeros=1)
Finalizes the matrix initialization.
Keep the diagonal value.
Definition: operator.hpp:50
DiagonalPolicy
Defines operator diagonal policy upon elimination of rows and/or columns.
Definition: operator.hpp:46
int dim
Definition: ex24.cpp:53
void AddDomainIntegrator(BilinearFormIntegrator *bfi)
Adds new Domain Integrator. Assumes ownership of bfi.
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
Arbitrary order H1-conforming (continuous) finite elements.
Definition: fe_coll.hpp:216
Base class for solvers.
Definition: operator.hpp:651
Abstract operator.
Definition: operator.hpp:24
virtual void Assemble(int skip_zeros=1)
Construct the internal matrix representation of the discrete linear operator.
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