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

Multirate implemented, must be debugged

parent 0386a196
No related branches found
No related tags found
No related merge requests found
Show whitespace changes
Inline Side-by-side
Solver/TESTCASES/rect.geo
+ 1
− 1
View file @ 90cfc361
......@@ -21,5 +21,5 @@ Physical Line("Top") = {6};
//Recombine Surface {9, 11};
Field[1] = MathEval;
Field[1].F = "0.01*1+0.01*1.7*y";
Field[1].F = "0.01*1+0.01*100*(y*y)";
Background Field = 1;
Solver/dgDofContainer.cpp
+ 7
− 0
View file @ 90cfc361
......@@ -160,6 +160,13 @@ void dgDofContainer::scale(double f)
_data->scale(f);
_ghostData->scale(f);
}
void dgDofContainer::scale(std::vector<dgGroupOfElements*>groupsVector, double f){
for(int i=0;i<groupsVector.size();i++){
dgGroupOfElements *g=groupsVector[i];
fullMatrix<double> &proxy=getGroupProxy(g);
proxy.scale(f);
}
}
void dgDofContainer::axpy(dgDofContainer &x, double a)
{
......
Solver/dgDofContainer.h
+ 1
− 0
View file @ 90cfc361
......@@ -26,6 +26,7 @@ private:
int _mshStep;
public:
void scale(double f);
void scale(std::vector<dgGroupOfElements*>groups, double f);
double norm();
void axpy(dgDofContainer &x, double a=1.);
void axpy(std::vector<dgGroupOfElements*>groups,dgDofContainer &x, double a=1.);
......
Solver/dgGroupOfElements.cpp
+ 10
− 8
View file @ 90cfc361
......@@ -920,8 +920,9 @@ void dgGroupCollection::buildGroupsOfInterfaces() {
// Split the groups of elements depending on their local time step
void dgGroupCollection::splitGroupsForMultirate(double dtRef,dgConservationLaw *claw, dgDofContainer *solution){
double dgGroupCollection::splitGroupsForMultirate(int maxLevels,dgConservationLaw *claw, dgDofContainer *solution){
Msg::Info("Splitting Groups for multirate time stepping");
maxLevels--;
int maxNumElems=getElementGroup(0)->getElement(0)->getGlobalNumber()+1;
std::vector<int>oldGroupIds;
oldGroupIds.resize(maxNumElems);
......@@ -1040,19 +1041,19 @@ void dgGroupCollection::splitGroupsForMultirate(double dtRef,dgConservationLaw *
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;
double 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++;
if(dtMaxExponent>maxLevels){
dtRef*=1.0/pow(2.0,(double)(dtMaxExponent-maxLevels));
dtMaxExponent=maxLevels;
}
_dtMaxExponent=dtMaxExponent;
std::vector<MElement *>currentNewGroup;
std::vector< dgGroupOfElements* >newGroups;// indexed by newGroupId
......@@ -1093,7 +1094,7 @@ void dgGroupCollection::splitGroupsForMultirate(double dtRef,dgConservationLaw *
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(localDt[el->getNum()]>=currentDt && (localDt[el->getNum()]<currentDt*2 || currentExponent==0)){
if(newGroupIds[el->getNum()]==-1){
mapNewGroups[iGroup].push_back(el);
newGroupIds[el->getNum()]=-2;
......@@ -1194,6 +1195,7 @@ void dgGroupCollection::splitGroupsForMultirate(double dtRef,dgConservationLaw *
Msg::Info("That makes a total of %d elements",count);
_elementGroups.clear();
_elementGroups=newGroups;
return dtRef;
}
......@@ -1243,7 +1245,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","solution",NULL);
cm->setArgNames("maxLevels","claw","solution",NULL);
cm = cb->addMethod("getNbElementGroups", &dgGroupCollection::getNbElementGroups);
cm->setDescription("return the number of dgGroupOfElements");
cm = cb->addMethod("getNbFaceGroups", &dgGroupCollection::getNbFaceGroups);
......
Solver/dgGroupOfElements.h
+ 1
− 1
View file @ 90cfc361
......@@ -245,7 +245,7 @@ class dgGroupCollection {
void buildGroupsOfElements (GModel *model,int dimension, int order);
void buildGroupsOfInterfaces ();
void splitGroupsForMultirate(double refDt,dgConservationLaw *claw, dgDofContainer *solution);
double splitGroupsForMultirate(int maxLevels,dgConservationLaw *claw, dgDofContainer *solution);
void find (MElement *elementToFind, int &iGroup, int &ithElementOfGroup);
......
Solver/dgRungeKutta.cpp
+ 71
− 29
View file @ 90cfc361
......@@ -4,6 +4,7 @@
#include "dgLimiter.h"
#include "dgResidual.h"
#include "dgGroupOfElements.h"
#include <stdio.h>
double dgRungeKutta::iterateEuler(const dgConservationLaw *claw, double dt, dgDofContainer *solution) {
double A[] = {0};
......@@ -423,7 +424,7 @@ dgRungeKuttaMultirate::dgRungeKuttaMultirate(dgGroupCollection* gc,dgConservatio
_residualVolume=new dgResidualVolume(*_law);
_residualInterface=new dgResidualInterface(*_law);
_K=new dgDofContainer*[10];
for(int i=0;i>10;i++){
for(int i=0;i<10;i++){
_K[i]=new dgDofContainer(gc,_law->getNbFields());
_K[i]->setAll(0.0);
}
......@@ -476,6 +477,7 @@ dgRungeKuttaMultirate::~dgRungeKuttaMultirate(){
}
void dgRungeKuttaMultirate::computeK(int iK,int exponent,bool isBuffer){
double localDt=_dt/pow(2.0,(double)exponent);
double _A[4][4]={
{0, 0, 0, 0},
{1.0/2.0, 0, 0 ,0},
......@@ -494,8 +496,6 @@ void dgRungeKuttaMultirate::computeK(int iK,int exponent,bool isBuffer){
{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]={
......@@ -510,29 +510,40 @@ void dgRungeKuttaMultirate::computeK(int iK,int exponent,bool isBuffer){
{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;
// Msg::Info("K%d, %d, %s",iK,exponent,isBuffer?"Buffer":"Bulk");
// 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;
std::vector<dgGroupOfElements *>&gVi=_innerBufferGroupsOfElements[exponent].first;
std::vector<dgGroupOfElements *>&gVo=_outerBufferGroupsOfElements[exponent].first;
_currentInput->scale(gVi,0.0);
_currentInput->scale(gVo,0.0);
_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]);
if(_AInner[iK][i]!=0.0)
_currentInput->axpy(gVi,*_K[i],_AInner[iK][i]*localDt);
if(_AOuter[iK][i]!=0.0)
_currentInput->axpy(gVo,*_K[i],_AOuter[iK][i]*localDt);
}
}
else{
std::vector<dgGroupOfElements *>gV=_bulkGroupsOfElements[exponent].first;
std::vector<dgGroupOfElements *> &gV=_bulkGroupsOfElements[exponent].first;
_currentInput->scale(gV,0.0);
_currentInput->axpy(gV,*_solution,1.0);
for(int i=0;i<iK;i++){
_currentInput->axpy(gV,*_K[i],_A[iK][i]);
if(_A[iK][i]!=0.0)
_currentInput->axpy(gV,*_K[i],_A[iK][i]*localDt);
}
}
computeResidual(exponent,isBuffer,_currentInput,_K[iK]);
/*
printf("%s%d\t",isBuffer?"B":"N",exponent);
for(int i=0;i<exponent;i++)
printf(" |");
printf(" K%d\n",iK);
*/
if(!isBuffer){
switch(iK){
case 0:
......@@ -555,28 +566,55 @@ void dgRungeKuttaMultirate::computeK(int iK,int exponent,bool isBuffer){
}
else{
if( (iK%5)<4 ){
computeK(iK%5,exponent+1,false);
computeK(iK%5,exponent,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){
void dgRungeKuttaMultirate::updateSolution(int exponent,bool isBuffer){
double localDt=_dt/pow(2.0,(double)exponent);
double b[4]={1.0/6.0, 1.0/3.0, 1.0/3.0, 1.0/6.0};
double c[4]={0, 1.0/2.0, 1.0/2.0, 1};
double bInner[10]={1.0/12.0, 1.0/6.0, 1.0/6.0,1.0/12.0,0,1.0/12.0, 1.0/6.0, 1.0/6.0,1.0/12.0,0};
double cInner[10]={0, 1.0/4.0, 1.0/4.0, 1.0/2.0, 1.0/2.0, 1.0/2.0, 3.0/4.0, 3.0/4.0, 1, 1 };
double bOuter[10]={1.0/6.0, 0 , 0 , 0 , 1.0/3.0,1.0/3.0, 0 , 0 , 0 , 1.0/6.0};
double cOuter[10]={0, 1.0/4.0, 1.0/4.0, 1.0/2.0, 1.0/2.0, 1.0/2.0, 3.0/4.0, 3.0/4.0, 1, 1 };
if(isBuffer){
std::vector<dgGroupOfElements *>*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);
std::vector<dgGroupOfElements *>&gVi=_innerBufferGroupsOfElements[exponent].first;
std::vector<dgGroupOfElements *>&gVo=_outerBufferGroupsOfElements[exponent].first;
for(int i=0;i<10;i++){
if(bInner[i]!=0.0)
_solution->axpy(gVi,*_K[i],bInner[i]*localDt);
if(bOuter[i]!=0.0)
_solution->axpy(gVi,*_K[i],bOuter[i]*localDt);
}
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;
else{
std::vector<dgGroupOfElements *>&gV=_bulkGroupsOfElements[exponent].first;
for(int i=0;i<4;i++){
if(b[i]!=0.0)
_solution->axpy(gV,*_K[i],b[i]*localDt);
}
}
}
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;
residual->scale(vei,0.0);
residual->scale(veo,0.0);
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());
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();
......@@ -592,13 +630,14 @@ void dgRungeKuttaMultirate::computeResidual(int exponent,bool isBuffer,dgDofCont
}
}
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];
std::vector<dgGroupOfElements *> &ve=_bulkGroupsOfElements[exponent].first;
residual->scale(ve,0.0);
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());
......@@ -615,6 +654,9 @@ void dgRungeKuttaMultirate::computeResidual(int exponent,bool isBuffer,dgDofCont
double dgRungeKuttaMultirate::iterate(double dt, dgDofContainer *solution){
_solution=solution;
_dt=dt;
_currentInput->scale(0.0);
_currentInput->axpy(*solution,1.0);
computeK(0,0,false);
computeK(1,0,false);
computeK(2,0,false);
......
Solver/dgRungeKutta.h
+ 1
− 0
View file @ 90cfc361
......@@ -34,6 +34,7 @@ class dgRungeKuttaMultirate: public dgRungeKutta{
private:
int _maxExponent;
int _minExponent;
double _dt;
dgConservationLaw *_law;
dgDofContainer *_solution;
dgDofContainer **_K;
......
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment