AMROC Main     Blockstructured Adaptive Mesh Refinement in object-oriented C++


Main Page   Class Hierarchy   Compound List   File List  

AMRFlaggingBase.h

Go to the documentation of this file.
00001 #ifndef AMROC_AMRFLAGGING_BASE_H
00002 #define AMROC_AMRFLAGGING_BASE_H
00003 
00011 #include "Flagging.h"
00012 #include "AMRFlaggingBaseInterface.h"
00013 
00014 #ifndef AMRFlaggingBaseName
00015 #define AMRFlaggingBase(dim)      name2(AMRFlaggingBase,dim)
00016 #define AMRFlaggingBaseName
00017 #endif
00018 
00027 template <class VectorType, class FixupType, class FlagType>
00028 class AMRFlaggingBase(DIM) : public Flagging(DIM)<VectorType,FixupType,FlagType>,
00029   public AMRFlaggingBaseInterface(DIM)<VectorType,FlagType> {
00030   typedef typename VectorType::InternalDataType DataType;
00031   typedef Solver(DIM)<VectorType,FixupType,FlagType> solver_type;    
00032   typedef GridData(DIM)<FlagType> flag_data_type;
00033   typedef GridData(DIM)<VectorType>   vec_grid_data_type;
00034  
00035 public:
00036   AMRFlaggingBase(DIM)(solver_type& solver) :
00037     Flagging(DIM)<VectorType,FixupType,FlagType>(solver), _PlotGrad(0), _PlotErrEst(0), 
00038     _BaseFlag(0) {
00039     Tolerance   = 0.0;
00040     ToleranceSp = 0.0;
00041     _NFlags = 2*NEquations();
00042   }
00043 
00044   virtual void register_at(ControlDevice& Ctrl) {}
00045   virtual void register_at(ControlDevice& Ctrl,const string& prefix) {
00046     LocCtrl = Ctrl.getSubDevice(prefix+"AMRFlagging");
00047     char ToleranceName[16];
00048     for (int d=0; d<NEquations(); d++) {
00049       sprintf(ToleranceName,"TolGrad(%d)",d+1);
00050       RegisterAt(LocCtrl,ToleranceName,ToleranceSp(d));
00051       sprintf(ToleranceName,"TolErrEst(%d)",d+1);
00052       RegisterAt(LocCtrl,ToleranceName,Tolerance(d));
00053       RegisterAt(LocCtrl,"PlotGrad",_PlotGrad);
00054       RegisterAt(LocCtrl,"PlotErrEst",_PlotErrEst);
00055     }    
00056   }
00057   virtual void update() {
00058     SetErrorEstimation(Tolerance>0.0 || ErrorEstimation());
00059   }
00060 
00061   virtual void SetFlags(const int Time, const int Level, double t, double dt) {
00062     FlagCellsByDifference(Time, Level);
00063     if (ErrorEstimation() && Time>0) {
00064       //---------------------------------------------------------------
00065       // Take a step on the U(). U()(Time,Level) has to remain 
00066       // unchanged in a fractional step method!
00067       //---------------------------------------------------------------
00068       END_WATCH_FLAGGING
00069         Solver_().IntegrateLevel(U(), Time, Level, t, dt, false, 0.0, 1); 
00070       START_WATCH
00071       FlagCellsByError(Time, Level);
00072     }
00073   }   
00074 
00075   virtual int NFlags() const { return _NFlags; }
00076 protected:  
00077   virtual void FlagCellsByDifference(const int Time, const int Level) { 
00078 
00079     if (Time < 0) return;
00080     int me = MY_PROC(GH); 
00081     int TStep = TimeStep(U(),Level);
00082     for (int k=0; k<NEquations(); k++) {
00083       if (ToleranceSp(k)<=0.0) continue;
00084       forall(U(),Time,Level,c)
00085         equals_from(Work()(Time,Level,c), U()(Time,Level,c), k);
00086         Work()(Time+TStep,Level,c) = 0.0;
00087       end_forall      
00088       FlagByDifference(Time, Level, Work(), ToleranceSp(k), k+1+BaseFlag());
00089 
00090       if (PlotGrad() && Solver_().LastOutputTime()==Time) {
00091         char name[DAGHBktGFNameWidth+16];
00092         sprintf(name,"grad_%d_%d",k+1,Time); 
00093         if (Level==0 && me==VizServer) 
00094           cout << " *** Writing " << name << endl;   
00095         Write(Work(),Time+TStep,Level,DAGH_Double,name); 
00096 #ifndef DAGH_NO_MPI
00097         MPI_Barrier(comm_service::comm()); 
00098 #endif
00099       }
00100     }
00101   }
00102 
00103   virtual void FlagCellsByError(const int Time, const int Level) { 
00104     
00105     if (Time <= 0) return;
00106     int me = MY_PROC(GH); 
00107     int TStep = TimeStep(U(),Level);
00108     int ShTStep = TimeStep(Ush(),Level);
00109     int myargc = sizeof(int); int myargs = NEquations();
00110     forall (U(),Time,Level,c)      
00111       //------------------------------------------------------------------------
00112       // First update Ush() from U()
00113       //------------------------------------------------------------------------
00114       f_restrict_amr(FA(DIM,U()(Time+TStep,Level,c)),
00115                      FA(DIM,Ush()(Time+ShTStep,Level,c)),
00116                      FBA(U().interiorbbox(Time+TStep,Level,c)), 
00117                      (char *)&myargs, &myargc); 
00118     end_forall
00119 
00120     //------------------------------------------------------------------------
00121     // Estimate error by Richardson-Extrapolation
00122     //------------------------------------------------------------------------
00123     DataType Order = pow(2.0,Integrator_().NMethodOrder()+1.0) - 2.0;
00124     for (int k=0; k<NEquations(); k++) {
00125       if (Tolerance(k)<=0.0) continue;
00126       forall(Ush(),Time,Level,c)
00127         equals_from(Worksh()(Time+ShTStep,Level,c), Ush()(Time+ShTStep,Level,c), k);
00128         minus_from(Worksh()(Time+ShTStep,Level,c), Ush()(Time,Level,c), k);
00129         Worksh()(Time+ShTStep,Level,c).divide(Order);
00130       end_forall      
00131       FlagByErrorEstimation(Time, Level, Worksh(), Tolerance(k), k+1+NEquations()+BaseFlag());
00132 
00133       if (PlotErrEst() && Solver_().LastOutputTime()==Time) {
00134         char name[DAGHBktGFNameWidth+16];
00135         sprintf(name,"errest_%d_%d",k+1,Time); 
00136         if (Level==0 && me==VizServer) 
00137           cout << " *** Writing " << name << endl;   
00138         Write(Worksh(),Time+ShTStep,Level,DAGH_Double,name); 
00139 #ifndef DAGH_NO_MPI
00140         MPI_Barrier(comm_service::comm()); 
00141 #endif
00142       }
00143     }
00144   } 
00145 
00146   inline const int& PlotGrad() const { return _PlotGrad; }
00147   inline const int& PlotErrEst() const { return _PlotErrEst; }
00148   inline const int& BaseFlag() const { return _BaseFlag; }
00149   void SetBaseFlag(const int bf) { _BaseFlag = bf; }
00150 
00151 protected:
00152   VectorType Tolerance, ToleranceSp;
00153   int _PlotGrad, _PlotErrEst;
00154   int _BaseFlag, _NFlags;
00155 };
00156 
00157 #endif


Quickstart     Users Guide     Programmers Reference     Installation      Examples     Download



AMROC Main      Home      Contact
last update: 06/01/04