AMROC - Source code documentation

00001 #ifndef AMROC_AMRFLAGGING_RELBASE_H
00002 #define AMROC_AMRFLAGGING_RELBASE_H
00003 
00010 #include "AMRFlaggingRelBaseInterface.h"
00011 
00012 #ifndef AMRFlaggingRelBaseName
00013 #define AMRFlaggingRelBase(dim)      name2(AMRFlaggingRelBase,dim)
00014 #define AMRFlaggingRelBaseName
00015 #endif
00016 
00027 template <class VectorType, class FixupType, class FlagType>
00028 class AMRFlaggingRelBase(DIM) : 
00029   public AMRFlagging(DIM)<VectorType,FixupType,FlagType>, 
00030   public AMRFlaggingRelBaseInterface(DIM)<VectorType,FlagType> {
00031   typedef AMRFlagging(DIM)<VectorType,FixupType,FlagType> flagging_type;    
00032   typedef typename VectorType::InternalDataType DataType;
00033   typedef Solver(DIM)<VectorType,FixupType,FlagType> solver_type;    
00034   typedef GridData(DIM)<FlagType> flag_data_type;
00035   typedef GridData(DIM)<VectorType>   vec_grid_data_type;
00036 
00037 public:
00038   AMRFlaggingRelBase(DIM)(solver_type& solver) : 
00039     flagging_type(solver), _BaseFlagRel(0) {
00040     TolGradRel = 0.0; TolErrEstRel = 0.0; ScGradRel = 0.0; ScErrEstRel = 0.0; 
00041     MaxRel = 0.0; ScalingGradRel = 0.0; ScalingErrEstRel = 0.0;
00042     _NFlagsRel = 2*NEquations();
00043     SetBaseFlagRel(flagging_type::NFlags());
00044   }
00045 
00046   virtual void register_at(ControlDevice& Ctrl) {}
00047   virtual void register_at(ControlDevice& Ctrl,const string& prefix) {
00048     flagging_type::register_at(Ctrl, prefix);
00049     char Name[64];
00050     for (int d=0; d<NEquations(); d++) {
00051       sprintf(Name,"Max(%d)",d+1);
00052       RegisterAt(LocCtrl,Name,MaxRel(d));
00053       sprintf(Name,"TolGradRel(%d)",d+1);
00054       RegisterAt(LocCtrl,Name,TolGradRel(d));
00055       sprintf(Name,"TolErrEstRel(%d)",d+1);
00056       RegisterAt(LocCtrl,Name,TolErrEstRel(d));
00057       sprintf(Name,"ScGradRel(%d)",d+1);     
00058       RegisterAt(LocCtrl,Name,ScGradRel(d));
00059       sprintf(Name,"ScErrEstRel(%d)",d+1);     
00060       RegisterAt(LocCtrl,Name,ScErrEstRel(d));
00061     }
00062   }
00063   virtual void update() {
00064     flagging_type::update();
00065     ScalingGradRel = ScGradRel*MaxRel;
00066     ScalingErrEstRel = ScErrEstRel*MaxRel;
00067     bool errest = ErrorEstimation();
00068     for (int d=0; d<NEquations(); d++) 
00069       errest = errest || (TolErrEstRel(d)>0.0 && ScErrEstRel(d)>0.0); 
00070     SetErrorEstimation(errest);
00071   }
00072 
00073   virtual void SetFlags(const int Time, const int Level, double t, double dt) {
00074     FlagCellsByDifference(Time, Level);
00075     FlagCellsByDifferenceRel(Time, Level);
00076     if (ErrorEstimation() && Time>0) {
00077       //---------------------------------------------------------------
00078       // Take a step on the U(). U()(Time,Level) has to remain 
00079       // unchanged in a fractional step method!
00080       //---------------------------------------------------------------
00081       END_WATCH_FLAGGING
00082         Solver_().IntegrateLevel(U(), Time, Level, t, dt, false, 0.0, 1); 
00083       START_WATCH
00084       FlagCellsByError(Time, Level);
00085       FlagCellsByErrorRel(Time, Level);
00086     }
00087   }   
00088 
00089   virtual int NFlags() const { return (flagging_type::NFlags()+_NFlagsRel); }
00090 protected:  
00091   DataType phi(DataType theta) {
00092     return (4.*theta/((1.+theta)*(1.+theta)));
00093   }
00094 //   DataType phi(DataType theta) {
00095 //     return (theta>1.0 ? 1.0/theta : theta);
00096 //   }
00097 
00098   virtual void FlagCellsByDifferenceRel(const int Time, const int Level) { 
00099 
00100     if (Time<0) return;
00101     int me = MY_PROC(GH); 
00102     int TStep = TimeStep(U(),Level);
00103     for (int k=0; k<NEquations(); k++) {
00104       if (ScalingGradRel(k)<=0.0 || TolGradRel(k)<=0.0) continue;
00105       forall(U(),Time,Level,c)
00106         equals_from(Work()(Time,Level,c), U()(Time,Level,c), k);
00107         Work()(Time+TStep,Level,c) = 1.0;
00108       end_forall      
00109       FlagByDifferenceRel(Time, Level, Work(), TolGradRel(k), ScalingGradRel(k), 
00110                           k+1+BaseFlagRel());
00111 
00112       if (PlotGrad() && Solver_().LastOutputTime()==Time) {
00113         char name[DAGHBktGFNameWidth+16];
00114         sprintf(name,"gradr_%d_%d",k+1,Time); 
00115         if (Level==0) {
00116           DataType MaxRel = MaxVal(Work(), Time, 0);
00117           if (me == VizServer) 
00118             cout << " *** Writing " << name << " Max on L0=" << MaxRel << endl;   
00119         }
00120         Write(Work(),Time+TStep,Level,DAGH_Double,name); 
00121 #ifndef DAGH_NO_MPI
00122         MPI_Barrier(comm_service::comm()); 
00123 #endif
00124       }
00125     }
00126   }
00127   
00128   virtual void FlagCellsByErrorRel(const int Time, const int Level) { 
00129     
00130     if (Time<=0) return;
00131     int me = MY_PROC(GH); 
00132     int TStep = TimeStep(U(),Level);
00133     int ShTStep = TimeStep(Ush(),Level);
00134     int myargc = sizeof(int); int myargs = NEquations();
00135     forall (U(),Time,Level,c)      
00136       //------------------------------------------------------------------------
00137       // First update Ush() from U()
00138       //------------------------------------------------------------------------
00139       f_restrict_amr(FA(DIM,U()(Time+TStep,Level,c)),
00140                      FA(DIM,Ush()(Time+ShTStep,Level,c)),
00141                      FBA(U().interiorbbox(Time+TStep,Level,c)), 
00142                      (char *)&myargs, &myargc);     
00143     end_forall
00144 
00145     //------------------------------------------------------------------------
00146     // Estimate error by Richardson-Extrapolation
00147     //------------------------------------------------------------------------
00148     DataType Order = pow(2.0,Integrator_().NMethodOrder()+1.0) - 2.0;
00149     for (int k=0; k<NEquations(); k++) {
00150       if (ScalingErrEstRel(k)<=0.0 || TolErrEstRel(k)<=0) continue;
00151       forall(Ush(),Time,Level,c)
00152         equals_from(Worksh()(Time,Level,c), Ush()(Time,Level,c), k);
00153         equals_from(Worksh()(Time+ShTStep,Level,c), Ush()(Time+ShTStep,Level,c), k);
00154         minus_from(Worksh()(Time+ShTStep,Level,c), Ush()(Time,Level,c), k);
00155         Worksh()(Time+ShTStep,Level,c).divide(Order);
00156       end_forall      
00157       FlagByErrorEstimationRel(Time, Level, Worksh(), TolErrEstRel(k), ScalingErrEstRel(k), 
00158                                k+1+NEquations()+BaseFlagRel());
00159 
00160       if (PlotErrEst() && Solver_().LastOutputTime()==Time) {
00161         char name[DAGHBktGFNameWidth+16];
00162         sprintf(name,"errestr_%d_%d",k+1,Time); 
00163         if (Level==0) {
00164           DataType MaxRel = MaxVal(Worksh(), Time, 0);
00165           if (me == VizServer) 
00166             cout << " *** Writing " << name << " Max on L0=" << MaxRel << endl;   
00167         }
00168         Write(Worksh(),Time+ShTStep,Level,DAGH_Double,name); 
00169 #ifndef DAGH_NO_MPI
00170         MPI_Barrier(comm_service::comm()); 
00171 #endif
00172       }
00173     }
00174   } 
00175 
00176   inline const int& BaseFlagRel() const { return _BaseFlagRel; }
00177   void SetBaseFlagRel(const int bf) { _BaseFlagRel = bf; }
00178 
00179 protected:
00180   VectorType TolGradRel, ScGradRel, TolErrEstRel, ScErrEstRel;
00181   VectorType MaxRel, ScalingGradRel, ScalingErrEstRel;
00182   int _BaseFlagRel, _NFlagsRel;
00183 };
00184 
00185 
00186 #endif
AMRFlaggingRelBase.h

last update: 06/01/04
