LduMatrix: Added support for both component-coupled and component-independent forms of PCG and PBiCG

applications/test/PisoFoam/PisoFoam.C

@@ -102,7 +102,7 @@ int main(int argc, char *argv[])
                         "{"
                         "    /*solver          SmoothSolver;*/"
                         "    smoother        GaussSeidel;"
-                        "    solver           PBiCG;"
+                        "    solver           PBiCCCG;"
                         "    preconditioner   DILU;"
                         "    tolerance        (1e-7 1e-7 1);"
                         "    relTol           (0 0 0);"

src/OpenFOAM/fields/Fields/Field/FieldFunctions.C

@@ -365,7 +365,7 @@ scalar sumProd(const UList<Type>& f1, const UList<Type>& f2)
     if (f1.size() && (f1.size() == f2.size()))
     {
         scalar SumProd = 0.0;
-        TFOR_ALL_S_OP_F_OP_F(scalar, SumProd, +=, Type, f1, *, Type, f2)
+        TFOR_ALL_S_OP_F_OP_F(scalar, SumProd, +=, Type, f1, &&, Type, f2)
         return SumProd;
     }
     else
@@ -498,7 +498,7 @@ template<class Type>
 scalar gSumProd(const UList<Type>& f1, const UList<Type>& f2)
 {
     scalar SumProd = sumProd(f1, f2);
-    reduce(SumProd, sumOp<Type>());
+    reduce(SumProd, sumOp<scalar>());
     return SumProd;
 }
 

src/OpenFOAM/fields/Fields/Field/FieldFunctions.H

@@ -23,8 +23,6 @@ License
 
 \*---------------------------------------------------------------------------*/
 
-#include "scalarField.H"
-
 #define TEMPLATE template<class Type>
 #include "FieldFunctionsM.H"
 
@@ -332,5 +330,6 @@ PRODUCT_OPERATOR(scalarProduct, &&, dotdot)
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
 #include "undefFieldFunctionsM.H"
+#include "scalarField.H"
 
 // ************************************************************************* //

src/OpenFOAM/fields/Fields/scalarField/scalarField.C

@@ -86,6 +86,24 @@ tmp<scalarField> stabilise(const tmp<scalarField>& tsf, const scalar s)
 }
 
 
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+template<>
+scalar sumProd(const UList<scalar>& f1, const UList<scalar>& f2)
+{
+    if (f1.size() && (f1.size() == f2.size()))
+    {
+        scalar SumProd = 0.0;
+        TFOR_ALL_S_OP_F_OP_F(scalar, SumProd, +=, scalar, f1, *, scalar, f2)
+        return SumProd;
+    }
+    else
+    {
+        return 0.0;
+    }
+}
+
+
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
 BINARY_TYPE_OPERATOR(scalar, scalar, scalar, +, add)

src/OpenFOAM/fields/Fields/scalarField/scalarField.H

@@ -72,6 +72,12 @@ tmp<scalarField> stabilise(const UList<scalar>&, const scalar s);
 tmp<scalarField> stabilise(const tmp<scalarField>&, const scalar s);
 
 
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+template<>
+scalar sumProd(const UList<scalar>& f1, const UList<scalar>& f2);
+
+
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
 BINARY_TYPE_OPERATOR(scalar, scalar, scalar, +, add)

src/OpenFOAM/matrices/LduMatrix/Solvers/PBiCCCG/PBiCCCG.C

@@ -23,12 +23,12 @@ License
 
 \*---------------------------------------------------------------------------*/
 
-#include "TPBiCG.H"
+#include "PBiCCCG.H"
 
 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 
 template<class Type, class DType, class LUType>
