Skip to content
Snippets Groups Projects
Commit a56a4507 authored by Jonathan Lambrechts's avatar Jonathan Lambrechts
Browse files

work on parallel

parent ad6994dc
Branches
Tags
No related merge requests found
Hide whitespace changes
Inline Side-by-side
Solver/TESTCASES/validation/WavePulse.lua
+ 2
− 4
View file @ a56a4507
......@@ -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)
......
Solver/dgAlgorithm.cpp
+ 59
− 18
View file @ a56a4507
......@@ -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);
......
Solver/dgAlgorithm.h
+ 7
− 1
View file @ a56a4507
......@@ -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);
};
......
Solver/dgGroupOfElements.cpp
+ 3
− 2
View file @ a56a4507
......@@ -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();
......
Solver/dgGroupOfElements.h
+ 3
− 1
View file @ a56a4507
......@@ -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 {
......
Solver/dgSystemOfEquations.cpp
+ 188
− 10
View file @ a56a4507
......@@ -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;
}
Solver/dgSystemOfEquations.h
+ 20
− 1
View file @ a56a4507
#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)
......
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment