From 0c44894171fe78be70a50647895109182470be20 Mon Sep 17 00:00:00 2001
From: Jonathan Lambrechts <jonathan.lambrechts@uclouvain.be>
Date: Fri, 15 Jan 2010 12:41:32 +0000
Subject: [PATCH] dg : first working parallel version
---
Numeric/fullMatrix.h | 4 +
Solver/CMakeLists.txt | 1 +
Solver/dgAlgorithm.cpp | 13 ++-
Solver/dgAlgorithm.h | 1 -
Solver/dgDofContainer.cpp | 139 +++++++++++++++++++++++++
Solver/dgDofContainer.h | 34 +++++++
Solver/dgGroupOfElements.cpp | 71 ++++++++++++-
Solver/dgGroupOfElements.h | 8 ++
Solver/dgSystemOfEquations.cpp | 181 +--------------------------------
Solver/dgSystemOfEquations.h | 31 +-----
10 files changed, 266 insertions(+), 217 deletions(-)
create mode 100644 Solver/dgDofContainer.cpp
create mode 100644 Solver/dgDofContainer.h
diff --git a/Numeric/fullMatrix.h b/Numeric/fullMatrix.h
index dcbb2416a9..dcf652c6ca 100644
--- a/Numeric/fullMatrix.h
+++ b/Numeric/fullMatrix.h
@@ -224,6 +224,10 @@ class fullMatrix
{
for(int i = 0; i < _r * _c; i++) _data[i] = m;
}
+ inline void setAll(const fullMatrix<scalar> &m)
+ {
+ for(int i = 0; i < _r * _c; i++) _data[i] = m._data[i];
+ }
inline void scale(const double s)
{
if(s == 0.)
diff --git a/Solver/CMakeLists.txt b/Solver/CMakeLists.txt
index aaea186caf..3392ddca02 100644
--- a/Solver/CMakeLists.txt
+++ b/Solver/CMakeLists.txt
@@ -20,6 +20,7 @@ set(SRC
dgConservationLawWaveEquation.cpp
dgConservationLawPerfectGas.cpp
dgSlopeLimiter.cpp
+ dgDofContainer.cpp
function.cpp
functionDerivator.cpp
luaFunction.cpp
diff --git a/Solver/dgAlgorithm.cpp b/Solver/dgAlgorithm.cpp
index 498f37b4bb..462ea79ade 100644
--- a/Solver/dgAlgorithm.cpp
+++ b/Solver/dgAlgorithm.cpp
@@ -259,7 +259,6 @@ void dgAlgorithm::rungeKutta (const dgConservationLaw &claw, // conservation l
double h, // time-step
dgDofContainer &sol,
dgDofContainer &resd,
- dgSystemOfEquations *syst,
dgLimiter *limiter,
int orderRK
) // order of RK integrator
@@ -281,8 +280,8 @@ void dgAlgorithm::rungeKutta (const dgConservationLaw &claw, // conservation l
int nbFields = claw.nbFields();
- dgDofContainer K (groups,claw);
- dgDofContainer Unp (groups,claw);
+ dgDofContainer K (groups,nbFields);
+ dgDofContainer Unp (groups,nbFields);
K._data->scale(0.0);
K._data->axpy(*(sol._data));
@@ -298,7 +297,7 @@ void dgAlgorithm::rungeKutta (const dgConservationLaw &claw, // conservation l
if (limiter){
limiter->apply(K, groups);
}
- syst->scatter(&K);
+ K.scatter();
this->residual(claw,groups,K._dataProxys,resd._dataProxys);
K._data->scale(0.);
for(int k=0; k < groups.getNbElementGroups(); k++) {
@@ -408,11 +407,11 @@ void dgAlgorithm::multirateRungeKutta (const dgConservationLaw &claw, // conse
dgDofContainer **K;
K=new dgDofContainer*[nStages];
for(int i=0;i<nStages;i++){
- K[i]=new dgDofContainer(groups,claw);
+ K[i]=new dgDofContainer(groups,nbFields);
K[i]->_data->scale(0.0);
}
- dgDofContainer Unp (groups,claw);
- dgDofContainer tmp (groups,claw);
+ dgDofContainer Unp (groups,nbFields);
+ dgDofContainer tmp (groups,nbFields);
Unp._data->scale(0.0);
Unp._data->axpy(*(sol._data));
diff --git a/Solver/dgAlgorithm.h b/Solver/dgAlgorithm.h
index c3dd4aa656..206bfcfd75 100644
--- a/Solver/dgAlgorithm.h
+++ b/Solver/dgAlgorithm.h
@@ -45,7 +45,6 @@ class dgAlgorithm {
double h,
dgDofContainer &residu,
dgDofContainer &sol,
- dgSystemOfEquations *syst,
dgLimiter *limiter=NULL,
int orderRK=4);
static void computeElementaryTimeSteps ( //dofManager &dof, // the DOF manager (maybe useless here)
diff --git a/Solver/dgDofContainer.cpp b/Solver/dgDofContainer.cpp
new file mode 100644
index 0000000000..9d20e7ba14
--- /dev/null
+++ b/Solver/dgDofContainer.cpp
@@ -0,0 +1,139 @@
+#include "GmshConfig.h"
+#include "dgDofContainer.h"
+#include "dgGroupOfElements.h"
+#ifdef HAVE_MPI
+#include "mpi.h"
+#endif
+dgDofContainer::dgDofContainer (dgGroupCollection &groups, int nbFields):
+ _groups(groups)
+{
+ _dataSize = 0;
+ _dataSizeGhost=0;
+ _totalNbElements = 0;
+ _parallelStructureExists = false;
+ _nbFields = nbFields;
+ for (int i=0; i< _groups.getNbElementGroups();i++){
+ int nbNodes = _groups.getElementGroup(i)->getNbNodes();
+ int nbElements = _groups.getElementGroup(i)->getNbElements();
+ _totalNbElements +=nbElements;
+ _dataSize += nbNodes*_nbFields*nbElements;
+ }
+
+ // allocate the big vectors
+ _data = new fullVector<double>(_dataSize);
+ // create proxys for each group
+ int offset = 0;
+ _dataProxys.resize(_groups.getNbElementGroups()+_groups.getNbGhostGroups());
+ for (int i=0;i<_groups.getNbElementGroups();i++){
+ dgGroupOfElements *group = _groups.getElementGroup(i);
+ int nbNodes = group->getNbNodes();
+ int nbElements = group->getNbElements();
+ _dataProxys[i] = new fullMatrix<double> (&(*_data)(offset),nbNodes, _nbFields*nbElements);
+ offset += nbNodes*_nbFields*nbElements;
+ }
+
+ //ghosts
+
+ int totalNbElementsGhost =0;
+ for (int i=0; i<_groups.getNbGhostGroups(); i++){
+ int nbNodes = _groups.getGhostGroup(i)->getNbNodes();
+ int nbElements = _groups.getGhostGroup(i)->getNbElements();
+ totalNbElementsGhost +=nbElements;
+ _dataSizeGhost += nbNodes*_nbFields*nbElements;
+ }
+
+ _dataProxys.resize(_groups.getNbElementGroups()+_groups.getNbGhostGroups());
+ _ghostData = new fullVector<double>(_dataSizeGhost);
+ offset=0;
+ for (int i=0;i<_groups.getNbGhostGroups();i++){
+ dgGroupOfElements *group = _groups.getGhostGroup(i);
+ int nbNodes = group->getNbNodes();
+ int nbElements = group->getNbElements();
+ _dataProxys[i+_groups.getNbElementGroups()] = new fullMatrix<double> (&(*_ghostData)(offset),nbNodes, _nbFields*nbElements);
+ offset += nbNodes*_nbFields*nbElements;
+ }
+
+}
+
+dgDofContainer::~dgDofContainer (){
+ if(_parallelStructureExists) {
+ delete []countSend;
+ delete []countRecv;
+ delete []shiftSend;
+ delete []shiftRecv;
+ delete []groupShiftRecv;
+ delete []sendBuf;
+ delete []recvBuf;
+ }
+ if (!_dataSize)return;
+ for (int i=0;i<_dataProxys.size();++i) delete _dataProxys[i];
+ delete _data;
+ delete _ghostData;
+}
+
+void dgDofContainer::buildParallelStructure(){
+ if (_parallelStructureExists)
+ return;
+
+ // MPI all2all buffers
+ countSend = new int[Msg::GetCommSize()];
+ shiftSend = new int[Msg::GetCommSize()];
+ countRecv = new int[Msg::GetCommSize()];
+ shiftRecv = new int[Msg::GetCommSize()];
+ groupShiftRecv = new int[_groups.getNbGhostGroups()];
+ for(int i=0;i<Msg::GetCommSize();i++){
+ countSend[i]=0;
+ countRecv[i]=0;
+ for(size_t j=0;j<_groups.getNbImageElementsOnPartition(i);j++){
+ countSend[i] += _groups.getElementGroup(_groups.getImageElementGroup(i,j))->getNbNodes()*_nbFields;
+ }
+ }
+ for(size_t i=0; i<_groups.getNbGhostGroups(); i++){
+ dgGroupOfElements *group = _groups.getGhostGroup(i);
+ groupShiftRecv[i] = countRecv[group->getGhostPartition()];
+ countRecv[group->getGhostPartition()]+=group->getNbElements()*group->getNbNodes()*_nbFields;
+ }
+ shiftSend[0] = shiftRecv[0]=0;
+ int totalSend = countSend[0];
+ int totalRecv = countRecv[0];
+ for(size_t i=1; i<Msg::GetCommSize(); i++){
+ shiftSend[i] = shiftSend[i-1]+countSend[i-1];
+ shiftRecv[i] = shiftRecv[i-1]+countRecv[i-1];
+ totalSend += countSend[i];
+ totalRecv += countRecv[i];
+ }
+ for(size_t i=0; i<_groups.getNbGhostGroups(); i++){
+ dgGroupOfElements *group = _groups.getGhostGroup(i);
+ groupShiftRecv[i] += shiftRecv[group->getGhostPartition()];
+ }
+ sendBuf = new double[totalSend];
+ recvBuf = new double[totalRecv];
+
+ _parallelStructureExists = true;
+}
+
+void dgDofContainer::scatter() {
+ if(!_parallelStructureExists)
+ buildParallelStructure();
+ //1) fill
+ int index=0;
+ for(int i=0;i<Msg::GetCommSize();i++){
+ for(size_t j=0;j<_groups.getNbImageElementsOnPartition(i);j++){
+ fullMatrix<double> &sol = getGroupProxy(_groups.getImageElementGroup(i,j));
+ int eid = _groups.getImageElementPositionInGroup(i,j);
+ for(int l=0;l<_nbFields;l++) {
+ for(int k=0;k<sol.size1();k++) {
+ sendBuf[index++] = sol(k, eid *_nbFields+l);
+ }
+ }
+ }
+ }
+ //2) send
+ MPI_Alltoallv(sendBuf,countSend,shiftSend,MPI_DOUBLE,recvBuf,countRecv,shiftRecv,MPI_DOUBLE,MPI_COMM_WORLD);
+ //3) distribute
+ for(int i=0; i< _groups.getNbGhostGroups();i++) {
+ fullMatrix<double> &sol = getGroupProxy(i+_groups.getNbElementGroups());
+ fullMatrix<double> recvProxy (recvBuf + groupShiftRecv[i], sol.size1(), sol.size2());
+ sol.setAll(recvProxy);
+ }
+}
diff --git a/Solver/dgDofContainer.h b/Solver/dgDofContainer.h
new file mode 100644
index 0000000000..7d4298507d
--- /dev/null
+++ b/Solver/dgDofContainer.h
@@ -0,0 +1,34 @@
+#ifndef _DG_DOF_CONTAINER_H_
+#define _DG_DOF_CONTAINER_H_
+
+#include <vector>
+#include "fullMatrix.h"
+class dgGroupCollection;
+
+struct dgDofContainer {
+private:
+ dgDofContainer (const dgDofContainer&);
+ int _totalNbElements;
+ int _nbFields;
+ dgGroupCollection &_groups;
+ int _dataSizeGhost;
+ fullVector<double> * _ghostData;
+ inline int getDataSize(){return _dataSize;}
+ // parallel
+ void buildParallelStructure();
+ bool _parallelStructureExists;
+ int *countSend,*countRecv,*shiftSend,*shiftRecv,*groupShiftRecv;
+ double *sendBuf, *recvBuf;
+public:
+ //todo move those 3 to private
+ fullVector<double> * _data; // the full data itself
+ std::vector<fullMatrix<double> *> _dataProxys; // proxys
+ int _dataSize; // the full data size i.e. concerning all groups (not ghost, see bellow)
+
+ inline fullMatrix<double> &getGroupProxy(int gId){ return *(_dataProxys[gId]); }
+ dgDofContainer (dgGroupCollection &groups, int nbFields);
+ ~dgDofContainer ();
+ int getNbElements() {return _totalNbElements;}
+ void scatter();
+};
+#endif
diff --git a/Solver/dgGroupOfElements.cpp b/Solver/dgGroupOfElements.cpp
index a1867f108e..16c6ac338b 100644
--- a/Solver/dgGroupOfElements.cpp
+++ b/Solver/dgGroupOfElements.cpp
@@ -1,3 +1,4 @@
+#include "GmshConfig.h"
#include "dgGroupOfElements.h"
#include "MElement.h"
#include "polynomialBasis.h"
@@ -7,6 +8,9 @@
#include "MLine.h"
#include "MPoint.h"
#include "GModel.h"
+#ifdef HAVE_MPI
+#include "mpi.h"
+#endif
static fullMatrix<double> * dgGetIntegrationRule (MElement *e, int p){
int npts;
@@ -515,7 +519,6 @@ dgGroupOfFaces::dgGroupOfFaces (const dgGroupOfElements &elGroup1, const dgGroup
}
}
- printf("groupe mixte %d %d elements %d faces\n",elGroup1.getNbElements(),elGroup2.getNbElements(),_faces.size());
break;
}
case 3 : {
@@ -681,8 +684,10 @@ void dgGroupCollection::buildGroups(GModel *model, int dim, int order)
}else{
std::multimap<MElement*, short>::iterator ghost=ghostsCells.lower_bound(el);
while(ghost!= ghostsCells.end() && ghost->first==el){
- nghosts+=(abs(ghost->second)-1==Msg::GetCommRank());
- ghostElements[el->getPartition()-1][el->getType()].push_back(el);
+ if(abs(ghost->second)-1==Msg::GetCommRank()){
+ ghostElements[el->getPartition()-1][el->getType()].push_back(el);
+ nghosts+=1;
+ }
ghost++;
}
}
@@ -750,7 +755,65 @@ void dgGroupCollection::buildGroups(GModel *model, int dim, int order)
delete gof;
}
}
- printf("nbGhostGroups : %i nbElementsGroups : %i\n",getNbGhostGroups(), getNbElementGroups());
+
+ // build the ghosts structure
+ int *nGhostElements = new int[Msg::GetCommSize()];
+ int *nParentElements = new int[Msg::GetCommSize()];
+ for (int i=0;i<Msg::GetCommSize();i++) {
+ nGhostElements[i]=0;
+ }
+ for (size_t i = 0; i< getNbGhostGroups(); i++){
+ dgGroupOfElements *group = getGhostGroup(i);
+ nGhostElements[group->getGhostPartition()] += group->getNbElements();
+ }
+ #ifdef HAVE_MPI
+ MPI_Alltoall(nGhostElements,1,MPI_INT,nParentElements,1,MPI_INT,MPI_COMM_WORLD);
+ #else
+ nParentElements[0]=nGhostElements[0];
+ #endif
+ int *shiftSend = new int[Msg::GetCommSize()];
+ int *shiftRecv = new int[Msg::GetCommSize()];
+ int *curShiftSend = new int[Msg::GetCommSize()];
+ int totalSend=0,totalRecv=0;
+ for(int i= 0; i<Msg::GetCommSize();i++){
+ shiftSend[i] = (i==0 ? 0 : shiftSend[i-1]+nGhostElements[i-1]);
+ curShiftSend[i] = shiftSend[i];
+ shiftRecv[i] = (i==0 ? 0 : shiftRecv[i-1]+nParentElements[i-1]);
+ totalSend += nGhostElements[i];
+ totalRecv += nParentElements[i];
+ }
+
+ int *idSend = new int[totalSend];
+ int *idRecv = new int[totalRecv];
+ for (int i = 0; i<getNbGhostGroups(); i++){
+ dgGroupOfElements *group = getGhostGroup(i);
+ int part = group->getGhostPartition();
+ for (int j=0; j< group->getNbElements() ; j++){
+ idSend[curShiftSend[part]++] = group->getElement(j)->getNum();
+ }
+ }
+ MPI_Alltoallv(idSend,nGhostElements,shiftSend,MPI_INT,idRecv,nParentElements,shiftRecv,MPI_INT,MPI_COMM_WORLD);
+ //create a Map elementNum :: group, position in group
+ std::map<int, std::pair<int,int> > elementMap;
+ for(size_t i = 0; i< getNbElementGroups(); i++) {
+ dgGroupOfElements *group = getElementGroup(i);
+ for(int j=0; j< group->getNbElements(); j++){
+ elementMap[group->getElement(j)->getNum()]=std::pair<int,int>(i,j);
+ }
+ }
+ _elementsToSend.resize(Msg::GetCommSize());
+ for(int i = 0; i<Msg::GetCommSize();i++){
+ _elementsToSend[i].resize(nParentElements[i]);
+ for(int j=0;j<nParentElements[i];j++){
+ int num = idRecv[shiftRecv[i]+j];
+ _elementsToSend[i][j]=elementMap[num];
+ }
+ }
+ delete []idRecv;
+ delete []idSend;
+ delete []curShiftSend;
+ delete []shiftSend;
+ delete []shiftRecv;
}
dgGroupCollection::~dgGroupCollection() {
diff --git a/Solver/dgGroupOfElements.h b/Solver/dgGroupOfElements.h
index 7c0f177772..55f742daa4 100644
--- a/Solver/dgGroupOfElements.h
+++ b/Solver/dgGroupOfElements.h
@@ -188,6 +188,10 @@ class dgGroupCollection {
std::vector<dgGroupOfFaces*> _faceGroups; //interface
std::vector<dgGroupOfFaces*> _boundaryGroups; //boundary
std::vector<dgGroupOfElements*> _ghostGroups; //ghost volume
+
+ //{group,id} of the elements to send to each partition for a scatter operation
+ std::vector< std::vector<std::pair<int,int> > >_elementsToSend;
+
public:
inline int getNbElementGroups() const {return _elementGroups.size();}
inline int getNbFaceGroups() const {return _faceGroups.size();}
@@ -198,6 +202,10 @@ class dgGroupCollection {
inline dgGroupOfFaces *getBoundaryGroup(int i) const {return _boundaryGroups[i];}
inline dgGroupOfElements *getGhostGroup(int i) const {return _ghostGroups[i];}
+ inline int getNbImageElementsOnPartition(int partId) const {return _elementsToSend[partId].size();}
+ inline int getImageElementGroup(int partId, int i) const {return _elementsToSend[partId][i].first;}
+ inline int getImageElementPositionInGroup(int partId, int i) const {return _elementsToSend[partId][i].second;}
+
void buildGroups (GModel *model,int dimension, int order);
~dgGroupCollection();
};
diff --git a/Solver/dgSystemOfEquations.cpp b/Solver/dgSystemOfEquations.cpp
index ddd220d488..405537ccc1 100644
--- a/Solver/dgSystemOfEquations.cpp
+++ b/Solver/dgSystemOfEquations.cpp
@@ -100,68 +100,13 @@ void dgSystemOfEquations::L2Projection (std::string functionName){
void dgSystemOfEquations::setup(){
if (!_claw) throw;
_groups.buildGroups(_gm,_dimension,_order);
- // build the ghosts structure
- int *nGhostElements = new int[Msg::GetCommSize()];
- int *nParentElements = new int[Msg::GetCommSize()];
- for (int i=0;i<Msg::GetCommSize();i++) {
- nGhostElements[i]=0;
- }
- for (size_t i = 0; i< _groups.getNbGhostGroups(); i++){
- dgGroupOfElements *group=_groups.getGhostGroup(i);
- nGhostElements[group->getGhostPartition()] += group->getNbElements();
- }
- #ifdef HAVE_MPI
- MPI_Alltoall(nGhostElements,1,MPI_INT,nParentElements,1,MPI_INT,MPI_COMM_WORLD);
- #else
- nParentElements[0]=nGhostElements[0];
- #endif
- shiftSend = new int[Msg::GetCommSize()];
- shiftRecv = new int[Msg::GetCommSize()];
- int *curShiftSend = new int[Msg::GetCommSize()];
- totalSend=0,totalRecv=0;
- for(int i= 0; i<Msg::GetCommSize();i++){
- shiftSend[i] = (i==0 ? 0 : shiftSend[i-1]+nGhostElements[i-1]);
- curShiftSend[i] = shiftSend[i];
- shiftRecv[i] = (i==0 ? 0 : shiftRecv[i-1]+nParentElements[i-1]);
- totalSend+=nGhostElements[i];
- totalRecv+=nParentElements[i];
- }
-
- int *idSend = new int[totalSend];
- int *idRecv = new int[totalRecv];
- for (int i = 0; i< _groups.getNbGhostGroups(); i++){
- dgGroupOfElements *group = _groups.getGhostGroup(i);
- int part = group->getGhostPartition();
- for (int j=0; j< group->getNbElements() ; j++){
- idSend[curShiftSend[part]++] = group->getElement(j)->getNum();
- }
- }
- MPI_Alltoallv(idSend,nGhostElements,shiftSend,MPI_INT,idRecv,nParentElements,shiftRecv,MPI_INT,MPI_COMM_WORLD);
- //create a Map elementNum :: group, position in group
- std::map<int, std::pair<int,int> > elementMap;
- for(size_t i = 0; i< _groups.getNbElementGroups(); i++) {
- dgGroupOfElements *group = _groups.getElementGroup(i);
- for(int j=0; j< group->getNbElements(); j++){
- elementMap[group->getElement(j)->getNum()]=std::pair<int,int>(i,j);
- }
- }
- _elementsToSend.resize(Msg::GetCommSize());
- for(int i = 0; i<Msg::GetCommSize();i++){
- _elementsToSend[i].resize(nParentElements[i]);
- for(int j=0;j<nParentElements[i];j++){
- int num = idRecv[shiftRecv[i]+j];
- _elementsToSend[i][j]=elementMap[num];
- }
- }
- // compute the total size of the soution
- delete curShiftSend;
- _solution = new dgDofContainer(_groups,*_claw);
- _rightHandSide = new dgDofContainer(_groups,*_claw);
+ _solution = new dgDofContainer(_groups,_claw->nbFields());
+ _rightHandSide = new dgDofContainer(_groups,_claw->nbFields());
}
double dgSystemOfEquations::RK44(double dt){
- _algo->rungeKutta(*_claw,_groups, dt, *_solution, *_rightHandSide, this);
+ _algo->rungeKutta(*_claw,_groups, dt, *_solution, *_rightHandSide);
return _solution->_data->norm();
}
@@ -184,13 +129,13 @@ double dgSystemOfEquations::computeInvSpectralRadius(){
double dgSystemOfEquations::RK44_limiter(double dt){
dgLimiter *sl = new dgSlopeLimiter(_claw);
- _algo->rungeKutta(*_claw,_groups, dt, *_solution, *_rightHandSide,this, sl);
+ _algo->rungeKutta(*_claw,_groups, dt, *_solution, *_rightHandSide, sl);
delete sl;
return _solution->_data->norm();
}
double dgSystemOfEquations::ForwardEuler(double dt){
- _algo->rungeKutta(*_claw, _groups, dt, *_solution, *_rightHandSide,this, NULL,1);
+ _algo->rungeKutta(*_claw, _groups, dt, *_solution, *_rightHandSide, NULL,1);
return _solution->_data->norm();
}
double dgSystemOfEquations::multirateRK43(double dt){
@@ -293,121 +238,5 @@ void dgSystemOfEquations::export_solution_as_is (const std::string &name){
delete pv;
}
-dgDofContainer::dgDofContainer (dgGroupCollection &groups, const dgConservationLaw &claw):
- _groups(groups)
-{
- _dataSize = 0;
- _dataSizeGhost=0;
- totalNbElements = 0;
- int totalNbElementsGhost =0;
- nbFields = claw.nbFields();
- for (int i=0; i< _groups.getNbElementGroups();i++){
- int nbNodes = _groups.getElementGroup(i)->getNbNodes();
- int nbElements = _groups.getElementGroup(i)->getNbElements();
- totalNbElements +=nbElements;
- _dataSize += nbNodes*nbFields*nbElements;
- }
- for (int i=0; i<_groups.getNbGhostGroups(); i++){
- int nbNodes = _groups.getGhostGroup(i)->getNbNodes();
- int nbElements = _groups.getGhostGroup(i)->getNbElements();
- totalNbElementsGhost +=nbElements;
- _dataSizeGhost += nbNodes*nbFields*nbElements;
- }
-
- // allocate the big vectors
- _data = new fullVector<double>(_dataSize);
- _ghostData = new fullVector<double>(_dataSizeGhost);
- // create proxys for each group
- int offset = 0;
- _dataProxys.resize(_groups.getNbElementGroups()+_groups.getNbGhostGroups());
- for (int i=0;i<_groups.getNbElementGroups();i++){
- dgGroupOfElements *group = _groups.getElementGroup(i);
- int nbNodes = group->getNbNodes();
- int nbElements = group->getNbElements();
- _dataProxys[i] = new fullMatrix<double> (&(*_data)(offset),nbNodes, nbFields*nbElements);
- offset += nbNodes*nbFields*nbElements;
- }
- offset=0;
- for (int i=0;i<_groups.getNbGhostGroups();i++){
- dgGroupOfElements *group = _groups.getGhostGroup(i);
- int nbNodes = group->getNbNodes();
- int nbElements = group->getNbElements();
- _dataProxys[i+_groups.getNbElementGroups()] = new fullMatrix<double> (&(*_ghostData)(offset),nbNodes, nbFields*nbElements);
- offset += nbNodes*nbFields*nbElements;
- }
-}
-
-dgDofContainer::~dgDofContainer (){
- if (!_dataSize)return;
- for (int i=0;i<_dataProxys.size();++i) delete _dataProxys[i];
- delete _data;
-}
-void dgSystemOfEquations::scatter(dgDofContainer *vector){
- //1) count
- int *countS=new int[Msg::GetCommSize()];
- int *shiftS=new int[Msg::GetCommSize()];
- int *countR=new int[Msg::GetCommSize()];
- int *shiftR=new int[Msg::GetCommSize()];
- for(int i=0;i<Msg::GetCommSize();i++){
- countS[i]=0;
- countR[i]=0;
- for(size_t j=0;j<_elementsToSend[i].size();j++){
- countS[i] += _groups.getElementGroup(_elementsToSend[i][j].first)->getNbNodes()*_claw->nbFields();
- }
- }
- for(size_t i=0; i<_groups.getNbGhostGroups(); i++){
- dgGroupOfElements *group = _groups.getGhostGroup(i);
- countR[group->getGhostPartition()]+=group->getNbElements()*group->getNbNodes()*_claw->nbFields();
- }
- shiftS[0]=shiftR[0]=0;
- int totalS=countS[0];
- int totalR=countR[0];
- for(size_t i=1; i<Msg::GetCommSize(); i++){
- shiftS[i]=shiftS[i-1]+countS[i-1];
- shiftR[i]=shiftR[i-1]+countR[i-1];
- totalS+=countS[i];
- totalR+=countR[i];
- }
- double *send=new double[totalS];
- double *recv=new double[totalR];
- //2) fill
- int index=0;
- for(int i=0;i<Msg::GetCommSize();i++){
- for(size_t j=0;j<_elementsToSend[i].size();j++){
- int gid = _elementsToSend[i][j].first;
- int eid = _elementsToSend[i][j].second;
- fullMatrix<double> &sol=vector->getGroupProxy(gid);
- for(int l=0;l<_claw->nbFields();l++){
- for(int k=0;k<sol.size1();k++){
- send[index++] = sol(k,eid*_claw->nbFields()+l);
- //send[index++] = _groups.getElementGroup(gid)->getElement(eid)->getNum()*9+l*3+k;
- }
- }
- }
- }
-
- //3) send
- MPI_Alltoallv(send,countS,shiftS,MPI_DOUBLE,recv,countR,shiftR,MPI_DOUBLE,MPI_COMM_WORLD);
- //4) distribute
- for(int i=0; i< _groups.getNbGhostGroups();i++){
- fullMatrix<double>&sol = vector->getGroupProxy(i+_groups.getNbElementGroups());
- int part = _groups.getGhostGroup(i)->getGhostPartition();
- int ndof = sol.size1()*sol.size2();
- for(int j=0;j<sol.size2();j++){
- for(int k=0;k<sol.size1();k++){
- sol(k,j)=recv[shiftR[part]++];
- //int a=(int)recv[shiftR[part]++];
- //int b= _groups.getGhostGroup(i)->getElement(j/3)->getNum()*9+3*(j%3)+k;
- ///*if(a!=b)*/printf("%i %i\n",a,b);
- }
- }
- }
- delete countS;
- delete countR;
- delete shiftS;
- delete shiftR;
- delete send;
- delete recv;
-}
diff --git a/Solver/dgSystemOfEquations.h b/Solver/dgSystemOfEquations.h
index 81182a3cfb..e8d89dff59 100644
--- a/Solver/dgSystemOfEquations.h
+++ b/Solver/dgSystemOfEquations.h
@@ -10,39 +10,13 @@
#include "dgConservationLaw.h"
#include "Gmsh.h"
#include "dgLimiter.h"
+#include "dgDofContainer.h"
-struct dgDofContainer {
-private:
- dgDofContainer (const dgDofContainer&);
- int totalNbElements;
- int nbFields;
- dgGroupCollection &_groups;
-public:
- int _dataSize; // the full data size i.e. concerning all groups (not ghost, see bellow)
- int _dataSizeGhost;
- std::vector<fullMatrix<double> *> _dataProxys; // proxys
- fullVector<double> * _data; // the full data itself
- fullVector<double> * _ghostData;
- inline int getDataSize(){return _dataSize;}
- inline fullMatrix<double> &getGroupProxy(int gId){ return *(_dataProxys[gId]); }
- dgDofContainer (dgGroupCollection &groups, const dgConservationLaw &claw);
- ~dgDofContainer ();
- int getNbElements() {return totalNbElements;}
-};
class binding;
class dgSystemOfEquations {
-//////////////
- //parallel section, should be moved to a new class
- int *shiftSend,*shiftRecv;
- int *nGhostElements,*nParentElements;
- int totalSend, totalRecv;
- public :
- void scatter(dgDofContainer *solution);
- std::vector< std::vector<std::pair<int,int> > >_elementsToSend; //{group,id} of the elements to send to each proc
-//////////////
private:
// the mesh and the model
GModel *_gm;
@@ -78,8 +52,7 @@ public:
static void registerBindings(binding *b);
inline const fullMatrix<double> getSolutionProxy (int iGroup, int iElement){
- return fullMatrix<double> ( *_solution->_dataProxys [iGroup] ,
- iElement * _claw->nbFields(),_claw->nbFields());
+ return fullMatrix<double> ( _solution->getGroupProxy(iGroup), iElement * _claw->nbFields(),_claw->nbFields());
}
void export_solution_as_is (const std::string &fileName);
void saveSolution (std::string fileName) ;
--
GitLab