Skip to content
Snippets Groups Projects
Commit 51ab7705 authored by Richard Comblen's avatar Richard Comblen
Browse files

Splitted multirate rk class

parent 75996ced
No related branches found
No related tags found
No related merge requests found
Hide whitespace changes
Inline Side-by-side
Solver/CMakeLists.txt
+ 1
− 0
View file @ 51ab7705
......@@ -14,6 +14,7 @@ set(SRC
dgGroupOfElements.cpp
dgAlgorithm.cpp
dgRungeKutta.cpp
dgRungeKuttaMultirate.cpp
dgConservationLaw.cpp
dgConservationLawAdvection.cpp
dgSystemOfEquations.cpp
......
Solver/TESTCASES/MultirateAdvection.lua
+ 1
− 1
View file @ 51ab7705
......@@ -78,7 +78,7 @@ LC = 0.1*.1
dt=dt
print('DT=',dt)
RK=dgRungeKutta()
multirateRK=dgRungeKuttaMultirate(GC,law)
multirateRK=dgRungeKuttaMultirate43(GC,law)
norm1=solution:norm()
norm2=solution2:norm()
print('Norm: ',norm1,norm2)
......
Solver/dgRungeKutta.cpp
+ 0
− 354
View file @ 51ab7705
......@@ -5,8 +5,6 @@
#include "dgLimiter.h"
#include "dgResidual.h"
#include "dgGroupOfElements.h"
#include <stdio.h>
#include <algorithm>
double dgRungeKutta::iterateEuler(const dgConservationLaw *claw, double dt, dgDofContainer *solution) {
double A[] = {0};
......@@ -245,355 +243,3 @@ void dgRungeKutta::registerBindings(binding *b) {
cm->setDescription("if a limiter is set, it is applied after each RK stage");
}
dgRungeKuttaMultirate::dgRungeKuttaMultirate(dgGroupCollection* gc,dgConservationLaw *law){
_law=law;
_residualVolume=new dgResidualVolume(*_law);
_residualInterface=new dgResidualInterface(*_law);
_residualBoundary=new dgResidualBoundary(*_law);
_K=new dgDofContainer*[10];
for(int i=0;i<10;i++){
_K[i]=new dgDofContainer(gc,_law->getNbFields());
_K[i]->setAll(0.0);
}
_currentInput=new dgDofContainer(gc,_law->getNbFields());
_currentInput->setAll(0.0);
_residual=new dgDofContainer(gc,_law->getNbFields());
_residual->setAll(0.0);
_minExponent=INT_MAX;
_maxExponent=0;
for(int iGroup=0;iGroup<gc->getNbElementGroups();iGroup++){
dgGroupOfElements *ge=gc->getElementGroup(iGroup);
_minExponent=std::min(_minExponent,ge->getMultirateExponent());
_maxExponent=std::max(_maxExponent,ge->getMultirateExponent());
_innerBufferGroupsOfElements.resize(_maxExponent+1);
_outerBufferGroupsOfElements.resize(_maxExponent+1);
_bulkGroupsOfElements.resize(_maxExponent+1);
if(ge->getIsInnerMultirateBuffer()){
_innerBufferGroupsOfElements[ge->getMultirateExponent()].first.push_back(ge);
}
else if(ge->getIsOuterMultirateBuffer()){
_outerBufferGroupsOfElements[ge->getMultirateExponent()].first.push_back(ge);
}
else{
_bulkGroupsOfElements[ge->getMultirateExponent()].first.push_back(ge);
}
}
// std::vector<std::set<dgGroupOfFaces*> >bulkFacesSet;
// std::vector<std::set<dgGroupOfFaces*> >innerBufferFacesSet;
// std::vector<std::set<dgGroupOfFaces*> >outerBufferFacesSet;
// bulkFacesSet.resize(_maxExponent+1);
// innerBufferFacesSet.resize(_maxExponent+1);
// outerBufferFacesSet.resize(_maxExponent+1);
for(int iGroup=0;iGroup<gc->getNbFaceGroups();iGroup++){
dgGroupOfFaces *gf=gc->getFaceGroup(iGroup);
const dgGroupOfElements *ge[2];
ge[0]=&gf->getGroupLeft();
ge[1]=&gf->getGroupRight();
for(int i=0;i<2;i++){
if(ge[i]->getIsInnerMultirateBuffer()){
_innerBufferGroupsOfElements[ge[i]->getMultirateExponent()].second.push_back(gf);
}
else if(ge[i]->getIsOuterMultirateBuffer()){
_outerBufferGroupsOfElements[ge[i]->getMultirateExponent()].second.push_back(gf);
}else{
_bulkGroupsOfElements[ge[i]->getMultirateExponent()].second.push_back(gf);
}
}
}
for(int iGroup=0;iGroup<gc->getNbBoundaryGroups();iGroup++){
dgGroupOfFaces *gf=gc->getBoundaryGroup(iGroup);
const dgGroupOfElements *ge[1];
ge[0]=&gf->getGroupLeft();
//ge[1]=&gf->getGroupRight();
for(int i=0;i<1;i++){
if(ge[i]->getIsInnerMultirateBuffer()){
_innerBufferGroupsOfElements[ge[i]->getMultirateExponent()].second.push_back(gf);
}
else if(ge[i]->getIsOuterMultirateBuffer()){
_outerBufferGroupsOfElements[ge[i]->getMultirateExponent()].second.push_back(gf);
}else{
_bulkGroupsOfElements[ge[i]->getMultirateExponent()].second.push_back(gf);
}
}
}
// Removing duplicate entries
for(int iExp=0;iExp<=_maxExponent;iExp++){
std::vector<dgGroupOfFaces*>*v[3];
v[0]=&_innerBufferGroupsOfElements[iExp].second;
v[1]=&_outerBufferGroupsOfElements[iExp].second;
v[2]=&_bulkGroupsOfElements[iExp].second;
for(int i=0;i<3;i++){
sort(v[i]->begin(),v[i]->end());
v[i]->erase(unique(v[i]->begin(),v[i]->end()),v[i]->end());
}
}
}
dgRungeKuttaMultirate::~dgRungeKuttaMultirate(){
for(int i=0;i>10;i++){
delete _K[i];
}
delete []_K;
delete _residual;
delete _currentInput;
delete _residualVolume;
delete _residualInterface;
delete _residualBoundary;
}
void dgRungeKuttaMultirate::computeInputForK(int iK,int exponent,bool isBuffer){
if(exponent>_maxExponent){
return;
}
//Msg::Info("Input for K%d at level %d %s",iK,exponent,isBuffer?"Buffer":"Bulk");
double localDt=_dt/pow(2.0,(double)exponent);
double _A[4][4]={
{0, 0, 0, 0},
{1.0/2.0, 0, 0 ,0},
{-1.0/6.0, 2.0/3.0, 0, 0},
{1.0/3.0, -1.0/3.0, 1, 0}
};
double _AInner[10][10]={
{0, 0, 0, 0, 0, 0, 0, 0, 0, 0},
{1.0/4.0 ,0, 0, 0, 0, 0, 0, 0, 0, 0},
{-1.0/12.0, 1.0/3.0, 0, 0, 0, 0, 0, 0, 0, 0},
{1.0/6.0, -1.0/6.0, 1.0/2.0, 0, 0, 0, 0, 0, 0, 0},
{1.0/12.0, 1.0/6.0, 1.0/6.0, 1.0/12.0, 0, 0, 0, 0, 0, 0},
{1.0/12.0, 1.0/6.0, 1.0/6.0, 1.0/12.0, 0, 0, 0, 0, 0, 0},
{1.0/12.0, 1.0/6.0, 1.0/6.0, 1.0/12.0, 0, 1.0/4.0, 0, 0, 0, 0},
{1.0/12.0, 1.0/6.0, 1.0/6.0, 1.0/12.0, 0, -1.0/12.0, 1.0/3.0, 0, 0, 0},
{1.0/12.0, 1.0/6.0, 1.0/6.0, 1.0/12.0, 0, 1.0/6.0, -1.0/6.0, 1.0/2.0, 0, 0},
{1.0/12.0, 1.0/6.0, 1.0/6.0, 1.0/12.0, 0, 1.0/12.0, 1.0/6.0, 1.0/6.0, 1.0/12.0, 0}
};
// Big step RK43
double _AOuter[10][10]={
{0, 0, 0, 0, 0, 0, 0, 0, 0, 0},
{1.0/4.0 , 0, 0, 0, 0, 0, 0, 0, 0, 0},
{1.0/4.0 , 0, 0, 0, 0, 0, 0, 0, 0, 0},
{1.0/2.0 , 0, 0, 0, 0, 0, 0, 0, 0, 0},
{1.0/2.0 , 0, 0, 0, 0, 0, 0, 0, 0, 0},
{-1.0/6.0, 0, 0, 0, 2.0/3.0, 0, 0, 0, 0, 0},
{1.0/12.0, 0, 0, 0, 1.0/6.0, 1.0/2.0, 0, 0, 0, 0},
{1.0/12.0, 0, 0, 0, 1.0/6.0, 1.0/2.0, 0, 0, 0, 0},
{1.0/3.0 , 0, 0, 0, -1.0/3.0, 1, 0, 0, 0, 0},
{1.0/3.0 , 0, 0, 0, -1.0/3.0, 1, 0, 0, 0, 0},
};
// Msg::Info("K%d, %d, %s",iK,exponent,isBuffer?"Buffer":"Bulk");
// compute K[iK] for this exponent and buffer
if(isBuffer){
std::vector<dgGroupOfElements *>&gVi=_innerBufferGroupsOfElements[exponent].first;
std::vector<dgGroupOfElements *>&gVo=_outerBufferGroupsOfElements[exponent].first;
_currentInput->scale(gVi,0.0);
_currentInput->scale(gVo,0.0);
_currentInput->axpy(gVi,*_solution,1.0);
_currentInput->axpy(gVo,*_solution,1.0);
for(int i=0;i<iK;i++){
if(_AInner[iK][i]!=0.0)
_currentInput->axpy(gVi,*_K[i],_AInner[iK][i]*localDt*2);
if(_AOuter[iK][i]!=0.0)
_currentInput->axpy(gVo,*_K[i],_AOuter[iK][i]*localDt*2);
}
}
else{
std::vector<dgGroupOfElements *> &gV=_bulkGroupsOfElements[exponent].first;
_currentInput->scale(gV,0.0);
_currentInput->axpy(gV,*_solution,1.0);
for(int i=0;i<iK;i++){
if(_A[iK][i]!=0.0){
_currentInput->axpy(gV,*_K[i],_A[iK][i]*localDt);
}
}
}
if(!isBuffer){
switch(iK){
case 0:
computeInputForK(0,exponent+1,true);
break;
case 1:
computeInputForK(4,exponent+1,true);
break;
case 2:
computeInputForK(5,exponent+1,true);
break;
case 3:
computeInputForK(9,exponent+1,true);
break;
}
}
else{
if( (iK%5)<4 ){
computeInputForK(iK%5,exponent,false);
}
}
computeK(iK,exponent,isBuffer);
// Msg::Info("Multirate %d %0.16e",iK,_K[iK]->norm());
if( (iK==3 && !isBuffer) || (iK==9 && isBuffer) ){
updateSolution(exponent,isBuffer);
}
if(!isBuffer){
switch(iK){
case 0:
for(int i=1;i<4;i++){
computeInputForK(i,exponent+1,true);
}
break;
case 2:
for(int i=6;i<9;i++){
computeInputForK(i,exponent+1,true);
}
break;
}
}
}
void dgRungeKuttaMultirate::updateSolution(int exponent,bool isBuffer){
//Msg::Info("Updating solution at level %d %s",exponent,isBuffer?"Buffer":"Bulk");
double localDt=_dt/pow(2.0,(double)exponent);
double b[4]={1.0/6.0, 1.0/3.0, 1.0/3.0, 1.0/6.0};
double c[4]={0, 1.0/2.0, 1.0/2.0, 1};
double bInner[10]={1.0/12.0, 1.0/6.0, 1.0/6.0,1.0/12.0,0,1.0/12.0, 1.0/6.0, 1.0/6.0,1.0/12.0,0};
double cInner[10]={0, 1.0/4.0, 1.0/4.0, 1.0/2.0, 1.0/2.0, 1.0/2.0, 3.0/4.0, 3.0/4.0, 1, 1 };
double bOuter[10]={1.0/6.0, 0 , 0 , 0 , 1.0/3.0,1.0/3.0, 0 , 0 , 0 , 1.0/6.0};
double cOuter[10]={0, 1.0/4.0, 1.0/4.0, 1.0/2.0, 1.0/2.0, 1.0/2.0, 3.0/4.0, 3.0/4.0, 1, 1 };
if(isBuffer){
std::vector<dgGroupOfElements *>&gVi=_innerBufferGroupsOfElements[exponent].first;
std::vector<dgGroupOfElements *>&gVo=_outerBufferGroupsOfElements[exponent].first;
for(int i=0;i<10;i++){
if(bInner[i]!=0.0)
_solution->axpy(gVi,*_K[i],bInner[i]*localDt*2);
if(bOuter[i]!=0.0)
_solution->axpy(gVo,*_K[i],bOuter[i]*localDt*2);
}
}
else{
std::vector<dgGroupOfElements *>&gV=_bulkGroupsOfElements[exponent].first;
for(int i=0;i<4;i++){
if(b[i]!=0.0)
_solution->axpy(gV,*_K[i],b[i]*localDt);
}
}
}
void dgRungeKuttaMultirate::computeK(int iK,int exponent,bool isBuffer){
//Msg::Info("Compute K%d at level %d %s",iK,exponent,isBuffer?"Buffer":"Bulk");
if(isBuffer){
std::vector<dgGroupOfElements *>&vei=_innerBufferGroupsOfElements[exponent].first;
std::vector<dgGroupOfElements *>&veo=_outerBufferGroupsOfElements[exponent].first;
_residual->scale(vei,0.0);
_residual->scale(veo,0.0);
for(int i=0;i<vei.size();i++){
_residualVolume->computeAndMap1Group(*vei[i], *_currentInput, *_residual);
}
for(int i=0;i<veo.size();i++){
_residualVolume->computeAndMap1Group(*veo[i], *_currentInput, *_residual);
}
std::vector<dgGroupOfFaces *>& vfi=_innerBufferGroupsOfElements[exponent].second;
std::vector<dgGroupOfFaces *>& vfo=_outerBufferGroupsOfElements[exponent].second;
std::set<dgGroupOfFaces *>gOFSet;
gOFSet.insert(vfo.begin(),vfo.end());
gOFSet.insert(vfi.begin(),vfi.end());
for(std::set<dgGroupOfFaces *>::iterator it=gOFSet.begin();it!=gOFSet.end();it++){
dgGroupOfFaces *faces=*it;
if(faces->isBoundary()){
_residualBoundary->computeAndMap1Group(*faces,*_currentInput,*_residual);
//Msg::Info("Buffer face group %p is boundary in multirate",faces);
}
else{
const dgGroupOfElements *gL=&faces->getGroupLeft();
const dgGroupOfElements *gR=&faces->getGroupRight();
fullMatrix<double> solInterface(faces->getNbNodes(),faces->getNbElements()*2*_law->getNbFields());
fullMatrix<double> residuInterface(faces->getNbNodes(),faces->getNbElements()*2*_law->getNbFields());
faces->mapToInterface(_law->getNbFields(), _currentInput->getGroupProxy(gL), _currentInput->getGroupProxy(gR), solInterface);
_residualInterface->compute1Group(*faces,solInterface,_currentInput->getGroupProxy(gL), _currentInput->getGroupProxy(gR),residuInterface);
//Msg::Info("Buffer face group %p is mapped left or right in multirate",faces);
if(gL->getMultirateExponent()==exponent && gL->getIsMultirateBuffer())
faces->mapLeftFromInterface(_law->getNbFields(), residuInterface,_residual->getGroupProxy(gL));
if(gR->getMultirateExponent()==exponent && gR->getIsMultirateBuffer())
faces->mapRightFromInterface(_law->getNbFields(), residuInterface,_residual->getGroupProxy(gR));
}
}
fullMatrix<double> iMassEl, KEl,residuEl;
for(int iGroup=0;iGroup<vei.size();iGroup++){
dgGroupOfElements *group=vei[iGroup];
int nbNodes = group->getNbNodes();
for(int iElem=0;iElem<group->getNbElements();iElem++){
residuEl.setAsProxy(_residual->getGroupProxy(group),iElem*_law->getNbFields(),_law->getNbFields());
KEl.setAsProxy(_K[iK]->getGroupProxy(group),iElem*_law->getNbFields(),_law->getNbFields());
iMassEl.setAsProxy(group->getInverseMassMatrix(),iElem*nbNodes,nbNodes);
iMassEl.mult(residuEl,KEl);
}
}
for(int iGroup=0;iGroup<veo.size();iGroup++){
dgGroupOfElements *group=veo[iGroup];
int nbNodes = group->getNbNodes();
for(int iElem=0;iElem<group->getNbElements();iElem++){
residuEl.setAsProxy(_residual->getGroupProxy(group),iElem*_law->getNbFields(),_law->getNbFields());
KEl.setAsProxy(_K[iK]->getGroupProxy(group),iElem*_law->getNbFields(),_law->getNbFields());
iMassEl.setAsProxy(group->getInverseMassMatrix(),iElem*nbNodes,nbNodes);
iMassEl.mult(residuEl,KEl);
}
}
}
else{
std::vector<dgGroupOfElements *> &ve=_bulkGroupsOfElements[exponent].first;
_residual->scale(ve,0.0);
for(int i=0;i<ve.size();i++){
_residualVolume->computeAndMap1Group(*ve[i], *_currentInput, *_residual);
}
std::vector<dgGroupOfFaces *> &vf=_bulkGroupsOfElements[exponent].second;
for(int i=0;i<vf.size();i++){
dgGroupOfFaces *faces=vf[i];
if(faces->isBoundary()){
_residualBoundary->computeAndMap1Group(*faces,*_currentInput,*_residual);
}
else{
const dgGroupOfElements *gL=&faces->getGroupLeft();
const dgGroupOfElements *gR=&faces->getGroupRight();
fullMatrix<double> solInterface(faces->getNbNodes(),faces->getNbElements()*2*_law->getNbFields());
fullMatrix<double> residuInterface(faces->getNbNodes(),faces->getNbElements()*2*_law->getNbFields());
faces->mapToInterface(_law->getNbFields(), _currentInput->getGroupProxy(gL), _currentInput->getGroupProxy(gR), solInterface);
_residualInterface->compute1Group(*faces,solInterface,_currentInput->getGroupProxy(gL), _currentInput->getGroupProxy(gR),residuInterface);
if(gL->getMultirateExponent()==exponent && !gL->getIsMultirateBuffer()){
faces->mapLeftFromInterface(_law->getNbFields(), residuInterface,_residual->getGroupProxy(gL));
}
if(gR->getMultirateExponent()==exponent && !gR->getIsMultirateBuffer()){
faces->mapRightFromInterface(_law->getNbFields(), residuInterface,_residual->getGroupProxy(gR));
}
}
}
fullMatrix<double> iMassEl, KEl,residuEl;
for(int iGroup=0;iGroup<ve.size();iGroup++){
dgGroupOfElements *group=ve[iGroup];
int nbNodes = group->getNbNodes();
for(int iElem=0;iElem<group->getNbElements();iElem++){
residuEl.setAsProxy(_residual->getGroupProxy(group),iElem*_law->getNbFields(),_law->getNbFields());
KEl.setAsProxy(_K[iK]->getGroupProxy(group),iElem*_law->getNbFields(),_law->getNbFields());
iMassEl.setAsProxy(group->getInverseMassMatrix(),iElem*nbNodes,nbNodes);
iMassEl.mult(residuEl,KEl);
}
}
}
}
double dgRungeKuttaMultirate::iterate(double dt, dgDofContainer *solution){
_solution=solution;
_dt=dt;
computeInputForK(0,0,false);
computeInputForK(1,0,false);
computeInputForK(2,0,false);
computeInputForK(3,0,false);
return solution->norm();
}
void dgRungeKuttaMultirate::registerBindings(binding *b) {
classBinding *cb = b->addClass<dgRungeKuttaMultirate>("dgRungeKuttaMultirate");
cb->setDescription("Multirate explicit Runge-Kutta");
methodBinding *cm;
cm = cb->setConstructor<dgRungeKuttaMultirate,dgGroupCollection *,dgConservationLaw*>();
cm->setArgNames("groupCollection","law",NULL);
cm->setDescription("A new multirate explicit runge kutta, pass parameters to the iterate function");
cm = cb->addMethod("iterate",&dgRungeKuttaMultirate::iterate);
cm->setArgNames("dt","solution",NULL);
cm->setDescription("update solution by doing a multirate RK43 (from Schlegel et al. Journal of Computational and Applied Mathematics, 2009) step of base time step dt for the conservation law");
cb->setParentClass<dgRungeKutta>();
}
Solver/dgRungeKutta.h
+ 0
− 28
View file @ 51ab7705
......@@ -30,32 +30,4 @@ class dgRungeKutta {
dgRungeKutta();
};
class dgRungeKuttaMultirate: public dgRungeKutta{
private:
int _maxExponent;
int _minExponent;
double _dt;
dgConservationLaw *_law;
dgDofContainer *_solution;
dgDofContainer **_K;
dgDofContainer *_currentInput;
dgDofContainer *_residual;
dgResidualVolume *_residualVolume;
dgResidualInterface *_residualInterface;
dgResidualBoundary *_residualBoundary;
std::vector<std::pair<std::vector<dgGroupOfElements*>,std::vector<dgGroupOfFaces*> > >_bulkGroupsOfElements;// int is the multirateExponent
std::vector<std::pair<std::vector<dgGroupOfElements*>,std::vector<dgGroupOfFaces*> > >_innerBufferGroupsOfElements;// int is the multirateExponent
std::vector<std::pair<std::vector<dgGroupOfElements*>,std::vector<dgGroupOfFaces*> > >_outerBufferGroupsOfElements;// int is the multirateExponent
void computeInputForK(int iK,int exponent,bool isBuffer);
void computeK(int iK,int exponent,bool isBuffer);
void updateSolution(int exponent,bool isBuffer);
void computeResidual(int exponent,bool isBuffer,dgDofContainer *solution, dgDofContainer *residual);
public:
dgRungeKuttaMultirate(dgGroupCollection *gc,dgConservationLaw* law);
~dgRungeKuttaMultirate();
double iterate(double dt, dgDofContainer *solution);
static void registerBindings(binding *b);
};
#endif
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment