src/petsc4Foam/Make/options

@@ -12,11 +12,15 @@ EXE_INC = \
     $(PFLAGS) $(PINC) \
     -Wno-old-style-cast \
     -I$(LIB_SRC)/finiteVolume/lnInclude \
+    -I$(LIB_SRC)/fileFormats/lnInclude \
+    -I$(LIB_SRC)/surfMesh/lnInclude \
     -I$(LIB_SRC)/meshTools/lnInclude \
     -I$(LIB_SRC)/Pstream/mpi \
     $(foreach dir,$(PETSC_INC_DIR),-I$(dir))
 
 LIB_LIBS = \
     -lfiniteVolume \
+    -lfileFormats \
+    -lsurfMesh \
     -lmeshTools \
     $(foreach dir,$(PETSC_LIB_DIR),-L$(dir)) -lpetsc

src/petsc4Foam/solvers/petscSolver.C

@@ -5,7 +5,7 @@
     \\  /    A nd           | www.openfoam.com
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
-    Copyright (C) 2018-2020 OpenCFD Ltd.
+    Copyright (C) 2018-2022 OpenCFD Ltd.
     Copyright (C) 2019-2020 Simone Bna
     Copyright (C) 2021 Stefano Zampini
 -------------------------------------------------------------------------------
@@ -33,15 +33,17 @@ License
 #include "PrecisionAdaptor.H"
 #include "cyclicLduInterface.H"
 #include "cyclicAMILduInterface.H"
-#include "addToRunTimeSelectionTable.H"
 #include "PstreamGlobals.H"
+#include "addToRunTimeSelectionTable.H"
 
 #include "petscSolver.H"
 #include "petscControls.H"
 #include "petscLinearSolverContexts.H"
 #include "petscUtils.H"
 #include "petscWrappedVector.H"
-#include <cstring>
+
+#include <algorithm>  // For std::max_element
+#include <cstring>    // For NULL
 
 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
 
@@ -156,14 +158,14 @@ Foam::solverPerformance Foam::petscSolver::scalarSolve
     Mat& Amat = ctx.Amat;
     KSP& ksp  = ctx.ksp;
 
