Include Conservative Multirate Schemes of Constantinescu, not yet implemented for order>2

Common/LuaBindings.cpp

@@ -361,6 +361,7 @@ binding::binding(){
   dgRungeKutta::registerBindings(this);
   dgRungeKuttaMultirate43::registerBindings(this);
   dgRungeKuttaMultirate22::registerBindings(this);
+  dgRungeKuttaMultirateConservative::registerBindings(this);
   dgSlopeLimiterRegisterBindings(this);
   dgSupraTransformNodalValueRegisterBindings(this);
   dgSystemOfEquations::registerBindings(this);

Solver/TESTCASES/ConservationTestMRAdvection.lua

@@ -65,86 +65,140 @@ GC=dgGroupCollection(myModel,2,order)
 solTmp=dgDofContainer(GC,1)
 solTmp:L2Projection(FS)
 dt=GC:splitGroupsForMultirate(2,1,2,law,solTmp)
+
+-- Example of Conservative Multirate scheme with Butcher tables as input (here RK2a)
+A=fullMatrix(2, 2)
+A:set(0, 0, 0)
+A:set(0, 1, 0)
+A:set(1, 0, 1)
+A:set(1, 1, 0)
+
+b=fullMatrix(1, 2)
+b:set(0, 0, 0.5)
+b:set(0, 1, 0.5)
+
+c=fullMatrix(1, 2)
+c:set(0, 0, 0)
+c:set(0, 1, 1)
+
+-- Building the RK objects
+
+multirateRK=dgRungeKuttaMultirateConservative.new2b(GC, law)
+multirateRK:printButcher()
+-- RK2a
+RK2a=dgRungeKutta()
+
+-- RK43
+RK43=dgRungeKutta()
+
+-- Multirate RK2a
+
+-- *** Default Multirate RK2a ***
+--multirateRK2a=dgRungeKuttaMultirate22(GC,law)
+
+-- *** Multirate RK2a with as input the butcher tables ***
+multirateRK2a=dgRungeKuttaMultirateConservative(GC, law, A, b, c)
+
+-- Multirate RK43
+multirateRK43=dgRungeKuttaMultirate43(GC,law)
+
+-- Multirate RK2b 
+multirateRK2b=dgRungeKuttaMultirateConservative.new2b(GC, law)
+-- multirateRK2b:printButcher()
+
 GC:buildGroupsOfInterfaces(myModel,2,order)
-solution=dgDofContainer(GC,1)
-solution:L2Projection(FS)
 solution0=dgDofContainer(GC,1)
+solution1=dgDofContainer(GC,1)
 solution2=dgDofContainer(GC,1)
 solution3=dgDofContainer(GC,1)
+solution4=dgDofContainer(GC,1)
 solution0:L2Projection(FS)
+solution1:L2Projection(FS)
 solution2:L2Projection(FS)
 solution3:L2Projection(FS)
-solution:exportGroupIdMsh()
+solution4:L2Projection(FS)
+solution0:exportGroupIdMsh()
 
 print'---- Groups splitted ----'
 
 print'*** setting the initial solution ***'
 
-solution:exportMsh(string.format("output/adv_sr43-%06d", 0)) 
-solution0:exportMsh(string.format("output/adv_sr22-%06d", 0)) 
-solution2:exportMsh(string.format("outputMultirate/adv_mr43-%06d", 0)) 
-solution3:exportMsh(string.format("outputMultirate/adv_mr22-%06d", 0))
+solution0:exportMsh(string.format("outputSinglerate/adv_sr2a-%06d", 0)) 
+solution1:exportMsh(string.format("outputSinglerate/adv_sr43-%06d", 0)) 
+solution2:exportMsh(string.format("outputMultirate/adv_mr2a-%06d", 0)) 
+solution3:exportMsh(string.format("outputMultirate/adv_mr43-%06d", 0))
+solution4:exportMsh(string.format("outputMultirate/adv_mr2b-%06d", 0))
 print'*** solve ***'
 
 LC = 0.1*.1
 dt=dt
 print('DT=',dt)
-RK22=dgRungeKutta()
-RK=dgRungeKutta()
-multirateRK=dgRungeKuttaMultirate43(GC,law)
-multirateRK22=dgRungeKuttaMultirate22(GC,law)
-norm1=solution:norm()
-norm0=solution:norm()
+
+norm0=solution0:norm()
+norm1=solution1:norm()
 norm2=solution2:norm()
 norm3=solution3:norm()
-print('Norm: ',norm0,norm1,norm2,norm3)
+norm4=solution4:norm()
+
+print('Norm: ', norm0, norm1, norm2, norm3, norm4)
 dt=dt
 time=0
 i=0
 integrator=dgFunctionIntegrator(functionLua(1, 'toIntegrate', {functionSolution.get()}))
        print('*** ITER-START ***')
-       v=fullMatrix(1,1);
-			 w=fullMatrix(1,1);
-			 u=fullMatrix(1,1);
-			 z=fullMatrix(1,1);
-       integrator:compute(solution,v)
-       integrator:compute(solution0,z)
-			 integrator:compute(solution2,w)
-			 integrator:compute(solution3,u)
-			 uref=u:get(0,0,0)
-			 vref=v:get(0,0,0)
-			 wref=w:get(0,0,0)
-			 zref=z:get(0,0,0)
+       w0=fullMatrix(1,1);
+			 w1=fullMatrix(1,1);
+			 w2=fullMatrix(1,1);
+			 w3=fullMatrix(1,1);
+			 w4=fullMatrix(1,1);
+       integrator:compute(solution0,w0)
+       integrator:compute(solution1,w1)
+			 integrator:compute(solution2,w2)
+			 integrator:compute(solution3,w3)
+			 integrator:compute(solution4,w4)
+			 w0ref=w0:get(0,0,0)
+			 w1ref=w1:get(0,0,0)
+			 w2ref=w2:get(0,0,0)
+			 w3ref=w3:get(0,0,0)
+			 w4ref=w4:get(0,0,0)
 			 print'*****************************************'
-       print('integral (SR43) : ',v:get(0,0,0))
-       print('integral (MR43) : ',w:get(0,0,0))
-       print('integral (SR22) : ',z:get(0,0,0))
-       print('integral (MR22) : ',u:get(0,0,0))
+       print('integral (SR2a) : ',w0:get(0,0,0))
+       print('integral (SR43) : ',w1:get(0,0,0))
+       print('integral (MR2a) : ',w2:get(0,0,0))
+       print('integral (MR43) : ',w3:get(0,0,0))
+       print('integral (MR2b) : ',w4:get(0,0,0))
 			 print'******************************************'
 while i<10 do
-		norm0=RK22:iterate22(law, dt, solution0)
-    norm1 = RK:iterate43SchlegelJCAM2009(law,dt,solution)
-    norm2 = multirateRK:iterate(dt,solution2)
-		norm3 = multirateRK22:iterate(dt, solution3)
+		norm0 = RK2a:iterate22(law, dt, solution0)
+    norm1 = RK43:iterate43SchlegelJCAM2009(law, dt, solution1)
+    norm2 = multirateRK2a:iterate(dt, solution2)
+		norm3 = multirateRK43:iterate(dt, solution3)
+		norm4 = multirateRK2b:iterate(dt, solution4)
     time=time+dt
     if (i % 1 == 0) then 
-       print('*** ITER ***',i)--,time,norm1, norm2,norm0, norm3)
+       print('*** ITER ***',i)--,time,norm0, norm1,norm2, norm3)
 			 print'-----------------------------------------'
