FEATURE: Updates to coupled p-U solver, including ddt term and performance tuning. Author: Hrvoje Jasak. Merge: Dominik Christ.

applications/solvers/coupled/pUCoupledFoam/UEqn.H

@@ -1,7 +1,8 @@
 // Momentum equation
 fvVectorMatrix UEqn
 (
-    fvm::div(phi, U)
+    fvm::ddt(U)
+  + fvm::div(phi, U)
   + turbulence->divDevReff(U)
 );
 

applications/solvers/coupled/pUCoupledFoam/createFields.H

@@ -45,7 +45,7 @@ volVector4Field Up
        runTime.timeName(),
        mesh,
        IOobject::NO_READ,
-       IOobject::AUTO_WRITE
+       IOobject::NO_WRITE
     ),
     mesh,
     dimensionedVector4("zero", dimless, vector4::zero)

src/foam/matrices/blockLduMatrix/BlockAmg/BlockMatrixCoarsening/BlockMatrixAgglomeration.C

@@ -213,8 +213,6 @@ void Foam::BlockMatrixAgglomeration<Type>::calcAgglomeration()
             }
         }
 
-        reduce(nSolo, sumOp<label>());
-
         if (nSolo > 0)
         {
             // Found solo equations
@@ -222,7 +220,7 @@ void Foam::BlockMatrixAgglomeration<Type>::calcAgglomeration()
 
             if (BlockLduMatrix<Type>::debug >= 2)
             {
-                Info<< "Found " << nSolo << " weakly connected equations."
+                Pout<< "Found " << nSolo << " weakly connected equations."
                     << endl;
             }
         }

src/foam/matrices/blockLduMatrix/BlockAmg/coarseBlockAmgLevel.C

@@ -37,7 +37,6 @@ Author
 #include "vector2D.H"
 #include "coeffFields.H"
 #include "BlockSolverPerformance.H"
-// #include "BlockGaussSeidelSolver.H"
 // #include "BlockBiCGStabSolver.H"
 // #include "BlockCGSolver.H"
 #include "BlockGMRESSolver.H"
@@ -78,7 +77,8 @@ Foam::coarseBlockAmgLevel<Type>::coarseBlockAmgLevel
             matrixPtr_,
             dict
         )
-    )
+    ),
+    Ax_()
 {}
 
 
@@ -221,12 +221,23 @@ void Foam::coarseBlockAmgLevel<Type>::solve
         return;
     }
 
+    // Switch of debug in top-level direct solve
+    label oldDebug = BlockLduMatrix<Type>::debug;
+
+    if (BlockLduMatrix<Type>::debug >= 3)
+    {
+        BlockLduMatrix<Type>::debug = 1;
+    }
+    else
+    {
+        BlockLduMatrix<Type>::debug = 0;
+    }
+
     if (matrixPtr_->symmetric())
     {
         topLevelDict.add("preconditioner", "Cholesky");
 
         coarseSolverPerf =
-        //BlockGaussSeidelSolver<Type>
 //         BlockCGSolver<Type>
         BlockGMRESSolver<Type>
         (
@@ -240,7 +251,6 @@ void Foam::coarseBlockAmgLevel<Type>::solve
         topLevelDict.add("preconditioner", "Cholesky");
 
         coarseSolverPerf =
-        //BlockGaussSeidelSolver<Type>
 //         BlockBiCGStabSolver<Type>
         BlockGMRESSolver<Type>
         (
@@ -250,6 +260,9 @@ void Foam::coarseBlockAmgLevel<Type>::solve
         ).solve(x, b);
     }
 
+    // Restore debug
+    BlockLduMatrix<Type>::debug = oldDebug;
+
     // Escape cases of top-level solver divergence
     if
     (
@@ -267,11 +280,6 @@ void Foam::coarseBlockAmgLevel<Type>::solve
         // Print top level correction failure as info for user
         coarseSolverPerf.print();
     }
-
-    if (BlockLduMatrix<Type>::debug >= 2)
-    {
-        coarseSolverPerf.print();
-    }
 }
 
 
