diff --git a/Solver/TESTCASES/ConvergenceMultirate.lua b/Solver/TESTCASES/ConvergenceMultirate.lua
index ad094513a8ac4ab756293ecc31c171ce588d431d..c2968487d690df7743f5a164f9c68bb94dda70bf 100644
--- a/Solver/TESTCASES/ConvergenceMultirate.lua
+++ b/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)
diff --git a/Solver/TESTCASES/MultirateAdvection.lua b/Solver/TESTCASES/MultirateAdvection.lua
index d1b9af574c86424a57c277083c1eb65df99a9931..4922582f73b3ad2ed8ec7b3349d78bba8adbf1b5 100644
--- a/Solver/TESTCASES/MultirateAdvection.lua
+++ b/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)
diff --git a/Solver/dgGroupOfElements.cpp b/Solver/dgGroupOfElements.cpp
index 618c1c4f8d8d5d57689047414ecd5ec0e09b18ac..f9c34e96527994b17b9a3e9f6a4a6959863a2607 100644
--- a/Solver/dgGroupOfElements.cpp
+++ b/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);
diff --git a/Solver/dgGroupOfElements.h b/Solver/dgGroupOfElements.h
index 9a8a25b0e12a3fe331dfafe8fd7c549a9721e2ce..42c8e1c218cc47d496275b340089284f242f3201 100644
--- a/Solver/dgGroupOfElements.h
+++ b/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);