00001
00002
00003
00004
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
00073
00074
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]==1) 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
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