Added: Allow using PETSc solvers with CGL2 version 2

src/ASM/GlbL2projector.C

@@ -21,6 +21,7 @@
 #include "IntegrandBase.h"
 #include "Function.h"
 #include "Profiler.h"
+#include "SparseMatrix.h"
 #ifdef HAS_PETSC
 #include "PETScMatrix.h"
 #include "LinSolParams.h"
@@ -106,37 +107,57 @@ private:
 };
 
 
-GlbL2::GlbL2 (IntegrandBase* p, size_t n) : A(SparseMatrix::SUPERLU)
+GlbL2::GlbL2 (IntegrandBase* p, size_t n) : problem(p)
 {
-  problem = p;
   nrhs = p->getNoFields(2);
-
-  A.redim(n,n);
-  B.redim(n*nrhs);
+  this->allocate(n);
 }
 
 
-GlbL2::GlbL2 (FunctionBase* f, size_t n) : A(SparseMatrix::SUPERLU)
+GlbL2::GlbL2 (FunctionBase* f, size_t n) : problem(nullptr), functions({f})
 {
-  problem = nullptr;
-  functions = { f };
   nrhs = f->dim();
-
-  A.redim(n,n);
-  B.redim(n*nrhs);
+  this->allocate(n);
 }
 
 
-GlbL2::GlbL2 (const FunctionVec& f, size_t n) : A(SparseMatrix::SUPERLU)
+GlbL2::GlbL2 (const FunctionVec& f, size_t n) : problem(nullptr), functions(f)
 {
-  problem = nullptr;
-  functions = f;
   nrhs = 0;
   for (FunctionBase* func : f)
     nrhs += func->dim();
+  this->allocate(n);
+}
 
-  A.redim(n,n);
-  B.redim(n*nrhs);
+
+GlbL2::~GlbL2()
+{
+  delete pA;
+  delete pB;
+#ifdef HAS_PETSC
+  delete adm;
+#endif
+}
+
+
+void GlbL2::allocate (size_t n)
+{
+#ifdef HAS_PETSC
+  adm = nullptr;
+  if (GlbL2::MatrixType == LinAlg::PETSC && GlbL2::SolverParams)
+  {
+    adm = new ProcessAdm();
+    pA = new PETScMatrix(*adm,*GlbL2::SolverParams);
+    pB = new PETScVector(*adm,n*nrhs);
+  }
+  else
+#endif
+  {
+    pA = new SparseMatrix(SparseMatrix::SUPERLU);
+    pB = new StdVector(n*nrhs);
+  }
+
+  pA->redim(n,n);
 }
 
 
