SolverControl.h

Go to the documentation of this file.

// -*- C++ -*-

// Copyright (C) 2002 Ralf Deiterding
// Brandenburgische Universitaet Cottbus
//
// Copyright (C) 2003-2007 California Institute of Technology
// Ralf Deiterding, ralf@amroc.net

#ifndef AMROC_SOLVER_CONTROL_H
#define AMROC_SOLVER_CONTROL_H

#if ! defined __PGI
# if defined(LINUX_387_DOUBLE) || defined(LINUX_387_EXTENDED)
#  include <fpu_control.h>
# endif
#endif

#ifdef HAVE_CONFIG_H
#include <config.h>        /* check config for unistd.h header */
#endif

#if defined(HAVE_UNISTD_H)
#include <unistd.h>        /* for chdir() declaration */
#else
extern "C" {
extern int chdir(const char *);
}
#endif

#define MaxStepControls  10
#define FACTOR 1.0e5

#include "Solver.h"

#ifndef DAGH_NO_MPI 
#include <mpi.h>
#endif 

#include <iostream>
#include <string>
#include <cstdio>
#include <cstdlib>
#include <cfloat>

class StepControl : public controlable {
public:
  StepControl() : LastTime(0.0), Outputs(0), OutputEvery(0), CheckPointEvery(0) {}

  virtual ~StepControl() {}

  virtual void register_at(ControlDevice& Ctrl) {}
  virtual void register_at(ControlDevice& Ctrl, const std::string& prefix) {
    LocCtrl = Ctrl.getSubDevice(prefix+"Control");
    RegisterAt(LocCtrl,"LastTime",LastTime);
    RegisterAt(LocCtrl,"Outputs",Outputs);
    RegisterAt(LocCtrl,"OutputEvery",OutputEvery);
    RegisterAt(LocCtrl,"CheckEvery",CheckPointEvery);
  }
  
public: 
  double LastTime;
  int Outputs, OutputEvery, CheckPointEvery;
  ControlDevice LocCtrl;     
};


class SolverControl : public controlable {
public:
  typedef Solver solver_type;

  SolverControl(solver_type& solver) : InitOutput(1), Restart(0), 
     _Solver(solver), NStepControls(1), allow_solve(false) {}
  
  virtual ~SolverControl() {}
  
  virtual void register_at(ControlDevice& Ctrl) {
    register_at(Ctrl,"");
  }
  virtual void register_at(ControlDevice& Ctrl, const std::string& prefix) {
    LocCtrl = Ctrl.getSubDevice(prefix+"SolverControl");
    Step[0].register_at(LocCtrl,"");    
    for (int cs=0; cs<MaxStepControls; cs++) {
      char text[5];
      std::sprintf(text,"%d-",cs+1);
      Step[cs].register_at(LocCtrl,std::string(text));
    }
    RegisterAt(LocCtrl,"Controls",NStepControls); 
    RegisterAt(LocCtrl,"Restart",Restart);  
    RegisterAt(LocCtrl,"InitOutput",InitOutput);      
    Solver_().register_at(Ctrl,"");
  }

  virtual void init() {     
    Solver_().init(); 
  }

  virtual void init(int argc, char *argv[]) {     
#if ! defined __PGI
# if defined(LINUX_387_DOUBLE) || defined(LINUX_387_EXTENDED)
    fpu_control_t cw;
    _FPU_GETCW(cw);
    cw &=  ~(_FPU_EXTENDED | _FPU_DOUBLE | _FPU_SINGLE);
  // cw &=  ~(_FPU_EXTENDED | _FPU_DOUBLE | _FPU_SINGLE | _FPU_MASK_ZM );
#  ifdef LINUX_387_DOUBLE
    cw |= _FPU_DOUBLE;
#  else
    cw |= _FPU_EXTENDED;
#  endif
    _FPU_SETCW(cw);
# endif
#endif
   
#ifndef DAGH_NO_MPI
    std::cout << "Initializing MPI\n" << std::flush;
    MPI_Init( &argc, &argv );
#endif  
    if (argc>1 && std::strlen(argv[1])>0)
      chdir((const char *)argv[1]);

    init();

    std::string sname = "solver.in";
    std::ifstream infile;
    infile.open(sname.c_str(), std::ios::in);
    if (!infile) {
      std::cerr << sname << " not found. Aborting." << std::endl;
      std::exit(-1);
    } 
    infile.close();

    MainCtrl = ControlDevice(GetFileControlDevice("solver.in",""));
    register_at(MainCtrl);
    MainCtrl.update();   

    update();

    allow_solve = Solver_().setup(); 
  }

  virtual void update() { 
    Solver_().update(); 
  }
  
  virtual void finish() {
    Solver_().finish();
#ifndef DAGH_NO_MPI
    std::cout << "Finalizing MPI\n" << std::flush;
    MPI_Finalize();
#endif  
  } 

  virtual void solve() {
    if (!allow_solve) return;

    double t=0., dt=0.;
    int CStepControl = 0;

    if (Restart) 
      Solver_().Restart(t,dt);

    if (!Restart || (t==0. && dt==0.))
      Solver_().Initialize(t,dt);     
    
    while (t>=Step[CStepControl].LastTime && CStepControl<NStepControls)
      CStepControl++;  
    
    double tout = (CStepControl>0 ? Step[CStepControl-1].LastTime : 0.0);
    double dtout = Step[CStepControl].LastTime - 
      (CStepControl>0 ? Step[CStepControl-1].LastTime : 0.0);
    if (Step[CStepControl].Outputs > 0)
      dtout /= Step[CStepControl].Outputs;
    while (t+DBL_EPSILON*FACTOR > tout+dtout)
      tout += dtout;
    
    if (InitOutput && (Step[CStepControl].Outputs>0 ||
                       Step[CStepControl].OutputEvery>0))
      Solver_().Output();       

    while (Step[NStepControls-1].LastTime-tout > DBL_EPSILON*FACTOR) {   
      while (tout+dtout-t > DBL_EPSILON*FACTOR) {   
        if (Step[CStepControl].OutputEvery==0)
          dt = std::min(dt,tout+dtout-t);
        if (t+dt > Step[NStepControls-1].LastTime) 
          dt = Step[NStepControls-1].LastTime-t;
        Solver_().Advance(t,dt);
        if (Step[CStepControl].OutputEvery > 0) 
          if (Solver_().NSteps() % Step[CStepControl].OutputEvery == 0) 
            Solver_().Output();
        if (Step[CStepControl].CheckPointEvery > 0) 
          if (Solver_().NSteps() % Step[CStepControl].CheckPointEvery == 0) 
            Solver_().Checkpointing();
      }
      Solver_().Output(); 
      tout += dtout;
      if (Step[CStepControl].LastTime-tout < DBL_EPSILON*FACTOR) {
        CStepControl++;     
        if (CStepControl<NStepControls) {   
          dtout = Step[CStepControl].LastTime - Step[CStepControl-1].LastTime;
          if (Step[CStepControl].Outputs > 0)
            dtout /= Step[CStepControl].Outputs;
        }
      }
    }
        
    Solver_().Output(); 
    Solver_().Checkpointing();
  } 

  inline solver_type& Solver_() { return _Solver; }
  inline const solver_type& Solver_() const { return _Solver; }
  
protected:
  int InitOutput;
  int Restart;  
  solver_type& _Solver; 
  int NStepControls;    
  bool allow_solve;
  StepControl Step[MaxStepControls];
  ControlDevice MainCtrl, LocCtrl;
};


#endif

---