diff --git a/Solver/TESTCASES/WavePulse.lua b/Solver/TESTCASES/WavePulse.lua
index 68bd5f0dc60430c8a556650b6542f96a7f44802f..5c679629b912b28c95a275d2e820c6b2dbc5980f 100644
--- a/Solver/TESTCASES/WavePulse.lua
+++ b/Solver/TESTCASES/WavePulse.lua
@@ -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
diff --git a/Solver/dgAlgorithm.cpp b/Solver/dgAlgorithm.cpp
index af63131bc0334574e59c71605440b37155cb4957..a376adddb593d06d722986741d61fee323d9698f 100644
--- a/Solver/dgAlgorithm.cpp
+++ b/Solver/dgAlgorithm.cpp
@@ -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,
diff --git a/Solver/dgAlgorithm.h b/Solver/dgAlgorithm.h
index a8ad3eb857a5aaf0eb6fd2f9d69e6709c5f21cef..823f9af3fbbfe4a1d2d361153c6e33795c3e4735 100644
--- a/Solver/dgAlgorithm.h
+++ b/Solver/dgAlgorithm.h
@@ -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,
diff --git a/Solver/dgSystemOfEquations.cpp b/Solver/dgSystemOfEquations.cpp
index 9035150742661bb60d94358b03570fdbe1dc349a..ae8357043d39a78b166aa010ba28a096a08f6d6a 100644
--- a/Solver/dgSystemOfEquations.cpp
+++ b/Solver/dgSystemOfEquations.cpp
@@ -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);
}
diff --git a/Solver/dgSystemOfEquations.h b/Solver/dgSystemOfEquations.h
index 7671f76d81383bb69f0da1959deef16d42e18084..e8cddcf7a547e23ab6e97a7c741bd90b2c08ad20 100644
--- a/Solver/dgSystemOfEquations.h
+++ b/Solver/dgSystemOfEquations.h
@@ -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[];