From 6113d9ae5dd5cf7acbec455665a2bc2b1cd7d8f5 Mon Sep 17 00:00:00 2001
From: akva <akva@e10b68d5-8a6e-419e-a041-bce267b0401d>
Date: Wed, 28 Nov 2012 17:23:54 +0000
Subject: [PATCH] added: explicit Runge-Kutta based time stepping.

- any tableaux can be added, currently support euler, heun, RK3, RK4.

example of usage in the AdvectionDiffusion application. the integrand
is expected to assemble the mass matrix as the system matrix,
while the load vector should be the value of the operator evaluation.

git-svn-id: http://svn.sintef.no/trondheim/IFEM/trunk@2174 e10b68d5-8a6e-419e-a041-bce267b0401d
---
 Apps/Common/SIMExplicitRK.h | 129 ++++++++++++++++++++++++++++++++++++
 Apps/Common/TimeIntUtils.h  |  10 +++
 2 files changed, 139 insertions(+)
 create mode 100644 Apps/Common/SIMExplicitRK.h

diff --git a/Apps/Common/SIMExplicitRK.h b/Apps/Common/SIMExplicitRK.h
new file mode 100644
index 00000000..ac68682a
--- /dev/null
+++ b/Apps/Common/SIMExplicitRK.h
@@ -0,0 +1,129 @@
+//==============================================================================
+//!
+//! \file SIMExplicitRK.h
+//!
+//! \date Nov 28 2012
+//!
+//! \author Arne Morten Kvarving / SINTEF
+//!
+//! \brief Explicit Runge-Kutta based time stepping for SIM classes
+//!
+//==============================================================================
+
+#ifndef SIM_EXPLICIT_RK_H_
+#define SIM_EXPLICIT_RK_H_
+
+#include "TimeIntUtils.h"
+
+namespace TimeIntegration {
+
+  //! \brief Template can be instanced over any SIM implementing ISolver,
+  //         and which derive from SIMbase
+  template<class Solver>
+class SIMExplicitRK
+{
+public:
+  SIMExplicitRK(Solver& solv, Method type) : 
+    solver(solv)
+  {
+    if (type == EULER) {
+      RK.order = 1;
+      RK.b.push_back(1.0);
+      RK.c.push_back(0.0);
+      RK.A.redim(1,1);
+    }
+    if (type == HEUN) {
+      RK.order = 2;
+      RK.b.push_back(0.5);
+      RK.b.push_back(0.5);
+      RK.c.push_back(0.0);
+      RK.c.push_back(1.0);
+      RK.A.redim(2,2);
+      RK.A(2,1) = 1.0;
+    }
+    if (type == RK3) {
+      RK.order = 3;
+      RK.b.push_back(1.0/6.0);
+      RK.b.push_back(2.0/3.0);
+      RK.b.push_back(1.0/6.0);
+      RK.c.push_back(0.0);
+      RK.c.push_back(0.5);
+      RK.c.push_back(1.0);
+      RK.A.redim(3,3);
+      RK.A(2,1) =  0.5;
+      RK.A(3,1) = -1.0;
+      RK.A(3,2) =  2.0;
+    }
+    if (type == RK4) {
+      RK.order = 4;
+      RK.b.push_back(1.0/6.0);
+      RK.b.push_back(1.0/3.0);
+      RK.b.push_back(1.0/3.0);
+      RK.b.push_back(1.0/6.0);
+      RK.c.push_back(0.0);
+      RK.c.push_back(0.5);
+      RK.c.push_back(0.5);
+      RK.c.push_back(1.0);
+      RK.A.redim(4,4);
+      RK.A(2,1) = 0.5;
+      RK.A(3,2) = 0.5;
+      RK.A(4,3) = 1.0;
+    }
+  }
+
+  bool solveStep(TimeStep& tp)
+  {
+    std::cout <<"\n  step = "<< tp.step <<"  time = "<< tp.time.t << std::endl;
+
+    std::vector<Vector> stages(RK.b.size());
+
+    TimeDomain time(tp.time);
+    Vector dum;
+    for (size_t i=0;i<stages.size();++i) {
+      Vector tmp(solver.getSolution());
+      for (size_t j=0;j<i;++j)
+        tmp.add(stages[j], tp.time.dt*RK.A(i+1,j+1));
+      time.t = tp.time.t+tp.time.dt*(RK.c[i]-1.0);
+      solver.updateDirichlet(time.t, &dum);
+      solver.applyDirichlet(tmp);
+      solver.assembleSystem(time, Vectors(1, tmp));
+      // solve Mu = Au + f
+      solver.solveSystem(stages[i]);
+    }
+
+    // finally construct solution as weighted stages
+    solver.updateDirichlet(tp.time.t, &dum);
+    for (int i=0;i<RK.order;++i)
+      solver.getSolution().add(stages[i], tp.time.dt*RK.b[i]);
+
+    solver.applyDirichlet(solver.getSolution());
+
+    solver.printSolutionSummary(solver.getSolution(), 
+                                0, solver.getProblem()->getField1Name(1));
+
+    return true;
+  }
+
+  bool advanceStep(TimeStep& tp)
+  {
+    return solver.advanceStep(tp);
+  }
+
+  bool saveModel(char* fileName, int& nBlock)
+  {
+    return solver.saveModel(fileName, nBlock);
+  }
+
+  bool saveStep(const TimeStep& tp, int& nBlock)
+  {
+    return solver.saveStep(tp, nBlock);
+  }
+
+protected:
+  Solver& solver;
+  RKTableaux RK;
+};
+
+}
+
+#endif
diff --git a/Apps/Common/TimeIntUtils.h b/Apps/Common/TimeIntUtils.h
index c0ba05f7..b1e80616 100644
--- a/Apps/Common/TimeIntUtils.h
+++ b/Apps/Common/TimeIntUtils.h
@@ -13,6 +13,9 @@
 #ifndef _TIME_INT_UTILS_H
 #define _TIME_INT_UTILS_H
 
+#include "DenseMatrix.h"
+#include <vector>
+
 
 namespace TimeIntegration //!< Time integration scope
 {
@@ -28,6 +31,13 @@ namespace TimeIntegration //!< Time integration scope
     BDF2  = 6   //!< Second order backward differencing, implicit
   };
 
+  //! \brief Struct holding a Runge-Kutta tableaux
+  typedef struct {
+    int order;             //!< order of Scheme
+    DenseMatrix A;         //!< Coefficient matrix
+    std::vector<double> b; //!< Stage weights
+    std::vector<double> c; //!< Stage levels
+  } RKTableaux;
 
   //! \brief Returns the temporal order of the given method
   //! \param[in] method The method