00001
00002
00003
00004
00005
00006 #ifndef AMROC_CLP_INTEGRATOR2_H
00007 #define AMROC_CLP_INTEGRATOR2_H
00008
00016 #include "ClpIntegratorBase.h"
00017
00018 #define f_step_2 FORTRAN_NAME(step2, STEP2)
00019
00020 extern "C" {
00021 void f_step_2();
00022 }
00023
00033 template <class VectorType, class AuxVectorType>
00034 class ClpIntegrator<VectorType,AuxVectorType,2> : public ClpIntegratorBase<VectorType,AuxVectorType,2> {
00035 typedef typename VectorType::InternalDataType DataType;
00036 typedef ClpIntegratorBase<VectorType,AuxVectorType,2> 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_2_func_type) ( const INTEGER& ixyz, const INTEGER& maxm,
00046 const INTEGER& meqn,
00047 const INTEGER& mwaves, const INTEGER& mbc, const INTEGER& mx,
00048 const DOUBLE ql[], const DOUBLE qr[],
00049 const INTEGER& maux, const DOUBLE auxl[], const DOUBLE auxr[],
00050 DOUBLE wave[], DOUBLE s[],
00051 DOUBLE amdq[], DOUBLE apdq[] );
00052 typedef void (*transverse_2_func_type) ( const INTEGER& ixyz, const INTEGER& maxm,
00053 const INTEGER& meqn,
00054 const INTEGER& mwaves, const INTEGER& mbc, const INTEGER& mx,
00055 const DOUBLE ql[], const DOUBLE qr[],
00056 const INTEGER& maux, const DOUBLE aux1[], const DOUBLE aux2[],
00057 const DOUBLE aux3[], const INTEGER& imp,
00058 DOUBLE asdq[], DOUBLE amdq[], DOUBLE apdq[] );
00059 typedef void (*step_2_func_type) (const INTEGER& mp, const INTEGER& maxm,
00060 const INTEGER& maxmx, const INTEGER& maxmy,
00061 const INTEGER& mvar, const INTEGER& meqn,
00062 const INTEGER& maux, const INTEGER& mwaves, const INTEGER& mbc,
00063 const INTEGER& mx,const INTEGER& my,
00064 VectorType qold[], const DOUBLE aux[],
00065 const DOUBLE& dx, const DOUBLE& dy,
00066 const DOUBLE& t, const DOUBLE& dt,
00067 INTEGER method[], INTEGER mthlim[], DOUBLE& cfl,
00068 VectorType fm[], VectorType fp[],
00069 VectorType gm[], VectorType gp[],
00070 DOUBLE faddm[], DOUBLE faddp[],
00071 DOUBLE gaddm[], DOUBLE gaddp[],
00072 DOUBLE q1d[], DOUBLE dtdx1d[], DOUBLE dtdy1d[],
00073 DOUBLE aux1[], DOUBLE aux2[], DOUBLE aux3[],
00074 DOUBLE work[], const INTEGER& mwork,
00075 normal_2_func_type rpn, transverse_2_func_type rpt);
00076
00077 ClpIntegrator(const int equ, const int wv, generic_func_type rpn, generic_func_type rpt) :
00078 base(equ,wv), f_step(f_step_2), f_rpn(rpn), f_rpt(rpt) {}
00079 ClpIntegrator(const int equ, const int wv, generic_func_type rpn, generic_func_type rpt,
00080 generic_func_type check, generic_func_type aux, generic_func_type source) :
00081 base(equ,wv,check,aux,source), f_step(f_step_2), f_rpn(rpn), f_rpt(rpt) {}
00082
00083
00084
00085
00086
00087
00088 virtual void ComputeRP(ld_vec_grid_data_type& LeftState, ld_aux_grid_data_type& auxl,
00089 ld_vec_grid_data_type& RightState, ld_aux_grid_data_type& auxr,
00090 ld_vec_grid_data_type& NewState, const int& direction) {
00091
00092 AllocError::SetTexts("ClpIntegrator<2>","ComputeRP(): work");
00093 int maxm = LeftState.extents(0) + 1;
00094 int mx = maxm-2*base::NGhosts();
00095 DataType* wave = new DataType[maxm*base::NEqUsed()*base::NWaves()];
00096 DataType* s = new DataType[maxm*base::NWaves()];
00097 DataType* LeftStateSl = new DataType[maxm*base::NEqUsed()];
00098 DataType* RightStateSl = new DataType[maxm*base::NEqUsed()];
00099 DataType* auxlSl = new DataType[maxm*base::NAux()];
00100 DataType* auxrSl = new DataType[maxm*base::NAux()];
00101
00102 AllocError::SetTexts("ClpIntegrator<2>","calculate_Riemann_Problem(): fluctuations");
00103 DataType* apdqSl = new DataType[maxm*base::NEqUsed()];
00104 DataType* amdqSl = new DataType[maxm*base::NEqUsed()];
00105
00106 copySlice(LeftStateSl, LeftState, base::NEqUsed(), 1);
00107 copySlice(RightStateSl, RightState, base::NEqUsed(), 0);
00108
00109 copyAuxSlice(auxlSl, auxl, base::NAux(), 1);
00110 copyAuxSlice(auxrSl, auxr, base::NAux(), 0);
00111
00112 ((normal_2_func_type) f_rpn)(direction, mx, base::NEqUsed(), base::NWaves(), base::NGhosts(), mx,
00113 LeftStateSl, RightStateSl,
00114 base::NAux(), auxlSl, auxrSl, wave, s, amdqSl, apdqSl);
00115
00116 addSlices(NewState, apdqSl, amdqSl, base::NEqUsed(), 1);
00117
00118 delete [] LeftStateSl; delete [] RightStateSl;
00119 delete [] auxlSl; delete [] auxrSl;
00120 delete [] apdqSl; delete [] amdqSl;
00121 delete [] wave;
00122 delete [] s;
00123 }
00124
00125 virtual double ComputeGrid(vec_grid_data_type& StateVec,
00126 vec_grid_data_type* Flux[], DataType* aux,
00127 const double& t, const double& dt, const int& mpass) {
00128
00129 assert(f_step && f_rpn && f_rpt);
00130 AllocError::SetTexts("ClpIntegrator<2>","calculate_time_step(): work");
00131 Coords ex = StateVec.extents();
00132 DCoords dx = base::GH().worldStep(StateVec.stepsize());
00133 int mx = ex(0)-2*base::NGhosts();
00134 int my = ex(1)-2*base::NGhosts();
00135 int maxm = Max(mx,my);
00136 int msl = 0, mssl = 0;
00137 #ifdef EXTENDED_INTEGRATOR
00138 msl = 4; if (base::method[0]==1) mssl = 4;
00139 #endif
00140
00141 int mwork = (mx + 2*base::NGhosts())*(my + 2*base::NGhosts())*mssl*base::NEquations() +
00142 (maxm + 2*base::NGhosts())*((12+msl+base::NWaves())*base::NEqUsed() +
00143 base::NWaves() + 3*base::NAux() + 2);
00144 DataType* work = new DataType[mwork];
00145
00146 int faddm = 0;
00147 int faddp = faddm + (maxm+2*base::NGhosts())*base::NEqUsed();
00148 int gaddm = faddp + (maxm+2*base::NGhosts())*base::NEqUsed();
00149 int gaddp = gaddm + 2*(maxm+2*base::NGhosts())*base::NEqUsed();
00150 int q1d = gaddp + 2*(maxm+2*base::NGhosts())*base::NEqUsed();
00151 int dtdx1d = q1d + (maxm+2*base::NGhosts())*base::NEqUsed();
00152 int dtdy1d = dtdx1d + (maxm+2*base::NGhosts());
00153 int aux1 = dtdy1d + (maxm+2*base::NGhosts());
00154 int aux2 = aux1 + (maxm+2*base::NGhosts())*base::NAux();
00155 int aux3 = aux2 + (maxm+2*base::NGhosts())*base::NAux();
00156 int next = aux3 + (maxm+2*base::NGhosts())*base::NAux();
00157 int mwork1 = mwork - next;
00158 double cfl = 0.0;
00159
00160
00161 ((step_2_func_type) f_step)(mpass,maxm,mx,my,base::NEquations(),base::NEqUsed(),
00162 base::NAux(),base::NWaves(),base::NGhosts(),mx,my,
00163 StateVec.data(),
00164 aux,dx(0),dx(1),t,dt,base::method,base::mthlim,cfl,
00165 Flux[0]->data(),Flux[1]->data(),
00166 Flux[2]->data(),Flux[3]->data(),
00167 &(work[faddm]),&(work[faddp]),
00168 &(work[gaddm]),&(work[gaddp]),
00169 &(work[q1d]),
00170 &(work[dtdx1d]),&(work[dtdy1d]),
00171 &(work[aux1]),&(work[aux2]),&(work[aux3]),
00172 &(work[next]),mwork1,
00173 ((normal_2_func_type) f_rpn),
00174 ((transverse_2_func_type) f_rpt));
00175
00176 delete [] work;
00177 return cfl;
00178 }
00179
00180 private:
00181 void copySlice(DataType* target, const ld_vec_grid_data_type &source,
00182 const int meqn, const int move) {
00183 const BBox& bbox = source.bbox();
00184 BeginFastIndex1(src, bbox, source.data(), const VectorType);
00185 for (int m=0; m<meqn; m++) {
00186 if (move) target++;
00187 BeginFastIndex1(tgt, bbox, target, DataType);
00188 for_1 (k, bbox, bbox.stepsize())
00189 FastIndex1(tgt, k) = FastIndex1(src, k)(m);
00190 end_for
00191 EndFastIndex1(tgt);
00192 if (!move) target++;
00193 target = &(target[bbox.size()]);
00194 }
00195 EndFastIndex1(src);
00196 }
00197
00198 void copyAuxSlice(DataType* target, const ld_aux_grid_data_type &source,
00199 const int meqn, const int move) {
00200 const BBox& bbox = source.bbox();
00201 BeginFastIndex1(src, bbox, source.data(), const AuxVectorType);
00202 for (int m=0; m<meqn; m++) {
00203 if (move) target++;
00204 BeginFastIndex1(tgt, bbox, target, DataType);
00205 for_1 (k, bbox, bbox.stepsize())
00206 FastIndex1(tgt, k) = FastIndex1(src, k)(m);
00207 end_for
00208 EndFastIndex1(tgt);
00209 if (!move) target++;
00210 target = &(target[bbox.size()]);
00211 }
00212 EndFastIndex1(src);
00213 }
00214
00215 void addSlices(ld_vec_grid_data_type &target,
00216 DataType* source1, DataType* source2,
00217 const int meqn, const int move) {
00218 const BBox& bbox = target.bbox();
00219 BeginFastIndex1(tgt, bbox, target.data(), VectorType);
00220 for (int m=0; m<meqn; m++) {
00221 if (move) { source1++; source2++; }
00222 BeginFastIndex1(src1, bbox, source1, const DataType);
00223 BeginFastIndex1(src2, bbox, source2, const DataType);
00224 for_1 (k, bbox, bbox.stepsize())
00225 FastIndex1(tgt, k)(m) = FastIndex1(src1, k) + FastIndex1(src2, k);
00226 end_for
00227 EndFastIndex1(src1);
00228 EndFastIndex1(src2);
00229 if (!move) { source1++; source2++; }
00230 source1 = &(source1[bbox.size()]);
00231 source2 = &(source2[bbox.size()]);
00232 }
00233 EndFastIndex1(tgt);
00234 }
00235
00236 public:
00237 inline void SetStepFunc(generic_func_type step) { f_step = step; }
00238 generic_func_type GetStepFunc() const { return f_step; }
00239 inline void SetRpnFunc(generic_func_type rpn) { f_rpn = rpn; }
00240 generic_func_type GetRpnFunc() const { return f_rpn; }
00241 inline void SetRptFunc(generic_func_type rpt) { f_rpt = rpt; }
00242 generic_func_type GetRptFunc() const { return f_rpt; }
00243
00244 protected:
00245 generic_func_type f_step, f_rpn, f_rpt;
00246 };
00247
00248
00249
00250 #endif