DG: Conservative multirate with optimal order, conservation needs to be verified

Solver/TESTCASES/ConvergenceMultirate.lua

@@ -62,7 +62,7 @@ FS = functionLua(1, 'initial_condition', {'XYZ'}):getName()
 GC=dgGroupCollection(myModel,2,order)
 solTmp=dgDofContainer(GC,1)
 solTmp:L2Projection(FS)
-dt=GC:splitGroupsForMultirate(20,1,law,solTmp)
+dt=GC:splitGroupsForMultirate(100,1,law,solTmp)
 GC:buildGroupsOfInterfaces(myModel,2,order)
 solution=dgDofContainer(GC,1)
 solutionRef=dgDofContainer(GC,1)
@@ -82,8 +82,8 @@ for i=1,nRefSteps do
   end
 end
 
-multirateRK=dgRungeKuttaMultirate43(GC,law)
---multirateRK=dgRungeKuttaMultirate22(GC,law)
+--multirateRK=dgRungeKuttaMultirate43(GC,law)
+multirateRK=dgRungeKuttaMultirate22(GC,law)
 integrator=dgFunctionIntegrator(functionLua(1, 'L2Norm', {'Solution'}):getName())
 
 oldnorm=1

Solver/dgRungeKuttaMultirate.cpp

@@ -479,16 +479,21 @@ dgRungeKuttaMultirate22::dgRungeKuttaMultirate22(dgGroupCollection *gc,dgConserv
 double dgRungeKuttaMultirate22::iterate(double dt, dgDofContainer *solution){
 	_solution=solution;
 	_dt=dt;
-	computeInputForK(0,0,false);
-	computeInputForK(1,0,false);
+	computeInputForK(0,0,false,-1);
+	computeInputForK(1,0,false,-1);
 	
 	return solution->norm();
 }
 
-void dgRungeKuttaMultirate22::computeInputForK(int iK,int exponent,bool isBuffer){
+void dgRungeKuttaMultirate22::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
 	
@@ -515,20 +520,49 @@ void dgRungeKuttaMultirate22::computeInputForK(int iK,int exponent,bool isBuffer
 			if(AOuter22[iK][i]!=0.0)
 				_currentInput->axpy(gVo,*_K[i],AOuter22[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(AOuter22[iK][i]!=0.0){
+          _currentInput->axpy(gVbigger,*_K[idOuter2Bulk[i]],AOuter22[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],b22[i]*localDt*2);
+        }
+        s+=2;
+      }
+      for(int i=0;i<iK;i++){
+        if(AOuter22[iK][i]!=0.0){
+          _currentInput->axpy(gViBigger,*_K[idOuter2Bulk[i]+s],AOuter22[iK][i]*localDt*2);
+        }
+      }
+    }
 	}
 	
 	if(!isBuffer){
 		switch(iK){
 			case 0:
-				computeInputForK(0,exponent+1,true);
+				computeInputForK(0,exponent+1,true,iK);
 				break;
 			case 1:
-				computeInputForK(3,exponent+1,true);
+				computeInputForK(3,exponent+1,true,iK);
 				break;
 		}
 	}
 	else{
-		computeInputForK(iK%2,exponent,false);
+		computeInputForK(iK%2,exponent,false,iK);
 	}
 	
 	if(exponent==0){
@@ -536,7 +570,7 @@ void dgRungeKuttaMultirate22::computeInputForK(int iK,int exponent,bool isBuffer
 		switch(iK){
 			case 0:
 				for(int i=1;i<3;i++){
-					computeInputForK(i,exponent+1,true);
+					computeInputForK(i,exponent+1,true,iK);
 				}
 				break;
 			case 1:
@@ -552,10 +586,10 @@ void dgRungeKuttaMultirate22::computeInputForK(int iK,int exponent,bool isBuffer
 			updateSolution(exponent, false);
 		if(iK==3)
 			updateSolution(exponent, true);
-		switch(iK){
+    switch(iK%2){
       case 0:
         for(int i=1;i<3;i++){
-			computeInputForK(i, exponent+1, true);
+          computeInputForK(i, exponent+1, true,iK);
         }
         break;
     }
@@ -563,7 +597,11 @@ void dgRungeKuttaMultirate22::computeInputForK(int iK,int exponent,bool isBuffer
 }
 
 void dgRungeKuttaMultirate22::updateSolution(int exponent,bool isBuffer){
-	//Msg::Info("Updating solution at level %d %s",exponent,isBuffer?"Buffer":"Bulk");
+#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;

Solver/dgRungeKuttaMultirate.h

@@ -17,7 +17,6 @@ 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
-  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);
@@ -40,7 +39,7 @@ class dgRungeKuttaMultirate43: public dgRungeKuttaMultirate{
 };
 
 class dgRungeKuttaMultirate22: public dgRungeKuttaMultirate{
-  void computeInputForK(int iK,int exponent,bool isBuffer);
+  void computeInputForK(int iK,int exponent,bool isBuffer,int upperLeveliK);
   void updateSolution(int exponent,bool isBuffer);
   public:
   double iterate(double dt, dgDofContainer *solution);