amroc/rim/RIMIntegrator.h

Go to the documentation of this file.
00001 // -*- C++ -*-
00002 // $Header: /home/proj/vtfng/vtf/amroc/rim/RIMIntegrator.h,v 1.6 2006/05/27 05:59:45 ralf Exp $
00003 
00004 #ifndef AMROC_RIM_INTEGRATOR_H
00005 #define AMROC_RIM_INTEGRATOR_H
00006 
00011 #include "Integrator.h"
00012 
00013 #include <string>
00014 
00015 #define f_init_common FORTRAN_NAME(combl, COMBL)
00016 
00017 extern "C" {
00018   void f_init_common();
00019 }
00020 
00026 template <class VectorType, int dim>
00027 class RIMIntegrator : public Integrator<VectorType,dim> {
00028   typedef typename VectorType::InternalDataType DataType;
00029   typedef Integrator<VectorType,dim> base;
00030 
00031 public:
00032   typedef typename base::vec_grid_fct_type vec_grid_fct_type;  
00033   typedef typename base::vec_grid_data_type vec_grid_data_type;
00034   typedef generic_fortran_func generic_func_type;
00035 
00036   typedef void (*check_1_func_type) ( FI(1,VectorType), BI, const INTEGER& meqn, 
00037                                       const INTEGER& mout, INTEGER& result );
00038   typedef void (*check_2_func_type) ( FI(2,VectorType), BI, const INTEGER& meqn, 
00039                                       const INTEGER& mout, INTEGER& result );
00040   typedef void (*check_3_func_type) ( FI(3,VectorType), BI, const INTEGER& meqn, 
00041                                       const INTEGER& mout, INTEGER& result );
00042 
00043   typedef void (*step_func_type) (const INTEGER& nx, const INTEGER& ny, const INTEGER& nbc, 
00044                                   const INTEGER& neqn, const DOUBLE& dx,
00045                                   const DOUBLE& dy, VectorType ux[], VectorType uxold[], 
00046                                   const DOUBLE& dt , DOUBLE& cfl, VectorType fx[], VectorType fy[]);
00047 
00048   RIMIntegrator(generic_func_type step) : 
00049     base(), f_step(step), f_chk(0),
00050     _order(2), _check(0) { 
00051     _name = ""; 
00052     FluxData = (VectorType*) 0;
00053   }
00054 
00055   RIMIntegrator(generic_func_type step, generic_func_type chk) : 
00056     base(), f_step(step), f_chk(chk),
00057     _order(2), _check(0) { 
00058     _name = ""; 
00059     FluxData = (VectorType*) 0;
00060   }
00061 
00062   ~RIMIntegrator() {
00063     if (FluxData) {
00064       delete [] FluxData;
00065       FluxData = (VectorType*) 0;
00066     }  
00067   }
00068 
00069   virtual void register_at(ControlDevice& Ctrl, const std::string& prefix) {
00070     ControlDevice LocCtrl = Ctrl.getSubDevice(prefix+"RIMIntegrator");
00071     _name = prefix+"RIMIntegrator";
00072     RegisterAt(LocCtrl,"Check",_check);
00073   }
00074   virtual void register_at(ControlDevice& Ctrl) {
00075     register_at(Ctrl, "");
00076   }
00077    
00078   virtual void SetupData(GridHierarchy* gh, const int& ghosts) {  
00079     base::SetupData(gh,ghosts);
00080     base::SetMaxIntegratorPasses(NMaxPass());
00081     base::_abort = (NCheck()>=10 ? DAGHTrue : DAGHFalse);
00082     f_init_common();
00083   }
00084 
00085   virtual double CalculateGrid(vec_grid_data_type& NewStateVec, 
00086                                vec_grid_data_type& OldStateVec,
00087                                vec_grid_data_type* Flux[],
00088                                const int& level, const double& t, 
00089                                const double& dt, const int& mpass) {
00090     double cfl = 0.0;
00091 
00092     int mx[DIM], d, per[DIM];
00093     Coords ex = NewStateVec.extents();
00094     for (d=0; d<DIM; d++) { 
00095       mx[d] = ex(d)-2*base::NGhosts();
00096       per[d] = base::GH().periodicboundary(d);
00097     }
00098     DCoords dx = base::GH().worldStep(NewStateVec.stepsize());
00099     DCoords lbc = base::GH().worldCoords(NewStateVec.lower(), NewStateVec.stepsize());
00100     DCoords ubc = base::GH().worldCoords(NewStateVec.upper(), NewStateVec.stepsize());
00101 
00102     if (NCheck()>0)
00103       base::CheckGrid(OldStateVec, level, OldStateVec.bbox(), t, "CalculateGrid::before f_step");
00104 
00105     // The ghostcells must remain untouched during integration!
00106     ((step_func_type) f_step)(mx[0],mx[1],base::NGhosts(),base::NEquations(),
00107                               dx(0),dx(1),NewStateVec.data(),OldStateVec.data(), 
00108                               dt,cfl,Flux[0]->data(),Flux[2]->data());
00109     
00110     if (NCheck()>=0) 
00111       // The negation captures cfl=NaN!
00112       if (ControlGrid(NewStateVec,level,NewStateVec.bbox(),t,0)==0 || !(cfl>0)) {
00113         ResetGridFluxes(Flux);
00114         NewStateVec.copy(OldStateVec);
00115         if (cfl>0) cfl = std::max(cfl, 2.0);
00116         return cfl;
00117       }
00118 
00119     if (NCheck()%10>1)
00120       base::CheckGrid(NewStateVec, level, NewStateVec.bbox(), t, "CalculateGrid::after f_step");
00121 
00122     return cfl;
00123   } 
00124 
00125   virtual void AllocGridFluxes(const BBox &bb, vec_grid_data_type**& Flux) {
00126     // The solver assumes a CONSECUTIVE array for the fluxes in ALL directions!
00127     // We need a member variable to store it between different member function calls.
00128     FluxData = new VectorType[base::Dim()*bb.size()];
00129     Flux = new vec_grid_data_type* [2*base::Dim()]; 
00130     for (register int d=0; d<base::Dim(); d++) {
00131       Flux[2*d] = new vec_grid_data_type(bb,&(FluxData[d*bb.size()])); 
00132       Flux[2*d+1] = Flux[2*d];
00133     }
00134   }
00135 
00136   virtual void DeAllocGridFluxes(vec_grid_data_type**& Flux) {
00137     if (Flux) {
00138       VectorType* dummy;
00139       for (register int d=0; d<base::Dim(); d++) 
00140         if (Flux[2*d]) {
00141           // Erease pointer to data
00142           Flux[2*d]->deallocate(dummy);
00143           delete Flux[2*d];
00144         }
00145       delete [] Flux;
00146     }
00147     if (FluxData) {
00148       delete [] FluxData;
00149       FluxData = (VectorType*) 0;
00150     }
00151   }
00152     
00153   virtual void ResetGridFluxes(vec_grid_data_type**& Flux) {
00154     if (Flux) {
00155       VectorType zero(0.0);
00156       for (register int d=0; d<base::Dim(); d++) 
00157         if (Flux[2*d]) 
00158           Flux[2*d]->equals(zero);
00159     }
00160   }
00161 
00162   virtual int ControlGrid(vec_grid_data_type& StateVec, const int& level, 
00163                           const BBox& where, const double& time, const int verbose) {
00164     int result = 1;
00165     if (!f_chk) return result;
00166     if (dim == 1) 
00167       ((check_1_func_type) f_chk)(FA(1,StateVec),BOUNDING_BOX(where),base::NEquations(),
00168                                   verbose,result);
00169     else if (dim == 2) 
00170       ((check_2_func_type) f_chk)(FA(2,StateVec),BOUNDING_BOX(where),base::NEquations(),
00171                                   verbose,result);
00172     else if (dim == 3) 
00173       ((check_3_func_type) f_chk)(FA(3,StateVec),BOUNDING_BOX(where),base::NEquations(),
00174                                   verbose,result);
00175     return result; 
00176   }  
00177 
00178   inline const int& NCheck() const { return _check; }
00179   virtual int NMethodOrder() const { return _order; }
00180   inline int NMaxPass() const { return 0; }
00181 
00182   inline void SetCheckFunc(generic_func_type check) { f_chk = check; }
00183   generic_func_type GetCheckFunc() const { return f_chk; }
00184   inline void SetStepFunc(generic_func_type step) { f_step = step; }
00185   generic_func_type GetStepFunc() const { return f_step; }
00186 
00187 protected:
00188   generic_func_type f_step, f_chk;
00189   std::string _name;
00190   int _order, _check;
00191   VectorType* FluxData;
00192 };
00193 
00194 #endif