Moukalled F., Mangani L., Darwish M. The finite volume method in computational fluid dynamics. An advanced introduction with OpenFOAM and Matlab (811443), страница 56
Текст из файла (страница 56)
Thus the lduMatrix.H filereads (Listing 10.1)class lduMatrix{// private data//- LDU mesh referenceconst lduMesh& lduMesh_;//- Coefficients (not including interfaces)scalarField *lowerPtr_, *diagPtr_, *upperPtr_;public://- Abstract base-class for lduMatrix solversclass solver{protected:class smoother{protected:class preconditioner{protected:Listing 10.1 The three base classes (solver, smoother, and preconditioner) defined with thelduMatrix class10.5Computational Pointers353So each smoother, solver, and preconditioner has to be derived from these threebase classes. For example, the DILU preconditioner is declared as shown in Listing10.2,class DILUPreconditioner:public lduMatrix::preconditioner{Listing 10.2 Syntax used to declare the DILU preconditionerwhile the Conjugate Gradient solver is declared as shown in Listing 10.3.class PCG:public lduMatrix::solverListing 10.3 Syntax used to declare the PCG solverIn either case the preconditioner or the solver is evidently derived from the baseclass defined under the lduMatrix class.Having clarified the basic concepts and organizational structure of algebraicsolvers in OpenFOAM®, an example investigating the details of implementing thepreconditioned CG method is now provided.
The files are located in the directory“$FOAM_SRC/OpenFOAM/matrices/lduMatrix/solvers/PCG”.The class derived from the lduMatrix::solver class, as shown in Listing 10.3,defines the main member function “solve” using the script in Listing 10.4.// Member Functions//- Solve the matrix with this solvervirtual solverPerformance solve(scalarField& psi,const scalarField& source,const direction cmpt=0) const;Listing 10.4 Script used to define the member function “solve”35410Solving the System of Algebraic EquationsThe function “solve” implements the solution algorithm of the chosen linearalgebraic solver in file “PCG.C”. Recalling the preconditioned conjugate gradientalgorithm, its sequence of events are given by1.
Calculate rð0Þ ¼ b A/ð0Þ2. Calculate dð0Þ ¼ P1 rð0Þ ðnÞ T 1 ðnÞrP rðnÞ3. Calculate a ¼ TdðnÞ AdðnÞ4. Calculate /ðnþ1Þ ¼ /ðnÞ þ aðnÞ dðnÞ5. Calculate rðnþ1Þ ¼ rðnÞ aðnÞ AdðnÞ6. If solution has converged stop ðnþ1Þ T 1 ðnþ1ÞrP rðnþ1Þ7. Calculate b¼T 1 ðnÞðnÞðr Þ P r8. Calculate dðnþ1Þ ¼ P1 rðnþ1Þ þ bðnþ1Þ dðnÞ9. Go to step 3The algorithm is directly implemented in “PCG.C” following the same procedure as described next.In step 1, depicted in Listing 10.5, the residual is evaluated and stored in the rAvariable while the wA variable stores the matrix-solution product A/ð0Þ .// --- Calculate A.psimatrix_.Amul(wA, psi, interfaceBouCoeffs_, interfaces_, cmpt);// --- Calculate initial residual fieldscalarField rA(source - wA);Listing 10.5 Script used to calculate the residualsStep 2 involves preconditioning.
Thus, first the type of preconditioner used isdefined with the object preconPtr as (Listing 10.6).// --- Select and construct the preconditionerautoPtr<lduMatrix::preconditioner> preconPtr =lduMatrix::preconditioner::New(*this,controlDict_);Listing 10.6 Defining the type of preconditioner usedThe constructor used is a generic one based on the base class and the “New”constructor. The preconditioner type is chosen from the dictionary at run time. Thenthe preconditioning operation, P1 rðnÞ , is applied to the residual rA (according tothe equation in step 2 of the algorithm).
The result is stored in the same variable wA10.5Computational Pointers355 Tto reduce memory usage and then used in the evaluation of rðnÞ P1 rðnÞ bysimply performing the scalar product of the two vectors wA and rA using thegSumProd function, as shown in Listing 10.7. Moreover the old value of wArAfrom the previous iteration ðn 1Þ is stored in the variable wArAold.wArAold = wArA;// --- Precondition residualpreconPtr->precondition(wA, rA, cmpt);// --- Update search directions:wArA = gSumProd(wA, rA, matrix().mesh().comm()); TListing 10.7 Syntax used to calculate rðnÞ P1 rðnÞ and to store its old valueFollowing the preconditioned conjugate gradient algorithm, steps 3, 4, and 5 areperformed as displayed in Listing 10.8.// --- Update preconditioned residualmatrix_.Amul(wA, pA, interfaceBouCoeffs_, interfaces_,cmpt);scalar wApA = gSumProd(wA, pA, matrix().mesh().comm());// --- Update solution and residual:scalar alpha = wArA/wApA;for (register label cell=0; cell<nCells; cell++){psiPtr[cell] += alpha*pAPtr[cell];rAPtr[cell] -= alpha*wAPtr[cell];}Listing 10.8 Script used to calculate a and to update the values of the dependent variable and theresidualsNow the variable wA stores the product AdðnÞ while pA represents the dðnÞ Tvector.
Again the product dðnÞ AdðnÞ is performed with the gSumProd functionand stored in the wApA variable. Once alpha is evaluated, the update of residualsand solution of steps 4 and 5 is performed in the for loop with the variables psiPtrand rAPtr shown in Listing 10.8 representing / and r, respectively.35610Solving the System of Algebraic EquationsThe iteration in the algorithm is completed by executing steps 7 and 8 using thescript in Listing 10.9.scalar beta = wArA/wArAold;for (register label cell=0; cell<nCells; cell++){pAPtr[cell] = wAPtr[cell] + beta*pAPtr[cell];}Listing 10.9 Script used for calculating new values for b and d to be used in the next iterationFor practical use in test cases, the definitions of the linear solver are made insidethe system directory in the fvSolution under the syntax solvers {} as shown inListing 10.10.solvers{T{solverpreconditionertolerancerelTolPCG;DIC;1e-06;0;}}Listing 10.10 Definition of the linear solverThe meaning of the various entries in Listing 10.10 are• solver: defines the type of solver (in this case PCG is the preconditioned conjugate gradient for symmetric matrices), with the various options being asfollows:– “PCG”: preconditioned conjugate gradient (for symmetric matrices only).– “PBiCG”: preconditioned biconjugate gradient (for asymmetric matricesonly).– “smoothSolver”: solver used only as a smoother to reduce residuals.– “GAMG”: generalised geometric algebraic multigrid.
It should be used forthe pressure equation on large grids.10.5Computational Pointers357• preconditioner: defines the type of the preconditioner to be used. The twoavailable options are– “DIC”: diagonal incomplete-cholesky preconditioner for symmetricmatrices.– “DILU”: diagonal incomplete-lower-upper preconditioner for asymmetricmatrices.• tolerance: the maximum allowable value of the absolute residual for the linearsolver to stop iterating.• relTol: the ratio between the initial residual and the actual residual for the linearsolver to stop iterating.Other examples are shown in Listing 10.11.solvers{T{solverpreconditionertolerancerelTolPBiCG;DILU;1e-06;0;}}solvers{T{solversmoothertolerancerelTolnSweepssmoothSolver;GaussSeidel;1e-8;0.1;1;}T{solverGAMG;tolerance1e-7;relTol0.01;smootherGaussSeidel;nPreSweeps0;nPostSweeps2;cacheAgglomeration on;agglomeratorfaceAreaPair;nCellsInCoarsestLevel 10;mergeLevels1;}}Listing 10.11 Examples of defining the linear solver35810Solving the System of Algebraic Equations10.6 ClosureThe chapter introduced the direct and iterative approaches for solving algebraicsystems of equations.
In each category a number of methods were described. Thealgebraic multigrid technique was also discussed. The next chapter will proceedwith the discretization of the conservation equation and will detail the discretizationof the convection term.10.7 ExercisesExercise 1In each of the following cases obtain the LU factorization of matrix A and use it tosolve the system of equation Ax = b by performing the backward and forwardsubstitution, i.e., Ly = b, Ux = y:0ðaÞ13:0B 2:0BA¼BB 4:0@ 3:04:00ðbÞ14:0B 7:0BA¼BB 2:0@ 4:03:00ðcÞ12:0B 4:0BA¼BB 3:0@ 6:04:05:012:02:03:06:06:06:015:05:007:03:04:09:05:015:04:0 CC5:0 CC5:0 A13:013:0B 12:0 CBCCb¼BB 13:0 C@ 21:0 A13:01:010:01:07:03:06:01:010:06:03:07:01:07:011:03:013:01:0 CC6:0 CC1:0 A14:0113:0B 13:0 CCBCb¼BB 16:0 C@ 13:0 A1:06:011:06:02:01:01:01:014:07:01:05:001:010:06:014:00 CC6:0 CC7:0 A12:0000113:0B 19:0 CBCCb¼BB 7:0 C@ 22:0 A30:0Exercise 2Using the Gauss-Seidel and Jacobi methods and starting with a zero initial guesssolve the following systems of equations while adopting as a stopping criteriamaxðkAx bkÞ\0:1:10.7ðaÞExercises06B0BA¼BB0@520ðbÞ36B0BB0A¼BB3B@35165503416433 CC0 CC0 A10054370124003600002136227B0BA¼@010270201265110 CC1 A26135922B 21 CB CCb¼BB 20 C@ 10 A22124B 14 CB CB 19 CCb¼BB 29 CB C@ 20 A80123 CC4 CC5 CC0 A4044002300ðcÞ0060210121B3 CCb¼B@ 10 A29Exercise 3For the systems of equations in Exercise 2 find the preconditioned matrix P for boththe Gauss-Seidel and Jacobi methods.
Use Eq. (10.81) with a zero initial guess toresolve the systems of equations subject to the same stopping criteria. Comparesolutions with those obtained in Exercise 2.Exercise 4Perform the ILU(0) factorization of the following M matrices:0ðaÞ6B0BB0BM¼BB0B3B@000ðbÞ6B3BB0BB0M¼BB0BB0B@40060000007000010006003500603000100601000060000032081300506010103CC0CC0CC0CC0A601100603500090340000060100CC0CC0CC0CC0CC0A6360100ðcÞ6B0BM¼BB0@30Solving the System of Algebraic Equations01010400650104 CC4 CC5 A1021065Exercise 5Based on the factorizations in Exercise 4 and knowing that the lower and upperparts of the factorization correspond to L and U, respectively, find for each of thecases the error matrix defined as R = M − P.Exercise 6Perform the DILU factorization of the following M matrices:0ðaÞ9B3BB0BM¼BB0B4B@040ðbÞ6B1BB1BB2M¼BB0BB1B@244751050260100200ðcÞ6B1BB3M¼BB4B@330010102000061002153600101CC4CC0CC1CC3A8525016220104131430084010410081103430280443410061263440308104033700504180101CC0CC0CC0A6100CC0CC0CC0CC1CC0A6Exercise 7Based on the diagonal factorizations in exercise 6 and using Eq.
![](https://s.studizba.com/z.php?f=/templates/si/i/noavatar.png&w=100&h=100&t=1"){kind=avatar}