-			 integrator:compute(solution,v)
-			 integrator:compute(solution2,w)
-			 integrator:compute(solution3, u)
-			 integrator:compute(solution0, z)
-       print('integral (SR43) : ',v:get(0,0,0))
-       print('integral (MR43) : ',w:get(0,0,0))
-       print('integral (SR22) : ',z:get(0,0,0))
-       print('integral (MR22) : ',u:get(0,0,0))
-			 print(string.format("Mass conservation relative error:\n SR43=%0.16e,\n MR43=%0.16e,\n SR22=%0.16e,\n MR22=%0.16e",math.abs((vref-v:get(0,0,0))/vref), math.abs((wref-w:get(0,0,0))/wref),math.abs((zref-z:get(0,0,0))/zref),math.abs((uref-u:get(0,0,0))/uref)));
+			 integrator:compute(solution0, w0)
+			 integrator:compute(solution1, w1)
+			 integrator:compute(solution2, w2)
+			 integrator:compute(solution3, w3)
+			 integrator:compute(solution4, w4)
+
+       print('integral (SR2a) : ',w0:get(0,0,0))
+       print('integral (SR43) : ',w1:get(0,0,0))
+       print('integral (MR2a) : ',w2:get(0,0,0))
+       print('integral (MR43) : ',w3:get(0,0,0))
+       print('integral (MR2b) : ',w4:get(0,0,0))
+			 
+			 print(string.format("Mass conservation relative error:\n SR2a=%0.16e,\n SR43=%0.16e,\n MR2a=%0.16e,\n MR43=%0.16e,\n MR2b=%0.16e",math.abs((w0ref-w0:get(0,0,0))/w0ref), math.abs((w1ref-w1:get(0,0,0))/w1ref),math.abs((w2ref-w2:get(0,0,0))/w2ref),math.abs((w3ref-w3:get(0,0,0))/w3ref), math.abs((w4ref-w4:get(0,0,0))/w4ref)));
 			 print'-----------------------------------------'
     end
     if (i % 10 == 0) then 
