From d39474160cb65b7401ac0e8cda2065359483349c Mon Sep 17 00:00:00 2001
From: Jonathan Lambrechts <jonathan.lambrechts@uclouvain.be>
Date: Tue, 27 Apr 2010 08:14:46 +0000
Subject: [PATCH] block linear solver for dg
---
Geo/GFaceCompound.cpp | 8 ++-
Mesh/BackgroundMesh.cpp | 6 +-
Mesh/highOrderSmoother.cpp | 15 +++--
Numeric/fullMatrix.h | 2 +
Plugin/Distance.cpp | 4 +-
Solver/CMakeLists.txt | 1 +
Solver/dofManager.h | 51 +++++++++++-----
Solver/elasticitySolver.cpp | 6 +-
Solver/function.cpp | 2 +-
Solver/linearSystem.cpp | 12 +---
Solver/linearSystemPETSc.cpp | 110 +++++++++++++++++++++++++++++++++++
Solver/linearSystemPETSc.h | 23 +++++---
Solver/multiscaleLaplace.cpp | 3 +-
13 files changed, 194 insertions(+), 49 deletions(-)
create mode 100644 Solver/linearSystemPETSc.cpp
diff --git a/Geo/GFaceCompound.cpp b/Geo/GFaceCompound.cpp
index 0cbbb498d8..eccdc6bf93 100644
--- a/Geo/GFaceCompound.cpp
+++ b/Geo/GFaceCompound.cpp
@@ -1100,7 +1100,8 @@ void GFaceCompound::parametrize(iterationStep step, typeOfMapping tom) const
for(std::set<MVertex *>::iterator itv = allNodes.begin(); itv !=allNodes.end() ; ++itv){
MVertex *v = *itv;
- double value = myAssembler.getDofValue(v, 0, 1);
+ double value;
+ myAssembler.getDofValue(value, v, 0, 1);
std::map<MVertex*,SPoint3>::iterator itf = coordinates.find(v);
if(itf == coordinates.end()){
SPoint3 p(0, 0, 0);
@@ -1348,8 +1349,9 @@ bool GFaceCompound::parametrize_conformal() const
for(std::set<MVertex *>::iterator itv = allNodes.begin(); itv !=allNodes.end() ; ++itv){
MVertex *v = *itv;
- double value1 = myAssembler.getDofValue(v,0,1);
- double value2 = myAssembler.getDofValue(v,0,2);
+ double value1, value2;
+ myAssembler.getDofValue(value1, v,0,1);
+ myAssembler.getDofValue(value2, v,0,2);
coordinates[v] = SPoint3(value1,value2,0.0);
}
diff --git a/Mesh/BackgroundMesh.cpp b/Mesh/BackgroundMesh.cpp
index 9deb77cb6b..7822b83d7e 100644
--- a/Mesh/BackgroundMesh.cpp
+++ b/Mesh/BackgroundMesh.cpp
@@ -371,8 +371,7 @@ void backgroundMesh::propagate1dMesh(GFace *_gf)
for ( ; itv2 != _2Dto3D.end(); ++itv2){
MVertex *v_2D = itv2->first;
MVertex *v_3D = itv2->second;
- double value = myAssembler.getDofValue(v_3D, 0, 1);
- _sizes[v_2D] = value;
+ myAssembler.getDofValue(_sizes[v_2D], v_3D, 0, 1);
}
delete _lsys;
#endif
@@ -473,8 +472,7 @@ void backgroundMesh::updateSizes(GFace *_gf)
for ( ; itv2 != _2Dto3D.end(); ++itv2){
MVertex *v_2D = itv2->first;
MVertex *v_3D = itv2->second;
- double value = myAssembler.getDofValue(v_3D, 0, 1);
- _sizes[v_2D] = value;
+ myAssembler.getDofValue(_sizes[v_2D], v_3D, 0, 1);
}
delete _lsys;
#endif
diff --git a/Mesh/highOrderSmoother.cpp b/Mesh/highOrderSmoother.cpp
index 95957bee56..281cde696c 100644
--- a/Mesh/highOrderSmoother.cpp
+++ b/Mesh/highOrderSmoother.cpp
@@ -599,8 +599,9 @@ double highOrderSmoother::smooth_metric_(std::vector<MElement*> & v,
if ((*it)->onWhat()->dim() == 2){
SPoint2 param;
reparamMeshVertexOnFace((*it), gf, param);
- SPoint2 dparam (myAssembler.getDofValue((*it), 0, getTag()),
- myAssembler.getDofValue((*it), 1, getTag()));
+ SPoint2 dparam;
+ myAssembler.getDofValue(dparam[0], (*it), 0, getTag());
+ myAssembler.getDofValue(dparam[1], (*it), 1, getTag());
SPoint2 newp = param+dparam;
dx += newp.x() * newp.x() + newp.y() * newp.y();
(*it)->setParameter(0, newp.x());
@@ -715,9 +716,13 @@ void highOrderSmoother::smooth(std::vector<MElement*> &all)
}
for (it = verticesToMove.begin(); it != verticesToMove.end(); ++it){
- it->first->x() += myAssembler.getDofValue(it->first, 0, getTag());
- it->first->y() += myAssembler.getDofValue(it->first, 1, getTag());
- it->first->z() += myAssembler.getDofValue(it->first, 2, getTag());
+ double ax, ay, az;
+ myAssembler.getDofValue(ax, it->first, 0, getTag());
+ myAssembler.getDofValue(ay, it->first, 1, getTag());
+ myAssembler.getDofValue(az, it->first, 2, getTag());
+ it->first->x() += ax;
+ it->first->y() += ay;
+ it->first->z() += az;
}
// delete matrices and vectors
diff --git a/Numeric/fullMatrix.h b/Numeric/fullMatrix.h
index 5d5a5f878e..04630c7ae5 100644
--- a/Numeric/fullMatrix.h
+++ b/Numeric/fullMatrix.h
@@ -63,6 +63,8 @@ class fullVector
{
if(s == 0.)
for(int i = 0; i < _r; ++i) _data[i] = 0.;
+ else if (s == -1.)
+ for(int i = 0; i < _r; ++i) _data[i] = -_data[i];
else
for(int i = 0; i < _r; ++i) _data[i] *= s;
}
diff --git a/Plugin/Distance.cpp b/Plugin/Distance.cpp
index a5cb8df70e..bb56344163 100644
--- a/Plugin/Distance.cpp
+++ b/Plugin/Distance.cpp
@@ -188,7 +188,9 @@ PView *GMSH_DistancePlugin::execute(PView *v)
for(std::map<MVertex*,double >::iterator itv = distance_map2.begin();
itv !=distance_map2.end() ; ++itv){
MVertex *v = itv->first;
- double value = std::min(0.9999, myAssembler.getDofValue(v, 0, 1));
+ double value;
+ myAssembler.getDofValue(value, v, 0, 1);
+ value = std::min(0.9999, value);
double dist = -mu * log(1. - value);
itv->second = dist;
}
diff --git a/Solver/CMakeLists.txt b/Solver/CMakeLists.txt
index fbeea2c22e..d297e54425 100644
--- a/Solver/CMakeLists.txt
+++ b/Solver/CMakeLists.txt
@@ -6,6 +6,7 @@
set(SRC
linearSystem.cpp
linearSystemCSR.cpp
+ linearSystemPETSc.cpp
groupOfElements.cpp
elasticityTerm.cpp
elasticitySolver.cpp
diff --git a/Solver/dofManager.h b/Solver/dofManager.h
index 9df7b57f88..31a1fb782b 100644
--- a/Solver/dofManager.h
+++ b/Solver/dofManager.h
@@ -48,12 +48,18 @@ template<class T> struct dofTraits
{
typedef T VecType;
typedef T MatType;
+ inline static void mult (VecType &r, const MatType &m, const VecType &v) { r = m*v;}
+ inline static void neg (VecType &r) { r = -r;}
+ inline static void setToZero (VecType &r) { r = 0;}
};
-template<class T> struct dofTraits<fullVector<T> >
+template<class T> struct dofTraits<fullMatrix<T> >
{
- typedef fullVector<T> VecType;
+ typedef fullMatrix<T> VecType;
typedef fullMatrix<T> MatType;
+ inline static void mult (VecType &r, const MatType &m, const VecType &v) { m.mult(v, r);}
+ inline static void neg (VecType &r) { r.scale(-1.);}
+ inline static void setToZero (VecType &r) { r.scale(0);}
};
/*
template<> struct dofTraits<fullVector<std::complex<double> > >
@@ -159,24 +165,30 @@ class dofManager{
return dataVec(0.0);
}
- inline dataVec getDofValue(int ent, int type) const
+ inline void getDofValue(dataVec &v, int ent, int type) const
{
Dof key(ent, type);
{
typename std::map<Dof, dataVec>::const_iterator it = fixed.find(key);
- if (it != fixed.end()) return it->second;
+ if (it != fixed.end()){
+ v = it->second;
+ return;
+ }
}
{
std::map<Dof, int>::const_iterator it = unknown.find(key);
- if (it != unknown.end())
- return _current->getFromSolution(it->second);
+ if (it != unknown.end()){
+ v = _current->getFromSolution(it->second);
+ return;
+ }
}
- return dataVec(0.0);
+ dofTraits<T>::setToZero(v);
}
- inline dataVec getDofValue(MVertex *v, int iComp, int iField) const
+ inline dataVec getDofValue(dataVec &value, MVertex *v, int iComp, int iField) const
{
- return getDofValue(v->getNum(), Dof::createTypeWithTwoInts(iComp, iField));
+ getDofValue(value, v->getNum(), Dof::createTypeWithTwoInts(iComp, iField));
}
+
inline void assemble(const Dof &R, const Dof &C, const dataMat &value)
{
if (!_current->isAllocated()) _current->allocate(unknown.size());
@@ -189,8 +201,12 @@ class dofManager{
}
else{
typename std::map<Dof, dataVec>::iterator itFixed = fixed.find(C);
- if (itFixed != fixed.end()){
- _current->addToRightHandSide(itR->second, -value * itFixed->second);
+ if (itFixed != fixed.end()) {
+ // tmp = -value * itFixed->second
+ dataVec tmp(itFixed->second);
+ dofTraits<T>::mult(tmp , value, itFixed->second);
+ dofTraits<T>::neg(tmp);
+ _current->addToRightHandSide(itR->second, tmp);
}
}
}
@@ -220,7 +236,11 @@ class dofManager{
else{
typename std::map<Dof, dataVec>::iterator itFixed = fixed.find(C[j]);
if (itFixed != fixed.end()){
- _current->addToRightHandSide(NR[i], -m(i,j) * itFixed->second);
+ // tmp = -m(i,j) * itFixed->second
+ dataVec tmp(itFixed->second);
+ dofTraits<T>::mult(tmp , m(i,j), itFixed->second);
+ dofTraits<T>::neg(tmp);
+ _current->addToRightHandSide(NR[i], tmp);
}
}
}
@@ -272,7 +292,11 @@ class dofManager{
{
typename std::map<Dof, dataVec>::iterator itFixed = fixed.find(R[j]);
if (itFixed != fixed.end()){
- _current->addToRightHandSide(NR[i], -m(i, j) * itFixed->second);
+ // tmp = -m(i,j) * itFixed->second
+ dataVec tmp(itFixed->second);
+ dofTraits<T>::mult(tmp , m(i,j), itFixed->second);
+ dofTraits<T>::neg(tmp);
+ _current->addToRightHandSide(NR[i], tmp);
}
}
}
@@ -336,5 +360,4 @@ class dofManager{
return 0;
}
};
-
#endif
diff --git a/Solver/elasticitySolver.cpp b/Solver/elasticitySolver.cpp
index 047e83fdbc..65540b388c 100644
--- a/Solver/elasticitySolver.cpp
+++ b/Solver/elasticitySolver.cpp
@@ -225,9 +225,9 @@ static double vonMises(dofManager<double> *a, MElement *e,
double valy[256];
double valz[256];
for (int k = 0; k < e->getNumVertices(); k++){
- valx[k] = a->getDofValue(e->getVertex(k), 0, _tag);
- valy[k] = a->getDofValue(e->getVertex(k), 1, _tag);
- valz[k] = a->getDofValue(e->getVertex(k), 2, _tag);
+ a->getDofValue(valx[k], e->getVertex(k), 0, _tag);
+ a->getDofValue(valy[k], e->getVertex(k), 1, _tag);
+ a->getDofValue(valz[k], e->getVertex(k), 2, _tag);
}
double gradux[3];
double graduy[3];
diff --git a/Solver/function.cpp b/Solver/function.cpp
index aadd3c0699..2d65bc040e 100644
--- a/Solver/function.cpp
+++ b/Solver/function.cpp
@@ -447,7 +447,7 @@ void function::registerBindings(binding *b){
mb->setDescription("return the unique instance of this class");
cb->setParentClass<function>();
- cb = b->addClass<functionSolution>("functionNormals");
+ cb = b->addClass<functionNormals>("functionNormals");
cb->setDescription("A function to access the face normals. This is a single-instance class, use the 'get' member to access the instance.");
mb = cb->addMethod("get",&functionNormals::get);
mb->setDescription("return the unique instance of this class");
diff --git a/Solver/linearSystem.cpp b/Solver/linearSystem.cpp
index 32d0022012..84193cd788 100644
--- a/Solver/linearSystem.cpp
+++ b/Solver/linearSystem.cpp
@@ -28,19 +28,13 @@ void linearSystem<double>::registerBindings(binding *b){
cb->setParentClass<linearSystem<double> >();
#endif
-#ifdef HAVE_PETSC
- cb = b->addClass<linearSystemPETSc<double> >("linearSystemPETSc");
- cb->setDescription("A linear system solver, based on PETSc");
- cm = cb->setConstructor<linearSystemPETSc<double> >();
- cm->setDescription ("A new PETSc<double> solver");
- cm->setArgNames(NULL);
- cb->setParentClass<linearSystem<double> >();
-#endif
-
cb = b->addClass<linearSystemFull<double> >("linearSystemFull");
cb->setDescription("A linear system solver, based on LAPACK (full matrices)");
cm = cb->setConstructor<linearSystemFull<double> >();
cm->setDescription ("A new Lapack based <double> solver");
cm->setArgNames(NULL);
cb->setParentClass<linearSystem<double> >();
+#ifdef HAVE_PETSC
+ linearSystemPETScRegisterBindings (b);
+#endif
}
diff --git a/Solver/linearSystemPETSc.cpp b/Solver/linearSystemPETSc.cpp
new file mode 100644
index 0000000000..c77af8d464
--- /dev/null
+++ b/Solver/linearSystemPETSc.cpp
@@ -0,0 +1,110 @@
+#include "GmshConfig.h"
+#if defined HAVE_PETSC
+#include "linearSystemPETSc.h"
+#include "fullMatrix.h"
+template <>
+void linearSystemPETSc<fullMatrix<double> >::addToMatrix(int row, int col, fullMatrix<double> val)
+{
+ for (int ii = 0; ii<val.size1(); ii++)
+ for (int jj = 0; jj < ii; jj ++) {
+ double buff = val(ii,jj);
+ val(ii,jj) = val (jj,ii);
+ val(jj,ii) = buff;
+ }
+ PetscInt i = row, j = col;
+ _try(MatSetValuesBlocked(_a, 1, &i, 1, &j, &val(0,0), ADD_VALUES));
+}
+template<>
+void linearSystemPETSc<fullMatrix<double> >::addToRightHandSide(int row, fullMatrix<double> val)
+{
+ for (int ii = 0; ii < _blockSize; ii ++) {
+ PetscInt i = row * _blockSize + ii;
+ _try(VecSetValues(_b, 1, &i, &val(ii,0), ADD_VALUES));
+ }
+}
+
+template<>
+fullMatrix<double> linearSystemPETSc<fullMatrix<double> >::getFromMatrix(int row, int col) const
+{
+ Msg::Error("getFromMatrix not implemented for PETSc");
+ return fullMatrix<double>(0,0);
+}
+
+template<>
+fullMatrix<double> linearSystemPETSc<fullMatrix<double> >::getFromRightHandSide(int row) const
+{
+ fullMatrix<double> val(_blockSize,1);
+ PetscScalar *tmp;
+ _try(VecGetArray(_b, &tmp));
+ for (int i=0; i<_blockSize; i++) {
+ PetscScalar s = tmp[row*_blockSize+i];
+#if defined(PETSC_USE_COMPLEX)
+ val(i,0) = s.real();
+#else
+ val(i,0) = s;
+#endif
+ }
+ _try(VecRestoreArray(_b, &tmp));
+ return val;
+}
+
+template<>
+fullMatrix<double> linearSystemPETSc<fullMatrix<double> >::getFromSolution(int row) const
+{
+ fullMatrix<double> val(_blockSize,1);
+ PetscScalar *tmp;
+ _try(VecGetArray(_x, &tmp));
+ for (int i=0; i<_blockSize; i++) {
+ PetscScalar s = tmp[row*_blockSize+i];
+#if defined(PETSC_USE_COMPLEX)
+ val(i,0) = s.real();
+#else
+ val(i,0) = s;
+#endif
+ }
+ _try(VecRestoreArray(_x, &tmp));
+ return val;
+}
+
+template<>
+void linearSystemPETSc<fullMatrix<double> >::allocate(int nbRows) {
+ clear();
+ _try(MatCreate(PETSC_COMM_WORLD, &_a));
+ _try(MatSetSizes(_a, PETSC_DECIDE, PETSC_DECIDE, nbRows*_blockSize, nbRows*_blockSize));
+ _try(MatSetType(_a, MATSEQBAIJ));
+ // override the default options with the ones from the option
+ // database (if any)
+ _try(MatSetFromOptions(_a));
+ _try(MatSeqBAIJSetPreallocation(_a, _blockSize, 4, NULL)); //todo preAllocate off-diagonal part
+ //_try(MatMPIBAIJSetPreallocation(_a, _blockSize, 4, NULL, 0, NULL)); //todo preAllocate off-diagonal part
+ _try(VecCreate(PETSC_COMM_WORLD, &_x));
+ _try(VecSetSizes(_x, PETSC_DECIDE, nbRows*_blockSize));
+ // override the default options with the ones from the option
+ // database (if any)
+ _try(VecSetFromOptions(_x));
+ _try(VecDuplicate(_x, &_b));
+ _isAllocated = true;
+}
+
+#include "Bindings.h"
+void linearSystemPETScRegisterBindings(binding *b) {
+ classBinding *cb;
+ methodBinding *cm;
+ cb = b->addClass<linearSystemPETSc<double> >("linearSystemPETSc");
+ cb->setDescription("A linear system solver, based on PETSc");
+ cm = cb->setConstructor<linearSystemPETSc<double> >();
+ cm->setDescription ("A new PETSc<double> solver");
+ cb->setParentClass<linearSystem<double> >();
+ cm->setArgNames(NULL);
+ cm = cb->addMethod("systemSolve", &linearSystem<fullMatrix<double> >::systemSolve);
+ cm->setDescription("compute x = A^{-1}b");
+
+ cb = b->addClass<linearSystemPETSc<fullMatrix<double> > >("linearSystemPETScBlock");
+ cb->setDescription("A linear system solver, based on PETSc");
+ cm = cb->setConstructor<linearSystemPETSc<fullMatrix<double> >, int>();
+ cm->setDescription ("A new PETScBlock<double> solver (we probably should get rid of the blockSize argument)");
+ cm->setArgNames("blockSize", NULL);
+ cm = cb->addMethod("systemSolve", &linearSystem<fullMatrix<double> >::systemSolve);
+ cm->setDescription("compute x = A^{-1}b");
+}
+#endif // HAVE_PETSC
diff --git a/Solver/linearSystemPETSc.h b/Solver/linearSystemPETSc.h
index acc4a040e0..d17a14c333 100644
--- a/Solver/linearSystemPETSc.h
+++ b/Solver/linearSystemPETSc.h
@@ -44,13 +44,13 @@
template <class scalar>
class linearSystemPETSc : public linearSystem<scalar> {
- private:
+ int _blockSize; // for block Matrix
bool _isAllocated;
Mat _a;
Vec _b, _x;
void _try(int ierr) const { CHKERRABORT(PETSC_COMM_WORLD, ierr); }
public:
- linearSystemPETSc() : _isAllocated(false) {}
+ linearSystemPETSc(int blockSize = -1) : _isAllocated(false) {_blockSize = blockSize;}
virtual ~linearSystemPETSc(){ clear(); }
virtual bool isAllocated() const { return _isAllocated; }
virtual void allocate(int nbRows)
@@ -83,12 +83,6 @@ class linearSystemPETSc : public linearSystem<scalar> {
}
_isAllocated = false;
}
- virtual void addToMatrix(int row, int col, scalar val)
- {
- PetscInt i = row, j = col;
- PetscScalar s = val;
- _try(MatSetValues(_a, 1, &i, 1, &j, &s, ADD_VALUES));
- }
virtual scalar getFromMatrix(int row, int col) const
{
Msg::Error("getFromMatrix not implemented for PETSc");
@@ -113,6 +107,12 @@ class linearSystemPETSc : public linearSystem<scalar> {
return s;
#endif
}
+ virtual void addToMatrix(int row, int col, scalar val)
+ {
+ PetscInt i = row, j = col;
+ PetscScalar s = val;
+ _try(MatSetValues(_a, 1, &i, 1, &j, &s, ADD_VALUES));
+ }
virtual scalar getFromSolution(int row) const
{
PetscScalar *tmp;
@@ -152,6 +152,9 @@ class linearSystemPETSc : public linearSystem<scalar> {
// override the default options with the ones from the option
// database (if any)
_try(KSPSetFromOptions(ksp));
+ PetscViewer view;
+ PetscViewerASCIIOpen(PETSC_COMM_WORLD,"mat.out", &view);
+ _try(MatView(_a, view));
_try(KSPSolve(ksp, _b, _x));
_try(KSPView(ksp, PETSC_VIEWER_STDOUT_SELF));
PetscInt its;
@@ -161,6 +164,10 @@ class linearSystemPETSc : public linearSystem<scalar> {
Mat &getMatrix(){ return _a; }
};
+class binding;
+void linearSystemPETScRegisterBindings(binding *b);
+
+
#else
template <class scalar>
diff --git a/Solver/multiscaleLaplace.cpp b/Solver/multiscaleLaplace.cpp
index 73637c6874..29acbb5a29 100644
--- a/Solver/multiscaleLaplace.cpp
+++ b/Solver/multiscaleLaplace.cpp
@@ -1199,7 +1199,8 @@ void multiscaleLaplace::parametrize_method (multiscaleLaplaceLevel & level,
int count = 0;
for(std::set<MVertex *>::iterator itv = allNodes.begin(); itv !=allNodes.end() ; ++itv){
MVertex *v = *itv;
- double value = myAssembler.getDofValue(v, 0, 1);
+ double value;
+ myAssembler.getDofValue(value, v, 0, 1);
if (step == 0)solution[v] = SPoint2(value,0);
else solution[v] = SPoint2(solution[v][0],value);
}
--
GitLab