amroc/amr/AMRGFMSolver.h

Go to the documentation of this file.
00001 // -*- C++ -*-
00002 
00003 // Copyright (C) 2003-2007 California Institute of Technology
00004 // Ralf Deiterding, ralf@cacr.caltech.edu
00005 
00006 #ifndef AMROC_GFMSOLVER_H
00007 #define AMROC_GFMSOLVER_H
00008 
00016 #include "AMRSolver.h"
00017 #include "GhostFluidMethod.h"
00018 #include "Criteria/RefinePhi.h"
00019 #include "Criteria/ResolvePhi.h"
00020 #include "Criteria/UnflagPhi.h"
00021 
00028 template <class VectorType, class FixupType, class FlagType, int dim>
00029 class AMRGFMSolver : public AMRSolver<VectorType,FixupType,FlagType,dim> {
00030   typedef typename VectorType::InternalDataType DataType;
00031   typedef AMRSolver<VectorType,FixupType,FlagType,dim> base;
00032   typedef AMRGFMSolver<VectorType,FixupType,FlagType,dim> self;
00033 public:
00034   typedef GhostFluidMethod<VectorType,dim> gfm_type;
00035   typedef typename base::vec_grid_fct_type vec_grid_fct_type;  
00036   typedef typename base::vec_grid_data_type vec_grid_data_type;  
00037   typedef typename base::grid_fct_type grid_fct_type;
00038   typedef typename base::grid_data_type grid_data_type;
00039 
00040   typedef typename base::integrator_type integrator_type;
00041   typedef typename base::initial_condition_type initial_condition_type;
00042   typedef typename base::boundary_conditions_type boundary_conditions_type;
00043   typedef typename base::leveltransfer_type leveltransfer_type;
00044   typedef typename base::flagging_type flagging_type;
00045   typedef typename base::fixup_type fixup_type;
00046 
00047   typedef GFRecomposeSpecificFunc<self,VectorType,dim> gfm_recompose_functor_type;
00048 
00049   typedef GridData<bool,dim>  bool_grid_data_type;
00050   typedef GridFunction<bool,dim> bool_grid_fct_type;  
00051 
00052   AMRGFMSolver(integrator_type& integ, 
00053                initial_condition_type& init,
00054                boundary_conditions_type& bc) : base(integ, init, bc), _MaxRecomposeLevel(-1),
00055                                                _RecoverExterior(0), _RecoverValue(0), _RecoveryValue(0.) {      
00056     _bf = (bool_grid_fct_type *) 0;
00057     _bf_sh = (bool_grid_fct_type *) 0;
00058     _GFMRecomposeFunc = (gfm_recompose_functor_type *) 0;
00059     _GFM = (gfm_type**) 0;
00060     _nGFM = 0;
00061   }
00062 
00063   virtual ~AMRGFMSolver() {    
00064     if (_bf) delete _bf;
00065     if (_bf_sh) delete _bf_sh;
00066     if (_GFMRecomposeFunc) delete _GFMRecomposeFunc;
00067   }
00068 
00069   virtual void init() {
00070     base::Flagging_().AddCriterion(new RefinePhi<VectorType,FixupType,FlagType,dim>(*this));
00071     base::Flagging_().AddCriterion(new UnflagPhi<VectorType,FixupType,FlagType,dim>(*this));    
00072     base::Flagging_().AddCriterion(new ResolvePhi<VectorType,FixupType,FlagType,dim>(*this,false));
00073     base::Flagging_().AddCriterion(new ResolvePhi<VectorType,FixupType,FlagType,dim>(*this,true));
00074     base::init();
00075     for (register int n=0; n<NGFM(); n++) 
00076       GFM(n).init();
00077   }
00078 
00079   virtual void register_at(ControlDevice& Ctrl, const std::string& prefix) {
00080     base::register_at(Ctrl,prefix);
00081     for (register int n=0; n<NGFM(); n++) {
00082       char text[5];
00083       std::sprintf(text,"%d-",n+1);
00084       GFM(n).register_at(base::LocCtrl, std::string(text));
00085     }
00086     RegisterAt(base::LocCtrl,"RecoverExterior",_RecoverExterior);
00087     RegisterAt(base::LocCtrl,"RecoverValue",_RecoverValue);
00088   } 
00089   virtual void register_at(ControlDevice& Ctrl) {
00090     register_at(Ctrl, "");
00091   }
00092 
00093   virtual void update() {
00094     base::update();
00095     for (register int n=0; n<NGFM(); n++) 
00096       GFM(n).update();
00097   }
00098   virtual void finish() {
00099     if (_bf) {
00100       delete _bf;
00101       _bf = (bool_grid_fct_type *) 0;
00102     }
00103     if (_bf_sh) {
00104       delete _bf_sh;
00105       _bf_sh = (bool_grid_fct_type *) 0;
00106     }
00107     base::finish();
00108     for (register int n=0; n<NGFM(); n++) 
00109       GFM(n).finish();
00110   }
00111   virtual void SetupData() {
00112     base::SetupGridHierarchy();
00113     for (register int n=0; n<NGFM(); n++) {
00114       GFM(n).SetErrorEstimation(base::ErrorEstimation()); 
00115       GFM(n).SetupData(base::PGH(), base::NGhosts());
00116     }
00117 
00118     int t_sten = 1;
00119     int s_sten = base::NGhosts();
00120     std::sprintf(BFName,"BFlags");
00121     _bf = new bool_grid_fct_type(BFName, t_sten, s_sten, base::GH(),
00122                                  DAGHCellCentered, DAGHNoComm, DAGHNoBoundary,
00123                                  DAGHNoAdaptBoundary, DAGHNoExternalGhost);
00124     SetCheckpointFlag(BF(), DAGHFalse);
00125     SetTimeAlias(BF(), -1, 1);
00126     if (base::ErrorEstimation()) {
00127       std::sprintf(BFNamesh,"BFlagssh");
00128       _bf_sh = new bool_grid_fct_type(BFNamesh, t_sten, s_sten, base::GH(),
00129                                       DAGHCellCentered, 0, DAGHShadowFactor, 
00130                                       DAGH_All, DAGHNoComm, DAGHNoBoundary,
00131                                       DAGHNoAdaptBoundary, DAGHNoExternalGhost);
00132       SetCheckpointFlag(BFsh(), DAGHFalse);
00133       SetTimeAlias(BFsh(), -1, 1);
00134     }
00135 
00136     base::SetupData();
00137     _GFMRecomposeFunc = new gfm_recompose_functor_type(this, &self::SetRecomposeBndry);
00138     base::U().GF_SetRecomposeFunc(_GFMRecomposeFunc);
00139     if (base::ErrorEstimation())  
00140       base::Ush().GF_SetRecomposeFunc(_GFMRecomposeFunc);
00141     for (register int n=0; n<NGFM(); n++) 
00142       GFM(n).SetGridFunctions(base::_u,base::_u_sh);
00143   }
00144 
00145   virtual void SetBndry(vec_grid_fct_type& u, const int Time, 
00146                         const int Level, double t) {    
00147     base::SetBndry(u,Time,Level,t);
00148     forall (BF(u),Time,Level,c)      
00149       BF(u)(Time,Level,c).equals(false);
00150     end_forall
00151     for (register int n=0; n<NGFM(); n++) {
00152       if (!GFM(n).Stationary() || Time==0) 
00153         GFM(n).SetLevelSet(u,Time,Level,t);     
00154       GFM(n).SetBndry(u,BF(u),Time,Level,t);
00155     }
00156   }
00157 
00158   virtual double IntegrateLevel(vec_grid_fct_type& u, const int Time, const int Level, 
00159                                 double t, double dt, bool DoFixup, double tc, 
00160                                 const int which) {
00161     double cfl = base::IntegrateLevel(u,Time,Level,t,dt,DoFixup,tc,which);
00162 
00163     int TStep = TimeStep(u,Level);     
00164     if (_RecoverExterior) {
00165       forall (u,Time,Level,c)           
00166         if (base::Dim() == 1) {
00167           BeginFastIndex1(u_old, u(Time,Level,c).bbox(),
00168                           u(Time,Level,c).data(), VectorType);
00169           BeginFastIndex1(u_new, u(Time+TStep,Level,c).bbox(),
00170                           u(Time+TStep,Level,c).data(), VectorType);
00171           BeginFastIndex1(bf, BF(u)(Time,Level,c).bbox(),
00172                           BF(u)(Time,Level,c).data(), bool);
00173           if (_RecoverValue) 
00174             for_1 (n, BF(u).interiorbbox(Time,Level,c), BF(u)(Time,Level,c).bbox().stepsize())
00175               if (FastIndex1(bf,n)) 
00176                 FastIndex1(u_new,n) = _RecoveryValue;
00177             end_for
00178           else
00179             for_1 (n, BF(u).interiorbbox(Time,Level,c), BF(u)(Time,Level,c).bbox().stepsize())
00180               if (FastIndex1(bf,n)) 
00181                 FastIndex1(u_new,n) = FastIndex1(u_old,n);
00182             end_for
00183           EndFastIndex1(bf);
00184           EndFastIndex1(u_new);
00185           EndFastIndex1(u_old);
00186         }      
00187         else if (base::Dim() == 2) {
00188           BeginFastIndex2(u_old, u(Time,Level,c).bbox(),
00189                           u(Time,Level,c).data(), VectorType);
00190           BeginFastIndex2(u_new, u(Time+TStep,Level,c).bbox(),
00191                           u(Time+TStep,Level,c).data(), VectorType);
00192           BeginFastIndex2(bf, BF(u)(Time,Level,c).bbox(),
00193                           BF(u)(Time,Level,c).data(), bool);
00194           if (_RecoverValue) 
00195             for_2 (n, m, BF(u).interiorbbox(Time,Level,c), BF(u)(Time,Level,c).bbox().stepsize())
00196               if (FastIndex2(bf,n,m)) 
00197                 FastIndex2(u_new,n,m) = _RecoveryValue;
00198             end_for
00199           else
00200             for_2 (n, m, BF(u).interiorbbox(Time,Level,c), BF(u)(Time,Level,c).bbox().stepsize())
00201               if (FastIndex2(bf,n,m)) 
00202                 FastIndex2(u_new,n,m) = FastIndex2(u_old,n,m);
00203             end_for
00204           EndFastIndex2(bf);
00205           EndFastIndex2(u_new);
00206           EndFastIndex2(u_old);
00207         }
00208         else if (base::Dim() == 3) {
00209           BeginFastIndex3(u_old, u(Time,Level,c).bbox(),
00210                           u(Time,Level,c).data(), VectorType);
00211           BeginFastIndex3(u_new, u(Time+TStep,Level,c).bbox(),
00212                           u(Time+TStep,Level,c).data(), VectorType);
00213           BeginFastIndex3(bf, BF(u)(Time,Level,c).bbox(),
00214                           BF(u)(Time,Level,c).data(), bool);
00215           if (_RecoverValue) 
00216             for_3 (n, m, l, BF(u).interiorbbox(Time,Level,c), BF(u)(Time,Level,c).bbox().stepsize())
00217               if (FastIndex3(bf,n,m,l)) 
00218                 FastIndex3(u_new,n,m,l) = _RecoveryValue;
00219             end_for
00220           else
00221             for_3 (n, m, l, BF(u).interiorbbox(Time,Level,c), BF(u)(Time,Level,c).bbox().stepsize())
00222               if (FastIndex3(bf,n,m,l)) 
00223                 FastIndex3(u_new,n,m,l) = FastIndex3(u_old,n,m,l);
00224             end_for
00225           EndFastIndex3(bf);
00226           EndFastIndex3(u_new);
00227           EndFastIndex3(u_old);
00228         }
00229       end_forall
00230     }
00231 
00232     // Revert GFM boundary cells back to values from interior.
00233     // This allows for moving boundaries and enables the creation
00234     // of new refinement regions along the GFM boundary just using
00235     // interior cell values.
00236     // Swapping time levels also for BF(u) ensures correct flags
00237     // also in (BF(u(Time+TStep)).
00238     for (register int n=0; n<NGFM(); n++) 
00239       if (Level<MaxLevel(base::GH()))
00240         GFM(n).SetBndry(u,BF(u),Time,Level,t,true);
00241 
00242     SwapTimeLevels(u,Level,Time,Time+TStep);
00243     SwapTimeLevels(BF(u),Level,Time,Time+TStep);
00244     for (register int n=0; n<NGFM(); n++) 
00245       if (Level<MaxLevel(base::GH()) || !GFM(n).Stationary())
00246         GFM(n).SetBndry(u,BF(u),Time,Level,t,true);
00247     SwapTimeLevels(u,Level,Time,Time+TStep);
00248     SwapTimeLevels(BF(u),Level,Time,Time+TStep);
00249 
00250     return cfl;
00251   }
00252 
00253   virtual void RecomposeGridHierarchy(const int Time, const int Level, 
00254                                       bool ShadowAllowed, bool DoFixup,
00255                                       bool RecomposeBaseLev, bool RecomposeHighLev) {
00256 
00257     // Recomposition of GridFunction requires LevelTransfer and a BoundaryConditions functor.
00258     // This is currently not available for Levelset functions.
00259     _MaxRecomposeLevel = DAGHNull;
00260     for (register int n=0; n<NGFM(); n++)
00261       GFM(n).SetMaxRecomposeLevel(_MaxRecomposeLevel);
00262     BF().GF_SetMaxRecomposeLevel(_MaxRecomposeLevel);
00263     if (base::ErrorEstimation())  
00264       BFsh().GF_SetMaxRecomposeLevel(_MaxRecomposeLevel);
00265     base::RecomposeGridHierarchy(Time,Level,ShadowAllowed,DoFixup,
00266                                  RecomposeBaseLev,RecomposeHighLev);
00267   }
00268 
00269   virtual void SetRecomposeBndry(vec_grid_fct_type& u, const int& Time, 
00270                                  const int& Level, const double& t) {
00271     if (Time>0 && CurrentTime(base::GH(),Level)==Time &&
00272         Level>_MaxRecomposeLevel) {
00273       forall (BF(u),Time,Level,c)      
00274         BF(u)(Time,Level,c).equals(false);
00275       end_forall
00276       for (register int n=0; n<NGFM(); n++) {
00277         GFM(n).SetLevelSet(u,Time,Level,t);     
00278         GFM(n).SetBndry(u,BF(u),Time,Level,t);
00279       }
00280     }
00281   }
00282 
00283   virtual void Initialize_(const double& dt_start) { 
00284     int* v = new int[NGFM()];
00285     for (register int n=0; n<NGFM(); n++) {
00286       v[n] = GFM(n).Boundary().GetVerbose();
00287       GFM(n).Boundary().SetVerbose(0);
00288     }
00289     base::Initialize_(dt_start);
00290     for (register int n=0; n<NGFM(); n++) 
00291       GFM(n).Boundary().SetVerbose(v[n]);
00292     delete [] v;
00293   }
00294 
00295   virtual void Output() {
00296     if (CurrentTime(base::GH(),0) == base::LastOutputTime()) return;
00297     if (base::_FileOutput) {
00298       START_WATCH 
00299         for (register int n=0; n<NGFM(); n++) 
00300           if (GFM(n).LevelSet().PlotPhi()) {
00301             double t = GetPhysicalTime(base::U(),CurrentTime(base::GH(),0),0);
00302             char name[DAGHBktGFNameWidth+16];
00303             std::sprintf(name,"phi_%d",n+1); 
00304             for (int l=0; l<=FineLevel(base::GH()); l++) {
00305               int Time = CurrentTime(base::GH(),l); 
00306               base::FileOutput_().WritePlain(GFM(n).Phi(), name, Time, l, t);  
00307             }
00308           }
00309       END_WATCH_OUPUT  
00310     }
00311     base::Output();
00312   }
00313 
00314   virtual bool Restart_(const char* CheckpointFile) {
00315     bool ret = base::Restart_(CheckpointFile);
00316     if (ret) {
00317       START_WATCH
00318         for (register int l=0; l<=FineLevel(base::GH()); l++) {
00319           int Time = CurrentTime(base::GH(),l);
00320           for (register int n=0; n<NGFM(); n++) {
00321             SetPhysicalTime(GFM(n).Phi(),Time,l,base::t[l]);
00322             if (base::ErrorEstimation()) {
00323               SetPhysicalTime(GFM(n).Phish(),Time,l,base::t[l]);
00324             }
00325           }
00326           _MaxRecomposeLevel = DAGHNull;
00327           SetRecomposeBndry(base::U(),Time,l,base::t[l]);
00328           if (base::ErrorEstimation()) 
00329             SetRecomposeBndry(base::Ush(),Time,l,base::t[l]);
00330         }
00331       END_WATCH_INITIALIZATION
00332     }
00333     return ret;
00334   }
00335 
00336   virtual void Checkpointing_(const char* CheckpointFile) {
00337     for (register int n=0; n<NGFM(); n++) {
00338       SetCheckpointFlag(GFM(n).Phi(), DAGHFalse);
00339       if (base::ErrorEstimation()) 
00340         SetCheckpointFlag(GFM(n).Phish(), DAGHFalse);
00341     }
00342 
00343     base::Checkpointing_(CheckpointFile);   
00344 
00345     for (register int n=0; n<NGFM(); n++) {
00346       SetCheckpointFlag(GFM(n).Phi(), DAGHTrue);
00347       if (base::ErrorEstimation()) 
00348         SetCheckpointFlag(GFM(n).Phish(), DAGHTrue);
00349     }
00350   }
00351   virtual void Restart_(std::stringstream& CheckpointStr) {
00352     base::Restart_(CheckpointStr);
00353   }
00354     
00355   virtual void Checkpointing_(std::stringstream& CheckpointStr) {
00356     base::Checkpointing_(CheckpointStr); 
00357   }
00358 
00359 
00360 
00361   void AddGFM(gfm_type* gfm) {
00362     if (gfm) {
00363       gfm_type** _GFM_new = new gfm_type*[_nGFM+1];
00364       if (_GFM) {
00365         for (register int n=0; n<_nGFM; n++)
00366           _GFM_new[n] = _GFM[n];
00367         delete [] _GFM;
00368       }
00369       _GFM = _GFM_new;
00370       _GFM[_nGFM] = gfm;
00371       _nGFM++;
00372     }
00373   }
00374 
00375   void EliminateGFM(gfm_type* gfm) {
00376     if (gfm) {
00377       int c=0;
00378       for (register int n=0; n<_nGFM; n++)
00379         if (_GFM[n] == gfm) c++;
00380       if (c>0) {
00381         gfm_type** _GFM_new = (gfm_type**) 0;
00382         if (_nGFM-c>0) {
00383           _GFM_new = new gfm_type*[_nGFM-c];
00384           for (register int n=0, m=0; n<_nGFM; n++)
00385             if (_GFM[n] != gfm) 
00386               _GFM_new[m++] = _GFM[n];
00387         }
00388         delete [] _GFM;
00389         _GFM = _GFM_new;
00390         _nGFM = _nGFM-c;
00391       }
00392     }
00393   }
00394 
00395   void DeleteGFM(gfm_type* gfm) {
00396     EliminateGFM(gfm);
00397     if (gfm) {
00398       delete gfm->GetBoundaryP();
00399       delete gfm->GetLevelSetP();
00400       delete gfm;
00401     }
00402   }
00403 
00404   inline bool_grid_fct_type* BFP() { return _bf; }
00405   inline bool_grid_fct_type& BF() { return *_bf; }
00406   inline const bool_grid_fct_type& BF() const { return *_bf; }
00407   inline bool_grid_fct_type& BFsh() { return *_bf_sh; }
00408   inline const bool_grid_fct_type& BFsh() const{ return *_bf_sh; }
00409   inline bool_grid_fct_type& BF(vec_grid_fct_type& u) 
00410   { return (&u==base::_u_sh ? *_bf_sh : *_bf); }
00411   inline const bool_grid_fct_type& BF(vec_grid_fct_type& u) const 
00412   { return (&u==base::_u_sh ? *_bf_sh : *_bf); }
00413 
00414   inline gfm_type* GFMP(const int n) { assert (n>=0 && n<_nGFM); return _GFM[n]; }
00415   inline gfm_type& GFM(const int n) { assert (n>=0 && n<_nGFM); return *(_GFM[n]); }
00416   inline const gfm_type& GFM(const int n) const { assert (n>=0 && n<_nGFM); return *(_GFM[n]); }
00417 
00418   inline const int& NGFM() const { return _nGFM; }
00419 
00420   inline void SetRecoveryValue(const VectorType& rec) { _RecoveryValue=rec; }
00421   inline VectorType& RecoveryValue() const { return _RecoveryValue; }
00422 
00423 protected:
00424   bool_grid_fct_type *_bf, *_bf_sh;
00425   gfm_recompose_functor_type* _GFMRecomposeFunc;
00426   gfm_type** _GFM;
00427   int _nGFM;
00428   int _MaxRecomposeLevel, _RecoverExterior, _RecoverValue;
00429   char BFName[DAGHBktGFNameWidth], BFNamesh[DAGHBktGFNameWidth];
00430   VectorType _RecoveryValue;
00431 };
00432 
00433 
00434 #endif
Generated on Thu Jun 30 02:19:12 2016 for AMROC Fluid-solver Framework by 1.6.3