-       solution:exportMsh(string.format("output/adv_sr43-%06d", i)) 
-       solution2:exportMsh(string.format("outputMultirate/adv_mr43-%06d", i)) 
-       solution3:exportMsh(string.format("outputMultirate/adv_mr22-%06d", i)) 
+       solution0:exportMsh(string.format("outputSinglerate/adv_sr2a-%06d", i)) 
+       solution1:exportMsh(string.format("outputSinglerate/adv_sr43-%06d", i)) 
+       solution2:exportMsh(string.format("outputMultirate/adv_mr2a-%06d", i)) 
+       solution3:exportMsh(string.format("outputMultirate/adv_mr43-%06d", i)) 
+       solution4:exportMsh(string.format("outputMultirate/adv_mr2b-%06d", i)) 
     end
     i=i+1
 end

Solver/dgRungeKuttaMultirate.cpp

@@ -77,24 +77,28 @@ static double cInner22[4]={0,  1.0/2.0, 1.0/2.0, 1.0};
 static double bOuter22[4]={1.0/4.0, 1.0/4.0, 1.0/4.0, 1.0/4.0};
 static double cOuter22[4]={0, 1.0, 0, 1.0};
 
-dgRungeKuttaMultirate::dgRungeKuttaMultirate(dgGroupCollection* gc,dgConservationLaw *law, int nStages){
+dgRungeKuttaMultirate::dgRungeKuttaMultirate(dgGroupCollection* gc,dgConservationLaw *law){
   _law=law;
   _residualVolume=new dgResidualVolume(*_law);
   _residualInterface=new dgResidualInterface(*_law);
   _gc=gc;
+	_init=false;
+}
+void dgRungeKuttaMultirate::initialize(int nStages){
+	if(_init==false){
   _K=new dgDofContainer*[nStages];
   for(int i=0;i<nStages;i++){
-    _K[i]=new dgDofContainer(gc,_law->getNbFields());
+    _K[i]=new dgDofContainer(_gc,_law->getNbFields());
     _K[i]->setAll(0.0);
   }
-  _currentInput=new dgDofContainer(gc,_law->getNbFields());
+  _currentInput=new dgDofContainer(_gc,_law->getNbFields());
   _currentInput->setAll(0.0);
-  _residual=new dgDofContainer(gc,_law->getNbFields());
+  _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);
+  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);
@@ -110,8 +114,8 @@ dgRungeKuttaMultirate::dgRungeKuttaMultirate(dgGroupCollection* gc,dgConservatio
       _bulkGroupsOfElements[ge->getMultirateExponent()].first.push_back(ge);
     }
   }
