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

dg : _data is private on dgDofContainer

parent b658a36f
No related branches found
No related tags found
No related merge requests found
Hide whitespace changes
Inline Side-by-side
Solver/dgAlgorithm.cpp
+ 35
− 43
View file @ f505c8a1
......@@ -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)
......
Solver/dgAlgorithm.h
+ 1
− 3
View file @ f505c8a1
......@@ -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,
......
Solver/dgDofContainer.cpp
+ 33
− 2
View file @ f505c8a1
......@@ -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);
}
Solver/dgDofContainer.h
+ 8
− 6
View file @ f505c8a1
......@@ -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
Solver/dgSystemOfEquations.cpp
+ 9
− 14
View file @ f505c8a1
......@@ -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){
......
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment