ENH: PCG: cache preconditioner

src/OpenFOAM/matrices/lduMatrix/lduMatrix/lduMatrix.H

@@ -575,6 +575,10 @@ public:
             {
                 NotImplemented;
             }
+
+            //- Signal end of solver
+            virtual void setFinished(const solverPerformance& perf) const
+            {}
     };
 
 

src/OpenFOAM/matrices/lduMatrix/solvers/FPCG/FPCG.C

@@ -159,12 +159,14 @@ Foam::solverPerformance Foam::FPCG::scalarSolve
     )
     {
         // --- Select and construct the preconditioner
-        autoPtr<lduMatrix::preconditioner> preconPtr =
+        if (!preconPtr_)
-            lduMatrix::preconditioner::New
+        {
+            preconPtr_ = lduMatrix::preconditioner::New
             (
                 *this,
                 controlDict_
             );
+        }
 
         FixedList<solveScalar, 2> globalSum;
 
@@ -175,7 +177,7 @@ Foam::solverPerformance Foam::FPCG::scalarSolve
             wArAold = wArA;
 
             // --- Precondition residual
-            preconPtr->precondition(wA, rA, cmpt);
+            preconPtr_->precondition(wA, rA, cmpt);
 
             // --- Update search directions and calculate residual:
             gSumMagProd(globalSum, wA, rA, matrix().mesh().comm());
@@ -237,6 +239,11 @@ Foam::solverPerformance Foam::FPCG::scalarSolve
         );
     }
 
+    if (preconPtr_)
+    {
+        preconPtr_->setFinished(solverPerf);
+    }
+
     matrix().setResidualField
     (
         ConstPrecisionAdaptor<scalar, solveScalar>(rA)(),

src/OpenFOAM/matrices/lduMatrix/solvers/FPCG/FPCG.H

@@ -61,6 +61,12 @@ class FPCG
 :
     public lduMatrix::solver
 {
+    // Private Member Data
+
+        //- Cached preconditioner
+        mutable autoPtr<lduMatrix::preconditioner> preconPtr_;
+
+
     // Private Member Functions
 
         //- Blocking version of sum(a*b), sum(mag(b))

src/OpenFOAM/matrices/lduMatrix/solvers/PBiCG/PBiCG.C

@@ -146,12 +146,14 @@ Foam::solverPerformance Foam::PBiCG::solve
         solveScalar wArT = 0;
 
         // --- Select and construct the preconditioner
-        autoPtr<lduMatrix::preconditioner> preconPtr =
+        if (!preconPtr_)
-        lduMatrix::preconditioner::New
+        {
+            preconPtr_ = lduMatrix::preconditioner::New
             (
                 *this,
                 controlDict_
             );
+        }
 
         // --- Solver iteration
         do
@@ -160,8 +162,8 @@ Foam::solverPerformance Foam::PBiCG::solve
             const solveScalar wArTold = wArT;
 
             // --- Precondition residuals
-            preconPtr->precondition(wA, rA, cmpt);
+            preconPtr_->precondition(wA, rA, cmpt);
-            preconPtr->preconditionT(wT, rT, cmpt);
+            preconPtr_->preconditionT(wT, rT, cmpt);
 
             // --- Update search directions:
             wArT = gSumProd(wA, rT, matrix().mesh().comm());
@@ -235,6 +237,11 @@ Foam::solverPerformance Foam::PBiCG::solve
             << exit(FatalError);
     }
 
+    if (preconPtr_)
+    {
+        preconPtr_->setFinished(solverPerf);
+    }
+
     matrix().setResidualField
     (
         ConstPrecisionAdaptor<scalar, solveScalar>(rA)(),

src/OpenFOAM/matrices/lduMatrix/solvers/PBiCG/PBiCG.H

@@ -54,6 +54,12 @@ class PBiCG
 :
     public lduMatrix::solver
 {
+    // Private Member Data
+
+        //- Cached preconditioner
+        mutable autoPtr<lduMatrix::preconditioner> preconPtr_;
+
+
     // Private Member Functions
 
         //- No copy construct

src/OpenFOAM/matrices/lduMatrix/solvers/PBiCGStab/PBiCGStab.C

@@ -149,12 +149,14 @@ Foam::solverPerformance Foam::PBiCGStab::scalarSolve
         solveScalar omega = 0;
 
         // --- Select and construct the preconditioner
-        autoPtr<lduMatrix::preconditioner> preconPtr =
+        if (!preconPtr_)
-        lduMatrix::preconditioner::New
+        {
+            preconPtr_ = lduMatrix::preconditioner::New
             (
                 *this,
                 controlDict_
             );
+        }
 
         // --- Solver iteration
         do
@@ -196,7 +198,7 @@ Foam::solverPerformance Foam::PBiCGStab::scalarSolve
             }
 
             // --- Precondition pA
-            preconPtr->precondition(yA, pA, cmpt);
+            preconPtr_->precondition(yA, pA, cmpt);
 
             // --- Calculate AyA
             matrix_.Amul(AyA, yA, interfaceBouCoeffs_, interfaces_, cmpt);
@@ -233,7 +235,7 @@ Foam::solverPerformance Foam::PBiCGStab::scalarSolve
             }
 
             // --- Precondition sA
-            preconPtr->precondition(zA, sA, cmpt);
+            preconPtr_->precondition(zA, sA, cmpt);
 
             // --- Calculate tA
             matrix_.Amul(tA, zA, interfaceBouCoeffs_, interfaces_, cmpt);
@@ -264,6 +266,11 @@ Foam::solverPerformance Foam::PBiCGStab::scalarSolve
         );
     }
 