@@ -283,18 +291,22 @@ void Foam::coarseBlockAmgLevel<Type>::scaleX
     Field<Type>& xBuffer
 ) const
 {
+    // KRJ: 2013-02-05: Removed subfield, creating a new field
+    // Initialise and size buffer to avoid multiple allocation.
+    // Buffer is created as private data of AMG level
+    // HJ, 26/Feb/2015
+    if (Ax_.empty())
+    {
+        Ax_.setSize(x.size());
+    }
 
-    // KRJ: 2013-02-05: Creating a new field not re-using
+    matrixPtr_->Amul(Ax_, x);
-    Field<Type> Ax(x.size());
 
-    matrixPtr_->Amul
+#if 0
-    (
-        reinterpret_cast<Field<Type>&>(Ax),
-        x
-    );
 
+    // Variant 1 (old): scale complete x with a single scaling factor
     scalar scalingFactorNum = sumProd(x, b);
-    scalar scalingFactorDenom = sumProd(x, Ax);
+    scalar scalingFactorDenom = sumProd(x, Ax_);
 
     vector2D scalingVector(scalingFactorNum, scalingFactorDenom);
     reduce(scalingVector, sumOp<vector2D>());
@@ -318,8 +330,50 @@ void Foam::coarseBlockAmgLevel<Type>::scaleX
     else
     {
         // Regular scaling
-        x *= scalingVector[0]/stabilise(scalingVector[1], SMALL);
+        x *= scalingVector[0]/stabilise(scalingVector[1], VSMALL);
     }
+
+#else
+
+    // Variant 2: scale each x individually
+    // HJ, 25/Feb/2015
+    Type scalingFactorNum = sumCmptProd(x, b);
+    Type scalingFactorDenom = sumCmptProd(x, Ax_);
+
+    Vector2D<Type> scalingVector(scalingFactorNum, scalingFactorDenom);
+    reduce(scalingVector, sumOp<Vector2D<Type> >());
+
+    Type scalingFactor = pTraits<Type>::one;
+
+    // Scale x
+    for (direction dir = 0; dir < pTraits<Type>::nComponents; dir++)
+    {
+        scalar num = component(scalingVector[0], dir);
+        scalar denom = component(scalingVector[1], dir);
+
+        if
+        (
+            mag(num) > GREAT || mag(denom) > GREAT
+         || num*denom <= 0 || mag(num) < mag(denom)
+        )
+        {
+            // Factor = 1.0, no scaling
+        }
+        else if (mag(num) > 2*mag(denom))
+        {
+            setComponent(scalingFactor, dir) = 2;
+        }
+        else
+        {
+            // Regular scaling
+            setComponent(scalingFactor, dir) = num/stabilise(denom, VSMALL);
+        }
+    }
+
+    // Scale X
+    cmptMultiply(x, x, scalingFactor);
+
+#endif
 }
 
 

src/foam/matrices/blockLduMatrix/BlockAmg/coarseBlockAmgLevel.H

@@ -76,6 +76,9 @@ class coarseBlockAmgLevel
         //- Smoother
         autoPtr<BlockLduSmoother<Type> > smootherPtr_;
 
+        //- Ax buffer
+        mutable Field<Type> Ax_;
+
 
     // Private Member Functions
 
@@ -91,6 +94,9 @@ public:
     //- Runtime type information
         TypeName("coarseBlockAmgLevel");
 
+
+    // Constructors
+
         //- Construct from components
         coarseBlockAmgLevel
         (

src/foam/matrices/blockLduMatrix/BlockAmg/fineBlockAmgLevel.C

@@ -74,7 +74,8 @@ Foam::fineBlockAmgLevel<Type>::fineBlockAmgLevel
             matrix,
             dict
         )
-    )
+    ),
+    Ax_()
 {}
 
 
@@ -184,7 +185,6 @@ void Foam::fineBlockAmgLevel<Type>::solve
 
     // Create artificial dictionary for finest solution
     dictionary finestDict;
-//     finestDict.add("nDirections", "5");
     finestDict.add("minIter", 1);
     finestDict.add("maxIter", 1000);
     finestDict.add("tolerance", tolerance);
@@ -201,11 +201,6 @@ void Foam::fineBlockAmgLevel<Type>::solve
                 matrix_,
                 finestDict
             ).solve(x, b);
-
-        if (BlockLduMatrix<Type>::debug >= 2)
-        {
-            coarseSolverPerf.print();
-        }
     }
     else
     {
@@ -218,11 +213,6 @@ void Foam::fineBlockAmgLevel<Type>::solve
                 matrix_,
                 finestDict
             ).solve(x, b);
