From 0759b72a32f7c98679ffead67d66bfb07c21d8c6 Mon Sep 17 00:00:00 2001
From: Laurent FACQ <laurent.facq@math.u-bordeaux.fr>
Date: Wed, 19 Jan 2022 14:45:50 +0100
Subject: [PATCH] Add solveAgainLaplacian() for multiple RHS solve with same
 matrix

---
 src/maths/Laplacian.hpp              |  1 +
 src/maths/solvers/LaplacianPetsc.cpp | 16 ++++++++++++----
 src/maths/solvers/LaplacianPetsc.hpp |  7 +++++++
 3 files changed, 20 insertions(+), 4 deletions(-)

diff --git a/src/maths/Laplacian.hpp b/src/maths/Laplacian.hpp
index 9b40f40..a8d9c68 100644
--- a/src/maths/Laplacian.hpp
+++ b/src/maths/Laplacian.hpp
@@ -64,6 +64,7 @@ virtual void addRHSValue(int row, double val) = 0;
 virtual double getRHSValue(const int row) = 0;
 virtual void zeroRHS() = 0;
 virtual double solveLaplacian() = 0;
+virtual double solveAgainLaplacian() = 0;
 virtual std::vector<double> getSolution() = 0;
 virtual std::vector<double> error() = 0;
 virtual void dumpMatrix() = 0;
diff --git a/src/maths/solvers/LaplacianPetsc.cpp b/src/maths/solvers/LaplacianPetsc.cpp
index b514841..57df2ba 100644
--- a/src/maths/solvers/LaplacianPetsc.cpp
+++ b/src/maths/solvers/LaplacianPetsc.cpp
@@ -37,8 +37,6 @@ void LaplacianPetsc::zeroMatrix()
 
 double LaplacianPetsc::solveLaplacian()
 {
-  PetscInt nbIteration;
-  PetscReal norm;
   PetscBool assembled;
 
   if (!MatAssembled(_LaplacianMatrix,&assembled)) {
@@ -46,8 +44,6 @@ double LaplacianPetsc::solveLaplacian()
     MatAssemblyEnd(_LaplacianMatrix, MAT_FINAL_ASSEMBLY);
   }
 
-  VecAssemblyBegin(_RHS);
-  VecAssemblyEnd(_RHS);
   VecDuplicate(_RHS, &_dum);
 
   KSPSetOperators(_ksp, _LaplacianMatrix, _LaplacianMatrix);
@@ -62,6 +58,18 @@ double LaplacianPetsc::solveLaplacian()
       MatSetNullSpace(_LaplacianMatrix, _nullsp);
     }
   }
+
+  return LaplacianPetsc::solveAgainLaplacian();
+}
+
+double LaplacianPetsc::solveAgainLaplacian()
+{
+  PetscInt nbIteration;
+  PetscReal norm;
+
+  VecAssemblyBegin(_RHS);
+  VecAssemblyEnd(_RHS);
+
   KSPSolve(_ksp, _RHS, _dum);
   KSPGetIterationNumber(_ksp, &nbIteration);
   KSPGetResidualNorm(_ksp, &norm);
diff --git a/src/maths/solvers/LaplacianPetsc.hpp b/src/maths/solvers/LaplacianPetsc.hpp
index 061f728..5006a07 100644
--- a/src/maths/solvers/LaplacianPetsc.hpp
+++ b/src/maths/solvers/LaplacianPetsc.hpp
@@ -111,6 +111,13 @@ void dumpRHS();
  */
 double solveLaplacian();
 
+/**
+ * @brief Solve again the LaplacianPetsc (same system, different RHS)
+ *
+ * @brief[out] return residual norm
+ */
+double solveAgainLaplacian();
+
 /**
  * @brief get Solution converted in std::vector
  */
-- 
GitLab