-Foam::TPBiCG<Type, DType, LUType>::TPBiCG
+Foam::PBiCCCG<Type, DType, LUType>::PBiCCCG
 (
     const word& fieldName,
     const LduMatrix<Type, DType, LUType>& matrix,
@@ -48,7 +48,7 @@ Foam::TPBiCG<Type, DType, LUType>::TPBiCG
 
 template<class Type, class DType, class LUType>
 typename Foam::LduMatrix<Type, DType, LUType>::solverPerformance
-Foam::TPBiCG<Type, DType, LUType>::solve
+Foam::PBiCCCG<Type, DType, LUType>::solve
 (
     Field<Type>& psi
 ) const
@@ -125,8 +125,7 @@ Foam::TPBiCG<Type, DType, LUType>::solve
             preconPtr->preconditionT(wT, rT);
 
             // --- Update search directions:
-            //wArT = gSumProd(wA, rT);
+            wArT = gSumProd(wA, rT);
-            wArT = gSum(wA && rT);
 
             if (solverPerf.nIterations() == 0)
             {
@@ -152,14 +151,14 @@ Foam::TPBiCG<Type, DType, LUType>::solve
             this->matrix_.Amul(wA, pA);
             this->matrix_.Tmul(wT, pT);
 
-            scalar wApT = gSum(wA && pT);
+            scalar wApT = gSumProd(wA, pT);
 
             // --- Test for singularity
             if
             (
                 solverPerf.checkSingularity
                 (
-                    cmptDivide(Type::one*mag(wApT), normFactor)
+                    cmptDivide(pTraits<Type>::one*mag(wApT), normFactor)
                 )
             )
             {

src/OpenFOAM/matrices/LduMatrix/Solvers/PBiCCCG/PBiCCCG.H

@@ -22,19 +22,19 @@ License
     along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
 
 Class
-    Foam::TPBiCG
+    Foam::PBiCCCG
 
 Description
     Preconditioned bi-conjugate gradient solver for asymmetric lduMatrices
     using a run-time selectable preconditiioner.
 
 SourceFiles
-    TPBiCG.C
+    PBiCCCG.C
 
 \*---------------------------------------------------------------------------*/
 
-#ifndef TPBiCG_H
+#ifndef PBiCCCG_H
-#define TPBiCG_H
+#define PBiCCCG_H
 
 #include "LduMatrix.H"
 
@@ -44,33 +44,33 @@ namespace Foam
 {
 
 /*---------------------------------------------------------------------------*\
-                           Class TPBiCG Declaration
+                           Class PBiCCCG Declaration
 \*---------------------------------------------------------------------------*/
 
 template<class Type, class DType, class LUType>
-class TPBiCG
+class PBiCCCG
 :
     public LduMatrix<Type, DType, LUType>::solver
 {
     // Private Member Functions
 
         //- Disallow default bitwise copy construct
-        TPBiCG(const TPBiCG&);
+        PBiCCCG(const PBiCCCG&);
 
         //- Disallow default bitwise assignment
-        void operator=(const TPBiCG&);
+        void operator=(const PBiCCCG&);
 
 
 public:
 
     //- Runtime type information
-    TypeName("PBiCG");
+    TypeName("PBiCCCG");
 
 
     // Constructors
 
         //- Construct from matrix components and solver data dictionary
-        TPBiCG
+        PBiCCCG
         (
             const word& fieldName,
             const LduMatrix<Type, DType, LUType>& matrix,
@@ -80,7 +80,7 @@ public:
 
     // Destructor
 
-        virtual ~TPBiCG()
+        virtual ~PBiCCCG()
         {}
 
 
@@ -101,7 +101,7 @@ public:
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
 #ifdef NoRepository
-#   include "TPBiCG.C"
+#   include "PBiCCCG.C"
 #endif
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

src/OpenFOAM/matrices/LduMatrix/Solvers/PBiCICG/PBiCICG.C

@@ -23,12 +23,12 @@ License
 
 \*---------------------------------------------------------------------------*/
 
-#include "TPBiCG.H"
+#include "PBiCICG.H"
 
 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 
 template<class Type, class DType, class LUType>
-Foam::TPBiCG<Type, DType, LUType>::TPBiCG
+Foam::PBiCICG<Type, DType, LUType>::PBiCICG
 (
     const word& fieldName,
     const LduMatrix<Type, DType, LUType>& matrix,
@@ -48,7 +48,7 @@ Foam::TPBiCG<Type, DType, LUType>::TPBiCG
 
 template<class Type, class DType, class LUType>
 typename Foam::LduMatrix<Type, DType, LUType>::solverPerformance
-Foam::TPBiCG<Type, DType, LUType>::solve(Field<Type>& psi) const
+Foam::PBiCICG<Type, DType, LUType>::solve(Field<Type>& psi) const
 {
     word preconditionerName(this->controlDict_.lookup("preconditioner"));
 

src/OpenFOAM/matrices/LduMatrix/Solvers/PBiCICG/PBiCICG.H

@@ -22,19 +22,19 @@ License
     along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
 
 Class
-    Foam::TPBiCG
+    Foam::PBiCICG
 
 Description
     Preconditioned bi-conjugate gradient solver for asymmetric lduMatrices
     using a run-time selectable preconditiioner.
 
 SourceFiles
-    TPBiCG.C
+    PBiCICG.C
 
 \*---------------------------------------------------------------------------*/
 
-#ifndef TPBiCG_H
+#ifndef PBiCICG_H
-#define TPBiCG_H
+#define PBiCICG_H
 
 #include "LduMatrix.H"
 
@@ -44,33 +44,33 @@ namespace Foam
 {
 
 /*---------------------------------------------------------------------------*\
-                           Class TPBiCG Declaration
+                           Class PBiCICG Declaration
 \*---------------------------------------------------------------------------*/
 
 template<class Type, class DType, class LUType>
-class TPBiCG
+class PBiCICG
 :
     public LduMatrix<Type, DType, LUType>::solver
 {
     // Private Member Functions
 
         //- Disallow default bitwise copy construct
-        TPBiCG(const TPBiCG&);
+        PBiCICG(const PBiCICG&);
 
         //- Disallow default bitwise assignment
-        void operator=(const TPBiCG&);
+        void operator=(const PBiCICG&);
 
 
 public:
 
     //- Runtime type information
-    TypeName("PBiCG");
+    TypeName("PBiCICG");
 
 
     // Constructors
 
         //- Construct from matrix components and solver data dictionary
-        TPBiCG
+        PBiCICG
         (
             const word& fieldName,
             const LduMatrix<Type, DType, LUType>& matrix,
@@ -80,7 +80,7 @@ public:
 
     // Destructor
 
-        virtual ~TPBiCG()
+        virtual ~PBiCICG()
         {}
 
 
@@ -101,7 +101,7 @@ public:
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
 #ifdef NoRepository
-#   include "PBiCGScalarAlpha.C"
+#   include "PBiCICG.C"
 #endif
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

src/OpenFOAM/matrices/LduMatrix/Solvers/PCICG/PCICG.C

@@ -23,12 +23,12 @@ License
 
 \*---------------------------------------------------------------------------*/
 
-#include "TPCG.H"
+#include "PCICG.H"
 
 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 
 template<class Type, class DType, class LUType>
-Foam::TPCG<Type, DType, LUType>::TPCG
+Foam::PCICG<Type, DType, LUType>::PCICG
 (
     const word& fieldName,
     const LduMatrix<Type, DType, LUType>& matrix,
@@ -48,7 +48,7 @@ Foam::TPCG<Type, DType, LUType>::TPCG
 
 template<class Type, class DType, class LUType>
 typename Foam::LduMatrix<Type, DType, LUType>::solverPerformance
-Foam::TPCG<Type, DType, LUType>::solve(Field<Type>& psi) const
+Foam::PCICG<Type, DType, LUType>::solve(Field<Type>& psi) const
 {
     word preconditionerName(this->controlDict_.lookup("preconditioner"));
 

src/OpenFOAM/matrices/LduMatrix/Solvers/PCICG/PCICG.H

@@ -22,19 +22,19 @@ License
     along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
 
 Class
-    Foam::TPCG
+    Foam::PCICG
 
 Description
     Preconditioned conjugate gradient solver for symmetric lduMatrices
     using a run-time selectable preconditiioner.
 
 SourceFiles
-    TPCG.C
+    PCICG.C
 
 \*---------------------------------------------------------------------------*/
 
-#ifndef TPCG_H
+#ifndef PCICG_H
-#define TPCG_H
+#define PCICG_H
 
 #include "LduMatrix.H"
 
@@ -44,33 +44,33 @@ namespace Foam
 {
 
 /*---------------------------------------------------------------------------*\
-                           Class TPCG Declaration
+                           Class PCICG Declaration
 \*---------------------------------------------------------------------------*/
 
 template<class Type, class DType, class LUType>
-class TPCG
+class PCICG
 :
     public LduMatrix<Type, DType, LUType>::solver
 {
     // Private Member Functions
 
         //- Disallow default bitwise copy construct
-        TPCG(const TPCG&);
+        PCICG(const PCICG&);
 
         //- Disallow default bitwise assignment
-        void operator=(const TPCG&);
+        void operator=(const PCICG&);
 
 
 public:
 
     //- Runtime type information
-    TypeName("PCG");
+    TypeName("PCICG");
 
 
     // Constructors
 
         //- Construct from matrix components and solver data dictionary
-        TPCG
+        PCICG
         (
             const word& fieldName,
             const LduMatrix<Type, DType, LUType>& matrix,
@@ -80,7 +80,7 @@ public:
 
     // Destructor
 
-        virtual ~TPCG()
+        virtual ~PCICG()
         {}
 
 
@@ -101,7 +101,7 @@ public:
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
 #ifdef NoRepository
-#   include "TPCG.C"
+#   include "PCICG.C"
 #endif
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

src/OpenFOAM/matrices/LduMatrix/Solvers/lduSolvers.C

@@ -23,9 +23,9 @@ License
 
 \*---------------------------------------------------------------------------*/
 
-#include "TPCG.H"
+#include "PCICG.H"
-//#include "TPBiCG.H"
+#include "PBiCCCG.H"
-#include "PBiCGScalarAlpha.H"
+#include "PBiCICG.H"
 #include "SmoothSolver.H"
 #include "fieldTypes.H"
 
@@ -35,11 +35,14 @@ License
     makeLduSymSolver(DiagonalSolver, Type, DType, LUType);                    \
     makeLduAsymSolver(DiagonalSolver, Type, DType, LUType);                   \
                                                                               \
-    makeLduSolver(TPCG, Type, DType, LUType);                                 \
+    makeLduSolver(PCICG, Type, DType, LUType);                                \
-    makeLduSymSolver(TPCG, Type, DType, LUType);                              \
+    makeLduSymSolver(PCICG, Type, DType, LUType);                             \
                                                                               \
-    makeLduSolver(TPBiCG, Type, DType, LUType);                               \
+    makeLduSolver(PBiCCCG, Type, DType, LUType);                              \
-    makeLduAsymSolver(TPBiCG, Type, DType, LUType);                           \
+    makeLduAsymSolver(PBiCCCG, Type, DType, LUType);                          \
+                                                                              \
+    makeLduSolver(PBiCICG, Type, DType, LUType);                              \
+    makeLduAsymSolver(PBiCICG, Type, DType, LUType);                          \
                                                                               \
     makeLduSolver(SmoothSolver, Type, DType, LUType);                         \
     makeLduSymSolver(SmoothSolver, Type, DType, LUType);                      \
@@ -47,7 +50,7 @@ License
 
 namespace Foam
 {
-    //makeLduSolvers(scalar, scalar, scalar);
+    makeLduSolvers(scalar, scalar, scalar);
     makeLduSolvers(vector, scalar, scalar);
     makeLduSolvers(sphericalTensor, scalar, scalar);
     makeLduSolvers(symmTensor, scalar, scalar);