diff --git a/src/petsc4Foam/solvers/petscSolver.C b/src/petsc4Foam/solvers/petscSolver.C
index b9f07abccc983aa61bb7ce75e32a02e9232c3b39..a55c35d8287db834ea7b04fcc61050f22d8ff6bb 100644
--- a/src/petsc4Foam/solvers/petscSolver.C
+++ b/src/petsc4Foam/solvers/petscSolver.C
@@ -166,14 +166,24 @@ Foam::solverPerformance Foam::petscSolver::scalarSolve
         ctx.initialized() = true;
         needsCacheUpdate = false;
 
+        PetscLogStageRegister
+        (
+            ("foam_" + eqName_ + "_mat").c_str(),
+            &ctx.matstage
+        );
+        PetscLogStageRegister
+        (
+            ("foam_" + eqName_ + "_ksp").c_str(),
+            &ctx.kspstage
+        );
         computeMatAllocation(Amat, lowNonZero, maxLowNonZeroPerRow);
-        buildMat(Amat, lowNonZero, maxLowNonZeroPerRow);
+        buildMat(Amat, lowNonZero, maxLowNonZeroPerRow, ctx.matstage);
         buildKsp(Amat, ksp);
     }
 
     if (ctx.caching.needsMatrixUpdate() && needsCacheUpdate)
     {
-        buildMat(Amat, lowNonZero, maxLowNonZeroPerRow);
+        buildMat(Amat, lowNonZero, maxLowNonZeroPerRow, ctx.matstage);
         updateKsp(ksp, Amat, !ctx.caching.needsPrecondUpdate());
     }
 
@@ -269,6 +279,8 @@ Foam::solverPerformance Foam::petscSolver::scalarSolve
        );
     }
 
+    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)
     KSPSetUp(ksp);
@@ -276,6 +288,8 @@ Foam::solverPerformance Foam::petscSolver::scalarSolve
     // Solve A x = b
     KSPSolve(ksp, petsc_source, petsc_psi);
 
+    PetscLogStagePop();
+
     // Set nIterations and final residual
     PetscInt nIters;
     KSPGetIterationNumber(ksp, &nIters);
@@ -527,9 +541,12 @@ void Foam::petscSolver::buildMat
 (
     Mat& Amat,
     List<label>& lowNonZero,
-    label& maxLowNonZeroPerRow
+    label& maxLowNonZeroPerRow,
+    PetscLogStage stage
 ) const
 {
+    PetscLogStagePush(stage);
+
     const lduAddressing& lduAddr = matrix_.mesh().lduAddr();
     const lduInterfacePtrsList interfaces(matrix_.mesh().interfaces());
 
@@ -710,6 +727,7 @@ void Foam::petscSolver::buildMat
 
     MatAssemblyBegin(Amat, MAT_FINAL_ASSEMBLY);
     MatAssemblyEnd(Amat, MAT_FINAL_ASSEMBLY);
+    PetscLogStagePop();
 }
 
 
diff --git a/src/petsc4Foam/solvers/petscSolver.H b/src/petsc4Foam/solvers/petscSolver.H
index fd07549b43ac25a0b2bcd7f0aa30dc7067922549..d0287a55c4e463661fa72e2aa1c099238bc088fa 100644
--- a/src/petsc4Foam/solvers/petscSolver.H
+++ b/src/petsc4Foam/solvers/petscSolver.H
@@ -105,7 +105,8 @@ class petscSolver
         (
             Mat& Amat,
             List<label>& lowNonZero,
-            label& maxLowNonZeroPerRow
+            label& maxLowNonZeroPerRow,
+            PetscLogStage stage
         ) const;
 
         //- Build the ksp_ krylov solver context
diff --git a/src/petsc4Foam/utils/petscLinearSolverContext.H b/src/petsc4Foam/utils/petscLinearSolverContext.H
index 75aaf75644bfc8b3b251d60536d469b2c7804949..7eb25b7e040e244509f701d7f3c36883241ca5fa 100644
--- a/src/petsc4Foam/utils/petscLinearSolverContext.H
+++ b/src/petsc4Foam/utils/petscLinearSolverContext.H
@@ -69,6 +69,9 @@ public:
         Mat Amat;
         KSP ksp;
 
+        PetscLogStage matstage;
+        PetscLogStage kspstage;
+
         solverPerformance performance;
         petscCacheManager caching;
 
@@ -143,7 +146,8 @@ public:
             }
         }
 
-        //- Compute the normFactor used in convergence testing, assumes ArowsSum has been already computed
+        //- Compute the normFactor used in convergence testing,
+        //- assumes ArowsSum has been already computed
         void computeNormFactor(Mat A, Vec psi, Vec source)
         {
             if (!init_aux_vectors_) return;