diff --git a/Solver/TESTCASES/ConvergenceMultirate.lua b/Solver/TESTCASES/ConvergenceMultirate.lua
index 8bfa3867a091bd27caa0d922f274d57363cc0b64..bae2734f50afa74534a04c8fcabaa807d94b2d2a 100644
--- a/Solver/TESTCASES/ConvergenceMultirate.lua
+++ b/Solver/TESTCASES/ConvergenceMultirate.lua
@@ -17,8 +17,8 @@ function velocity( XYZ, FCT )
for i=0,XYZ:size1()-1 do
X = XYZ:get(i,0)
Y = XYZ:get(i,1)
- FCT:set(i,0,1)
- FCT:set(i,1,0)
+ FCT:set(i,0,0)
+ FCT:set(i,1,1)
FCT:set(i,2,0)
end
end
@@ -62,12 +62,13 @@ FS = functionLua(1, 'initial_condition', {'XYZ'}):getName()
GC=dgGroupCollection(myModel,2,order)
solTmp=dgDofContainer(GC,1)
solTmp:L2Projection(FS)
-dt=GC:splitGroupsForMultirate(3,law,solTmp)
+dt=GC:splitGroupsForMultirate(20,1,law,solTmp)
GC:buildGroupsOfInterfaces(myModel,2,order)
solution=dgDofContainer(GC,1)
solutionRef=dgDofContainer(GC,1)
solutionRef:L2Projection(FS)
solution:exportGroupIdMsh()
+solution:exportMultirateGroupMsh()
solutionRef:exportMsh("solution")
nRefSteps=100
diff --git a/Solver/TESTCASES/MultirateAdvection.lua b/Solver/TESTCASES/MultirateAdvection.lua
index 4bd85943d51fef3b74f62f40ed6ae0e37b399b79..d1b9af574c86424a57c277083c1eb65df99a9931 100644
--- a/Solver/TESTCASES/MultirateAdvection.lua
+++ b/Solver/TESTCASES/MultirateAdvection.lua
@@ -17,6 +17,7 @@ function velocity( XYZ, FCT )
Y = XYZ:get(i,1)
FCT:set(i,0,0)
FCT:set(i,1,1)
+ FCT:set(i,2,0)
end
end
@@ -48,7 +49,7 @@ elseif (order == 5) then
end
print'*** Create a dg solver ***'
-velF = functionLua(2, 'velocity', {'XYZ'}):getName()
+velF = functionLua(3, 'velocity', {'XYZ'}):getName()
law=dgConservationLawAdvectionDiffusion(velF,"")
zero=functionConstant({0.}):getName();
@@ -59,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(4,law,solTmp)
+dt=GC:splitGroupsForMultirate(2,10,law,solTmp)
GC:buildGroupsOfInterfaces(myModel,2,order)
solution=dgDofContainer(GC,1)
solution:L2Projection(FS)
@@ -95,7 +96,7 @@ time=0
--for i=1,1000
i=0
integrator=dgFunctionIntegrator(functionLua(1, 'toIntegrate', {'XYZ'}):getName())
-while time<0.2 do
+while i<1 do
-- norm1 = RK:iterate43SchlegelJCAM2009(law,dt,solution)
-- TEST with Explicit Euler multirate !!!
norm2 = multirateRK:iterate(dt,solution2)
diff --git a/Solver/dgDofContainer.cpp b/Solver/dgDofContainer.cpp
index c6bce12be812f7d517c0ad84007728a272f86cd0..cbf6afa0af39791c8142f6544b84ae34ac80c3d6 100644
--- a/Solver/dgDofContainer.cpp
+++ b/Solver/dgDofContainer.cpp
@@ -192,6 +192,19 @@ double dgDofContainer::norm() {
return localNorm;
#endif
}
+double dgDofContainer::norm(std::vector<dgGroupOfElements *>gV){
+ double norm=0;
+ for(int i=0;i<gV.size();i++){
+ dgGroupOfElements* g=gV[i];
+ fullMatrix<double>&groupData=getGroupProxy(g);
+ for(int r=0;r<groupData.size1();r++){
+ for(int c=0;c<groupData.size1();c++){
+ norm+=groupData(r,c)*groupData(r,c);
+ }
+ }
+ }
+ return norm;
+}
void dgDofContainer::save(const std::string name) {
FILE *f = fopen (name.c_str(),"wb");
_data->binarySave(f);
@@ -233,7 +246,7 @@ void dgDofContainer::exportGroupIdMsh(){
// the elementnodedata::export does not work !!
std::ostringstream name_oss;
- name_oss<<"groups.msh";
+ name_oss<<"groupsId.msh";
if(Msg::GetCommSize()>1)
name_oss<<"_"<<Msg::GetCommRank();
FILE *f = fopen (name_oss.str().c_str(),"w");
@@ -294,6 +307,76 @@ void dgDofContainer::exportGroupIdMsh(){
*/
}
+void dgDofContainer::exportMultirateGroupMsh(){
+ // the elementnodedata::export does not work !!
+
+ std::ostringstream name_oss;
+ name_oss<<"groupsMultirateType.msh";
+ if(Msg::GetCommSize()>1)
+ name_oss<<"_"<<Msg::GetCommRank();
+ FILE *f = fopen (name_oss.str().c_str(),"w");
+ int COUNT = 0;
+ for (int i=0;i < _groups.getNbElementGroups() ;i++){
+ COUNT += _groups.getElementGroup(i)->getNbElements();
+ }
+ fprintf(f,"$MeshFormat\n2.1 0 8\n$EndMeshFormat\n");
+ fprintf(f,"$ElementNodeData\n");
+ fprintf(f,"1\n");
+ fprintf(f,"\"%s\"\n","GroupsIds");
+ fprintf(f,"1\n");
+ fprintf(f,"0.0\n");
+ fprintf(f,"%d\n", Msg::GetCommSize() > 1 ? 4 : 3);
+ fprintf(f,"0\n 1\n %d\n",COUNT);
+ if(Msg::GetCommSize() > 1) fprintf(f,"%d\n", Msg::GetCommRank());
+ for (int i=0;i < _groups.getNbElementGroups() ;i++){
+ dgGroupOfElements *group = _groups.getElementGroup(i);
+ int groupType=0;//-3*group->getMultirateExponent();
+ if(group->getIsInnerMultirateBuffer())
+ groupType+=1;
+ if(group->getIsOuterMultirateBuffer())
+ groupType+=2;
+ for (int iElement=0 ; iElement< group->getNbElements() ;++iElement) {
+ MElement *e = group->getElement(iElement);
+ int num = e->getNum();
+ fullMatrix<double> sol(getGroupProxy(i), iElement*_nbFields,_nbFields);
+ fprintf(f,"%d %d",num,sol.size1());
+ for (int k=0;k<sol.size1();++k) {
+ fprintf(f," %.16E ",groupType*1.0);
+ }
+ fprintf(f,"\n");
+ }
+ }
+ fprintf(f,"$EndElementNodeData\n");
+ fclose(f);
+ return;
+ // we should discuss that : we do a copy of the solution, this should
+ // be avoided !
+
+ /*std::map<int, std::vector<double> > data;
+
+ for (int i=0;i < _groups.getNbElementGroups() ;i++){
+ dgGroupOfElements *group = _groups.getElementGroup(i);
+ for (int iElement=0 ; iElement< group->getNbElements() ;++iElement) {
+ MElement *e = group->getElement(iElement);
+ int num = e->getNum();
+ fullMatrix<double> sol(getGroupProxy(i), iElement*_nbFields,_nbFields);
+ std::vector<double> val;
+ // val.resize(sol.size2()*sol.size1());
+ val.resize(sol.size1());
+ int counter = 0;
+ // for (int iC=0;iC<sol.size2();iC++)
+ printf("%g %g %g\n",sol(0,0),sol(1,0),sol(2,0));
+ for (int iR=0;iR<sol.size1();iR++)val[counter++] = sol(iR,0);
+ data[num] = val;
+ }
+ }
+
+ PView *pv = new PView (name, "ElementNodeData", _gm, data, 0.0, 1);
+ pv->getData()->writeMSH(name+".msh", false);
+ delete pv;
+ */
+}
+
void dgDofContainer::exportMsh(const std::string name)
{
// the elementnodedata::export does not work !!
@@ -386,7 +469,9 @@ void dgDofContainer::registerBindings(binding *b){
cm->setArgNames("filename", NULL);
cm = cb->addMethod("exportGroupIdMsh",&dgDofContainer::exportGroupIdMsh);
cm->setDescription("Export the group ids for gmsh visualization");
- cm = cb->addMethod("norm",&dgDofContainer::norm);
+ cm = cb->addMethod("exportMultirateGroupMsh",&dgDofContainer::exportMultirateGroupMsh);
+ cm->setDescription("Export the group Multirate properties for gmsh visualization");
+ cm = cb->addMethod("norm",(double (dgDofContainer::*)())&dgDofContainer::norm);
cm->setDescription("Returns the norm of the vector");
cm = cb->addMethod("scale",(void (dgDofContainer::*)(double))&dgDofContainer::scale);
cm->setArgNames("factor",NULL);
diff --git a/Solver/dgDofContainer.h b/Solver/dgDofContainer.h
index c1ce0989748ade9f861ac5d9b6d3161f60972f71..ad926cd1ef44d70dad7215c00b556a5552b066f7 100644
--- a/Solver/dgDofContainer.h
+++ b/Solver/dgDofContainer.h
@@ -28,6 +28,7 @@ public:
void scale(double f);
void scale(std::vector<dgGroupOfElements*>groups, double f);
double norm();
+ double norm(std::vector<dgGroupOfElements*>groups);
void axpy(dgDofContainer &x, double a=1.);
void axpy(std::vector<dgGroupOfElements*>groups,dgDofContainer &x, double a=1.);
inline fullMatrix<double> &getGroupProxy(int gId){ return *(_dataProxys[gId]); }
@@ -44,6 +45,7 @@ public:
void Mesh2Mesh_L2Projection(dgDofContainer &other);
void exportMsh(const std::string name);
void exportGroupIdMsh();
+ void exportMultirateGroupMsh();
static void registerBindings(binding *b);
inline dgGroupCollection *getGroups() { return &_groups; }
diff --git a/Solver/dgGroupOfElements.cpp b/Solver/dgGroupOfElements.cpp
index 158aca05e69c9e228f30280edd39cb250fe2c8d6..dfdb3fe8f066bca608b79155bb2e438acbe45cfc 100644
--- a/Solver/dgGroupOfElements.cpp
+++ b/Solver/dgGroupOfElements.cpp
@@ -19,6 +19,7 @@
#include <string.h>
#endif
+
static fullMatrix<double> * dgGetIntegrationRule (MElement *e, int p){
int npts;
IntPt *pts;
@@ -737,7 +738,7 @@ void dgGroupCollection::buildParallelStructure()
}
// Split the groups of elements depending on their local time step
-double dgGroupCollection::splitGroupsForMultirate(int maxLevels,dgConservationLaw *claw, dgDofContainer *solution) {
+double dgGroupCollection::splitGroupsForMultirate(int maxLevels,int bufferSize,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
@@ -845,32 +846,33 @@ double dgGroupCollection::splitGroupsForMultirate(int maxLevels,dgConservationLa
else{
// Add the neighbors elements to the new groups
// For buffer AND non buffer layers
- int _lowerLevelGroupIdStart=lowerLevelGroupIdStart;
- int _lowerLevelGroupIdEnd=lowerLevelGroupIdEnd;
lowerLevelGroupIdStart=lowerLevelGroupIdEnd;
- for(int iInterface=0;iInterface<miniInterfaceV->size();iInterface++){
- dgMiniInterface &interface=miniInterfaceV->at(iInterface);
- bool toAdd=false;
- // if one of the elements adjacent to the interface is mapped to the previous level
- // we check all elements adjacent to this interface, and add those that does not
- // already have a new group to the current new group
- for(int iConn=0;iConn<interface.connections.size();iConn++){
- int gId=interface.connections[iConn].iGroup;
- int eId=interface.connections[iConn].iElement;
- int newGroupId=newGroupIds[gId][eId];
- if(newGroupId>=0 /*newGroupId >= _lowerLevelGroupIdStart && newGroupId<_lowerLevelGroupIdEnd*/){
- toAdd=true;
- continue;
- }
- }
- if(toAdd){
+ // We add bufferSize elements (most of the time, bufferSize=1 is enough)
+ for(int iLoop=0;iLoop<bufferSize;iLoop++){
+ for(int iInterface=0;iInterface<miniInterfaceV->size();iInterface++){
+ dgMiniInterface &interface=miniInterfaceV->at(iInterface);
+ bool toAdd=false;
+ // if one of the elements adjacent to the interface is mapped to the previous level
+ // we check all elements adjacent to this interface, and add those that does not
+ // already have a new group to the current new group
for(int iConn=0;iConn<interface.connections.size();iConn++){
int gId=interface.connections[iConn].iGroup;
int eId=interface.connections[iConn].iElement;
int newGroupId=newGroupIds[gId][eId];
- if(newGroupId==-1){
- mapNewGroups[gId].push_back(eId);
- newGroupIds[gId][eId]=-2;
+ if(newGroupId>=0 || ( newGroupId>-2-iLoop && newGroupId<-1) ){
+ toAdd=true;
+ continue;
+ }
+ }
+ if(toAdd){
+ for(int iConn=0;iConn<interface.connections.size();iConn++){
+ int gId=interface.connections[iConn].iGroup;
+ int eId=interface.connections[iConn].iElement;
+ int newGroupId=newGroupIds[gId][eId];
+ if(newGroupId==-1){
+ mapNewGroups[gId].push_back(eId);
+ newGroupIds[gId][eId]=-2-iLoop;
+ }
}
}
}
@@ -892,22 +894,41 @@ double dgGroupCollection::splitGroupsForMultirate(int maxLevels,dgConservationLa
// We split those elements into bulk and innerBuffer (i.e. those who have a neighbor with a bigger time step)
lowerLevelGroupIdStart=currentNewGroupId;
+
+#define MAXBUFFERSIZE 100000
+ if(currentExponent>0){
+ for(int iLoop=0;iLoop<bufferSize;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;
+ int oldGId=it->first;
+ int oldEId=it->second[i];
+ // We only consider elements in the map and not in the inner buffer for a different iLoop
+ if(newGroupIds[oldGId][oldEId]>-2 || newGroupIds[oldGId][oldEId]<-(MAXBUFFERSIZE-1) )
+ continue;
+ for(int iNeighbor=0;iNeighbor<elementToNeighbors[oldGId][oldEId].size();iNeighbor++){
+ std::pair<int,int> neighIds=elementToNeighbors[oldGId][oldEId][iNeighbor];
+ if(newGroupIds[neighIds.first][neighIds.second]==-1-iLoop*MAXBUFFERSIZE){
+ inInnerBuffer=true;
+ continue;
+ }
+ }
+ if(inInnerBuffer){
+ newGroupIds[oldGId][oldEId]=-1-(iLoop+1)*MAXBUFFERSIZE;
+ }
+ }
+ }
+ }
+ }
+
for (std::map<int, std::vector<int> >::iterator it = mapNewGroups.begin(); it !=mapNewGroups.end() ; ++it){
if(!it->second.empty()){
std::vector<int>forBulk;
std::vector<int>forInnerBuffer;
for(int i=0;i<it->second.size();i++){
- bool inInnerBuffer=false;
int oldGId=it->first;
int oldEId=it->second[i];
- for(int iNeighbor=0;iNeighbor<elementToNeighbors[oldGId][oldEId].size();iNeighbor++){
- std::pair<int,int> neighIds=elementToNeighbors[oldGId][oldEId][iNeighbor];
- if(newGroupIds[neighIds.first][neighIds.second]==-1){
- inInnerBuffer=true;
- continue;
- }
- }
- if(inInnerBuffer){
+ if( newGroupIds[oldGId][oldEId]>-2 || newGroupIds[oldGId][oldEId]<-(MAXBUFFERSIZE-1) ){
forInnerBuffer.push_back(oldEId);
}
else{
@@ -994,7 +1015,7 @@ double dgGroupCollection::splitGroupsForMultirate(int maxLevels,dgConservationLa
// Some stats
int count=0;
for(int i=0;i<newGroups.size();i++){
- Msg::Info("old: %d, level %d, New group # %d has %d elements",oldGroupIds[i],newGroups[i]->getMultirateExponent(),i,newGroups[i]->getNbElements());
+ Msg::Info("Old group #%d, exponent %d, New group #%d has %d elements",oldGroupIds[i],newGroups[i]->getMultirateExponent(),i,newGroups[i]->getNbElements());
if(newGroups[i]->getIsInnerMultirateBuffer())
Msg::Info("InnerBuffer");
else if(newGroups[i]->getIsOuterMultirateBuffer())
@@ -1182,7 +1203,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","claw","solution",NULL);
+ cm->setArgNames("maxLevels","bufferSize","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 39bc4c2b9b0f8d4079fc3ae27d6d9a45845b4f6c..0190ef5ccdf3d8921a3b95a55754f27186c6d29f 100644
--- a/Solver/dgGroupOfElements.h
+++ b/Solver/dgGroupOfElements.h
@@ -247,7 +247,7 @@ class dgGroupCollection {
void buildGroupsOfElements (GModel *model,int dimension, int order);
void buildGroupsOfInterfaces ();
- double splitGroupsForMultirate(int maxLevels,dgConservationLaw *claw, dgDofContainer *solution);
+ double splitGroupsForMultirate(int maxLevels,int bufferSize,dgConservationLaw *claw, dgDofContainer *solution);
void splitGroupsByVerticalLayer(std::vector<std::string> topLevelTags);
void find (MElement *elementToFind, int &iGroup, int &ithElementOfGroup);
diff --git a/Solver/dgRungeKuttaMultirate.cpp b/Solver/dgRungeKuttaMultirate.cpp
index 0770696a2c9b79bb8ea3877d962a41c5b7149bba..ad35b10a625d8df976041896b00771124c65a8fa 100644
--- a/Solver/dgRungeKuttaMultirate.cpp
+++ b/Solver/dgRungeKuttaMultirate.cpp
@@ -6,12 +6,56 @@
#include "dgGroupOfElements.h"
#include <algorithm>
#include <limits.h>
+#include <stdio.h>
+//#define MULTIRATEVERBOSE
+
+
+static double A43[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}
+ };
+static double AInner43[10][10]={
+ {0, 0, 0, 0, 0, 0, 0, 0, 0, 0},
+ {1.0/4.0 ,0, 0, 0, 0, 0, 0, 0, 0, 0},
+ {-1.0/12.0, 1.0/3.0, 0, 0, 0, 0, 0, 0, 0, 0},
+ {1.0/6.0, -1.0/6.0, 1.0/2.0, 0, 0, 0, 0, 0, 0, 0},
+ {1.0/12.0, 1.0/6.0, 1.0/6.0, 1.0/12.0, 0, 0, 0, 0, 0, 0},
+ {1.0/12.0, 1.0/6.0, 1.0/6.0, 1.0/12.0, 0, 0, 0, 0, 0, 0},
+ {1.0/12.0, 1.0/6.0, 1.0/6.0, 1.0/12.0, 0, 1.0/4.0, 0, 0, 0, 0},
+ {1.0/12.0, 1.0/6.0, 1.0/6.0, 1.0/12.0, 0, -1.0/12.0, 1.0/3.0, 0, 0, 0},
+ {1.0/12.0, 1.0/6.0, 1.0/6.0, 1.0/12.0, 0, 1.0/6.0, -1.0/6.0, 1.0/2.0, 0, 0},
+ {1.0/12.0, 1.0/6.0, 1.0/6.0, 1.0/12.0, 0, 1.0/12.0, 1.0/6.0, 1.0/6.0, 1.0/12.0, 0}
+ };
+
+ // Big step RK43
+static double AOuter43[10][10]={
+ {0, 0, 0, 0, 0, 0, 0, 0, 0, 0},
+ {1.0/4.0 , 0, 0, 0, 0, 0, 0, 0, 0, 0},
+ {1.0/4.0 , 0, 0, 0, 0, 0, 0, 0, 0, 0},
+ {1.0/2.0 , 0, 0, 0, 0, 0, 0, 0, 0, 0},
+ {1.0/2.0 , 0, 0, 0, 0, 0, 0, 0, 0, 0},
+ {-1.0/6.0, 0, 0, 0, 2.0/3.0, 0, 0, 0, 0, 0},
+ {1.0/12.0, 0, 0, 0, 1.0/6.0, 1.0/2.0, 0, 0, 0, 0},
+ {1.0/12.0, 0, 0, 0, 1.0/6.0, 1.0/2.0, 0, 0, 0, 0},
+ {1.0/3.0 , 0, 0, 0, -1.0/3.0, 1, 0, 0, 0, 0},
+ {1.0/3.0 , 0, 0, 0, -1.0/3.0, 1, 0, 0, 0, 0},
+ };
+
+static double b[4]={1.0/6.0, 1.0/3.0, 1.0/3.0, 1.0/6.0};
+static double c[4]={0, 1.0/2.0, 1.0/2.0, 1};
+static double bInner[10]={1.0/12.0, 1.0/6.0, 1.0/6.0,1.0/12.0,0,1.0/12.0, 1.0/6.0, 1.0/6.0,1.0/12.0,0};
+static double cInner[10]={0, 1.0/4.0, 1.0/4.0, 1.0/2.0, 1.0/2.0, 1.0/2.0, 3.0/4.0, 3.0/4.0, 1, 1 };
+static double bOuter[10]={1.0/6.0, 0 , 0 , 0 , 1.0/3.0,1.0/3.0, 0 , 0 , 0 , 1.0/6.0};
+static double cOuter[10]={0, 1.0/4.0, 1.0/4.0, 1.0/2.0, 1.0/2.0, 1.0/2.0, 3.0/4.0, 3.0/4.0, 1, 1 };
dgRungeKuttaMultirate::dgRungeKuttaMultirate(dgGroupCollection* gc,dgConservationLaw *law, int nStages){
_law=law;
_residualVolume=new dgResidualVolume(*_law);
_residualInterface=new dgResidualInterface(*_law);
_residualBoundary=new dgResidualBoundary(*_law);
+ _gc=gc;
_K=new dgDofContainer*[nStages];
for(int i=0;i<nStages;i++){
_K[i]=new dgDofContainer(gc,_law->getNbFields());
@@ -40,12 +84,6 @@ dgRungeKuttaMultirate::dgRungeKuttaMultirate(dgGroupCollection* gc,dgConservatio
_bulkGroupsOfElements[ge->getMultirateExponent()].first.push_back(ge);
}
}
-// std::vector<std::set<dgGroupOfFaces*> >bulkFacesSet;
-// std::vector<std::set<dgGroupOfFaces*> >innerBufferFacesSet;
-// std::vector<std::set<dgGroupOfFaces*> >outerBufferFacesSet;
-// bulkFacesSet.resize(_maxExponent+1);
-// innerBufferFacesSet.resize(_maxExponent+1);
-// outerBufferFacesSet.resize(_maxExponent+1);
for(int iGroup=0;iGroup<gc->getNbFaceGroups();iGroup++){
dgGroupOfFaces *gf=gc->getFaceGroup(iGroup);
const dgGroupOfElements *ge[2];
@@ -66,7 +104,6 @@ dgRungeKuttaMultirate::dgRungeKuttaMultirate(dgGroupCollection* gc,dgConservatio
dgGroupOfFaces *gf=gc->getBoundaryGroup(iGroup);
const dgGroupOfElements *ge[1];
ge[0]=&gf->getGroupOfConnections(0).getGroupOfElements();
- //ge[1]=&gf->getGroupRight();
for(int i=0;i<1;i++){
if(ge[i]->getIsInnerMultirateBuffer()){
_innerBufferGroupsOfElements[ge[i]->getMultirateExponent()].second.push_back(gf);
@@ -103,7 +140,11 @@ dgRungeKuttaMultirate::~dgRungeKuttaMultirate(){
}
void dgRungeKuttaMultirate::computeK(int iK,int exponent,bool isBuffer){
- //Msg::Info("Compute K%d at level %d %s",iK,exponent,isBuffer?"Buffer":"Bulk");
+#ifdef MULTIRATEVERBOSE
+ for(int i=0;i<exponent;i++)
+ printf("\t");
+ printf("Exponent %d, %s, compute K%d\n",exponent,isBuffer?" Buffer":"Not buffer",iK);
+#endif
if(isBuffer){
std::vector<dgGroupOfElements *>&vei=_innerBufferGroupsOfElements[exponent].first;
std::vector<dgGroupOfElements *>&veo=_outerBufferGroupsOfElements[exponent].first;
@@ -124,7 +165,6 @@ void dgRungeKuttaMultirate::computeK(int iK,int exponent,bool isBuffer){
dgGroupOfFaces *faces=*it;
if(faces->getNbGroupOfConnections()==1){
_residualBoundary->computeAndMap1Group(*faces,*_currentInput,*_residual);
- //Msg::Info("Buffer face group %p is boundary in multirate",faces);
}
else{
const dgGroupOfElements *gL = &faces->getGroupOfConnections(0).getGroupOfElements();
@@ -133,11 +173,12 @@ void dgRungeKuttaMultirate::computeK(int iK,int exponent,bool isBuffer){
fullMatrix<double> residuInterface(faces->getNbNodes(),faces->getNbElements()*2*_law->getNbFields());
faces->mapToInterface(_law->getNbFields(), _currentInput->getGroupProxy(gL), _currentInput->getGroupProxy(gR), solInterface);
_residualInterface->compute1Group(*faces,solInterface,_currentInput->getGroupProxy(gL), _currentInput->getGroupProxy(gR),residuInterface);
- //Msg::Info("Buffer face group %p is mapped left or right in multirate",faces);
- if(gL->getMultirateExponent()==exponent && gL->getIsMultirateBuffer())
+ if(gL->getMultirateExponent()==exponent && gL->getIsMultirateBuffer()){
faces->mapLeftFromInterface(_law->getNbFields(), residuInterface,_residual->getGroupProxy(gL));
- if(gR->getMultirateExponent()==exponent && gR->getIsMultirateBuffer())
+ }
+ if(gR->getMultirateExponent()==exponent && gR->getIsMultirateBuffer()){
faces->mapRightFromInterface(_law->getNbFields(), residuInterface,_residual->getGroupProxy(gR));
+ }
}
}
fullMatrix<double> iMassEl, KEl,residuEl;
@@ -236,43 +277,25 @@ void dgRungeKuttaMultirate43::computeInputForK(int iK,int exponent,bool isBuffer
if(exponent>_maxExponent){
return;
}
- //Msg::Info("Input for K%d at level %d %s",iK,exponent,isBuffer?"Buffer":"Bulk");
+#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);
- double _A[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 _AInner[10][10]={
- {0, 0, 0, 0, 0, 0, 0, 0, 0, 0},
- {1.0/4.0 ,0, 0, 0, 0, 0, 0, 0, 0, 0},
- {-1.0/12.0, 1.0/3.0, 0, 0, 0, 0, 0, 0, 0, 0},
- {1.0/6.0, -1.0/6.0, 1.0/2.0, 0, 0, 0, 0, 0, 0, 0},
- {1.0/12.0, 1.0/6.0, 1.0/6.0, 1.0/12.0, 0, 0, 0, 0, 0, 0},
- {1.0/12.0, 1.0/6.0, 1.0/6.0, 1.0/12.0, 0, 0, 0, 0, 0, 0},
- {1.0/12.0, 1.0/6.0, 1.0/6.0, 1.0/12.0, 0, 1.0/4.0, 0, 0, 0, 0},
- {1.0/12.0, 1.0/6.0, 1.0/6.0, 1.0/12.0, 0, -1.0/12.0, 1.0/3.0, 0, 0, 0},
- {1.0/12.0, 1.0/6.0, 1.0/6.0, 1.0/12.0, 0, 1.0/6.0, -1.0/6.0, 1.0/2.0, 0, 0},
- {1.0/12.0, 1.0/6.0, 1.0/6.0, 1.0/12.0, 0, 1.0/12.0, 1.0/6.0, 1.0/6.0, 1.0/12.0, 0}
- };
- // Big step RK43
- double _AOuter[10][10]={
- {0, 0, 0, 0, 0, 0, 0, 0, 0, 0},
- {1.0/4.0 , 0, 0, 0, 0, 0, 0, 0, 0, 0},
- {1.0/4.0 , 0, 0, 0, 0, 0, 0, 0, 0, 0},
- {1.0/2.0 , 0, 0, 0, 0, 0, 0, 0, 0, 0},
- {1.0/2.0 , 0, 0, 0, 0, 0, 0, 0, 0, 0},
- {-1.0/6.0, 0, 0, 0, 2.0/3.0, 0, 0, 0, 0, 0},
- {1.0/12.0, 0, 0, 0, 1.0/6.0, 1.0/2.0, 0, 0, 0, 0},
- {1.0/12.0, 0, 0, 0, 1.0/6.0, 1.0/2.0, 0, 0, 0, 0},
- {1.0/3.0 , 0, 0, 0, -1.0/3.0, 1, 0, 0, 0, 0},
- {1.0/3.0 , 0, 0, 0, -1.0/3.0, 1, 0, 0, 0, 0},
- };
-// Msg::Info("K%d, %d, %s",iK,exponent,isBuffer?"Buffer":"Bulk");
- // compute K[iK] for this exponent and buffer
- if(isBuffer){
+ 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(A43[iK][i]!=0.0){
+ _currentInput->axpy(gV,*_K[i],A43[iK][i]*localDt);
+ }
+ }
+ }
+ else{
std::vector<dgGroupOfElements *>&gVi=_innerBufferGroupsOfElements[exponent].first;
std::vector<dgGroupOfElements *>&gVo=_outerBufferGroupsOfElements[exponent].first;
_currentInput->scale(gVi,0.0);
@@ -280,21 +303,44 @@ void dgRungeKuttaMultirate43::computeInputForK(int iK,int exponent,bool isBuffer
_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);
+ if(AInner43[iK][i]!=0.0)
+ _currentInput->axpy(gVi,*_K[i],AInner43[iK][i]*localDt*2);
+ if(AOuter43[iK][i]!=0.0)
+ _currentInput->axpy(gVo,*_K[i],AOuter43[iK][i]*localDt*2);
}
- }
- else{
- 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);
+ /*
+ NOT NEEDED
+ // 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<4) || (iK>5 && iK<9) ){
+ 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,0,1,2,2,2,2,3};
+ for(int i=0;i<iK;i++){
+ if(AOuter43[iK][i]!=0.0){
+ // ON DIRAIT QUE CETTE LIGNE DE CODE NE FAIT RIEN ...
+ _currentInput->axpy(gVbigger,*_K[idOuter2Bulk[i]],AOuter43[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>=5){
+ for(int i=0;i<4;i++){
+ _currentInput->axpy(gViBigger,*_K[i],b[i]*localDt*2);
+ }
+ s+=5;
+ }
+ for(int i=0;i<iK;i++){
+ if(AOuter43[iK][i]!=0.0){
+ _currentInput->axpy(gViBigger,*_K[idOuter2Bulk[i]+s],AOuter43[iK][i]*localDt*2);
+ }
}
}
+ */
}
if(!isBuffer){
@@ -318,12 +364,10 @@ void dgRungeKuttaMultirate43::computeInputForK(int iK,int exponent,bool isBuffer
computeInputForK(iK%5,exponent,false);
}
}
- computeK(iK,exponent,isBuffer);
-// Msg::Info("Multirate %d %0.16e",iK,_K[iK]->norm());
- if( (iK==3 && !isBuffer) || (iK==9 && isBuffer) ){
- updateSolution(exponent,isBuffer);
- }
- if(!isBuffer){
+ if(exponent==0){
+ computeK(iK,exponent,false);
+ if(iK==3)
+ updateSolution(exponent,false);
switch(iK){
case 0:
for(int i=1;i<4;i++){
@@ -337,17 +381,37 @@ void dgRungeKuttaMultirate43::computeInputForK(int iK,int exponent,bool isBuffer
break;
}
}
+ if(isBuffer && exponent>0){
+ if( (iK%5)<4)
+ computeK(iK%5,exponent,false);
+ computeK(iK,exponent,true);
+ if( (iK%5)==3)
+ updateSolution(exponent,false);
+ if(iK==9)
+ updateSolution(exponent,true);
+ switch(iK%5){
+ case 0:
+ for(int i=1;i<4;i++){
+ computeInputForK(i,exponent+1,true);
+ }
+ break;
+ case 2:
+ for(int i=6;i<9;i++){
+ computeInputForK(i,exponent+1,true);
+ }
+ break;
+ }
+ }
}
void dgRungeKuttaMultirate43::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);
- double b[4]={1.0/6.0, 1.0/3.0, 1.0/3.0, 1.0/6.0};
- double c[4]={0, 1.0/2.0, 1.0/2.0, 1};
- double bInner[10]={1.0/12.0, 1.0/6.0, 1.0/6.0,1.0/12.0,0,1.0/12.0, 1.0/6.0, 1.0/6.0,1.0/12.0,0};
- double cInner[10]={0, 1.0/4.0, 1.0/4.0, 1.0/2.0, 1.0/2.0, 1.0/2.0, 3.0/4.0, 3.0/4.0, 1, 1 };
- double bOuter[10]={1.0/6.0, 0 , 0 , 0 , 1.0/3.0,1.0/3.0, 0 , 0 , 0 , 1.0/6.0};
- double cOuter[10]={0, 1.0/4.0, 1.0/4.0, 1.0/2.0, 1.0/2.0, 1.0/2.0, 3.0/4.0, 3.0/4.0, 1, 1 };
if(isBuffer){
std::vector<dgGroupOfElements *>&gVi=_innerBufferGroupsOfElements[exponent].first;
std::vector<dgGroupOfElements *>&gVo=_outerBufferGroupsOfElements[exponent].first;
@@ -364,6 +428,8 @@ void dgRungeKuttaMultirate43::updateSolution(int exponent,bool isBuffer){
if(b[i]!=0.0)
_solution->axpy(gV,*_K[i],b[i]*localDt);
}
+ _currentInput->scale(gV,0.0);
+ _currentInput->axpy(gV,*_solution,1.0);
}
}
diff --git a/Solver/dgRungeKuttaMultirate.h b/Solver/dgRungeKuttaMultirate.h
index ddeb70c79d7c3ee9bbd5a95d26fafc1cc0ef330f..f7a98fa322403950882548e48954835356b2f785 100644
--- a/Solver/dgRungeKuttaMultirate.h
+++ b/Solver/dgRungeKuttaMultirate.h
@@ -13,6 +13,7 @@ class dgRungeKuttaMultirate: public dgRungeKutta{
dgResidualVolume *_residualVolume;
dgResidualInterface *_residualInterface;
dgResidualBoundary *_residualBoundary;
+ dgGroupCollection *_gc;
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