diff --git a/src/petsc4Foam/solvers/petscSolver.C b/src/petsc4Foam/solvers/petscSolver.C index c22d93b8eaa63278e342fd6e77fb31c295c1362c..e0eb9399a2ea3e04da3bb1b8365257fe6da9e7d0 100644 --- a/src/petsc4Foam/solvers/petscSolver.C +++ b/src/petsc4Foam/solvers/petscSolver.C @@ -248,18 +248,24 @@ Foam::solverPerformance Foam::petscSolver::scalarSolve PetscWrappedVector petsc_source(source, Amat); // Initialize data to compute L1 residual or to monitor it while solving - if (!usePetscResidualNorm || monitorFoamResidualNorm) + ctx.performance.initialResidual() = 0; + KSPNormType normType; + KSPGetNormType(ksp, &normType); + if (!usePetscResidualNorm || monitorFoamResidualNorm || normType == KSP_NORM_NONE) { - ctx.createAuxVecs(petsc_psi); - ctx.initArowsSumVec(Amat); - ctx.computeNormFactor(Amat,petsc_psi,petsc_source); + ctx.createAuxVecs(); + ctx.initArowsSumVec(); + ctx.computeNormFactor(petsc_psi,petsc_source); + if (normType == KSP_NORM_NONE && !usePetscResidualNorm) + { + ctx.performance.initialResidual() = ctx.getResidualNormL1(petsc_source); + } } // Convergence testing if (usePetscResidualNorm) { - // Add monitor to record initial residual - ctx.performance.initialResidual() = 0; + // Add monitor to record initial residual computed by PETSc KSPMonitorSet ( ksp, @@ -335,20 +341,31 @@ Foam::solverPerformance Foam::petscSolver::scalarSolve KSPSolve(ksp, petsc_source, petsc_psi); PetscLogStagePop(); + ctx.caching.eventEnd(); + // Set nIterations and final residual PetscInt nIters; KSPGetIterationNumber(ksp, &nIters); ctx.performance.nIterations() = nIters; + // If using the PETSc infrastructure to monitor convergence, + // we need to report back the final residual if (usePetscResidualNorm) { PetscReal rnorm; - KSPGetResidualNorm(ksp, &rnorm); + KSPNormType normType; + KSPGetNormType(ksp, &normType); + if (normType == KSP_NORM_NONE) // 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 + { + KSPGetResidualNorm(ksp, &rnorm); + } ctx.performance.finalResidual() = rnorm; } - ctx.caching.eventEnd(); - // Return solver performance to OpenFOAM return ctx.performance; } diff --git a/src/petsc4Foam/utils/petscLinearSolverContext.H b/src/petsc4Foam/utils/petscLinearSolverContext.H index ef6145a2494f1375e8419efea4c77d29389df452..08a491668bd842d70e50292385c4f19e578c5f59 100644 --- a/src/petsc4Foam/utils/petscLinearSolverContext.H +++ b/src/petsc4Foam/utils/petscLinearSolverContext.H @@ -125,32 +125,33 @@ public: } //- Create auxiliary rows for calculation purposes - void createAuxVecs(Vec psi) + void createAuxVecs() { + if (!Amat) return; if (!init_aux_vectors_) { + MatCreateVecs(Amat, NULL, &ArowsSum); + VecDuplicateVecs(ArowsSum, 2, &res_l1_w); init_aux_vectors_ = true; - VecDuplicate(psi, &ArowsSum); - VecDuplicateVecs(psi, 2, &res_l1_w); } } //- Sum the rows of A, placing the result in ArowsSum - void initArowsSumVec(Mat A) + void initArowsSumVec() { if (init_aux_vectors_) { - MatGetRowSum(A, ArowsSum); + MatGetRowSum(Amat, ArowsSum); } } //- Compute the normFactor used in convergence testing, //- assumes ArowsSum has been already computed - void computeNormFactor(Mat A, Vec psi, Vec source) + void computeNormFactor(Vec psi, Vec source) { if (!init_aux_vectors_) return; - MatMult(A, psi, res_l1_w[1]); + MatMult(Amat, psi, res_l1_w[1]); normFactor = Foam::PetscUtils::normFactor @@ -161,6 +162,32 @@ public: ArowsSum ); } + + //- Compute residual (L1) norm + //- assumes normFactor has been already computed + PetscReal getResidualNormL1(Vec psi, Vec source) + { + if (!init_aux_vectors_) return -1.; + + PetscReal rnorm; + MatMult(Amat, psi, res_l1_w[1]); + VecAXPY(res_l1_w[1], -1., source); + VecNorm(res_l1_w[1], NORM_1, &rnorm); + return rnorm/normFactor; + } + + //- Compute residual (L1) norm + //- assumes normFactor has been already computed + //- as well as Adotpsi is stored in res_l1_w[1] + PetscReal getResidualNormL1(Vec source) + { + if (!init_aux_vectors_) return -1.; + + PetscReal rnorm; + VecAXPY(res_l1_w[1], -1., source); + VecNorm(res_l1_w[1], NORM_1, &rnorm); + return rnorm/normFactor; + } }; // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //