FreeFOAM The Cross-Platform CFD Toolkit
Hosted by SourceForge:
Get FreeFOAM at SourceForge.net.
            Fast, secure and Free Open Source software downloads

lduMatrixATmul.C

Go to the documentation of this file.
00001 /*---------------------------------------------------------------------------*\
00002   =========                 |
00003   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
00004    \\    /   O peration     |
00005     \\  /    A nd           | Copyright (C) 1991-2010 OpenCFD Ltd.
00006      \\/     M anipulation  |
00007 -------------------------------------------------------------------------------
00008 License
00009     This file is part of OpenFOAM.
00010 
00011     OpenFOAM is free software: you can redistribute it and/or modify it
00012     under the terms of the GNU General Public License as published by
00013     the Free Software Foundation, either version 3 of the License, or
00014     (at your option) any later version.
00015 
00016     OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
00017     ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
00018     FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
00019     for more details.
00020 
00021     You should have received a copy of the GNU General Public License
00022     along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
00023 
00024 Description
00025     Multiply a given vector (second argument) by the matrix or its transpose
00026     and return the result in the first argument.
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     // Initialise the update of interfaced interfaces
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     // Update interface interfaces
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     // Initialise the update of interfaced interfaces
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     // Update interface interfaces
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     // Add the interface internal coefficients to diagonal
00186     // and the interface boundary coefficients to the sum-off-diagonal
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     // Parallel boundary initialisation.
00226     // Note: there is a change of sign in the coupled
00227     // interface update.  The reason for this is that the
00228     // internal coefficients are all located at the l.h.s. of
00229     // the matrix whereas the "implicit" coefficients on the
00230     // coupled boundaries are all created as if the
00231     // coefficient contribution is of a source-kind (i.e. they
00232     // have a sign as if they are on the r.h.s. of the matrix.
00233     // To compensate for this, it is necessary to turn the
00234     // sign of the contribution.
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     // Initialise the update of interfaced interfaces
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     // Update interface interfaces
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 // ************************ vim: set sw=4 sts=4 et: ************************ //
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines