dg implicit : newton + euler implicit + ok with more than one field + add options on linear system

Solver/function.cpp

@@ -631,3 +631,6 @@ functionConstant *function::getDT() {
     _dtFunction = functionConstantNew(0.);
   return _dtFunction;
 }
+function *getFunctionCoordinates(){
+  return functionCoordinates::get();
+}

Solver/function.h

@@ -242,4 +242,5 @@ class functionConstant : public function {
   }
   void set(double val); 
 };
+function *getFunctionCoordinates();
 #endif

Solver/linearSystem.cpp

@@ -43,3 +43,9 @@ void linearSystem<double>::registerBindings(binding *b){
   linearSystemPETScRegisterBindings (b);
 #endif
 }
+
+void linearSystemBase::setParameter (std::string key, std::string value) {
+  if (isAllocated())
+    Msg::Error("this system is already allocated, parameters cannot be set");
+  _parameters[key] = value;
+}

Solver/linearSystem.h

@@ -5,27 +5,35 @@
 
 #ifndef _LINEAR_SYSTEM_H_
 #define _LINEAR_SYSTEM_H_
+#include <map>
+#include <string>
 
 // A class that encapsulates a linear system solver interface :
 // building a sparse matrix, solving a linear system
 class binding;
+class linearSystemBase {
+  protected:
+  std::map<std::string, std::string> _parameters;
+  public :
+  virtual bool isAllocated() const = 0;
+  virtual void allocate(int nbRows) = 0;
+  virtual void clear() = 0;
+  virtual void zeroMatrix() = 0;
+  virtual void zeroRightHandSide() = 0;
+  virtual int systemSolve() = 0;
+  void setParameter (std::string key, std::string value);
+};
 
 template <class scalar>
-class linearSystem {
+class linearSystem : public linearSystemBase {
  public :
   linearSystem (){}
   virtual ~linearSystem (){}
-  virtual bool isAllocated() const = 0;
-  virtual void allocate(int nbRows) = 0;
-  virtual void clear() = 0;
   virtual void addToMatrix(int _row, int _col, const scalar &val) = 0;
   virtual void getFromMatrix(int _row, int _col, scalar &val) const = 0;
   virtual void addToRightHandSide(int _row, const scalar &val) = 0;
   virtual void getFromRightHandSide(int _row, scalar &val) const = 0;
   virtual void getFromSolution(int _row, scalar &val) const = 0;
-  virtual void zeroMatrix() = 0;
-  virtual void zeroRightHandSide() = 0;
-  virtual int systemSolve() = 0;
   static void registerBindings (binding*);
   virtual double normInfRightHandSide() const = 0;
 };

Solver/linearSystemPETSc.cpp

@@ -3,6 +3,8 @@
 #if defined(HAVE_PETSC)
 #include "linearSystemPETSc.h"
 #include "fullMatrix.h"
+#include <stdlib.h>
+#include "GmshMessage.h"
 
 template <>
 void linearSystemPETSc<fullMatrix<PetscScalar> >::addToMatrix(int row, int col, const fullMatrix<PetscScalar> &val)
@@ -56,6 +58,7 @@ void linearSystemPETSc<fullMatrix<PetscScalar> >::getFromRightHandSide(int row,
   }
   _try(VecRestoreArray(_b, &tmp));
 }
