Go to the documentation of this file.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 #include "DILUSmoother.H"
00027 #include <OpenFOAM/DILUPreconditioner.H>
00028
00029
00030
00031 namespace Foam
00032 {
00033 defineTypeNameAndDebug(DILUSmoother, 0);
00034
00035 lduMatrix::smoother::addasymMatrixConstructorToTable<DILUSmoother>
00036 addDILUSmootherAsymMatrixConstructorToTable_;
00037 }
00038
00039
00040
00041
00042 Foam::DILUSmoother::DILUSmoother
00043 (
00044 const word& fieldName,
00045 const lduMatrix& matrix,
00046 const FieldField<Field, scalar>& interfaceBouCoeffs,
00047 const FieldField<Field, scalar>& interfaceIntCoeffs,
00048 const lduInterfaceFieldPtrsList& interfaces
00049 )
00050 :
00051 lduMatrix::smoother
00052 (
00053 fieldName,
00054 matrix,
00055 interfaceBouCoeffs,
00056 interfaceIntCoeffs,
00057 interfaces
00058 ),
00059 rD_(matrix_.diag())
00060 {
00061 DILUPreconditioner::calcReciprocalD(rD_, matrix_);
00062 }
00063
00064
00065
00066
00067 void Foam::DILUSmoother::smooth
00068 (
00069 scalarField& psi,
00070 const scalarField& source,
00071 const direction cmpt,
00072 const label nSweeps
00073 ) const
00074 {
00075 const scalar* const __restrict__ rDPtr = rD_.begin();
00076
00077 const label* const __restrict__ uPtr =
00078 matrix_.lduAddr().upperAddr().begin();
00079 const label* const __restrict__ lPtr =
00080 matrix_.lduAddr().lowerAddr().begin();
00081
00082 const scalar* const __restrict__ upperPtr = matrix_.upper().begin();
00083 const scalar* const __restrict__ lowerPtr = matrix_.lower().begin();
00084
00085
00086 scalarField rA(rD_.size());
00087 scalar* __restrict__ rAPtr = rA.begin();
00088
00089 for (label sweep=0; sweep<nSweeps; sweep++)
00090 {
00091 matrix_.residual
00092 (
00093 rA,
00094 psi,
00095 source,
00096 interfaceBouCoeffs_,
00097 interfaces_,
00098 cmpt
00099 );
00100
00101 rA *= rD_;
00102
00103 register label nFaces = matrix_.upper().size();
00104 for (register label face=0; face<nFaces; face++)
00105 {
00106 register label u = uPtr[face];
00107 rAPtr[u] -= rDPtr[u]*lowerPtr[face]*rAPtr[lPtr[face]];
00108 }
00109
00110 register label nFacesM1 = nFaces - 1;
00111 for (register label face=nFacesM1; face>=0; face--)
00112 {
00113 register label l = lPtr[face];
00114 rAPtr[l] -= rDPtr[l]*upperPtr[face]*rAPtr[uPtr[face]];
00115 }
00116
00117 psi += rA;
00118 }
00119 }
00120
00121
00122