Blockstructured Adaptive Mesh Refinement in object-oriented C++

---

       

---

AMRDQFlagging.h

Go to the documentation of this file.

#ifndef AMROC_AMRDQ_FLAGGING_H
#define AMROC_AMRDQ_FLAGGING_H

#ifndef AMRDQFlaggingName
#define AMRDQFlagging(dim)      name2(AMRDQFlagging,dim)
#define AMRDQFlaggingName
#endif

template <class VectorType, class FixupType, class FlagType>
class AMRDQFlagging(DIM) : 
  public AMRFlaggingRel(DIM)<VectorType,FixupType,FlagType> {
  typedef typename VectorType::InternalDataType DataType;
  typedef GridData(DIM)<VectorType>   vec_grid_data_type;
 
  typedef AMRFlaggingRel(DIM)<VectorType,FixupType,FlagType> flagging_type;    
  typedef Solver(DIM)<VectorType,FixupType,FlagType> solver_type;    

public:
  AMRDQFlagging(DIM)(solver_type& solver, const int quan) :
    flagging_type(solver), _Quantities(quan) {
    DQTolerance   = new DataType[NQuantities()];
    DQToleranceSp = new DataType[NQuantities()];
    for (int k=0; k<NQuantities(); k++) {
      DQTolerance[k]   = 0.0;
      DQToleranceSp[k] = 0.0;
    }
    _DQNFlags = 2*NQuantities(); 
    SetDQBaseFlag(flagging_type::NFlags());
  }

  ~AMRDQFlagging(DIM)() {
    delete [] DQTolerance;
    delete [] DQToleranceSp;
  }

  virtual void register_at(ControlDevice& Ctrl) {}
  virtual void register_at(ControlDevice& Ctrl,const string& prefix) {
    flagging_type::register_at(Ctrl,prefix);
    int d;
    char ToleranceName[16];
    for (d=0; d<NQuantities(); d++) {
      sprintf(ToleranceName,"DQTolGrad(%d)",d+1);
      RegisterAt(LocCtrl,ToleranceName,DQToleranceSp[d]);
      sprintf(ToleranceName,"DQTolErrEst(%d)",d+1);
      RegisterAt(LocCtrl,ToleranceName,DQTolerance[d]);
    }    
  }
  virtual void update() {
    flagging_type::update();
    bool est = ErrorEstimation();
    for (int k=0; k<NQuantities(); k++) 
      est = est || (DQTolerance[k] > 0.0);
    SetErrorEstimation(est);
  }

  virtual void SetFlags(const int Time, const int Level, double t, double dt) {
    FlagCellsByDifference(Time, Level);
    FlagCellsByDifferenceRel(Time, Level);
    FlagCellsDQByDifference(Time, Level, dt);
    if (ErrorEstimation() && Time>0) {
      //---------------------------------------------------------------
      // Take a step on the U(). U()(Time,Level) has to remain 
      // unchanged in a fractional step method!
      //---------------------------------------------------------------
      END_WATCH_FLAGGING
        Solver_().IntegrateLevel(U(), Time, Level, t, dt, false, 0.0, 1); 
      START_WATCH
      FlagCellsByError(Time, Level);
      FlagCellsByErrorRel(Time, Level);
      FlagCellsDQByError(Time, Level, dt);
    }
  }   

  virtual int NFlags() const { return (flagging_type::NFlags()+_DQNFlags); }
protected:  
  void FlagCellsDQByDifference(const int Time, const int Level, double dt) { 
    
    if (Time < 0) return;
    int me = MY_PROC(GH); 
    int TStep = TimeStep(U(),Level);
    for (int k=1; k<=NQuantities(); k++) {
      if (DQToleranceSp[k-1] <= 0.0) continue;
      forall (U(),Time,Level,c)    
        //------------------------------------------------------------------------
        // Compute derived quantities from vector of state
        //------------------------------------------------------------------------
        f_dq_flag(FA(DIM,U()(Time,Level,c)), 
                  FA(DIM,Work()(Time,Level,c)),
                  FBA(Work()(Time,Level,c)),NEquations(),k,dt);
        Work()(Time+TStep,Level,c) = 0.0;       
      end_forall
      FlagByDifference(Time, Level, Work(), DQToleranceSp[k-1], k+DQBaseFlag());

      if (PlotGrad() && Solver_().LastOutputTime()==Time) {
        char name[DAGHBktGFNameWidth+16];
        sprintf(name,"dqgrad_%d_%d",k,Time); 
        if (Level==0 && me==VizServer) 
          cout << " *** Writing " << name << endl;   
        Write(Work(),Time+TStep,Level,DAGH_Double,name); 
#ifndef DAGH_NO_MPI
        MPI_Barrier(comm_service::comm()); 
#endif
      }
    } 
  }

  void FlagCellsDQByError(const int Time, const int Level, double dt) { 
    
    if (Time <= 0) return;
    int me = MY_PROC(GH); 
    int TStep = TimeStep(U(),Level);
    int ShTStep = TimeStep(Ush(),Level);
    int myargc = sizeof(int); int myargs = 1;
    DataType Order = pow(2.0, 
       ((flagging_type::integrator_type&)Integrator_()).NMethodOrder()+1.0) - 2.0;
    for (int k=1; k<=NQuantities(); k++) {
      if (DQTolerance[k-1] <= 0.0) continue;
      forall (U(),Time,Level,c)    
        //------------------------------------------------------------------------
        // Compute derived quantities from vector of state
        //------------------------------------------------------------------------
        f_dq_flag(FA(DIM,U()(Time+TStep,Level,c)), 
                  FA(DIM,Work()(Time+TStep,Level,c)),
                  FBA(Work().interiorbbox(Time+TStep,Level,c)),
                  NEquations(),k,dt);
        f_dq_flag(FA(DIM,Ush()(Time,Level,c)), 
                  FA(DIM,Worksh()(Time,Level,c)),
                  FBA(Worksh().interiorbbox(Time,Level,c)),
                  NEquations(),k,dt);
        //------------------------------------------------------------------------
        //       // First update Worksh() from Work()
        //------------------------------------------------------------------------
        f_restrict_amr(FA(DIM,(vec_grid_data_type &)Work()(Time+TStep,Level,c)),
                       FA(DIM,(vec_grid_data_type &)Worksh()(Time+ShTStep,Level,c)),
                       FBA(Work().interiorbbox(Time+TStep,Level,c)), 
                       (char *)&myargs, &myargc); 
        //------------------------------------------------------------------------
        // Estimate error by Richardson-Extrapolation
        //------------------------------------------------------------------------
        Worksh()(Time+ShTStep,Level,c).minus(Worksh()(Time,Level,c));
        Worksh()(Time+ShTStep,Level,c).divide(Order);
      end_forall
      FlagByErrorEstimation(Time, Level, Worksh(), DQTolerance[k-1], 
                            k+NQuantities()+DQBaseFlag());

      if (PlotErrEst() && Solver_().LastOutputTime()==Time) {
        char name[DAGHBktGFNameWidth+16];
        sprintf(name,"dqerrest_%d_%d",k,Time); 
        if (Level==0 && me==VizServer) 
          cout << " *** Writing " << name << endl;   
        Write(Worksh(),Time+ShTStep,Level,DAGH_Double,name); 
#ifndef DAGH_NO_MPI
        MPI_Barrier(comm_service::comm()); 
#endif
      }
    } 
  }

  inline const int& NQuantities() const { return _Quantities; }
  inline const int& DQBaseFlag() const { return _DQBaseFlag; }
  void SetDQBaseFlag(const int bf) { _DQBaseFlag = bf; }

protected:
  DataType *DQTolerance, *DQToleranceSp;
  int _Quantities;
  int _DQBaseFlag;
  int _DQNFlags;
};



#endif

---

Quickstart     Users Guide     Programmers Reference     Installation      Examples     Download

---

AMROC Main      Home      Contact

last update: 06/01/04