diff --git a/Solver/TESTCASES/ConvergenceMultirate.lua b/Solver/TESTCASES/ConvergenceMultirate.lua
index bae2734f50afa74534a04c8fcabaa807d94b2d2a..71ac0ad30ddc68453fa16f87d614e2e2553f3fb0 100644
--- a/Solver/TESTCASES/ConvergenceMultirate.lua
+++ b/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
diff --git a/Solver/dgRungeKuttaMultirate.cpp b/Solver/dgRungeKuttaMultirate.cpp
index ee67970e8db6fe585404ebe0d4b2b7a342265323..9f8238e04418894ee554c1a0226208fcd75e3610 100644
--- a/Solver/dgRungeKuttaMultirate.cpp
+++ b/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,18 +586,22 @@ void dgRungeKuttaMultirate22::computeInputForK(int iK,int exponent,bool isBuffer
updateSolution(exponent, false);
if(iK==3)
updateSolution(exponent, true);
- switch(iK){
- case 0:
- for(int i=1;i<3;i++){
- computeInputForK(i, exponent+1, true);
- }
- break;
- }
+ switch(iK%2){
+ case 0:
+ for(int i=1;i<3;i++){
+ computeInputForK(i, exponent+1, true,iK);
+ }
+ break;
+ }
}
}
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;
diff --git a/Solver/dgRungeKuttaMultirate.h b/Solver/dgRungeKuttaMultirate.h
index f7a98fa322403950882548e48954835356b2f785..79004a329a66297b6b02b696a6d6f2a90789264f 100644
--- a/Solver/dgRungeKuttaMultirate.h
+++ b/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);