00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021
00022
00023
00024
00025
00026
00027
00028
00029
00030 #include <OpenFOAM/lduMatrix.H>
00031
00032
00033
00034 void Foam::lduMatrix::Amul
00035 (
00036 scalarField& Apsi,
00037 const tmp<scalarField>& tpsi,
00038 const FieldField<Field, scalar>& interfaceBouCoeffs,
00039 const lduInterfaceFieldPtrsList& interfaces,
00040 const direction cmpt
00041 ) const
00042 {
00043 scalar* __restrict__ ApsiPtr = Apsi.begin();
00044
00045 const scalarField& psi = tpsi();
00046 const scalar* const __restrict__ psiPtr = psi.begin();
00047
00048 const scalar* const __restrict__ diagPtr = diag().begin();
00049
00050 const label* const __restrict__ uPtr = lduAddr().upperAddr().begin();
00051 const label* const __restrict__ lPtr = lduAddr().lowerAddr().begin();
00052
00053 const scalar* const __restrict__ upperPtr = upper().begin();
00054 const scalar* const __restrict__ lowerPtr = lower().begin();
00055
00056
00057 initMatrixInterfaces
00058 (
00059 interfaceBouCoeffs,
00060 interfaces,
00061 psi,
00062 Apsi,
00063 cmpt
00064 );
00065
00066 register const label nCells = diag().size();
00067 for (register label cell=0; cell<nCells; cell++)
00068 {
00069 ApsiPtr[cell] = diagPtr[cell]*psiPtr[cell];
00070 }
00071
00072
00073 register const label nFaces = upper().size();
00074
00075 for (register label face=0; face<nFaces; face++)
00076 {
00077 ApsiPtr[uPtr[face]] += lowerPtr[face]*psiPtr[lPtr[face]];
00078 ApsiPtr[lPtr[face]] += upperPtr[face]*psiPtr[uPtr[face]];
00079 }
00080
00081
00082 updateMatrixInterfaces
00083 (
00084 interfaceBouCoeffs,
00085 interfaces,
00086 psi,
00087 Apsi,
00088 cmpt
00089 );
00090
00091 tpsi.clear();
00092 }
00093
00094
00095 void Foam::lduMatrix::Tmul
00096 (
00097 scalarField& Tpsi,
00098 const tmp<scalarField>& tpsi,
00099 const FieldField<Field, scalar>& interfaceIntCoeffs,
00100 const lduInterfaceFieldPtrsList& interfaces,
00101 const direction cmpt
00102 ) const
00103 {
00104 scalar* __restrict__ TpsiPtr = Tpsi.begin();
00105
00106 const scalarField& psi = tpsi();
00107 const scalar* const __restrict__ psiPtr = psi.begin();
00108
00109 const scalar* const __restrict__ diagPtr = diag().begin();
00110
00111 const label* const __restrict__ uPtr = lduAddr().upperAddr().begin();
00112 const label* const __restrict__ lPtr = lduAddr().lowerAddr().begin();
00113
00114 const scalar* const __restrict__ lowerPtr = lower().begin();
00115 const scalar* const __restrict__ upperPtr = upper().begin();
00116
00117
00118 initMatrixInterfaces
00119 (
00120 interfaceIntCoeffs,
00121 interfaces,
00122 psi,
00123 Tpsi,
00124 cmpt
00125 );
00126
00127 register const label nCells = diag().size();
00128 for (register label cell=0; cell<nCells; cell++)
00129 {
00130 TpsiPtr[cell] = diagPtr[cell]*psiPtr[cell];
00131 }
00132
00133 register const label nFaces = upper().size();
00134 for (register label face=0; face<nFaces; face++)
00135 {
00136 TpsiPtr[uPtr[face]] += upperPtr[face]*psiPtr[lPtr[face]];
00137 TpsiPtr[lPtr[face]] += lowerPtr[face]*psiPtr[uPtr[face]];
00138 }
00139
00140
00141 updateMatrixInterfaces
00142 (
00143 interfaceIntCoeffs,
00144 interfaces,
00145 psi,
00146 Tpsi,
00147 cmpt
00148 );
00149
00150 tpsi.clear();
00151 }
00152
00153
00154 void Foam::lduMatrix::sumA
00155 (
00156 scalarField& sumA,
00157 const FieldField<Field, scalar>& interfaceBouCoeffs,
00158 const lduInterfaceFieldPtrsList& interfaces
00159 ) const
00160 {
00161 scalar* __restrict__ sumAPtr = sumA.begin();
00162
00163 const scalar* __restrict__ diagPtr = diag().begin();
00164
00165 const label* __restrict__ uPtr = lduAddr().upperAddr().begin();
00166 const label* __restrict__ lPtr = lduAddr().lowerAddr().begin();
00167
00168 const scalar* __restrict__ lowerPtr = lower().begin();
00169 const scalar* __restrict__ upperPtr = upper().begin();
00170
00171 register const label nCells = diag().size();
00172 register const label nFaces = upper().size();
00173
00174 for (register label cell=0; cell<nCells; cell++)
00175 {
00176 sumAPtr[cell] = diagPtr[cell];
00177 }
00178
00179 for (register label face=0; face<nFaces; face++)
00180 {
00181 sumAPtr[uPtr[face]] += lowerPtr[face];
00182 sumAPtr[lPtr[face]] += upperPtr[face];
00183 }
00184
00185
00186
00187 forAll(interfaces, patchI)
00188 {
00189 if (interfaces.set(patchI))
00190 {
00191 const unallocLabelList& pa = lduAddr().patchAddr(patchI);
00192 const scalarField& pCoeffs = interfaceBouCoeffs[patchI];
00193
00194 forAll(pa, face)
00195 {
00196 sumAPtr[pa[face]] -= pCoeffs[face];
00197 }
00198 }
00199 }
00200 }
00201
00202
00203 void Foam::lduMatrix::residual
00204 (
00205 scalarField& rA,
00206 const scalarField& psi,
00207 const scalarField& source,
00208 const FieldField<Field, scalar>& interfaceBouCoeffs,
00209 const lduInterfaceFieldPtrsList& interfaces,
00210 const direction cmpt
00211 ) const
00212 {
00213 scalar* __restrict__ rAPtr = rA.begin();
00214
00215 const scalar* const __restrict__ psiPtr = psi.begin();
00216 const scalar* const __restrict__ diagPtr = diag().begin();
00217 const scalar* const __restrict__ sourcePtr = source.begin();
00218
00219 const label* const __restrict__ uPtr = lduAddr().upperAddr().begin();
00220 const label* const __restrict__ lPtr = lduAddr().lowerAddr().begin();
00221
00222 const scalar* const __restrict__ upperPtr = upper().begin();
00223 const scalar* const __restrict__ lowerPtr = lower().begin();
00224
00225
00226
00227
00228
00229
00230
00231
00232
00233
00234
00235
00236 FieldField<Field, scalar> mBouCoeffs(interfaceBouCoeffs.size());
00237
00238 forAll(mBouCoeffs, patchi)
00239 {
00240 if (interfaces.set(patchi))
00241 {
00242 mBouCoeffs.set(patchi, -interfaceBouCoeffs[patchi]);
00243 }
00244 }
00245
00246
00247 initMatrixInterfaces
00248 (
00249 mBouCoeffs,
00250 interfaces,
00251 psi,
00252 rA,
00253 cmpt
00254 );
00255
00256 register const label nCells = diag().size();
00257 for (register label cell=0; cell<nCells; cell++)
00258 {
00259 rAPtr[cell] = sourcePtr[cell] - diagPtr[cell]*psiPtr[cell];
00260 }
00261
00262
00263 register const label nFaces = upper().size();
00264
00265 for (register label face=0; face<nFaces; face++)
00266 {
00267 rAPtr[uPtr[face]] -= lowerPtr[face]*psiPtr[lPtr[face]];
00268 rAPtr[lPtr[face]] -= upperPtr[face]*psiPtr[uPtr[face]];
00269 }
00270
00271
00272 updateMatrixInterfaces
00273 (
00274 mBouCoeffs,
00275 interfaces,
00276 psi,
00277 rA,
00278 cmpt
00279 );
00280 }
00281
00282
00283 Foam::tmp<Foam::scalarField> Foam::lduMatrix::residual
00284 (
00285 const scalarField& psi,
00286 const scalarField& source,
00287 const FieldField<Field, scalar>& interfaceBouCoeffs,
00288 const lduInterfaceFieldPtrsList& interfaces,
00289 const direction cmpt
00290 ) const
00291 {
00292 tmp<scalarField> trA(new scalarField(psi.size()));
00293 residual(trA(), psi, source, interfaceBouCoeffs, interfaces, cmpt);
00294 return trA;
00295 }
00296
00297
00298