AMROC - Source code documentation

00001 #ifndef AMROC_CLP_INTEGRATORBASE_H
00002 #define AMROC_CLP_INTEGRATORBASE_H
00003 
00011 #include "Integrator.h"
00012 #include "ClpIntegratorInterface.h"
00013 #include "Timing.h"
00014 #include <float.h>
00015 
00016 #define FACTOR 1.0e5
00017 
00018 #ifndef ClpIntegratorBaseName
00019 #define ClpIntegratorBase(dim)      name2(ClpIntegratorBase,dim)
00020 #define ClpIntegratorBaseName
00021 #endif
00022 
00031 template <class VectorType, class AuxVectorType>
00032 class ClpIntegratorBase(DIM) : public Integrator(DIM)<VectorType>,
00033     public ClpIntegratorInterface(DIM)<VectorType, AuxVectorType> {
00034   typedef typename VectorType::InternalDataType DataType;
00035   typedef GridData(DIM)<VectorType> vec_grid_data_type;
00036   typedef GridFunction(DIM)<VectorType> vec_grid_fct_type;  
00037   typedef GridData(DIM_1)<VectorType>    ld_vec_grid_data_type;
00038   typedef GridData(DIM_1)<AuxVectorType> ld_aux_grid_data_type;
00039 
00040   typedef TSSrcIntegrator(DIM)<VectorType> src_integrator_type;
00041   typedef ClpIntegratorBase(DIM)<VectorType,AuxVectorType> self;
00042   typedef ClpTSSrcIntegratorBase(DIM)<VectorType,AuxVectorType> src_type;
00043   typedef Integrator(DIM)<VectorType> integrator_base_type;
00044   typedef ClpSetAux(DIM)<VectorType,AuxVectorType> set_aux_type;
00045 
00046 public:
00047   ClpIntegratorBase(DIM)(const int equ, const int wv, const int gh) : 
00048     integrator_base_type(gh), _EqUsed(equ), _Waves(wv) {
00049     normal_flux_func = f_normalflux;  
00050     limiter = new int[NWaves()];
00051     mthlim = new int[NWaves()];
00052     for (int i=0; i<NWaves(); i++) {
00053       limiter[i] = -1; mthlim[i] = -1;
00054     }
00055     _name = "";
00056     _set_aux = new set_aux_type(*this);
00057     for (int d=0; d<7; d++) method[d] = 0;
00058   }
00059 
00060   ~ClpIntegratorBase(DIM)() {
00061     if (limiter) delete [] limiter;
00062     if (mthlim) delete [] mthlim;
00063     delete _set_aux;
00064   }
00065 
00066   virtual void register_at(ControlDevice& Ctrl, const string& prefix) {
00067     ControlDevice LocCtrl = Ctrl.getSubDevice(prefix+"ClawpackIntegrator");
00068     _name = prefix+"ClawpackIntegrator";
00069     char MethodName[16];
00070     register int d;
00071     for (d=0; d<7; d++) {
00072       sprintf(MethodName,"Method(%d)",d+1);
00073       RegisterAt(LocCtrl,MethodName,method[d]);
00074     }
00075     char LimiterName[16];
00076     for (d=0; d<NWaves(); d++) {
00077       sprintf(LimiterName,"Limiter(%d)",d+1);
00078       RegisterAt(LocCtrl,LimiterName,limiter[d]);
00079     }
00080     if (_SrcIntegrator)
00081       SrcIntegrator().register_at(LocCtrl,prefix);
00082   }
00083   virtual void register_at(ControlDevice& Ctrl) {
00084     register_at(Ctrl, "");
00085   }
00086    
00087   virtual void SetupData() {
00088     for (int i=0; i<NWaves(); i++) 
00089       if (limiter[i] > 0) mthlim[i] = limiter[i];
00090       else mthlim[i] = limiter[0];
00091     if (NAux()!=0 && NAux()!=AuxVectorType::Length())
00092       method[6] = 0;
00093     if (!_SrcIntegrator)
00094       method[4] = 0;    
00095     SetSrcTreatment(method[4]);
00096     SetMaxIntegratorPasses(NMaxPass());
00097     _abort = (NCheck()>=10 ? DAGHTrue : DAGHFalse);
00098     integrator_base_type::SetupData();
00099   }
00100 
00101   virtual void finish() {
00102     if (limiter) {
00103       delete [] limiter;
00104       limiter = (int *) 0;
00105     }
00106     if (mthlim) {
00107       delete [] mthlim;
00108       mthlim = (int *) 0;
00109     }
00110     integrator_base_type::finish();
00111   } 
00112 
00113   void CalculateRP(ld_vec_grid_data_type &LeftState, BBox& bbl, 
00114                    double tl, double dtl, DCoords& dxl,
00115                    ld_vec_grid_data_type &RightState, BBox& bbr, 
00116                    double tr, double dtr, DCoords& dxr,
00117                    ld_vec_grid_data_type &NewState, 
00118                    const int& direction) {
00119     
00120     //------------------------------------------------------------------                      
00121     // Calculate AuxArray and source term on left and right state.
00122     //------------------------------------------------------------------                      
00123     vec_grid_data_type left(bbl, LeftState.data());
00124     vec_grid_data_type right(bbr, RightState.data());
00125     AllocError::SetTexts("ClpIntegrator","calculate_Riemann_Problem(): aux arrays");
00126     DataType* auxldata = new DataType[LeftState.bbox().size()*NAux()];
00127     DataType* auxrdata = new DataType[RightState.bbox().size()*NAux()];
00128 
00129     // Set aux arrays.
00130     if (NAux()>0) {
00131       AuxArray().Set(auxldata,left,tl,dtl,dxl);
00132       AuxArray().Set(auxrdata,right,tr,dtr,dxr);
00133     }      
00134 
00135     // Apply the source term, if necessary.
00136     START_INTERMEDIATE_WATCH
00137     START_WATCH
00138       if (NSrc() > 0) {
00139         // Integrate source term on left and right data to the same time level.
00140         int srcbnd = 1;
00141         if (tl < tr) {
00142           ((src_type &)SrcIntegrator()).SetAuxArray(auxldata,NAux());
00143           ((src_type &)SrcIntegrator()).CalculateSrcOnGrid(left,tl,tr-tl,dxl,srcbnd);
00144           tl += tr-tl;
00145         }
00146         if (tl > tr) {
00147           ((src_type &)SrcIntegrator()).SetAuxArray(auxrdata,NAux());        
00148           ((src_type &)SrcIntegrator()).CalculateSrcOnGrid(right,tr,tl-tr,dxr,srcbnd);
00149           tr += tl-tr;
00150         }
00151 
00152         if (NSrc() > 1) {
00153           double dt = Min(dtl,dtr);
00154           if (NSrc()==2) {
00155             ((src_type &)SrcIntegrator()).SetAuxArray(auxldata,NAux());
00156             ((src_type &)SrcIntegrator()).CalculateSrcOnGrid(left,tl,dt/2.0,dxl,srcbnd);
00157             ((src_type &)SrcIntegrator()).SetAuxArray(auxrdata,NAux());      
00158             ((src_type &)SrcIntegrator()).CalculateSrcOnGrid(right,tr,dt/2.0,dxr,srcbnd);
00159           }
00160           if (NSrc()==3) {
00161             ((src_type &)SrcIntegrator()).SetAuxArray(auxldata,NAux());      
00162             ((src_type &)SrcIntegrator()).CalculateSrcOnGrid(left,tl,dt,dxl,srcbnd);
00163             ((src_type &)SrcIntegrator()).SetAuxArray(auxrdata,NAux());      
00164             ((src_type &)SrcIntegrator()).CalculateSrcOnGrid(right,tr,dt,dxr,srcbnd);
00165           }
00166           if (NAux()>0) {
00167             AuxArray().Set(auxldata,left,tl,dt,dxl);
00168             AuxArray().Set(auxrdata,right,tr,dt,dxr);
00169           }      
00170         }         
00171       }
00172     END_WATCH_SOURCE_INTEGRATION
00173     END_INTERMEDIATE_WATCH
00174     VectorType* data_dummy;
00175     left.deallocate(data_dummy);
00176     right.deallocate(data_dummy);
00177     ld_aux_grid_data_type auxl(LeftState.bbox(), (AuxVectorType*) auxldata);
00178     ld_aux_grid_data_type auxr(RightState.bbox(), (AuxVectorType*) auxrdata);
00179 
00180     // Solve the Riemann problems. 
00181     ComputeRP(LeftState, auxl, RightState, auxr, NewState, direction);
00182 
00183     AuxVectorType* aux_data_dummy;
00184     auxl.deallocate(aux_data_dummy);
00185     auxr.deallocate(aux_data_dummy);
00186     delete [] auxldata; delete [] auxrdata;
00187   }
00188 
00189   virtual double CalculateGrid(vec_grid_data_type& NewStateVec, 
00190                                vec_grid_data_type& OldStateVec,
00191                                vec_grid_data_type* Flux[],
00192                                double t, double dt, DCoords& dx,
00193                                const int& mpass) {
00194     double cfl = 0.0;
00195 #ifdef DEBUG_PRINT_INTEGRATOR
00196     ( comm_service::log() << "on: " << OldStateVec.bbox() ).flush();
00197 #endif
00198 
00199     if (NMaxPass()==0 || mpass==1) 
00200       NewStateVec.copy(OldStateVec);
00201 
00202     if (dt >= DBL_EPSILON*FACTOR) {
00203       AllocError::SetTexts("ClpIntegrator","calculate_time_step(): aux");
00204       DataType* aux = new DataType[OldStateVec.bbox().size()*NAux()];     
00205       if (NAux() > 0) AuxArray().Set(aux,NewStateVec,t,dt,dx);
00206       if (NSrc() > 0) 
00207         ((src_type &)SrcIntegrator()).SetAuxArray(aux,NAux());  
00208 
00209 #ifdef DEBUG_AMR
00210       if (NCheck()%10>1 && (NSrc()==2 || NSrc()==3))
00211         CheckGrid(NewStateVec, NewStateVec.bbox(), "CalculateGrid::before SrcInt");
00212 #endif
00213 
00214       // The source term has to be applied also on ghostcells, otherwise boundary conditions
00215       // would not be consistent. See CLAWPACK-Manual Note #5 "Source Terms".
00216       START_INTERMEDIATE_WATCH
00217       START_WATCH
00218         if (NSrc()>1 && (NMaxPass()==0 || mpass==1)) {
00219           int srcbnd = 1;
00220           if (NSrc()==2)
00221             SrcIntegrator().CalculateSrcOnGrid(NewStateVec, t, dt/2.0, dx, srcbnd);
00222           if (NSrc()==3) 
00223             SrcIntegrator().CalculateSrcOnGrid(NewStateVec, t, dt, dx, srcbnd);
00224           if (NAux()>0) AuxArray().Set(aux,NewStateVec,t,dt,dx);
00225         }
00226       END_WATCH_SOURCE_INTEGRATION
00227       END_INTERMEDIATE_WATCH
00228 
00229       if (NCheck() && (NMaxPass()==0 || mpass==1))
00230         CheckGrid(NewStateVec, NewStateVec.bbox(), "CalculateGrid::before f_step");
00231 
00232       // Update NewStateVec. 
00233       f_init_rcommon(t, dt, AA(DIM,dx()), mpass); 
00234       cfl = ComputeGrid(NewStateVec, Flux, aux, dt, dx);
00235 
00236 #ifdef DEBUG_AMR
00237       if (NCheck()%10>0 && (NMaxPass()==0 || mpass==NMaxPass()))
00238         CheckGrid(NewStateVec, NewStateVec.bbox(), "CalculateGrid::after f_step");
00239 #endif
00240 
00241       // The source term is not applied to ghostcells, because they are overwritten anyway by
00242       // the application of boundary conditions in the next step.
00243       START_INTERMEDIATE_WATCH
00244       START_WATCH
00245         if (NSrc()>0 &&  NSrc()<3 && (NMaxPass()==0 || mpass==NMaxPass())) {
00246           int srcbnd = 0;
00247           if (NAux()>0) AuxArray().Set(aux,NewStateVec,t,dt,dx);
00248           if (NSrc()==2)
00249             SrcIntegrator().CalculateSrcOnGrid(NewStateVec, t+dt/2.0,dt/2.0, dx, srcbnd);
00250           if (NSrc()==1)
00251             SrcIntegrator().CalculateSrcOnGrid(NewStateVec, t, dt, dx, srcbnd);
00252         }
00253       END_WATCH_SOURCE_INTEGRATION
00254       END_INTERMEDIATE_WATCH
00255 
00256 #ifdef DEBUG_AMR
00257       if (NCheck()%10>2 && (NSrc()==2 || NSrc()==1) && (mpass==0 || mpass==NMaxPass()))
00258         CheckGrid(NewStateVec, NewStateVec.bbox(), "CalculateGrid::after SrcInt");
00259 #endif
00260 
00261       delete [] aux;
00262     }
00263     else {
00264       VectorType zero(0.0);
00265       for (int d=0; d<DIM; d++)
00266         Flux[d]->equals(zero);
00267 #ifdef DEBUG_PRINT_INTEGRATOR
00268       ( comm_service::log() <<  "dt < eps*1.0e5  " ).flush();
00269 #endif
00270     }
00271 
00272 #ifdef DEBUG_PRINT_INTEGRATOR
00273       ( comm_service::log() << "  CFL = " << cfl <<
00274         "   t = " << t << "   dt = " << dt << "\n" ).flush();
00275 #endif
00276     
00277     return cfl;
00278   }
00279 
00280   virtual int ControlGrid(vec_grid_data_type& StateVec, const BBox& where) {
00281     int result;
00282     f_check(FA(DIM,StateVec),BOUNDING_BOX(where),
00283             NEquations(),result);
00284     return result; 
00285   }  
00286 
00287   virtual int NMethodOrder() { return method[1]; }
00288   inline const int& NCheck() const { return method[3]; }
00289   inline int NSrc() const { return (method[4]%10); }
00290   inline int NMaxPass() const { return (method[2]<-1 && DIM>1 ? DIM : 0); }
00291   inline const int& NAux() const { return method[6]; }
00292   inline const int& NEqUsed() const { return _EqUsed; }
00293   inline const int& NWaves() const { return _Waves; }  
00294   inline set_aux_type& AuxArray() { return *_set_aux; }
00295   inline const set_aux_type& AuxArray() const { return *_set_aux; }
00296 
00297   ostream& debug(ostream& out, const self& P) {
00298     out << "ClpIntegrator registered at " << P._name << "\n";
00299     out << "   Methods:";
00300     for (int i=0; i<7; i++)
00301       out << " " << P.method[i];
00302     out << "   Limiter:";
00303     for (int j=0; j<P.NWaves(); j++)
00304       out << " " << P.limiter[j];
00305     out << "\n";
00306     return out;
00307   }
00308 
00309 protected:
00310   int _EqUsed, _Waves;
00311   normal_func_type normal_flux_func;
00312   int method[7];
00313   int *limiter, *mthlim;
00314   string _name;
00315   set_aux_type* _set_aux;
00316 };
00317 
00318 
00319 
00320 #endif
ClpIntegratorBase.h

last update: 06/01/04