+    if (preconPtr_)
+    {
+        preconPtr_->setFinished(solverPerf);
+    }
+
     matrix().setResidualField
     (
         ConstPrecisionAdaptor<scalar, solveScalar>(rA)(),

src/OpenFOAM/matrices/lduMatrix/solvers/PBiCGStab/PBiCGStab.H

@@ -69,6 +69,12 @@ class PBiCGStab
 :
     public lduMatrix::solver
 {
+    // Private Member Data
+
+        //- Cached preconditioner
+        mutable autoPtr<lduMatrix::preconditioner> preconPtr_;
+
+
     // Private Member Functions
 
         //- No copy construct

src/OpenFOAM/matrices/lduMatrix/solvers/PCG/PCG.C

@@ -6,7 +6,7 @@
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
     Copyright (C) 2011-2017 OpenFOAM Foundation
-    Copyright (C) 2019-2021 OpenCFD Ltd.
+    Copyright (C) 2019-2023 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -129,12 +129,14 @@ Foam::solverPerformance Foam::PCG::scalarSolve
     )
     {
         // --- Select and construct the preconditioner
-        autoPtr<lduMatrix::preconditioner> preconPtr =
+        if (!preconPtr_)
-            lduMatrix::preconditioner::New
+        {
+            preconPtr_ = lduMatrix::preconditioner::New
             (
                 *this,
                 controlDict_
             );
+        }
 
         // --- Solver iteration
         do
@@ -143,7 +145,7 @@ Foam::solverPerformance Foam::PCG::scalarSolve
             wArAold = wArA;
 
             // --- Precondition residual
-            preconPtr->precondition(wA, rA, cmpt);
+            preconPtr_->precondition(wA, rA, cmpt);
 
             // --- Update search directions:
             wArA = gSumProd(wA, rA, matrix().mesh().comm());
@@ -199,6 +201,11 @@ Foam::solverPerformance Foam::PCG::scalarSolve
         );
     }
 
+    if (preconPtr_)
+    {
+        preconPtr_->setFinished(solverPerf);
+    }
+
     matrix().setResidualField
     (
         ConstPrecisionAdaptor<scalar, solveScalar>(rA)(),

src/OpenFOAM/matrices/lduMatrix/solvers/PCG/PCG.H

@@ -57,6 +57,12 @@ class PCG
 :
     public lduMatrix::solver
 {
+    // Private Member Data
+
+        //- Cached preconditioner
+        mutable autoPtr<lduMatrix::preconditioner> preconPtr_;
+
+
     // Private Member Functions
 
         //- No copy construct

src/OpenFOAM/matrices/lduMatrix/solvers/PPCG/PPCG.C

@@ -113,16 +113,18 @@ Foam::solverPerformance Foam::PPCG::scalarSolveCG
     }
 
     // --- Select and construct the preconditioner
-    autoPtr<lduMatrix::preconditioner> preconPtr =
+    if (!preconPtr_)
-        lduMatrix::preconditioner::New
+    {
+        preconPtr_ = lduMatrix::preconditioner::New
         (
             *this,
             controlDict_
         );
+    }
 
     // --- Precondition residual (= u0)
     solveScalarField u(nCells);
-    preconPtr->precondition(u, r, cmpt);
+    preconPtr_->precondition(u, r, cmpt);
 
     // --- Calculate A*u - reuse w
     matrix_.Amul(w, u, interfaceBouCoeffs_, interfaces_, cmpt);
@@ -143,12 +145,12 @@ Foam::solverPerformance Foam::PPCG::scalarSolveCG
         gSumMagProd(globalSum, u, r, w, r, outstandingRequest, comm);
 
         // --- Precondition residual
-        preconPtr->precondition(m, w, cmpt);
+        preconPtr_->precondition(m, w, cmpt);
     }
     else
     {
         // --- Precondition residual
-        preconPtr->precondition(m, w, cmpt);
+        preconPtr_->precondition(m, w, cmpt);
 
         // --- Start global reductions for inner products
         gSumMagProd(globalSum, w, u, m, r, outstandingRequest, comm);
@@ -229,12 +231,12 @@ Foam::solverPerformance Foam::PPCG::scalarSolveCG
             gSumMagProd(globalSum, u, r, w, r, outstandingRequest, comm);
 
             // --- Precondition residual
-            preconPtr->precondition(m, w, cmpt);
+            preconPtr_->precondition(m, w, cmpt);
         }
         else
         {
             // --- Precondition residual
-            preconPtr->precondition(m, w, cmpt);
+            preconPtr_->precondition(m, w, cmpt);
 
             // --- Start global reductions for inner products
             gSumMagProd(globalSum, w, u, m, r, outstandingRequest, comm);
@@ -247,11 +249,23 @@ Foam::solverPerformance Foam::PPCG::scalarSolveCG
     // Cleanup any outstanding requests
     outstandingRequest.wait();
 
+    if (preconPtr_)
+    {
+        preconPtr_->setFinished(solverPerf);
+    }
+
+    //TBD
+    //matrix().setResidualField
+    //(
+    //    ConstPrecisionAdaptor<scalar, solveScalar>(rA)(),
+    //    fieldName_,
+    //    false
+    //);
+
     return solverPerf;
 }
 
 
