FreeFOAM The Cross-Platform CFD Toolkit
GAMGSolverSolve.C
Go to the documentation of this file.
1 /*---------------------------------------------------------------------------*\
2  ========= |
3  \\ / F ield | OpenFOAM: The Open Source CFD Toolbox
4  \\ / O peration |
5  \\ / A nd | Copyright (C) 1991-2010 OpenCFD Ltd.
6  \\/ M anipulation |
7 -------------------------------------------------------------------------------
8 License
9  This file is part of OpenFOAM.
10 
11  OpenFOAM is free software: you can redistribute it and/or modify it
12  under the terms of the GNU General Public License as published by
13  the Free Software Foundation, either version 3 of the License, or
14  (at your option) any later version.
15 
16  OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
17  ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
18  FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
19  for more details.
20 
21  You should have received a copy of the GNU General Public License
22  along with OpenFOAM. If not, see <http://www.gnu.org/licenses/>.
23 
24 \*---------------------------------------------------------------------------*/
25 
26 #include "GAMGSolver.H"
27 #include <OpenFOAM/ICCG.H>
28 #include <OpenFOAM/BICCG.H>
29 #include <OpenFOAM/SubField.H>
30 
31 // * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
32 
34 (
36  const scalarField& source,
37  const direction cmpt
38 ) const
39 {
40  // Setup class containing solver performance data
41  lduMatrix::solverPerformance solverPerf(typeName, fieldName_);
42 
43  // Calculate A.psi used to calculate the initial residual
44  scalarField Apsi(psi.size());
45  matrix_.Amul(Apsi, psi, interfaceBouCoeffs_, interfaces_, cmpt);
46 
47  // Create the storage for the finestCorrection which may be used as a
48  // temporary in normFactor
49  scalarField finestCorrection(psi.size());
50 
51  // Calculate normalisation factor
52  scalar normFactor = this->normFactor(psi, source, Apsi, finestCorrection);
53 
54  if (debug >= 2)
55  {
56  Pout<< " Normalisation factor = " << normFactor << endl;
57  }
58 
59  // Calculate initial finest-grid residual field
60  scalarField finestResidual(source - Apsi);
61 
62  // Calculate normalised residual for convergence test
63  solverPerf.initialResidual() = gSumMag(finestResidual)/normFactor;
64  solverPerf.finalResidual() = solverPerf.initialResidual();
65 
66 
67  // Check convergence, solve if not converged
68  if (!solverPerf.checkConvergence(tolerance_, relTol_))
69  {
70  // Create coarse grid correction fields
71  PtrList<scalarField> coarseCorrFields;
72 
73  // Create coarse grid sources
74  PtrList<scalarField> coarseSources;
75 
76  // Create the smoothers for all levels
78 
79  // Initialise the above data structures
80  initVcycle(coarseCorrFields, coarseSources, smoothers);
81 
82  do
83  {
84  Vcycle
85  (
86  smoothers,
87  psi,
88  source,
89  Apsi,
90  finestCorrection,
91  finestResidual,
92  coarseCorrFields,
93  coarseSources,
94  cmpt
95  );
96 
97  // Calculate finest level residual field
98  matrix_.Amul(Apsi, psi, interfaceBouCoeffs_, interfaces_, cmpt);
99  finestResidual = source;
100  finestResidual -= Apsi;
101 
102  solverPerf.finalResidual() = gSumMag(finestResidual)/normFactor;
103 
104  if (debug >= 2)
105  {
106  solverPerf.print();
107  }
108  } while
109  (
110  ++solverPerf.nIterations() < maxIter_
111  && !(solverPerf.checkConvergence(tolerance_, relTol_))
112  );
113  }
114 
115  return solverPerf;
116 }
117 
118 
119 void Foam::GAMGSolver::Vcycle
120 (
121  const PtrList<lduMatrix::smoother>& smoothers,
122  scalarField& psi,
123  const scalarField& source,
124  scalarField& Apsi,
125  scalarField& finestCorrection,
126  scalarField& finestResidual,
127  PtrList<scalarField>& coarseCorrFields,
128  PtrList<scalarField>& coarseSources,
129  const direction cmpt
130 ) const
131 {
132  //debug = 2;
133 
134  const label coarsestLevel = matrixLevels_.size() - 1;
135 
136  // Restrict finest grid residual for the next level up
137  agglomeration_.restrictField(coarseSources[0], finestResidual, 0);
138 
139  if (debug >= 2 && nPreSweeps_)
140  {
141  Pout<< "Pre-smoothing scaling factors: ";
142  }
143 
144 
145  // Residual restriction (going to coarser levels)
146  for (label leveli = 0; leveli < coarsestLevel; leveli++)
147  {
148  // If the optional pre-smoothing sweeps are selected
149  // smooth the coarse-grid field for the restriced source
150  if (nPreSweeps_)
151  {
152  coarseCorrFields[leveli] = 0.0;
153 
154  smoothers[leveli + 1].smooth
155  (
156  coarseCorrFields[leveli],
157  coarseSources[leveli],
158  cmpt,
159  nPreSweeps_ + leveli
160  );
161 
163  (
164  Apsi,
165  coarseCorrFields[leveli].size()
166  );
167 
168  // Scale coarse-grid correction field
169  // but not on the coarsest level because it evaluates to 1
170  if (scaleCorrection_ && leveli < coarsestLevel - 1)
171  {
172  scalar sf = scalingFactor
173  (
174  const_cast<scalarField&>(ACf.operator const scalarField&()),
175  matrixLevels_[leveli],
176  coarseCorrFields[leveli],
177  interfaceLevelsBouCoeffs_[leveli],
178  interfaceLevels_[leveli],
179  coarseSources[leveli],
180  cmpt
181  );
182 
183  if (debug >= 2)
184  {
185  Pout<< sf << " ";
186  }
187 
188  coarseCorrFields[leveli] *= sf;
189  }
190 
191  // Correct the residual with the new solution
192  matrixLevels_[leveli].Amul
193  (
194  const_cast<scalarField&>(ACf.operator const scalarField&()),
195  coarseCorrFields[leveli],
196  interfaceLevelsBouCoeffs_[leveli],
197  interfaceLevels_[leveli],
198  cmpt
199  );
200 
201  coarseSources[leveli] -= ACf;
202  }
203 
204  // Residual is equal to source
205  agglomeration_.restrictField
206  (
207  coarseSources[leveli + 1],
208  coarseSources[leveli],
209  leveli + 1
210  );
211  }
212 
213  if (debug >= 2 && nPreSweeps_)
214  {
215  Pout<< endl;
216  }
217 
218 
219  // Solve Coarsest level with either an iterative or direct solver
220  solveCoarsestLevel
221  (
222  coarseCorrFields[coarsestLevel],
223  coarseSources[coarsestLevel]
224  );
225 
226 
227  if (debug >= 2)
228  {
229  Pout<< "Post-smoothing scaling factors: ";
230  }
231 
232  // Smoothing and prolongation of the coarse correction fields
233  // (going to finer levels)
234  for (label leveli = coarsestLevel - 1; leveli >= 0; leveli--)
235  {
236  // Create a field for the pre-smoothed correction field
237  // as a sub-field of the finestCorrection which is not
238  // currently being used
239  scalarField::subField preSmoothedCoarseCorrField
240  (
241  finestCorrection,
242  coarseCorrFields[leveli].size()
243  );
244 
245  // Only store the preSmoothedCoarseCorrField is pre-smoothing is used
246  if (nPreSweeps_)
247  {
248  preSmoothedCoarseCorrField.assign(coarseCorrFields[leveli]);
249  }
250 
251  agglomeration_.prolongField
252  (
253  coarseCorrFields[leveli],
254  coarseCorrFields[leveli + 1],
255  leveli + 1
256  );
257 
258  // Scale coarse-grid correction field
259  // but not on the coarsest level because it evaluates to 1
260  if (scaleCorrection_ && leveli < coarsestLevel - 1)
261  {
262  // Create A.psi for this coarse level as a sub-field of Apsi
264  (
265  Apsi,
266  coarseCorrFields[leveli].size()
267  );
268 
269  scalar sf = scalingFactor
270  (
271  const_cast<scalarField&>(ACf.operator const scalarField&()),
272  matrixLevels_[leveli],
273  coarseCorrFields[leveli],
274  interfaceLevelsBouCoeffs_[leveli],
275  interfaceLevels_[leveli],
276  coarseSources[leveli],
277  cmpt
278  );
279 
280 
281  if (debug >= 2)
282  {
283  Pout<< sf << " ";
284  }
285 
286  coarseCorrFields[leveli] *= sf;
287  }
288 
289  // Only add the preSmoothedCoarseCorrField is pre-smoothing is used
290  if (nPreSweeps_)
291  {
292  coarseCorrFields[leveli] += preSmoothedCoarseCorrField;
293  }
294 
295  smoothers[leveli + 1].smooth
296  (
297  coarseCorrFields[leveli],
298  coarseSources[leveli],
299  cmpt,
300  nPostSweeps_ + leveli
301  );
302  }
303 
304  // Prolong the finest level correction
305  agglomeration_.prolongField
306  (
307  finestCorrection,
308  coarseCorrFields[0],
309  0
310  );
311 
312  if (scaleCorrection_)
313  {
314  // Calculate finest level scaling factor
315  scalar fsf = scalingFactor
316  (
317  Apsi,
318  matrix_,
319  finestCorrection,
321  interfaces_,
322  finestResidual,
323  cmpt
324  );
325 
326  if (debug >= 2)
327  {
328  Pout<< fsf << endl;
329  }
330 
331  forAll(psi, i)
332  {
333  psi[i] += fsf*finestCorrection[i];
334  }
335  }
336  else
337  {
338  forAll(psi, i)
339  {
340  psi[i] += finestCorrection[i];
341  }
342  }
343 
344  smoothers[0].smooth
345  (
346  psi,
347  source,
348  cmpt,
349  nFinestSweeps_
350  );
351 }
352 
353 
354 void Foam::GAMGSolver::initVcycle
355 (
356  PtrList<scalarField>& coarseCorrFields,
357  PtrList<scalarField>& coarseSources,
358  PtrList<lduMatrix::smoother>& smoothers
359 ) const
360 {
361  coarseCorrFields.setSize(matrixLevels_.size());
362  coarseSources.setSize(matrixLevels_.size());
363  smoothers.setSize(matrixLevels_.size() + 1);
364 
365  // Create the smoother for the finest level
366  smoothers.set
367  (
368  0,
370  (
371  fieldName_,
372  matrix_,
373  interfaceBouCoeffs_,
374  interfaceIntCoeffs_,
375  interfaces_,
376  controlDict_
377  )
378  );
379 
380  forAll (matrixLevels_, leveli)
381  {
382  coarseCorrFields.set
383  (
384  leveli,
385  new scalarField
386  (
387  agglomeration_.meshLevel(leveli + 1).lduAddr().size()
388  )
389  );
390 
391  coarseSources.set
392  (
393  leveli,
394  new scalarField
395  (
396  agglomeration_.meshLevel(leveli + 1).lduAddr().size()
397  )
398  );
399 
400  smoothers.set
401  (
402  leveli + 1,
404  (
405  fieldName_,
406  matrixLevels_[leveli],
407  interfaceLevelsBouCoeffs_[leveli],
408  interfaceLevelsIntCoeffs_[leveli],
409  interfaceLevels_[leveli],
410  controlDict_
411  )
412  );
413  }
414 }
415 
416 
417 void Foam::GAMGSolver::solveCoarsestLevel
418 (
419  scalarField& coarsestCorrField,
420  const scalarField& coarsestSource
421 ) const
422 {
423  if (directSolveCoarsest_)
424  {
425  coarsestCorrField = coarsestSource;
426  coarsestLUMatrixPtr_->solve(coarsestCorrField);
427  }
428  else
429  {
430  const label coarsestLevel = matrixLevels_.size() - 1;
431  coarsestCorrField = 0;
432  lduMatrix::solverPerformance coarseSolverPerf;
433 
434  if (matrixLevels_[coarsestLevel].asymmetric())
435  {
436  coarseSolverPerf = BICCG
437  (
438  "coarsestLevelCorr",
439  matrixLevels_[coarsestLevel],
440  interfaceLevelsBouCoeffs_[coarsestLevel],
441  interfaceLevelsIntCoeffs_[coarsestLevel],
442  interfaceLevels_[coarsestLevel],
443  tolerance_,
444  relTol_
445  ).solve
446  (
447  coarsestCorrField,
448  coarsestSource
449  );
450  }
451  else
452  {
453  coarseSolverPerf = ICCG
454  (
455  "coarsestLevelCorr",
456  matrixLevels_[coarsestLevel],
457  interfaceLevelsBouCoeffs_[coarsestLevel],
458  interfaceLevelsIntCoeffs_[coarsestLevel],
459  interfaceLevels_[coarsestLevel],
460  tolerance_,
461  relTol_
462  ).solve
463  (
464  coarsestCorrField,
465  coarsestSource
466  );
467  }
468 
469  if (debug >= 2)
470  {
471  coarseSolverPerf.print();
472  }
473  }
474 }
475 
476 
477 // ************************ vim: set sw=4 sts=4 et: ************************ //