MFEM  v4.4.0
Finite element discretization library
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Pages
multigrid.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 "multigrid.hpp"
13 
14 namespace mfem
15 {
16 
18  : cycleType(CycleType::VCYCLE), preSmoothingSteps(1), postSmoothingSteps(1)
19 {}
20 
22  const Array<Solver*>& smoothers_,
23  const Array<Operator*>& prolongations_,
24  const Array<bool>& ownedOperators_,
25  const Array<bool>& ownedSmoothers_,
26  const Array<bool>& ownedProlongations_)
27  : Solver(operators_.Last()->NumRows()), cycleType(CycleType::VCYCLE),
28  preSmoothingSteps(1), postSmoothingSteps(1),
29  X(operators_.Size()), Y(X.Size()), R(X.Size()), Z(X.Size())
30 {
31  operators_.Copy(operators);
32  smoothers_.Copy(smoothers);
33  prolongations_.Copy(prolongations);
34  ownedOperators_.Copy(ownedOperators);
35  ownedSmoothers_.Copy(ownedSmoothers);
36  ownedProlongations_.Copy(ownedProlongations);
37 
38  for (int level = 0; level < operators.Size(); ++level)
39  {
40  X[level] = new Vector(operators[level]->NumRows());
41  *X[level] = 0.0;
42  Y[level] = new Vector(operators[level]->NumRows());
43  *Y[level] = 0.0;
44  R[level] = new Vector(operators[level]->NumRows());
45  *R[level] = 0.0;
46  Z[level] = new Vector(operators[level]->NumRows());
47  *Z[level] = 0.0;
48  }
49 }
50 
52 {
53  for (int i = 0; i < operators.Size(); ++i)
54  {
55  if (ownedOperators[i])
56  {
57  delete operators[i];
58  }
59  if (ownedSmoothers[i])
60  {
61  delete smoothers[i];
62  }
63  delete X[i];
64  delete Y[i];
65  delete R[i];
66  delete Z[i];
67  }
68 
69  for (int i = 0; i < prolongations.Size(); ++i)
70  {
71  if (ownedProlongations[i])
72  {
73  delete prolongations[i];
74  }
75  }
76 
77  operators.DeleteAll();
78  smoothers.DeleteAll();
79  prolongations.DeleteAll();
80  X.DeleteAll();
81  Y.DeleteAll();
82  R.DeleteAll();
83  Z.DeleteAll();
84 }
85 
86 void Multigrid::AddLevel(Operator* opr, Solver* smoother, bool ownOperator,
87  bool ownSmoother)
88 {
89  operators.Append(opr);
90  smoothers.Append(smoother);
91  ownedOperators.Append(ownOperator);
92  ownedSmoothers.Append(ownSmoother);
93  width = opr->Width();
94  height = opr->Height();
95 
96  X.Append(new Vector(height));
97  *X.Last() = 0.0;
98  Y.Append(new Vector(height));
99  *Y.Last() = 0.0;
100  R.Append(new Vector(height));
101  *R.Last() = 0.0;
102  Z.Append(new Vector(height));
103  *Z.Last() = 0.0;
104 }
105 
106 int Multigrid::NumLevels() const { return operators.Size(); }
107 
108 int Multigrid::GetFinestLevelIndex() const { return NumLevels() - 1; }
109 
111 {
112  return operators[level];
113 }
114 
116 {
117  return operators[level];
118 }
119 
121 {
122  return GetOperatorAtLevel(operators.Size() - 1);
123 }
124 
126 {
127  return GetOperatorAtLevel(operators.Size() - 1);
128 }
129 
131 {
132  return smoothers[level];
133 }
134 
136 {
137  return smoothers[level];
138 }
139 
140 void Multigrid::SetCycleType(CycleType cycleType_, int preSmoothingSteps_,
141  int postSmoothingSteps_)
142 {
143  cycleType = cycleType_;
144  preSmoothingSteps = preSmoothingSteps_;
145  postSmoothingSteps = postSmoothingSteps_;
146 }
147 
148 void Multigrid::Mult(const Vector& x, Vector& y) const
149 {
150  MFEM_ASSERT(NumLevels() > 0, "");
151  *X.Last() = x;
152  *Y.Last() = 0.0;
153  Cycle(GetFinestLevelIndex());
154  y = *Y.Last();
155 }
156 
158 {
159  MFEM_ABORT("SetOperator not supported in Multigrid");
160 }
161 
162 void Multigrid::SmoothingStep(int level, bool transpose) const
163 {
164  GetOperatorAtLevel(level)->Mult(*Y[level], *R[level]); // r = A x
165  subtract(*X[level], *R[level], *R[level]); // r = b - A x
166  if (transpose)
167  {
168  GetSmootherAtLevel(level)->MultTranspose(*R[level], *Z[level]); // z = S r
169  }
170  else
171  {
172  GetSmootherAtLevel(level)->Mult(*R[level], *Z[level]); // z = S r
173  }
174  add(*Y[level], 1.0, *Z[level], *Y[level]); // x = x + S (b - A x)
175 }
176 
177 void Multigrid::Cycle(int level) const
178 {
179  if (level == 0)
180  {
181  GetSmootherAtLevel(level)->Mult(*X[level], *Y[level]);
182  return;
183  }
184 
185  for (int i = 0; i < preSmoothingSteps; i++)
186  {
187  SmoothingStep(level, false);
188  }
189 
190  // Compute residual
191  GetOperatorAtLevel(level)->Mult(*Y[level], *R[level]);
192  subtract(*X[level], *R[level], *R[level]);
193 
194  // Restrict residual
195  GetProlongationAtLevel(level - 1)->MultTranspose(*R[level], *X[level - 1]);
196 
197  // Init zeros
198  *Y[level - 1] = 0.0;
199 
200  // Corrections
201  int corrections = 1;
203  {
204  corrections = 2;
205  }
206  for (int correction = 0; correction < corrections; ++correction)
207  {
208  Cycle(level - 1);
209  }
210 
211  // Prolongate
212  GetProlongationAtLevel(level - 1)->Mult(*Y[level - 1], *R[level]);
213 
214  // Add update
215  *Y[level] += *R[level];
216 
217  // Post-smooth
218  for (int i = 0; i < postSmoothingSteps; i++)
219  {
220  SmoothingStep(level, true);
221  }
222 }
223 
224 const Operator* Multigrid::GetProlongationAtLevel(int level) const
225 {
226  return prolongations[level];
227 }
228 
230 {
231  for (int i = 0; i < bfs.Size(); ++i)
232  {
233  delete bfs[i];
234  }
235 
236  bfs.DeleteAll();
237 
238  for (int i = 0; i < essentialTrueDofs.Size(); ++i)
239  {
240  delete essentialTrueDofs[i];
241  }
242 
243  essentialTrueDofs.DeleteAll();
244 }
245 
247  OperatorHandle& A,
248  Vector& X, Vector& B)
249 {
250  bfs.Last()->FormLinearSystem(*essentialTrueDofs.Last(), x, b, A, X, B);
251 }
252 
254  const Vector& b, Vector& x)
255 {
256  bfs.Last()->RecoverFEMSolution(X, b, x);
257 }
258 
259 const Operator* GeometricMultigrid::GetProlongationAtLevel(int level) const
260 {
261  return fespaces.GetProlongationAtLevel(level);
262 }
263 
264 } // namespace mfem
Array< Vector * > Z
Definition: multigrid.hpp:50
Array< Vector * > X
Definition: multigrid.hpp:47
int Width() const
Get the width (size of input) of the Operator. Synonym with NumCols().
Definition: operator.hpp:72
virtual ~GeometricMultigrid()
Destructor.
Definition: multigrid.cpp:229
Pointer to an Operator of a specified type.
Definition: handle.hpp:33
void Copy(Array &copy) const
Create a copy of the internal array to the provided copy.
Definition: array.hpp:856
Array< BilinearForm * > bfs
Definition: multigrid.hpp:124
Array< Vector * > Y
Definition: multigrid.hpp:48
virtual void Mult(const Vector &x, Vector &y) const =0
Operator application: y=A(x).
void add(const Vector &v1, const Vector &v2, Vector &v)
Definition: vector.cpp:300
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
Array< Array< int > * > essentialTrueDofs
Definition: multigrid.hpp:123
Array< Operator * > prolongations
Definition: multigrid.hpp:37
CycleType cycleType
Definition: multigrid.hpp:43
int Height() const
Get the height (size of output) of the Operator. Synonym with NumRows().
Definition: operator.hpp:66
int Append(const T &el)
Append element &#39;el&#39; to array, resize if necessary.
Definition: array.hpp:751
double b
Definition: lissajous.cpp:42
Array< Vector * > R
Definition: multigrid.hpp:49
const FiniteElementSpaceHierarchy & fespaces
Definition: multigrid.hpp:122
Array< bool > ownedProlongations
Definition: multigrid.hpp:41
virtual void SetOperator(const Operator &op) override
Not supported for multigrid.
Definition: multigrid.cpp:157
int NumRows() const
Get the number of rows (size of output) of the Operator. Synonym with Height().
Definition: operator.hpp:69
void subtract(const Vector &x, const Vector &y, Vector &z)
Definition: vector.cpp:447
virtual ~Multigrid()
Destructor.
Definition: multigrid.cpp:51
Array< Solver * > smoothers
Definition: multigrid.hpp:36
int height
Dimension of the output / number of rows in the matrix.
Definition: operator.hpp:27
void AddLevel(Operator *opr, Solver *smoother, bool ownOperator, bool ownSmoother)
Adds a level to the multigrid operator hierarchy.
Definition: multigrid.cpp:86
Solver * GetSmootherAtLevel(int level) const
Returns smoother at given level.
Definition: multigrid.cpp:130
const Operator * GetOperatorAtFinestLevel() const
Returns operator at finest level.
Definition: multigrid.cpp:120
virtual void Mult(const Vector &x, Vector &y) const override
Application of the multigrid as a preconditioner.
Definition: multigrid.cpp:148
void RecoverFineFEMSolution(const Vector &X, const Vector &b, Vector &x)
Recover the solution of a linear system formed with FormFineLinearSystem()
Definition: multigrid.cpp:253
int GetFinestLevelIndex() const
Returns the index of the finest level.
Definition: multigrid.cpp:108
void FormFineLinearSystem(Vector &x, Vector &b, OperatorHandle &A, Vector &X, Vector &B)
Definition: multigrid.cpp:246
int NumLevels() const
Returns the number of levels.
Definition: multigrid.cpp:106
const Operator * GetOperatorAtLevel(int level) const
Returns operator at given level.
Definition: multigrid.cpp:110
Array< Operator * > operators
Definition: multigrid.hpp:35
Vector data type.
Definition: vector.hpp:60
Multigrid()
Constructs an empty multigrid hierarchy.
Definition: multigrid.cpp:17
int postSmoothingSteps
Definition: multigrid.hpp:45
Array< bool > ownedSmoothers
Definition: multigrid.hpp:40
Base class for solvers.
Definition: operator.hpp:651
Operator * GetProlongationAtLevel(int level) const
Returns the prolongation operator from the finite element space at level to the finite element space ...
Abstract operator.
Definition: operator.hpp:24
Array< bool > ownedOperators
Definition: multigrid.hpp:39
int width
Dimension of the input / number of columns in the matrix.
Definition: operator.hpp:28
void SetCycleType(CycleType cycleType_, int preSmoothingSteps_, int postSmoothingSteps_)
Set cycle type and number of pre- and post-smoothing steps used by Mult.
Definition: multigrid.cpp:140