Skip to content
Snippets Groups Projects
Commit e885723d authored by Bruno Seny's avatar Bruno Seny
Browse files

First work on multirate Runge-Kutta

parent 9f9b5d32
No related branches found
No related tags found
No related merge requests found
Hide whitespace changes
Inline Side-by-side
Solver/TESTCASES/WavePulse.lua
+ 2
− 2
View file @ e885723d
......@@ -43,10 +43,10 @@ print'*** solve ***'
dt = 0.00125;
N = 1000;
for i=1,N do
norm = DG:RK44(dt)
norm = DG:multirateRK43(dt)
print('*** ITER ***',i,norm)
if (i % 10 == 0) then
DG:exportSolution(string.format("output/solution-%03d", i))
DG:exportSolution(string.format("output/solution-%04d", i))
end
end
......
Solver/dgAlgorithm.cpp
+ 134
− 0
View file @ e885723d
......@@ -319,6 +319,140 @@ void dgAlgorithm::rungeKutta (const dgConservationLaw &claw, // conservation l
}
}
void dgAlgorithm::multirateRungeKutta (const dgConservationLaw &claw, // conservation law
std::vector<dgGroupOfElements*> &eGroups, // group of elements
std::vector<dgGroupOfFaces*> &fGroups, // group of interfacs
std::vector<dgGroupOfFaces*> &bGroups, // group of boundaries
double h, // time-step
dgDofContainer &sol,
dgDofContainer &resd,
int orderRK) // order of RK integrator
{
// U_{n+1}=U_n+h*(SUM_i a_i*K_i)
// K_i=M^(-1)R(U_n+b_i*K_{i-1})
int nStages=10;
// classical RK44
// double A[4][4]={
// {0, 0, 0, 0},
// {1.0/2.0, 0, 0 ,0},
// {0, 1.0/2.0, 0, 0},
// {0, 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};
// 3/8 RK44
// double A[4][4]={
// {0, 0, 0, 0},
// {1.0/3.0, 0, 0 ,0},
// {-1.0/3.0, 1.0, 0, 0},
// {1, -1, 1, 0}
// };
// double b[4]={1.0/8.0, 3.0/8.0, 3.0/8.0, 1.0/8.0};
// double c[4]={0, 1.0/3.0, 2.0/3.0, 1};
// RK43 from Schlegel et al. JCAM 2009
// 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};
// Small step RK43
double A[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 b[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 c[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 A[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 b[10]={1.0/6.0, 0 , 0 , 0 , 1.0/3.0,1.0/3.0, 0 , 0 , 0 , 1.0/6.0};
// double c[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> K(sol);
// Current updated solution
// fullMatrix<double> Unp(sol);
fullMatrix<double> residuEl, KEl;
fullMatrix<double> iMassEl;
int nbFields = claw.nbFields();
dgDofContainer **K;
K=new dgDofContainer*[nStages];
for(int i=0;i<nStages;i++){
K[i]=new dgDofContainer(eGroups,claw);
K[i]->_data->scale(0.0);
}
dgDofContainer Unp (eGroups,claw);
dgDofContainer tmp (eGroups,claw);
Unp._data->scale(0.0);
Unp._data->axpy(*(sol._data));
for(int j=0; j<nStages;j++){
tmp._data->scale(0.0);
tmp._data->axpy(*(sol._data));
for(int k=0;k < eGroups.size();k++) {
for(int i=0;i<j;i++){
if(fabs(A[j][i])>1e-12){
tmp.getGroupProxy(k).add(K[i]->getGroupProxy(k),h*A[j][i]);
}
}
}
this->residual(claw,eGroups,fGroups,bGroups,tmp._dataProxys,resd._dataProxys);
for(int k=0;k < eGroups.size();k++) {
int nbNodes = eGroups[k]->getNbNodes();
for(int i=0;i<eGroups[k]->getNbElements();i++) {
residuEl.setAsProxy(*(resd._dataProxys[k]),i*nbFields,nbFields);
KEl.setAsProxy(*(K[j]->_dataProxys[k]),i*nbFields,nbFields);
iMassEl.setAsProxy(eGroups[k]->getInverseMassMatrix(),i*nbNodes,nbNodes);
iMassEl.mult(residuEl,KEl);
}
Unp.getGroupProxy(k).add(K[j]->getGroupProxy(k),h*b[j]);
}
}
for (int i=0;i<sol._dataSize;i++){
// printf("tempSol[%d] = %g\n",i,(*tempSol._data)(i));
// memcp
(*sol._data)(i)=(*Unp._data)(i);
}
for(int i=0;i<nStages;i++){
delete K[i];
}
delete []K;
}
void dgAlgorithm::residualBoundary ( //dofManager &dof, // the DOF manager (maybe useless here)
const dgConservationLaw &claw, // the conservation law
dgGroupOfFaces &group,
......
Solver/dgAlgorithm.h
+ 9
− 0
View file @ e885723d
......@@ -48,6 +48,15 @@ class dgAlgorithm {
dgDofContainer &residu,
dgDofContainer &sol,
int orderRK=4);
void multirateRungeKutta (
const dgConservationLaw &claw,
std::vector<dgGroupOfElements*> &eGroups, //group of elements
std::vector<dgGroupOfFaces*> &fGroups, // group of interfacs
std::vector<dgGroupOfFaces*> &bGroups, // group of boundaries
double h,
dgDofContainer &residu,
dgDofContainer &sol,
int orderRK=3);
static void multAddInverseMassMatrix ( /*dofManager &dof,*/
const dgGroupOfElements & group,
fullMatrix<double> &residu,
......
Solver/dgSystemOfEquations.cpp
+ 7
− 0
View file @ e885723d
......@@ -49,6 +49,7 @@ methodBinding *dgSystemOfEquations::methods[]={
new methodBindingTemplate<dgSystemOfEquations,void,std::string>("exportSolution",&dgSystemOfEquations::exportSolution),
new methodBindingTemplate<dgSystemOfEquations,void,std::string>("L2Projection",&dgSystemOfEquations::L2Projection),
new methodBindingTemplate<dgSystemOfEquations,double,double>("RK44",&dgSystemOfEquations::RK44),
new methodBindingTemplate<dgSystemOfEquations,double,double>("multirateRK43",&dgSystemOfEquations::multirateRK43),
0};
// do a L2 projection
......@@ -80,6 +81,12 @@ double dgSystemOfEquations::RK44(double dt){
return _solution->_data->norm();
}
double dgSystemOfEquations::multirateRK43(double dt){
_algo->multirateRungeKutta(*_claw, _elementGroups, _faceGroups, _boundaryGroups, dt, *_solution, *_rightHandSide);
return _solution->_data->norm();
}
void dgSystemOfEquations::exportSolution(std::string outputFile){
export_solution_as_is(outputFile);
}
......
Solver/dgSystemOfEquations.h
+ 4
− 1
View file @ e885723d
......@@ -13,8 +13,10 @@ private:
dgDofContainer (const dgDofContainer&);
public:
int _dataSize; // the full data size i.e. concerning all groups
fullVector<double> * _data; // the full data itself
std::vector<fullMatrix<double> *> _dataProxys; // proxys
fullVector<double> * _data; // the full data itself
inline int getDataSize(){return _dataSize;}
inline fullMatrix<double> &getGroupProxy(int gId){ return *(_dataProxys[gId]); }
dgDofContainer (std::vector<dgGroupOfElements*> &groups, const dgConservationLaw &claw);
~dgDofContainer ();
};
......@@ -51,6 +53,7 @@ public:
void setup (); // setup the groups and allocate
void exportSolution (std::string filename); // export the solution
double RK44 (double dt);
double multirateRK43 (double dt);
void L2Projection (std::string functionName); // assign the solution to a given function
static const char className[];
......
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment