AMROC - Source code documentation

00001 #ifndef AMROC_AMRDQ_FLAGGINGREL_H
00002 #define AMROC_AMRDQ_FLAGGINGREL_H
00003 
00011 #ifndef AMRDQFlaggingRelName
00012 #define AMRDQFlaggingRel(dim)      name2(AMRDQFlaggingRel,dim)
00013 #define AMRDQFlaggingRelName
00014 #endif
00015 
00026 template <class VectorType, class FixupType, class FlagType>
00027 class AMRDQFlaggingRel(DIM) : 
00028   public AMRDQFlagging(DIM)<VectorType,FixupType,FlagType> {
00029   typedef AMRDQFlagging(DIM)<VectorType,FixupType,FlagType> flagging_type;    
00030   typedef Solver(DIM)<VectorType,FixupType,FlagType> solver_type;    
00031   typedef typename VectorType::InternalDataType DataType;
00032   typedef GridData(DIM)<VectorType>   vec_grid_data_type;
00033 
00034 public:
00035   AMRDQFlaggingRel(DIM)(solver_type& solver, const int quan) :
00036     flagging_type(solver, quan), _Quantities(quan) {
00037     DQTolGradRel = new DataType[NQuantities()];
00038     DQTolErrEstRel = new DataType[NQuantities()];
00039     DQScGradRel = new DataType[NQuantities()];
00040     DQScErrEstRel = new DataType[NQuantities()];
00041     DQMax = new DataType[NQuantities()];
00042     DQScalingGradRel = new DataType[NQuantities()];
00043     DQScalingErrEstRel = new DataType[NQuantities()];
00044     for (int k=0; k<NQuantities(); k++) {
00045       DQTolGradRel[k] = 0.0; DQTolErrEstRel[k] = 0.0; 
00046       DQScGradRel[k] = 0.0; DQScErrEstRel[k] = 0.0;
00047       DQMax[k] = 0.0; DQScalingGradRel[k] = 0.0; DQScalingErrEstRel[k] = 0.0;
00048     }
00049     _DQNFlagsRel = 2*NQuantities(); 
00050     SetDQBaseFlagRel(flagging_type::NFlags());
00051   }
00052 
00053   ~AMRDQFlaggingRel(DIM)() {
00054     delete [] DQTolGradRel; delete [] DQTolErrEstRel;
00055     delete [] DQScGradRel; delete [] DQScErrEstRel;
00056     delete [] DQMax; 
00057     delete [] DQScalingGradRel; delete [] DQScalingErrEstRel;
00058   }
00059 
00060   virtual void register_at(ControlDevice& Ctrl) {}
00061   virtual void register_at(ControlDevice& Ctrl,const string& prefix) {
00062     flagging_type::register_at(Ctrl,prefix);
00063     char Name[64];
00064     int d;
00065     for (d=0; d<NQuantities(); d++) {
00066       sprintf(Name,"DQMax(%d)",d+1);
00067       RegisterAt(LocCtrl,Name,DQMax[d]);
00068       sprintf(Name,"DQTolGradRel(%d)",d+1);
00069       RegisterAt(LocCtrl,Name,DQTolGradRel[d]);
00070       sprintf(Name,"DQTolErrEstRel(%d)",d+1);
00071       RegisterAt(LocCtrl,Name,DQTolErrEstRel[d]);
00072       sprintf(Name,"DQScGradRel(%d)",d+1);     
00073       RegisterAt(LocCtrl,Name,DQScGradRel[d]);
00074       sprintf(Name,"DQScErrEstRel(%d)",d+1);     
00075       RegisterAt(LocCtrl,Name,DQScErrEstRel[d]);
00076     }    
00077   }
00078   virtual void update() {
00079     flagging_type::update();
00080     bool errest = ErrorEstimation();
00081     int k;
00082     for (k=0; k<NQuantities(); k++) 
00083       if (DQTolErrEstRel[k]>0.0 && DQScErrEstRel[k]>0.0) {
00084         errest = true; break;
00085       }
00086     SetErrorEstimation(errest);
00087     for (k=0; k<NQuantities(); k++) {
00088       DQScalingGradRel[k] = DQScGradRel[k]*DQMax[k]; 
00089       DQScalingErrEstRel[k] = DQScErrEstRel[k]*DQMax[k]; 
00090     }      
00091   }
00092 
00093   virtual void SetFlags(const int Time, const int Level, double t, double dt) {
00094     FlagCellsByDifference(Time, Level);
00095     FlagCellsByDifferenceRel(Time, Level);
00096     FlagCellsDQByDifference(Time, Level, dt);
00097     FlagCellsDQByDifferenceRel(Time, Level, dt);
00098     if (ErrorEstimation() && Time>0) {
00099       //---------------------------------------------------------------
00100       // Take a step on the U(). U()(Time,Level) has to remain 
00101       // unchanged in a fractional step method!
00102       //---------------------------------------------------------------
00103       END_WATCH_FLAGGING
00104         Solver_().IntegrateLevel(U(), Time, Level, t, dt, false, 0.0, 1); 
00105       START_WATCH
00106       FlagCellsByError(Time, Level);
00107       FlagCellsByErrorRel(Time, Level);
00108       FlagCellsDQByError(Time, Level, dt);
00109       FlagCellsDQByErrorRel(Time, Level, dt);
00110     }
00111   }   
00112 
00113   virtual int NFlags() const { return (flagging_type::NFlags()+_DQNFlagsRel); }
00114 protected:  
00115   void FlagCellsDQByDifferenceRel(const int Time, const int Level, double dt) { 
00116     
00117     if (Time<0) return;
00118     int me = MY_PROC(GH); 
00119     int TStep = TimeStep(U(),Level);
00120     for (int k=1; k<=NQuantities(); k++) {
00121       if (DQScalingGradRel[k-1]<=0.0 || DQTolGradRel[k-1]<=0.0) continue;
00122       forall (U(),Time,Level,c)    
00123         //------------------------------------------------------------------------
00124         // Compute derived quantities from vector of state
00125         //------------------------------------------------------------------------
00126         f_dq_flag(FA(DIM,U()(Time,Level,c)), 
00127                   FA(DIM,Work()(Time,Level,c)),
00128                   FBA(Work()(Time,Level,c)),NEquations(),k,dt);
00129         Work()(Time+TStep,Level,c) = 1.0;
00130       end_forall
00131       FlagByDifferenceRel(Time, Level, Work(), DQTolGradRel[k-1], DQScalingGradRel[k-1], 
00132                           k+DQBaseFlagRel());
00133   
00134       if (PlotGrad() && Solver_().LastOutputTime()==Time) {
00135         char name[DAGHBktGFNameWidth+16];
00136         sprintf(name,"dqgradr_%d_%d",k,Time); 
00137         if (Level==0) {
00138           DataType MaxRel = MaxVal(Work(), Time, 0);
00139           if (me == VizServer) 
00140             cout << " *** Writing " << name << " Max on L0=" << MaxRel << endl;   
00141         }
00142         Write(Work(),Time+TStep,Level,DAGH_Double,name); 
00143 #ifndef DAGH_NO_MPI
00144         MPI_Barrier(comm_service::comm()); 
00145 #endif
00146       }
00147     } 
00148   }
00149 
00150   void FlagCellsDQByErrorRel(const int Time, const int Level, double dt) { 
00151     
00152     if (Time<=0) return;
00153     int me = MY_PROC(GH); 
00154     int TStep = TimeStep(U(),Level);
00155     int ShTStep = TimeStep(Ush(),Level);
00156     int myargc = sizeof(int); int myargs = 1;
00157     DataType Order = pow(2.0, 
00158        ((flagging_type::integrator_type&)Integrator_()).NMethodOrder()+1.0) - 2.0;
00159     for (int k=1; k<=NQuantities(); k++) {
00160       if (DQScalingErrEstRel[k-1]<=0.0 || DQTolErrEstRel[k-1]<=0.0) continue;
00161       forall (U(),Time,Level,c)    
00162         //------------------------------------------------------------------------
00163         // Compute derived quantities from vector of state
00164         //------------------------------------------------------------------------
00165         f_dq_flag(FA(DIM,U()(Time+TStep,Level,c)), 
00166                   FA(DIM,Work()(Time+TStep,Level,c)),
00167                   FBA(Work().interiorbbox(Time+TStep,Level,c)),
00168                   NEquations(),k,dt);
00169         f_dq_flag(FA(DIM,Ush()(Time,Level,c)), 
00170                   FA(DIM,Worksh()(Time,Level,c)),
00171                   FBA(Worksh().interiorbbox(Time,Level,c)),
00172                   NEquations(),k,dt);
00173         //------------------------------------------------------------------------
00174         // First update Worksh() from Work() 
00175         //------------------------------------------------------------------------
00176         f_restrict_amr(FA(DIM,(vec_grid_data_type &)Work()(Time+TStep,Level,c)),
00177                        FA(DIM,(vec_grid_data_type &)Worksh()(Time+ShTStep,Level,c)),
00178                        FBA(Work().interiorbbox(Time+TStep,Level,c)), 
00179                        (char *)&myargs, &myargc); 
00180         //------------------------------------------------------------------------
00181         // Estimate error by Richardson-Extrapolation
00182         //------------------------------------------------------------------------
00183         Worksh()(Time+ShTStep,Level,c).minus(Worksh()(Time,Level,c));
00184         Worksh()(Time+ShTStep,Level,c).divide(Order);
00185       end_forall
00186       FlagByErrorEstimationRel(Time, Level, Worksh(), DQTolErrEstRel[k-1], 
00187                                DQScalingErrEstRel[k-1], k+NQuantities()+DQBaseFlagRel());
00188 
00189       if (PlotErrEst() && Solver_().LastOutputTime()==Time) {
00190         char name[DAGHBktGFNameWidth+16];
00191         sprintf(name,"dqerrestr_%d_%d",k,Time); 
00192         if (Level==0) {
00193           DataType MaxRel = MaxVal(Worksh(), Time, 0);
00194           if (me == VizServer) 
00195             cout << " *** Writing " << name << " Max on L0=" << MaxRel << endl;   
00196         }
00197         Write(Worksh(),Time+ShTStep,Level,DAGH_Double,name); 
00198 #ifndef DAGH_NO_MPI
00199         MPI_Barrier(comm_service::comm()); 
00200 #endif
00201       }
00202     }
00203   }
00204 
00205   inline const int& NQuantities() const { return _Quantities; }
00206   inline const int& DQBaseFlagRel() const { return _DQBaseFlagRel; }
00207   void SetDQBaseFlagRel(const int bf) { _DQBaseFlagRel = bf; }
00208 
00209 protected:
00210   DataType *DQTolGradRel, *DQTolErrEstRel, *DQScGradRel, *DQScErrEstRel; 
00211   DataType *DQMax, *DQScalingGradRel, *DQScalingErrEstRel;
00212   int _Quantities;
00213   int _DQBaseFlagRel, _DQNFlagsRel;
00214 };
00215 
00216 
00217 
00218 #endif
AMRDQFlaggingRel.h

last update: 06/01/04
