diagnostic.h

Go to the documentation of this file.

#ifndef DIAGNOSTIC_H
#define DIAGNOSTIC_H
#include "normalization.h"
#include "dimensionality.h"
#include "task.h"
#include "parameter.h"
#include "simulation.h"
#include "numeric.h"
#include <vector>
#include <fstream>
//---------------------------------------------------------------------------------------------------------------------
class Species;

#define ID_DIAGPORT "DiagPort"
#define ID_DIAGNOSTIC "Diagnostic"



#define GID_1DGRID "1d-grid"
#define GID_VALUE "scalar"
#define GID_1DPOS "1d-vec"
#define GID_1DPS "4d-vec"


#define DOUBLEVEC std::vector<double>

class DiagPort : public Task {
 protected:
  bool Series;
  Simulation* Sim;
 public:
  DiagPort () : Sim(NULL) {};
  virtual ~DiagPort () {};
  std::string GetClassName () const {return ID_DIAGPORT;};
  
  void Init ();
  
  bool Execute () {return true;}; 
  
  virtual void BeginStep (const std::string& strType = "", const std::string& strName = "") = 0;
  virtual void EndStep () = 0;
  virtual void Out (const double d) = 0;
  virtual void Out (const Position &pos) = 0;
  virtual void Out (const DOUBLEVEC& vec) = 0;
  virtual void Out (const Velocity& V, double dF = 1.0) = 0;
  virtual void Out (const grid<double>& g, double dF = 1.0) = 0;
  virtual void Out (const grid<Position>& g, double dF = 1.0) = 0;
  virtual void Out (const grid<Velocity>& g, double dF = 1.0) = 0;
  

 protected:
  PARAMETERMAP* MakeParamMap (PARAMETERMAP* pm);
};

class File : public DiagPort {
 protected:
  std::string FileName;
  std::ofstream outfile;
 public:
  File () {};
  virtual ~File () { if (outfile) outfile.close(); };
  
  virtual void Init ();
  
  virtual void BeginStep (const std::string& strType = "", const std::string& strName = "");
  virtual void EndStep ();
  
  virtual void Out (double d);
  virtual void Out (const Position &pos);
  
  virtual void Out (const DOUBLEVEC& vec);
  
  virtual void Out (const Velocity& V, double dF = 1.0);

  virtual void Out (const grid<double>& g, double dF = 1.0);
  
  virtual void Out (const grid<Position>& g, double dF = 1.0);
  
  virtual void Out (const grid<Velocity>& g, double dF = 1.0);
 protected:
  std::string ParseFileName (std::string strFile) const;
  // Serialization
 protected:
  virtual PARAMETERMAP* MakeParamMap (PARAMETERMAP* pm = NULL);
};

class Diagnostic : public Task {
 private:
  std::string strDiagPort;
 protected:
  Simulation* Sim;
  DiagPort* DP;
 public:
  Diagnostic () : DP(NULL), Sim(NULL) {}
  virtual ~Diagnostic () {}
  virtual std::string GetClassName () const { return ID_DIAGNOSTIC; }
  
  virtual void Init ();
 protected:
  virtual PARAMETERMAP* MakeParamMap (PARAMETERMAP* pm = NULL);
};

class Trajectory : public Diagnostic {
 protected:
  Species* pS;
  int ParticleNo;
 public:
  Trajectory () : pS(NULL) {};
  virtual ~Trajectory () {};
  
  virtual void Init ();
    
  virtual bool Execute ();
 protected:
  virtual PARAMETERMAP* MakeParamMap (PARAMETERMAP* pm = NULL);
};

// Diagnostic to write out the complete phase space of a Species
class PhaseSpace : public Diagnostic {
 protected:
  Species* pS;
 public:
  PhaseSpace () : pS(NULL) {};
  virtual ~PhaseSpace () {};
  
  virtual void Init ();
  
  virtual bool Execute ();
}; 

template<class TYPE>
class DiagValue : public Diagnostic {
 private:
  int Mean, F;
  TYPE Val;
  TYPE* pObj;
 protected:
  int LastOne;
 private:
  DiagValue();
 public:       
  DiagValue (TYPE* pobj);
  virtual ~DiagValue ();
  
  virtual bool Execute ();
 protected:
  virtual PARAMETERMAP* MakeParamMap (PARAMETERMAP* pm = NULL);
};

template<class TYPE>
class DiagGrid : public Diagnostic {
 private:
  grid<TYPE> G;
  grid<TYPE>* pG;
 protected:
  int LastOne;
  int Mean;
  int F;
  double Scale;
 private:
  DiagGrid ();        
 public:
  DiagGrid (grid<TYPE>* apG);
  virtual ~DiagGrid ();
  virtual void Init ();
  virtual bool Execute ();
 protected:
  virtual PARAMETERMAP* MakeParamMap (PARAMETERMAP* pm = NULL);
};