-    List<label>& lowNonZero = ctx.lowNonZero;
-    label& maxLowNonZeroPerRow = ctx.maxLowNonZeroPerRow;
+    List<PetscInt>& lowNonZero = ctx.lowNonZero;
+    PetscInt& maxLowNonZeroPerRow = ctx.maxLowNonZeroPerRow;
 
     if (firsttimein)
     {
         DebugInfo<< "Initializing PETSc Linear Solver " << eqName_ << nl;
 
-        ctx.initialized() = true;
+        ctx.initialized(true);
 
         PetscLogStageRegister
         (
@@ -190,6 +192,7 @@ Foam::solverPerformance Foam::petscSolver::scalarSolve
     const bool pcup = ctx.caching.needsPrecondUpdate();
     if (matup)
     {
+        DebugInfo<< "Update PETSc Linear Solver Matrix for: " << eqName_ << nl;
         PetscLogStagePush(ctx.matstage);
         updateMat(Amat, lowNonZero, maxLowNonZeroPerRow);
         PetscLogStagePop();
@@ -255,10 +258,11 @@ Foam::solverPerformance Foam::petscSolver::scalarSolve
     {
         ctx.createAuxVecs();
         ctx.initArowsSumVec();
-        ctx.computeNormFactor(petsc_psi,petsc_source);
+        ctx.computeNormFactor(petsc_psi, petsc_source);
         if (normType == KSP_NORM_NONE && !usePetscResidualNorm)
         {
-           ctx.performance.initialResidual() = ctx.getResidualNormL1(petsc_source);
+            ctx.performance.initialResidual() =
+                ctx.getResidualNormL1(petsc_source);
         }
     }
 
@@ -360,13 +364,16 @@ Foam::solverPerformance Foam::petscSolver::scalarSolve
         PetscReal rnorm;
         KSPNormType normType;
         KSPGetNormType(ksp, &normType);
-        if (normType == KSP_NORM_NONE) // We run with PETSc norm KSP_NORM_NONE -> report L1 norm
+        if (normType == KSP_NORM_NONE)
         {
-           rnorm = ctx.getResidualNormL1(petsc_psi, petsc_source);
+            // We run with PETSc norm KSP_NORM_NONE -> report L1 norm
+            rnorm = ctx.getResidualNormL1(petsc_psi, petsc_source);
         }
-        else // report final PETSc norm since users explicitly ask to run with PETSc norm
+        else
         {
-           KSPGetResidualNorm(ksp, &rnorm);
+            // Report final PETSc norm since users explicitly ask to run
+            // with PETSc norm
+            KSPGetResidualNorm(ksp, &rnorm);
         }
         ctx.performance.finalResidual() = rnorm;
     }
@@ -460,8 +467,8 @@ void Foam::petscSolver::buildKsp
 void Foam::petscSolver::buildMat
 (
     Mat& Amat,
-    List<label>& lowNonZero,
-    label& maxLowNonZeroPerRow
+    List<PetscInt>& lowNonZero,
+    PetscInt& maxLowNonZeroPerRow
 ) const
 {
     // communicator for the matrix
@@ -631,7 +638,20 @@ void Foam::petscSolver::buildMat
             ++lowNonZero[celli];
         }
 
-        maxLowNonZeroPerRow = max(lowNonZero);
+
+        // PetscInt may be an unusually sized type.
+        // Use max_element instead of Foam::max(List..)
+
+        maxLowNonZeroPerRow = 0;
+        {
+            auto maxElem =
+                std::max_element(lowNonZero.cbegin(), lowNonZero.cend());
+
+            if (maxElem != lowNonZero.cend())  // Extra safety
+            {
+                maxLowNonZeroPerRow = *maxElem;
+            }
+        }
 
         // Off-processor non-zero entries, per row.
 
@@ -721,8 +741,8 @@ void Foam::petscSolver::buildMat
 void Foam::petscSolver::updateMat
 (
     Mat& Amat,
-    List<label>& lowNonZero,
-    label& maxLowNonZeroPerRow
+    const List<PetscInt>& lowNonZero,
+    const PetscInt maxLowNonZeroPerRow
 ) const
 {
     const lduAddressing& lduAddr = matrix_.mesh().lduAddr();
@@ -794,7 +814,9 @@ void Foam::petscSolver::updateMat
        (*f)(Amat,coo_v,INSERT_VALUES);
 #endif
        PetscFree(coo_v);
-    } else {
+    }
+    else
+    {
        // Number of internal faces (connectivity)
        const label nIntFaces_ = upp.size();
 
@@ -832,7 +854,7 @@ void Foam::petscSolver::updateMat
            {
                globRow = globalNumbering_.toGlobal(low[kidx]);
 
-               for (label jidx=0; jidx<lowNonZero[idx]; ++jidx)
+               for (label jidx=0; jidx < lowNonZero[idx]; ++jidx)
                {
                    // The row is sorted, col is unsorted
                    globCols[jidx] = globalNumbering_.toGlobal(upp[kidx]);

src/petsc4Foam/solvers/petscSolver.H

@@ -5,7 +5,7 @@
     \\  /    A nd           | www.openfoam.com
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
-    Copyright (C) 2018-2020 OpenCFD Ltd.
+    Copyright (C) 2018-2022 OpenCFD Ltd.
     Copyright (C) 2019-2020 Simone Bna
 -------------------------------------------------------------------------------
 License
@@ -55,8 +55,8 @@ SourceFiles
 
 \*---------------------------------------------------------------------------*/
 
-#ifndef petscFoamSolver_H
-#define petscFoamSolver_H
+#ifndef Foam_module_petscSolver_H
+#define Foam_module_petscSolver_H
 
 #include "lduMatrix.H"
 
@@ -96,16 +96,16 @@ class petscSolver
         void buildMat
         (
             Mat& Amat,
-            List<label>& lowNonZero,
-            label& maxLowNonZeroPerRow
+            List<PetscInt>& lowNonZero,
+            PetscInt& maxLowNonZeroPerRow
         ) const;
 
         //- Update Amat_ matrix by inserting new values
         void updateMat
         (
             Mat& Amat,
-            List<label>& lowNonZero,
-            label& maxLowNonZeroPerRow
+            const List<PetscInt>& lowNonZero,
+            const PetscInt maxLowNonZeroPerRow
         ) const;
 
         //- Build the ksp_ krylov solver context

src/petsc4Foam/utils/petscCacheManager.H

@@ -61,8 +61,8 @@ SourceFiles
 
 \*---------------------------------------------------------------------------*/
 
-#ifndef petscFoamCacheManager_H
-#define petscFoamCacheManager_H
+#ifndef Foam_module_petscUtils_Caching_H
+#define Foam_module_petscUtils_Caching_H
 
 #include "Enum.H"
 #include "clockValue.H"

src/petsc4Foam/utils/petscControls.H

@@ -36,8 +36,8 @@ SourceFiles
 
 \*---------------------------------------------------------------------------*/
 
-#ifndef petscFoamControls_H
-#define petscFoamControls_H
+#ifndef Foam_module_petscControls_H
+#define Foam_module_petscControls_H
 
 #include "IOdictionary.H"
 #include "MeshObject.H"

src/petsc4Foam/utils/petscLinearSolverContext.H

@@ -5,7 +5,7 @@
     \\  /    A nd           | www.openfoam.com
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
-    Copyright (C) 2019 OpenCFD Ltd.
+    Copyright (C) 2019-2022 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -34,8 +34,8 @@ SourceFiles
 
 \*---------------------------------------------------------------------------*/
 
-#ifndef petscFoamLinearSolverContext_H
-#define petscFoamLinearSolverContext_H
+#ifndef Foam_module_petscLinearSolverContext_H
+#define Foam_module_petscLinearSolverContext_H
 
 #include "solverPerformance.H"
 #include "petscCacheManager.H"
@@ -66,10 +66,12 @@ public:
 
     // Public Members
 
+        //- Using OpenFOAM norm for convergence testing
+        bool useFoamTest;
+
         Mat Amat;
         KSP ksp;
 
-        bool useFoamTest;
         Vec ArowsSum;
         Vec *res_l1_w;
 
@@ -77,7 +79,7 @@ public:
         petscCacheManager caching;
 
         List<PetscInt> lowNonZero;
-        label maxLowNonZeroPerRow;
+        PetscInt maxLowNonZeroPerRow = 0;
 
         PetscScalar normFactor;
 
@@ -93,9 +95,9 @@ public:
         :
             init_(false),
             init_aux_vectors_(false),
+            useFoamTest(false),
             Amat(nullptr),
             ksp(nullptr),
-            useFoamTest(false),
             ArowsSum(nullptr),
             res_l1_w(nullptr)
         {}
@@ -114,15 +116,17 @@ public:
     // Member Functions
 
         //- Return value of initialized
-        bool initialized() const
+        bool initialized() const noexcept
         {
             return init_;
         }
 
-        //- Return value of initialized
-        bool& initialized()
+        //- Change state of initialized, return previous value
+        bool initialized(const bool yes) noexcept
         {
-            return init_;
+            bool old(init_);
+            init_ = yes;
+            return old;
         }
 
         //- Create auxiliary rows for calculation purposes
@@ -168,7 +172,7 @@ public:
         //- assumes normFactor has been already computed
         PetscReal getResidualNormL1(Vec psi, Vec source)
         {
-            if (!init_aux_vectors_) return -1.;
+            if (!init_aux_vectors_) return static_cast<PetscReal>(-1);
 
             PetscReal rnorm;
             MatMult(Amat, psi, res_l1_w[1]);
@@ -182,15 +186,25 @@ public:
         //- as well as Adotpsi is stored in res_l1_w[1]
         PetscReal getResidualNormL1(Vec source)
         {
-            if (!init_aux_vectors_) return -1.;
+            if (!init_aux_vectors_) return static_cast<PetscReal>(-1);
 
             PetscReal rnorm;
             VecAXPY(res_l1_w[1], -1., source);
             VecNorm(res_l1_w[1], NORM_1, &rnorm);
             return rnorm/normFactor;
         }
+
+
+    // Housekeeping
+
+        //- Non-const access to state of initialized
+        bool& initialized() noexcept
+        {
+            return init_;
+        }
 };
 
+
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
 } // End namespace Foam

src/petsc4Foam/utils/petscLinearSolverContexts.H

@@ -34,8 +34,8 @@ SourceFiles
 
 \*---------------------------------------------------------------------------*/
 
-#ifndef petscFoamLinearSolverContexts_H
-#define petscFoamLinearSolverContexts_H
+#ifndef Foam_module_petscLinearSolverContexts_H
+#define Foam_module_petscLinearSolverContexts_H
 
 #include "fvMesh.H"
 #include "MeshObject.H"

src/petsc4Foam/utils/petscUtils.H

@@ -35,8 +35,8 @@ SourceFiles
 
 \*---------------------------------------------------------------------------*/
 
-#ifndef petscFoamUtils_H
-#define petscFoamUtils_H
+#ifndef Foam_module_petscUtils_H
+#define Foam_module_petscUtils_H
 
 // OpenFOAM
 #include "scalarField.H"

src/petsc4Foam/utils/petscWrappedVector.H

@@ -33,8 +33,8 @@ SourceFiles
 
 \*---------------------------------------------------------------------------*/
 
-#ifndef petscFoamWrappedVector_H
-#define petscFoamWrappedVector_H
+#ifndef Foam_module_petscWrappedVector_H
+#define Foam_module_petscWrappedVector_H
 
 // OpenFOAM
 #include "List.H"