SIMImplicitLMM.h

Go to the documentation of this file.

//==============================================================================
//==============================================================================
 
#ifndef SIM_IMPLICIT_LMM_H_
#define SIM_IMPLICIT_LMM_H_
 
#include "NonLinSIM.h"
#include "SystemMatrix.h"
#include "SIMenums.h"
#include "TimeIntUtils.h"
#include "TimeStep.h"
 
class DataExporter;
 
 
namespace TimeIntegration {
 
  template<class Solver>
class SIMImplicitLMM
{
public:
  SIMImplicitLMM(Solver& solv, Method type, bool = false,
                 const std::string& solField = "") :
    solver(solv), nSim(solver,loads), fieldName(solField)
  {
    if (type == AM2)
      order = 2;
    else if (type == AM3)
      order = 3;
    else if (type == AM4)
      order = 4;
    else
      order = 1;
 
    loads.resize(order, nullptr);
 
    Solver::msgLevel = 1; // prints primary solution summary only
  }
 
  ~SIMImplicitLMM()
  {
    for (SystemVector* v : loads)
      delete v;
  }
 
  const ProcessAdm& getProcessAdm() const { return solver.getProcessAdm(); }
 
  bool solveStep(TimeStep& tp)
  {
    const std::vector<std::vector<double>> AM_coefs =
      {{1.0},
       {0.5, 0.5},
       {-1.0/12.0, 2.0/3.0, 5.0/12.0},
       {1.0/24.0, -5.0/24, 19.0/24.0, 9.0/24.0},
       {-19.0/720.0, 106.0/720.0, -264.0/720.0, 646.0/720.0, 251.0/720.0}};
 
    const int c_order = hasICs ? order-1 : std::min(order-1, tp.step-1);
    nSim.setCoefs(AM_coefs[c_order]);
 
    nSim.initSol(2);
    nSim.theSolutions()[0] = solver.getSolution();
    nSim.theSolutions()[1] = solver.getSolution();
    solver.setTimeScale(AM_coefs[c_order].back()*tp.time.dt);
    if (nSim.solveStep(tp, SIM::DYNAMIC) != SIM::CONVERGED)
      return false;
 
    solver.getSolution() = nSim.getSolution(0);
 
    solver.setMode(SIM::RHS_ONLY);
    solver.setTimeScale(1.0);
    if (!solver.assembleSystem(tp.time, Vectors(1, solver.getSolution()), false))
      return false;
 
    loads[0] = solver.getRHSvector(0, true);
 
    return true;
  }
 
  bool advanceStep(TimeStep& tp)
  {
    // Evaluate fluxes for initial conditions.
    if (tp.step == 1) {
      hasICs = true;
      for (int j = 2; j <= order && hasICs; ++j)
        if (std::string fName = fieldName + std::to_string(j);
            solver.hasIC(fName)) {
          TimeDomain time(tp.time);
          time.t = tp.time.t - j*tp.time.dt;
          if (!solver.assembleSystem(time, Vectors(1, solver.getSolution(j-1))))
            return false;
 
          loads[j-2] = solver.getRHSvector(0, true);
        }
        else
          hasICs = false;
    }
 
    delete loads.back();
    for (int j = order-2; j >= 0; --j)
      loads[j+1] = loads[j];
 
    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);
  }
 
protected:
  class LMMNonLinSIM : public NonLinSIM
  {
  public:
    LMMNonLinSIM(SIMbase& sim, std::vector<SystemVector*>& load)
      : NonLinSIM(sim), loads(load) {}
 
    void setCoefs(const std::vector<double>& coef) { coefs = coef; }
 
  protected:
    bool assembleSystem(const TimeDomain& time, const Vectors& pSol,
                        bool newLHSmatrix = true, bool poorConvg = false)
    {
      if (!model.assembleSystem(time, pSol, newLHSmatrix, poorConvg))
        return false;
 
      for (size_t i = 0; i < coefs.size()-1; ++i)
        model.addToRHSvector(0, *loads[coefs.size()-i-1], coefs[i]*time.dt);
 
      return true;
    }
 
  private:
    std::vector<SystemVector*>& loads; 
    std::vector<double>         coefs; 
  };
 
  Solver& solver; 
  LMMNonLinSIM nSim; 
  std::vector<SystemVector*> loads; 
  int order; 
  const std::string fieldName; 
  bool hasICs = false; 
};
 
}
 
#endif