ExactSolution.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_EXACT_SOLUTION_H
#define AMROC_EXACT_SOLUTION_H

template <class VectorType, int dim>
class ExactSolution : public AMRBase<VectorType,dim> {
  typedef AMRBase<VectorType,dim> base;
  typedef typename VectorType::InternalDataType DataType;

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

  ExactSolution() : base(), ENorm(0), ENormCutOut(1), ENormAll(1), ENormComp(1) {} 

  virtual ~ExactSolution() {}

  virtual void SetGrid(vec_grid_data_type& gd, const int& level, double t) {} 

  virtual void register_at(ControlDevice& Ctrl, const std::string& prefix) {
    LocCtrl = Ctrl.getSubDevice(prefix+"Error");    
    RegisterAt(LocCtrl,"Norm",ENorm);    
    RegisterAt(LocCtrl,"NoRestrict",ENormCutOut);
    RegisterAt(LocCtrl,"CombineLevels",ENormAll);
    RegisterAt(LocCtrl,"ComponentWise",ENormComp);
  }
  virtual void register_at(ControlDevice& Ctrl) { register_at(Ctrl, ""); }
  virtual void init() {}
  virtual void update() {}
  virtual void finish() {}

  //---------------------------------------------------------------------- 
  // Calculate Error Norm for conservative variables 
  //---------------------------------------------------------------------- 
  virtual void ErrorNorm(vec_grid_fct_type& u, grid_fct_type& work, 
                         const double& t) {
    if (ENorm<1 || ENorm>3) 
      return;
    int me = MY_PROC;
    if (me == VizServer) 
      std::cout << "                         Error:";
    for (int lev=0; lev<=FineLevel(base::GH()); lev++) {
      int Time = CurrentTime(base::GH(),lev); 
      int TStep = TimeStep(u,lev);
      forall (u,Time+TStep,lev,c) 
        SetGrid(u(Time+TStep,lev,c), lev, t);
        u(Time+TStep,lev,c).minus(u(Time,lev,c));
      end_forall      
    }

    DataType *Error = (DataType *) 0;
    if (ENormComp) {
      for (int i=0; i<base::NEquations(); i++) { 
        for (int lev=0; lev<=FineLevel(base::GH()); lev++) {
          int Time = CurrentTime(base::GH(),lev); 
          int TStep = TimeStep(u,lev);
          forall (u,Time+TStep,lev,c) 
            equals_from(work(Time+TStep,lev,c), u(Time+TStep,lev,c), i);
          end_forall 
        }       
        CalculateNorm(work, Error);
        if (ENormAll) std::cout << "   " << Error[0];
        else for (int lev=0; lev<=FineLevel(base::GH()); lev++) 
          std::cout << "   L(" << i << "," << lev << ")=" << Error[lev];
      }
    }
    else {
      CalculateNorm(u, Error);
      if (ENormAll) std::cout << "   " << Error[0];
      else for (int lev=0; lev<=FineLevel(base::GH()); lev++) 
        std::cout << "   L(" << lev << ")=" << Error[lev];
    }
    if (me == VizServer) 
      std::cout << std::endl;
    if (Error) 
      delete [] Error;
  }  

  virtual void CalculateNorm(vec_grid_fct_type& u, double*& Error) {
    if (ENormAll)
      Error = new double[1];
    else
      Error = new double[FineLevel(base::GH())];

    Error[0] = 0.0;
    for (int lev=0; lev<=FineLevel(base::GH()); lev++) {
      int Time = CurrentTime(base::GH(),lev); 
      int TStep = TimeStep(u,lev), FTStep=0;
      if (lev<FineLevel(base::GH()))
        FTStep = TimeStep(u,lev+1);
      
      forall (u,Time+TStep,lev,c) 
        if (lev<FineLevel(base::GH()) && ENormCutOut) {
          forall (u,Time+FTStep,lev+1,c2) 
            BBox bb(coarsen(u.interiorbbox(Time+FTStep,lev+1,c2),base::GH().refinefactor(lev)));
            u(Time+TStep,lev,c).equals(VectorType(0.0),bb);
          end_forall      
        }
      end_forall      

      double ErrorLev = Norm(u,Time+TStep,lev,ENorm);
      if (ENormAll) {
        switch (ENorm) {
        case DAGHNormL1:
        case DAGHNormL2:
          Error[0] += ErrorLev;
          break;
        case DAGHNormMax: 
          Error[0] = Max(Error[0], ErrorLev);
          break;
        default:
          break;
        }
      }
      else 
        Error[lev] = ErrorLev;
    }
  }  

  virtual void CalculateNorm(grid_fct_type& u, double*& Error) {
    if (ENormAll)
      Error = new double[1];
    else
      Error = new double[FineLevel(base::GH())];

    Error[0] = 0.0;
    for (int lev=0; lev<=FineLevel(base::GH()); lev++) {
      int Time = CurrentTime(base::GH(),lev); 
      int TStep = TimeStep(u,lev), FTStep=0;
      if (lev<FineLevel(base::GH()))
        FTStep = TimeStep(u,lev+1);
      
      forall (u,Time+TStep,lev,c) 
        if (lev<FineLevel(base::GH()) && ENormCutOut) {
          forall (u,Time+FTStep,lev+1,c2) 
            BBox bb(coarsen(u.interiorbbox(Time+FTStep,lev+1,c2),base::GH().refinefactor(lev)));
            u(Time+TStep,lev,c).equals(DataType(0.0),bb);
          end_forall      
        }
      end_forall      

      double ErrorLev = Norm(u,Time+TStep,lev,ENorm);
      if (ENormAll) {
        switch (ENorm) {
        case DAGHNormL1:
        case DAGHNormL2:
          Error[0] += ErrorLev;
          break;
        case DAGHNormMax: 
          Error[0] = Max(Error[0], ErrorLev);
          break;
        default:
          break;
        }
      }
      else 
        Error[lev] = ErrorLev;
    }
  }  

 protected:
  int ENorm, ENormCutOut, ENormAll, ENormComp;
  ControlDevice LocCtrl;  
};


#endif

---