00001
00002
00003
00004
00005
00006
00007
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
00092
00093
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
00122
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]==1) 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
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