Add initial set of comments

src/petsc/csrMatrix/solvers/petscSolver.C

@@ -113,6 +113,8 @@ Foam::solverPerformance Foam::petscSolver::solve
 
     bool needsCacheUpdate = true;
 
+// SZ logic about this?
+// SZ remeshing sets initialized to false?
     if (!ctx.initialized())
     {
         ctx.initialized() = true;
@@ -126,11 +128,24 @@ Foam::solverPerformance Foam::petscSolver::solve
         Info<< "Initializing PETSc Linear Solver " << eqName_ << nl;
     }
 
+// SZ what 'caching' actually flags?
     const bool caching
     (
         petscDict_.lookupOrDefault("caching", true)
     );
 
+// SZ update matrix coefficients only?
+// SZ Is openfoam solving a linear system without updating the linear operator?
+// PETSc has support for updating the operators but not setup the new preconditioner
+// via KSPSetReusePreconditioner
+// SZ: PROPOSAL
+// At this point, we don't know if the values stored in the lduMatrix has been updated
+// or not. If we reuse the Mat object (and we want to do this, unless the sparsity pattern
+// has changed, e.g. via remeshing), then it would be nice to know whether we need to
+// reset the values in Amat or not. A possible solution is to add a int64  variable to the lduMatrix
+// that gets incremented every time its values are updated. This solver class will store that value
+// and check here if the value stored is the same of the current one in the lduMatrix.
+// If so, don't set values in Amat; otherwise set the new values and store the lduMatrix counter.
     if (!caching && needsCacheUpdate)
     {
         buildMat(Amat);
@@ -139,6 +154,10 @@ Foam::solverPerformance Foam::petscSolver::solve
         Info<< "Update cached PETSc Linear Solver " << eqName_ << nl;
     }
 
+// SZ It would be nice to add a prefix to the solver
+// I imagine different solvers may need different PETSc prefixes to experiment
+// via rcfile for quick and fine tuning, e.g.
+// -sub_foamsolverprefix_pc_type xxx
 
     // Solver name from petsc dictionary
     const word solverName
@@ -151,11 +170,14 @@ Foam::solverPerformance Foam::petscSolver::solve
         petscDict_.lookupOrDefault("verbose", false)
     );
 
+// SZ this should be removed
     const bool ksp_view
     (
         petscDict_.lookupOrDefault("ksp_view", false)
     );
 
