amroc/clawpack/ClpIntegrator3.h

Go to the documentation of this file.
00001 // -*- C++ -*-
00002 
00003 // Copyright (C) 2002 Ralf Deiterding
00004 // Brandenburgische Universitaet Cottbus
00005 //
00006 // Copyright (C) 2003-2007 California Institute of Technology
00007 // Ralf Deiterding, ralf@cacr.caltech.edu
00008 
00009 #ifndef AMROC_CLP_INTEGRATOR3_H
00010 #define AMROC_CLP_INTEGRATOR3_H
00011 
00019 #include "ClpIntegratorBase.h"
00020 
00021 #define f_step_3 FORTRAN_NAME(step3, STEP3)
00022 
00023 extern "C" {
00024   void f_step_3();
00025 }
00026 
00036 template <class VectorType, class AuxVectorType>
00037 class ClpIntegrator<VectorType,AuxVectorType,3> : public ClpIntegratorBase<VectorType,AuxVectorType,3> {
00038   typedef typename VectorType::InternalDataType DataType;
00039   typedef ClpIntegratorBase<VectorType,AuxVectorType,3> base;
00040 
00041 public:
00042   typedef typename base::vec_grid_fct_type vec_grid_fct_type;  
00043   typedef typename base::vec_grid_data_type vec_grid_data_type;
00044   typedef typename base::ld_vec_grid_data_type ld_vec_grid_data_type;
00045   typedef typename base::ld_aux_grid_data_type ld_aux_grid_data_type;
00046   typedef typename base::generic_func_type generic_func_type;
00047 
00048   typedef void (*normal_3_func_type) ( const INTEGER& ixyz, const INTEGER& maxm, 
00049                                        const INTEGER& meqn,
00050                                        const INTEGER& mwaves, const INTEGER& mbc, const INTEGER& mx,
00051                                        const DOUBLE ql[], const DOUBLE qr[],
00052                                        const INTEGER& maux, const DOUBLE auxl[], const DOUBLE auxr[],
00053                                        DOUBLE wave[], DOUBLE s[],
00054                                        DOUBLE amdq[], DOUBLE apdq[] );
00055   typedef void (*transverse_3_func_type) ( const INTEGER& ixyz, const INTEGER& icoor,
00056                                            const INTEGER& maxm, const INTEGER& meqn,
00057                                            const INTEGER& mwaves, const INTEGER& mbc, const INTEGER& mx,
00058                                            const DOUBLE ql[], const DOUBLE qr[],
00059                                            const INTEGER& maux, const DOUBLE aux1[],
00060                                            const DOUBLE aux2[], const DOUBLE aux3[],
00061                                            const INTEGER& ilr, DOUBLE asdq[],
00062                                            DOUBLE bmasdq[], DOUBLE bpasdq[] );
00063   typedef void (*step_3_func_type) ( const INTEGER& mp, const INTEGER& maxm, 
00064                                      const INTEGER& maxmx, const INTEGER& maxmy, const INTEGER& maxmz,
00065                                      const INTEGER& mvar, const INTEGER& meqn, 
00066                                      const INTEGER& maux, const INTEGER& mwaves, const INTEGER& mbc,
00067                                      const INTEGER& mx, const INTEGER& my, const INTEGER& mz,
00068                                      VectorType qold[], const DOUBLE aux[],          
00069                                      const DOUBLE& dx, const DOUBLE& dy, const DOUBLE& dz,
00070                                      const DOUBLE& t, const DOUBLE& dt, INTEGER method[], 
00071                                      INTEGER mthlim[], DOUBLE& cfl, 
00072                                      VectorType fm[], VectorType fp[], 
00073                                      VectorType gm[], VectorType gp[], 
00074                                      VectorType hm[], VectorType hp[], 
00075                                      DOUBLE faddm[], DOUBLE faddp[], 
00076                                      DOUBLE gaddm[], DOUBLE gaddp[], 
00077                                      DOUBLE haddm[], DOUBLE haddp[],
00078                                      DOUBLE q1d[], 
00079                                      DOUBLE dtdx1d[], DOUBLE dtdy1d[], DOUBLE dtdz1d[],
00080                                      DOUBLE aux1[], DOUBLE aux2[], DOUBLE aux3[], 
00081                                      DOUBLE work[], const INTEGER& mwork, 
00082                                      normal_3_func_type rpn, transverse_3_func_type rpt );
00083 
00084   ClpIntegrator(const int equ, const int wv, generic_func_type rpn, generic_func_type rpt) : 
00085     base(equ,wv), f_step(f_step_3), f_rpn(rpn), f_rpt(rpt) {}
00086   ClpIntegrator(const int equ, const int wv, generic_func_type rpn, generic_func_type rpt,
00087                 generic_func_type check, generic_func_type aux, generic_func_type source) : 
00088     base(equ,wv,check,aux,source), f_step(f_step_3), f_rpn(rpn), f_rpt(rpt) {}
00089 
00090   //------------------------------------------------------------------
00091   // Solve Riemann problems. Copy left and right data into appropriate
00092   // format for normal_flux_func(). Note that left state has to shifted
00093   // one element. This is done in copySlice().
00094   //------------------------------------------------------------------                      
00095   virtual void ComputeRP(ld_vec_grid_data_type& LeftState, ld_aux_grid_data_type& auxl, 
00096                          ld_vec_grid_data_type& RightState, ld_aux_grid_data_type& auxr, 
00097                          ld_vec_grid_data_type& NewState, const int& direction) {
00098 
00099     AllocError::SetTexts("ClpIntegrator<3>","ComputeRP(): work");
00100     int maxm = LeftState.extents(0) + 1;
00101     int mx = maxm-2*base::NGhosts();
00102     DataType* wave = new DataType[maxm*base::NEqUsed()*base::NWaves()];
00103     DataType* s = new DataType[maxm*base::NWaves()];
00104 
00105     AllocError::SetTexts("ClpIntegrator<3>","calculate_Riemann_Problem(): slices");
00106     DataType* LeftStateSl  = new DataType[maxm*base::NEqUsed()];
00107     DataType* RightStateSl = new DataType[maxm*base::NEqUsed()];
00108     DataType* apdqSl = new DataType[maxm*base::NEqUsed()];
00109     DataType* amdqSl = new DataType[maxm*base::NEqUsed()];
00110     DataType* auxlSl = new DataType[maxm*3*base::NAux()];
00111     DataType* auxrSl = new DataType[maxm*3*base::NAux()];
00112 
00113     BBox bb(LeftState.bbox());
00114     BBox bbSl(bb);   
00115     AlignBBoxToSlice(bbSl);
00116     for (register int y=bb.lower(1); y<=bb.upper(1); y+=bb.stepsize(1)) {
00117 
00118       copySlice(LeftStateSl, bbSl, LeftState, y, base::NEqUsed(), 1);
00119       copySlice(RightStateSl, bbSl, RightState, y, base::NEqUsed(), 0);
00120         
00121       // auxiliary data along neighboring slices generally aren't needed in 
00122       // the rpn3 routine.
00123       copyAuxSlice(&(auxlSl[maxm*base::NAux()]), bbSl, auxl, y, base::NAux(), 1);
00124       copyAuxSlice(&(auxrSl[maxm*base::NAux()]), bbSl, auxr, y, base::NAux(), 0);
00125 
00126       ((normal_3_func_type) f_rpn)(direction, mx, base::NEqUsed(), base::NWaves(), base::NGhosts(), mx,
00127                                    LeftStateSl, RightStateSl,  
00128                                    base::NAux(), auxlSl, auxrSl, wave, s, amdqSl, apdqSl);
00129     
00130       addSlices(NewState, apdqSl, amdqSl, bbSl, y, base::NEqUsed(), 1);
00131     }
00132     
00133     delete [] LeftStateSl; delete [] RightStateSl;
00134     delete [] auxlSl; delete [] auxrSl;
00135     delete [] apdqSl; delete [] amdqSl;
00136     delete [] wave;
00137     delete [] s;
00138   }
00139 
00140   virtual double ComputeGrid(vec_grid_data_type& StateVec, 
00141                              vec_grid_data_type* Flux[], DataType* aux,
00142                              const double& t, const double& dt, const int& mpass) {
00143 
00144     assert(f_step && f_rpn && f_rpt);
00145     AllocError::SetTexts("ClpIntegrator<3>","calculate_time_step(): work");
00146     Coords ex = StateVec.extents();
00147     DCoords dx = base::GH().worldStep(StateVec.stepsize());
00148     int mx = ex(0)-2*base::NGhosts();
00149     int my = ex(1)-2*base::NGhosts();      
00150     int mz = ex(2)-2*base::NGhosts();      
00151     int maxm = Max(Max(mx,my),mz);
00152     int msl = 0, mssl = 0;
00153 #ifdef EXTENDED_INTEGRATOR
00154     msl = 4; if (base::method[0]>0) mssl = 6; 
00155 #endif  
00156     int mwork = (mx + 2*base::NGhosts())*(my + 2*base::NGhosts())*
00157       (mz + 2*base::NGhosts())*mssl*base::NEquations() + 
00158       (maxm + 2*base::NGhosts())*((54+msl+base::NWaves())*base::NEqUsed() + 
00159                                   base::NWaves() + 9*base::NAux() + 3);      
00160     DataType* work = new DataType[mwork];   
00161      
00162     int faddm = 0;
00163     int faddp = faddm + (maxm+2*base::NGhosts())*base::NEqUsed();
00164     int gaddm = faddp + (maxm+2*base::NGhosts())*base::NEqUsed();
00165     int gaddp = gaddm + 6*(maxm+2*base::NGhosts())*base::NEqUsed();
00166     int haddm = gaddp + 6*(maxm+2*base::NGhosts())*base::NEqUsed();
00167     int haddp = haddm + 6*(maxm+2*base::NGhosts())*base::NEqUsed();
00168     int q1d   = haddp + 6*(maxm+2*base::NGhosts())*base::NEqUsed();
00169     int dtdx1d = q1d + (maxm+2*base::NGhosts())*base::NEqUsed();
00170     int dtdy1d = dtdx1d + (maxm+2*base::NGhosts());
00171     int dtdz1d = dtdy1d + (maxm+2*base::NGhosts());
00172     int aux1 = dtdz1d + (maxm+2*base::NGhosts());
00173     int aux2 = aux1 + 3*(maxm+2*base::NGhosts())*base::NAux();
00174     int aux3 = aux2 + 3*(maxm+2*base::NGhosts())*base::NAux();
00175     int next = aux3 + 3*(maxm+2*base::NGhosts())*base::NAux();
00176     int mwork1 = mwork - next;
00177     double cfl = 0.0;
00178 
00179     // The ghostcells must remain untouched during integration!
00180     ((step_3_func_type) f_step)(mpass,maxm,mx,my,mz,base::NEquations(),base::NEqUsed(),
00181                                 base::NAux(),base::NWaves(),base::NGhosts(),mx,my,mz,
00182                                 StateVec.data(),
00183                                 aux,dx(0),dx(1),dx(2),t,dt,base::method,base::mthlim,cfl, 
00184                                 Flux[0]->data(),Flux[1]->data(),
00185                                 Flux[2]->data(),Flux[3]->data(), 
00186                                 Flux[4]->data(),Flux[5]->data(), 
00187                                 &(work[faddm]),&(work[faddp]),
00188                                 &(work[gaddm]),&(work[gaddp]),
00189                                 &(work[haddm]),&(work[haddp]),
00190                                 &(work[q1d]),
00191                                 &(work[dtdx1d]),&(work[dtdy1d]),&(work[dtdz1d]),
00192                                 &(work[aux1]),&(work[aux2]),&(work[aux3]),
00193                                 &(work[next]),mwork1,
00194                                 ((normal_3_func_type) f_rpn), 
00195                                 ((transverse_3_func_type) f_rpt));
00196     
00197     delete [] work;
00198     return cfl;
00199   }
00200 
00201 private:
00202   void AlignBBoxToSlice(BBox &bb) {
00203     gdbAlignBBox(1, bb, DAGH_X);
00204     bb.rank = 1;
00205     bb.lower().rank = 1;
00206     bb.upper().rank = 1;
00207     bb.stepsize().rank = 1;
00208   } 
00209   
00210   void copySlice(DataType* target, const BBox &targetbbox, 
00211                  const ld_vec_grid_data_type &source, 
00212                  const int y, const int meqn, const int move) {
00213     BBox bbox = targetbbox * source.bbox();
00214     BeginFastIndex2(src, source.bbox(), source.data(), const VectorType);
00215     for (int m=0; m<meqn; m++) {
00216       if (move) target++;
00217       BeginFastIndex1(tgt, targetbbox, target, DataType);
00218       for_1 (k, bbox, bbox.stepsize()) 
00219         FastIndex1(tgt, k) = FastIndex2(src, k, y)(m);
00220       end_for
00221       EndFastIndex1(tgt);
00222       if (!move) target++;
00223       target = &(target[targetbbox.size()]);
00224     }
00225     EndFastIndex2(src);
00226   }
00227 
00228   void copyAuxSlice(DataType* target, const BBox &targetbbox, 
00229                     const ld_aux_grid_data_type &source, 
00230                     const int y, const int meqn, const int move) {
00231     BBox bbox = targetbbox * source.bbox();
00232     BeginFastIndex2(src, source.bbox(), source.data(), const AuxVectorType);
00233     for (int m=0; m<meqn; m++) {
00234       if (move) target++;
00235       BeginFastIndex1(tgt, targetbbox, target, DataType);
00236       for_1 (k, bbox, bbox.stepsize()) 
00237         FastIndex1(tgt, k) = FastIndex2(src, k, y)(m);
00238       end_for
00239       EndFastIndex1(tgt);
00240       if (!move) target++;
00241       target = &(target[targetbbox.size()]);
00242     }
00243     EndFastIndex2(src);
00244   }
00245 
00246   void addSlices(ld_vec_grid_data_type &target,
00247                  DataType* source1, DataType* source2, const BBox &sourcebbox, 
00248                  const int y, const int meqn, const int move) {
00249     BBox bbox = target.bbox() * sourcebbox;
00250     BeginFastIndex2(tgt, target.bbox(), target.data(), VectorType);
00251     for (int m=0; m<meqn; m++) {
00252       if (move) { source1++; source2++; } 
00253       BeginFastIndex1(src1, sourcebbox, source1, const DataType);
00254       BeginFastIndex1(src2, sourcebbox, source2, const DataType);
00255       for_1 (k, bbox, bbox.stepsize()) 
00256         FastIndex2(tgt, k, y)(m) = FastIndex1(src1, k) + FastIndex1(src2, k);
00257       end_for
00258       EndFastIndex1(src1);
00259       EndFastIndex1(src2);
00260       if (!move) { source1++; source2++; } 
00261       source1 = &(source1[sourcebbox.size()]);
00262       source2 = &(source2[sourcebbox.size()]);
00263     }
00264     EndFastIndex2(tgt);
00265   }
00266 
00267 public:
00268   inline void SetStepFunc(generic_func_type step) { f_step = step; }
00269   generic_func_type GetStepFunc() const { return f_step; }
00270   inline void SetRpnFunc(generic_func_type rpn) { f_rpn = rpn; }
00271   generic_func_type GetRpnFunc() const { return f_rpn; }
00272   inline void SetRptFunc(generic_func_type rpt) { f_rpt = rpt; }
00273   generic_func_type GetRptFunc() const { return f_rpt; }
00274 
00275 protected:
00276   generic_func_type f_step, f_rpn, f_rpt;
00277 };
00278 
00279 
00280 #endif