From 36f25f99b394bfa0e9afab9e390077646726a3a3 Mon Sep 17 00:00:00 2001
From: Richard Comblen <richard.comblen@uclouvain.be>
Date: Wed, 17 Feb 2010 08:09:35 +0000
Subject: [PATCH] Added many explicit Runge-Kutta
---
Solver/dgLimiter.cpp | 2 +-
Solver/dgLimiter.h | 2 +-
Solver/dgRungeKutta.cpp | 153 ++++++++++++++++++++++++++++++++++++++--
Solver/dgRungeKutta.h | 7 +-
4 files changed, 157 insertions(+), 7 deletions(-)
diff --git a/Solver/dgLimiter.cpp b/Solver/dgLimiter.cpp
index 518674751d..45c7b55e99 100644
--- a/Solver/dgLimiter.cpp
+++ b/Solver/dgLimiter.cpp
@@ -132,7 +132,7 @@ void dgLimiter::registerBindings(binding *b) {
cm->setDescription("apply the limiter on the solution");
}
-void dgSlopeLimiter::registerBindings(binding *b) {
+void dgSlopeLimiterRegisterBindings(binding *b) {
classBinding *cb = b->addClass<dgSlopeLimiter>("dgSlopeLimiter");
cb->setDescription("The usual DG slope limiter: keep the mean, reduces uniformly the slope until it does not overshoot the neighbors' mean");
methodBinding *cm;
diff --git a/Solver/dgLimiter.h b/Solver/dgLimiter.h
index 5e4072beae..d8630beb45 100644
--- a/Solver/dgLimiter.h
+++ b/Solver/dgLimiter.h
@@ -21,7 +21,7 @@ class dgSlopeLimiter : public dgLimiter{
public :
dgSlopeLimiter (dgConservationLaw *claw) : dgLimiter (claw) {}
virtual int apply ( dgDofContainer *solution);
- static void registerBindings(binding *b);
};
+void dgSlopeLimiterRegisterBindings(binding *b);
#endif
diff --git a/Solver/dgRungeKutta.cpp b/Solver/dgRungeKutta.cpp
index ba3572a990..3480c5f5cc 100644
--- a/Solver/dgRungeKutta.cpp
+++ b/Solver/dgRungeKutta.cpp
@@ -11,18 +11,87 @@ double dgRungeKutta::iterateEuler(const dgConservationLaw *claw, double dt, dgDo
return diagonalRK(claw, dt, solution, 1, A, b);
}
+double dgRungeKutta::iterate22(const dgConservationLaw *claw, double dt, dgDofContainer *solution) {
+ double A[] = {0, 1};
+ double b[] = {1./2, 1./2};
+ return diagonalRK(claw, dt, solution, 2, A, b);
+}
+double dgRungeKutta::iterate33(const dgConservationLaw *claw, double dt, dgDofContainer *solution) {
+ double A[] = {0, 1./3, 2./3};
+ double b[] = {1./4,0, 3./4};
+ return diagonalRK(claw, dt, solution, 3, A, b);
+}
double dgRungeKutta::iterate44(const dgConservationLaw *claw, double dt, dgDofContainer *solution) {
double A[] = {0, 0.5, 0.5, 1};
double b[] = {1./6, 1./3, 1./3, 1./6};
return diagonalRK(claw, dt, solution, 4, A, b);
}
-double dgRungeKutta::iterate22(const dgConservationLaw *claw, double dt, dgDofContainer *solution) {
- double A[] = {0, 1};
- double b[] = {1./2, 1./2};
- return diagonalRK(claw, dt, solution, 2, A, b);
+
+double dgRungeKutta::iterate44ThreeEight(const dgConservationLaw *claw, double dt, dgDofContainer *solution) {
+ // 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}
+ };
+ fullMatrix<double> Afm(4,4);
+ for(int i=0;i<4;i++){
+ for(int j=0;j<4;j++){
+ Afm(i,j)=A[i][j];
+ }
+ }
+ 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};
+ return nonDiagonalRK(claw, dt, solution, 4, Afm, b, c);
}
+double dgRungeKutta::iterateRKFehlberg45(const dgConservationLaw *claw, double dt, dgDofContainer *solution){
+ double A[6][6]={
+ {0,0,0,0,0,0},
+ {1.0/4,0,0,0,0,0},
+ {3.0/32,9.0/32,0,0,0,0},
+ {1932.0/2197,-7200.0/2197,7296.0/2197,0,0,0},
+ {439.0/216,-8,3680.0/513,-845.0/4104,0,0},
+ {8.0/27,2,-3544.0/2565,-1859.0/4104,-11.0/40,0},
+ };
+ double c[6];
+ fullMatrix<double> Afm(6,6);
+ for(int i=0;i<6;i++){
+ c[i]=0;
+ for(int j=0;j<6;j++){
+ Afm(i,j)=A[i][j];
+ c[i]+=A[i][j];
+ }
+ }
+ double b[6]={16.0/235,6656.0/12825,28561.0/56430,-9.0/50,2.0/55};
+ return nonDiagonalRK(claw, dt, solution, 6, Afm, b, c);
+}
+
+
+double dgRungeKutta::iterate43SchlegelJCAM2009(const dgConservationLaw *claw, double dt, dgDofContainer *solution) {
+ // 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}
+ };
+ fullMatrix<double> Afm(4,4);
+ for(int i=0;i<4;i++){
+ for(int j=0;j<4;j++){
+ Afm(i,j)=A[i][j];
+ }
+ }
+ 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};
+ return nonDiagonalRK(claw, dt, solution, 4, Afm, b, c);
+}
+
+
+
+
double dgRungeKutta::diagonalRK (const dgConservationLaw *claw,
double dt, // time-step
dgDofContainer *sol,
@@ -70,6 +139,70 @@ double dgRungeKutta::diagonalRK (const dgConservationLaw *claw,
return sol->norm();
}
+
+
+double dgRungeKutta::nonDiagonalRK(const dgConservationLaw *claw,
+ double dt,
+ dgDofContainer *solution,
+ int nStages,
+ fullMatrix<double>&A, // c | A
+ double *b, // Standard Butcher tableau notation :___|__
+ double *c // |b^T
+ )
+{
+ 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) {}
@@ -87,9 +220,21 @@ void dgRungeKutta::registerBindings(binding *b) {
cm = cb->addMethod("iterate22",&dgRungeKutta::iterate22);
cm->setArgNames("law","dt","solution",NULL);
cm->setDescription("update solution by doing Heun second order runge-kutta step of time step dt for the conservation law");
+ cm = cb->addMethod("iterate33",&dgRungeKutta::iterate33);
+ cm->setArgNames("law","dt","solution",NULL);
+ cm->setDescription("update solution by doing classical third order runge-kutta step of time step dt for the conservation law");
cm = cb->addMethod("iterate44",&dgRungeKutta::iterate44);
cm->setArgNames("law","dt","solution",NULL);
cm->setDescription("update solution by doing classical fourth order runge-kutta step of time step dt for the conservation law");
+ cm = cb->addMethod("iterate44ThreeEight",&dgRungeKutta::iterate44ThreeEight);
+ cm->setArgNames("law","dt","solution",NULL);
+ cm->setDescription("update solution by doing a fourth order runge-kutta (3/8 rule, not the classical one, see Hairer & Wanner book) step of time step dt for the conservation law");
+ cm = cb->addMethod("iterate43SchlegelJCAM2009",&dgRungeKutta::iterate43SchlegelJCAM2009);
+ cm->setArgNames("law","dt","solution",NULL);
+ cm->setDescription("update solution by doing a four stage third order runge-kutta (from Schlegel et al. Journal of Computational and Applied Mathematics, 2009) step of time step dt for the conservation law");
+ 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("setLimiter",&dgRungeKutta::setLimiter);
cm->setArgNames("limiter",NULL);
cm->setDescription("if a limiter is set, it is applied after each RK stage");
diff --git a/Solver/dgRungeKutta.h b/Solver/dgRungeKutta.h
index c4a2495dc9..09a19cd2a9 100644
--- a/Solver/dgRungeKutta.h
+++ b/Solver/dgRungeKutta.h
@@ -4,15 +4,20 @@ class dgConservationLaw;
class dgDofContainer;
class dgLimiter;
class binding;
+#include "fullMatrix.h"
class dgRungeKutta {
double diagonalRK(const dgConservationLaw *claw, double dt, dgDofContainer *solution, int nStages, double *A, double *b); // c == A
- //static double nonDiagonalRK();
+ double nonDiagonalRK(const dgConservationLaw *claw, double dt, dgDofContainer *solution, int nStages, fullMatrix<double>&A, double *b, double *c);
dgLimiter *_limiter;
public:
void setLimiter(dgLimiter *limiter) { _limiter = limiter; }
double iterateEuler(const dgConservationLaw *claw, double dt, dgDofContainer *solution);
double iterate22(const dgConservationLaw *claw, double dt, dgDofContainer *solution);
+ double iterate33(const dgConservationLaw *claw, double dt, dgDofContainer *solution);
double iterate44(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 iterateRKFehlberg45(const dgConservationLaw *claw, double dt, dgDofContainer *solution);
static void registerBindings (binding *b);
dgRungeKutta();
};
--
GitLab