SIMExplicitLMM.h

Go to the documentation of this file.

//==============================================================================
//==============================================================================
 
#ifndef SIM_EXPLICIT_LMM_H_
#define SIM_EXPLICIT_LMM_H_
 
#include "SystemMatrix.h"
#include "SIMenums.h"
#include "TimeIntUtils.h"
#include "TimeStep.h"
 
#include <memory>
#include <sstream>
#include <utility>
 
class DataExporter;
 
 
namespace TimeIntegration {
 
template<class Solver>
class SIMExplicitLMM
{
public:
  SIMExplicitLMM(Solver& solv, Method type, bool standalone = true,
                 const std::string& solField = "") :
    solver(solv), alone(standalone), fieldName(solField)
  {
    if (type == AB2)
      order = 2;
    else if (type == AB3)
      order = 3;
    else if (type == AB4)
      order = 4;
    else if (type == AB5)
      order = 5;
    else
      order = 1;
 
    loads.resize(order);
 
    Solver::msgLevel = 1; // prints primary solution summary only
  }
 
  const ProcessAdm& getProcessAdm() const { return solver.getProcessAdm(); }
 
  bool solveStep(TimeStep& tp)
  {
    if (alone)
      solver.getProcessAdm().cout <<"\n  step = "<< tp.step <<"  time = "<< tp.time.t << std::endl;
 
    // Initialize time-dependent Dirichlet boundary conditions at start of step.
    if (!solver.initDirichlet(tp.time.t))
      return false;
 
    TimeDomain time(tp.time);
    time.t = tp.time.t - tp.time.dt;
 
    // Set up constraints at stage time for assembly.
    // Note: we don't apply absolute BCs before assembly like RK does;
    // the assembly process handles the constraint formulation.
    Vector zeroRef(solver.getSolution(1).size());
    if (!solver.updateDirichlet(time.t, &zeroRef, true))
      return false;
 
    if (!solver.assembleSystem(time, Vectors(1, solver.getSolution(1)),
                               !linear || (tp.step == 1)))
      return false;
 
    loads[0].reset(solver.getRHSvector(0, true));
 
    const std::vector<std::vector<double>> AB_coefs = {
       {1.0},
       {-0.5, 1.5},
       {5.0/12.0, -16.0/12.0, 23.0/12.0},
       {-9.0/24.0, 37.0/24, -59.0/24.0, 55.0/24.0},
       {251.0/720.0, -1274.0/720.0, 2616.0/720.0, -2774.0/720.0, 1901.0/720.0}
    };
 
    const int c_order = hasICs ? order-1 : std::min(order-1, tp.step-1);
    const std::vector<double>& AB_coef = AB_coefs[c_order];
 
    SystemVector* rhs = solver.getRHSvector(0, false);
    rhs->mult(AB_coef.back() * tp.time.dt);
 
    for (size_t j = 0; j < AB_coef.size()-1; ++j)
      solver.addToRHSvector(0, *loads[c_order-j], AB_coef[j]*tp.time.dt);
 
    if (!solver.solveSystem(solver.getSolution()))
      return false;
 
    if (linear)
      solver.setMode(SIM::RHS_ONLY);
 
    solver.getSolution() += solver.getSolution(1);
 
    // Enforce absolute boundary values at the new time level.
    Vector dum;
    if (!solver.updateDirichlet(tp.time.t, &dum))
      return false;
    if (!solver.applyDirichlet(solver.getSolution()))
      return false;
 
    if (alone)
      solver.printSolutionSummary(solver.getSolution(), 0,
                                  solver.getProblem()->getField1Name(1).c_str());
 
    return true;
  }
 
  bool advanceStep(TimeStep& tp)
  {
    // Evaluate fluxes for initial conditions.
    if (tp.step == 1) {
      hasICs = true;
      for (int j = 2; j <= order; ++j) {
        std::stringstream str;
        str << fieldName << j;
        if (solver.hasIC(str.str())) {
          TimeDomain time(tp.time);
          time.t = tp.time.t - j*tp.time.dt;
 
          Vector zeroRef(solver.getSolution(j-1).size());
          if (!solver.updateDirichlet(time.t, &zeroRef, true))
            return false;
 
          if (!solver.assembleSystem(time, Vectors(1, solver.getSolution(j-1))))
            return false;
 
          loads[j-2].reset(solver.getRHSvector(0, true));
        } else {
          hasICs = false;
          break;
        }
      }
    }
 
    loads.back().reset();
    for (int j = order-2; j >= 0; --j)
      loads[j+1] = std::move(loads[j]);
    loads[0].reset();
 
    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; 
  std::vector<std::unique_ptr<SystemVector>> loads; 
  int order; 
  bool alone; 
  const std::string fieldName; 
  bool hasICs = false; 
  bool linear = false; 
};
 
}
 
#endif