@@ -247,12 +268,15 @@ bool GlbL2::evalIntMx (LocalIntegral& elmInt,
 
 void GlbL2::preAssemble (const std::vector<IntVec>& MMNPC, size_t nel)
 {
-  A.preAssemble(MMNPC,nel);
+  pA->preAssemble(MMNPC,nel);
 }
 
 
 bool GlbL2::solve (Matrix& sField)
 {
+  SparseMatrix& A = *pA;
+  StdVector&    B = *pB;
+
   // Insert a 1.0 value on the diagonal for equations with no contributions.
   // Needed in immersed boundary calculations with "totally outside" elements.
   size_t i, j, nnod = A.dim();
@@ -289,6 +313,9 @@ bool GlbL2::solve (const std::vector<Matrix*>& sField)
     return false;
   }
 
+  SparseMatrix& A = *pA;
+  StdVector&    B = *pB;
+
   // Insert a 1.0 value on the diagonal for equations with no contributions.
   // Needed in immersed boundary calculations with "totally outside" elements.
   size_t i, j, nnod = A.dim();
@@ -330,7 +357,7 @@ bool ASMbase::L2projection (Matrix& sField,
   PROFILE2("ASMbase::L2projection");
 
   GlbL2 gl2(integrand,this->getNoNodes(1));
-  L2GlobalInt dummy(gl2.A,gl2.B);
+  L2GlobalInt dummy(*gl2.pA,*gl2.pB);
 
   gl2.preAssemble(MNPC,this->getNoElms(true));
   return this->integrate(gl2,dummy,time) && gl2.solve(sField);
@@ -342,7 +369,7 @@ bool ASMbase::L2projection (Matrix& sField, FunctionBase* function, double t)
   PROFILE2("ASMbase::L2projection");
 
   GlbL2 gl2(function,this->getNoNodes(1));
-  L2GlobalInt dummy(gl2.A,gl2.B);
+  L2GlobalInt dummy(*gl2.pA,*gl2.pB);
   TimeDomain time; time.t = t;
 
   gl2.preAssemble(MNPC,this->getNoElms(true));
@@ -356,7 +383,7 @@ bool ASMbase::L2projection (const std::vector<Matrix*>& sField,
   PROFILE2("ASMbase::L2projection");
 
   GlbL2 gl2(function,this->getNoNodes(1));
-  L2GlobalInt dummy(gl2.A,gl2.B);
+  L2GlobalInt dummy(*gl2.pA,*gl2.pB);
   TimeDomain time; time.t = t;
 
   gl2.preAssemble(MNPC,this->getNoElms(true));

src/ASM/GlbL2projector.h

@@ -15,12 +15,17 @@
 #define _GLB_L2_PROJECTOR_H
 
 #include "Integrand.h"
-#include "SparseMatrix.h"
+#include "LinAlgenums.h"
+#include "MatVec.h"
 
 class IntegrandBase;
 class FunctionBase;
 class LinSolParams;
+class SparseMatrix;
+class StdVector;
+class ProcessAdm;
 
+typedef std::vector<int>           IntVec;      //!< Vector of integers
 typedef std::vector<size_t>        uIntVec;     //!< Vector of unsigned integers
 typedef std::vector<FunctionBase*> FunctionVec; //!< Vector of functions
 
@@ -47,8 +52,8 @@ public:
   //! \param[in] f The functions to do L2-projection on
   //! \param[in] n Dimension of the L2-projection matrices (number of nodes)
   GlbL2(const FunctionVec& f, size_t n);
-  //! \brief Empty destructor.
-  virtual ~GlbL2() {}
+  //! \brief The destructor frees the system matrix and system vector.
+  virtual ~GlbL2();
 
   //! \brief Defines which FE quantities are needed by the integrand.
   virtual int getIntegrandType() const;
@@ -116,14 +121,21 @@ public:
   //! \param[out] sField Nodal/control-point values of the projected results.
   bool solve(const std::vector<Matrix*>& sField);
 
+private:
+  //! \brief Allocates the system L2-projection matrices.
+  void allocate(size_t n);
+
 public:
-  mutable SparseMatrix A; //!< Left-hand-side matrix of the L2-projection
-  mutable StdVector    B; //!< Right-hand-side vectors of the L2-projection
+  mutable SparseMatrix* pA; //!< Left-hand-side matrix of the L2-projection
+  mutable StdVector*    pB; //!< Right-hand-side vectors of the L2-projection
 
 private:
   IntegrandBase* problem; //!< The main problem integrand
   FunctionVec  functions; //!< Explicit functions to L2-project
   size_t            nrhs; //!< Number of right-hand-size vectors
+#ifdef HAS_PETSC
+  ProcessAdm* adm; //!< Process administrator for PETSc
+#endif
 };
 
 #endif

src/LinAlg/PETScMatrix.h

@@ -105,9 +105,6 @@ public:
   //! \brief Returns the matrix type.
   virtual LinAlg::MatrixType getType() const { return LinAlg::PETSC; }
 
-  //! \brief Returns the dimension of the system matrix.
-  virtual size_t dim(int = 1) const { return 0; }
-
   //! \brief Initializes the element assembly process.
   //! \details Must be called once before the element assembly loop.
   //! The PETSc data structures are initialized and the all symbolic operations