amroc/clawpack/ClpIntegrator1.h

Go to the documentation of this file.
00001 // -*- C++ -*-
00002 
00003 // Copyright (C) 2002 Ralf Deiterding
00004 // Brandenburgische Universitaet Cottbus
00005 
00006 #ifndef AMROC_CLP_INTEGRATOR1_H
00007 #define AMROC_CLP_INTEGRATOR1_H
00008 
00016 #include "ClpIntegratorBase.h"
00017 
00018 #define f_step_1 FORTRAN_NAME(step1, STEP1)
00019 
00020 extern "C" {
00021   void f_step_1();
00022 }
00023 
00033 template <class VectorType, class AuxVectorType>
00034 class ClpIntegrator<VectorType,AuxVectorType,1> : public ClpIntegratorBase<VectorType,AuxVectorType,1> {
00035   typedef typename VectorType::InternalDataType DataType;
00036   typedef ClpIntegratorBase<VectorType,AuxVectorType,1> base;
00037 
00038 public:
00039   typedef typename base::vec_grid_fct_type vec_grid_fct_type;  
00040   typedef typename base::vec_grid_data_type vec_grid_data_type;
00041   typedef typename base::ld_vec_grid_data_type ld_vec_grid_data_type;
00042   typedef typename base::ld_aux_grid_data_type ld_aux_grid_data_type;
00043   typedef typename base::generic_func_type generic_func_type;
00044 
00045   typedef void (*normal_1_func_type) ( const INTEGER& maxmx, const INTEGER& meqn,
00046                                        const INTEGER& mwaves, const INTEGER& mbc, const INTEGER& mx,
00047                                        const DOUBLE ql[], const DOUBLE qr[],
00048                                        const INTEGER& maux, const DOUBLE auxl[], const DOUBLE auxr[],
00049                                        DOUBLE wave[], DOUBLE s[],
00050                                        DOUBLE amdq[], DOUBLE apdq[] );
00051   typedef void (*step_1_func_type) (const INTEGER& maxmx, const INTEGER& mvar, const INTEGER& meqn, 
00052                                     const INTEGER& maux, const INTEGER& mwaves, const INTEGER& mbc,
00053                                     const INTEGER& mx,
00054                                     VectorType qold[], const DOUBLE aux[],           
00055                                     const DOUBLE& dx, const DOUBLE& t, const DOUBLE& dt, 
00056                                     INTEGER method[], INTEGER mthlim[], DOUBLE& cfl, 
00057                                     VectorType fm[], VectorType fp[], 
00058                                     DOUBLE wave[], DOUBLE s[], 
00059                                     DOUBLE amdq[], DOUBLE apdq[], 
00060                                     DOUBLE dtdx[], DOUBLE aux1[], DOUBLE q1[], 
00061                                     DOUBLE work[], const INTEGER& mwork, 
00062                                     normal_1_func_type rpn);
00063 
00064   ClpIntegrator(const int equ, const int wv, generic_func_type rpn) : 
00065     base(equ,wv), f_step(f_step_1), f_rpn(rpn) {}
00066   ClpIntegrator(const int equ, const int wv, generic_func_type rpn,
00067                 generic_func_type check, generic_func_type aux, 
00068                 generic_func_type source) : 
00069     base(equ,wv,check,aux,source), f_step(f_step_1), f_rpn(rpn) {}
00070 
00071   //------------------------------------------------------------------
00072   // Solve Riemann problems. Copy left and right data into appropriate
00073   // format for normal_flux_func(). Note that left state has to shifted
00074   // one element. This is done in copySlice().
00075   //------------------------------------------------------------------                      
00076   virtual void ComputeRP(ld_vec_grid_data_type& LeftState, ld_aux_grid_data_type& auxl, 
00077                          ld_vec_grid_data_type& RightState, ld_aux_grid_data_type& auxr, 
00078                          ld_vec_grid_data_type& NewState, const int& direction) {
00079 
00080     AllocError::SetTexts("ClpIntegrator<1>","calculate_Riemann_Problem(): work");
00081     Coords ex = LeftState.extents();
00082     int maxm = LeftState.extents(0) + 1;
00083     int nghosts = 0;
00084     DataType* wave = new DataType[maxm*base::NEqUsed()*base::NWaves()];
00085     DataType* s = new DataType[maxm*base::NWaves()];
00086     DataType* LeftStateSl  = new DataType[maxm*base::NEqUsed()];
00087     DataType* RightStateSl = new DataType[maxm*base::NEqUsed()];
00088     DataType* auxlSl = new DataType[maxm*base::NAux()];
00089     DataType* auxrSl = new DataType[maxm*base::NAux()];
00090 
00091     AllocError::SetTexts("ClpIntegrator<1>","calculate_Riemann_Problem(): fluctuations");
00092     DataType* apdqSl = new DataType[(NewState.bbox().size()+1)*base::NEqUsed()];
00093     DataType* amdqSl = new DataType[(NewState.bbox().size()+1)*base::NEqUsed()];
00094 
00095     copySlice(LeftStateSl, LeftState, base::NEqUsed(), 1);
00096     copySlice(RightStateSl, RightState, base::NEqUsed(), 0);
00097         
00098     copyAuxSlice(auxlSl, auxl, base::NAux(), 1);
00099     copyAuxSlice(auxrSl, auxr, base::NAux(), 0);
00100 
00101     ((normal_1_func_type) f_rpn)(maxm, base::NEqUsed(), base::NWaves(), nghosts, maxm,
00102                                  LeftStateSl, RightStateSl,
00103                                  base::NAux(), auxlSl, auxrSl, wave, s, amdqSl, apdqSl);
00104 
00105     addSlices(NewState, apdqSl, amdqSl, base::NEqUsed(), 1);
00106     
00107     delete [] LeftStateSl; delete [] RightStateSl;
00108     delete [] auxlSl; delete [] auxrSl;
00109     delete [] apdqSl; delete [] amdqSl;
00110     delete [] wave;
00111     delete [] s;
00112   }
00113 
00114   virtual double ComputeGrid(vec_grid_data_type& StateVec, 
00115                              vec_grid_data_type* Flux[], DataType* aux,
00116                              const double& t, const double& dt, const int& mpass) {
00117     
00118     assert(f_step && f_rpn);
00119     AllocError::SetTexts("ClpIntegrator<1>","calculate_time_step(): work");
00120     Coords ex = StateVec.extents();
00121     DCoords dx = base::GH().worldStep(StateVec.stepsize());
00122     int mx = ex(0)-2*base::NGhosts();
00123     int msl = 0, mssl = 0;
00124 #ifdef EXTENDED_INTEGRATOR 
00125     msl = 4; if (base::method[0]>0) mssl = 2;
00126 #endif  
00127     int mwork = (mx + 2*base::NGhosts()) * ((base::NEqUsed()+1)*base::NWaves() + 
00128                                             (3+msl)*base::NEqUsed()+ 
00129                                             base::NAux() + 1 + mssl*base::NEquations());      
00130     DataType* work = new DataType[mwork];
00131      
00132     int wave = 0;
00133     int s    = wave + (mx+2*base::NGhosts())*base::NEqUsed()*base::NWaves();
00134     int amdq =    s + (mx+2*base::NGhosts())*base::NWaves();
00135     int apdq = amdq + (mx+2*base::NGhosts())*base::NEqUsed();
00136     int dtdx = apdq + (mx+2*base::NGhosts())*base::NEqUsed();
00137     int aux1 = dtdx + (mx+2*base::NGhosts());
00138     int q1   = aux1 + (mx+2*base::NGhosts())*base::NAux();
00139     int next = q1 + (mx+2*base::NGhosts())*base::NEqUsed();
00140     int mwork1 = mwork - next;
00141     double cfl = 0.0;
00142 
00143     // The ghostcells must remain untouched during integration!
00144     ((step_1_func_type) f_step)(mx,base::NEquations(),base::NEqUsed(),
00145                                 base::NAux(),base::NWaves(),base::NGhosts(),mx,
00146                                 StateVec.data(),
00147                                 aux,dx(0),t,dt,base::method,base::mthlim,cfl, 
00148                                 Flux[0]->data(),Flux[1]->data(),
00149                                 &(work[wave]),&(work[s]),
00150                                 &(work[amdq]),&(work[apdq]),
00151                                 &(work[dtdx]),&(work[aux1]),&(work[q1]),
00152                                 &(work[next]),mwork1,
00153                                 ((normal_1_func_type) f_rpn));
00154 
00155     delete [] work;    
00156     return cfl;
00157   }
00158 
00159 private:
00160   void copySlice(DataType* target, const ld_vec_grid_data_type &source, 
00161                  const int meqn, const int move) {
00162     const BBox& bbox = source.bbox();
00163     for (int m=0; m<meqn; m++) {
00164       if (move) target++;
00165       BeginFastIndex1(tgt, bbox, target, DataType);
00166       for_1 (k, bbox, bbox.stepsize()) 
00167         FastIndex1(tgt, k) = source(k)(m);
00168       end_for
00169       EndFastIndex1(tgt);
00170       if (!move) target++;
00171       target = &(target[bbox.size()]);
00172     }
00173   }
00174 
00175   void copyAuxSlice(DataType* target, const ld_aux_grid_data_type &source, 
00176                     const int meqn, const int move) {
00177     const BBox& bbox = source.bbox();
00178     for (int m=0; m<meqn; m++) {
00179       if (move) target++;
00180       BeginFastIndex1(tgt, bbox, target, DataType);
00181       for_1 (k, bbox, bbox.stepsize()) 
00182         FastIndex1(tgt, k) = source(k)(m);
00183       end_for
00184       EndFastIndex1(tgt);
00185       if (!move) target++;
00186       target = &(target[bbox.size()]);
00187     }
00188   }
00189 
00190   void addSlices(ld_vec_grid_data_type &target,
00191                  DataType* source1, DataType* source2, 
00192                  const int meqn, const int move) {
00193     const BBox& bbox = target.bbox();
00194     for (int m=0; m<meqn; m++) {
00195       if (move) { source1++; source2++; } 
00196       BeginFastIndex1(src1, bbox, source1, const DataType);
00197       BeginFastIndex1(src2, bbox, source2, const DataType);
00198       for_1 (k, bbox, bbox.stepsize()) 
00199         target(k)(m) = FastIndex1(src1, k) + FastIndex1(src2, k);
00200       end_for
00201       EndFastIndex1(src1);
00202       EndFastIndex1(src2);
00203       if (!move) { source1++; source2++; } 
00204       source1 = &(source1[bbox.size()]);
00205       source2 = &(source2[bbox.size()]);
00206     }
00207   }
00208   
00209 public:
00210   inline void SetStepFunc(generic_func_type step) { f_step = step; }
00211   generic_func_type GetStepFunc() const { return f_step; }
00212   inline void SetRpnFunc(generic_func_type rpn) { f_rpn = rpn; }
00213   generic_func_type GetRpnFunc() const { return f_rpn; }
00214 
00215 protected:
00216   generic_func_type f_step, f_rpn;
00217 };
00218 
00219 
00220 #endif