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

@@ -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

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