diff --git a/Solver/TESTCASES/MultirateAdvection.lua b/Solver/TESTCASES/MultirateAdvection.lua
index 34dfde0c54fdda6796c9583fe67055ed5cf3028b..e876a9b2bf5b9f4740ed7afd855bbeb1e96d5a21 100644
--- a/Solver/TESTCASES/MultirateAdvection.lua
+++ b/Solver/TESTCASES/MultirateAdvection.lua
@@ -24,7 +24,7 @@ end
Example of a lua program driving the DG code
--]]
-order = 1
+order = 2
print'*** Loading the mesh and the model ***'
myModel = GModel ()
myModel:load ('rect.geo')
@@ -53,7 +53,7 @@ FS = functionLua(1, 'initial_condition', {'XYZ'}):getName()
GC=dgGroupCollection(myModel,2,order)
solTmp=dgDofContainer(GC,1)
solTmp:L2Projection(FS)
-dt=GC:splitGroupsForMultirate(10,law,solTmp)
+dt=GC:splitGroupsForMultirate(1,law,solTmp)
GC:buildGroupsOfInterfaces(myModel,2,order)
solution=dgDofContainer(GC,1)
solution:L2Projection(FS)
@@ -87,7 +87,7 @@ time=0
-- multirateRK:setLimiter(limiter)
--for i=1,1000
i=0
-while time<0.5 do
+while time<0.1 do
-- norm1 = RK:iterate43SchlegelJCAM2009(law,dt,solution)
-- TEST with Explicit Euler multirate !!!
norm2 = multirateRK:iterate(dt,solution2)
@@ -95,10 +95,10 @@ while time<0.5 do
if (i % 10 == 0) then
print('*** ITER ***',i,time,norm2)
end
- if (i % 10 == 0) then
- solution:exportMsh(string.format("output/rt-%06d", i))
- solution2:exportMsh(string.format("outputMultirate/rt-%06d", i))
- end
+ -- if (i % 10 == 0) then
+ -- solution:exportMsh(string.format("output/rt-%06d", i))
+ -- solution2:exportMsh(string.format("outputMultirate/rt-%06d", i))
+ -- end
i=i+1
end
diff --git a/Solver/dgRungeKuttaMultirate.cpp b/Solver/dgRungeKuttaMultirate.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..b5686bdecb34ad90310fc1c43843c986a171c357
--- /dev/null
+++ b/Solver/dgRungeKuttaMultirate.cpp
@@ -0,0 +1,380 @@
+#include "dgRungeKuttaMultirate.h"
+#include "Bindings.h"
+#include "dgResidual.h"
+#include "dgConservationLaw.h"
+#include "dgDofContainer.h"
+#include "dgGroupOfElements.h"
+#include <algorithm>
+
+dgRungeKuttaMultirate::dgRungeKuttaMultirate(dgGroupCollection* gc,dgConservationLaw *law, int nStages){
+ _law=law;
+ _residualVolume=new dgResidualVolume(*_law);
+ _residualInterface=new dgResidualInterface(*_law);
+ _residualBoundary=new dgResidualBoundary(*_law);
+ _K=new dgDofContainer*[nStages];
+ for(int i=0;i<nStages;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::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);
+ }
+ }
+ }
+}
+
+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>();
+ */
+}
+
+dgRungeKuttaMultirate43::dgRungeKuttaMultirate43(dgGroupCollection *gc,dgConservationLaw* law):dgRungeKuttaMultirate(gc,law,10){}
+
+
+double dgRungeKuttaMultirate43::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 dgRungeKuttaMultirate43::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 dgRungeKuttaMultirate43::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 dgRungeKuttaMultirate43::registerBindings(binding *b) {
+ classBinding *cb = b->addClass<dgRungeKuttaMultirate43>("dgRungeKuttaMultirate43");
+ cb->setDescription("Multirate explicit Runge-Kutta with SchlegelJCAM2009 4 stages 3rd order method");
+ methodBinding *cm;
+ cm = cb->setConstructor<dgRungeKuttaMultirate43,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",&dgRungeKuttaMultirate43::iterate);
+ cm->setArgNames("dt","solution",NULL);
+ cm->setDescription("update the 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");
+}
diff --git a/Solver/dgRungeKuttaMultirate.h b/Solver/dgRungeKuttaMultirate.h
new file mode 100644
index 0000000000000000000000000000000000000000..c65f02690a25411eb5c02b7a5ec9f9639e0a8a7d
--- /dev/null
+++ b/Solver/dgRungeKuttaMultirate.h
@@ -0,0 +1,39 @@
+#include "dgRungeKutta.h"
+
+class dgRungeKuttaMultirate: public dgRungeKutta{
+ protected:
+ 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
+ virtual void computeInputForK(int iK,int exponent,bool isBuffer)=0;
+ void computeK(int iK,int exponent,bool isBuffer);
+ virtual void updateSolution(int exponent,bool isBuffer)=0;
+ void computeResidual(int exponent,bool isBuffer,dgDofContainer *solution, dgDofContainer *residual);
+ public:
+ dgRungeKuttaMultirate(dgGroupCollection *gc,dgConservationLaw* law, int nStages);
+ virtual ~dgRungeKuttaMultirate();
+
+
+ virtual double iterate(double dt, dgDofContainer *solution)=0;
+ static void registerBindings(binding *b);
+};
+
+class dgRungeKuttaMultirate43: public dgRungeKuttaMultirate{
+ void computeInputForK(int iK,int exponent,bool isBuffer);
+ void updateSolution(int exponent,bool isBuffer);
+ public:
+ double iterate(double dt, dgDofContainer *solution);
+ dgRungeKuttaMultirate43(dgGroupCollection *gc,dgConservationLaw *law);
+ static void registerBindings(binding *b);
+};