SIMSolver.h

Go to the documentation of this file.

// $Id$
//==============================================================================
//==============================================================================
 
#ifndef _SIM_SOLVER_H_
#define _SIM_SOLVER_H_
 
#include "IFEM.h"
#include "SIMadmin.h"
#include "TimeStep.h"
#include "HDF5Restart.h"
#include "HDF5Writer.h"
#include "tinyxml2.h"
 
 
template<class T1> class SIMSolverStat : public SIMadmin
{
public:
  explicit SIMSolverStat(T1& s1, const char* head = nullptr) : SIMadmin(head), S1(s1)
  {
    exporter = nullptr;
    startExpLevel = 0;
  }
 
  virtual ~SIMSolverStat() { delete exporter; }
 
  virtual bool read(const char*) { return true; }
 
  void handleDataOutput(const std::string& hdf5file,
                        const ProcessAdm& modelAdm,
                        int saveInterval = 1)
  {
    if (IFEM::getOptions().discretization == ASM::Spectral && !hdf5file.empty())
      IFEM::cout <<"\n  ** HDF5 output is not available for spectral"
                 <<" discretization. Deactivating...\n"<< std::endl;
    else
    {
      exporter = new DataExporter(true,saveInterval,startExpLevel);
      exporter->registerWriter(new HDF5Writer(hdf5file,modelAdm));
      S1.registerFields(*exporter);
      IFEM::registerCallback(*exporter);
    }
  }
 
protected:
  void printHeading(const char* heading)
  {
    if (heading)
    {
      std::string myHeading(heading);
      size_t n = myHeading.find_last_of('\n');
      if (n+1 < myHeading.size()) n = myHeading.size()-n;
      IFEM::cout <<"\n\n"<< myHeading <<"\n";
      for (size_t i = 0; i < n && i < myHeading.size(); i++) IFEM::cout <<'=';
      IFEM::cout << std::endl;
    }
  }
 
public:
  virtual int solveProblem(char* infile, const char* heading = nullptr)
  {
    // Save FE model to VTF for visualization
    int geoBlk = 0, nBlock = 0;
    if (!S1.saveModel(infile,geoBlk,nBlock))
      return 1;
 
    this->printHeading(heading);
 
    // Solve the stationary problem
    TimeStep dummy;
    if (!S1.solveStep(dummy))
      return 2;
 
    // Save the results
    if (!S1.saveStep(dummy,nBlock))
      return 4;
 
    if (exporter && !exporter->dumpTimeLevel())
      return 5;
 
    return 0;
  }
 
protected:
  T1& S1; 
 
  DataExporter* exporter; 
  int      startExpLevel; 
};
 
 
template<class T1> class SIMSolver : public SIMSolverStat<T1>
{
public:
  explicit SIMSolver(T1& s1) : SIMSolverStat<T1>(s1,"Time integration driver")
  {
    saveDivergedSol = dumpLog = false;
    restartAdm = nullptr;
    restartProcAdm = nullptr;
    startRstLevel = 0;
  }
 
  virtual ~SIMSolver() { delete restartAdm; delete restartProcAdm; }
 
  virtual bool read(const char* file) { return this->SIMadmin::read(file); }
 
  const TimeStep& getTimePrm() const { return tp; }
 
  bool advanceStep() { return tp.increment() && this->S1.advanceStep(tp); }
 
  virtual int solveProblem(char* infile, const char* heading = nullptr)
  {
    // Save FE model to VTF and HDF5 for visualization.
    // Save the initial configuration also if multi-step simulation.
    int geoBlk = 0, nBlock = 0;
    if (!this->saveState(geoBlk,nBlock,true,infile,tp.multiSteps()))
      return 2;
 
    this->printHeading(heading);
 
    // Solve for each time step up to final time
    while (this->advanceStep())
      if (!this->S1.solveStep(tp))
        return saveDivergedSol && !this->S1.saveStep(tp,nBlock) ? 4 : 3;
      else if (!this->saveState(geoBlk,nBlock))
        return 4;
      else
        IFEM::pollControllerFifo();
 
    return 0;
  }
 
  void handleDataOutput(const std::string& hdf5file,
                        const ProcessAdm& modelAdm,
                        int saveInterval = 1,
                        int restartInterval = 0)
  {
    this->SIMSolverStat<T1>::handleDataOutput(hdf5file, modelAdm, saveInterval);
    if (restartInterval < 1) return; // No restart output
 
    const ProcessAdm* adm = nullptr;
    if (!modelAdm.dd.isPartitioned())
      adm = &modelAdm;
    else if (modelAdm.getProcId() == 0)
      adm = restartProcAdm = new ProcessAdm();
    else
      return; // Restart output on processor 0 only
 
    restartAdm = new HDF5Restart(hdf5file+"_restart", *adm, restartInterval,
                                 startRstLevel);
  }
 
  virtual bool serialize(HDF5Restart::SerializeData& data)
  {
    return tp.serialize(data) && this->S1.serialize(data);
  }
 
  virtual bool deSerialize(const HDF5Restart::SerializeData& data)
  {
    return tp.deSerialize(data) && this->S1.deSerialize(data);
  }
 
protected:
  virtual bool parse(char* keyw, std::istream& is) { return tp.parse(keyw,is); }
 
  virtual bool parse(const tinyxml2::XMLElement* elem)
  {
    if (strcasecmp(elem->Value(),"postprocessing"))
      return tp.parse(elem);
 
    const tinyxml2::XMLElement* child = elem->FirstChildElement();
    for (; child; child = child->NextSiblingElement())
      if (!strncasecmp(child->Value(),"savediverg",10))
        saveDivergedSol = true;
      else if (!strncasecmp(child->Value(),"dumplog",7))
        dumpLog = true;
 
    return true;
  }
 
  bool saveState(int& geoBlk, int& nBlock, bool newMesh = false,
                 char* infile = nullptr, bool saveRes = true)
  {
    if (newMesh && !this->S1.saveModel(infile,geoBlk,nBlock))
      return false; // saving new geometry to VTF failed
    else if (!saveRes)
      return true;  // no result saving
    else if (!this->S1.saveStep(tp,nBlock))
      return false; // saving results to VTF failed
 
    DataExporter* exportAdm = SIMSolverStat<T1>::exporter;
    if (exportAdm && !exportAdm->dumpTimeLevel(&tp,newMesh,dumpLog))
      return false; // saving results to HDF5 failed
 
    if (!restartAdm || !restartAdm->dumpStep(tp))
      return true;  // no restart output for this step
 
    if (dumpLog)
      IFEM::cout <<"  Serializing simulator state for restart (time level "
                 << restartAdm->getTimeLevel()+1 <<")"<< std::endl;
 
    HDF5Restart::SerializeData data;
    if (exportAdm)
      data["TimeLevel"] = std::to_string(exportAdm->getTimeLevel());
    return this->serialize(data) ? restartAdm->writeData(data) : true;
  }
 
public:
  int restart(const std::string& restartFile, int restartStep)
  {
    if (restartFile.empty()) return 0;
 
    ProcessAdm oneAdm;
    const ProcessAdm& adm = this->S1.getProcessAdm();
    HDF5Restart hdf(restartFile, adm.dd.isPartitioned() ? oneAdm : adm);
 
    HDF5Restart::SerializeData data;
    if ((restartStep = hdf.readData(data,restartStep-1)) >= 0)
    {
      IFEM::cout <<"\n === Restarting from a serialized state ==="
                 <<"\n     file = "<< restartFile
                 <<"\n     step = "<< restartStep+1;
      if (this->deSerialize(data))
      {
        // Record time levels in case we are saving results in the restart also
        startRstLevel = restartStep;
        HDF5Restart::SerializeData::const_iterator it = data.find("TimeLevel");
        if (it != data.end())
        {
          SIMSolverStat<T1>::startExpLevel = atoi(it->second.c_str());
          IFEM::cout <<" "<< SIMSolverStat<T1>::startExpLevel+1;
        }
        else
          SIMSolverStat<T1>::startExpLevel = startRstLevel;
      }
      else
        restartStep = -2;
      IFEM::cout << std::endl;
    }
 
    if (restartStep < 0)
      std::cerr <<" *** SIMSolver: Failed to read restart data."<< std::endl;
    return restartStep;
  }
 
private:
  bool saveDivergedSol; 
 
protected:
  bool dumpLog; 
  TimeStep tp;  
  HDF5Restart*    restartAdm; 
  ProcessAdm* restartProcAdm; 
  int          startRstLevel; 
};
 
#endif