OpenFOAM logo
Open Source CFD Toolkit

lduMatrix.H

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-2005 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 the
00013     Free Software Foundation; either version 2 of the License, or (at your
00014     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, write to the Free Software Foundation,
00023     Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
00024 
00025 Class
00026     lduMatrix
00027 
00028 Description
00029     lduMatrix is a general matrix class in which the coefficients are
00030     stored as three arrays, one for the upper triangle, one for the
00031     lower triangle and a third for the diagonal.  Addressing arrays must
00032     be supplied for the upper and lower triangles.  Diagonal, ICCG and BCG
00033     solvers are supplied and selected automatically by the solve member
00034     function according the the particular matrix structure.
00035 
00036     It would be better if this class were organised as a hierachy starting
00037     from an empty matrix, then deriving diagonal, symmetric and asymmetric
00038     matrices.  A job for the future.
00039 
00040 SourceFiles
00041     lduMatrixATmul.C
00042     lduMatrix.C
00043     lduMatrixHoperations.C
00044     lduMatrixOperations.C
00045     lduMatrixSolver.C
00046     lduMatrixTests.C
00047     lduMatrixUpdateCoupledInterfaces.C
00048 
00049 \*---------------------------------------------------------------------------*/
00050 
00051 #ifndef lduMatrix_H
00052 #define lduMatrix_H
00053 
00054 #include "primitiveFields.H"
00055 #include "FieldField.H"
00056 #include "lduAddressing.H"
00057 #include "processorTopology.H"
00058 #include "lduCoupledInterfacePtrsList.H"
00059 #include "typeInfo.H"
00060 #include "autoPtr.H"
00061 #include "runTimeSelectionTables.H"
00062 
00063 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
00064 
00065 namespace Foam
00066 {
00067 
00068 /*---------------------------------------------------------------------------*\
00069                            Class lduMatrix Declaration
00070 \*---------------------------------------------------------------------------*/
00071 
00072 class lduMatrix
00073 {
00074 
00075 public:
00076 
00077     //- Class returned by the solver
00078     //  containing performance statistics
00079     class solverPerformance
00080     {
00081         word   solverName_;
00082         word   fieldName_;
00083         scalar initialResidual_;
00084         scalar finalResidual_;
00085         label  noIterations_;
00086         bool   converged_;
00087         bool   singular_;
00088 
00089 
00090     public:
00091 
00092         // Constructors
00093 
00094             solverPerformance()
00095             :
00096                 initialResidual_(0),
00097                 finalResidual_(0),
00098                 noIterations_(0),
00099                 converged_(false),
00100                 singular_(false)
00101             {}
00102 
00103 
00104             solverPerformance
00105             (
00106                 const word&  solverName,
00107                 const word&  fieldName,
00108                 const scalar iRes = 0,
00109                 const scalar fRes = 0,
00110                 const label  nIter = 0,
00111                 const bool   converged = false,
00112                 const bool   singular = false
00113             )
00114             :
00115                 solverName_(solverName),
00116                 fieldName_(fieldName),
00117                 initialResidual_(iRes),
00118                 finalResidual_(fRes),
00119                 noIterations_(nIter),
00120                 converged_(converged),
00121                 singular_(singular)
00122             {}
00123 
00124 
00125         // Member functions
00126 
00127             //- Return solver name
00128             const word& solverName() const
00129             {
00130                 return solverName_;
00131             }
00132 
00133             //- Return initial residual
00134             scalar initialResidual() const
00135             {
00136                 return initialResidual_;
00137             }
00138 
00139             //- Return initial residual
00140             scalar& initialResidual()
00141             {
00142                 return initialResidual_;
00143             }
00144 
00145 
00146             //- Return final residual
00147             scalar finalResidual() const
00148             {
00149                 return finalResidual_;
00150             }
00151 
00152             //- Return final residual
00153             scalar& finalResidual()
00154             {
00155                 return finalResidual_;
00156             }
00157 
00158 
00159             //- Return number of iterations
00160             label nIterations() const
00161             {
00162                 return noIterations_;
00163             }
00164 
00165             //- Return number of iterations
00166             label& nIterations()
00167             {
00168                 return noIterations_;
00169             }
00170 
00171 
00172             //- Has the solver converged?
00173             bool converged() const
00174             {
00175                 return converged_;
00176             }
00177 
00178             //- Is the matrix singular?
00179             bool singular() const
00180             {
00181                 return singular_;
00182             }
00183 
00184             //- Convergence test
00185             bool checkConvergence
00186             (
00187                 const scalar tolerance,
00188                 const scalar relTolerance
00189             );
00190 
00191             //- Singularity test
00192             bool checkSingularity(const scalar residual);
00193 
00194             //- Print summary of solver performance
00195             void print() const;
00196     };
00197 
00198 
00199     class solver
00200     {
00201 
00202     protected:
00203 
00204         // Protected data
00205 
00206             word fieldName_;
00207             scalarField& psi_;
00208             const lduMatrix& matrix_;
00209             const scalarField& source_;
00210             const FieldField<Field, scalar>& coupleBouCoeffs_;
00211             const FieldField<Field, scalar>& coupleIntCoeffs_;
00212             const lduCoupledInterfacePtrsList& interfaces_;
00213             const direction cmpt_;
00214 
00215 
00216     public:
00217 
00218         //- Runtime type information
00219         virtual const word& type() const = 0;
00220 
00221 
00222         // Declare run-time constructor selection tables
00223 
00224             declareRunTimeSelectionTable
00225             (
00226                 autoPtr,
00227                 solver,
00228                 symMatrix,
00229                 (
00230                     const word& fieldName,
00231                     scalarField& psi,
00232                     lduMatrix& matrix,
00233                     const scalarField& source,
00234                     const FieldField<Field, scalar>& coupleBouCoeffs,
00235                     const FieldField<Field, scalar>& coupleIntCoeffs,
00236                     const lduCoupledInterfacePtrsList& interfaces,
00237                     const direction cmpt,
00238                     Istream& solverData
00239                 ),
00240                 (
00241                     fieldName,
00242                     psi,
00243                     matrix,
00244                     source,
00245                     coupleBouCoeffs,
00246                     coupleIntCoeffs,
00247                     interfaces,
00248                     cmpt,
00249                     solverData
00250                 )
00251             );
00252 
00253             declareRunTimeSelectionTable
00254             (
00255                 autoPtr,
00256                 solver,
00257                 asymMatrix,
00258                 (
00259                     const word& fieldName,
00260                     scalarField& psi,
00261                     lduMatrix& matrix,
00262                     const scalarField& source,
00263                     const FieldField<Field, scalar>& coupleBouCoeffs,
00264                     const FieldField<Field, scalar>& coupleIntCoeffs,
00265                     const lduCoupledInterfacePtrsList& interfaces,
00266                     const direction cmpt,
00267                     Istream& solverData
00268                 ),
00269                 (
00270                     fieldName,
00271                     psi,
00272                     matrix,
00273                     source,
00274                     coupleBouCoeffs,
00275                     coupleIntCoeffs,
00276                     interfaces,
00277                     cmpt,
00278                     solverData
00279                 )
00280             );
00281 
00282 
00283         // Constructors
00284 
00285             solver
00286             (
00287                 const word& fieldName,
00288                 scalarField& psi,
00289                 const lduMatrix& matrix,
00290                 const scalarField& source,
00291                 const FieldField<Field, scalar>& coupleBouCoeffs,
00292                 const FieldField<Field, scalar>& coupleIntCoeffs,
00293                 const lduCoupledInterfacePtrsList& interfaces,
00294                 const direction cmpt
00295             )
00296             :
00297                 fieldName_(fieldName),
00298                 psi_(psi),
00299                 matrix_(matrix),
00300                 source_(source),
00301                 coupleBouCoeffs_(coupleBouCoeffs),
00302                 coupleIntCoeffs_(coupleIntCoeffs),
00303                 interfaces_(interfaces),
00304                 cmpt_(cmpt)
00305             {}
00306 
00307 
00308         // Selectors
00309 
00310             //- Return a new solver
00311             static autoPtr<solver> New
00312             (
00313                 const word& fieldName,
00314                 scalarField& psi,
00315                 lduMatrix& matrix,
00316                 const scalarField& source,
00317                 const FieldField<Field, scalar>& coupleBouCoeffs,
00318                 const FieldField<Field, scalar>& coupleIntCoeffs,
00319                 const lduCoupledInterfacePtrsList& interfaces,
00320                 const direction cmpt,
00321                 Istream& solverData
00322             );
00323 
00324 
00325         // Destructor
00326 
00327             virtual ~solver()
00328             {}
00329 
00330 
00331         // Member functions
00332 
00333             virtual solverPerformance solve() = 0;
00334     };
00335 
00336 
00337 private:
00338 
00339     // private data
00340 
00341         // L-U addressing
00342         const lduAddressing& lduAddr_;
00343 
00344         // lduCoupledInterface evaluation schedule
00345         const patchScheduleList& lduCoupledInterfaceSchedule_;
00346 
00347         //- Face matrix elements
00348         scalarField *lowerPtr_, *diagPtr_, *upperPtr_;
00349 
00350 
00351 public:
00352 
00353     // Static data
00354 
00355         // Declare name of the class and it's debug switch
00356         ClassName("lduMatrix");
00357 
00358         //- Large scalar for the use in solvers
00359         static const scalar great_;
00360 
00361         //- Small scalar for the use in solvers
00362         static const scalar small_;
00363 
00364 
00365     // Constructors
00366 
00367         //- Construct given a field to solve, L-U matrix addressing
00368         //  and tolerence
00369         lduMatrix
00370         (
00371             const lduAddressing& ldu,
00372             const patchScheduleList& lduCoupledInterfaceSchedule
00373         );
00374 
00375         //- Construct as copy
00376         lduMatrix(const lduMatrix&);
00377 
00378         //- Construct as copy or re-use as specified.
00379         lduMatrix(lduMatrix&, bool reUse);
00380 
00381         //- Construct from Istream
00382         lduMatrix
00383         (
00384             const lduAddressing& ldu,
00385             const patchScheduleList& lduCoupledInterfaceSchedule,
00386             Istream&
00387         );
00388 
00389 
00390     // Destructor
00391 
00392         ~lduMatrix();
00393 
00394 
00395     // Member functions
00396 
00397         // Access
00398 
00399             scalarField& lower();
00400             scalarField& diag();
00401             scalarField& upper();
00402 
00403             const scalarField& lower() const;
00404             const scalarField& diag() const;
00405             const scalarField& upper() const;
00406 
00407             const lduAddressing& lduAddr() const
00408             {
00409                 return lduAddr_;
00410             }
00411 
00412             const patchScheduleList& lduCoupledInterfaceSchedule() const
00413             {
00414                 return lduCoupledInterfaceSchedule_;
00415             }
00416 
00417             bool hasDiag() const
00418             {
00419                 return (diagPtr_);
00420             }
00421 
00422             bool hasUpper() const
00423             {
00424                 return (upperPtr_);
00425             }
00426 
00427             bool hasLower() const
00428             {
00429                 return (lowerPtr_);
00430             }
00431 
00432             bool diagonal() const
00433             {
00434                 return (diagPtr_ && !lowerPtr_ && !upperPtr_);
00435             }
00436 
00437             bool symmetric() const
00438             {
00439                 return (diagPtr_ && (!lowerPtr_ && upperPtr_));
00440             }
00441 
00442             bool asymmetric() const
00443             {
00444                 return (diagPtr_ && lowerPtr_ && upperPtr_);
00445             }
00446 
00447 
00448         // operations
00449 
00450             void sumDiag();
00451             void negSumDiag();
00452 
00453             //- Relax matrix (for steady-state solution)
00454             //  alpha=1 : diagonally equal
00455             //  alpha<1 :    ,,      dominant
00456             void relax
00457             (
00458                 const FieldField<Field, scalar>& intCoeffsCmptAvg,
00459                 const FieldField<Field, scalar>& magCoupleBouCoeffs,
00460                 const lduCoupledInterfacePtrsList& interfaces,
00461                 const scalar alpha=1.0
00462             );
00463 
00464 
00465             //- Matrix multiplication with updated coupled interfaces.
00466             void Amul
00467             (
00468                 scalarField&,
00469                 const tmp<scalarField>&,
00470                 const FieldField<Field, scalar>&,
00471                 const lduCoupledInterfacePtrsList&,
00472                 const direction cmpt
00473             ) const;
00474     
00475             //- Matrix transpose multiplication with updated coupled interfaces.
00476             void Tmul
00477             (
00478                 scalarField&,
00479                 const tmp<scalarField>&,
00480                 const FieldField<Field, scalar>&,
00481                 const lduCoupledInterfacePtrsList&,
00482                 const direction cmpt
00483             ) const;
00484 
00485 
00486             //- Initialise the update of coupled interfaces
00487             //  for matrix operations
00488             void initMatrixInterfaces
00489             (
00490                 const FieldField<Field, scalar>& coupleCoeffs,
00491                 const lduCoupledInterfacePtrsList& interfaces,
00492                 const scalarField& psiif,
00493                 scalarField& result,
00494                 const direction cmpt
00495             ) const;
00496     
00497             //- Update coupled interfaces for matrix operations
00498             void updateMatrixInterfaces
00499             (
00500                 const FieldField<Field, scalar>& coupleCoeffs,
00501                 const lduCoupledInterfacePtrsList& interfaces,
00502                 const scalarField& psiif,
00503                 scalarField& result,
00504                 const direction cmpt
00505             ) const;
00506     
00507 
00508             tmp<scalarField> residual
00509             (
00510                 const scalarField& psi,
00511                 const FieldField<Field, scalar>& coupleBouCoeffs,
00512                 const lduCoupledInterfacePtrsList& interfaces,
00513                 const direction cmpt
00514             ) const;
00515 
00516             tmp<scalarField> residual
00517             (
00518                 const scalarField& psi,
00519                 const scalarField& source,
00520                 const FieldField<Field, scalar>& coupleBouCoeffs,
00521                 const lduCoupledInterfacePtrsList& interfaces,
00522                 const direction cmpt
00523             ) const;
00524 
00525 
00526             template<class Type>
00527             tmp<Field<Type> > H(const Field<Type>&) const;
00528 
00529             template<class Type>
00530             tmp<Field<Type> > H(const tmp<Field<Type> >&) const;
00531 
00532             template<class Type>
00533             tmp<Field<Type> > faceH(const Field<Type>&) const;
00534 
00535             template<class Type>
00536             tmp<Field<Type> > faceH(const tmp<Field<Type> >&) const;
00537 
00538 
00539     // Member operators
00540 
00541         void operator=(const lduMatrix&);
00542 
00543         void negate();
00544 
00545         void operator+=(const lduMatrix&);
00546         void operator-=(const lduMatrix&);
00547 
00548         void operator*=(const scalarField&);
00549         void operator*=(scalar);
00550 
00551 
00552     // Ostream operator
00553 
00554         friend Ostream& operator<<(Ostream&, const lduMatrix&);
00555 };
00556 
00557 
00558 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
00559 
00560 } // End namespace Foam
00561 
00562 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
00563 
00564 #ifdef NoRepository
00565 #   include "lduMatrixHoperations.C"
00566 #endif
00567 
00568 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
00569 
00570 #endif
00571 
00572 // ************************************************************************* //
For further information go to www.openfoam.org