-
-        if (BlockLduMatrix<Type>::debug >= 2)
-        {
-            coarseSolverPerf.print();
-        }
     }
 }
 
@@ -235,26 +225,31 @@ void Foam::fineBlockAmgLevel<Type>::scaleX
     Field<Type>& xBuffer
 ) const
 {
-
     // KRJ: 2013-02-05: Removed subfield, creating a new field
-    Field<Type> Ax(x.size());
+    // Initialise and size buffer to avoid multiple allocation.
+    // Buffer is created as private data of AMG level
+    // HJ, 26/Feb/2015
+    if (Ax_.empty())
+    {
+        Ax_.setSize(x.size());
+    }
 
-    matrix_.Amul
+    matrix_.Amul(Ax_, x);
-    (
+
-        reinterpret_cast<Field<Type>&>(Ax),
+#if 0
-        x
-    );
 
     scalar scalingFactorNum = sumProd(x, b);
-    scalar scalingFactorDenom = sumProd(x, Ax);
+    scalar scalingFactorDenom = sumProd(x, Ax_);
 
-    vector scalingVector(scalingFactorNum, scalingFactorDenom, 0);
+    vector2D scalingVector(scalingFactorNum, scalingFactorDenom);
-    reduce(scalingVector, sumOp<vector>());
+    reduce(scalingVector, sumOp<vector2D>());
 
     // Scale x
     if
     (
-        scalingVector[0]*scalingVector[1] <= 0
+        mag(scalingVector[0]) > GREAT
+     || mag(scalingVector[1]) > GREAT
+     || scalingVector[0]*scalingVector[1] <= 0
      || mag(scalingVector[0]) < mag(scalingVector[1])
     )
     {
@@ -268,8 +263,48 @@ void Foam::fineBlockAmgLevel<Type>::scaleX
     else
     {
         // Regular scaling
-        x *= scalingVector[0]/stabilise(scalingVector[1], SMALL);
+        x *= scalingVector[0]/stabilise(scalingVector[1], VSMALL);
     }
+
+#else
+
+    Type scalingFactorNum = sumCmptProd(x, b);
+    Type scalingFactorDenom = sumCmptProd(x, Ax_);
+
+    Vector2D<Type> scalingVector(scalingFactorNum, scalingFactorDenom);
+    reduce(scalingVector, sumOp<Vector2D<Type> >());
+
+    Type scalingFactor = pTraits<Type>::one;
+
+    // Scale x
+    for (direction dir = 0; dir < pTraits<Type>::nComponents; dir++)
+    {
+        scalar num = component(scalingVector[0], dir);
+        scalar denom = component(scalingVector[1], dir);
+
+        if
+        (
+            mag(num) > GREAT || mag(denom) > GREAT
+         || num*denom <= 0 || mag(num) < mag(denom)
+        )
+        {
+            // Factor = 1.0, no scaling
+        }
+        else if (mag(num) > 2*mag(denom))
+        {
+            setComponent(scalingFactor, dir) = 2;
+        }
+        else
+        {
+            // Regular scaling
+            setComponent(scalingFactor, dir) = num/stabilise(denom, VSMALL);
+        }
+    }
+
+    // Scale X
+    cmptMultiply(x, x, scalingFactor);
+
+#endif
 }
 
 

src/foam/matrices/blockLduMatrix/BlockAmg/fineBlockAmgLevel.H

@@ -80,6 +80,9 @@ class fineBlockAmgLevel
         //- Smoother
         autoPtr<BlockLduSmoother<Type> > smootherPtr_;
 
+        //- Ax buffer
+        mutable Field<Type> Ax_;
+
 
     // Private Member Functions
 
@@ -96,6 +99,8 @@ public:
     TypeName("fineBlockAmgLevel");
 
 
+    // Constructors
+
         //- Construct from components
         fineBlockAmgLevel
         (

src/foam/matrices/blockLduMatrix/BlockLduSolvers/BlockGMRES/BlockGMRESSolver.C

@@ -113,9 +113,6 @@ Foam::BlockGMRESSolver<Type>::solve
 
     scalar norm = this->normFactor(x, b);
 
-    // Multiplication helper
-    typename BlockCoeff<Type>::multiply mult;
-
     Field<Type> wA(x.size());
 
     // Calculate initial residual