Skip to content
Snippets Groups Projects
Commit 0386a196 authored by Richard Comblen's avatar Richard Comblen
Browse files

Multirate on the way, almost finished, test needed

parent c58fd0ba
No related branches found
No related tags found
No related merge requests found
...@@ -343,6 +343,7 @@ binding::binding(){ ...@@ -343,6 +343,7 @@ binding::binding(){
dgPerfectGasLaw2dRegisterBindings(this); dgPerfectGasLaw2dRegisterBindings(this);
dgResidual::registerBindings(this); dgResidual::registerBindings(this);
dgRungeKutta::registerBindings(this); dgRungeKutta::registerBindings(this);
dgRungeKuttaMultirate::registerBindings(this);
dgSlopeLimiterRegisterBindings(this); dgSlopeLimiterRegisterBindings(this);
dgSystemOfEquations::registerBindings(this); dgSystemOfEquations::registerBindings(this);
fullMatrix<double>::registerBindings(this); fullMatrix<double>::registerBindings(this);
......
...@@ -996,7 +996,7 @@ static void checkHighOrderTetrahedron(const char* cc, GModel *m, ...@@ -996,7 +996,7 @@ static void checkHighOrderTetrahedron(const char* cc, GModel *m,
minGGlob, avg / (count ? count : 1)); 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) static void printJacobians(GModel *m, const char *nm)
{ {
......
...@@ -18,8 +18,8 @@ Plane Surface(11) = {10}; ...@@ -18,8 +18,8 @@ Plane Surface(11) = {10};
Physical Surface("sprut") = {11, 9}; Physical Surface("sprut") = {11, 9};
Physical Line("Walls") = {5, 7, 3, 4, 1}; Physical Line("Walls") = {5, 7, 3, 4, 1};
Physical Line("Top") = {6}; Physical Line("Top") = {6};
Transfinite Line {5, 7, 1, 3} = 50 Using Progression 1; //Recombine Surface {9, 11};
Transfinite Line {6, 4, 2} = 25 Using Progression 1;
Transfinite Surface {9} Alternated; Field[1] = MathEval;
Transfinite Surface {11} Alternated; Field[1].F = "0.01*1+0.01*1.7*y";
Recombine Surface {9, 11}; Background Field = 1;
...@@ -19,6 +19,64 @@ ...@@ -19,6 +19,64 @@
// works for any number of groups // 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) void dgAlgorithm::computeElementaryTimeSteps ( //dofManager &dof, // the DOF manager (maybe useless here)
const dgConservationLaw &claw, // the conservation law const dgConservationLaw &claw, // the conservation law
const dgGroupOfElements & group, const dgGroupOfElements & group,
......
...@@ -166,6 +166,14 @@ void dgDofContainer::axpy(dgDofContainer &x, double a) ...@@ -166,6 +166,14 @@ void dgDofContainer::axpy(dgDofContainer &x, double a)
_data->axpy(*x._data,a); _data->axpy(*x._data,a);
_ghostData->axpy(*x._ghostData,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 dgDofContainer::norm() {
double localNorm = _data->norm(); double localNorm = _data->norm();
...@@ -213,6 +221,70 @@ void dgDofContainer::L2Projection(std::string functionName){ ...@@ -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) void dgDofContainer::exportMsh(const std::string name)
{ {
...@@ -236,7 +308,7 @@ void dgDofContainer::exportMsh(const std::string name) ...@@ -236,7 +308,7 @@ void dgDofContainer::exportMsh(const std::string name)
fprintf(f,"1\n"); fprintf(f,"1\n");
fprintf(f,"%d\n", _mshStep); // should print actual time here fprintf(f,"%d\n", _mshStep); // should print actual time here
fprintf(f,"%d\n", Msg::GetCommSize() > 1 ? 4 : 3); 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()); if(Msg::GetCommSize() > 1) fprintf(f,"%d\n", Msg::GetCommRank());
for (int i=0;i < _groups.getNbElementGroups() ;i++){ for (int i=0;i < _groups.getNbElementGroups() ;i++){
dgGroupOfElements *group = _groups.getElementGroup(i); dgGroupOfElements *group = _groups.getElementGroup(i);
...@@ -300,4 +372,6 @@ void dgDofContainer::registerBindings(binding *b){ ...@@ -300,4 +372,6 @@ void dgDofContainer::registerBindings(binding *b){
cm = cb->addMethod("exportMsh",&dgDofContainer::exportMsh); cm = cb->addMethod("exportMsh",&dgDofContainer::exportMsh);
cm->setDescription("Export the dof for gmsh visualization"); cm->setDescription("Export the dof for gmsh visualization");
cm->setArgNames("filename", NULL); cm->setArgNames("filename", NULL);
cm = cb->addMethod("exportGroupIdMsh",&dgDofContainer::exportGroupIdMsh);
cm->setDescription("Export the group ids for gmsh visualization");
} }
...@@ -28,6 +28,7 @@ public: ...@@ -28,6 +28,7 @@ public:
void scale(double f); void scale(double f);
double norm(); double norm();
void axpy(dgDofContainer &x, double a=1.); 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(int gId){ return *(_dataProxys[gId]); }
inline fullMatrix<double> &getGroupProxy(const dgGroupOfElements* g){ return *(_dataProxys[_groupId[g]]); } inline fullMatrix<double> &getGroupProxy(const dgGroupOfElements* g){ return *(_dataProxys[_groupId[g]]); }
dgDofContainer (dgGroupCollection *groups, int nbFields); dgDofContainer (dgGroupCollection *groups, int nbFields);
...@@ -41,6 +42,7 @@ public: ...@@ -41,6 +42,7 @@ public:
void L2Projection(std::string functionName); void L2Projection(std::string functionName);
void Mesh2Mesh_L2Projection(dgDofContainer &other); void Mesh2Mesh_L2Projection(dgDofContainer &other);
void exportMsh(const std::string name); void exportMsh(const std::string name);
void exportGroupIdMsh();
static void registerBindings(binding *b); static void registerBindings(binding *b);
inline dgGroupCollection *getGroups() { return &_groups; } inline dgGroupCollection *getGroups() { return &_groups; }
......
#include <iostream>
#include "GmshConfig.h" #include "GmshConfig.h"
#include "GmshMessage.h"
#include "dgGroupOfElements.h" #include "dgGroupOfElements.h"
#include "dgConservationLaw.h" #include "dgConservationLaw.h"
#include "dgDofContainer.h" #include "dgDofContainer.h"
...@@ -131,6 +131,10 @@ dgGroupOfElements::dgGroupOfElements(const std::vector<MElement*> &e, ...@@ -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(){ dgGroupOfElements::~dgGroupOfElements(){
delete _integration; delete _integration;
delete _redistributionFluxes[0]; delete _redistributionFluxes[0];
...@@ -622,6 +626,36 @@ void dgGroupOfFaces::mapFromInterface ( int nFields, ...@@ -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, void dgAlgorithm::buildGroupsOfFaces(GModel *model, int dim, int order,
...@@ -675,8 +709,11 @@ void dgGroupCollection::buildGroupsOfElements(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 // 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){ 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() { ...@@ -745,7 +782,7 @@ void dgGroupCollection::buildGroupsOfInterfaces() {
for (int i=0;i<_elementGroups.size();i++){ 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) { switch(dim) {
case 1 : { case 1 : {
std::map<const std::string, std::set<MVertex*> >::iterator mapIt; std::map<const std::string, std::set<MVertex*> >::iterator mapIt;
...@@ -782,7 +819,11 @@ void dgGroupCollection::buildGroupsOfInterfaces() { ...@@ -782,7 +819,11 @@ void dgGroupCollection::buildGroupsOfInterfaces() {
} }
} }
// create a group of faces for every face that is common to elements of the same group // 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 // create a groupe of d
for (int j=i+1;j<_elementGroups.size();j++){ for (int j=i+1;j<_elementGroups.size();j++){
dgGroupOfFaces *gof = new dgGroupOfFaces(*_elementGroups[i],*_elementGroups[j],order); dgGroupOfFaces *gof = new dgGroupOfFaces(*_elementGroups[i],*_elementGroups[j],order);
...@@ -795,7 +836,11 @@ void dgGroupCollection::buildGroupsOfInterfaces() { ...@@ -795,7 +836,11 @@ void dgGroupCollection::buildGroupsOfInterfaces() {
//create ghost groups //create ghost groups
for(int i=0;i<Msg::GetCommSize();i++){ 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){ 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 //create face group for ghostGroups
...@@ -875,29 +920,282 @@ void dgGroupCollection::buildGroupsOfInterfaces() { ...@@ -875,29 +920,282 @@ void dgGroupCollection::buildGroupsOfInterfaces() {
// Split the groups of elements depending on their local time step // Split the groups of elements depending on their local time step
void dgGroupCollection::splitGroupsForMultirate(double refDt,dgConservationLaw *claw){ void dgGroupCollection::splitGroupsForMultirate(double dtRef,dgConservationLaw *claw, dgDofContainer *solution){
// 1) find the element with the smallest stable time step Msg::Info("Splitting Groups for multirate time stepping");
double sr = DBL_MAX; int maxNumElems=getElementGroup(0)->getElement(0)->getGlobalNumber()+1;
dgDofContainer *solution = new dgDofContainer((this),claw->getNbFields()); std::vector<int>oldGroupIds;
solution->scale(0.); 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++){
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]);
}
}
std::vector<double>localDt;
localDt.resize(maxNumElems);
for (int i=0;i<getNbElementGroups();i++){ for (int i=0;i<getNbElementGroups();i++){
std::vector<double> DTS; dgGroupOfElements *elGroup=getElementGroup(i);
dgAlgorithm::computeElementaryTimeSteps(*claw, *getElementGroup(i), solution->getGroupProxy(i), DTS); for(int j=0;j<DTS[i].size();j++){
for (int k=0;k<DTS.size();k++){ MElement *el=elGroup->getElement(j);
std::cout << "SR: " << DTS[k] << std::endl; localDt[el->getNum()]=DTS[i][j];
sr = std::min(sr,DTS[k]);
} }
} }
delete solution; delete []DTS;
#ifdef HAVE_MPI #ifdef HAVE_MPI
double sr_min; double dtMin_min;
MPI_Allreduce((void *)&sr, &sr_min, 1, MPI_DOUBLE, MPI_MIN, MPI_COMM_WORLD); MPI_Allreduce((void *)&dtMin, &dtMin_min, 1, MPI_DOUBLE, MPI_MIN, MPI_COMM_WORLD);
sr=sr_min; dtMin=dtMin_min;
double dtMax_max;
MPI_Allreduce((void *)&dtMax, &dtMax_max, 1, MPI_DOUBLE, MPI_MAX, MPI_COMM_WORLD);
dtMax=dtMax_max;
#endif #endif
// sr is the time step for the most constrained element. // dtMin is the time step for the most constrained element.
std::cout << "SR: " << sr << std::endl; 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() dgGroupCollection::dgGroupCollection()
{ {
...@@ -945,7 +1243,7 @@ void dgGroupCollection::registerBindings(binding *b){ ...@@ -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->setDescription("Build the group of interfaces, i.e. boundary interfaces and inter-element interfaces");
cm = cb->addMethod("splitGroupsForMultirate",&dgGroupCollection::splitGroupsForMultirate); cm = cb->addMethod("splitGroupsForMultirate",&dgGroupCollection::splitGroupsForMultirate);
cm->setDescription("Split the groups according to their own stable time step"); 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 = cb->addMethod("getNbElementGroups", &dgGroupCollection::getNbElementGroups);
cm->setDescription("return the number of dgGroupOfElements"); cm->setDescription("return the number of dgGroupOfElements");
cm = cb->addMethod("getNbFaceGroups", &dgGroupCollection::getNbFaceGroups); cm = cb->addMethod("getNbFaceGroups", &dgGroupCollection::getNbFaceGroups);
......
...@@ -22,6 +22,7 @@ class GModel; ...@@ -22,6 +22,7 @@ class GModel;
class binding; class binding;
class dgConservationLaw; class dgConservationLaw;
class dgDofContainer;
class dgElement { class dgElement {
...@@ -40,10 +41,13 @@ public: ...@@ -40,10 +41,13 @@ public:
MElement *element() const { return _element;} MElement *element() const { return _element;}
}; };
class dgGroupOfFaces;
// store topological and geometrical data for 1 group for 1 discretisation // store topological and geometrical data for 1 group for 1 discretisation
class dgGroupOfElements { class dgGroupOfElements {
int _ghostPartition; // -1 : this is not a ghosted group, otherwise the id of the parent partition int _ghostPartition; // -1 : this is not a ghosted group, otherwise the id of the parent partition
// N elements in the group // N elements in the group
std::vector<MElement*> _elements; std::vector<MElement*> _elements;
// inverse map that gives the index of an element in the group // inverse map that gives the index of an element in the group
...@@ -79,7 +83,16 @@ class dgGroupOfElements { ...@@ -79,7 +83,16 @@ class dgGroupOfElements {
// forbid the copy // forbid the copy
// dgGroupOfElements (const dgGroupOfElements &e, int order) {} // dgGroupOfElements (const dgGroupOfElements &e, int order) {}
// dgGroupOfElements & operator = (const dgGroupOfElements &e) {} // dgGroupOfElements & operator = (const dgGroupOfElements &e) {}
protected:
int _multirateExponent;
bool _multirateInnerBuffer;
bool _multirateOuterBuffer;
public: 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); dgGroupOfElements (const std::vector<MElement*> &e, int pOrder,int ghostPartition=-1);
virtual ~dgGroupOfElements (); virtual ~dgGroupOfElements ();
inline int getNbElements() const {return _elements.size();} inline int getNbElements() const {return _elements.size();}
...@@ -101,11 +114,13 @@ public: ...@@ -101,11 +114,13 @@ public:
inline const fullMatrix<double> getMapping (int iElement) const {return fullMatrix<double>(*_mapping, iElement, 1);} inline const fullMatrix<double> getMapping (int iElement) const {return fullMatrix<double>(*_mapping, iElement, 1);}
inline int getOrder() const {return _order;} inline int getOrder() const {return _order;}
inline int getGhostPartition() const {return _ghostPartition;} inline int getGhostPartition() const {return _ghostPartition;}
void copyPrivateDataFrom(const dgGroupOfElements *from);
inline int getIndexOfElement (MElement *e) const { inline int getIndexOfElement (MElement *e) const {
std::map<MElement*,int>::const_iterator it = element_to_index.find(e); std::map<MElement*,int>::const_iterator it = element_to_index.find(e);
if (it == element_to_index.end())return -1; if (it == element_to_index.end())return -1;
return it->second; return it->second;
} }
friend class dgGroupCollection;
}; };
class dgGroupOfFaces { class dgGroupOfFaces {
...@@ -191,6 +206,8 @@ public: ...@@ -191,6 +206,8 @@ public:
//keep this outside the Algorithm because this is the only place where data overlap //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 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 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;} const polynomialBasis * getPolynomialBasis() const {return _fsFace;}
}; };
...@@ -205,7 +222,12 @@ class dgGroupCollection { ...@@ -205,7 +222,12 @@ class dgGroupCollection {
std::vector< std::vector<std::pair<int,int> > >_elementsToSend; std::vector< std::vector<std::pair<int,int> > >_elementsToSend;
bool _groupsOfElementsBuilt; bool _groupsOfElementsBuilt;
bool _groupsOfInterfacesBuilt; bool _groupsOfInterfacesBuilt;
// Multirate stuff
int _dtMaxExponent;
public: public:
inline int getDtMaxExponent() const {return _dtMaxExponent;}
inline GModel* getModel() {return _model;} inline GModel* getModel() {return _model;}
inline int getNbElementGroups() const {return _elementGroups.size();} inline int getNbElementGroups() const {return _elementGroups.size();}
inline int getNbFaceGroups() const {return _faceGroups.size();} inline int getNbFaceGroups() const {return _faceGroups.size();}
...@@ -223,7 +245,7 @@ class dgGroupCollection { ...@@ -223,7 +245,7 @@ class dgGroupCollection {
void buildGroupsOfElements (GModel *model,int dimension, int order); void buildGroupsOfElements (GModel *model,int dimension, int order);
void buildGroupsOfInterfaces (); void buildGroupsOfInterfaces ();
void splitGroupsForMultirate(double refDt,dgConservationLaw *claw); void splitGroupsForMultirate(double refDt,dgConservationLaw *claw, dgDofContainer *solution);
void find (MElement *elementToFind, int &iGroup, int &ithElementOfGroup); void find (MElement *elementToFind, int &iGroup, int &ithElementOfGroup);
......
...@@ -203,6 +203,181 @@ double dgRungeKutta::nonDiagonalRK(const dgConservationLaw *claw, ...@@ -203,6 +203,181 @@ double dgRungeKutta::nonDiagonalRK(const dgConservationLaw *claw,
return solution->norm(); 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) {} dgRungeKutta::dgRungeKutta():_limiter(NULL) {}
...@@ -235,7 +410,229 @@ void dgRungeKutta::registerBindings(binding *b) { ...@@ -235,7 +410,229 @@ void dgRungeKutta::registerBindings(binding *b) {
cm = cb->addMethod("iterateRKFehlberg45",&dgRungeKutta::iterateRKFehlberg45); cm = cb->addMethod("iterateRKFehlberg45",&dgRungeKutta::iterateRKFehlberg45);
cm->setArgNames("law","dt","solution",NULL); 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->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 = cb->addMethod("setLimiter",&dgRungeKutta::setLimiter);
cm->setArgNames("limiter",NULL); cm->setArgNames("limiter",NULL);
cm->setDescription("if a limiter is set, it is applied after each RK stage"); 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>();
}
#ifndef _DG_RUNGE_KUTTA_H_ #ifndef _DG_RUNGE_KUTTA_H_
#define _DG_RUNGE_KUTTA_H_ #define _DG_RUNGE_KUTTA_H_
#include <vector>
#include <map>
class dgConservationLaw; class dgConservationLaw;
class dgDofContainer; class dgDofContainer;
class dgLimiter; class dgLimiter;
class dgGroupCollection;
class dgGroupOfElements;
class dgGroupOfFaces;
class dgResidualVolume;
class dgResidualInterface;
class binding; class binding;
#include "fullMatrix.h" #include "fullMatrix.h"
class dgRungeKutta { class dgRungeKutta {
...@@ -18,7 +25,33 @@ class dgRungeKutta { ...@@ -18,7 +25,33 @@ class dgRungeKutta {
double iterate44ThreeEight(const dgConservationLaw *claw, double dt, dgDofContainer *solution); double iterate44ThreeEight(const dgConservationLaw *claw, double dt, dgDofContainer *solution);
double iterate43SchlegelJCAM2009(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 iterateRKFehlberg45(const dgConservationLaw *claw, double dt, dgDofContainer *solution);
double iterate43Multirate(const dgConservationLaw *claw, double dt, dgDofContainer *solution);
static void registerBindings (binding *b); static void registerBindings (binding *b);
dgRungeKutta(); 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 #endif
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment