doc/README.md

@@ -14,7 +14,7 @@ make PETSC_DIR=... PETSC_ARCH=... install
 ```
 
 Note that this approach does not account for mixing of mpi
-distributions or versions.
+distributions or versions. Also, it will configure PETSc in debug mode, with an impact on performances.
 
 Export the appropriate PETSC_ARCH_PATH location and/or adjust config
 files(s):
@@ -24,20 +24,35 @@ files(s):
 
 
 
-Another possible configuration
+Another possible configuration for an optimized build of PETSc
 ```
 ./configure
     --with-64-bit-indices=0 \
     --with-precision=double \
+    --with-debugging=0 \
+    --COPTFLAGS=-O3 \
+    --CXXOPTFLAGS=-O3 \
+    --FOPTFLAGS=-O3 \
     --prefix=$WM_THIRD_PARTY_DIR/platforms/$WM_ARCH$WM_COMPILER$WM_PRECISION_OPTION$WM_LABEL_OPTION/petsc-git \
     PETSC_ARCH=$WM_OPTIONS \
-    --download-f2cblaslapack \
+    --download-openblas \
     --with-fc=0 \
     --with-mpi-dir=$MPI_ARCH_PATH
 ```
 
 The ThirdParty [makePETSC][makePETSC] script may be of some use.
 
+For performant runs, we strongly advise you to specify the location of your favorite BLAS/LAPACK libraries either using
+```
+--with-blaslapack-lib=_comma_separated_list_of_libraries_needed_to_resolve_blas/lapack_symbols
+```
+or
+```
+--with-blaslapack-dir=_prefix_location_where_to_find_include_and_libraries
+```
+
+For additional information on the PETSc configure process, see [here](https://www.mcs.anl.gov/petsc/documentation/installation.html).
+
 # Running
 
 Requires LD_LIBRARY_PATH to include PETSc information.

src/petsc4Foam/solvers/petscSolver.C

@@ -157,14 +157,11 @@ Foam::solverPerformance Foam::petscSolver::scalarSolve
     List<label>& lowNonZero = ctx.lowNonZero;
     label& maxLowNonZeroPerRow = ctx.maxLowNonZeroPerRow;
 
-    bool needsCacheUpdate = true;
-
     if (!ctx.initialized())
     {
         DebugInfo<< "Initializing PETSc Linear Solver " << eqName_ << nl;
 
         ctx.initialized() = true;
-        needsCacheUpdate = false;
 
         PetscLogStageRegister
         (
@@ -172,20 +169,30 @@ Foam::solverPerformance Foam::petscSolver::scalarSolve
             &ctx.matstage
         );
         PetscLogStageRegister
+        (
+            ("foam_" + eqName_ + "_pc").c_str(),
+            &ctx.pcstage
+        );
+        PetscLogStageRegister
         (
             ("foam_" + eqName_ + "_ksp").c_str(),
             &ctx.kspstage
         );
-        computeMatAllocation(Amat, lowNonZero, maxLowNonZeroPerRow);
-        buildMat(Amat, lowNonZero, maxLowNonZeroPerRow, ctx.matstage);
+        PetscLogStagePush(ctx.matstage);
+        buildMat(Amat, lowNonZero, maxLowNonZeroPerRow);
+        PetscLogStagePop();
         buildKsp(Amat, ksp);
     }
 
-    if (ctx.caching.needsMatrixUpdate() && needsCacheUpdate)
+    const bool matup = ctx.caching.needsMatrixUpdate();
+    const bool pcup = ctx.caching.needsPrecondUpdate();
+    if (matup)
     {
-        buildMat(Amat, lowNonZero, maxLowNonZeroPerRow, ctx.matstage);
-        updateKsp(ksp, Amat, !ctx.caching.needsPrecondUpdate());
+        PetscLogStagePush(ctx.matstage);
+        updateMat(Amat, lowNonZero, maxLowNonZeroPerRow);
+        PetscLogStagePop();
     }
+    updateKsp(ksp, Amat, !pcup);
 
     // Use built-in residual norm computation
     const bool usePetscResidualNorm
@@ -273,24 +280,24 @@ Foam::solverPerformance Foam::petscSolver::scalarSolve
     // Monitoring (print to stdout) residual reduction in OpenFOAM norm
     if (monitorFoamResidualNorm)
     {
-       KSPMonitorSet
-       (
-           ksp,
-           PetscUtils::foamKSPMonitorFoam,
-           &ctx,
-           NULL
-       );
+        KSPMonitorSet
+        (
+            ksp,
+            PetscUtils::foamKSPMonitorFoam,
+            &ctx,
+            NULL
+        );
     }
 
-    PetscLogStagePush(ctx.kspstage);
-
     // Setup KSP (this is not explicitly needed, but we call it to separate
     // PCSetUp from KSPSolve timings when requesting -log_view from PETSc)
+    PetscLogStagePush(ctx.pcstage);
     KSPSetUp(ksp);
+    PetscLogStagePop();
 
     // Solve A x = b
+    PetscLogStagePush(ctx.kspstage);
     KSPSolve(ksp, petsc_source, petsc_psi);
-
     PetscLogStagePop();
 
     // Set nIterations and final residual
@@ -336,15 +343,15 @@ void Foam::petscSolver::updateKsp
     const bool reuse
 ) const
 {
-    DebugInfo<< "Cache-Hit: updating PETSc Amat " << eqName_ << nl;
-
     if (reuse)
     {
+        DebugInfo<< "Cache-Hit: reuse preconditioner " << eqName_ << nl;
+
         KSPSetReusePreconditioner(ksp, PETSC_TRUE);
     }
     else
     {
-        DebugInfo<< "Cache-Hit: updating PETSc Pmat " << eqName_ << nl;
+        DebugInfo<< "Cache-Hit: rebuild preconditioner " << eqName_ << nl;
 
         KSPSetReusePreconditioner(ksp, PETSC_FALSE);
     }
@@ -393,7 +400,7 @@ void Foam::petscSolver::buildKsp
 }
 
 
-void Foam::petscSolver::computeMatAllocation
+void Foam::petscSolver::buildMat
 (
     Mat& Amat,
     List<label>& lowNonZero,
@@ -505,7 +512,7 @@ void Foam::petscSolver::computeMatAllocation
         nonZeroNeigh_.data()
     );
 
-    // The general AIJ preallocator do not handle SELL
+    // The general AIJ preallocator does not handle SELL
     MatSeqSELLSetPreallocation
     (
         Amat,
@@ -537,16 +544,13 @@ void Foam::petscSolver::computeMatAllocation
 }
 
 
-void Foam::petscSolver::buildMat
+void Foam::petscSolver::updateMat
 (
     Mat& Amat,
     List<label>& lowNonZero,
-    label& maxLowNonZeroPerRow,
-    PetscLogStage stage
+    label& maxLowNonZeroPerRow
 ) const
 {
-    PetscLogStagePush(stage);
-
     const lduAddressing& lduAddr = matrix_.mesh().lduAddr();
     const lduInterfacePtrsList interfaces(matrix_.mesh().interfaces());
 
@@ -727,7 +731,6 @@ void Foam::petscSolver::buildMat
 
     MatAssemblyBegin(Amat, MAT_FINAL_ASSEMBLY);
     MatAssemblyEnd(Amat, MAT_FINAL_ASSEMBLY);
-    PetscLogStagePop();
 }
 
 

src/petsc4Foam/solvers/petscSolver.H

@@ -93,20 +93,19 @@ class petscSolver
 
         //- Compute allocation (proc and off-proc)
         //- for the petsc matrix Amat_
-        void computeMatAllocation
+        void buildMat
         (
             Mat& Amat,
             List<label>& lowNonZero,
             label& maxLowNonZeroPerRow
         ) const;
 
-        //- Build Amat_ matrix by inserting new values
-        void buildMat
+        //- Update Amat_ matrix by inserting new values
+        void updateMat
         (
             Mat& Amat,
             List<label>& lowNonZero,
-            label& maxLowNonZeroPerRow,
-            PetscLogStage stage
+            label& maxLowNonZeroPerRow
         ) const;
 
         //- Build the ksp_ krylov solver context

src/petsc4Foam/utils/petscCacheManager.H

@@ -162,8 +162,11 @@ class petscCacheManager
 {
     // Private Data
 
-        //- The current (relative) iteration
-        label iter;
+        //- The current (relative) iterations for matrix
+        label miter;
+
+        //- The current (relative) iterations for preconditioner
+        label piter;
 
         PetscUtils::Caching matrixCaching;
         PetscUtils::Caching precondCaching;
@@ -176,7 +179,8 @@ public:
         //- Default construct
         petscCacheManager()
         :
-            iter{0}
+            miter{0},
+            piter{0}
         {}
 
 
@@ -195,12 +199,12 @@ public:
 
     bool needsMatrixUpdate()
     {
-        return needsUpdate(matrixCaching, iter);
+        return needsUpdate(matrixCaching, miter);
     }
 
     bool needsPrecondUpdate()
     {
-        return needsUpdate(precondCaching, iter);
+        return needsUpdate(precondCaching, piter);
     }
 
     void eventBegin()
@@ -219,11 +223,11 @@ public:
             // Elapsed timing interval (s)
             precondCaching.timeIter = precondCaching.timer_.elapsed();
 
-            if (iter == 0)
+            if (piter == 0)
             {
                 precondCaching.time0 = precondCaching.timeIter;
             }
-            else if (iter == 1)
+            else if (piter == 1)
             {
                 precondCaching.time1 = precondCaching.timeIter;
             }
@@ -235,7 +239,14 @@ public:
          || precondCaching.updateType_ == PetscUtils::Caching::Adaptive
         )
         {
-            ++iter;
+            ++piter;
+        }
+        if
+        (
+            matrixCaching.updateType_ == PetscUtils::Caching::Periodic
+        )
+        {
+            ++miter;
         }
     }
 
@@ -244,17 +255,19 @@ private:
 
     bool needsUpdate(const PetscUtils::Caching& caching, label& iter) const
     {
+        // Default: Always update
+        bool need = true;
+
         switch (caching.updateType_)
         {
             case PetscUtils::Caching::Never:
             {
-                return false;
+                need = false;
                 break;
             }
 
             case PetscUtils::Caching::Always:
             {
-                return true;
                 break;
             }
 
@@ -267,7 +280,10 @@ private:
                 )
                 {
                     iter = 0;
-                    return true;
+                }
+                else
+                {
+                    need = false;
                 }
                 break;
             }
@@ -277,10 +293,10 @@ private:
                 if (iter > 3) // we need at least three times
                 {
                     const double ratio0 =
-                        precondCaching.time0 / precondCaching.timeIter;
+                        caching.time0 / caching.timeIter;
 
                     const double ratio1 =
-                        precondCaching.time1 / precondCaching.timeIter;
+                        caching.time1 / caching.timeIter;
 
                     const int nsteps =
                         min(1e5, ratio0 * (1. / mag(1. - ratio1 + 1e-6)));
@@ -288,16 +304,17 @@ private:
                     if (iter >= nsteps)
                     {
                         iter = 0;
-                        return true;
+                    }
+                    else
+                    {
+                        need = false;
                     }
                 }
-
                 break;
             }
         }
 
-        // Default: Always update
-        return true;
+        return need;
     }
 };
 

src/petsc4Foam/utils/petscLinearSolverContext.H

@@ -69,8 +69,9 @@ public:
         Mat Amat;
         KSP ksp;
 
-        PetscLogStage matstage;
-        PetscLogStage kspstage;
+        // Auxiliary vectors for the OpenFOAM-PETSc L1-norm
+        Vec ArowsSum;
+        Vec *res_l1_w;
 
         solverPerformance performance;
         petscCacheManager caching;
@@ -78,12 +79,12 @@ public:
         List<PetscInt> lowNonZero;
         label maxLowNonZeroPerRow;
 
-        // Auxiliary vectors for the OpenFOAM-PETSc L1-norm
-        Vec ArowsSum;
-        Vec *res_l1_w;
-
         PetscScalar normFactor;
 
+        PetscLogStage matstage;
+        PetscLogStage pcstage;
+        PetscLogStage kspstage;
+
 
     // Constructors
 
@@ -91,24 +92,21 @@ public:
         petscLinearSolverContext()
         :
             init_(false),
-            init_aux_vectors_(false)
+            init_aux_vectors_(false),
+            Amat(nullptr),
+            ksp(nullptr),
+            ArowsSum(nullptr),
+            res_l1_w(nullptr)
         {}
 
 
     //- Destructor
     virtual ~petscLinearSolverContext()
     {
-        if (init_)
-        {
-            MatDestroy(&Amat);
-            KSPDestroy(&ksp);
-        }
-
-        if (init_aux_vectors_)
-        {
-            VecDestroy(&ArowsSum);
-            VecDestroyVecs(2,&res_l1_w);
-        }
+        MatDestroy(&Amat);
+        KSPDestroy(&ksp);
+        VecDestroy(&ArowsSum);
+        if (res_l1_w) VecDestroyVecs(2,&res_l1_w);
     }
 
 
@@ -131,9 +129,9 @@ public:
         {
             if (!init_aux_vectors_)
             {
+                init_aux_vectors_ = true;
                 VecDuplicate(psi, &ArowsSum);
                 VecDuplicateVecs(psi, 2, &res_l1_w);
-                init_aux_vectors_ = true;
             }
         }
 
@@ -151,6 +149,7 @@ public:
         void computeNormFactor(Mat A, Vec psi, Vec source)
         {
             if (!init_aux_vectors_) return;
+
             MatMult(A, psi, res_l1_w[1]);
 
             normFactor =