amroc/balans/BalansIntegrator.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_BALANS_INTEGRATOR_H
00007 #define AMROC_BALANS_INTEGRATOR_H
00008 
00013 #include "Integrator.h"
00014 
00015 #include <string>
00016 
00017 #define f_init_common FORTRAN_NAME(combl, COMBL)
00018 
00019 extern "C" {
00020   void f_init_common();
00021 }
00022 
00026 template <class VectorType, int dim>
00027 class BalansIntegrator : 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& mx, const INTEGER& my, 
00044                                   const INTEGER& mbc, const INTEGER& meqn,
00045                                   const DOUBLE x[], const DOUBLE y[], 
00046                                   const DOUBLE& dx, const DOUBLE& dy, 
00047                                   VectorType q[], 
00048                                   const DOUBLE& dt , DOUBLE& cfl, 
00049                                   VectorType fx[], VectorType fy[]);
00050 
00051   BalansIntegrator(generic_func_type step) : 
00052     base(), f_step(step), f_chk(0),
00053     _order(2), _check(-1) { 
00054     _name = ""; 
00055     FluxData = (VectorType*) 0;
00056   }
00057 
00058   BalansIntegrator(generic_func_type step, generic_func_type chk) : 
00059     base(), f_step(step), f_chk(chk),
00060     _order(2), _check(0) { 
00061     _name = ""; 
00062     FluxData = (VectorType*) 0;
00063   }
00064 
00065   ~BalansIntegrator() {
00066     if (FluxData) {
00067       delete [] FluxData;
00068       FluxData = (VectorType*) 0;
00069     }  
00070   }
00071 
00072   virtual void register_at(ControlDevice& Ctrl, const std::string& prefix) {
00073     ControlDevice LocCtrl = Ctrl.getSubDevice(prefix+"BalansIntegrator");
00074     _name = prefix+"BalansIntegrator";
00075     RegisterAt(LocCtrl,"Check",_check);
00076   }
00077   virtual void register_at(ControlDevice& Ctrl) {
00078     register_at(Ctrl, "");
00079   }
00080    
00081   virtual void SetupData(GridHierarchy* gh, const int& ghosts) {  
00082     base::SetupData(gh,ghosts);
00083     base::SetMaxIntegratorPasses(NMaxPass());
00084     base::_abort = (NCheck()>=10 ? DAGHTrue : DAGHFalse);
00085     f_init_common();
00086   }
00087 
00088   virtual double CalculateGrid(vec_grid_data_type& NewStateVec, 
00089                                vec_grid_data_type& OldStateVec,
00090                                vec_grid_data_type* Flux[],
00091                                const int& level, const double& t, 
00092                                const double& dt, const int& mpass) {
00093     double cfl = 0.0;
00094     NewStateVec.copy(OldStateVec);
00095 
00096     Coords ex = OldStateVec.extents();
00097     DCoords lbcorner = base::GH().worldCoords(OldStateVec.lower(), OldStateVec.stepsize());
00098     DCoords dx = base::GH().worldStep(OldStateVec.stepsize());
00099     int maxm[3], mx[3], d;
00100     DataType* x[3];
00101     for (d=0; d<dim; d++) {
00102       maxm[d] = ex(d);
00103       mx[d] = ex(d)-2*base::NGhosts();
00104       x[d] = new DataType[maxm[d]];
00105       for (int i=0; i<maxm[d]; i++) x[d][i] = (i+0.5)*dx(d)+lbcorner(d);
00106     }
00107 
00108     if (NCheck()>0)
00109       base::CheckGrid(OldStateVec, level, OldStateVec.bbox(), t, "CalculateGrid::before f_step");
00110 
00111     // The ghostcells must remain untouched during integration!
00112     ((step_func_type) f_step)(mx[0],mx[1],base::NGhosts(),base::NEquations(),
00113                               x[0],x[1],dx(0),dx(1),NewStateVec.data(), 
00114                               dt,cfl,Flux[0]->data(),Flux[2]->data());
00115     
00116     if (NCheck()>=0) 
00117       // The negation captures cfl=NaN!
00118       if (ControlGrid(NewStateVec,level,NewStateVec.bbox(),t,0)==0 || !(cfl>0)) {
00119         ResetGridFluxes(Flux);
00120         NewStateVec.copy(OldStateVec);
00121         if (cfl>0) cfl = std::max(cfl, 2.0);
00122         return cfl;
00123       }
00124 
00125     if (NCheck()%10>1)
00126       base::CheckGrid(NewStateVec, level, NewStateVec.bbox(), t, "CalculateGrid::after f_step");
00127 
00128     return cfl;
00129   } 
00130 
00131   virtual void AllocGridFluxes(const BBox &bb, vec_grid_data_type**& Flux) {
00132     // The solver assumes a CONSECUTIVE array for the fluxes in ALL directions!
00133     // We need a member variable to store it between different member function calls.
00134     FluxData = new VectorType[base::Dim()*bb.size()];
00135     Flux = new vec_grid_data_type* [2*base::Dim()]; 
00136     for (register int d=0; d<base::Dim(); d++) {
00137       Flux[2*d] = new vec_grid_data_type(bb,&(FluxData[d*bb.size()])); 
00138       Flux[2*d+1] = Flux[2*d];
00139     }
00140   }
00141 
00142   virtual void DeAllocGridFluxes(vec_grid_data_type**& Flux) {
00143     if (Flux) {
00144       VectorType* dummy;
00145       for (register int d=0; d<base::Dim(); d++) 
00146         if (Flux[2*d]) {
00147           // Erease pointer to data
00148           Flux[2*d]->deallocate(dummy);
00149           delete Flux[2*d];
00150         }
00151       delete [] Flux;
00152     }
00153     if (FluxData) {
00154       delete [] FluxData;
00155       FluxData = (VectorType*) 0;
00156     }
00157   }
00158     
00159   virtual void ResetGridFluxes(vec_grid_data_type**& Flux) {
00160     if (Flux) {
00161       VectorType zero(0.0);
00162       for (register int d=0; d<base::Dim(); d++) 
00163         if (Flux[2*d]) 
00164           Flux[2*d]->equals(zero);
00165     }
00166   }
00167 
00168   virtual int ControlGrid(vec_grid_data_type& StateVec, const int& level, 
00169                           const BBox& where, const double& time, const int verbose) {
00170     int result = 1;
00171     if (!f_chk) return result;
00172     if (dim == 1) 
00173       ((check_1_func_type) f_chk)(FA(1,StateVec),BOUNDING_BOX(where),base::NEquations(),
00174                                   verbose,result);
00175     else if (dim == 2) 
00176       ((check_2_func_type) f_chk)(FA(2,StateVec),BOUNDING_BOX(where),base::NEquations(),
00177                                   verbose,result);
00178     else if (dim == 3) 
00179       ((check_3_func_type) f_chk)(FA(3,StateVec),BOUNDING_BOX(where),base::NEquations(),
00180                                   verbose,result);
00181     return result; 
00182   }  
00183 
00184   inline const int& NCheck() const { return _check; }
00185   virtual int NMethodOrder() const { return _order; }
00186   inline int NMaxPass() const { return 0; }
00187 
00188   inline void SetCheckFunc(generic_func_type check) { f_chk = check; }
00189   generic_func_type GetCheckFunc() const { return f_chk; }
00190   inline void SetStepFunc(generic_func_type step) { f_step = step; }
00191   generic_func_type GetStepFunc() const { return f_step; }
00192 
00193 protected:
00194   generic_func_type f_step, f_chk;
00195   std::string _name;
00196   int _order, _check;
00197   VectorType* FluxData;
00198 };
00199 
00200 #endif