diff --git a/Solver/dgAlgorithm.cpp b/Solver/dgAlgorithm.cpp
index 462ea79ade56dd6b6799c4d881618c1f73b705d7..fd3ad00617e947c8b78a57515dfdbef683b60c03 100644
--- a/Solver/dgAlgorithm.cpp
+++ b/Solver/dgAlgorithm.cpp
@@ -282,44 +282,39 @@ void dgAlgorithm::rungeKutta (const dgConservationLaw &claw, // conservation l
dgDofContainer K (groups,nbFields);
dgDofContainer Unp (groups,nbFields);
-
- K._data->scale(0.0);
- K._data->axpy(*(sol._data));
- Unp._data->scale(0.0);
- Unp._data->axpy(*(sol._data));
+ K.scale(0.);
+ K.axpy(sol);
+ Unp.scale(0.);
+ Unp.axpy(sol);
for(int j=0; j<orderRK;j++){
if(j){
- K._data->scale(b[j]);
- K._data->axpy(*(sol._data));
+ K.scale(b[j]);
+ K.axpy(sol);
}
if (limiter){
limiter->apply(K, groups);
}
- K.scatter();
- this->residual(claw,groups,K._dataProxys,resd._dataProxys);
- K._data->scale(0.);
+ this->residual(claw,groups,K,resd);
+ K.scale(0.);
for(int k=0; k < groups.getNbElementGroups(); k++) {
dgGroupOfElements *group = groups.getElementGroup(k);
int nbNodes = group->getNbNodes();
for(int i=0;i<group->getNbElements();i++) {
- residuEl.setAsProxy(*(resd._dataProxys[k]),i*nbFields,nbFields);
- KEl.setAsProxy(*(K._dataProxys[k]),i*nbFields,nbFields);
+ residuEl.setAsProxy(resd.getGroupProxy(k),i*nbFields,nbFields);
+ KEl.setAsProxy(K.getGroupProxy(k),i*nbFields,nbFields);
iMassEl.setAsProxy(group->getInverseMassMatrix(),i*nbNodes,nbNodes);
iMassEl.mult(residuEl,KEl);
}
}
- Unp._data->axpy(*(K._data),a[j]);
+ Unp.axpy(K,a[j]);
}
if (limiter) limiter->apply(Unp, groups);
- for (int i=0;i<sol._dataSize;i++){
- (*sol._data)(i)=(*Unp._data)(i);
- }
+ sol.scale(0.);
+ sol.axpy(Unp);
}
-
-
void dgAlgorithm::multirateRungeKutta (const dgConservationLaw &claw, // conservation law
dgGroupCollection &groups,
double h, // time-step
@@ -408,17 +403,17 @@ void dgAlgorithm::multirateRungeKutta (const dgConservationLaw &claw, // conse
K=new dgDofContainer*[nStages];
for(int i=0;i<nStages;i++){
K[i]=new dgDofContainer(groups,nbFields);
- K[i]->_data->scale(0.0);
+ K[i]->scale(0.);
}
dgDofContainer Unp (groups,nbFields);
dgDofContainer tmp (groups,nbFields);
- Unp._data->scale(0.0);
- Unp._data->axpy(*(sol._data));
+ Unp.scale(0.0);
+ Unp.axpy(sol);
for(int j=0; j<nStages;j++){
- tmp._data->scale(0.0);
- tmp._data->axpy(*(sol._data));
+ tmp.scale(0.0);
+ tmp.axpy(sol);
for(int k=0;k < groups.getNbElementGroups();k++) {
for(int i=0;i<j;i++){
if(fabs(A[j][i])>1e-12){
@@ -426,23 +421,22 @@ void dgAlgorithm::multirateRungeKutta (const dgConservationLaw &claw, // conse
}
}
}
- this->residual(claw,groups,tmp._dataProxys,resd._dataProxys);
+ this->residual(claw,groups,tmp,resd);
for(int k=0;k < groups.getNbElementGroups(); k++) {
dgGroupOfElements *group = groups.getElementGroup(k);
int nbNodes = group->getNbNodes();
for(int i=0;i<group->getNbElements();i++) {
- residuEl.setAsProxy(*(resd._dataProxys[k]),i*nbFields,nbFields);
- KEl.setAsProxy(*(K[j]->_dataProxys[k]),i*nbFields,nbFields);
+ residuEl.setAsProxy(resd.getGroupProxy(k),i*nbFields,nbFields);
+ KEl.setAsProxy(K[j]->getGroupProxy(k),i*nbFields,nbFields);
iMassEl.setAsProxy(group->getInverseMassMatrix(),i*nbNodes,nbNodes);
iMassEl.mult(residuEl,KEl);
}
Unp.getGroupProxy(k).add(K[j]->getGroupProxy(k),h*b[j]);
}
}
-
- for (int i=0;i<sol._dataSize;i++){
- (*sol._data)(i)=(*Unp._data)(i);
- }
+ sol.scale(0.);
+ sol.axpy(Unp);
+
for(int i=0;i<nStages;i++){
delete K[i];
}
@@ -536,14 +530,15 @@ void dgAlgorithm::residualBoundary ( //dofManager &dof, // the DOF manager (mayb
// works for any number of groups
void dgAlgorithm::residual( const dgConservationLaw &claw,
dgGroupCollection &groups,
- std::vector<fullMatrix<double> *> &solution, // solution
- std::vector<fullMatrix<double> *> &residu) // residual
+ dgDofContainer &solution,
+ dgDofContainer &residu)
{
+ solution.scatter();
int nbFields=claw.nbFields();
//volume term
for(size_t i=0;i<groups.getNbElementGroups() ; i++) {
- residu[i]->scale(0);
- residualVolume(claw,*groups.getElementGroup(i),*solution[i],*residu[i]);
+ residu.getGroupProxy(i).scale(0);
+ residualVolume(claw,*groups.getElementGroup(i),solution.getGroupProxy(i),residu.getGroupProxy(i));
}
// residu[0]->print("Volume");
//interface term
@@ -563,11 +558,10 @@ void dgAlgorithm::residual( const dgConservationLaw &claw,
}
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);
- residualInterface(claw,faces,solInterface,*solution[iGroupLeft], *solution[iGroupRight],residuInterface);
- faces.mapFromInterface(nbFields, residuInterface, *residu[iGroupLeft], *residu[iGroupRight]);
+ faces.mapToInterface(nbFields, solution.getGroupProxy(iGroupLeft), solution.getGroupProxy(iGroupRight), solInterface);
+ residualInterface(claw,faces,solInterface,solution.getGroupProxy(iGroupLeft), solution.getGroupProxy(iGroupRight),residuInterface);
+ faces.mapFromInterface(nbFields, residuInterface,residu.getGroupProxy(iGroupLeft), residu.getGroupProxy(iGroupRight));
}
- // residu[0]->print("Interfaces");
//boundaries
for(size_t i=0;i<groups.getNbBoundaryGroups() ; i++) {
dgGroupOfFaces &faces = *groups.getBoundaryGroup(i);
@@ -579,12 +573,10 @@ void dgAlgorithm::residual( const dgConservationLaw &claw,
}
fullMatrix<double> solInterface(faces.getNbNodes(),faces.getNbElements()*nbFields);
fullMatrix<double> residuInterface(faces.getNbNodes(),faces.getNbElements()*nbFields);
- faces.mapToInterface(nbFields, *solution[iGroupLeft], *solution[iGroupRight], solInterface);
- residualBoundary(claw,faces,solInterface,*solution[iGroupLeft],residuInterface);
- faces.mapFromInterface(nbFields, residuInterface, *residu[iGroupLeft], *residu[iGroupRight]);
+ faces.mapToInterface(nbFields, solution.getGroupProxy(iGroupLeft), solution.getGroupProxy(iGroupRight), solInterface);
+ residualBoundary(claw,faces,solInterface,solution.getGroupProxy(iGroupLeft),residuInterface);
+ faces.mapFromInterface(nbFields, residuInterface, residu.getGroupProxy(iGroupLeft), residu.getGroupProxy(iGroupRight));
}
-
- // residu[0]->print("Boundaries");
}
void dgAlgorithm::computeElementaryTimeSteps ( //dofManager &dof, // the DOF manager (maybe useless here)
diff --git a/Solver/dgAlgorithm.h b/Solver/dgAlgorithm.h
index 206bfcfd75d5ee332aac2c6a288b7dd86e122890..25a37420ee7a09bad5ae1b97c4ebd9b36546770f 100644
--- a/Solver/dgAlgorithm.h
+++ b/Solver/dgAlgorithm.h
@@ -36,9 +36,7 @@ class dgAlgorithm {
);
// works only if there is only 1 group of element
static void residual( const dgConservationLaw &claw, dgGroupCollection &groups,
- std::vector<fullMatrix<double> *> &solution, // solution !! at faces nodes
- std::vector<fullMatrix<double> *> &residual // residual !! at faces nodes
- );
+ dgDofContainer &solution, dgDofContainer &residual);
void rungeKutta (
const dgConservationLaw &claw,
dgGroupCollection &groups,
diff --git a/Solver/dgDofContainer.cpp b/Solver/dgDofContainer.cpp
index 4464e03d37d4e7f05e661a128d200ef69dc3fb6c..34aa06db1c578717e1a69db2e56c15c5ec2365e1 100644
--- a/Solver/dgDofContainer.cpp
+++ b/Solver/dgDofContainer.cpp
@@ -9,7 +9,7 @@
dgDofContainer::dgDofContainer (dgGroupCollection &groups, int nbFields):
_groups(groups)
{
- _dataSize = 0;
+ int _dataSize = 0;
_dataSizeGhost=0;
_totalNbElements = 0;
_parallelStructureExists = false;
@@ -67,7 +67,6 @@ dgDofContainer::~dgDofContainer (){
delete []sendBuf;
delete []recvBuf;
}
- if (!_dataSize)return;
for (int i=0;i<_dataProxys.size();++i) delete _dataProxys[i];
delete _data;
delete _ghostData;
@@ -143,3 +142,35 @@ void dgDofContainer::scatter() {
sol.setAll(recvProxy);
}
}
+void dgDofContainer::scale(double f)
+{
+ _data->scale(f);
+ _ghostData->scale(f);
+}
+
+void dgDofContainer::axpy(dgDofContainer &x, double a)
+{
+ _data->axpy(*x._data,a);
+ _ghostData->axpy(*x._ghostData,a);
+}
+
+double dgDofContainer::norm() {
+ double localNorm = _data->norm();
+ #ifdef HAVE_MPI
+ double globalNorm;
+ MPI_Allreduce(&localNorm, &globalNorm, 1, MPI_DOUBLE, MPI_SUM, MPI_COMM_WORLD);
+ return globalNorm;
+ #else
+ return localNorm;
+ #endif
+}
+void dgDofContainer::save(const std::string name) {
+ FILE *f = fopen (name.c_str(),"rb");
+ _data->binaryLoad(f);
+ fclose(f);
+}
+void dgDofContainer::load(const std::string name) {
+ FILE *f = fopen (name.c_str(),"rb");
+ _data->binaryLoad(f);
+ fclose(f);
+}
diff --git a/Solver/dgDofContainer.h b/Solver/dgDofContainer.h
index 7d4298507daca491276719d003c2c3da58c76ac5..8eaf22b0e31e4f80cc10c94cc3e8fdf728fd4ccf 100644
--- a/Solver/dgDofContainer.h
+++ b/Solver/dgDofContainer.h
@@ -12,23 +12,25 @@ private:
int _nbFields;
dgGroupCollection &_groups;
int _dataSizeGhost;
+ fullVector<double> * _data; // the full data itself
fullVector<double> * _ghostData;
- inline int getDataSize(){return _dataSize;}
+ //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)
-
+public:
+ void scale(double f);
+ double norm();
+ void axpy(dgDofContainer &x, double a=1.);
inline fullMatrix<double> &getGroupProxy(int gId){ return *(_dataProxys[gId]); }
dgDofContainer (dgGroupCollection &groups, int nbFields);
~dgDofContainer ();
int getNbElements() {return _totalNbElements;}
void scatter();
+ void save(const std::string name);
+ void load(const std::string name);
};
#endif
diff --git a/Solver/dgSystemOfEquations.cpp b/Solver/dgSystemOfEquations.cpp
index d369c32f9bb15b99b335a79d20d1fd7b34961d8f..d036a13ea85246e5d22f2d99cf1a6acf05c751cd 100644
--- a/Solver/dgSystemOfEquations.cpp
+++ b/Solver/dgSystemOfEquations.cpp
@@ -30,7 +30,6 @@ dgSystemOfEquations::dgSystemOfEquations(GModel *gm){
_solution = 0;
}
-
// set the order of interpolation
void dgSystemOfEquations::setOrder(int o){
_order = o;
@@ -91,8 +90,8 @@ void dgSystemOfEquations::L2Projection (std::string functionName){
dgConservationLawL2Projection Law(functionName,*_claw);
for (int i=0;i<_groups.getNbElementGroups();i++){
dgGroupOfElements *group = _groups.getElementGroup(i);
- _algo->residualVolume(Law,*group,*_solution->_dataProxys[i],*_rightHandSide->_dataProxys[i]);
- _algo->multAddInverseMassMatrix(*group,*_rightHandSide->_dataProxys[i],*_solution->_dataProxys[i]);
+ _algo->residualVolume(Law,*group,_solution->getGroupProxy(i),_rightHandSide->getGroupProxy(i));
+ _algo->multAddInverseMassMatrix(*group,_rightHandSide->getGroupProxy(i),_solution->getGroupProxy(i));
}
}
@@ -107,14 +106,14 @@ void dgSystemOfEquations::setup(){
double dgSystemOfEquations::RK44(double dt){
_algo->rungeKutta(*_claw,_groups, dt, *_solution, *_rightHandSide);
- return _solution->_data->norm();
+ return _solution->norm();
}
double dgSystemOfEquations::computeInvSpectralRadius(){
double sr = 1.e22;
for (int i=0;i<_groups.getNbElementGroups();i++){
std::vector<double> DTS;
- _algo->computeElementaryTimeSteps(*_claw, *_groups.getElementGroup(i), *_solution->_dataProxys[i], DTS);
+ _algo->computeElementaryTimeSteps(*_claw, *_groups.getElementGroup(i), _solution->getGroupProxy(i), DTS);
for (int k=0;k<DTS.size();k++) sr = std::min(sr,DTS[k]);
}
#ifdef HAVE_MPI
@@ -130,16 +129,16 @@ double dgSystemOfEquations::RK44_limiter(double dt){
dgLimiter *sl = new dgSlopeLimiter(_claw);
_algo->rungeKutta(*_claw,_groups, dt, *_solution, *_rightHandSide, sl);
delete sl;
- return _solution->_data->norm();
+ return _solution->norm();
}
double dgSystemOfEquations::ForwardEuler(double dt){
_algo->rungeKutta(*_claw, _groups, dt, *_solution, *_rightHandSide, NULL,1);
- return _solution->_data->norm();
+ return _solution->norm();
}
double dgSystemOfEquations::multirateRK43(double dt){
_algo->multirateRungeKutta(*_claw, _groups, dt, *_solution, *_rightHandSide);
- return _solution->_data->norm();
+ return _solution->norm();
}
void dgSystemOfEquations::exportSolution(std::string outputFile){
@@ -161,15 +160,11 @@ dgSystemOfEquations::~dgSystemOfEquations(){
}
void dgSystemOfEquations::saveSolution (std::string name) {
- FILE *f = fopen (name.c_str(),"wb");
- _solution->_data->binarySave(f);
- fclose(f);
+ _solution->save(name);
}
void dgSystemOfEquations::loadSolution (std::string name){
- FILE *f = fopen (name.c_str(),"rb");
- _solution->_data->binaryLoad(f);
- fclose(f);
+ _solution->load(name);
}
void dgSystemOfEquations::export_solution_as_is (const std::string &name){