DG: conservative multirate is conservative

Solver/TESTCASES/ConvergenceMultirate.lua

@@ -29,6 +29,12 @@ function L2Norm( sol, FCT )
     FCT:set(i,0,s*s) 
   end
 end
+function solF( sol, FCT )
+  for i=0,sol:size1()-1 do
+    s = sol:get(i,0)
+    FCT:set(i,0,s) 
+  end
+end
 --[[ 
      Example of a lua program driving the DG code
 --]]
@@ -63,7 +69,7 @@ FS = functionLua(1, 'initial_condition', {xyzF})
 GC=dgGroupCollection(myModel,2,order)
 solTmp=dgDofContainer(GC,1)
 solTmp:L2Projection(FS)
-dt=GC:splitGroupsForMultirate(100,1,law,solTmp)
+dt=GC:splitGroupsForMultirate(20,1,2,law,solTmp)
 GC:buildGroupsOfInterfaces(myModel,2,order)
 solution=dgDofContainer(GC,1)
 solutionRef=dgDofContainer(GC,1)
@@ -72,7 +78,7 @@ solution:exportGroupIdMsh()
 solution:exportMultirateGroupMsh()
 solutionRef:exportMsh("solution") 
 
-nRefSteps=100
+nRefSteps=40
 dtRef=dt/nRefSteps
 RK=dgRungeKutta()
 print'Compute reference solution'
@@ -86,15 +92,23 @@ end
 --multirateRK=dgRungeKuttaMultirate43(GC,law)
 multirateRK=dgRungeKuttaMultirate22(GC,law)
 integrator=dgFunctionIntegrator(functionLua(1, 'L2Norm', {functionSolution.get()}))
+intSol=dgFunctionIntegrator(functionLua(1, 'solF', {functionSolution.get()}))
 
 oldnorm=1
 for n=0,3 do
   solution:L2Projection(FS)
   print('solnorm ',solution:norm())
+  sN1=fullMatrix(1,1)
+  intSol:compute(solution,sN1)
+  print(string.format("solution integral before: %0.16e",sN1:get(0,0)))
   for i=1,2^n do
     multirateRK:iterate(dt/(2^n),solution)
     print('it ',i)
   end
+  sN2=fullMatrix(1,1)
+  intSol:compute(solution,sN2)
+  print(string.format("solution integral after: %0.16e",sN2:get(0,0)))
+  print(string.format("Mass conservation relative error: %0.16e",(sN1:get(0,0)-sN2:get(0,0))/sN1:get(0,0) ))
   solution:axpy(solutionRef,-1)
   lsq=fullMatrix(1,1)
   integrator:compute(solution,lsq)

Solver/TESTCASES/MultirateAdvection.lua

@@ -60,7 +60,7 @@ FS = functionLua(1, 'initial_condition', {'XYZ'}):getName()
 GC=dgGroupCollection(myModel,2,order)
 solTmp=dgDofContainer(GC,1)
 solTmp:L2Projection(FS)
-dt=GC:splitGroupsForMultirate(2,10,law,solTmp)
+dt=GC:splitGroupsForMultirate(2,1,1,law,solTmp)
 GC:buildGroupsOfInterfaces(myModel,2,order)
 solution=dgDofContainer(GC,1)
 solution:L2Projection(FS)

Solver/dgGroupOfElements.cpp

@@ -18,6 +18,7 @@
 #else
 #include <string.h>
 #endif
+//#define TEST
 
 
 static fullMatrix<double> * dgGetIntegrationRule (MElement *e, int p){
@@ -758,7 +759,7 @@ void dgGroupCollection::buildParallelStructure()
 }
 
 // Split the groups of elements depending on their local time step
