SIMExplicitRK.h

Go to the documentation of this file.

//==============================================================================
//==============================================================================
 
#ifndef SIM_EXPLICIT_RK_H_
#define SIM_EXPLICIT_RK_H_
 
#include "SIMenums.h"
#include "TimeIntUtils.h"
#include "TimeStep.h"
 
class DataExporter;
 
 
namespace TimeIntegration {
 
template<class Solver>
class SIMExplicitRK
{
public:
  SIMExplicitRK(Solver& solv, Method type, bool standalone = true) :
    solver(solv), alone(standalone)
  {
    if (type == EULER) {
      RK.order = 1;
      RK.b.push_back(1.0);
      RK.c.push_back(0.0);
      RK.A.resize(1,1);
    }
    else 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.resize(2,2);
      RK.A(2,1) = 1.0;
    }
    else 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.resize(3,3);
      RK.A(2,1) =  0.5;
      RK.A(3,1) = -1.0;
      RK.A(3,2) =  2.0;
    }
    else 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.resize(4,4);
      RK.A(2,1) = 0.5;
      RK.A(3,2) = 0.5;
      RK.A(4,3) = 1.0;
    }
    else
      RK.order = 0;
 
    Solver::msgLevel = 1; // prints primary solution summary only
  }
 
  const ProcessAdm& getProcessAdm() const { return solver.getProcessAdm(); }
 
  virtual bool solveStep(TimeStep& tp)
  {
    if (alone)
      solver.getProcessAdm().cout <<"\n  step = "<< tp.step <<"  time = "<< tp.time.t << std::endl;
 
    Vectors stages;
    return this->solveRK(stages, tp);
  }
 
  bool solveRK(Vectors& stages, const TimeStep& tp)
  {
    TimeDomain time(tp.time);
    Vector dum;
 
    stages.resize(RK.b.size());
 
    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);
      if (!solver.assembleSystem(time, Vectors(1, tmp),
                                 !linear || (tp.step == 1 && i == 0)))
        return false;
 
      // solve Mu = Au + f
      if (!solver.solveSystem(stages[i]))
        return false;
 
      if (linear)
        solver.setMode(SIM::RHS_ONLY);
    }
 
    // finally construct solution as weighted stages
    solver.updateDirichlet(tp.time.t, &dum);
    for (size_t i=0;i<RK.b.size();++i)
      solver.getSolution().add(stages[i], tp.time.dt*RK.b[i]);
    solver.applyDirichlet(solver.getSolution());
 
    if (alone)
      solver.printSolutionSummary(solver.getSolution(), 0,
                                  solver.getProblem()->getField1Name(1).c_str());
 
    return true;
  }
 
  bool advanceStep(TimeStep& tp)
  {
    return solver.advanceStep(tp);
  }
 
  bool saveModel(char* fileName, int& geoBlk, int& nBlock)
  {
    return solver.saveModel(fileName, geoBlk, nBlock);
  }
 
  bool saveStep(const TimeStep& tp, int& nBlock)
  {
    return solver.saveStep(tp, nBlock);
  }
 
  void registerFields(DataExporter& exporter)
  {
    solver.registerFields(exporter);
  }
 
  bool serialize(std::map<std::string,std::string>& data)
  {
    return solver.serialize(data);
  }
 
  bool deSerialize(const std::map<std::string,std::string>& data)
  {
    return solver.deSerialize(data);
  }
 
  void setLinear(bool enable) { linear = enable; }
 
protected:
  Solver& solver; 
  RKTableaux RK;  
  bool alone; 
  bool linear = false; 
};
 
}
 
#endif