+
 template<>
 void linearSystemPETSc<fullMatrix<PetscScalar> >::getFromSolution(int row, fullMatrix<PetscScalar> &val) const
 {
@@ -75,6 +78,9 @@ void linearSystemPETSc<fullMatrix<PetscScalar> >::getFromSolution(int row, fullM
 template<>
 void linearSystemPETSc<fullMatrix<PetscScalar> >::allocate(int nbRows) 
 {
+  _blockSize = strtol (_parameters["blockSize"].c_str(), NULL, 10);
+  if (_blockSize == 0)
+    Msg::Error ("'blockSize' parameters must be set for linearSystemPETScBlock");
   clear();
   _try(MatCreate(PETSC_COMM_WORLD, &_a)); 
   _try(MatSetSizes(_a, PETSC_DECIDE, PETSC_DECIDE, nbRows * _blockSize, nbRows * _blockSize));
@@ -108,10 +114,9 @@ void linearSystemPETScRegisterBindings(binding *b)
   cm->setArgNames(NULL);
   cb = b->addClass<linearSystemPETSc<fullMatrix<PetscScalar> > >("linearSystemPETScBlock");
   cb->setDescription("A linear system solver, based on PETSc");
-  cm = cb->setConstructor<linearSystemPETSc<fullMatrix<PetscScalar> >, int>();
-  cm->setDescription ("A new PETScBlock<PetscScalar> solver (we probably should get rid of the blockSize argument)");
+  cm = cb->setConstructor<linearSystemPETSc<fullMatrix<PetscScalar> > >();
+  cm->setDescription ("A new PETScBlock<PetscScalar> solver");
   cb->setParentClass<linearSystem<fullMatrix<PetscScalar> > >();
-  cm->setArgNames("blockSize", NULL);
 #endif // FIXME
 
 }

Solver/linearSystemPETSc.h

@@ -45,13 +45,32 @@
 template <class scalar>
 class linearSystemPETSc : public linearSystem<scalar> {
   int _blockSize; // for block Matrix
-  bool _isAllocated;
+  bool _isAllocated, _kspAllocated;
   Mat _a;
   Vec _b, _x;
+  KSP _ksp;
   void _try(int ierr) const { CHKERRABORT(PETSC_COMM_WORLD, ierr); }
+  void _kspCreate() {
+    _try(KSPCreate(PETSC_COMM_WORLD, &_ksp));
+    PC pc;
+    _try(KSPGetPC(_ksp, &pc));
+    // set some default options
+    _try(PCSetType(pc, PCILU));
+    _try(PCFactorSetMatOrderingType(pc, MATORDERING_RCM));
+    _try(PCFactorSetLevels(pc, 1));
+    _try(KSPSetTolerances(_ksp, 1.e-8, PETSC_DEFAULT, PETSC_DEFAULT, PETSC_DEFAULT));
+    // override the default options with the ones from the option
+    // database (if any)
+    _try(KSPSetFromOptions(_ksp));
+    _kspAllocated = true;
+  }
  public:
-  linearSystemPETSc(int blockSize = -1) : _isAllocated(false) {_blockSize = blockSize;}
-  virtual ~linearSystemPETSc(){ clear(); }
+  linearSystemPETSc() : _isAllocated(false) {_blockSize = 0; _kspAllocated = false;}
+  virtual ~linearSystemPETSc(){ 
+    clear(); 
+    if (_kspAllocated)
+      _try (KSPDestroy (_ksp));
+  }
   virtual bool isAllocated() const { return _isAllocated; }
   virtual void allocate(int nbRows)
   {
@@ -132,36 +151,39 @@ class linearSystemPETSc : public linearSystem<scalar> {
   }
   virtual void zeroMatrix()
   {
+    if (_isAllocated) {
+      _try(MatAssemblyBegin(_a, MAT_FINAL_ASSEMBLY));
+      _try(MatAssemblyEnd(_a, MAT_FINAL_ASSEMBLY));
       _try(MatZeroEntries(_a));
     }
+  }
   virtual void zeroRightHandSide()
   {
+    if (_isAllocated) {
+      _try(VecAssemblyBegin(_b));
+      _try(VecAssemblyEnd(_b));
       _try(VecZeroEntries(_b));
     }
-  virtual int systemSolve()
+  }
+  int systemSolve()
   {
+    if (!_kspAllocated)
+      _kspCreate();
+    if (linearSystem<scalar>::_parameters["matrix_reuse"] == "same_sparsity")
+      _try(KSPSetOperators(_ksp, _a, _a, SAME_NONZERO_PATTERN));
+    else if (linearSystem<scalar>::_parameters["matrix_reuse"] == "same_matrix")
+      _try(KSPSetOperators(_ksp, _a, _a, SAME_PRECONDITIONER));
+    else 
+      _try(KSPSetOperators(_ksp, _a, _a, DIFFERENT_NONZERO_PATTERN));
     _try(MatAssemblyBegin(_a, MAT_FINAL_ASSEMBLY));
     _try(MatAssemblyEnd(_a, MAT_FINAL_ASSEMBLY));
     _try(VecAssemblyBegin(_b));
     _try(VecAssemblyEnd(_b));
-    KSP ksp;
-    _try(KSPCreate(PETSC_COMM_WORLD, &ksp));
-    _try(KSPSetOperators(ksp, _a, _a, DIFFERENT_NONZERO_PATTERN));
-    PC pc;
-    _try(KSPGetPC(ksp, &pc));
-    // set some default options
-    _try(PCSetType(pc, PCILU));
-    _try(PCFactorSetMatOrderingType(pc, MATORDERING_RCM));
-    _try(PCFactorSetLevels(pc, 1));
-    _try(KSPSetTolerances(ksp, 1.e-8, PETSC_DEFAULT, PETSC_DEFAULT, PETSC_DEFAULT));
-    // override the default options with the ones from the option
-    // database (if any)
-    _try(KSPSetFromOptions(ksp));
-    _try(KSPSolve(ksp, _b, _x));
-    _try(KSPView(ksp, PETSC_VIEWER_STDOUT_SELF));
-    PetscInt its;
-    _try(KSPGetIterationNumber(ksp, &its));
-    Msg::Info("%d iterations", its);
+    _try(KSPSolve(_ksp, _b, _x));
+    //_try(KSPView(ksp, PETSC_VIEWER_STDOUT_SELF));
+    //PetscInt its;
+    //_try(KSPGetIterationNumber(ksp, &its));
+    //Msg::Info("%d iterations", its);
     return 1;
   }
   Mat &getMatrix(){ return _a; }