#ifndef ONE_DIMENSIONAL

template<class TYPE>
class DiagGridLine : public Diagnostic {
 protected:
  int line,dir;
  grid<TYPE>* pG;
 private:
  DiagGridLine ();
 public:
  DiagGridLine (grid<TYPE>* apG);
  virtual ~DiagGridLine ();
  virtual void Init ();
  virtual bool Execute ();
 protected:
  virtual PARAMETERMAP* MakeParamMap (PARAMETERMAP* pm = NULL);
};
#endif
template<class TYPE>
DiagValue<TYPE>::DiagValue (TYPE* pobj) : pObj(pobj) {
  F = 0; 
  Val = 0; 
  LastOne = 0; 
}
template<class TYPE>
DiagValue<TYPE>::~DiagValue () {}
template<class TYPE>
bool DiagValue<TYPE>::Execute () {
  if (nEnd<0) return true;
  int nCycle = pProcess->GetCycle();
  // Accumulate average
  if (nCycle>LastOne+nStep-Mean) {
    F++;
    Val += (*pObj);
  }
  if (!DoNow()) return false;
  DP->BeginStep(GID_VALUE, strName);
  DP->Out( Val/double(F) );
  DP->EndStep();
  // reset values for next calculation of average
  F = 0;
  Val = 0;
  LastOne = nCycle;
  return false;
};

template<class TYPE>
PARAMETERMAP* DiagValue<TYPE>::MakeParamMap (PARAMETERMAP* pm) {
  pm = Diagnostic::MakeParamMap(pm);
  (*pm)["mean"] = WParameter(new ParameterValue<int>(&Mean, 1));
  return pm;
};
//--diagnostic.t-------------------------------------------------------------------------------------------------------
template<class TYPE>
DiagGrid<TYPE>::DiagGrid (grid<TYPE>* apG) : pG(apG) {}
template<class TYPE>
DiagGrid<TYPE>::~DiagGrid () {}
template<class TYPE>
void DiagGrid<TYPE>::Init () {
  Diagnostic::Init();
  GridPosition Nx = pG->Size();
  G.Resize(Nx);
  G = 0.0;
  F = 0;
  LastOne = nStart;
};
template<class TYPE>
bool DiagGrid<TYPE>::Execute () {
  if (nEnd < 0) return true;
  int nCycle = pProcess->GetCycle();
  if (nCycle>LastOne+nStep-Mean) {
    F++;
    // Accumulate values
    G += *pG;
  }
  if (!DoNow()) return false;
  DP->BeginStep(GID_1DGRID, strName);
  DP->Out(G, 1 / (F*Scale));
  DP->EndStep();
  // Reset averaging grid and count
  G = 0.0;
  F = 0;
  LastOne = nCycle;
  return false;
};
template<class TYPE>
//---------------------------------------------------------------------------------------------------------------------
PARAMETERMAP* DiagGrid<TYPE>::MakeParamMap(PARAMETERMAP* pm) {
  pm = Diagnostic::MakeParamMap(pm);
  (*pm)["mean"] = WParameter(new ParameterValue<int>(&Mean, 1));
  (*pm)["scale"] = WParameter(new ParameterValue<double>(&Scale, 1.0));
  return pm;
};
//---------------------------------------------------------------------------------------------------------------------
// DiagGridLine
#ifndef ONE_DIMENSIONAL
template<class TYPE>
DiagGridLine<TYPE>::DiagGridLine(grid<TYPE>* apG) : pG(apG) {}
template<class TYPE>
DiagGridLine<TYPE>::~DiagGridLine () {}
template<class TYPE>
void DiagGridLine<TYPE>::Init () {
  Diagnostic::Init();
  // error if the direction is not within the dimensional range
  if ((dir<0) || (dir>=Dimension)) {
    std::cerr << "Invalid direction " << dir << " in " << GetName() << std::endl;
    dir = 0;
  }
};
template<class TYPE>
bool DiagGridLine<TYPE>::Execute () {
  if (nEnd < 0) return true;
  int nCycle = pProcess->GetCycle();
  if (!DoNow()) return false;
  DP->BeginStep(GID_1DGRID, strName);
  // set all components to line
  GridPosition i =line;
  GridPosition NX = pG->Size();
  // and iterate through the component specified by direction
  for (i[dir]=0; i[dir]<=NX[dir]; i[dir]++) DP->Out((*pG)[i]);    
  DP->EndStep();
  return false;
}
template<class TYPE>
PARAMETERMAP* DiagGridLine<TYPE>::MakeParamMap (PARAMETERMAP* pm) {
  pm = Diagnostic::MakeParamMap(pm);
  (*pm)["line"] = WParameter(new ParameterValue<int>(&line, 1));
  (*pm)["direction"] = WParameter(new ParameterValue<int>(&dir, 0));
  return pm;
}
#endif
#endif // DIAGNOSTIC_H