Fix for the GPU case + performance improvement for CG and foam L1 norm computation

For the GPU case: if we do not call VecResetArray, the placed CPU array will never be updated

The residual is now accessed by KSPBuildResidual which, for CG, will consist in a VecCopy call, see this MR in PETSc https://gitlab.com/petsc/petsc/-/merge_requests/2860 We no longer need a patchy PETSc to compute the L1 norm inside CG

Added support to monitor the L1 norm / normFactor convergence of the solver It also prints out the correct initial residual when using PETSc convergence test

Removed some no longer needed functions

src/petsc4Foam/solvers/petscSolver.C

@@ -175,17 +175,38 @@ Foam::solverPerformance Foam::petscSolver::solve
         updateKsp(ksp, Amat, !ctx.caching.needsPrecondUpdate());
     }
 
-    // Solver name from petsc
-    KSPType ksptype;
-    KSPGetType(ksp, &ksptype);
-    const word solverName(ksptype);
-
     // Use built-in residual norm computation
     const bool usePetscResidualNorm
     (
         petscDict_.getOrDefault("use_petsc_residual_norm", false)
     );
 
+    // Monitor built-in residual
+    const bool monitorFoamResidualNorm
+    (
+        petscDict_.getOrDefault("monitor_foam_residual_norm", false)
+    );
+
+
+    // Disable KSP default computation of residual norm if we are
+    // monitoring convergence a-la OpenFOAM
+    if (!usePetscResidualNorm)
+    {
+        KSPSetNormType(ksp, KSP_NORM_NONE);
+    }
+    else
+    {
+        KSPSetNormType(ksp, KSP_NORM_DEFAULT);
+    }
+
+    // ksp set options from db (may change norm type here if needed)
+    KSPSetFromOptions(ksp);
+
+    // Solver name from petsc (may have been changed from option database)
+    KSPType ksptype;
+    KSPGetType(ksp, &ksptype);
+    const word solverName(ksptype);
+
     // Setup class containing solver performance data
     solverPerformance solverPerf
     (
@@ -195,38 +216,66 @@ Foam::solverPerformance Foam::petscSolver::solve
 
     ctx.solverPerf = &solverPerf;
 
-    // Create solution vector
+    // Create solution and rhs vectors for PETSc
+    // We could create these once, and call PlaceArray/ResetArray instead
     PetscWrappedVector petsc_psi(psi, Amat);
     PetscWrappedVector petsc_source(source, Amat);
 
-    // Initial residual (scaled)
+    // Initialize data to compute L1 residual or to monitor it while solving
+    if (!usePetscResidualNorm || monitorFoamResidualNorm)
+    {
+        ctx.createAuxVecs(petsc_psi);
+        ctx.initArowsSumVec(Amat);
+        ctx.computeNormFactor(Amat,petsc_psi,petsc_source);
+    }
+
+    // Convergence testing
     if (usePetscResidualNorm)
     {
+        // Add monitor to record initial residual
         solverPerf.initialResidual() = 0;
+        KSPMonitorSet
+        (
+            ksp,
+            PetscUtils::foamKSPMonitorRecordInit,
+            &ctx,
+            NULL
+        );
     }
     else
     {
-        ctx.createAuxVecs(petsc_psi);
-        ctx.initArowsSumVec(Amat);
-
         KSPSetConvergenceTest
         (
             ksp,
             PetscUtils::foamKSPConverge,
             &ctx,
-            nullptr
+            NULL
         );
     }
 
+    // Monitoring (print to stdout) residual reduction in OpenFOAM norm
+    if (monitorFoamResidualNorm)
+    {
+       KSPMonitorSet
+       (
+           ksp,
+           PetscUtils::foamKSPMonitorFoam,
+           &ctx,
+           NULL
+       );
+    }
+
+    // Setup KSP (this is not explicitly needed, but we call it to separate
+    // PCSetUp from KSPSolve timings when requesting -log_view from PETSc)
+    KSPSetUp(ksp);
+
     // Solve A x = b
     KSPSolve(ksp, petsc_source, petsc_psi);
 
-    // Set nIterations and compute final residual
-    {
-        PetscInt nIters;
-        KSPGetIterationNumber(ksp, &nIters);
-        solverPerf.nIterations() = nIters;
-    }
+    // Set nIterations and final residual
+    PetscInt nIters;
+    KSPGetIterationNumber(ksp, &nIters);
+    solverPerf.nIterations() = nIters;
 
     if (usePetscResidualNorm)
     {
@@ -282,15 +331,6 @@ void Foam::petscSolver::buildKsp
     KSP& ksp
 ) const
 {
-    // Redirected petsc preconditioner and solver names
-    // In the future may wish to support solver config via dictionary
-    // instead of a "Coeffs" dictionary
-
-    const bool verbose
-    (
-        petscDict_.getOrDefault("verbose", false)
-    );
-
     // Create parallel solver context
     KSPCreate(PETSC_COMM_WORLD, &ksp);
 
@@ -300,10 +340,6 @@ void Foam::petscSolver::buildKsp
     // ksp set operator and preconditioner
     KSPSetOperators(ksp, Amat, Amat);
 
-    // Set the preconditioner
-    PC pc;
-    KSPGetPC(ksp, &pc);
-
     // OpenFOAM relative tolerance -> tolerance_
     KSPSetTolerances
     (
@@ -314,20 +350,8 @@ void Foam::petscSolver::buildKsp
         maxIter_
     );
 
-    // ksp set options from db
-    KSPSetFromOptions(ksp);
-
-    KSPType ksptype;
-    KSPGetType(ksp, &ksptype);
-    const word solverName(ksptype);
-
-    // Use solution from the previous timestep if the solver is not 'preonly'
-    if (solverName != "preonly")
-    {
-        KSPSetInitialGuessNonzero(ksp, PETSC_TRUE);
-    }
-
-    KSPSetUp(ksp);
+    // Use solution from the previous timestep as initial guess
+    KSPSetInitialGuessNonzero(ksp, PETSC_TRUE);
 }
 
 
@@ -431,6 +455,14 @@ void Foam::petscSolver::computeMatAllocation
         nonZero_.data(),
         nonZeroNeigh_.data()
     );
+
+    // The general AIJ preallocator do not handle SELL
+    MatSeqSELLSetPreallocation
+    (
+        Amat,
+        0, nonZero_.data()        // on-processor
+    );
+
     MatMPISELLSetPreallocation
     (
         Amat,
@@ -450,6 +482,9 @@ void Foam::petscSolver::computeMatAllocation
     }
 
     MatSetOption(Amat, MAT_KEEP_NONZERO_PATTERN, PETSC_TRUE);
+
+    // Allow general setup of other matrices instead of erroring
+    MatSetUp(Amat);
 }