diff --git a/Solver/TESTCASES/validation/WavePulse.lua b/Solver/TESTCASES/validation/WavePulse.lua
index b22218354a043625a1a9b3d3d401bcbb560892b3..abf15298e4c4b6917b57aaa5f6432b9d35e86e51 100644
--- a/Solver/TESTCASES/validation/WavePulse.lua
+++ b/Solver/TESTCASES/validation/WavePulse.lua
@@ -22,10 +22,8 @@ myModel = GModel ()
--myModel:load('box.geo')
--myModel:load('box.msh')
--myModel:load('square_quads.msh')
-myModel:load('square.geo')
-myModel:load('test.msh')
-
-
+myModel:load('square_part.msh')
+myModel:mesh(2)
print'*** Create a dg solver ***'
DG = dgSystemOfEquations (myModel)
diff --git a/Solver/dgAlgorithm.cpp b/Solver/dgAlgorithm.cpp
index 7f76ede9df3262e6f3f6c53058e3d9069acfe966..c05c3e8571196b43214699ce246121feb3550323 100644
--- a/Solver/dgAlgorithm.cpp
+++ b/Solver/dgAlgorithm.cpp
@@ -258,11 +258,14 @@ void dgAlgorithm::rungeKutta (const dgConservationLaw &claw, // conservation l
std::vector<dgGroupOfElements*> &eGroups, // group of elements
std::vector<dgGroupOfFaces*> &fGroups, // group of interfacs
std::vector<dgGroupOfFaces*> &bGroups, // group of boundaries
+ std::vector<dgGroupOfElements*> &ghostGroups, // group of boundaries
double h, // time-step
dgDofContainer &sol,
dgDofContainer &resd,
+ dgSystemOfEquations *syst,
dgLimiter *limiter,
- int orderRK) // order of RK integrator
+ int orderRK
+ ) // order of RK integrator
{
// U_{n+1}=U_n+h*(SUM_i a_i*K_i)
@@ -275,16 +278,14 @@ void dgAlgorithm::rungeKutta (const dgConservationLaw &claw, // conservation l
a[0] = h;
}
- // fullMatrix<double> K(sol);
// Current updated solution
- // fullMatrix<double> Unp(sol);
fullMatrix<double> residuEl, KEl;
fullMatrix<double> iMassEl;
int nbFields = claw.nbFields();
- dgDofContainer K (eGroups,claw);
- dgDofContainer Unp (eGroups,claw);
+ dgDofContainer K (eGroups,ghostGroups,claw);
+ dgDofContainer Unp (eGroups,ghostGroups,claw);
K._data->scale(0.0);
K._data->axpy(*(sol._data));
@@ -300,7 +301,8 @@ void dgAlgorithm::rungeKutta (const dgConservationLaw &claw, // conservation l
if (limiter){
limiter->apply(K, eGroups, fGroups);
}
- this->residual(claw,eGroups,fGroups,bGroups,K._dataProxys,resd._dataProxys);
+ syst->scatter(&K);
+ this->residual(claw,eGroups,fGroups,bGroups,ghostGroups,K._dataProxys,resd._dataProxys);
K._data->scale(0.);
for(int k=0;k < eGroups.size();k++) {
int nbNodes = eGroups[k]->getNbNodes();
@@ -325,6 +327,7 @@ void dgAlgorithm::multirateRungeKutta (const dgConservationLaw &claw, // conse
std::vector<dgGroupOfElements*> &eGroups, // group of elements
std::vector<dgGroupOfFaces*> &fGroups, // group of interfacs
std::vector<dgGroupOfFaces*> &bGroups, // group of boundaries
+ std::vector<dgGroupOfElements*> &ghostGroups, // group of boundaries
double h, // time-step
dgDofContainer &sol,
dgDofContainer &resd,
@@ -429,7 +432,7 @@ void dgAlgorithm::multirateRungeKutta (const dgConservationLaw &claw, // conse
}
}
}
- this->residual(claw,eGroups,fGroups,bGroups,tmp._dataProxys,resd._dataProxys);
+ this->residual(claw,eGroups,fGroups,bGroups,ghostGroups,tmp._dataProxys,resd._dataProxys);
for(int k=0;k < eGroups.size();k++) {
int nbNodes = eGroups[k]->getNbNodes();
for(int i=0;i<eGroups[k]->getNbElements();i++) {
@@ -568,14 +571,20 @@ void dgAlgorithm::partitionGroup(dgGroupOfElements &eGroup,
void dgAlgorithm::buildGroups(GModel *model, int dim, int order,
std::vector<dgGroupOfElements*> &eGroups,
std::vector<dgGroupOfFaces*> &fGroups,
- std::vector<dgGroupOfFaces*> &bGroups)
+ std::vector<dgGroupOfFaces*> &bGroups,
+ std::vector<dgGroupOfElements*> &ghostGroups
+ )
{
std::map<const std::string,std::set<MVertex*> > boundaryVertices;
std::map<const std::string,std::set<MEdge, Less_Edge> > boundaryEdges;
std::map<const std::string,std::set<MFace, Less_Face> > boundaryFaces;
std::vector<GEntity*> entities;
model->getEntities(entities);
- std::map<int, std::vector<MElement *> >allElements;
+ std::map<int, std::vector<MElement *> >localElements;
+ std::vector<std::map<int, std::vector<MElement *> > >ghostElements(Msg::GetCommSize()); // [partition][elementType]
+ int nlocal=0;
+ int nghosts=0;
+ std::multimap<MElement*, short> &ghostsCells = model->getGhostCells();
for(unsigned int i = 0; i < entities.size(); i++){
GEntity *entity = entities[i];
if(entity->dim() == dim-1){
@@ -599,16 +608,27 @@ void dgAlgorithm::buildGroups(GModel *model, int dim, int order,
}
}
}else if(entity->dim() == dim){
- for (int iel=0; iel<entity->getNumMeshElements(); iel++)
- allElements[entity->getMeshElement(iel)->getType()].push_back(entity->getMeshElement(iel));
+ for (int iel=0; iel<entity->getNumMeshElements(); iel++){
+ MElement *el=entity->getMeshElement(iel);
+ if(el->getPartition()==Msg::GetCommRank()+1 || el->getPartition()==0){
+ localElements[el->getType()].push_back(el);
+ nlocal++;
+ }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);
+ ghost++;
+ }
+ }
+ }
}
}
-
- // do a group of elements for every element type that is present in the mesh
- Msg::Info("%d groups of elements",allElements.size());
- for (std::map<int, std::vector<MElement *> >::iterator it = allElements.begin(); it !=allElements.end() ; ++it){
-
+
+ Msg::Info("%d groups of elements",localElements.size());
+ // do a group of elements for every element type in the mesh
+ for (std::map<int, std::vector<MElement *> >::iterator it = localElements.begin(); it !=localElements.end() ; ++it){
eGroups.push_back(new dgGroupOfElements(it->second,order));
}
@@ -644,9 +664,25 @@ void dgAlgorithm::buildGroups(GModel *model, int dim, int order,
for (int j=i+1;j<eGroups.size();j++){
dgGroupOfFaces *gof = new dgGroupOfFaces(*eGroups[i],*eGroups[j],order);
if (gof->getNbElements())
- fGroups.push_back(gof);
+ fGroups.push_back(gof);
+ else
+ delete gof;
+ }
+ }
+ //create ghost groups
+ for(int i=0;i<Msg::GetCommSize();i++){
+ for (std::map<int, std::vector<MElement *> >::iterator it = ghostElements[i].begin(); it !=ghostElements[i].end() ; ++it){
+ ghostGroups.push_back(new dgGroupOfElements(it->second,order,i));
+ }
+ }
+ //create face group for ghostGroups
+ for (int i=0; i<ghostGroups.size(); i++){
+ for (int j=0;j<eGroups.size();j++){
+ dgGroupOfFaces *gof = new dgGroupOfFaces(*ghostGroups[i],*eGroups[j],order);
+ if (gof->getNbElements())
+ fGroups.push_back(gof);
else
- delete gof;
+ delete gof;
}
}
}
@@ -655,6 +691,7 @@ void dgAlgorithm::residual( const dgConservationLaw &claw,
std::vector<dgGroupOfElements*> &eGroups, //group of elements
std::vector<dgGroupOfFaces*> &fGroups, // group of interfacs
std::vector<dgGroupOfFaces*> &bGroups, // group of boundaries
+ std::vector<dgGroupOfElements*> &ghostGroups, // group of boundaries
std::vector<fullMatrix<double> *> &solution, // solution
std::vector<fullMatrix<double> *> &residu) // residual
{
@@ -674,6 +711,10 @@ void dgAlgorithm::residual( const dgConservationLaw &claw,
if (eGroups[j] == &faces.getGroupLeft())iGroupLeft = j;
if (eGroups[j] == &faces.getGroupRight())iGroupRight= j;
}
+ for(size_t j=0;j<ghostGroups.size() ; j++) {
+ if (ghostGroups[j] == &faces.getGroupLeft())iGroupLeft = j+eGroups.size();
+ if (ghostGroups[j] == &faces.getGroupRight())iGroupRight= j+eGroups.size();
+ }
fullMatrix<double> solInterface(faces.getNbNodes(),faces.getNbElements()*2*nbFields);
fullMatrix<double> residuInterface(faces.getNbNodes(),faces.getNbElements()*2*nbFields);
faces.mapToInterface(nbFields, *solution[iGroupLeft], *solution[iGroupRight], solInterface);
diff --git a/Solver/dgAlgorithm.h b/Solver/dgAlgorithm.h
index 2402b2a6093ee6ea6889fb91a7b1531faae5ddc7..2ca505a71d7d60d1f192669220d908ebca0a7231 100644
--- a/Solver/dgAlgorithm.h
+++ b/Solver/dgAlgorithm.h
@@ -10,6 +10,7 @@ class dgConservationLaw;
class dgDofContainer;
class dgTerm;
class dgLimiter;
+class dgSystemOfEquations;
class dgAlgorithm {
public :
@@ -37,6 +38,7 @@ class dgAlgorithm {
std::vector<dgGroupOfElements*> &eGroups, //group of elements
std::vector<dgGroupOfFaces*> &fGroups, // group of interfacs
std::vector<dgGroupOfFaces*> &bGroups, // group of boundaries
+ std::vector<dgGroupOfElements*> &ghostGroups, // group of boundaries
std::vector<fullMatrix<double> *> &solution, // solution !! at faces nodes
std::vector<fullMatrix<double> *> &residual // residual !! at faces nodes
);
@@ -45,9 +47,11 @@ class dgAlgorithm {
std::vector<dgGroupOfElements*> &eGroups, //group of elements
std::vector<dgGroupOfFaces*> &fGroups, // group of interfacs
std::vector<dgGroupOfFaces*> &bGroups, // group of boundaries
+ std::vector<dgGroupOfElements*> &ghostGroups, // group of boundaries
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)
@@ -61,6 +65,7 @@ class dgAlgorithm {
std::vector<dgGroupOfElements*> &eGroups, //group of elements
std::vector<dgGroupOfFaces*> &fGroups, // group of interfacs
std::vector<dgGroupOfFaces*> &bGroups, // group of boundaries
+ std::vector<dgGroupOfElements*> &ghostGroups, // group of boundaries
double h,
dgDofContainer &residu,
dgDofContainer &sol,
@@ -73,7 +78,8 @@ class dgAlgorithm {
static void buildGroups(GModel *model, int dim, int order,
std::vector<dgGroupOfElements*> &eGroups,
std::vector<dgGroupOfFaces*> &fGroups,
- std::vector<dgGroupOfFaces*> &bGroups);
+ std::vector<dgGroupOfFaces*> &bGroups,
+ std::vector<dgGroupOfElements*> &ghostGroups);
};
diff --git a/Solver/dgGroupOfElements.cpp b/Solver/dgGroupOfElements.cpp
index 5b746696452c4ef8e7272440ae1c8ef2e01d6dc7..fc5080aad90b0d4b6fed6138e12279a9cd3aeb52 100644
--- a/Solver/dgGroupOfElements.cpp
+++ b/Solver/dgGroupOfElements.cpp
@@ -43,10 +43,11 @@ static fullMatrix<double> dgGetFaceIntegrationRuleOnElement (
}
-dgGroupOfElements::dgGroupOfElements(const std::vector<MElement*> &e, int polyOrder)
+dgGroupOfElements::dgGroupOfElements(const std::vector<MElement*> &e, int polyOrder,int ghostPartition)
: _elements(e),
_fs(*_elements[0]->getFunctionSpace(polyOrder)),
- _integration(dgGetIntegrationRule (_elements[0], polyOrder))
+ _integration(dgGetIntegrationRule (_elements[0], polyOrder)),
+ _ghostPartition(ghostPartition)
{
_order=polyOrder;
_dimUVW = _dimXYZ = e[0]->getDim();
diff --git a/Solver/dgGroupOfElements.h b/Solver/dgGroupOfElements.h
index c1f7de0c26fa2e8cc963dc7c24feadbe8ccfb2fe..61c9925384972a46564b28df85da4132d69ce5be 100644
--- a/Solver/dgGroupOfElements.h
+++ b/Solver/dgGroupOfElements.h
@@ -38,6 +38,7 @@ public:
// store topological and geometrical data for 1 group for 1 discretisation
class dgGroupOfElements {
+ int _ghostPartition; // -1 : this is not a ghosted group, otherwise the id of the parent partition
// N elements in the group
std::vector<MElement*> _elements;
// the ONLY function space that is used to
@@ -72,7 +73,7 @@ class dgGroupOfElements {
// dgGroupOfElements (const dgGroupOfElements &e, int order) {}
// dgGroupOfElements & operator = (const dgGroupOfElements &e) {}
public:
- dgGroupOfElements (const std::vector<MElement*> &e, int pOrder);
+ dgGroupOfElements (const std::vector<MElement*> &e, int pOrder,int ghostPartition=-1);
virtual ~dgGroupOfElements ();
inline int getNbElements() const {return _elements.size();}
inline int getNbNodes() const {return _collocation->size2();}
@@ -91,6 +92,7 @@ public:
inline fullMatrix<double> &getInverseMassMatrix () const {return *_imass;}
inline const fullMatrix<double> getMapping (int iElement) const {return fullMatrix<double>(*_mapping, iElement, 1);}
inline int getOrder() const {return _order;}
+ inline int getGhostPartition() const {return _ghostPartition;}
};
class dgGroupOfFaces {
diff --git a/Solver/dgSystemOfEquations.cpp b/Solver/dgSystemOfEquations.cpp
index f1598de057abdc1f83b8d8a150f97830157618e2..45f7ec2b170a19f6d1bc2809c2093bca6b63aeb2 100644
--- a/Solver/dgSystemOfEquations.cpp
+++ b/Solver/dgSystemOfEquations.cpp
@@ -7,6 +7,9 @@
#include "PView.h"
#include "PViewData.h"
#include "dgLimiter.h"
+#ifdef HAVE_MPI
+#include "mpi.h"
+#endif
class dgConservationLawL2Projection : public dgConservationLaw {
std::string _functionName;
@@ -100,15 +103,68 @@ void dgSystemOfEquations::setup(){
_order,
_elementGroups,
_faceGroups,
- _boundaryGroups);
+ _boundaryGroups,
+ _ghostGroups
+ );
+ // 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< _ghostGroups.size(); i++){
+ nGhostElements[_ghostGroups[i]->getGhostPartition()] += _ghostGroups[i]->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 (size_t i = 0; i< _ghostGroups.size(); i++){
+ int part = _ghostGroups[i]->getGhostPartition();
+ for (int j=0; j< _ghostGroups[i]->getNbElements() ; j++){
+ idSend[curShiftSend[part]++] = _ghostGroups[i]->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<_elementGroups.size(); i++) {
+ for(int j=0; j<_elementGroups[i]->getNbElements(); j++){
+ elementMap[_elementGroups[i]->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
- _solution = new dgDofContainer(_elementGroups,*_claw);
- _rightHandSide = new dgDofContainer(_elementGroups,*_claw);
+ delete curShiftSend;
+ _solution = new dgDofContainer(_elementGroups,_ghostGroups,*_claw);
+ _rightHandSide = new dgDofContainer(_elementGroups,_ghostGroups,*_claw);
}
double dgSystemOfEquations::RK44(double dt){
- _algo->rungeKutta(*_claw, _elementGroups, _faceGroups, _boundaryGroups, dt, *_solution, *_rightHandSide, NULL);
+ _algo->rungeKutta(*_claw, _elementGroups, _faceGroups, _boundaryGroups,_ghostGroups, dt, *_solution, *_rightHandSide, this);
return _solution->_data->norm();
}
@@ -119,22 +175,29 @@ double dgSystemOfEquations::computeInvSpectralRadius(){
_algo->computeElementaryTimeSteps(*_claw, *_elementGroups[i], *_solution->_dataProxys[i], DTS);
for (int k=0;k<DTS.size();k++) sr = std::min(sr,DTS[k]);
}
- return sr;
+ #ifdef HAVE_MPI
+ double sr_min;
+ MPI_Allreduce((void *)&sr, &sr_min, 1, MPI_DOUBLE, MPI_MIN,
+ MPI_COMM_WORLD);
+ return sr_min;
+ #else
+ return sr
+ #endif
}
double dgSystemOfEquations::RK44_limiter(double dt){
dgLimiter *sl = new dgSlopeLimiter(_claw);
- _algo->rungeKutta(*_claw, _elementGroups, _faceGroups, _boundaryGroups, dt, *_solution, *_rightHandSide, sl);
+ _algo->rungeKutta(*_claw, _elementGroups, _faceGroups, _boundaryGroups,_ghostGroups, dt, *_solution, *_rightHandSide,this, sl);
delete sl;
return _solution->_data->norm();
}
double dgSystemOfEquations::ForwardEuler(double dt){
- _algo->rungeKutta(*_claw, _elementGroups, _faceGroups, _boundaryGroups, dt, *_solution, *_rightHandSide, NULL,1);
+ _algo->rungeKutta(*_claw, _elementGroups, _faceGroups, _boundaryGroups,_ghostGroups, dt, *_solution, *_rightHandSide,this, NULL,1);
return _solution->_data->norm();
}
double dgSystemOfEquations::multirateRK43(double dt){
- _algo->multirateRungeKutta(*_claw, _elementGroups, _faceGroups, _boundaryGroups, dt, *_solution, *_rightHandSide);
+ _algo->multirateRungeKutta(*_claw, _elementGroups, _faceGroups, _boundaryGroups,_ghostGroups, dt, *_solution, *_rightHandSide);
return _solution->_data->norm();
}
@@ -177,6 +240,8 @@ void dgSystemOfEquations::export_solution_as_is (const std::string &name){
for (int ICOMP = 0; ICOMP<_claw->nbFields();++ICOMP){
std::ostringstream name_oss;
name_oss<<name<<"_COMP_"<<ICOMP<<".msh";
+ if(Msg::GetCommSize()>1)
+ name_oss<<"_"<<Msg::GetCommRank();
FILE *f = fopen (name_oss.str().c_str(),"w");
int COUNT = 0;
for (int i=0;i < _elementGroups.size() ;i++){
@@ -232,9 +297,11 @@ void dgSystemOfEquations::export_solution_as_is (const std::string &name){
delete pv;
}
-dgDofContainer::dgDofContainer (std::vector<dgGroupOfElements*> &elementGroups, const dgConservationLaw &claw){
+dgDofContainer::dgDofContainer (std::vector<dgGroupOfElements*> &elementGroups,std::vector<dgGroupOfElements*> ghostGroups, const dgConservationLaw &claw){
_dataSize = 0;
+ _dataSizeGhost=0;
totalNbElements = 0;
+ int totalNbElementsGhost =0;
nbFields = claw.nbFields();
for (int i=0;i<elementGroups.size();i++){
int nbNodes = elementGroups[i]->getNbNodes();
@@ -242,15 +309,57 @@ dgDofContainer::dgDofContainer (std::vector<dgGroupOfElements*> &elementGroups,
totalNbElements +=nbElements;
_dataSize += nbNodes*nbFields*nbElements;
}
+ for (int i=0;i<ghostGroups.size();i++){
+ int nbNodes = ghostGroups[i]->getNbNodes();
+ int nbElements = ghostGroups[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(elementGroups.size());
+ _dataProxys.resize(elementGroups.size()+ghostGroups.size());
for (int i=0;i<elementGroups.size();i++){
+ dgGroupOfElements *group=elementGroups[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<ghostGroups.size();i++){
+ dgGroupOfElements *group=ghostGroups[i];
+ int nbNodes = group->getNbNodes();
+ int nbElements = group->getNbElements();
+ _dataProxys[i+elementGroups.size()] = new fullMatrix<double> (&(*_ghostData)(offset),nbNodes, nbFields*nbElements);
+ offset += nbNodes*nbFields*nbElements;
+ }
+ // printf("datasize = %d\n",_dataSize);
+}
+
+dgDofContainer::dgDofContainer (std::vector<dgGroupOfElements*> &elementGroups, const dgConservationLaw &claw){
+ _dataSize = 0;
+ totalNbElements = 0;
+ nbFields = claw.nbFields();
+ for (int i=0;i<elementGroups.size();i++){
int nbNodes = elementGroups[i]->getNbNodes();
int nbElements = elementGroups[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(elementGroups.size());
+ for (int i=0;i<elementGroups.size();i++){
+ dgGroupOfElements *group=elementGroups[i];
+ int nbNodes = group->getNbNodes();
+ int nbElements = group->getNbElements();
_dataProxys[i] = new fullMatrix<double> (&(*_data)(offset),nbNodes, nbFields*nbElements);
offset += nbNodes*nbFields*nbElements;
}
@@ -262,3 +371,72 @@ dgDofContainer::~dgDofContainer (){
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] += _elementGroups[_elementsToSend[i][j].first]->getNbNodes()*_claw->nbFields();
+ }
+ }
+ for(size_t i=0; i<_ghostGroups.size(); i++){
+ dgGroupOfElements *group=_ghostGroups[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++] = _elementGroups[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< _ghostGroups.size();i++){
+ fullMatrix<double>&sol = vector->getGroupProxy(i+_elementGroups.size());
+ int part = _ghostGroups[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= _ghostGroups[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 f92ae8b44e75d4224788e468387373808bab85b3..e8f4bcd583240a09b19f8ba8a9b774a8e9bdadc5 100644
--- a/Solver/dgSystemOfEquations.h
+++ b/Solver/dgSystemOfEquations.h
@@ -1,5 +1,7 @@
#ifndef _DG_SYSTEM_OF_EQUATIONS_
#define _DG_SYSTEM_OF_EQUATIONS_
+#include <vector>
+#include <utility>
#include "GmshConfig.h"
#include "GModel.h"
#include "dgAlgorithm.h"
@@ -15,19 +17,33 @@ private:
int totalNbElements;
int nbFields;
public:
- int _dataSize; // the full data size i.e. concerning all groups
+ 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 (std::vector<dgGroupOfElements*> &groups,std::vector<dgGroupOfElements*> ghostGroups, const dgConservationLaw &claw);
dgDofContainer (std::vector<dgGroupOfElements*> &groups, const dgConservationLaw &claw);
~dgDofContainer ();
int getNbElements() {return totalNbElements;}
+ //collective, should be called on all proc
+ void setGhostedGroups (std::vector<dgGroupOfElements*> &ghostGroups);
};
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);
+//////////////
+ private:
// the mesh and the model
GModel *_gm;
// the algorithm that computes DG operators
@@ -44,6 +60,9 @@ class dgSystemOfEquations {
dgDofContainer *_rightHandSide;
// groups of elements (volume terms)
std::vector<dgGroupOfElements*> _elementGroups;
+ //ghost structure
+ std::vector<dgGroupOfElements*> _ghostGroups;
+ std::vector< std::vector<std::pair<int,int> > >_elementsToSend; //{group,id} of the elements to send to each proc
// groups of faces (interface terms)
std::vector<dgGroupOfFaces*> _faceGroups;
// groups of faces (boundary conditions)