diff --git a/Common/LuaBindings.cpp b/Common/LuaBindings.cpp
index 0011274ee36142c3784c687412ef70059be0c106..347a7b6784f0206275c623357fe5003b3ba7e532 100644
--- a/Common/LuaBindings.cpp
+++ b/Common/LuaBindings.cpp
@@ -343,6 +343,7 @@ binding::binding(){
dgPerfectGasLaw2dRegisterBindings(this);
dgResidual::registerBindings(this);
dgRungeKutta::registerBindings(this);
+ dgRungeKuttaMultirate::registerBindings(this);
dgSlopeLimiterRegisterBindings(this);
dgSystemOfEquations::registerBindings(this);
fullMatrix<double>::registerBindings(this);
diff --git a/Mesh/HighOrder.cpp b/Mesh/HighOrder.cpp
index 73f512ba6d916413f231c0dda2a529a8890926d9..61a846627520c9d6c6fd3cba1a8d3e1983169a55 100644
--- a/Mesh/HighOrder.cpp
+++ b/Mesh/HighOrder.cpp
@@ -996,7 +996,7 @@ static void checkHighOrderTetrahedron(const char* cc, GModel *m,
minGGlob, avg / (count ? count : 1));
}
-extern double mesh_functional_distorsion(MTriangle *t, double u, double v);
+extern double mesh_functional_distorsion(MElement *t, double u, double v);
static void printJacobians(GModel *m, const char *nm)
{
diff --git a/Solver/TESTCASES/rect.geo b/Solver/TESTCASES/rect.geo
index 9ab82d5364eb20af385cf148029a3365bec7b0da..6f581ea1d64723210792af9058832d22aad94c69 100644
--- a/Solver/TESTCASES/rect.geo
+++ b/Solver/TESTCASES/rect.geo
@@ -18,8 +18,8 @@ Plane Surface(11) = {10};
Physical Surface("sprut") = {11, 9};
Physical Line("Walls") = {5, 7, 3, 4, 1};
Physical Line("Top") = {6};
-Transfinite Line {5, 7, 1, 3} = 50 Using Progression 1;
-Transfinite Line {6, 4, 2} = 25 Using Progression 1;
-Transfinite Surface {9} Alternated;
-Transfinite Surface {11} Alternated;
-Recombine Surface {9, 11};
+//Recombine Surface {9, 11};
+
+Field[1] = MathEval;
+Field[1].F = "0.01*1+0.01*1.7*y";
+Background Field = 1;
diff --git a/Solver/dgAlgorithm.cpp b/Solver/dgAlgorithm.cpp
index 4246e56c1ca33ed786a51a7cf05d205795b088e0..0366c576e59555cacb87affd418b3590f4f9f180 100644
--- a/Solver/dgAlgorithm.cpp
+++ b/Solver/dgAlgorithm.cpp
@@ -19,6 +19,64 @@
// works for any number of groups
+/*
+void dgAlgorithm::residualForSomeGroups( const dgConservationLaw &claw,
+ std::vector<dgGroupOfElements *>groups,
+ dgGroupCollection &groupCollection,
+ dgDofContainer &solution,
+ dgDofContainer &residu)
+{
+ solution.scatter();
+ int nbFields=claw.getNbFields();
+ //volume term
+ std::set<dgGroupOfFaces *>groupOfInternalFacesToUpdate;
+ std::vector<dgGroupOfFaces *>groupOfBoundaryFacesToUpdate;
+ for(size_t i=0;i<groups.size() ; i++) {
+ int groupId=groups[i]->getId();
+ residu.getGroupProxy(groupId).scale(0);
+ residualVolume(claw,*groups[i],solution.getGroupProxy(groupId),residu.getGroupProxy(groupId));
+ groupOfInternalFacesToUpdate.insert(groups[i]->getGroupsOfNeighboringFaces()->begin(),groups[i]->getGroupsOfNeighboringFaces()->end());
+ groupOfBoundaryFacesToUpdate.insert(groupOfBoundaryFacesToUpdate.end(),groups[i]->getGroupsOfBoundaryFaces()->begin(),groups[i]->getGroupsOfBoundaryFaces()->end());
+ groupOfInternalFacesToUpdate.insert(groups[i]->getGroupOfInsideFaces());
+ }
+ for(std::set<dgGroupOfFaces *>::iterator it=groupOfInternalFacesToUpdate.begin();it!=groupOfInternalFacesToUpdate.end();it++) {
+ dgGroupOfFaces *faces = *it;
+ int iGroupLeft = -1, iGroupRight = -1;
+ for(size_t j=0;j<groupCollection.getNbElementGroups(); j++) {
+ dgGroupOfElements *groupj = groupCollection.getElementGroup(j);
+ if (groupj == &faces->getGroupLeft()) iGroupLeft = j;
+ if (groupj == &faces->getGroupRight()) iGroupRight= j;
+ }
+ for(size_t j=0;j<groupCollection.getNbGhostGroups(); j++) {
+ dgGroupOfElements *groupj = groupCollection.getGhostGroup(j);
+ if (groupj == &faces->getGroupLeft()) iGroupLeft = j;
+ if (groupj == &faces->getGroupLeft())iGroupLeft = j + groupCollection.getNbElementGroups();
+ if (groupj == &faces->getGroupRight())iGroupRight= j + groupCollection.getNbElementGroups();
+ }
+ fullMatrix<double> solInterface(faces->getNbNodes(),faces->getNbElements()*2*nbFields);
+ fullMatrix<double> residuInterface(faces->getNbNodes(),faces->getNbElements()*2*nbFields);
+ faces->mapToInterface(nbFields, solution.getGroupProxy(iGroupLeft), solution.getGroupProxy(iGroupRight), solInterface);
+ residualInterface(claw,*faces,solInterface,solution.getGroupProxy(iGroupLeft), solution.getGroupProxy(iGroupRight),residuInterface);
+ faces->mapFromInterface(nbFields, residuInterface,residu.getGroupProxy(iGroupLeft), residu.getGroupProxy(iGroupRight));
+ }
+ //boundaries
+ for(size_t i=0;i<groupOfBoundaryFacesToUpdate.size() ; i++) {
+ dgGroupOfFaces *faces = groupOfBoundaryFacesToUpdate[i];
+ int iGroupLeft = -1, iGroupRight = -1;
+ for(size_t j=0;j<groupCollection.getNbElementGroups(); j++) {
+ dgGroupOfElements *groupj = groupCollection.getElementGroup(j);
+ if (groupj == &faces->getGroupLeft())iGroupLeft = j;
+ if (groupj == &faces->getGroupRight())iGroupRight= j;
+ }
+ fullMatrix<double> solInterface(faces->getNbNodes(),faces->getNbElements()*nbFields);
+ fullMatrix<double> residuInterface(faces->getNbNodes(),faces->getNbElements()*nbFields);
+ faces->mapToInterface(nbFields, solution.getGroupProxy(iGroupLeft), solution.getGroupProxy(iGroupRight), solInterface);
+ residualBoundary(claw,*faces,solInterface,solution.getGroupProxy(iGroupLeft),residuInterface);
+ faces->mapFromInterface(nbFields, residuInterface, residu.getGroupProxy(iGroupLeft), residu.getGroupProxy(iGroupRight));
+ }
+}
+*/
+
void dgAlgorithm::computeElementaryTimeSteps ( //dofManager &dof, // the DOF manager (maybe useless here)
const dgConservationLaw &claw, // the conservation law
const dgGroupOfElements & group,
diff --git a/Solver/dgDofContainer.cpp b/Solver/dgDofContainer.cpp
index 7a2b0d0cf6a1da41a2ca8908a1e77b1e1838a917..2f63ec8fd30a3828ea68b9ba1dfd6bb40ff42a65 100644
--- a/Solver/dgDofContainer.cpp
+++ b/Solver/dgDofContainer.cpp
@@ -166,6 +166,14 @@ void dgDofContainer::axpy(dgDofContainer &x, double a)
_data->axpy(*x._data,a);
_ghostData->axpy(*x._ghostData,a);
}
+void dgDofContainer::axpy(std::vector<dgGroupOfElements*>groupsVector,dgDofContainer &x, double a){
+ for(int i=0;i<groupsVector.size();i++){
+ dgGroupOfElements *g=groupsVector[i];
+ fullMatrix<double> &proxy=getGroupProxy(g);
+ fullMatrix<double> &xProxy=x.getGroupProxy(g);
+ proxy.add(xProxy,a);
+ }
+}
double dgDofContainer::norm() {
double localNorm = _data->norm();
@@ -213,6 +221,70 @@ void dgDofContainer::L2Projection(std::string functionName){
}
}
+void dgDofContainer::exportGroupIdMsh(){
+ // the elementnodedata::export does not work !!
+
+ std::ostringstream name_oss;
+ name_oss<<"groups.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);
+ 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 ",i*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)
{
@@ -236,7 +308,7 @@ void dgDofContainer::exportMsh(const std::string name)
fprintf(f,"1\n");
fprintf(f,"%d\n", _mshStep); // should print actual time here
fprintf(f,"%d\n", Msg::GetCommSize() > 1 ? 4 : 3);
- fprintf(f,"%d\n 1\n %d\n", _mshStep, COUNT);
+ fprintf(f,"%d\n 1\n %d\n", 0/*_mshStep*/, COUNT);
if(Msg::GetCommSize() > 1) fprintf(f,"%d\n", Msg::GetCommRank());
for (int i=0;i < _groups.getNbElementGroups() ;i++){
dgGroupOfElements *group = _groups.getElementGroup(i);
@@ -300,4 +372,6 @@ void dgDofContainer::registerBindings(binding *b){
cm = cb->addMethod("exportMsh",&dgDofContainer::exportMsh);
cm->setDescription("Export the dof for gmsh visualization");
cm->setArgNames("filename", NULL);
+ cm = cb->addMethod("exportGroupIdMsh",&dgDofContainer::exportGroupIdMsh);
+ cm->setDescription("Export the group ids for gmsh visualization");
}
diff --git a/Solver/dgDofContainer.h b/Solver/dgDofContainer.h
index 12094bf1ad11d47429c7f8686bcd1b818a7d6c50..6c3f15ac14f305b70a4ff8b3e8973921ce461a09 100644
--- a/Solver/dgDofContainer.h
+++ b/Solver/dgDofContainer.h
@@ -28,6 +28,7 @@ public:
void scale(double f);
double norm();
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]); }
inline fullMatrix<double> &getGroupProxy(const dgGroupOfElements* g){ return *(_dataProxys[_groupId[g]]); }
dgDofContainer (dgGroupCollection *groups, int nbFields);
@@ -41,6 +42,7 @@ public:
void L2Projection(std::string functionName);
void Mesh2Mesh_L2Projection(dgDofContainer &other);
void exportMsh(const std::string name);
+ void exportGroupIdMsh();
static void registerBindings(binding *b);
inline dgGroupCollection *getGroups() { return &_groups; }
diff --git a/Solver/dgGroupOfElements.cpp b/Solver/dgGroupOfElements.cpp
index c21ccaf325f7af74111c244e99859be4fb4a60ad..48b36559f1970c48d63340d4ebcbbb90fa32718d 100644
--- a/Solver/dgGroupOfElements.cpp
+++ b/Solver/dgGroupOfElements.cpp
@@ -1,5 +1,5 @@
-#include <iostream>
#include "GmshConfig.h"
+#include "GmshMessage.h"
#include "dgGroupOfElements.h"
#include "dgConservationLaw.h"
#include "dgDofContainer.h"
@@ -131,6 +131,10 @@ dgGroupOfElements::dgGroupOfElements(const std::vector<MElement*> &e,
}
}
+void dgGroupOfElements::copyPrivateDataFrom(const dgGroupOfElements *from){
+// Msg::Error("dgGroupOfElements::copyPrivateDataFrom not yet implemented, ask Richard");
+}
+
dgGroupOfElements::~dgGroupOfElements(){
delete _integration;
delete _redistributionFluxes[0];
@@ -622,6 +626,36 @@ void dgGroupOfFaces::mapFromInterface ( int nFields,
}
}
}
+void dgGroupOfFaces::mapLeftFromInterface ( int nFields,
+ const fullMatrix<double> &v,
+ fullMatrix<double> &vLeft
+ )
+{
+ for(int i=0; i<getNbElements(); i++) {
+ const std::vector<int> &closureRight = getClosureRight(i);
+ const std::vector<int> &closureLeft = getClosureLeft(i);
+ for (int iField=0; iField<nFields; iField++){
+ for(size_t j =0; j < closureLeft.size(); j++)
+ vLeft(closureLeft[j], _left[i]*nFields + iField) += v(j, i*2*nFields + iField);
+ }
+ }
+}
+void dgGroupOfFaces::mapRightFromInterface ( int nFields,
+ const fullMatrix<double> &v,
+ fullMatrix<double> &vRight
+ )
+{
+ for(int i=0; i<getNbElements(); i++) {
+ const std::vector<int> &closureRight = getClosureRight(i);
+ const std::vector<int> &closureLeft = getClosureLeft(i);
+ for (int iField=0; iField<nFields; iField++){
+ for(size_t j =0; j < closureRight.size(); j++)
+ vRight(closureRight[j], _right[i]*nFields + iField) += v(j, (i*2+1)*nFields + iField);
+ }
+ }
+}
+
+
/*
void dgAlgorithm::buildGroupsOfFaces(GModel *model, int dim, int order,
@@ -675,8 +709,11 @@ void dgGroupCollection::buildGroupsOfElements(GModel *model, int dim, int order)
}
// do a group of elements for every element type in the mesh
+ int id=_elementGroups.size();
for (std::map<int, std::vector<MElement *> >::iterator it = localElements.begin(); it !=localElements.end() ; ++it){
- _elementGroups.push_back(new dgGroupOfElements(it->second,order));
+ dgGroupOfElements *newGroup=new dgGroupOfElements(it->second,order);
+ _elementGroups.push_back(newGroup);
+ id++;
}
}
@@ -745,7 +782,7 @@ void dgGroupCollection::buildGroupsOfInterfaces() {
for (int i=0;i<_elementGroups.size();i++){
- // create a group of faces on the boundary for every group ef elements
+ // create a group of faces on the boundary for every group of elements
switch(dim) {
case 1 : {
std::map<const std::string, std::set<MVertex*> >::iterator mapIt;
@@ -782,7 +819,11 @@ void dgGroupCollection::buildGroupsOfInterfaces() {
}
}
// create a group of faces for every face that is common to elements of the same group
- _faceGroups.push_back(new dgGroupOfFaces(*_elementGroups[i],order));
+ dgGroupOfFaces *gof=new dgGroupOfFaces(*_elementGroups[i],order);
+ if(gof->getNbElements())
+ _faceGroups.push_back(gof);
+ else
+ delete gof;
// create a groupe of d
for (int j=i+1;j<_elementGroups.size();j++){
dgGroupOfFaces *gof = new dgGroupOfFaces(*_elementGroups[i],*_elementGroups[j],order);
@@ -795,7 +836,11 @@ void dgGroupCollection::buildGroupsOfInterfaces() {
//create ghost groups
for(int i=0;i<Msg::GetCommSize();i++){
for (std::map<int, std::vector<MElement *> >::iterator it = ghostElements[i].begin(); it !=ghostElements[i].end() ; ++it){
- _ghostGroups.push_back(new dgGroupOfElements(it->second,order,i));
+ dgGroupOfElements *gof=new dgGroupOfElements(it->second,order,i);
+ if (gof->getNbElements())
+ _ghostGroups.push_back(gof);
+ else
+ delete gof;
}
}
//create face group for ghostGroups
@@ -875,30 +920,283 @@ void dgGroupCollection::buildGroupsOfInterfaces() {
// Split the groups of elements depending on their local time step
-void dgGroupCollection::splitGroupsForMultirate(double refDt,dgConservationLaw *claw){
- // 1) find the element with the smallest stable time step
- double sr = DBL_MAX;
- dgDofContainer *solution = new dgDofContainer((this),claw->getNbFields());
- solution->scale(0.);
+void dgGroupCollection::splitGroupsForMultirate(double dtRef,dgConservationLaw *claw, dgDofContainer *solution){
+ Msg::Info("Splitting Groups for multirate time stepping");
+ int maxNumElems=getElementGroup(0)->getElement(0)->getGlobalNumber()+1;
+ std::vector<int>oldGroupIds;
+ oldGroupIds.resize(maxNumElems);
+ std::vector<MElement*>allElements;
+ allElements.resize(maxNumElems);
+ for(int iGroup=0;iGroup<getNbElementGroups();iGroup++){
+ dgGroupOfElements *elGroup=getElementGroup(iGroup);
+ for(int iElement=0;iElement<elGroup->getNbElements();iElement++){
+ MElement *el=elGroup->getElement(iElement);
+ oldGroupIds[el->getNum()]=iGroup;
+ allElements[el->getNum()]=el;
+ }
+ }
+
+ std::vector<int>newGroupIds;
+ newGroupIds.assign(maxNumElems,-1);
+
+ std::vector<std::vector<int> >elementToNeighbors;
+ elementToNeighbors.resize(maxNumElems);
+
+ switch(getElementGroup(0)->getElement(0)->getDim()){
+ case 1:
+ {
+ std::map<MVertex*,int> vertexMap;
+ for(int iGroup=0;iGroup<getNbElementGroups();iGroup++){
+ dgGroupOfElements *elGroup=getElementGroup(iGroup);
+ for(int iElement=0; iElement<elGroup->getNbElements(); iElement++){
+ MElement *el = elGroup->getElement(iElement);
+ for (int iVertex=0; iVertex<el->getNumVertices(); iVertex++){
+ MVertex* vertex = el->getVertex(iVertex);
+ if(vertexMap.find(vertex) == vertexMap.end()){
+ vertexMap[vertex] = el->getNum();
+ }else{
+ elementToNeighbors[vertexMap[vertex]].push_back(el->getNum());
+ elementToNeighbors[el->getNum()].push_back(vertexMap[vertex]);
+ }
+ }
+ }
+ }
+ vertexMap.clear();
+ }
+ break;
+ case 2:
+ {
+ std::map<MEdge,int,Less_Edge> edgeMap;
+ for(int iGroup=0;iGroup<getNbElementGroups();iGroup++){
+ dgGroupOfElements *elGroup=getElementGroup(iGroup);
+ for(int iElement=0; iElement<elGroup->getNbElements(); iElement++){
+ MElement *el = elGroup->getElement(iElement);
+ for (int iEdge=0; iEdge<el->getNumEdges(); iEdge++){
+ MEdge edge = el->getEdge(iEdge);
+ if(edgeMap.find(edge) == edgeMap.end()){
+ edgeMap[edge] = el->getNum();
+ }else{
+ elementToNeighbors[edgeMap[edge]].push_back(el->getNum());
+ elementToNeighbors[el->getNum()].push_back(edgeMap[edge]);
+ }
+ }
+ }
+ }
+ edgeMap.clear();
+ }
+ break;
+ case 3:
+ {
+ std::map<MFace,int,Less_Face> faceMap;
+ for(int iGroup=0;iGroup<getNbElementGroups();iGroup++){
+ dgGroupOfElements *elGroup=getElementGroup(iGroup);
+ for(int iElement=0; iElement<elGroup->getNbElements(); iElement++){
+ MElement *el = elGroup->getElement(iElement);
+ for (int iFace=0; iFace<el->getNumFaces(); iFace++){
+ MFace face = el->getFace(iFace);
+ if(faceMap.find(face) == faceMap.end()){
+ faceMap[face] = el->getNum();
+ }else{
+ elementToNeighbors[faceMap[face]].push_back(el->getNum());
+ elementToNeighbors[el->getNum()].push_back(faceMap[face]);
+ }
+ }
+ }
+ }
+ faceMap.clear();
+ }
+ break;
+ default:
+ break;
+ }
+
+ // find the range of time steps
+ double dtMin = DBL_MAX;
+ double dtMax = 0;
+ std::vector<double> *DTS=new std::vector<double>[getNbElementGroups()];
for (int i=0;i<getNbElementGroups();i++){
- std::vector<double> DTS;
- dgAlgorithm::computeElementaryTimeSteps(*claw, *getElementGroup(i), solution->getGroupProxy(i), DTS);
- for (int k=0;k<DTS.size();k++){
- std::cout << "SR: " << DTS[k] << std::endl;
- sr = std::min(sr,DTS[k]);
+ dgAlgorithm::computeElementaryTimeSteps(*claw, *getElementGroup(i), solution->getGroupProxy(i), DTS[i]);
+ for (int k=0;k<DTS[i].size();k++){
+ dtMin = std::min(dtMin,DTS[i][k]);
+ dtMax = std::max(dtMax,DTS[i][k]);
}
}
- delete solution;
- #ifdef HAVE_MPI
- double sr_min;
- MPI_Allreduce((void *)&sr, &sr_min, 1, MPI_DOUBLE, MPI_MIN, MPI_COMM_WORLD);
- sr=sr_min;
- #endif
- // sr is the time step for the most constrained element.
- std::cout << "SR: " << sr << std::endl;
+ std::vector<double>localDt;
+ localDt.resize(maxNumElems);
+ for (int i=0;i<getNbElementGroups();i++){
+ dgGroupOfElements *elGroup=getElementGroup(i);
+ for(int j=0;j<DTS[i].size();j++){
+ MElement *el=elGroup->getElement(j);
+ localDt[el->getNum()]=DTS[i][j];
+ }
+ }
+ delete []DTS;
+#ifdef HAVE_MPI
+ double dtMin_min;
+ MPI_Allreduce((void *)&dtMin, &dtMin_min, 1, MPI_DOUBLE, MPI_MIN, MPI_COMM_WORLD);
+ dtMin=dtMin_min;
+ double dtMax_max;
+ MPI_Allreduce((void *)&dtMax, &dtMax_max, 1, MPI_DOUBLE, MPI_MAX, MPI_COMM_WORLD);
+ dtMax=dtMax_max;
+#endif
+ // dtMin is the time step for the most constrained element.
+ if(dtRef<=dtMin){
+ Msg::Info("No multirate, the reference time step is stable for all elements.");
+ return;
+ }
+
+ Msg::Info("Time step for standard RK should be %e",dtMin);
+ Msg::Info("Multirate base time step should be %e",dtMax);
+
+ dtRef=dtMax*0.8;
+ // time steps are dtRef*2^(-dtExponent), with dtExponent ranging in [0:dtMaxExponent]
+ int dtMaxExponent=0;
+ while(dtRef/pow(2.0,(double)dtMaxExponent)>dtMin)
+ dtMaxExponent++;
+ _dtMaxExponent=dtMaxExponent;
+ std::vector<MElement *>currentNewGroup;
+ std::vector< dgGroupOfElements* >newGroups;// indexed by newGroupId
+ std::map<int,int>oldToNewGroupsMap;// oldGroupId to newGroupId
+
+ int lowerLevelGroupIdStart=-1;
+ int lowerLevelGroupIdEnd=-1;
+
+ bool isOuterBufferLayer=false;
+ int currentNewGroupId=0;
+ int loopId=0;
+ for(int currentExponent=dtMaxExponent;currentExponent>=0;(!isOuterBufferLayer)?currentExponent--:currentExponent=currentExponent){
+ double currentDt=dtRef/pow(2.0,(double)currentExponent);
+ std::map<int,std::vector<MElement*> >mapNewGroups;
+ if(lowerLevelGroupIdStart==-1){
+ lowerLevelGroupIdStart=0;
+ }
+ else{
+ // Add the neighbors elements to the new groups
+ int _lowerLevelGroupIdStart=lowerLevelGroupIdStart;
+ int _lowerLevelGroupIdEnd=lowerLevelGroupIdEnd;
+ lowerLevelGroupIdStart=lowerLevelGroupIdEnd;
+ for(int iNewGroup=_lowerLevelGroupIdStart;iNewGroup<_lowerLevelGroupIdEnd;iNewGroup++){
+ for(int iElement=0; iElement<newGroups[iNewGroup]->getNbElements(); iElement++){
+ MElement *el = newGroups[iNewGroup]->getElement(iElement);
+ for(int iNeighbor=0;iNeighbor<elementToNeighbors[el->getNum()].size();iNeighbor++){
+ int neighId=elementToNeighbors[el->getNum()][iNeighbor];
+ if(newGroupIds[neighId]==-1){
+ mapNewGroups[oldGroupIds[neighId]].push_back(allElements[neighId]);
+ newGroupIds[neighId]=-2;
+ }
+ }
+ }
+ }
+ }
+ if(!isOuterBufferLayer){
+ for(int iGroup=0;iGroup<getNbElementGroups();iGroup++){
+ dgGroupOfElements *elGroup=getElementGroup(iGroup);
+ for(int iElement=0;iElement<elGroup->getNbElements();iElement++){
+ MElement *el=elGroup->getElement(iElement);
+ if(localDt[el->getNum()]>=currentDt && localDt[el->getNum()]<currentDt*2){
+ if(newGroupIds[el->getNum()]==-1){
+ mapNewGroups[iGroup].push_back(el);
+ newGroupIds[el->getNum()]=-2;
+ }
+ }
+ }
+ }
+ lowerLevelGroupIdStart=currentNewGroupId;
+ for (std::map<int, std::vector<MElement *> >::iterator it = mapNewGroups.begin(); it !=mapNewGroups.end() ; ++it){
+ if(!it->second.empty()){
+ std::vector<MElement *>forBulk;
+ std::vector<MElement *>forInnerBuffer;
+ for(int i=0;i<it->second.size();i++){
+ MElement* el=it->second[i];
+ bool inInnerBuffer=false;
+ for(int iNeighbor=0;iNeighbor<elementToNeighbors[el->getNum()].size();iNeighbor++){
+ int neighId=elementToNeighbors[el->getNum()][iNeighbor];
+ if(newGroupIds[neighId]==-1){
+ inInnerBuffer=true;
+ continue;
+ }
+ }
+ if(inInnerBuffer){
+ forInnerBuffer.push_back(el);
+ }
+ else{
+ forBulk.push_back(el);
+ }
+ }
+ for(int i=0;i<forBulk.size();i++){
+ newGroupIds[it->second[i]->getNum()]=currentNewGroupId;
+ }
+ dgGroupOfElements *oldGroup=getElementGroup(it->first);
+ dgGroupOfElements *newGroup;
+ if(!forBulk.empty()){
+ newGroup=new dgGroupOfElements(forBulk,oldGroup->getOrder(),oldGroup->getGhostPartition());
+ newGroup->copyPrivateDataFrom(oldGroup);
+ newGroup->_multirateExponent=currentExponent;
+ newGroup->_multirateOuterBuffer=false;
+ newGroup->_multirateInnerBuffer=false;
+ newGroups.resize(currentNewGroupId+1);
+ newGroups[currentNewGroupId]=newGroup;
+ oldToNewGroupsMap[it->first]=currentNewGroupId;
+ currentNewGroupId++;
+ }
+ if(!forInnerBuffer.empty()){
+ newGroup=new dgGroupOfElements(forInnerBuffer,oldGroup->getOrder(),oldGroup->getGhostPartition());
+ newGroup->copyPrivateDataFrom(oldGroup);
+ newGroup->_multirateExponent=currentExponent;
+ newGroup->_multirateOuterBuffer=false;
+ newGroup->_multirateInnerBuffer=true;
+ newGroups.resize(currentNewGroupId+1);
+ newGroups[currentNewGroupId]=newGroup;
+ oldToNewGroupsMap[it->first]=currentNewGroupId;
+ currentNewGroupId++;
+ }
+ }
+ else
+ Msg::Info("Empty Group !!!!!!!!!!!!!!!!!!");
+ }
+ }
+ else{
+ for (std::map<int, std::vector<MElement *> >::iterator it = mapNewGroups.begin(); it !=mapNewGroups.end() ; ++it){
+ if(!it->second.empty()){
+ for(int i=0;i<it->second.size();i++){
+ newGroupIds[it->second[i]->getNum()]=currentNewGroupId;
+ }
+ dgGroupOfElements *oldGroup=getElementGroup(it->first);
+ dgGroupOfElements *newGroup=new dgGroupOfElements(it->second,oldGroup->getOrder(),oldGroup->getGhostPartition());
+ newGroup->copyPrivateDataFrom(oldGroup);
+ newGroup->_multirateExponent=currentExponent;
+ newGroup->_multirateOuterBuffer=true;
+ newGroup->_multirateInnerBuffer=false;
+ newGroups.resize(currentNewGroupId+1);
+ newGroups[currentNewGroupId]=newGroup;
+ oldToNewGroupsMap[it->first]=currentNewGroupId;
+ currentNewGroupId++;
+ }
+ else
+ Msg::Info("Empty Group !!!!!!!!!!!!!!!!!!");
+ }
+ }
+ lowerLevelGroupIdEnd=currentNewGroupId;
+ isOuterBufferLayer=!isOuterBufferLayer;
+ }
+ int count=0;
+ for(int i=0;i<newGroups.size();i++){
+ Msg::Info("%d New group # %d has %d elements",newGroups[i]->getMultirateExponent(),i,newGroups[i]->getNbElements());
+ if(newGroups[i]->getIsInnerMultirateBuffer())
+ Msg::Info("InnerBuffer");
+ else if(newGroups[i]->getIsOuterMultirateBuffer())
+ Msg::Info("OuterBuffer");
+ else
+ Msg::Info("Not Buffer");
+ count+=newGroups[i]->getNbElements();
+ }
+ Msg::Info("That makes a total of %d elements",count);
+ _elementGroups.clear();
+ _elementGroups=newGroups;
}
+
dgGroupCollection::dgGroupCollection()
{
_groupsOfElementsBuilt=false;_groupsOfInterfacesBuilt=false;
@@ -945,7 +1243,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("refDt","claw",NULL);
+ cm->setArgNames("refDt","claw","solution",NULL);
cm = cb->addMethod("getNbElementGroups", &dgGroupCollection::getNbElementGroups);
cm->setDescription("return the number of dgGroupOfElements");
cm = cb->addMethod("getNbFaceGroups", &dgGroupCollection::getNbFaceGroups);
diff --git a/Solver/dgGroupOfElements.h b/Solver/dgGroupOfElements.h
index 41c6f46450d61d9d12ed016a5ae609f363453bf8..f108b15c06dbf89e40f098d002084e4ecff9844b 100644
--- a/Solver/dgGroupOfElements.h
+++ b/Solver/dgGroupOfElements.h
@@ -22,6 +22,7 @@ class GModel;
class binding;
class dgConservationLaw;
+class dgDofContainer;
class dgElement {
@@ -40,10 +41,13 @@ public:
MElement *element() const { return _element;}
};
+class dgGroupOfFaces;
+
// store topological and geometrical data for 1 group for 1 discretisation
class dgGroupOfElements {
int _ghostPartition; // -1 : this is not a ghosted group, otherwise the id of the parent partition
+
// N elements in the group
std::vector<MElement*> _elements;
// inverse map that gives the index of an element in the group
@@ -79,7 +83,16 @@ class dgGroupOfElements {
// forbid the copy
// dgGroupOfElements (const dgGroupOfElements &e, int order) {}
// dgGroupOfElements & operator = (const dgGroupOfElements &e) {}
+protected:
+ int _multirateExponent;
+ bool _multirateInnerBuffer;
+ bool _multirateOuterBuffer;
public:
+
+ inline int getMultirateExponent() const {return _multirateExponent;}
+ inline int getIsInnerMultirateBuffer() const {return _multirateInnerBuffer;}
+ inline int getIsOuterMultirateBuffer() const {return _multirateOuterBuffer;}
+ inline int getIsMultirateBuffer()const {return _multirateInnerBuffer || _multirateOuterBuffer;}
dgGroupOfElements (const std::vector<MElement*> &e, int pOrder,int ghostPartition=-1);
virtual ~dgGroupOfElements ();
inline int getNbElements() const {return _elements.size();}
@@ -101,11 +114,13 @@ public:
inline const fullMatrix<double> getMapping (int iElement) const {return fullMatrix<double>(*_mapping, iElement, 1);}
inline int getOrder() const {return _order;}
inline int getGhostPartition() const {return _ghostPartition;}
+ void copyPrivateDataFrom(const dgGroupOfElements *from);
inline int getIndexOfElement (MElement *e) const {
std::map<MElement*,int>::const_iterator it = element_to_index.find(e);
if (it == element_to_index.end())return -1;
return it->second;
}
+ friend class dgGroupCollection;
};
class dgGroupOfFaces {
@@ -191,6 +206,8 @@ public:
//keep this outside the Algorithm because this is the only place where data overlap
void mapToInterface(int nFields, const fullMatrix<double> &vLeft, const fullMatrix<double> &vRight, fullMatrix<double> &v);
void mapFromInterface(int nFields, const fullMatrix<double> &v, fullMatrix<double> &vLeft, fullMatrix<double> &vRight);
+ void mapLeftFromInterface(int nFields, const fullMatrix<double> &v, fullMatrix<double> &vLeft);
+ void mapRightFromInterface(int nFields, const fullMatrix<double> &v, fullMatrix<double> &vRight);
const polynomialBasis * getPolynomialBasis() const {return _fsFace;}
};
@@ -205,7 +222,12 @@ class dgGroupCollection {
std::vector< std::vector<std::pair<int,int> > >_elementsToSend;
bool _groupsOfElementsBuilt;
bool _groupsOfInterfacesBuilt;
+
+ // Multirate stuff
+ int _dtMaxExponent;
+
public:
+ inline int getDtMaxExponent() const {return _dtMaxExponent;}
inline GModel* getModel() {return _model;}
inline int getNbElementGroups() const {return _elementGroups.size();}
inline int getNbFaceGroups() const {return _faceGroups.size();}
@@ -223,7 +245,7 @@ class dgGroupCollection {
void buildGroupsOfElements (GModel *model,int dimension, int order);
void buildGroupsOfInterfaces ();
- void splitGroupsForMultirate(double refDt,dgConservationLaw *claw);
+ void splitGroupsForMultirate(double refDt,dgConservationLaw *claw, dgDofContainer *solution);
void find (MElement *elementToFind, int &iGroup, int &ithElementOfGroup);
diff --git a/Solver/dgRungeKutta.cpp b/Solver/dgRungeKutta.cpp
index 2e6898738a1b3b466359ad04369c93af5b45a56d..a1b03cbdc3714b53888bd6adff1fcf517c0bfaea 100644
--- a/Solver/dgRungeKutta.cpp
+++ b/Solver/dgRungeKutta.cpp
@@ -203,6 +203,181 @@ double dgRungeKutta::nonDiagonalRK(const dgConservationLaw *claw,
return solution->norm();
}
+
+double dgRungeKutta::iterate43Multirate(const dgConservationLaw *claw, double dt, dgDofContainer *solution){
+ 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 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};
+
+ dgGroupCollection *groups=solution->getGroups();
+ int dtMaxExponent=groups->getDtMaxExponent();
+ int nSmallSteps=1;
+ {
+ int i=0;
+ while(i<dtMaxExponent){
+ nSmallSteps*=2;
+ i++;
+ }
+ }
+
+ int nbFields = claw->getNbFields();
+ dgDofContainer **K;
+ K=new dgDofContainer*[4];
+ for(int i=0;i<4;i++){
+ K[i]=new dgDofContainer(groups,nbFields);
+ K[i]->scale(0.);
+ }
+ dgDofContainer Unp (groups,nbFields);
+ dgDofContainer tmp (groups,nbFields);
+ dgDofContainer resd(groups,nbFields);
+
+// axpy for each group to build the right input vector for the residual
+// At each stage, update the residual for a std::vector of groups, using residualForSomeGroups
+
+// K0 for everyone
+// K1 K2 K3 for the smallestDt group and its corresponding buffer
+// smallestDt can be updated for the first smallset time step
+// K4 for the smallest Dt and the buffer and K1 for the larger Dt
+// K0 for the smallest Dt, K5 for the buffer, and K2 for the larger Dt
+// K1 K2 K3 for the smallestDt and K6 K7 K8 for the buffer
+// K4 for the smallest Dt and K9 for the buffer, K3 for the larger Dt
+// update for everyone
+
+
+
+/*
+ for(int iStep=0;iStep<nSmallSteps;iStep++){
+ for(int iGroup=0; iGroup < groups->getNbElementGroups(); iGroup++) {
+ for(int iStage=0;iStage<4;iStage++){
+ dgGroupOfElements *group=groups->getElementGroup(iGroup);
+ int groupDtExponent=group->getMultirateExponent();
+ }
+ tmp.scale(0.0);
+ tmp.axpy(*solution);
+ for(int j=0;j<i;j++){
+ if(fabs(A(i,j))>1e-12){
+ tmp.axpy(*K[j],dt*A(i,j));
+ }
+ }
+ if (_limiter) _limiter->apply(&tmp);
+ dgAlgorithm::residual(*claw,*groups,tmp,resd);
+
+ // Multiply the spatial residual by the inverse of the mass matrix
+ int nbNodes = group->getNbNodes();
+ for(int j=0;j<group->getNbElements();j++) {
+ residuEl.setAsProxy(resd.getGroupProxy(k),j*nbFields,nbFields);
+ KEl.setAsProxy(K[i]->getGroupProxy(k),j*nbFields,nbFields);
+ iMassEl.setAsProxy(group->getInverseMassMatrix(),j*nbNodes,nbNodes);
+ iMassEl.mult(residuEl,KEl);
+ }
+ }
+ Unp.axpy(*K[i],dt*b[i]);
+ }
+ }
+ */
+/*
+ int nStages=10;
+ // Small step RK43
+ double _A1[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}
+ };
+ double b1[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 c1[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 };
+
+ // Big step RK43
+ double _A2[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},
+ };
+ double b2[10]={1.0/6.0, 0 , 0 , 0 , 1.0/3.0,1.0/3.0, 0 , 0 , 0 , 1.0/6.0};
+ double c2[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 };
+ fullMatrix<double> A1(10,10);
+ fullMatrix<double> A2(10,10);
+ for(int i=0;i<10;i++){
+ for(int j=0;j<10;j++){
+ A1(i,j)=_A1[i][j];
+ A2(i,j)=_A2[i][j];
+ }
+ }
+
+
+ dgGroupCollection *groups=solution->getGroups();
+ fullMatrix<double> residuEl, KEl;
+ fullMatrix<double> iMassEl;
+
+ int nbFields = claw->getNbFields();
+
+ dgDofContainer **K;
+ K=new dgDofContainer*[nStages];
+ for(int i=0;i<nStages;i++){
+ K[i]=new dgDofContainer(groups,nbFields);
+ K[i]->scale(0.);
+ }
+ dgDofContainer Unp (groups,nbFields);
+ dgDofContainer tmp (groups,nbFields);
+ dgDofContainer resd(groups,nbFields);
+
+ Unp.scale(0.0);
+ Unp.axpy(*solution);
+ for(int i=0; i<nStages;i++){
+ tmp.scale(0.0);
+ tmp.axpy(*solution);
+ for(int j=0;j<i;j++){
+ if(fabs(A(i,j))>1e-12){
+ tmp.axpy(*K[j],dt*A(i,j));
+ }
+ }
+ if (_limiter) _limiter->apply(&tmp);
+ dgAlgorithm::residual(*claw,*groups,tmp,resd);
+
+ // Multiply the spatial residual by the inverse of the mass matrix
+ for(int k=0; k < groups->getNbElementGroups(); k++) {
+ dgGroupOfElements *group = groups->getElementGroup(k);
+ int nbNodes = group->getNbNodes();
+ for(int j=0;j<group->getNbElements();j++) {
+ residuEl.setAsProxy(resd.getGroupProxy(k),j*nbFields,nbFields);
+ KEl.setAsProxy(K[i]->getGroupProxy(k),j*nbFields,nbFields);
+ iMassEl.setAsProxy(group->getInverseMassMatrix(),j*nbNodes,nbNodes);
+ iMassEl.mult(residuEl,KEl);
+ }
+ }
+ Unp.axpy(*K[i],dt*b[i]);
+ }
+ if (_limiter) _limiter->apply(&Unp);
+ solution->scale(0.);
+ solution->axpy(Unp);
+
+ for(int i=0;i<nStages;i++){
+ delete K[i];
+ }
+ delete []K;
+ return solution->norm();
+*/
+}
+
+
dgRungeKutta::dgRungeKutta():_limiter(NULL) {}
@@ -235,7 +410,229 @@ void dgRungeKutta::registerBindings(binding *b) {
cm = cb->addMethod("iterateRKFehlberg45",&dgRungeKutta::iterateRKFehlberg45);
cm->setArgNames("law","dt","solution",NULL);
cm->setDescription("update solution by doing a six stage fifth order Runge-Kutta-Fehlberg step of time step dt for the conservation law");
+ cm = cb->addMethod("iterate43Multirate",&dgRungeKutta::iterate43Multirate);
+ cm->setArgNames("law","dt","solution",NULL);
+ cm->setDescription("update solution by doing a multirate RK43 (from Schlegel et al. Journal of Computational and Applied Mathematics, 2009) step of base time step dt for the conservation law");
cm = cb->addMethod("setLimiter",&dgRungeKutta::setLimiter);
cm->setArgNames("limiter",NULL);
cm->setDescription("if a limiter is set, it is applied after each RK stage");
}
+
+dgRungeKuttaMultirate::dgRungeKuttaMultirate(dgGroupCollection* gc,dgConservationLaw *law){
+ _law=law;
+ _residualVolume=new dgResidualVolume(*_law);
+ _residualInterface=new dgResidualInterface(*_law);
+ _K=new dgDofContainer*[10];
+ for(int i=0;i>10;i++){
+ _K[i]=new dgDofContainer(gc,_law->getNbFields());
+ _K[i]->setAll(0.0);
+ }
+ _currentInput=new dgDofContainer(gc,_law->getNbFields());
+ _currentInput->setAll(0.0);
+ _minExponent=INT_MAX;
+ _maxExponent=0;
+ for(int iGroup=0;iGroup<gc->getNbElementGroups();iGroup++){
+ dgGroupOfElements *ge=gc->getElementGroup(iGroup);
+ _minExponent=std::min(_minExponent,ge->getMultirateExponent());
+ _maxExponent=std::max(_maxExponent,ge->getMultirateExponent());
+ _innerBufferGroupsOfElements.resize(_maxExponent+1);
+ _outerBufferGroupsOfElements.resize(_maxExponent+1);
+ _bulkGroupsOfElements.resize(_maxExponent+1);
+ if(ge->getIsInnerMultirateBuffer()){
+ _innerBufferGroupsOfElements[ge->getMultirateExponent()].first.push_back(ge);
+ }
+ else if(ge->getIsOuterMultirateBuffer()){
+ _outerBufferGroupsOfElements[ge->getMultirateExponent()].first.push_back(ge);
+ }
+ else{
+ _bulkGroupsOfElements[ge->getMultirateExponent()].first.push_back(ge);
+ }
+ }
+ for(int iGroup=0;iGroup<gc->getNbFaceGroups();iGroup++){
+ dgGroupOfFaces *gf=gc->getFaceGroup(iGroup);
+ const dgGroupOfElements *ge[2];
+ ge[0]=&gf->getGroupLeft();
+ ge[1]=&gf->getGroupRight();
+ if(ge[0]==ge[1]){
+ _bulkGroupsOfElements[ge[0]->getMultirateExponent()].second.push_back(gf);
+ }
+ for(int i=0;i<2;i++){
+ if(ge[i]->getIsInnerMultirateBuffer()){
+ _innerBufferGroupsOfElements[ge[i]->getMultirateExponent()].second.push_back(gf);
+ }
+ else if(ge[i]->getIsOuterMultirateBuffer()){
+ _outerBufferGroupsOfElements[ge[i]->getMultirateExponent()].second.push_back(gf);
+ }else{
+ _bulkGroupsOfElements[ge[i]->getMultirateExponent()].second.push_back(gf);
+ }
+ }
+ }
+}
+dgRungeKuttaMultirate::~dgRungeKuttaMultirate(){
+ for(int i=0;i>10;i++){
+ delete _K[i];
+ }
+ delete []_K;
+}
+
+void dgRungeKuttaMultirate::computeK(int iK,int exponent,bool isBuffer){
+ 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}
+ };
+ 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 };
+
+ // 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},
+ };
+ 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(exponent>_maxExponent)
+ return;
+ // compute K[iK] for this exponent and buffer
+ _currentInput->scale(0.0);
+ if(isBuffer){
+ std::vector<dgGroupOfElements *>gVi=_innerBufferGroupsOfElements[exponent].first;
+ std::vector<dgGroupOfElements *>gVo=_outerBufferGroupsOfElements[exponent].first;
+ _currentInput->axpy(gVi,*_solution,1.0);
+ _currentInput->axpy(gVo,*_solution,1.0);
+ for(int i=0;i<iK;i++){
+ _currentInput->axpy(gVi,*_K[i],_AInner[iK][i]);
+ _currentInput->axpy(gVo,*_K[i],_AOuter[iK][i]);
+ }
+ }
+ else{
+ std::vector<dgGroupOfElements *>gV=_bulkGroupsOfElements[exponent].first;
+ _currentInput->axpy(gV,*_solution,1.0);
+ for(int i=0;i<iK;i++){
+ _currentInput->axpy(gV,*_K[i],_A[iK][i]);
+ }
+ }
+ if(!isBuffer){
+ switch(iK){
+ case 0:
+ for(int i=0;i<4;i++){
+ computeK(i,exponent+1,true);
+ }
+ break;
+ case 1:
+ computeK(4,exponent+1,true);
+ break;
+ case 2:
+ for(int i=5;i<9;i++){
+ computeK(i,exponent+1,true);
+ }
+ break;
+ case 3:
+ computeK(9,exponent+1,true);
+ break;
+ }
+ }
+ else{
+ if( (iK%5)<4 ){
+ computeK(iK%5,exponent+1,false);
+ }
+ }
+ if( (iK==3 && !isBuffer) || (iK==9 && isBuffer) ){
+ updateSolution(exponent,isBuffer);
+ }
+}
+void dgRungeKuttaMultirate::updateSolution(int exponent,bool isBuffer){}
+void dgRungeKuttaMultirate::computeResidual(int exponent,bool isBuffer,dgDofContainer *solution, dgDofContainer *residual){
+ if(isBuffer){
+ std::vector<dgGroupOfElements *>*vei=&_innerBufferGroupsOfElements[exponent].first;
+ std::vector<dgGroupOfElements *>*veo=&_outerBufferGroupsOfElements[exponent].first;
+ for(int i=0;i<vei->size();i++){
+ _residualVolume->computeAndMap1Group(*(*vei)[i], *solution, *residual);
+ }
+ for(int i=0;i<veo->size();i++){
+ _residualVolume->computeAndMap1Group(*(*veo)[i], *solution, *residual);
+ }
+ std::vector<dgGroupOfFaces *>* vfi=&_innerBufferGroupsOfElements[exponent].second;
+ std::vector<dgGroupOfFaces *>* vfo=&_outerBufferGroupsOfElements[exponent].second;
+ std::set<dgGroupOfFaces *>gOFSet;
+ gOFSet.insert(vfo->begin(),vfo->end());
+ for(std::set<dgGroupOfFaces *>::iterator it=gOFSet.begin();it!=gOFSet.end();it++){
+ dgGroupOfFaces *faces=*it;
+ const dgGroupOfElements *gL=&faces->getGroupLeft();
+ const dgGroupOfElements *gR=&faces->getGroupRight();
+ fullMatrix<double> solInterface(faces->getNbNodes(),faces->getNbElements()*2*_law->getNbFields());
+ fullMatrix<double> residuInterface(faces->getNbNodes(),faces->getNbElements()*2*_law->getNbFields());
+ faces->mapToInterface(_law->getNbFields(), solution->getGroupProxy(gL), solution->getGroupProxy(gR), solInterface);
+ _residualInterface->compute1Group(*faces,solInterface,solution->getGroupProxy(gL), solution->getGroupProxy(gR),residuInterface);
+ if(gL->getMultirateExponent()==exponent && gL->getIsMultirateBuffer())
+ faces->mapLeftFromInterface(_law->getNbFields(), residuInterface,residual->getGroupProxy(gL));
+ if(gL->getMultirateExponent()==exponent && gL->getIsMultirateBuffer())
+ faces->mapRightFromInterface(_law->getNbFields(), residuInterface,residual->getGroupProxy(gR));
+ }
+ }
+ else{
+ std::vector<dgGroupOfElements *> *ve=&_bulkGroupsOfElements[exponent].first;
+ for(int i=0;i<ve->size();i++){
+ _residualVolume->computeAndMap1Group(*(*ve)[i], *solution, *residual);
+ }
+ std::vector<dgGroupOfFaces *> *vf=&_bulkGroupsOfElements[exponent].second;
+ for(int i=0;i<vf->size();i++){
+ dgGroupOfFaces *faces=(*vf)[i];
+ const dgGroupOfElements *gL=&faces->getGroupLeft();
+ const dgGroupOfElements *gR=&faces->getGroupRight();
+ fullMatrix<double> solInterface(faces->getNbNodes(),faces->getNbElements()*2*_law->getNbFields());
+ fullMatrix<double> residuInterface(faces->getNbNodes(),faces->getNbElements()*2*_law->getNbFields());
+ faces->mapToInterface(_law->getNbFields(), solution->getGroupProxy(gL), solution->getGroupProxy(gR), solInterface);
+ _residualInterface->compute1Group(*faces,solInterface,solution->getGroupProxy(gL), solution->getGroupProxy(gR),residuInterface);
+ if(gL->getMultirateExponent()==exponent && !gL->getIsMultirateBuffer())
+ faces->mapLeftFromInterface(_law->getNbFields(), residuInterface,residual->getGroupProxy(gL));
+ if(gL->getMultirateExponent()==exponent && !gL->getIsMultirateBuffer())
+ faces->mapRightFromInterface(_law->getNbFields(), residuInterface,residual->getGroupProxy(gR));
+ }
+ }
+}
+
+double dgRungeKuttaMultirate::iterate(double dt, dgDofContainer *solution){
+ _solution=solution;
+ computeK(0,0,false);
+ computeK(1,0,false);
+ computeK(2,0,false);
+ computeK(3,0,false);
+ updateSolution(0,false);
+
+ return solution->norm();
+}
+
+void dgRungeKuttaMultirate::registerBindings(binding *b) {
+ classBinding *cb = b->addClass<dgRungeKuttaMultirate>("dgRungeKuttaMultirate");
+ cb->setDescription("Multirate explicit Runge-Kutta");
+ methodBinding *cm;
+ cm = cb->setConstructor<dgRungeKuttaMultirate,dgGroupCollection *,dgConservationLaw*>();
+ cm->setArgNames("groupCollection","law",NULL);
+ cm->setDescription("A new multirate explicit runge kutta, pass parameters to the iterate function");
+ cm = cb->addMethod("iterate",&dgRungeKuttaMultirate::iterate);
+ cm->setArgNames("dt","solution",NULL);
+ cm->setDescription("update solution by doing a multirate RK43 (from Schlegel et al. Journal of Computational and Applied Mathematics, 2009) step of base time step dt for the conservation law");
+ cb->setParentClass<dgRungeKutta>();
+}
diff --git a/Solver/dgRungeKutta.h b/Solver/dgRungeKutta.h
index 09a19cd2a963ec216a8286b098c4724462a7262e..4168dd6ad0fb61a6c0d4629969c25b5f4bfcc62e 100644
--- a/Solver/dgRungeKutta.h
+++ b/Solver/dgRungeKutta.h
@@ -1,8 +1,15 @@
#ifndef _DG_RUNGE_KUTTA_H_
#define _DG_RUNGE_KUTTA_H_
+#include <vector>
+#include <map>
class dgConservationLaw;
class dgDofContainer;
class dgLimiter;
+class dgGroupCollection;
+class dgGroupOfElements;
+class dgGroupOfFaces;
+class dgResidualVolume;
+class dgResidualInterface;
class binding;
#include "fullMatrix.h"
class dgRungeKutta {
@@ -18,7 +25,33 @@ class dgRungeKutta {
double iterate44ThreeEight(const dgConservationLaw *claw, double dt, dgDofContainer *solution);
double iterate43SchlegelJCAM2009(const dgConservationLaw *claw, double dt, dgDofContainer *solution);
double iterateRKFehlberg45(const dgConservationLaw *claw, double dt, dgDofContainer *solution);
+ double iterate43Multirate(const dgConservationLaw *claw, double dt, dgDofContainer *solution);
static void registerBindings (binding *b);
dgRungeKutta();
};
+
+class dgRungeKuttaMultirate: public dgRungeKutta{
+ private:
+ int _maxExponent;
+ int _minExponent;
+ dgConservationLaw *_law;
+ dgDofContainer *_solution;
+ dgDofContainer **_K;
+ dgDofContainer *_currentInput;
+ dgResidualVolume *_residualVolume;
+ dgResidualInterface *_residualInterface;
+ 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
+ void computeK(int iK,int exponent,bool isBuffer);
+ void updateSolution(int exponent,bool isBuffer);
+ void computeResidual(int exponent,bool isBuffer,dgDofContainer *solution, dgDofContainer *residual);
+ public:
+ dgRungeKuttaMultirate(dgGroupCollection *gc,dgConservationLaw* law);
+ ~dgRungeKuttaMultirate();
+
+
+ double iterate(double dt, dgDofContainer *solution);
+ static void registerBindings(binding *b);
+};
#endif