00001
00002
00003
00004
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
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
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
00133
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
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