-  for(int iGroup=0;iGroup<gc->getNbFaceGroups();iGroup++){
-    dgGroupOfFaces *gf=gc->getFaceGroup(iGroup);
+  for(int iGroup=0;iGroup<_gc->getNbFaceGroups();iGroup++){
+    dgGroupOfFaces *gf=_gc->getFaceGroup(iGroup);
     for(int i=0;i<gf->getNbGroupOfConnections();i++){
       const dgGroupOfElements *ge = &gf->getGroupOfConnections(i).getGroupOfElements();
       if(ge->getIsInnerMultirateBuffer()){
@@ -135,6 +139,11 @@ dgRungeKuttaMultirate::dgRungeKuttaMultirate(dgGroupCollection* gc,dgConservatio
       v[i]->erase(unique(v[i]->begin(),v[i]->end()),v[i]->end());
     }
   }
+	_init=true;
+	}
+	else
+		return;
+
 }
 dgRungeKuttaMultirate::~dgRungeKuttaMultirate(){
   for(int i=0;i>10;i++){
@@ -273,10 +282,11 @@ void dgRungeKuttaMultirate::registerBindings(binding *b) {
   */
 }
 
-dgRungeKuttaMultirate43::dgRungeKuttaMultirate43(dgGroupCollection *gc,dgConservationLaw* law):dgRungeKuttaMultirate(gc,law,10){}
+dgRungeKuttaMultirate43::dgRungeKuttaMultirate43(dgGroupCollection *gc,dgConservationLaw* law): dgRungeKuttaMultirate(gc,law){}
 
 
 double dgRungeKuttaMultirate43::iterate(double dt, dgDofContainer *solution){
+  initialize(10);
   _solution=solution;
   _dt=dt;
   computeInputForK(0,0,false);
@@ -460,10 +470,11 @@ void dgRungeKuttaMultirate43::registerBindings(binding *b) {
   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");
 }
 
-dgRungeKuttaMultirate22::dgRungeKuttaMultirate22(dgGroupCollection *gc,dgConservationLaw* law):dgRungeKuttaMultirate(gc,law,4){}
+dgRungeKuttaMultirate22::dgRungeKuttaMultirate22(dgGroupCollection *gc,dgConservationLaw* law): dgRungeKuttaMultirate(gc,law){}
 
 
 double dgRungeKuttaMultirate22::iterate(double dt, dgDofContainer *solution){
+	initialize(4);
 	_solution=solution;
 	_dt=dt;
 	computeInputForK(0,0,false,-1);
@@ -620,3 +631,329 @@ void dgRungeKuttaMultirate22::registerBindings(binding *b) {
 	cm->setArgNames("dt","solution",NULL);
 	cm->setDescription("update the solution by doing a multirate RK2a (from Constantinescu and Sandu,  'Update on Multirate Timestepping Methods for Hyperbolic Conservation Laws', Computer Sciance Technical Report,  2007) step of base time step dt for the conservation law");
 }
+
+dgRungeKuttaMultirateConservative::dgRungeKuttaMultirateConservative(dgGroupCollection *gc, dgConservationLaw *law,  fullMatrix<double> *A,  fullMatrix<double> *b, fullMatrix<double> *c):  _A(A), _b(b), _c(c), dgRungeKuttaMultirate(gc, law){
+	if(A->size1()!=A->size2() || A->size1()!=b->size2() || A->size1()!=c->size2())
+		Msg::Error("The base method's Butcher tables have incompatible sizes");
+	//Number of stages of the base method
+	int bStages = A->size1();
+	
+	//Compute Butcher Table for Outer Buffer
+	_AOuter=new fullMatrix<double> (2*bStages,2*bStages);
+	_bOuter=new fullMatrix<double> (1,2*bStages);
+	_cOuter=new fullMatrix<double> (1,2*bStages);
+
+	//Compute Butcher Table for Inner Buffer
+	_AInner=new fullMatrix<double> (2*bStages,2*bStages);
+	_bInner=new fullMatrix<double> (1,2*bStages);
+	_cInner=new fullMatrix<double> (1,2*bStages);
+	for(int i=0;i<bStages;i++){
+		for(int j=0;j<bStages;j++){
+			
+			(*_AOuter)(i, j)=(*_A)(i, j);
+			(*_AOuter)(bStages+i,bStages+j)=(*_A)(i, j);
+			(*_bOuter)(0, i)=0.5*(*_b)(0, i);
+			(*_bOuter)(0,bStages+i)=0.5*(*_b)(0, i);
+			(*_cOuter)(0, i)=(*_c)(0, i);
+			(*_cOuter)(0,bStages+i)=(*_c)(0, i);
+			
+			(*_AInner)(i, j)=0.5*(*_A)(i, j);
+			(*_AInner)(bStages+i,bStages+j)=0.5*(*_A)(i, j);
+			(*_AInner)(bStages+i,j)=0.5*(*_b)(0, j);
+			(*_bInner)(0, i)=0.5*(*_b)(0, i);
+			(*_bInner)(0,bStages+i)=0.5*(*_b)(0, i);
+			(*_cInner)(0, i)=0.5*(*_c)(0, i);
+			(*_cInner)(0,bStages+i)=0.5+0.5*(*_c)(0, i);
+
+		}
+	}
+}
+void dgRungeKuttaMultirateConservative::printButcher(){
+	Msg::Info("------- Printing Inner Butcher Tables --------");
+	printf("\n ************************************************** \n");
+	_AInner->print("A_inner");
+	printf("\n ********************* \n");
+	_bInner->print("b_inner");
+	printf("\n ********************* \n");
+	_cInner->print("c_inner");
+	Msg::Info("------- Printing Outer Butcher Tables --------");
+	printf("\n ************************************************** \n");
+	_AOuter->print("A_outer");
+	printf("\n ********************* \n");
+	_bOuter->print("b_outer");
+	printf("\n ********************* \n");
+	_cOuter->print("c_outer");
+
+
+}
+double dgRungeKuttaMultirateConservative::iterate(double dt, dgDofContainer *solution){
+	initialize(4);
+	_solution=solution;
+	_dt=dt;
+	computeInputForK(0,0,false,-1);
+	computeInputForK(1,0,false,-1);
+	
+	return solution->norm();
+}
+
+void dgRungeKuttaMultirateConservative::computeInputForK(int iK,int exponent,bool isBuffer,int upperLeveliK){
+	if(exponent>_maxExponent){
+		return;
+	}
+#ifdef MULTIRATEVERBOSE
+  for(int i=0;i<exponent;i++)
+    printf("\t");
+  printf("Exponent %d, %s, input   K%d\n",exponent,isBuffer?"    Buffer":"Not buffer",iK);
+#endif
+	double localDt=_dt/pow(2.0,(double)exponent);
+	// compute K[iK] for this exponent and buffer
+	
+	if(!isBuffer){
+		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);
+			}
+		}
+	}
+	else{
+		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);
+		}
+    // We need to update input for the neighboring elements with bigger time step.
+    // if there is no corresponding K to be computed
+    if( (iK>0 && iK<3) ){
+      std::vector<dgGroupOfElements *>&gVbigger=_bulkGroupsOfElements[exponent-1].first;
+      _currentInput->scale(gVbigger,0.0);
+      _currentInput->axpy(gVbigger,*_solution,1.0);
+      int idOuter2Bulk[10]={0,0,0,1};
+      for(int i=0;i<iK;i++){
+        if((*_AOuter)(iK, i)!=0.0){
+          _currentInput->axpy(gVbigger,*_K[idOuter2Bulk[i]],(*_AOuter)(iK, i)*localDt*2);
+        }
+      }
+      std::vector<dgGroupOfElements *>&gViBigger=_innerBufferGroupsOfElements[exponent-1].first;
+      _currentInput->scale(gViBigger,0.0);
+      _currentInput->axpy(gViBigger,*_solution,1.0);
+      // shift
+      int s=0;
+      if(upperLeveliK>=2){
+        for(int i=0;i<2;i++){
+          _currentInput->axpy(gViBigger,*_K[i],(*_b)(0, i)*localDt*2);
+        }
+        s+=2;
+      }
+      for(int i=0;i<iK;i++){
+        if((*_AOuter)(iK, i)!=0.0){
+          _currentInput->axpy(gViBigger,*_K[idOuter2Bulk[i]+s],(*_AOuter)(iK, i)*localDt*2);
+        }
+      }
+    }
+	}
+	
+	if(!isBuffer){
+		switch(iK){
+			case 0:
+				computeInputForK(0,exponent+1,true,iK);
+				break;
+			case 1:
+				computeInputForK(3,exponent+1,true,iK);
+				break;
+		}
+	}
+	else{
+		computeInputForK(iK%2,exponent,false,iK);
+	}
+	
+	if(exponent==0){
+		computeK(iK, exponent, false);
+		switch(iK){
+			case 0:
+				for(int i=1;i<3;i++){
+					computeInputForK(i,exponent+1,true,iK);
+				}
+				break;
+			case 1:
+				updateSolution(exponent, false);
+				break;
+		}
+	}
+
+	if(isBuffer && exponent>0){
+		computeK(iK%2, exponent, false);
+		computeK(iK, exponent, true);
+		if(iK%2==1)
+			updateSolution(exponent, false);
+		if(iK==3)
+			updateSolution(exponent, true);
+    switch(iK%2){
+      case 0:
+        for(int i=1;i<3;i++){
+          computeInputForK(i, exponent+1, true,iK);
+        }
+        break;
+    }
+ }
+}
+
+void dgRungeKuttaMultirateConservative::updateSolution(int exponent,bool isBuffer){
+#ifdef MULTIRATEVERBOSE
+  for(int i=0;i<exponent;i++)
+    printf("\t");
+  printf("Updating solution at level %d %s\n",exponent,isBuffer?"Buffer":"Bulk");
+#endif
+	double localDt=_dt/pow(2.0,(double)exponent);
+	if(isBuffer){
+		std::vector<dgGroupOfElements *>&gVi=_innerBufferGroupsOfElements[exponent].first;
+		std::vector<dgGroupOfElements *>&gVo=_outerBufferGroupsOfElements[exponent].first;
+		for(int i=0;i<4;i++){
+			if((*_bInner)(0, i)!=0.0)
+				_solution->axpy(gVi,*_K[i],(*_bInner)(0, i)*localDt*2);
+			if((*_bOuter)(0, i)!=0.0)
+				_solution->axpy(gVo,*_K[i],(*_bOuter)(0, i)*localDt*2);
+		}
+	}
+	else{
+		std::vector<dgGroupOfElements *>&gV=_bulkGroupsOfElements[exponent].first;
+		for(int i=0;i<2;i++){
+			if((*_b)(0, i)!=0.0)
+				_solution->axpy(gV,*_K[i],(*_b)(0, i)*localDt);
+		}
+	}
+}
+
+dgRungeKuttaMultirateConservative* dgRungeKuttaMultirateConservative::new44(dgGroupCollection *gc,  dgConservationLaw *law){
+	Msg::Error("---- Generic Multirate of Constantinescu doesn't work yet for this standard Butcher tables --------");
+	fullMatrix<double> *A=new fullMatrix<double>(4, 4);
+	fullMatrix<double> *b=new fullMatrix<double>(1, 4);
+	fullMatrix<double> *c=new fullMatrix<double>(1, 4);
+
+  double At[4][4]={
+		{0, 0, 0, 0}, 
+		{1.0/2.0, 0, 0, 0}, 
+		{0, 1.0/2.0, 0, 0}, 
+		{0, 0, 1, 0}
+	};
+  double bt[4]={1.0/6.0, 1.0/3.0, 1.0/3.0, 1.0/6.0};
+	double ct[4]={0, 1.0/2.0, 1.0/2.0, 1};
+
+	for(int i=0; i<4; i++){
+		for (int j=0; j<4; j++){
+			(*A)(i, j)=At[i][j];
+		}
+		(*b)(0, i)=bt[i];
+		(*c)(0, i)=ct[i];
+	}
+
+	return new dgRungeKuttaMultirateConservative(gc, law, A, b, c);
+}
+dgRungeKuttaMultirateConservative* dgRungeKuttaMultirateConservative::new43(dgGroupCollection *gc,  dgConservationLaw *law){
+	Msg::Error("---- Generic Multirate of Constantinescu doesn't work yet for this standard Butcher tables --------");
+	fullMatrix<double> *A=new fullMatrix<double>(4, 4);
+	fullMatrix<double> *b=new fullMatrix<double>(1, 4);
+	fullMatrix<double> *c=new fullMatrix<double>(1, 4);
+
+  double At[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 bt[4]={1.0/6.0, 1.0/3.0, 1.0/3.0, 1.0/6.0};
+	double ct[4]={0, 1.0/2.0, 1.0/2.0, 1};
+
+	for(int i=0; i<4; i++){
+		for (int j=0; j<4; j++){
+			(*A)(i, j)=At[i][j];
+		}
+		(*b)(0, i)=bt[i];
+		(*c)(0, i)=ct[i];
+	}
+
+	return new dgRungeKuttaMultirateConservative(gc, law, A, b, c);
+}
+dgRungeKuttaMultirateConservative* dgRungeKuttaMultirateConservative::new2a(dgGroupCollection *gc,  dgConservationLaw *law){
+	
+	fullMatrix<double> *A=new fullMatrix<double>(2, 2);
+	fullMatrix<double> *b=new fullMatrix<double>(1, 2);
+	fullMatrix<double> *c=new fullMatrix<double>(1, 2);
+
+  double At[2][2]={
+		{0, 0}, 
+		{1, 0}
+	};
+  double bt[2]={0.5, 0.5};
+	double ct[2]={0, 1};
+
+	for(int i=0; i<2; i++){
+		for (int j=0; j<2; j++){
+			(*A)(i, j)=At[i][j];
+		}
+		(*b)(0, i)=bt[i];
+		(*c)(0, i)=ct[i];
+	}
+
+	return new dgRungeKuttaMultirateConservative(gc, law, A, b, c);
+}
+dgRungeKuttaMultirateConservative* dgRungeKuttaMultirateConservative::new2b(dgGroupCollection *gc,  dgConservationLaw *law){
+	
+	fullMatrix<double> *A=new fullMatrix<double>(2, 2);
+	fullMatrix<double> *b=new fullMatrix<double>(1, 2);
+	fullMatrix<double> *c=new fullMatrix<double>(1, 2);
+
+  double At[2][2]={
+		{0, 0}, 
+		{1.0/2.0, 0}
+	};
+  double bt[2]={0, 1};
+	double ct[2]={0, 1.0/2.0};
+
+	for(int i=0; i<2; i++){
+		for (int j=0; j<2; j++){
+			(*A)(i, j)=At[i][j];
+		}
+		(*b)(0, i)=bt[i];
+		(*c)(0, i)=ct[i];
+	}
+
+	return new dgRungeKuttaMultirateConservative(gc, law, A, b, c);
+}
+
+void dgRungeKuttaMultirateConservative::registerBindings(binding *b) {
+	classBinding *cb = b->addClass<dgRungeKuttaMultirateConservative>("dgRungeKuttaMultirateConservative");
+	cb->setDescription("Multirate explicit Runge-Kutta with Constantinescu n (base method) stages 2nd order method");
+	methodBinding *cm;
+	cm = cb->setConstructor<dgRungeKuttaMultirateConservative,dgGroupCollection*,dgConservationLaw* , fullMatrix<double>*, fullMatrix<double>*, fullMatrix<double>* >();
+	cm->setArgNames("groupCollection","law","A","b","c",NULL);
+	cm->setDescription("A new multirate explicit runge kutta, pass parameters to the iterate function");
+	cm = cb->addMethod("iterate",&dgRungeKuttaMultirateConservative::iterate);
+	cm->setArgNames("dt","solution",NULL);
+	cm->setDescription("update the solution by doing a multirate RK constructed from the given Butcher tables (from Constantinescu and Sandu, 'Update on Multirate Timestepping Methods for Hyperbolic Conservation Laws', Computer Sciance Technical Report,  2007) step of base time step dt for the conservation law");
+  cm = cb->addMethod("new44", &dgRungeKuttaMultirateConservative::new43);
+	cm->setDescription("Creates a new Conservative Runge-Kutta scheme based on the base method RK44");
+	cm->setArgNames("groupCollection", "law", NULL);
+  cm = cb->addMethod("new43", &dgRungeKuttaMultirateConservative::new43);
+	cm->setDescription("Creates a new Conservative Runge-Kutta scheme based on the base method RK43");
+	cm->setArgNames("groupCollection", "law", NULL);
+  cm = cb->addMethod("new2a", &dgRungeKuttaMultirateConservative::new2a);
+	cm->setDescription("Creates a new Conservative Runge-Kutta scheme based on the base method RK2a");
+	cm->setArgNames("groupCollection", "law", NULL);
+  cm = cb->addMethod("new2b", &dgRungeKuttaMultirateConservative::new2b);
+	cm->setDescription("Creates a new Conservative Runge-Kutta scheme based on the base method RK2b");
+	cm->setArgNames("groupCollection", "law", NULL);
+  cm = cb->addMethod("printButcher", &dgRungeKuttaMultirateConservative::printButcher);
+	cm->setDescription("Print the Butcher Tables for Buffers");
+	cm->setArgNames(NULL);
+}
+

Solver/dgRungeKuttaMultirate.h

@@ -2,6 +2,7 @@
 
 class dgRungeKuttaMultirate: public dgRungeKutta{
   protected:
+  bool _init;	
   int _maxExponent;
   int _minExponent;
   double _dt;
@@ -16,14 +17,13 @@ class dgRungeKuttaMultirate: public dgRungeKutta{
   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 initialize(int nStages);
   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);
+  dgRungeKuttaMultirate(dgGroupCollection *gc,dgConservationLaw* law);
   virtual ~dgRungeKuttaMultirate();
-
-
   virtual double iterate(double dt, dgDofContainer *solution)=0;
   static void registerBindings(binding *b);
 };
@@ -45,3 +45,27 @@ class dgRungeKuttaMultirate22: public dgRungeKuttaMultirate{
   dgRungeKuttaMultirate22(dgGroupCollection *gc,dgConservationLaw *law);
   static void registerBindings(binding *b);
 };
+class dgRungeKuttaMultirateConservative: public dgRungeKuttaMultirate{
+	private:
+	fullMatrix<double> *_A;
+	fullMatrix<double> *_b;
+	fullMatrix<double> *_c;
+	fullMatrix<double> *_AOuter;
+	fullMatrix<double> *_bOuter;
+	fullMatrix<double> *_cOuter;
+	fullMatrix<double> *_AInner;
+	fullMatrix<double> *_bInner;
+	fullMatrix<double> *_cInner;
+
+	void computeInputForK(int iK,int exponent,bool isBuffer,int upperLeveliK);
+	void updateSolution(int exponent,bool isBuffer);
+	public:
+	double iterate(double dt, dgDofContainer *solution);
+	dgRungeKuttaMultirateConservative(dgGroupCollection *gc,dgConservationLaw *law, fullMatrix<double> *A, fullMatrix<double> *b,fullMatrix<double> *c);
+	void printButcher();
+	static dgRungeKuttaMultirateConservative* new44(dgGroupCollection *gc, dgConservationLaw *law);
+	static dgRungeKuttaMultirateConservative* new43(dgGroupCollection *gc, dgConservationLaw *law);
+	static dgRungeKuttaMultirateConservative* new2a(dgGroupCollection *gc, dgConservationLaw *law);
+	static dgRungeKuttaMultirateConservative* new2b(dgGroupCollection *gc, dgConservationLaw *law);
+	static void registerBindings(binding *b);
+};