FreeFOAM The Cross-Platform CFD Toolkit
DICSmoother.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 "DICSmoother.H"
28 
29 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
30 
31 namespace Foam
32 {
33  defineTypeNameAndDebug(DICSmoother, 0);
34 
35  lduMatrix::smoother::addsymMatrixConstructorToTable<DICSmoother>
37 }
38 
39 
40 // * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
41 
43 (
44  const word& fieldName,
45  const lduMatrix& matrix,
46  const FieldField<Field, scalar>& interfaceBouCoeffs,
47  const FieldField<Field, scalar>& interfaceIntCoeffs,
48  const lduInterfaceFieldPtrsList& interfaces
49 )
50 :
52  (
53  fieldName,
54  matrix,
55  interfaceBouCoeffs,
56  interfaceIntCoeffs,
57  interfaces
58  ),
59  rD_(matrix_.diag())
60 {
61  DICPreconditioner::calcReciprocalD(rD_, matrix_);
62 }
63 
64 
65 // * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
66 
68 (
70  const scalarField& source,
71  const direction cmpt,
72  const label nSweeps
73 ) const
74 {
75  const scalar* const __restrict__ rDPtr = rD_.begin();
76  const scalar* const __restrict__ upperPtr = matrix_.upper().begin();
77  const label* const __restrict__ uPtr = matrix_.lduAddr().upperAddr().begin();
78  const label* const __restrict__ lPtr = matrix_.lduAddr().lowerAddr().begin();
79 
80  // Temporary storage for the residual
81  scalarField rA(rD_.size());
82  scalar* __restrict__ rAPtr = rA.begin();
83 
84  for (label sweep=0; sweep<nSweeps; sweep++)
85  {
86  matrix_.residual
87  (
88  rA,
89  psi,
90  source,
91  interfaceBouCoeffs_,
92  interfaces_,
93  cmpt
94  );
95 
96  rA *= rD_;
97 
98  register label nFaces = matrix_.upper().size();
99  for (register label face=0; face<nFaces; face++)
100  {
101  register label u = uPtr[face];
102  rAPtr[u] -= rDPtr[u]*upperPtr[face]*rAPtr[lPtr[face]];
103  }
104 
105  register label nFacesM1 = nFaces - 1;
106  for (register label face=nFacesM1; face>=0; face--)
107  {
108  register label l = lPtr[face];
109  rAPtr[l] -= rDPtr[l]*upperPtr[face]*rAPtr[uPtr[face]];
110  }
111 
112  psi += rA;
113  }
114 }
115 
116 
117 // ************************ vim: set sw=4 sts=4 et: ************************ //