-
 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 
 Foam::PPCG::PPCG

src/OpenFOAM/matrices/lduMatrix/solvers/PPCG/PPCG.H

@@ -66,6 +66,12 @@ class PPCG
 :
     public lduMatrix::solver
 {
+    // Private Member Data
+
+        //- Cached preconditioner
+        mutable autoPtr<lduMatrix::preconditioner> preconPtr_;
+
+
     // Private Member Functions
 
         //- Non-blocking version of sum(a*b), sum(a*c), sum(mag(sumMag))

src/meshTools/matrices/lduMatrix/preconditioners/distributedDILUPreconditioner/distributedDILUPreconditioner.C

@@ -248,10 +248,18 @@ void Foam::distributedDILUPreconditioner::send
 
 void Foam::distributedDILUPreconditioner::wait
 (
-    DynamicList<UPstream::Request>& requests
+    DynamicList<UPstream::Request>& requests,
+    const bool cancel
 ) const
 {
+    if (cancel)
+    {
+        UPstream::cancelRequests(requests);
+    }
+    else
+    {
         UPstream::waitRequests(requests);
+    }
     requests.clear();
 }
 
@@ -757,12 +765,20 @@ Foam::distributedDILUPreconditioner::distributedDILUPreconditioner
 
 Foam::distributedDILUPreconditioner::~distributedDILUPreconditioner()
 {
+    DebugPout<< "~distributedDILUPreconditioner()" << endl;
+
     // Wait on all requests before storage is being taken down
-    // (could rely on construction order?)
+
     wait(lowerSendRequests_);
     wait(lowerRecvRequests_);
     wait(higherSendRequests_);
     wait(higherRecvRequests_);
+
+    // TBD: cancel/ignore outstanding messages
+    //wait(lowerSendRequests_, true);
+    //wait(lowerRecvRequests_, true);
+    //wait(higherSendRequests_, true);
+    //wait(higherRecvRequests_, true);
 }
 
 
@@ -894,4 +910,26 @@ void Foam::distributedDILUPreconditioner::precondition
 }
 
 
+void Foam::distributedDILUPreconditioner::setFinished
+(
+    const solverPerformance& s
+) const
+{
+    DebugPout<< "setFinished fieldName:" << s.fieldName() << endl;
+
+    // Wait on all requests before storage is being taken down
+    // (could rely on construction order?)
+    wait(lowerSendRequests_);
+    wait(lowerRecvRequests_);
+    wait(higherSendRequests_);
+    wait(higherRecvRequests_);
+
+    // TBD: cancel/ignore outstanding messages
+    //wait(lowerSendRequests_, true);
+    //wait(lowerRecvRequests_, true);
+    //wait(higherSendRequests_, true);
+    //wait(higherRecvRequests_, true);
+}
+
+
 // ************************************************************************* //

src/meshTools/matrices/lduMatrix/preconditioners/distributedDILUPreconditioner/distributedDILUPreconditioner.H

@@ -62,8 +62,8 @@ SourceFiles
 
 \*---------------------------------------------------------------------------*/
 
-#ifndef distributedDILUPreconditioner_H
+#ifndef Foam_distributedDILUPreconditioner_H
-#define distributedDILUPreconditioner_H
+#define Foam_distributedDILUPreconditioner_H
 
 #include "lduMatrix.H"
 
@@ -82,7 +82,7 @@ class distributedDILUPreconditioner
 {
 protected:
 
-    // Protected data
+    // Protected Data
 
         //- Precondition across global coupled bc
         const bool coupled_;
@@ -187,8 +187,12 @@ protected:
             DynamicList<UPstream::Request>& requests
         ) const;
 
-        //- Wait for explicit requests
+        //- Wait for requests or cancel/free requests
-        void wait(DynamicList<UPstream::Request>& requests) const;
+        void wait
+        (
+            DynamicList<UPstream::Request>& requests,
+            const bool cancel = false
+        ) const;
 
         //- Update diagonal for interface
         virtual void addInterfaceDiag
@@ -260,6 +264,9 @@ public:
             const solveScalarField& rA,
             const direction cmpt=0
         ) const;
+
+        //- Signal end of solver
+        virtual void setFinished(const solverPerformance& perf) const;
 };