WENOStatistics.h

Go to the documentation of this file.

// -*- C++ -*-

// Copyright (C) California Institute of Technology
// Carlos Pantano

#ifndef AMROC_WENOSTATISTICS_H
#define AMROC_WENOSTATISTICS_H

#include <iostream>
#include <fstream>
#include <string>
#include <stdlib.h>
#include "DAGH.h"
#include "AMRBase.h"
#include "Interpolation.h"

#include "StatParser/StatParser.h"

//#define DEBUG_PRINT_STATISTICS

template <class VectorType, class InterpolationType, class OutputType, int dim>
class WENOStatistics : public AMRBase<VectorType,dim>, 
                       private StatParser<typename InterpolationType::point_type,
                                          typename VectorType::InternalDataType,dim> {
  typedef AMRBase<VectorType,dim> base;
  typedef typename VectorType::InternalDataType DataType;
  typedef typename InterpolationType::point_type point_type;
  typedef StatParser<point_type,DataType,dim> stparser;
  typedef typename stparser::symrec symrec;

  // bring members from templated base class into scope
  using StatParser<point_type,DataType,dim>::group_table;
  using StatParser<point_type,DataType,dim>::putsym;
  using StatParser<point_type,DataType,dim>::unlinksym;
  using StatParser<point_type,DataType,dim>::vectorize;
  using StatParser<point_type,DataType,dim>::buildcoords;
  using StatParser<point_type,DataType,dim>::clearcoords;

public:
  typedef typename base::vec_grid_fct_type vec_grid_fct_type;  
  typedef GridFunction<DataType,dim> grid_fct_type;

public:
  WENOStatistics(InterpolationType *int_ref, OutputType* out_ref) : base() {
    _Step = 1;
    _FileOutput = out_ref;
    _Interpolation = int_ref;
    _OutputFile = "stats";
    _DumpMeshes = 0;
  }

  void Setup(GridHierarchy *gh, const int& ghosts, int* shape, double* geom) {
    base::SetupData(gh, ghosts);

    int MaxLev;
    int nx, ny, nz;
    char s[36];

    MaxLev = MaxLevel(base::GH());

    nx = shape[0];
    putsym(TOK_INTEGER, STAT_NX)->vdata.ivalue = nx;
    // add level dependent point number
    sprintf(s, "%s(%d)", STAT_NX, 1);
    putsym(TOK_INTEGER, s)->vdata.ivalue = nx;
    for ( int i = 0; i < MaxLev-1 ; i++ ) {
        sprintf(s, "%s(%d)", STAT_NX, i+2);
        nx *= RefineFactor(base::GH(),i);
        putsym(TOK_INTEGER, s)->vdata.ivalue = nx;
    }
    putsym(TOK_INTEGER, STAT_NXF)->vdata.ivalue = nx;
    putsym(TOK_DOUBLE, STAT_XMIN)->vdata.dvalue = geom[0];
    putsym(TOK_DOUBLE, STAT_XMAX)->vdata.dvalue = geom[1];

    if ( dim > 1 ) {
        ny = shape[1];
        putsym(TOK_INTEGER, STAT_NY)->vdata.ivalue = ny;
        // add level dependent point number
        sprintf(s, "%s(%d)", STAT_NY, 1);
        putsym(TOK_INTEGER, s)->vdata.ivalue = ny;
        for ( int i = 0; i < MaxLev-1 ; i++ ) {
            sprintf(s, "%s(%d)", STAT_NY, i+2);
            ny *= RefineFactor(base::GH(),i);;
            putsym(TOK_INTEGER, s)->vdata.ivalue = ny;
        }
        putsym(TOK_INTEGER, STAT_NYF)->vdata.ivalue = ny;
        putsym(TOK_DOUBLE, STAT_YMIN)->vdata.dvalue = geom[2];
        putsym(TOK_DOUBLE, STAT_YMAX)->vdata.dvalue = geom[3];
    }

    if ( dim > 2 ) {
        nz = shape[2];
        putsym(TOK_INTEGER, STAT_NZ)->vdata.ivalue = nz;
        // add level dependent point number
        sprintf(s, "%s(%d)", STAT_NZ, 1);
        putsym(TOK_INTEGER, s)->vdata.ivalue = nz;
        for ( int i = 0; i < MaxLev-1 ; i++ ) {
            sprintf(s, "%s(%d)", STAT_NZ, i+2);
            nz *= RefineFactor(base::GH(),i);;
            putsym(TOK_INTEGER, s)->vdata.ivalue = nz;
        }
        putsym(TOK_INTEGER, STAT_NZF)->vdata.ivalue = nz;
        putsym(TOK_DOUBLE, STAT_ZMIN)->vdata.dvalue = geom[4];
        putsym(TOK_DOUBLE, STAT_ZMAX)->vdata.dvalue = geom[5];
    }

    /* simulation variables (dimension dependent) */
    // by convention I take positive iptr's as conservative
    // variables and negative iptr's as primitive variables.
    symrec *ptr;

    putsym(TOK_VARIABLE, "rho")->vdata.iptr = 1;
    putsym(TOK_VARIABLE, "r")->vdata.iptr = 1;
    (ptr = putsym(TOK_VARIABLE, "u"))->vdata.iptr = 2; ptr->src = STAT_DERIVED;
    if ( dim > 1 ) {
      (ptr = putsym(TOK_VARIABLE, "v"))->vdata.iptr = 3; ptr->src = STAT_DERIVED;
    }
    if ( dim > 2 ) {
      (ptr = putsym(TOK_VARIABLE, "w"))->vdata.iptr = 4; ptr->src = STAT_DERIVED;
    }
    putsym(TOK_VARIABLE, "E")->vdata.iptr = 5;
    for ( int i = 0; i < NVARS - 5 ; i++ ) {
      char name[12];
      sprintf(name, "Y%d", i+1);
      (ptr = putsym(TOK_VARIABLE, name))->vdata.iptr = 
        ((dim == 1 ? 7 : (dim == 2 ? 8 : 9))+i); ptr->src = STAT_DERIVED;
    }
    if ( dim == 3 ) {
      (ptr = putsym(TOK_VARIABLE, "sgske"))->vdata.iptr = (NVARS+5); ptr->src = STAT_DERIVED;
    }
    putsym(TOK_VARIABLE, "T")->vdata.iptr = NVARS+1;
    (ptr = putsym(TOK_VARIABLE, "p"))->vdata.iptr = 
      (dim == 1 ? 5 : (dim == 2 ? 6 : 7)); ptr->src = STAT_DERIVED;
    (ptr = putsym(TOK_VARIABLE, "gamma"))->vdata.iptr = 
      (dim == 1 ? 6 : (dim == 2 ? 7 : 8)); ptr->src = STAT_DERIVED;
    putsym(TOK_VARINARRAY, "q")->vdata.iptr = 0;
    (ptr = putsym(TOK_VARINARRAY, "qout"))->vdata.iptr = 0; ptr->src = STAT_DERIVED;

    // Process the rack definition file
    if ( _DefFile.length() ) { 
      stparser::parse(_DefFile.c_str());
      // if any group did not specify step sampling rate, then set it to default common value
      for ( typename stparser::grpvec::const_iterator g = group_table.begin() ; 
            g != group_table.end(); g++) 
        if ( (*g)->step == 0 ) (*g)->step = _Step;
      
      // write probe information: coordinates
      int me = MY_PROC; 
      if (me == VizServer) {
        if ( _DumpMeshes ) {
          int ig = 1;
          for ( typename stparser::grpvec::const_iterator g = group_table.begin() ; 
                g != group_table.end(); g++, ig++) {
            // for each probe
            int ip = 1;
            for ( typename stparser::provec::const_iterator p = (*g)->probe_table.begin() ;
                  p != (*g)->probe_table.end() ; p ++, ip++ )
              {
                if ( buildcoords(*p) != RETURN_OK ) return;
                
                int i, j;
                FILE * fp; 
                char name[128];
                
                sprintf(name, "group%d-probe%d.dat", ig, ip);
                fp = fopen(name, "w");
                if ( (*p)->type == TOK_THREAD ) {
                  for ( i = 0 ; i < (*p)->options.thread.na ; i++ ) {
                    fprintf(fp, "%12.8lf %12.8lf %12.8lf\n", (*p)->coords[i](0), 
                            (*p)->coords[i](1), 
                            (*p)->coords[i](2));
                  } 
                } else {
                  int na, nb;
                  if ( (*p)->type == TOK_SURFACE ) {
                    na = (*p)->options.surface.na;
                    nb = (*p)->options.surface.nb;
                  } else if ((*p)->type == TOK_PLANES ) {
                    na = (*p)->options.plane.na;
                    nb = (*p)->options.plane.nb;
                  }
                  fprintf(fp, "VARIABLES = \"X\", \"Y\", \"Z\"\n");
                  fprintf(fp, "ZONE I=%d, J=%d, F=POINT\n", na, nb);
                  for ( j = 0 ; j < nb ; j++ )
                    for ( i = 0 ; i < na ; i++ ) {
                      int l = i+j*na;
                      fprintf(fp, "%12.8lf %12.8lf %12.8lf\n", (*p)->coords[l](0),
                              (*p)->coords[l](1), (*p)->coords[l](2));
                    }
                }
                fclose(fp);
                
                clearcoords(*p);
              }
          }
        }
      }
    }
  }

  virtual void register_at(ControlDevice& Ctrl) { register_at(Ctrl, ""); }
  virtual void register_at(ControlDevice& Ctrl, const std::string& prefix) {
    LocCtrl = Ctrl.getSubDevice(prefix+"Statistics");
    
    RegisterAt(LocCtrl,"DefinitionFile",_DefFile);
    RegisterAt(LocCtrl,"OutputFile",_OutputFile);
    RegisterAt(LocCtrl,"Step",_Step);
    RegisterAt(LocCtrl,"WriteMeshes", _DumpMeshes);
  }

  ~WENOStatistics() {
    _FileOutput = NULL;
    _Interpolation = NULL;
  }

 public:

  void Evaluate(double * t, vec_grid_fct_type& U, grid_fct_type& Work)
    {
      int me ; 
      int Npoints;
      
      std::ofstream outfile;
      std::ostream* out;

      if ( group_table.size() == 0 ) return;

      //int TimeC = CurrentTime(base::GH(),0);
      int TimeC = StepsTaken(base::GH(),0);

      // check step sampling rate
      int gcheck = 0;
      for ( typename stparser::grpvec::const_iterator g = group_table.begin() ; 
            g != group_table.end(); g++) {
        if ( TimeC%((*g)->step) == 0 ) {
          gcheck = 1;
          break;
        }
      }

      //if ( TimeC%_Step != 0 ) return;
      if ( gcheck == 0 ) return;

      me = MY_PROC; 

      // open output file
      if (me == VizServer) {
        char fname[256];
        int time = CurrentTime(base::GH(), 0);
        std::sprintf(fname,"%s%d.dat",_OutputFile.c_str(),time);
        outfile.open(fname, std::ios::out );
        out = new std::ostream(outfile.rdbuf());
        // write time header
        *out << "# time " << t[0] << std::endl;
      }
      
      // for each group
      int ig = 1, iptr;
      for ( typename stparser::grpvec::const_iterator g = group_table.begin() ; 
            g != group_table.end(); g++, ig++ ) {

        if ( TimeC%((*g)->step) == 0 ) { 

        if (me == VizServer) 
          *out << "# group " << ig << std::endl;
        // for each probe
        int ip = 1;
        for ( typename stparser::provec::const_iterator p = (*g)->probe_table.begin() ; 
              p != (*g)->probe_table.end() ; p ++, ip++ )
          {
            if ( buildcoords(*p) != RETURN_OK ) return;
            
            // write probe type
            if (me == VizServer) 
              *out << "# probe " << token((*p)->type) << " " << ip << std::endl;
            
            // compute derived information
            switch ( (*p)->type ) {
            case TOK_THREAD:
              Npoints = (*p)->options.thread.na;
              break;
            case TOK_SURFACE:
              Npoints = (*p)->options.surface.na*(*p)->options.surface.nb;
              break;
            case TOK_PLANES:
              Npoints = (*p)->options.plane.na*(*p)->options.plane.nb;
              break;
            case TOK_VOLUME:
              /* not implemented yet */
              break;
            }
            
            symrec *ptr, *ptr_ref;
  
            /* first, make a copy of the stack */
            ptr = clonestack((*g)->sym_table);
  
            /* traverse the stack and evaluate */
            while ( ptr ) {
#ifdef DEBUG_PRINT_STATISTICS
                (comm_service::log() << token(ptr->type) << "\n").flush();
#endif
              // vectorize the data array that will be used for the result
              vectorize(ptr, Npoints);
              // the arguments of the operator do not need to be vectorized
              // for all processors since they are not used. Only VizServer
              // performs operations with them.
              
              switch( ptr->type ) {
              case TOK_OP_PLUS:
                  if (me == VizServer) {
                      vectorize(ptr->prev, Npoints);
                      vectorize(ptr->prev->prev, Npoints);
                      for ( int l = 0 ; l < Npoints ; l++ )
                          ptr->pdata[l] = ptr->prev->prev->pdata[l] + ptr->prev->pdata[l];
                  }
                  unlinksym(ptr->prev->prev);
                  unlinksym(ptr->prev);
                break;
              case TOK_OP_MINUS:
                if (me == VizServer) {
                    vectorize(ptr->prev, Npoints);
                    vectorize(ptr->prev->prev, Npoints);
                    for ( int l = 0 ; l < Npoints ; l++ )
                        ptr->pdata[l] = ptr->prev->prev->pdata[l] - ptr->prev->pdata[l];
                }
                unlinksym(ptr->prev->prev);
                unlinksym(ptr->prev);
                break;
              case TOK_OP_PROD:
                if (me == VizServer) {
                    vectorize(ptr->prev, Npoints);
                    vectorize(ptr->prev->prev, Npoints);
                    for ( int l = 0 ; l < Npoints ; l++ )
                        ptr->pdata[l] = ptr->prev->prev->pdata[l] * ptr->prev->pdata[l];
                }
                unlinksym(ptr->prev->prev);
                unlinksym(ptr->prev);
                break;
              case TOK_OP_DIV:
                if (me == VizServer) {
                    vectorize(ptr->prev, Npoints);
                    vectorize(ptr->prev->prev, Npoints);
                    for ( int l = 0 ; l < Npoints ; l++ )
                        ptr->pdata[l] = ptr->prev->prev->pdata[l] / ptr->prev->pdata[l];
                }
                unlinksym(ptr->prev->prev);
                unlinksym(ptr->prev);
                break;
              case TOK_OP_NEG:
                if (me == VizServer) {
                    vectorize(ptr->prev, Npoints);
                    for ( int l = 0 ; l < Npoints ; l++ )
                        ptr->pdata[l] = - ptr->prev->pdata[l];
                }
                unlinksym(ptr->prev);
                break;
              case TOK_OP_POWER:
                if (me == VizServer) {
                    vectorize(ptr->prev, Npoints);
                    vectorize(ptr->prev->prev, Npoints);
                    for ( int l = 0 ; l < Npoints ; l++ )
                        ptr->pdata[l] = pow(ptr->prev->prev->pdata[l], ptr->prev->pdata[l]);
                }
                unlinksym(ptr->prev->prev);
                unlinksym(ptr->prev);
                break;
              case TOK_POINTER:
                ptr_ref = ptr->ref;
                switch(ptr_ref->type) {
                case TOK_FUNCTION1:
                  if (me == VizServer) {
                      vectorize(ptr->prev, Npoints);
                      for ( int l = 0 ; l < Npoints ; l++ )
                          ptr->pdata[l] = ((func_t1)ptr_ref->vdata.ptr)(ptr->prev->pdata[l]);
                  }
                  unlinksym(ptr->prev);
                  break;
                case TOK_FUNCTION2:
                  if (me == VizServer) {
                      vectorize(ptr->prev, Npoints);
                      vectorize(ptr->prev->prev, Npoints);
                      for ( int l = 0 ; l < Npoints ; l++ )
                          ptr->pdata[l] = ((func_t2)ptr_ref->vdata.ptr)(ptr->prev->prev->pdata[l], 
                                                                        ptr->prev->pdata[l]);
                  }
                  unlinksym(ptr->prev->prev);
                  unlinksym(ptr->prev);
                  break;
                case TOK_VARINARRAY:
                case TOK_VARIABLE:
                  
                    // for variables in one of the AMR arrays, the pointer is located in the
                    // current token, not the defered one (ref). For variables in regular
                    // variables, the pointer is stored in the pointed variable.
                  if ( ptr_ref->type == TOK_VARINARRAY ) 
                    iptr = ptr->vdata.iptr;
                  else 
                    iptr = ptr_ref->vdata.iptr;
                  
                  // process the variables here
                  if ( ptr_ref->src == STAT_DERIVED ) {
                    // output variable type case
                    for (register int l=0; l<=FineLevel(base::GH()); l++) {
                      int Time = CurrentTime(base::GH(),l);
                      // if derived variable (call Transform)
                      _FileOutput->Transform(U, Work, Time, l, iptr, t[l]);
                    }
#ifdef DEBUG_PRINT_STATISTICS
                    (comm_service::log() << " convert output\n" ).flush();
#endif
                  } else {
                    int _dim = dim;
                    // conservative vector of state variable case
                    if ( _dim == 2 ) { 
                      for ( int Level = FineLevel(base::GH()) ; Level >= 0 ; Level --) {
                        int Time = CurrentTime(base::GH(),Level);
                        forall (U,Time,Level,c)
                          Coords ss(U(Time,Level,c).bbox().stepsize());
                          BBox bbi(U.boundingbbox(Time,Level,c));
                          BeginFastIndex2(U, U(Time,Level,c).bbox(), 
                                          U(Time,Level,c).data(), VectorType);
                          BeginFastIndex2(Work, Work(Time,Level,c).bbox(), 
                                          Work(Time,Level,c).data(), DataType);
                          for_2 (i, j, bbi, ss)
                            FastIndex2(Work,i,j) = FastIndex2(U,i,j)(iptr-1);
                          end_for
                          EndFastIndex2(U);
                          EndFastIndex2(Work);
                        end_forall
                        } 
                    } else if ( _dim == 3 ) {
                      for ( int Level = FineLevel(base::GH()) ; Level >= 0 ; Level --) {
                        int Time = CurrentTime(base::GH(),Level);
                        forall (U,Time,Level,c)
                          Coords ss(U(Time,Level,c).bbox().stepsize());
                          BBox bbi(U.boundingbbox(Time,Level,c));
                          BeginFastIndex3(U, U(Time,Level,c).bbox(),  
                                          U(Time,Level,c).data(), VectorType); 
                          BeginFastIndex3(Work, Work(Time,Level,c).bbox(), 
                                          Work(Time,Level,c).data(), DataType);
                          for_3 (i, j, k, bbi, ss)
                            FastIndex3(Work,i,j,k) = FastIndex3(U,i,j,k)(iptr-1);
                          end_for
                          EndFastIndex3(U);
                          EndFastIndex3(Work); 
                          end_forall
                      } 
                    }
#ifdef DEBUG_PRINT_STATISTICS
                    (comm_service::log() << " convert array\n" ).flush();
#endif
                  }
                  // note that we use data so the pointer token will be converted.
                  _Interpolation->PointsValues(Work,t[0],Npoints,(*p)->coords,ptr->pdata);
#ifdef DEBUG_PRINT_STATISTICS
                  (comm_service::log() << "before (" << iptr << ") ..." ).flush();
#endif
                  _Interpolation->ArrayCombine(VizServer,Npoints,ptr->pdata);
#ifdef DEBUG_PRINT_STATISTICS
                  (comm_service::log() << "after\n" ).flush();
#endif
#ifndef DAGH_NO_MPI
                  MPI_Barrier(comm_service::comm()); 
#endif
                  break;
                case TOK_COORDINATE:
                  if (me == VizServer) {
                    int iptr = (ptr_ref->vdata.iptr)-1;
                    if ( iptr >= dim ) {
                      staterror("invalid coordinate in this stack");
                      return;
                    }
                    for ( int l = 0 ; l < Npoints ; l++ )
                      ptr->pdata[l] = (*p)->coords[l](iptr);
                  }
                  break;
                  
                default:
                  staterror("there should not be any other pointer type left in this stack");
                  return;
                }
                break;
              case TOK_SEPARATOR:
                // output
                if (me == VizServer) {
                    int na, nb, avea, aveb;
                    vectorize(ptr->prev, Npoints);
                  // write probe data
                  switch((*p)->type) {
                  case TOK_THREAD:
                    if ( (*p)->options.thread.ave == STAT_AVERAGED ) {
                      double sum = 0.0;
                      for ( int ip = 0 ; ip < Npoints ; ip ++ )
                        sum += ptr->prev->pdata[ip];
                      *out << (sum/Npoints) << std::endl;
                    } else {
                      for ( int ip = 0 ; ip < Npoints ; ip ++ )
                        *out << ptr->prev->pdata[ip] << " ";
                      *out << std::endl;
                    }
                    break;
                  case TOK_SURFACE:
                      avea = (*p)->options.surface.avea;
                      aveb = (*p)->options.surface.aveb;
                      na = (*p)->options.surface.na;
                      nb = (*p)->options.surface.nb;
                  case TOK_PLANES:
                      if ( (*p)->type == TOK_PLANES ) {
                          avea = STAT_AVERAGED;
                          aveb = STAT_AVERAGED;
                          na = (*p)->options.plane.na;
                          nb = (*p)->options.plane.nb;
                      }
                      
                    if ( avea == STAT_AVERAGED && aveb == STAT_AVERAGED ) {
                      double sum = 0.0;
                      for ( int ip = 0 ; ip < Npoints ; ip ++ )
                        sum += ptr->prev->pdata[ip];
                      *out << (sum/Npoints) << std::endl;
                    } else if ( avea == STAT_RAW && aveb == STAT_AVERAGED ) {
                      for ( int i = 0 ; i < na ; i++ ) {
                        double sum=0.0;
                        for ( int j = 0 ; j < nb ; j++ ) {
                          sum += ptr->prev->pdata[i+j*na];
                          *out << (sum/nb) << " ";
                        }
                      }
                      *out << std::endl;
                    }  else if ( avea == STAT_AVERAGED && aveb == STAT_RAW ) {
                      for ( int j = 0 ; j < nb ; j++ ) {
                        double sum=0.0;
                        for ( int i = 0 ; i < na ; i++ ) {
                          sum += ptr->prev->pdata[i+j*na];
                          *out << (sum/na) << " ";
                        }
                      }
                      *out << std::endl;
                    } else {
                      for ( int ip = 0 ; ip < Npoints ; ip ++ )
                        *out << ptr->prev->pdata[ip] << " ";
                      *out << std::endl;
                    }
                    break;
                  case TOK_VOLUME:
                    
                    break;
                  }
                }
                unlinksym(ptr->prev);
                break;
              case TOK_FUNCTION1:
              case TOK_FUNCTION2:
              case TOK_VARINARRAY:
              case TOK_VARIABLE:
              case TOK_COORDINATE:
                staterror("there should not be any pointer type tokens"
                          " left in the stack");
                return;
              default:
                  if ( ptr->type != TOK_DOUBLE ) {
                      staterror("unhandled stack token");
                      return;
                  }
                break;
              }
              // need special code to handle the separator token, to flush the queue
              if ( ptr->type == TOK_SEPARATOR ) {
                  ptr_ref = ptr->next;
                  unlinksym(ptr);
                  ptr = ptr_ref;
              } else {
                  ptr->type = TOK_DOUBLE;
                  ptr = ptr->next;
              }
              
#ifdef DEBUG_PRINT_STATISTICS
              (comm_service::log() << " end of token switch\n" ).flush();
#endif
            }
          clearcoords(*p);  
          }
        }
      }
      // close output file
      if (me == VizServer) {
        outfile.close();
        delete out;
      }
    }
  
private:
  int _DumpMeshes;
  int _Step;
  ControlDevice LocCtrl;
  std::string _OutputFile;
  std::string _DefFile;
  
  OutputType* _FileOutput;
  InterpolationType* _Interpolation;
};

#endif