-double dgGroupCollection::splitGroupsForMultirate(int maxLevels,int bufferSize,dgConservationLaw *claw, dgDofContainer *solution) {
+double dgGroupCollection::splitGroupsForMultirate(int maxLevels,int innerBufferSize,int outerBufferSize,dgConservationLaw *claw, dgDofContainer *solution) {
   // What are the levels/layers:
   // bulk: elements that are time stepped using the "normal" 4 stage Runge-Kutta
   // innerBuffer: elements that use the small time step but talks to elements using the big time step
@@ -856,6 +857,7 @@ double dgGroupCollection::splitGroupsForMultirate(int maxLevels,int bufferSize,d
 
   // We pass two times by each exponent level: one for the normal(bulk) + innerBuffer groups, and one for the outerBuffer groups
   for(int currentExponent=dtMaxExponent;currentExponent>=0;(!isOuterBufferLayer)?currentExponent--:currentExponent=currentExponent){
+    int currentBufferSize=isOuterBufferLayer?outerBufferSize:innerBufferSize;
     double currentDt=dtRef/pow(2.0,(double)currentExponent);
     // element for the new groups at this level:
     // mapNewGroups[oldGroupId][oldElementId]
@@ -868,7 +870,7 @@ double dgGroupCollection::splitGroupsForMultirate(int maxLevels,int bufferSize,d
       // For buffer AND non buffer layers
       lowerLevelGroupIdStart=lowerLevelGroupIdEnd;
       // We add bufferSize elements (most of the time, bufferSize=1 is enough)
-      for(int iLoop=0;iLoop<bufferSize;iLoop++){
+      for(int iLoop=0;iLoop<currentBufferSize;iLoop++){
         for(int iInterface=0;iInterface<miniInterfaceV->size();iInterface++){
           dgMiniInterface &interface=miniInterfaceV->at(iInterface);
           bool toAdd=false;
@@ -917,7 +919,7 @@ double dgGroupCollection::splitGroupsForMultirate(int maxLevels,int bufferSize,d
 
 #define MAXBUFFERSIZE 100000
       if(currentExponent>0){
-        for(int iLoop=0;iLoop<bufferSize;iLoop++){
+        for(int iLoop=0;iLoop<currentBufferSize;iLoop++){
           for (std::map<int, std::vector<int> >::iterator it = mapNewGroups.begin(); it !=mapNewGroups.end() ; ++it){
             for(int i=0;i<it->second.size();i++){
               bool inInnerBuffer=false;
@@ -969,9 +971,22 @@ double dgGroupCollection::splitGroupsForMultirate(int maxLevels,int bufferSize,d
             oldGroupIds[currentNewGroupId]=it->first;
             newGroup=new dgGroupOfElements(forBulkV,oldGroup->getOrder(),oldGroup->getGhostPartition());
             newGroup->copyPrivateDataFrom(oldGroup);
+            #ifdef TEST
+            if(currentExponent==0){
+              newGroup->_multirateExponent=currentExponent+1;
+              newGroup->_multirateOuterBuffer=true;
+              newGroup->_multirateInnerBuffer=false;
+            }
+            else{
+              newGroup->_multirateExponent=currentExponent;
+              newGroup->_multirateOuterBuffer=false;
+              newGroup->_multirateInnerBuffer=true;
+            }
+            #else
             newGroup->_multirateExponent=currentExponent;
             newGroup->_multirateOuterBuffer=false;
             newGroup->_multirateInnerBuffer=false;
+            #endif
             newGroups.resize(currentNewGroupId+1);
             newGroups[currentNewGroupId]=newGroup;
             currentNewGroupId++;
@@ -1237,7 +1252,7 @@ void dgGroupCollection::registerBindings(binding *b)
   cm->setDescription("Build the group of interfaces, i.e. boundary interfaces and inter-element interfaces");
   cm = cb->addMethod("splitGroupsForMultirate",&dgGroupCollection::splitGroupsForMultirate);
   cm->setDescription("Split the groups according to their own stable time step");
-  cm->setArgNames("maxLevels","bufferSize","claw","solution",NULL);
+  cm->setArgNames("maxLevels","innerBufferSize","outerBufferSize","claw","solution",NULL);
   cm = cb->addMethod("splitGroupsByVerticalLayer",&dgGroupCollection::splitGroupsByVerticalLayer);
   cm->setDescription("Split the groups according vertical layer structure. The first is defined by the topLevelTags.");
   cm->setArgNames("topLevelTags",NULL);

Solver/dgGroupOfElements.h

@@ -250,7 +250,7 @@ class dgGroupCollection {
   void buildGroupsOfElements (GModel *model,int dimension, int order, std::vector<std::string>* physicalTags);
   void buildGroupsOfInterfaces ();
 
-  double splitGroupsForMultirate(int maxLevels,int bufferSize,dgConservationLaw *claw, dgDofContainer *solution);
+  double splitGroupsForMultirate(int maxLevels,int innerBufferSize,int outerBufferSize,dgConservationLaw *claw, dgDofContainer *solution);
   void splitGroupsByVerticalLayer(std::vector<std::string> topLevelTags);
 
   void find (MElement *elementToFind, int &iGroup, int &ithElementOfGroup);