amroc/amr/Interfaces/SchemeExactSolution.h

Go to the documentation of this file.

// -*- C++ -*-

// Copyright (C) 2013 Ralf Deiterding
// German Aerospace Center (DLR)

#ifndef SCHEME_EXACT_SOLUTION_H
#define SCHEME_EXACT_SOLUTION_H

#include "ExactSolution.h"

template <class SchemeType, int dim>
class SchemeExactSolution : public ExactSolution<typename SchemeType::VectorType,dim> {
  typedef ExactSolution<typename SchemeType::VectorType,dim> base;
public:
  typedef typename SchemeType::VectorType VectorType;
  typedef typename VectorType::InternalDataType DataType;
  typedef typename base::vec_grid_data_type vec_grid_data_type;
  typedef typename base::grid_data_type grid_data_type;
  typedef typename base::vec_grid_fct_type vec_grid_fct_type;  
  typedef typename base::grid_fct_type grid_fct_type;

  SchemeExactSolution(SchemeType &scheme) : base(), _scheme(scheme) {
    OutputFileName = "-";
    OutputIndices = "-";
  }
  virtual ~SchemeExactSolution() {}
  
  virtual void SetGrid(vec_grid_data_type& fvec, grid_data_type& workvec, 
                       const int& level, const double& t) {}

  virtual void register_at(ControlDevice& Ctrl, const std::string& prefix) {
    base::register_at(Ctrl,prefix);
    RegisterAt(base::LocCtrl,"OutputFileName",OutputFileName);
    RegisterAt(base::LocCtrl,"OutputIndices",OutputIndices);
  }
  virtual void register_at(ControlDevice& Ctrl) { register_at(Ctrl, ""); }

  virtual void ErrorNorm(vec_grid_fct_type& u, grid_fct_type& work, 
                         const double& t) {

    if (base::ENorm<1 || base::ENorm>3) 
      return;
    if (!(this->Set(u,work,t)))
      return;

    if (OutputIndices.c_str()[0] != '-') {
      int me = MY_PROC;
      std::ofstream outfile;
      std::ostream* out = 0;
      if (me==VizServer && OutputFileName.c_str()[0] != '-') {
        outfile.open(OutputFileName.c_str(), std::ios::out | std::ios::app);
        out = new std::ostream(outfile.rdbuf());
      }
      if (me == VizServer) {
        std::cout << "                         Output Error:";
        if (out) *out << t;
      }
      
      DataType *Error = (DataType *) 0;
      std::string delimiter = ",";
      std::string s = OutputIndices+delimiter;
      size_t pos = 0;
      std::string token;
      while ((pos = s.find(delimiter)) != std::string::npos) {
        token = s.substr(0, pos);
        int cnt = std::atoi( token.c_str() );
        s.erase(0, pos + delimiter.length());
        
        for (int lev=0; lev<=FineLevel(base::GH()); lev++) {
          int Time = CurrentTime(base::GH(),lev); 
          int TStep = TimeStep(u,lev);
          forall (u,Time,lev,c) 
            Scheme().Output(u(Time,lev,c), work(Time,lev,c), cnt);
          Scheme().Output(u(Time+TStep,lev,c), work(Time+TStep,lev,c), cnt);
          work(Time+TStep,lev,c).minus(work(Time,lev,c));
          end_forall            
        }
        
        this->CalculateNorm(work, Error);
        
        if (base::ENormAll) {
          std::cout << "   " << Error[0];
          if (out) *out << " " << Error[0];
        }
        else 
          for (int lev=0; lev<=FineLevel(base::GH()); lev++) {
            std::cout << "   L(" << lev << ")=" << Error[lev];
            if (out) *out << " " << Error[lev];   
          }
      }
      if (me == VizServer) {
        std::cout << std::endl;
        if (out) *out << std::endl;       
        outfile.close();
        delete out;
      }
      if (Error) delete [] Error;
    }
    base::ErrorNorm(u,work,t);
  }  

  inline SchemeType& Scheme() { return _scheme; }
  inline const SchemeType& Scheme() const { return _scheme; }

protected:
  SchemeType &_scheme;
  std::string OutputFileName, OutputIndices;
};


#endif

---