+// SZ it appears the norm used inside KSP is always the 2 norm, not the openfoam one
+// SZ only the initial/final values of the openfoam norm is computed 
     // Use built-in residual norm computation
     const bool usePetscResidualNorm
     (
@@ -192,6 +214,7 @@ Foam::solverPerformance Foam::petscSolver::solve
     // Initial residual (scaled)
     if (usePetscResidualNorm)
     {
+// SZ why zero?
         solverPerf.initialResidual() = 0.;
     }
     else
@@ -203,6 +226,10 @@ Foam::solverPerformance Foam::petscSolver::solve
     // Solve A x = b
     KSPSolve(ksp, petsc_source, petsc_psi);
 
+// SZ why explicitly managing it here? Also, the View just give an idea of the operators involved
+// If you want to do this, then this should be moved after the first setup time of the KSP
+// PETSc's ksp_view instead already prints after each solve
+// REMOVE PETSC_VIEWER_STDOUT_SELF, USE NULL 
     // View the ksp
     if (ksp_view)
     {
@@ -249,7 +276,7 @@ Foam::solverPerformance Foam::petscSolver::solve
     return solverPerf;
 }
 
-
+// SZ does this happen with different matrix size?
 void Foam::petscSolver::updateKsp
 (
     Mat& Amat,
@@ -300,9 +327,12 @@ void Foam::petscSolver::buildKsp
     dictionary subPrecondDict = petscDict_.subOrEmptyDict(subPrecondName + "Coeffs");
     dictionary solverDict  = petscDict_.subOrEmptyDict(solverName + "Coeffs");
 
+
+// SZ Usage of PETSC_COMM_WORLD should be avoided, can we get the comm from the mesh?
     // Create parallel solver context
     KSPCreate(PETSC_COMM_WORLD, &ksp);
 
+// SZ since we have an updateksp method, we should use that for consistency
     // ksp set operator and preconditioner
     KSPSetOperators(ksp, Amat, Amat);
 
@@ -310,6 +340,19 @@ void Foam::petscSolver::buildKsp
     PC pc;
     KSPGetPC(ksp, &pc);
     PCSetType(pc, precondName.c_str());
+// SZ If we use the prefix (e.g the name of the solver component, like 'p' or 'U' etc), we need some changes here
+// We need to add KSPSetOptionsPrefix(ksp,prefix) after the creation of the KSP object
+// My proposal is to have explicit options set in the dictionaries precondDict and solverDict
+// and get rid of subPreconditionerDict (which is specific to PCBJACOBI only)
+// and factorDict (specific to PCFACTOR only)
+// As an example, for the pressure (prefix=p) solver in the preconditioner coeffs for bjacobi, one will explicitly set
+// sub_p_pc_type=lu
+// sub_p_pc_factor_mat_solver_type=mumps
+// sub_p_ksp_type=bcgs
+// and those flags are set directly with one call 
+// PetscUtils::setSolverOptions(precondDict, verbose);
+// the user will have much more control on the various solvers, as the options are passed directly to PETSc
+// and we don't need to handle corner cases ourselves.
     PetscUtils::setFlag("-sub_pc_type", subPrecondName, verbose);
 
     PetscUtils::setSolverFlags("-pc_" + precondName, precondDict, verbose);
@@ -330,7 +373,10 @@ void Foam::petscSolver::buildKsp
         maxIter_
     );
 
+// SZ if we want to use openfoam residual computation in KSP, we should use the KSPSetConvergenceTest
+
     // Use the solution from the previous timestep if the solver is not preonly
+// SZ PETSc has also support to compute initial guesses using POD when matrices are fixed for few time steps
     if (solverName != "preonly")
     {
         KSPSetInitialGuessNonzero(ksp, PETSC_TRUE);
@@ -399,11 +445,18 @@ void Foam::petscSolver::computeMatAllocation
 
 
     // Create and allocate petsc matrix
+// SZ PETSC_COMM_WORLD
     MatCreate(PETSC_COMM_WORLD, &Amat);
     MatSetType(Amat, MATMPIAIJ);
     MatSetFromOptions(Amat);
     MatSetSizes(Amat, nrows_, nrows_, PETSC_DECIDE, PETSC_DECIDE);
-    MatMPIAIJSetPreallocation
+// SZ do we want to support MATIS? Needs l2gmap to be computed (needs hash table with visited global neighs
+
+// SZ others should be added (e.g. MATHYPRE)
+// I would use MatXAIJSetPreallocation here
+// Also, it seems the vector may use components, but the matrix will not
+// block size should be set here
+     MatMPIAIJSetPreallocation
     (
         Amat,
         0, nonZero_.data(),      // on-processor
@@ -490,6 +543,7 @@ void Foam::petscSolver::buildMat
         PetscInt globCol = globalNumbering_.toGlobal(upp[idx]);
         PetscScalar val = uppVal[idx];
 
+// SZ is there a way to set full rows here? or this information is totally lost at this point?
         MatSetValues
         (
             Amat,
@@ -508,6 +562,7 @@ void Foam::petscSolver::buildMat
         PetscInt globCol = globalNumbering_.toGlobal(low[idx]);
         PetscScalar val = lowVal[idx];
 
+// SZ is there a way to set full rows here? or this information is totally lost at this point?
         MatSetValues
         (
             Amat,
@@ -577,6 +632,7 @@ void Foam::petscSolver::buildMat
                     PetscInt globCol = nbrCells[i];  // <- already global
                     PetscScalar val = -bCoeffs[i];
 
+// SZ is there a way to set full rows here? or this information is totally lost at this point?
                     MatSetValues
                     (
                         Amat,

src/petsc/csrMatrix/solvers/petscSolver.H

@@ -61,6 +61,7 @@ SourceFiles
 #define petscSolver_H
 
 #include "lduMatrix.H"
+// SZ forward declare PETSc objects instead of including petsc.h?
 #include "petsc.h"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

src/petsc/utils/PetscUtils.C

@@ -122,6 +122,7 @@ void Foam::PetscUtils::copyVector
     VecAssemblyEnd(output);
 }
 
+// SZ this function should be rewritten by using MatMultAdd and not allocate  
 Foam::scalar Foam::PetscUtils::calcResidual
 (
     Mat A,
@@ -152,6 +153,8 @@ Foam::scalar Foam::PetscUtils::calcResidual
     return resid;
 }
 
+// SZ this function should be rewritten since it uses a lot of separate blas-1 kernels that should be fused together
+// usage of VecDuplicate should be avoided
 Foam::scalar Foam::PetscUtils::normFactor
 (
     Mat A,
@@ -207,6 +210,7 @@ Foam::scalar Foam::PetscUtils::normFactor
     return normFactor;
 }
 
+// this function does not seem to be used
 PetscErrorCode Foam::PetscUtils::foamKSPConverge
 (
     KSP ksp,

src/petsc/utils/PetscWrappedVector.H

@@ -66,6 +66,8 @@ class PetscWrappedVector
         template<class T>
         void createImpl(const UList<T>& list)
         {
+// SZ thinking about GPUs: this constructor should be replaced with a more standard
+// VecCreate + VecSetSize + VecSetType + VecPlaceArray sequence
             VecCreateMPIWithArray
             (
                 PETSC_COMM_WORLD,

src/petsc/utils/petscControls.C

@@ -42,6 +42,10 @@ bool Foam::petscControls::loaded_ = false;
 
 void Foam::petscControls::start(const fileName& optionsFile)
 {
+// SZ: what if PETSc is already initialized (not by openFOAM)?
+// this should be replaced by
+// loaded_ = !PetscInitializeCalled
+// Also ,is openfoam always using MPI_COMM_WORLD?
     if (!loaded_)
     {
         #if OPENFOAM < 1904
@@ -92,6 +96,7 @@ void Foam::petscControls::stop()
     }
     else
     {
+// SZ change the comment: no need to finalize PETSc
         Info<< "PETSc already finalized" << nl;
     }
 }

src/petsc/utils/petscLinearSolverContext.H

@@ -36,6 +36,8 @@ SourceFiles
 #define petscLinearSolverContext_H
 
 // PETSc
+// SZ: I prefer forward declaring KSP and Mat instead of including the whole PETSc here.
+// This should be done consistently for all public headers
 #include "petsc.h"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

[Mark OLESEN]

No worries about this. The petscSolver.H will never be included by any user coding. The petscSolver is called using the OpenFOAM run-time selection mechanism (see lines 14-20 in the petscSolver.C file).