amroc/weno/WENOIntegrator.h

Go to the documentation of this file.
00001 // -*- C++ -*-
00002 
00003 // Copyright (C) 2003-2007 California Institute of Technology
00004 // Ralf Deiterding, Carlos Pantano, David Hill 
00005 
00006 #ifndef AMROC_WENO_INTEGRATOR_H
00007 #define AMROC_WENO_INTEGRATOR_H
00008 
00016 #include "Integrator.h"
00017 #include "CLESLog.h"
00018 
00019 #include <string>
00020 
00021 #ifndef f_init_common
00022 #define f_init_common FORTRAN_NAME(combl, COMBL)
00023 #endif
00024 
00025 #define f_init_weno FORTRAN_NAME(initweno, INITWENO)
00026 #define f_finish_weno FORTRAN_NAME(finishweno, FINISHWENO)
00027 
00028 extern "C" {
00029   void f_init_weno(const INTEGER& dim, const INTEGER& meqn, const INTEGER& nvars, 
00030                    const INTEGER& mghosts, 
00031                    const INTEGER& order, const INTEGER& optimized,
00032                    const INTEGER& use_carbfix, const INTEGER& method, 
00033                    const INTEGER& viscous, const INTEGER& les, 
00034                    const INTEGER& usrc, 
00035                    const INTEGER& noTimeRefine, const DOUBLE& filter_alpha);
00036   void f_finish_weno();
00037   void f_init_common();
00038 }
00039 
00040 
00049 template <class VectorType, int dim>
00050 class WENOIntegrator : public Integrator<VectorType,dim> {
00051   typedef typename VectorType::InternalDataType DataType;
00052   typedef Integrator<VectorType,dim> base;
00053 
00054 public:
00055   typedef typename base::vec_grid_fct_type vec_grid_fct_type;  
00056   typedef typename base::vec_grid_data_type vec_grid_data_type;
00057   typedef generic_fortran_func generic_func_type;
00058 
00059   typedef int DCFlagDataType;
00060   typedef GridData<DCFlagDataType,dim> flag_grid_data_type;  
00061 
00062   typedef void (*check_1_func_type) ( FI(1,VectorType), BI, const INTEGER& meqn, 
00063                                       const INTEGER& mout, INTEGER& result );
00064   typedef void (*check_2_func_type) ( FI(2,VectorType), BI, const INTEGER& meqn, 
00065                                       const INTEGER& mout, INTEGER& result );
00066   typedef void (*check_3_func_type) ( FI(3,VectorType), BI, const INTEGER& meqn, 
00067                                       const INTEGER& mout, INTEGER& result );
00068 
00069   typedef void (*step_func_type) (const INTEGER& rk, VectorType ux[], 
00070                                   VectorType uxold[], const DOUBLE& dt , 
00071                                   DOUBLE& cfl, VectorType fxi[], 
00072                                   INTEGER dcflag[], INTEGER& iFilter);
00073   typedef void (*bnds_func_type) (INTEGER ix[], DOUBLE lbc[], DOUBLE ubc[], 
00074                                   DOUBLE dx[], const DOUBLE* bnd, 
00075                                   const INTEGER& mb, const INTEGER per[]);
00076 
00077   WENOIntegrator(generic_func_type step, generic_func_type bnds, generic_func_type check) : 
00078     base(), f_step(step), f_bnds(bnds), f_chk(check),
00079     _order(2), _optimized(0), _use_carbfix(1), _method(0),
00080     _visc(0), _les(0),_usrc(0), _check(0), _DCFlagAdaptBndry(1),
00081     _noTimeRefine(0), _FilterStep(-1), _FilterStrength(0.0)
00082   { 
00083     _name = ""; 
00084     FluxData = (VectorType*) 0;
00085   }
00086 
00087   ~WENOIntegrator() {
00088     if (FluxData) {
00089       delete [] FluxData;
00090       FluxData = (VectorType*) 0;
00091     }  
00092   }
00093 
00094   virtual double PassTimeStepFraction(int mpass) { 
00095     switch ( mpass ) {
00096     case 1:
00097       return 0.0;
00098     case 2:
00099       return 1.0;
00100     case 3:
00101       return 0.5;
00102     default:
00103       return 0.0;
00104     }
00105   }
00106 
00107   virtual void register_at(ControlDevice& Ctrl, const std::string& prefix) {
00108     ControlDevice LocCtrl = Ctrl.getSubDevice(prefix+"WENOIntegrator");
00109     _name = prefix+"WENOIntegrator";
00110     RegisterAt(LocCtrl,"Order",_order);
00111     RegisterAt(LocCtrl,"Optimized",_optimized);
00112     RegisterAt(LocCtrl,"CarbuncleFix",_use_carbfix);
00113     RegisterAt(LocCtrl,"Method",_method);
00114     RegisterAt(LocCtrl,"UseViscous",_visc);
00115     RegisterAt(LocCtrl,"UseLES",_les);
00116     RegisterAt(LocCtrl,"UseSource",_usrc);
00117     RegisterAt(LocCtrl,"Check",_check);
00118     RegisterAt(LocCtrl,"DCFlagAdaptBndry",_DCFlagAdaptBndry);
00119     
00120     RegisterAt(LocCtrl,"FilterStep",_FilterStep);
00121     RegisterAt(LocCtrl,"FilterStrength",_FilterStrength);
00122 
00123     log.register_at(LocCtrl, ""); 
00124   }
00125   virtual void register_at(ControlDevice& Ctrl) {
00126     register_at(Ctrl, "");
00127   }
00128    
00129   virtual void SetupData(GridHierarchy* gh, const int& ghosts) {  
00130     base::SetupData(gh,ghosts);
00131     base::SetMaxIntegratorPasses(NMaxPass());
00132     base::_abort = (NCheck()>=10 ? DAGHTrue : DAGHFalse);
00133     int _dim = dim;
00134 
00135     log.initialize();
00136 
00137     f_init_weno(_dim, base::NEquations(), NVARS, base::NGhosts(),
00138                 _order, _optimized, _use_carbfix, _method, _visc, _les,
00139                 _usrc, _noTimeRefine, _FilterStrength);
00140     f_init_common();
00141   }
00142 
00143   virtual void finish() { f_finish_weno(); log.finalize(); } 
00144 
00145   virtual double CalculateGrid(vec_grid_data_type& NewStateVec, 
00146                                vec_grid_data_type& OldStateVec,
00147                                vec_grid_data_type* Flux[],
00148                                const int& level, const double& t, 
00149                                const double& dt, const int& mpass) {
00150     double cfl = 0.0;
00151 #ifdef DEBUG_PRINT_INTEGRATOR
00152     ( comm_service::log() << "on: " << OldStateVec.bbox() ).flush();
00153 #endif
00154 
00155     if (mpass==1 || mpass==0) NewStateVec.copy(OldStateVec);
00156 
00157     int mx[dim], d, per[dim];
00158     Coords ex = NewStateVec.extents();
00159     for (d=0; d<dim; d++) { 
00160       mx[d] = ex(d)-2*base::NGhosts();
00161       per[d] = base::GH().periodicboundary(d);
00162     }
00163 
00164     flag_grid_data_type** DCFlag;
00165     DCFlagDataType* DCFlagData;
00166     AllocDCFlag(NewStateVec.bbox(), DCFlag, DCFlagData);
00167     if (mpass==1 || mpass==0) SetDCFlag(NewStateVec.bbox(), DCFlag, level);
00168 
00169     DCoords dx = base::GH().worldStep(NewStateVec.stepsize());
00170     DCoords lbc = base::GH().worldCoords(NewStateVec.lower(), NewStateVec.stepsize());
00171     DCoords ubc = base::GH().worldCoords(NewStateVec.upper(), NewStateVec.stepsize());
00172 
00173     ((bnds_func_type) f_bnds)(mx,lbc(),ubc(),dx(),base::GH().wholebndry(),
00174                               base::GH().nbndry(),per);
00175     
00176     if (NCheck()>0 && (mpass==1 || mpass==0))
00177       base::CheckGrid(NewStateVec, level, NewStateVec.bbox(), t, "CalculateGrid::before f_step");
00178 
00179     int iFilter = 0;
00180     int Time = CurrentTime(base::GH(), 0 );
00181     if ( _FilterStep>0 && Time%(_FilterStep*StepSize(base::GH(),0)) == 0 ) iFilter = 1;
00182 
00183     // The ghostcells must remain untouched during integration!
00184     ((step_func_type) f_step)(mpass,NewStateVec.data(),OldStateVec.data(), 
00185                               dt,cfl,Flux[0]->data(),DCFlag[0]->data(),iFilter);
00186 
00187     if (NCheck()>=0) 
00188       // The negation captures cfl=NaN!
00189       if (ControlGrid(NewStateVec,level,NewStateVec.bbox(),t,0)==0 || !(cfl>=0) ||
00190           (!(cfl>0) && (mpass==NMaxPass()))) {  
00191         ResetGridFluxes(Flux);
00192         NewStateVec.copy(OldStateVec);
00193         if (cfl>0) cfl = std::max(cfl, 2.0);
00194       }
00195     
00196     DeAllocDCFlag(DCFlag,DCFlagData);
00197 
00198 #ifdef DEBUG_PRINT_INTEGRATOR
00199       ( comm_service::log() << "  CFL = " << cfl <<
00200         "   t = " << t << "   dt = " << dt << "\n" ).flush();
00201 #endif
00202     
00203     return cfl;
00204   } 
00205 
00206   void SetDCFlag(const BBox &databbox, flag_grid_data_type** DCFlag, const int& l) {
00207     int d;
00208     for (d=0; d<dim; d++) DCFlag[d]->equals(0);
00209 
00210     if ( _DCFlagAdaptBndry==0 ) return;
00211 
00212     BBox interiorbbox = grow(databbox, -base::NGhosts());
00213 
00214     // check if need to flag fine side.
00215     if (base::U().prolongbndrylevel(l)!=0 ) {
00216       for (int s=0; s<2*dim; s++) { 
00217         d = s/2;
00218         BBox localbbndry(interiorbbox);
00219         if (s%2 == 0) {
00220           localbbndry.growlower(d,-1); 
00221           localbbndry.upper(d) = localbbndry.lower(d);
00222         }
00223         else {
00224           localbbndry.growupper(d,1); 
00225           localbbndry.lower(d) = localbbndry.upper(d);
00226         }
00227         BBoxList localbndrylist = *base::U().prolongbndrylevel(l);
00228         localbndrylist *= localbbndry;
00229         for (register BBox* bb=localbndrylist.first();bb;
00230              bb=localbndrylist.next()) {
00231           if (s%2 == 0) {
00232             bb->lower(d) += bb->stepsize(d);
00233             bb->upper(d) = bb->lower(d);
00234           } 
00235           DCFlag[d]->equals((-1+2*(s%2)),*bb);
00236         }
00237       }
00238     }
00239 
00240     if ( _DCFlagAdaptBndry==1 ) return;
00241     
00242     // check if we need to flag coarse side
00243     if ( l < FineLevel(base::GH()) ) {
00244       int by = StepSize(base::GH(),l)/StepSize(base::GH(),l+1);
00245 #ifdef _DEBUG_DCFLAG
00246       printf("\nGhost cells %d\n", base::NGhosts());
00247       printf("father(%d) (%d %d) (%d %d) [%d %d]\n",l,
00248              interiorbbox.lower(0), interiorbbox.upper(0), 
00249              interiorbbox.lower(1), interiorbbox.upper(1), 
00250              interiorbbox.stepsize(0), interiorbbox.stepsize(1));
00251 #endif
00252       GridBoxList* gbchild = base::GH().ggbl(l+1);
00253       if ( gbchild )
00254         for (register GridBox *g=gbchild->first();g;g=gbchild->next())
00255           {
00256             //BBox bb = grow(g->gbBBox(), base::NGhosts());
00257             BBox bb(g->gbBBox());
00258             bb.coarsen(by);
00259 #ifdef _DEBUG_DCFLAG
00260             printf("child(%d) (%d %d) (%d %d) [%d %d]\n",l+1,
00261                    bb.lower(0), bb.upper(0), 
00262                    bb.lower(1), bb.upper(1), 
00263                    bb.stepsize(0), bb.stepsize(1));
00264 #endif
00265             
00266             for (int s=0; s<2*dim; s++) { 
00267               d = s/2;
00268               BBox localbbndry(bb);
00269               if (s%2 == 0) {
00270                 localbbndry.growlower(d,-1); 
00271                 localbbndry.upper(d) = localbbndry.lower(d);
00272               }
00273               else {
00274                 localbbndry.growupper(d,1); 
00275                 localbbndry.lower(d) = localbbndry.upper(d);
00276               }
00277 #ifdef _DEBUG_DCFLAG
00278               printf("to intersect[%d] (%d %d) (%d %d) [%d %d]\n", s,
00279                      localbbndry.lower(0), localbbndry.upper(0), 
00280                      localbbndry.lower(1), localbbndry.upper(1), 
00281                      localbbndry.stepsize(0), localbbndry.stepsize(1));
00282 #endif
00283               localbbndry *= interiorbbox;
00284               if ( ! localbbndry.empty() ) {
00285                 if ( s%2 == 0 ) {
00286                   localbbndry.growlower(d,1); 
00287                   localbbndry.upper(d) = localbbndry.lower(d);
00288                 }
00289 #ifdef _DEBUG_DCFLAG
00290                 printf("intersection (%d %d) (%d %d) [%d %d]\n",
00291                        localbbndry.lower(0), localbbndry.upper(0), 
00292                        localbbndry.lower(1), localbbndry.upper(1), 
00293                        localbbndry.stepsize(0), localbbndry.stepsize(1));
00294 #endif
00295                 
00296                 DCFlag[d]->equals((1-2*(s%2)),localbbndry);
00297               }
00298             }
00299           }
00300     }
00301   }
00302     
00303   void AllocDCFlag(const BBox &databbox, flag_grid_data_type**& DCFlag, 
00304                    DCFlagDataType*& DCFlagData) {
00305     BBox fluxbbox = grow(databbox, -base::NGhosts());
00306     fluxbbox.growupper(1);
00307     // Create a CONSECUTIVE array for the flags in ALL directions!
00308     DCFlagData = new DCFlagDataType[base::Dim()*fluxbbox.size()];
00309     DCFlag = new flag_grid_data_type* [base::Dim()]; 
00310     for (register int d=0; d<base::Dim(); d++)
00311       DCFlag[d] = new flag_grid_data_type(fluxbbox,&(DCFlagData[d*fluxbbox.size()])); 
00312   }
00313 
00314   void DeAllocDCFlag(flag_grid_data_type**& DCFlag, DCFlagDataType*& DCFlagData) {
00315     if (DCFlag) {
00316       DCFlagDataType* dummy;
00317       for (register int d=0; d<base::Dim(); d++) 
00318         if (DCFlag[d]) {
00319           // Erease pointer to data
00320           DCFlag[d]->deallocate(dummy);
00321           delete DCFlag[d];
00322         }
00323       delete [] DCFlag;
00324     }
00325     if (DCFlagData) {
00326       delete [] DCFlagData;
00327       DCFlagData = (DCFlagDataType*) 0;
00328     }
00329   }
00330     
00331   virtual void AllocGridFluxes(const BBox &bb, vec_grid_data_type**& Flux) {
00332     // The solver assumes a CONSECUTIVE array for the fluxes in ALL directions!
00333     // We need a member variable to store it between different member function calls.
00334     FluxData = new VectorType[base::Dim()*bb.size()];
00335     Flux = new vec_grid_data_type* [2*base::Dim()]; 
00336     for (register int d=0; d<base::Dim(); d++) {
00337       Flux[2*d] = new vec_grid_data_type(bb,&(FluxData[d*bb.size()])); 
00338       Flux[2*d+1] = Flux[2*d];
00339     }
00340   }
00341 
00342   virtual void DeAllocGridFluxes(vec_grid_data_type**& Flux) {
00343     if (Flux) {
00344       VectorType* dummy;
00345       for (register int d=0; d<base::Dim(); d++) 
00346         if (Flux[2*d]) {
00347           // Erease pointer to data
00348           Flux[2*d]->deallocate(dummy);
00349           delete Flux[2*d];
00350         }
00351       delete [] Flux;
00352     }
00353     if (FluxData) {
00354       delete [] FluxData;
00355       FluxData = (VectorType*) 0;
00356     }
00357   }
00358     
00359   virtual void ResetGridFluxes(vec_grid_data_type**& Flux) {
00360     if (Flux) {
00361       VectorType zero(0.0);
00362       for (register int d=0; d<base::Dim(); d++) 
00363         if (Flux[2*d]) 
00364           Flux[2*d]->equals(zero);
00365     }
00366   }
00367 
00368   virtual int ControlGrid(vec_grid_data_type& StateVec, const int& level, 
00369                           const BBox& where, const double& time, const int verbose) {
00370     int result = 1;
00371     if (!f_chk) return result;
00372     if (dim == 1) 
00373       ((check_1_func_type) f_chk)(FA(1,StateVec),BOUNDING_BOX(where),base::NEquations(),
00374                                   verbose,result);
00375     else if (dim == 2) 
00376       ((check_2_func_type) f_chk)(FA(2,StateVec),BOUNDING_BOX(where),base::NEquations(),
00377                                   verbose,result);
00378     else if (dim == 3) 
00379       ((check_3_func_type) f_chk)(FA(3,StateVec),BOUNDING_BOX(where),base::NEquations(),
00380                                   verbose,result);
00381     return result; 
00382   }  
00383 
00384   inline const int& NCheck() const { return _check; }
00385   virtual int NMethodOrder() const { return _order; }
00386   inline int NMaxPass() const { return 3; }
00387   virtual void SetNoTimeRefine(int flag) { _noTimeRefine = flag; }
00388 
00389   inline void SetCheckFunc(generic_func_type check) { f_chk = check; }
00390   generic_func_type GetCheckFunc() const { return f_chk; }
00391   inline void SetStepFunc(generic_func_type step) { f_step = step; }
00392   generic_func_type GetStepFunc() const { return f_step; }
00393   inline void SetBndsFunc(generic_func_type bnds) { f_bnds = bnds; }
00394   generic_func_type GetBndsFunc() const { return f_bnds; }
00395 
00396 protected:
00397   generic_func_type f_step, f_bnds, f_chk;
00398   std::string _name;
00399   int _order,_optimized, _use_carbfix;
00400   int _method, _visc, _les, _usrc, _check;
00401   int _DCFlagAdaptBndry;
00402   int _noTimeRefine;
00403   VectorType* FluxData;
00404   CLESLog log;
00405   int _FilterStep;
00406   double _FilterStrength;
00407 };
00408 
00409 
00410 
00411 #endif