00001
00002
00003
00004
00005
00006 #ifndef EGLM2D_H
00007 #define EGLM2D_H
00008
00015 #include "Interfaces/SchemeBase.h"
00016
00023 template <class DataType>
00024 class EGLM2D : public SchemeBase<Vector<DataType,9>, 2> {
00025 typedef SchemeBase<Vector<DataType,9>, 2> base;
00026 public:
00027 typedef typename base::vec_grid_data_type vec_grid_data_type;
00028 typedef typename base::grid_data_type grid_data_type;
00029 typedef typename base::VectorType VectorType;
00030 typedef typename base::SideName SideName;
00031 typedef typename base::point_type point_type;
00032
00033 enum ICPredefined { Constant, DiscontinuityX, DiscontinuityY, DiscontinuityXY };
00034 enum BCPredefined { Symmetry, SlipWall, Inlet, Outlet };
00035 enum GFMPredefined { GFMExtrapolation, GFMSlipWall };
00036
00037 EGLM2D() : base() {
00038 SchemeNb = 1;
00039 no_clean = 0;
00040 gamma = DataType(5./3.);
00041 epsilon = DataType(1.e-8);
00042 cfl_clean = 1.;
00043 cr = DataType(0.18);
00044 om = DataType(0.);
00045 Limiter = 0;
00046 }
00047
00048 virtual ~EGLM2D() {}
00049
00050 virtual void register_at(ControlDevice& Ctrl, const std::string& prefix) {
00051 base::LocCtrl = Ctrl.getSubDevice(prefix+"EGLM2D");
00052 RegisterAt(base::LocCtrl,"Scheme",SchemeNb);
00053 RegisterAt(base::LocCtrl,"Limiter",Limiter);
00054 RegisterAt(base::LocCtrl,"Omega",om);
00055 RegisterAt(base::LocCtrl,"Gamma",gamma);
00056 RegisterAt(base::LocCtrl,"CFLClean",cfl_clean);
00057 RegisterAt(base::LocCtrl,"NoClean",no_clean);
00058 RegisterAt(base::LocCtrl,"cr",cr);
00059 }
00060
00061 virtual void SetupData(GridHierarchy* gh, const int& ghosts) {
00062 base::SetupData(gh,ghosts);
00063 WriteInit();
00064 }
00065
00066 virtual double Step(vec_grid_data_type& vec, vec_grid_data_type& oldvec, vec_grid_data_type* Flux[],
00067 const double& t, const double& dt, const int& mpass) const {
00068
00069 grid_data_type rE(vec.bbox());
00070 DCoords dxm = base::GH().worldStep(vec.stepsize());
00071 double dx = dxm(0), dy = dxm(1);
00072
00073
00074 double ch = cfl_clean*std::min(dx,dy)/dt;
00075 double cfl_grid = 0.;
00076
00077 if (mpass<=1) {
00078 Flux[0]->equals(VectorType(0.)); Flux[2]->equals(VectorType(0.));
00079 DataType xSlopeMax = DataType(0.), ySlopeMax = DataType(0.);
00080
00081 conservative2primitive(vec,rE);
00082
00083 if (Limiter<=0) {
00084 vec_grid_data_type Fx(vec.bbox()), Fy(vec.bbox());
00085
00086 computeFluxX(vec,Fx,rE);
00087 NumericalFluxX(vec,vec,Fx,Fx,rE,rE,*Flux[0],xSlopeMax,ch);
00088
00089 computeFluxY(vec,Fy,rE);
00090 NumericalFluxY(vec,vec,Fy,Fy,rE,rE,*Flux[2],ySlopeMax,ch);
00091 }
00092
00093 else {
00094 vec_grid_data_type ql(vec.bbox()), qr(vec.bbox());
00095 vec_grid_data_type Fl(vec.bbox()), Fr(vec.bbox());
00096 grid_data_type rEl(vec.bbox()), rEr(vec.bbox());
00097
00098 MUSCLX(vec,qr,ql,Fr,Fl,rEr,rEl,dt,dx,ch);
00099 computeFluxX(ql,Fl,rEl); computeFluxX(qr,Fr,rEr);
00100 NumericalFluxX(qr,ql,Fr,Fl,rEr,rEl,*Flux[0],xSlopeMax,ch);
00101
00102 MUSCLY(vec,qr,ql,Fr,Fl,rEr,rEl,dt,dy,ch);
00103 computeFluxY(ql,Fl,rEl); computeFluxY(qr,Fr,rEr);
00104 NumericalFluxY(qr,ql,Fr,Fl,rEr,rEl,*Flux[2],ySlopeMax,ch);
00105 }
00106
00107
00108 DataType slopeMax = std::max(xSlopeMax/dx,ySlopeMax/dy);
00109 cfl_grid = dt*slopeMax;
00110
00111
00112 conservationStep(vec,rE,*Flux[0],*Flux[2],dt,dx,dy);
00113 }
00114
00115 else if (mpass==2 && !no_clean) {
00116 conservative2primitive(vec,rE);
00117
00118
00119 computeCorrection(vec,rE,ch,dt,dx,dy);
00120 }
00121
00122 return cfl_grid;
00123 }
00124
00125 virtual void ICStandard(vec_grid_data_type &vec, const int type,
00126 DataType* aux=0, const int naux=0, const int scaling=0) const {
00127 int Nx = vec.extents(0), Ny = vec.extents(1);
00128 int mxs = base::NGhosts(), mys = base::NGhosts();
00129 int mxe = Nx-base::NGhosts()-1, mye = Ny-base::NGhosts()-1;
00130 VectorType *mesh = (VectorType *)vec.databuffer();
00131
00132 DCoords lbcorner = base::GH().worldCoords(vec.lower(), vec.stepsize());
00133 DCoords dx = base::GH().worldStep(vec.stepsize());
00134
00135 if (type==Constant) {
00136 VectorType g;
00137 for (register int i=0; i<base::NEquations(); i++)
00138 g(i) = aux[i];
00139 for (register int j=mys; j<=mye; j++)
00140 for (register int i=mxs; i<=mxe; i++)
00141 mesh[base::idx(i,j,Nx)] = g;
00142 }
00143 else if (type==DiscontinuityX) {
00144 VectorType gl, gr;
00145 for (register int i=0; i<base::NEquations(); i++) {
00146 gl(i) = aux[i]; gr(i) = aux[base::NEquations()+i];
00147 }
00148 for (register int j=mys; j<=mye; j++)
00149 for (register int i=mxs; i<=mxe; i++)
00150 if ((i+0.5)*dx(0)+lbcorner(0) < aux[2*base::NEquations()])
00151 mesh[base::idx(i,j,Nx)] = gl;
00152 else
00153 mesh[base::idx(i,j,Nx)] = gr;
00154 }
00155 else if (type==DiscontinuityY) {
00156 VectorType gl, gr;
00157 for (register int i=0; i<base::NEquations(); i++) {
00158 gl(i) = aux[i]; gr(i) = aux[base::NEquations()+i];
00159 }
00160 for (register int j=mys; j<=mye; j++)
00161 for (register int i=mxs; i<=mxe; i++) {
00162 if ((j+0.5)*dx(1)+lbcorner(1) < aux[2*base::NEquations()])
00163 mesh[base::idx(i,j,Nx)] = gl;
00164 else
00165 mesh[base::idx(i,j,Nx)] = gr;
00166 }
00167 }
00168 else if (type==DiscontinuityXY) {
00169 VectorType g1, g2, g3, g4;
00170 for (register int i=0; i<base::NEquations(); i++) {
00171 g1(i) = aux[i];
00172 g2(i) = aux[base::NEquations()+i];
00173 g3(i) = aux[2*base::NEquations()+i];
00174 g4(i) = aux[3*base::NEquations()+i];
00175 }
00176 for (register int j=mys; j<=mye; j++)
00177 for (register int i=mxs; i<=mxe; i++) {
00178 double xc = (i+0.5)*dx(0)+lbcorner(0);
00179 double yc = (j+0.5)*dx(1)+lbcorner(1);
00180 register int pos = base::idx(i,j,Nx);
00181
00182 if (xc>=aux[4*base::NEquations()] && yc>=aux[4*base::NEquations()+1])
00183 mesh[pos] = g1;
00184 else if (xc<aux[4*base::NEquations()] && yc>=aux[4*base::NEquations()+1])
00185 mesh[pos] = g2;
00186 else if (xc<aux[4*base::NEquations()] && yc<aux[4*base::NEquations()+1])
00187 mesh[pos] = g3;
00188 else
00189 mesh[pos] = g4;
00190 }
00191 }
00192
00193 if (scaling==base::Primitive) {
00194 grid_data_type rhoE(vec.bbox());
00195 computeRhoE(vec,rhoE);
00196 primitive2conservative(vec,rhoE);
00197 }
00198
00199 }
00200
00201
00202 virtual void BCStandard(vec_grid_data_type &vec, const BBox &bb, const int type, const int side,
00203 DataType* aux=0, const int naux=0, const int scaling=0) const {
00204 int Nx = vec.extents(0);
00205 int b = vec.bottom();
00206 assert (vec.stepsize(0)==bb.stepsize(0) && vec.stepsize(1)==bb.stepsize(1));
00207 int mxs = std::max(bb.lower(0)/bb.stepsize(0),vec.lower(0)/vec.stepsize(0));
00208 int mys = std::max(bb.lower(1)/bb.stepsize(1),vec.lower(1)/vec.stepsize(1));
00209 int mxe = std::min(bb.upper(0)/bb.stepsize(0),vec.upper(0)/vec.stepsize(0));
00210 int mye = std::min(bb.upper(1)/bb.stepsize(1),vec.upper(1)/vec.stepsize(1));
00211 VectorType *mesh = (VectorType *)vec.databuffer();
00212
00213 if (type==Symmetry || type==SlipWall) {
00214 switch (side) {
00215 case base::Left:
00216 for (register int i=mxe; i>=mxs; i--) {
00217 register int isym = 2*mxe+1-i;
00218 for (register int j=mys; j<=mye; j++) {
00219 mesh[b+base::idx(i,j,Nx)] = mesh[b+base::idx(isym,j,Nx)];
00220 mesh[b+base::idx(i,j,Nx)](3) = -mesh[b+base::idx(isym,j,Nx)](3);
00221 if (type==Symmetry)
00222 mesh[b+base::idx(i,j,Nx)](6) = -mesh[b+base::idx(isym,j,Nx)](6);
00223 }
00224 }
00225 break;
00226 case base::Right:
00227 for (register int i=mxs; i<=mxe; i++) {
00228 register int isym = 2*mxs-1-i;
00229 for (register int j=mys; j<=mye; j++) {
00230 mesh[b+base::idx(i,j,Nx)] = mesh[b+base::idx(isym,j,Nx)];
00231 mesh[b+base::idx(i,j,Nx)](3) = -mesh[b+base::idx(isym,j,Nx)](3);
00232 if (type==Symmetry)
00233 mesh[b+base::idx(i,j,Nx)](6) = -mesh[b+base::idx(isym,j,Nx)](6);
00234 }
00235 }
00236 break;
00237 case base::Bottom:
00238 for (register int j=mye; j>=mys; j--) {
00239 register int jsym = 2*mye+1-j;
00240 for (register int i=mxs; i<=mxe; i++) {
00241 mesh[b+base::idx(i,j,Nx)] = mesh[b+base::idx(i,jsym,Nx)];
00242 mesh[b+base::idx(i,j,Nx)](4) = -mesh[b+base::idx(i,jsym,Nx)](4);
00243 if (type==Symmetry)
00244 mesh[b+base::idx(i,j,Nx)](7) = -mesh[b+base::idx(i,jsym,Nx)](7);
00245 }
00246 }
00247 break;
00248 case base::Top:
00249 for (register int j=mys; j<=mye; j++) {
00250 register int jsym = 2*mys-1-j;
00251 for (register int i=mxs; i<=mxe; i++) {
00252 mesh[b+base::idx(i,j,Nx)] = mesh[b+base::idx(i,jsym,Nx)];
00253 mesh[b+base::idx(i,j,Nx)](4) = -mesh[b+base::idx(i,jsym,Nx)](4);
00254 if (type==Symmetry)
00255 mesh[b+base::idx(i,j,Nx)](7) = -mesh[b+base::idx(i,jsym,Nx)](7);
00256 }
00257 }
00258 break;
00259 }
00260 }
00261
00262 else if (type==Inlet) {
00263 if (aux==0 || naux<9) return;
00264 VectorType g;
00265 for (register int i=0; i<base::NEquations(); i++)
00266 g(i) = aux[i];
00267 if (scaling==base::Primitive) {
00268 DataType B2 = DataType(0.5)*(g(6)*g(6)+g(7)*g(7)+g(8)*g(8));
00269 DataType v2 = DataType(0.5)*(g(3)*g(3)+g(4)*g(4)+g(5)*g(5));
00270 g(1) = g(0)*v2 + g(1)/(gamma-DataType(1.))+B2;
00271 g(3) *= g(0); g(4) *= g(0); g(5) *= g(0);
00272 }
00273 switch (side) {
00274 case base::Left:
00275 for (register int i=mxe; i>=mxs; i--)
00276 for (register int j=mys; j<=mye; j++)
00277 mesh[b+base::idx(i,j,Nx)] = g;
00278 break;
00279 case base::Right:
00280 for (register int i=mxs; i<=mxe; i++)
00281 for (register int j=mys; j<=mye; j++)
00282 mesh[b+base::idx(i,j,Nx)] = g;
00283 break;
00284 case base::Bottom:
00285 for (register int j=mye; j>=mys; j--)
00286 for (register int i=mxs; i<=mxe; i++)
00287 mesh[b+base::idx(i,j,Nx)] = g;
00288 break;
00289 case base::Top:
00290 for (register int j=mys; j<=mye; j++)
00291 for (register int i=mxs; i<=mxe; i++)
00292 mesh[b+base::idx(i,j,Nx)] = g;
00293 break;
00294 }
00295 }
00296
00297 else if (type==Outlet) {
00298 switch (side) {
00299 case base::Left:
00300 for (register int i=mxe; i>=mxs; i--)
00301 for (register int j=mys; j<=mye; j++)
00302 mesh[b+base::idx(i,j,Nx)] = mesh[b+base::idx(i+1,j,Nx)];
00303 break;
00304 case base::Right:
00305 for (register int i=mxs; i<=mxe; i++)
00306 for (register int j=mys; j<=mye; j++)
00307 mesh[b+base::idx(i,j,Nx)] = mesh[b+base::idx(i-1,j,Nx)];
00308 break;
00309 case base::Bottom:
00310 for (register int j=mye; j>=mys; j--)
00311 for (register int i=mxs; i<=mxe; i++)
00312 mesh[b+base::idx(i,j,Nx)] = mesh[b+base::idx(i,j+1,Nx)];
00313 break;
00314 case base::Top:
00315 for (register int j=mys; j<=mye; j++)
00316 for (register int i=mxs; i<=mxe; i++)
00317 mesh[b+base::idx(i,j,Nx)] = mesh[b+base::idx(i,j-1,Nx)];
00318 break;
00319 }
00320 }
00321 }
00322
00323 virtual bool GFMUseTransform() const { return true; }
00324 virtual void GFMTransform(vec_grid_data_type &vec, vec_grid_data_type &vechelp) const {
00325 grid_data_type rE(vec.bbox());
00326 vechelp.copy(vec);
00327 conservative2primitive(vechelp,rE);
00328 }
00329
00330 virtual void GFMBCStandard(vec_grid_data_type &vec, const int type, const int& nc, const int* idx,
00331 const VectorType* f, const point_type* xc, const DataType* distance,
00332 const point_type* normal, DataType* aux=0, const int naux=0,
00333 const int scaling=0) const {
00334
00335
00336
00337 if (type==GFMSlipWall) {
00338 for (int n=0; n<nc; n++) {
00339 DataType rho = f[n](0), pre = f[n](1), psi = f[n](2);
00340 DataType u = -f[n](3), v = -f[n](4), w = f[n](5);
00341 DataType Bx = f[n](6), By = f[n](7), Bz = f[n](8);
00342
00343
00344 if (naux>=2) {
00345 u += aux[naux*n];
00346 v += aux[naux*n+1];
00347 }
00348
00349 DataType vl = DataType(2.)*(u*normal[n](0)+v*normal[n](1));
00350 u = f[n](3) + vl*normal[n](0);
00351 v = f[n](4) + vl*normal[n](1);
00352
00353 VectorType &mesh = vec(Coords(base::Dim(),idx+base::Dim()*n));
00354 mesh(0) = rho;
00355 DataType B2 = DataType(0.5)*(Bx*Bx+By*By+Bz*Bz);
00356 DataType v2 = DataType(0.5)*(u*u+v*v+w*w);
00357 mesh(1) = rho*v2 + pre/(gamma-1)+B2;
00358 mesh(2) = psi;
00359 mesh(3) = u*rho; mesh(4) = v*rho; mesh(5) = w*rho;
00360 mesh(6) = Bx; mesh(7) = By; mesh(8) = Bz;
00361 }
00362 }
00363
00364
00365 else if (type==GFMExtrapolation) {
00366 for (int n=0; n<nc; n++) {
00367 DataType rho = f[n](0), pre = f[n](1), psi = f[n](2);
00368 DataType u = -f[n](3), v = -f[n](4), w = f[n](5);
00369 DataType Bx = f[n](6), By = f[n](7), Bz = f[n](8);
00370
00371
00372 DataType vl = u*normal[n](0)+v*normal[n](1);
00373 u = vl*normal[n](0);
00374 v = vl*normal[n](1);
00375
00376
00377 DataType Bl = Bx*normal[n](0)+By*normal[n](1);
00378 Bx = Bl*normal[n](0);
00379 By = Bl*normal[n](1);
00380
00381 VectorType &mesh = vec(Coords(base::Dim(),idx+base::Dim()*n));
00382 mesh(0) = rho;
00383 DataType B2 = DataType(0.5)*(Bx*Bx+By*By+Bz*Bz);
00384 DataType v2 = DataType(0.5)*(u*u+v*v+w*w);
00385 mesh(1) = rho*v2 + pre/(gamma-1)+B2;
00386 mesh(2) = psi;
00387 mesh(3) = u*rho; mesh(4) = v*rho; mesh(5) = w*rho;
00388 mesh(6) = Bx; mesh(7) = By; mesh(8) = Bz;
00389 }
00390 }
00391 }
00392
00393 virtual void Input(vec_grid_data_type& fvec, grid_data_type& work, const int cnt,
00394 const int skip_ghosts=1) const {}
00395
00396 virtual void Output(vec_grid_data_type& vec, grid_data_type& workvec, const int cnt,
00397 const int skip_ghosts=1) const {
00398
00399 int Nx = vec.extents(0), Ny = vec.extents(1);
00400 int mxs = base::NGhosts()*skip_ghosts, mys = base::NGhosts()*skip_ghosts;
00401 int mxe = Nx-base::NGhosts()*skip_ghosts-1, mye = Ny-base::NGhosts()*skip_ghosts-1;
00402 VectorType *mesh = (VectorType *)vec.databuffer();
00403 DataType *work = (DataType *)workvec.databuffer();
00404 DCoords dx = base::GH().worldStep(vec.stepsize());
00405 if (skip_ghosts==0) {
00406 if (cnt==12 || cnt==14) { mxs+=1; mxe-=1; }
00407 if (cnt==12 || cnt==15) { mys+=1; mye-=1; }
00408 }
00409
00410 for (register int j=mys; j<=mye; j++)
00411 for (register int i=mxs; i<=mxe; i++) {
00412 register int pos = base::idx(i,j,Nx);
00413
00414 switch(cnt) {
00415 case 1:
00416 work[pos]=mesh[pos](0);
00417 break;
00418 case 2:
00419 work[pos]=mesh[pos](3)/mesh[pos](0);
00420 break;
00421 case 3:
00422 work[pos]=mesh[pos](4)/mesh[pos](0);
00423 break;
00424 case 4:
00425 work[pos]=mesh[pos](1);
00426 break;
00427 case 6:
00428 work[pos]=(mesh[pos](1)-
00429 (mesh[pos](3)*mesh[pos](3)+mesh[pos](4)*mesh[pos](4)+
00430 mesh[pos](5)*mesh[pos](5))/mesh[pos](0)/DataType(2.)-
00431 (mesh[pos](6)*mesh[pos](6)+mesh[pos](7)*mesh[pos](7)+
00432 mesh[pos](8)*mesh[pos](8))/DataType(2.))*(gamma-DataType(1.));
00433 break;
00434 case 8:
00435 work[pos]=mesh[pos](5)/mesh[pos](0);
00436 break;
00437 case 9:
00438 work[pos]=mesh[pos](6);
00439 break;
00440 case 10:
00441 work[pos]=mesh[pos](7);
00442 break;
00443 case 11:
00444 work[pos]=mesh[pos](8);
00445 break;
00446 case 12:
00447 work[pos]=(mesh[base::idx(i+1,j,Nx)](6)-mesh[base::idx(i-1,j,Nx)](6))/(2.*dx(0)) +
00448 (mesh[base::idx(i,j+1,Nx)](7)-mesh[base::idx(i,j-1,Nx)](7))/(2.*dx(1));
00449 break;
00450 case 13:
00451 work[pos]=mesh[pos](2);
00452 break;
00453 case 14:
00454 work[pos]=(mesh[base::idx(i+1,j,Nx)](6)-mesh[base::idx(i-1,j,Nx)](6))/(2.*dx(0));
00455 break;
00456 case 15:
00457 work[pos]=(mesh[base::idx(i,j+1,Nx)](7)-mesh[base::idx(i,j-1,Nx)](7))/(2.*dx(1));
00458 break;
00459 case 16: {
00460 DataType u = mesh[pos](3)/mesh[pos](0);
00461 DataType v = mesh[pos](4)/mesh[pos](0);
00462 DataType w = mesh[pos](5)/mesh[pos](0);
00463 work[pos]=mesh[pos](6)*(mesh[pos](8)*v-mesh[pos](7)*w)+
00464 mesh[pos](7)*(mesh[pos](6)*w-mesh[pos](8)*u)+
00465 mesh[pos](8)*(mesh[pos](7)*u-mesh[pos](6)*v);
00466 break;
00467 }
00468 }
00469 }
00470 }
00471
00472 virtual int Check(vec_grid_data_type& vec, const BBox &bb, const double& time,
00473 const int verbose) const {
00474 int Nx = vec.extents(0);
00475 int b = vec.bottom();
00476 assert (vec.stepsize(0)==bb.stepsize(0) && vec.stepsize(1)==bb.stepsize(1));
00477 BBox bbmax = grow(vec.bbox(),-base::NGhosts());
00478 int mxs = std::max(bb.lower(0)/bb.stepsize(0),bbmax.lower(0)/bbmax.stepsize(0));
00479 int mys = std::max(bb.lower(1)/bb.stepsize(1),bbmax.lower(1)/bbmax.stepsize(1));
00480 int mxe = std::min(bb.upper(0)/bb.stepsize(0),bbmax.upper(0)/bbmax.stepsize(0));
00481 int mye = std::min(bb.upper(1)/bb.stepsize(1),bbmax.upper(1)/bbmax.stepsize(1));
00482 VectorType *mesh = (VectorType *)vec.databuffer();
00483
00484 int result = 1;
00485 for (register int j=mys; j<=mye; j++)
00486 for (register int i=mxs; i<=mxe; i++) {
00487 register int pos = b+base::idx(i,j,Nx);
00488 if (!(mesh[pos](0)>-1.e37)) {
00489 result = 0;
00490 if (verbose)
00491 std::cerr << "EGLM2D-Check(): mesh(" << i << "," << j << ")(0)="
00492 << mesh[pos](0) << " " << vec.bbox() << std::endl;
00493 }
00494 if (!(mesh[pos](1)>-1.e37)) {
00495 result = 0;
00496 if (verbose)
00497 std::cerr << "EGLM2D-Check(): mesh(" << i << "," << j << ")(1)="
00498 << mesh[pos](1) << " " << vec.bbox() << std::endl;
00499 }
00500 DataType p = (mesh[pos](1)-
00501 (mesh[pos](3)*mesh[pos](3)+mesh[pos](4)*mesh[pos](4)+
00502 mesh[pos](5)*mesh[pos](5))/mesh[pos](0)/DataType(2.)-
00503 (mesh[pos](6)*mesh[pos](6)+mesh[pos](7)*mesh[pos](7)+
00504 mesh[pos](8)*mesh[pos](8))/DataType(2.))*(gamma-DataType(1.));
00505
00506 if (!(p>-1.e37)) {
00507 result = 0;
00508 if (verbose)
00509 std::cerr << "EGLM2D-Check(): p(" << i << "," << j << ")=" << p << " " << vec.bbox() << std::endl;
00510 }
00511 }
00512 return result;
00513 }
00514
00515 virtual int NMethodOrder() const { return 1; }
00516 virtual int NMaxPass() const { return 2; }
00517
00518 protected:
00519 void WriteInit() const {
00520 int me = MY_PROC;
00521 if (me == VizServer) {
00522 std::ofstream ofs("chem.dat", std::ios::out);
00523 ofs << "RU 1.0" << std::endl;
00524 ofs << "PA 1.0" << std::endl;
00525 ofs << "Sp(01) w" << std::endl;
00526 ofs << "W(01) 1.0" << std::endl;
00527 ofs << "Sp(02) Bx" << std::endl;
00528 ofs << "W(02) 1.0" << std::endl;
00529 ofs << "Sp(03) By" << std::endl;
00530 ofs << "W(03) 1.0" << std::endl;
00531 ofs << "Sp(04) Bz" << std::endl;
00532 ofs << "W(04) 1.0" << std::endl;
00533 ofs << "Sp(05) Divergence" << std::endl;
00534 ofs << "W(05) 1.0" << std::endl;
00535 ofs << "Sp(06) Psi" << std::endl;
00536 ofs << "W(06) 1.0" << std::endl;
00537 ofs << "Sp(07) Div_Bx-comp" << std::endl;
00538 ofs << "W(07) 1.0" << std::endl;
00539 ofs << "Sp(08) Div_By-comp" << std::endl;
00540 ofs << "W(08) 1.0" << std::endl;
00541 ofs << "Sp(09) Helicity" << std::endl;
00542 ofs << "W(09) 1.0" << std::endl;
00543 ofs.close();
00544 }
00545 }
00546
00547
00548
00549 void computeFluxX(vec_grid_data_type &vec, vec_grid_data_type &fluxx, grid_data_type &rE) const
00550 {
00551 DataType half = DataType(0.5);
00552 DataType rho,pre,vx,vy,vz,vx2,Bx,By,Bz,Bx2,By2,Bz2,B2;
00553 int Nx = vec.extents(0), Ny = vec.extents(1);
00554 VectorType *mesh = (VectorType *)vec.databuffer();
00555 VectorType *Fx = (VectorType *)fluxx.databuffer();
00556 DataType *rhoE = (DataType *)rE.databuffer();
00557
00558 for (register int j=0; j<Ny; j++)
00559 for (register int i=0; i<Nx; i++) {
00560 register int pos = base::idx(i,j,Nx);
00561 rho = mesh[pos](0);
00562 pre = mesh[pos](1);
00563 vx = mesh[pos](3);
00564 vy = mesh[pos](4);
00565 vz = mesh[pos](5);
00566 Bx = mesh[pos](6);
00567 By = mesh[pos](7);
00568 Bz = mesh[pos](8);
00569
00570 Bx2 = Bx*Bx; By2 = By*By; Bz2 = Bz*Bz;
00571 B2 = half*(Bz2+Bx2+By2);
00572 vx2 = vx*vx;
00573
00574 Fx[pos](0) = rho*vx;
00575 Fx[pos](3) = rho*vx2 + pre + half*(Bz2+By2-Bx2);
00576 Fx[pos](4) = rho*vx*vy - Bx*By;
00577 Fx[pos](5) = rho*vx*vz - Bx*Bz;
00578 Fx[pos](1) = (rhoE[pos] + pre + B2)*vx - Bx*(vx*Bx + vy*By + vz*Bz);
00579 Fx[pos](6) = 0.0;
00580 Fx[pos](7) = vx*By - vy*Bx;
00581 Fx[pos](8) = vx*Bz - vz*Bx;
00582 Fx[pos](2) = 0.0;
00583 }
00584 }
00585
00586
00587
00588 void computeFluxY(vec_grid_data_type &vec, vec_grid_data_type &fluxy, grid_data_type &rE) const
00589 {
00590 DataType half = DataType(0.5);
00591 DataType rho,pre,vx,vy,vz,vy2,Bx,By,Bz,Bx2,By2,Bz2,B2;
00592 int Nx = vec.extents(0), Ny = vec.extents(1);
00593 VectorType *mesh = (VectorType *)vec.databuffer();
00594 VectorType *Fy = (VectorType *)fluxy.databuffer();
00595 DataType *rhoE = (DataType *)rE.databuffer();
00596
00597 for (register int j=0; j<Ny; j++)
00598 for (register int i=0; i<Nx; i++) {
00599 register int pos = base::idx(i,j,Nx);
00600 rho = mesh[pos](0);
00601 pre = mesh[pos](1);
00602 vx = mesh[pos](3);
00603 vy = mesh[pos](4);
00604 vz = mesh[pos](5);
00605 Bx = mesh[pos](6);
00606 By = mesh[pos](7);
00607 Bz = mesh[pos](8);
00608
00609 Bx2 = Bx*Bx; By2 = By*By; Bz2 = Bz*Bz;
00610 B2 = half*(Bz2+Bx2+By2);
00611 vy2 = vy*vy;
00612
00613 Fy[pos](0) = rho*vy;
00614 Fy[pos](3) = rho*vx*vy - Bx*By;
00615 Fy[pos](4) = rho*vy2 + pre + half*(Bx2+Bz2-By2);
00616 Fy[pos](5) = rho*vy*vz - By*Bz;
00617 Fy[pos](1) = (rhoE[pos] + pre + B2)*vy - By*(vx*Bx + vy*By + vz*Bz);
00618 Fy[pos](6) = vy*Bx - vx*By;
00619 Fy[pos](7) = 0.0;
00620 Fy[pos](8) = vy*Bz - vz*By;
00621 Fy[pos](2) = 0.0;
00622 }
00623 }
00624
00625 void NumericalFluxX(vec_grid_data_type &qrt, vec_grid_data_type &qlt, vec_grid_data_type &Frt, vec_grid_data_type &Flt,
00626 grid_data_type& rErt, grid_data_type& rElt, vec_grid_data_type &FluxX,
00627 DataType& xSlopeMax, DataType ch) const {
00628
00629 VectorType *ql = (VectorType *)qlt.databuffer(), *qr = (VectorType *)qrt.databuffer();
00630 VectorType *Fl = (VectorType *)Flt.databuffer(), *Fr = (VectorType *)Frt.databuffer();
00631 DataType *rEl = (DataType *)rElt.databuffer(), *rEr = (DataType *)rErt.databuffer();
00632 VectorType *F = (VectorType *)FluxX.databuffer();
00633
00634 int Nx = qrt.extents(0), Ny = qrt.extents(1);
00635 int mxs = base::NGhosts(), mys = base::NGhosts();
00636 int mxe = Nx-base::NGhosts()-1, mye = Ny-base::NGhosts()-1;
00637
00638 DataType auxRhoE[2];
00639 VectorType vc[2], auxF[2];
00640
00641 for(register int j=mys; j<=mye; j++)
00642 for(register int i=mxs; i<=mxe+1;i++) {
00643 vc[0] = qr[base::idx(i-1,j,Nx)];
00644 vc[1] = ql[base::idx(i,j,Nx)];
00645 auxF[0] = Fr[base::idx(i-1,j,Nx)];
00646 auxF[1] = Fl[base::idx(i,j,Nx)];
00647 auxRhoE[0] = rEr[base::idx(i-1,j,Nx)];
00648 auxRhoE[1] = rEl[base::idx(i,j,Nx)];
00649
00650 if(SchemeNb ==1) RiemannSolver_HLL(vc,auxRhoE,auxF,3,xSlopeMax,6);
00651 else RiemannSolver_HLLD(vc,auxRhoE,auxF,3,xSlopeMax,6);
00652
00653 F[base::idx(i,j,Nx)] = auxF[0];
00654 }
00655 if (!no_clean) fluxCorrectionX(qrt,qlt,FluxX,ch);
00656 }
00657
00658
00659 void NumericalFluxY(vec_grid_data_type &qrt, vec_grid_data_type &qlt, vec_grid_data_type &Frt, vec_grid_data_type &Flt,
00660 grid_data_type& rErt, grid_data_type& rElt, vec_grid_data_type &FluxY,
00661 DataType& ySlopeMax, DataType ch) const {
00662
00663 VectorType *ql = (VectorType *)qlt.databuffer(), *qr = (VectorType *)qrt.databuffer();
00664 VectorType *Fl = (VectorType *)Flt.databuffer(), *Fr = (VectorType *)Frt.databuffer();
00665 DataType *rEl = (DataType *)rElt.databuffer(), *rEr = (DataType *)rErt.databuffer();
00666 VectorType *G = (VectorType *)FluxY.databuffer();
00667
00668 int Nx = qrt.extents(0), Ny = qrt.extents(1);
00669 int mxs = base::NGhosts(), mys = base::NGhosts();
00670 int mxe = Nx-base::NGhosts()-1, mye = Ny-base::NGhosts()-1;
00671
00672 DataType auxRhoE[2];
00673 VectorType vc[2], auxF[2];
00674
00675 for(register int i=mxs; i<=mxe;i++)
00676 for(register int j=mys; j<=mye+1;j++){
00677 vc[0] = qr[base::idx(i,j-1,Nx)];
00678 vc[1] = ql[base::idx(i,j,Nx)];
00679 auxF[0] = Fr[base::idx(i,j-1,Nx)];
00680 auxF[1] = Fl[base::idx(i,j,Nx)];
00681 auxRhoE[0] = rEr[base::idx(i,j-1,Nx)];
00682 auxRhoE[1] = rEl[base::idx(i,j,Nx)];
00683
00684 if(SchemeNb ==1) RiemannSolver_HLL(vc,auxRhoE,auxF,4,ySlopeMax,7);
00685 else RiemannSolver_HLLD(vc,auxRhoE,auxF,4,ySlopeMax,7);
00686
00687 G[base::idx(i,j,Nx)] = auxF[0];
00688 }
00689
00690 if (!no_clean) fluxCorrectionY(qrt,qlt,FluxY,ch);
00691 }
00692
00693 void MUSCLX(vec_grid_data_type &vec, vec_grid_data_type &qrt, vec_grid_data_type &qlt,
00694 vec_grid_data_type &Frt, vec_grid_data_type &Flt, grid_data_type& rErt, grid_data_type& rElt,
00695 DataType dt, DataType dx, DataType ch) const {
00696
00697 VectorType *ql = (VectorType *)qlt.databuffer(), *qr = (VectorType *)qrt.databuffer();
00698 VectorType *Fl = (VectorType *)Flt.databuffer(), *Fr = (VectorType *)Frt.databuffer();
00699 VectorType *primitive = (VectorType *)vec.databuffer();
00700
00701 int Nx = vec.extents(0), Ny = vec.extents(1);
00702 int mxs = base::NGhosts(), mys = base::NGhosts();
00703 int mxe = Nx-base::NGhosts()-1, mye = Ny-base::NGhosts()-1;
00704
00705 qlt.copy(vec); qrt.copy(vec);
00706 for(register int j=mys; j<=mye; j++)
00707 for(register int i=mxs-1; i<=mxe+1; i++){
00708 register int pos = base::idx(i,j,Nx), posm = base::idx(i-1,j,Nx), posp = base::idx(i+1,j,Nx);
00709 for (register int n=0; n<base::NEquations(); n++)
00710 Reconstruct(primitive[pos](n),primitive[posm](n),primitive[posp](n),ql[pos](n),qr[pos](n));
00711 }
00712
00713 computeRhoE(qlt,rElt); computeRhoE(qrt,rErt);
00714 computeFluxX(qlt,Flt,rElt); computeFluxX(qrt,Frt,rErt);
00715 if (!no_clean) { fluxCorrectionX(qlt,qlt,Flt,ch); fluxCorrectionX(qrt,qrt,Frt,ch); }
00716 primitive2conservative(qlt,rElt); primitive2conservative(qrt,rErt);
00717 DataType hdtx = DataType(0.5)*dt/dx;
00718 for(register int j=mys;j<=mye;j++)
00719 for(register int i=mxs-1;i<=mxe+1;i++) {
00720 ql[base::idx(i,j,Nx)] -= hdtx*(Fr[base::idx(i,j,Nx)]-Fl[base::idx(i,j,Nx)]);
00721 qr[base::idx(i,j,Nx)] -= hdtx*(Fr[base::idx(i,j,Nx)]-Fl[base::idx(i,j,Nx)]);
00722 }
00723 conservative2primitive(qlt,rElt); conservative2primitive(qrt,rErt);
00724 }
00725
00726 void MUSCLY(vec_grid_data_type &vec, vec_grid_data_type &qrt, vec_grid_data_type &qlt,
00727 vec_grid_data_type &Frt, vec_grid_data_type &Flt, grid_data_type& rErt, grid_data_type& rElt,
00728 DataType dt, DataType dy, DataType ch) const {
00729
00730 VectorType *ql = (VectorType *)qlt.databuffer(), *qr = (VectorType *)qrt.databuffer();
00731 VectorType *Fl = (VectorType *)Flt.databuffer(), *Fr = (VectorType *)Frt.databuffer();
00732 VectorType *primitive = (VectorType *)vec.databuffer();
00733
00734 int Nx = vec.extents(0), Ny = vec.extents(1);
00735 int mxs = base::NGhosts(), mys = base::NGhosts();
00736 int mxe = Nx-base::NGhosts()-1, mye = Ny-base::NGhosts()-1;
00737
00738 qlt.copy(vec); qrt.copy(vec);
00739 for (register int i=mxs; i<=mxe; i++)
00740 for(register int j=mys-1; j<=mye+1; j++) {
00741 register int pos = base::idx(i,j,Nx), posm = base::idx(i,j-1,Nx), posp = base::idx(i,j+1,Nx);
00742 for (register int n=0; n<base::NEquations(); n++)
00743 Reconstruct(primitive[pos](n),primitive[posm](n),primitive[posp](n),ql[pos](n),qr[pos](n));
00744 }
00745
00746 computeRhoE(qlt,rElt); computeRhoE(qrt,rErt);
00747 computeFluxY(qlt,Flt,rElt); computeFluxY(qrt,Frt,rErt);
00748 if (!no_clean) { fluxCorrectionY(qlt,qlt,Flt,ch); fluxCorrectionY(qrt,qrt,Frt,ch); }
00749 primitive2conservative(qlt,rElt); primitive2conservative(qrt,rErt);
00750 DataType hdty = DataType(0.5)*dt/dy;
00751 for(register int i=mxs;i<=mxe;i++)
00752 for(register int j=mys-1;j<=mye+1;j++) {
00753 ql[base::idx(i,j,Nx)] -= hdty*(Fr[base::idx(i,j,Nx)]-Fl[base::idx(i,j,Nx)]);
00754 qr[base::idx(i,j,Nx)] -= hdty*(Fr[base::idx(i,j,Nx)]-Fl[base::idx(i,j,Nx)]);
00755 }
00756 conservative2primitive(qlt,rElt); conservative2primitive(qrt,rErt);
00757 }
00758
00759
00760
00761
00762
00763
00764 void RiemannSolver_HLLD(VectorType vec[2], DataType rhoE[2], VectorType F[2], int position, DataType& slopeMax, int BFieldPosition) const
00765 {
00766 DataType slopeLeft, slopeRight,slopeLeftStar,slopeRightStar,slopeM;
00767 VectorType Ustar[2],Ustarstar[2];
00768
00769 slope_HLLD (vec,position,slopeLeft,slopeRight,slopeMax,BFieldPosition);
00770
00771 computeStateUstarsHLLD (vec,rhoE,Ustar,Ustarstar,position,BFieldPosition,slopeLeft,slopeRight,slopeLeftStar,slopeRightStar,slopeM);
00772
00773
00774
00775 vec [0][1] = rhoE[0]; vec [1][1] = rhoE[1];
00776 vec [0][3] = vec[0][0]*vec[0][3]; vec [1][3] = vec[1][0]*vec[1][3];
00777 vec [0][4] = vec[0][0]*vec[0][4]; vec [1][4] = vec[1][0]*vec[1][4];
00778 vec [0][5] = vec[0][0]*vec[0][5]; vec [1][5] = vec[1][0]*vec[1][5];
00779
00780
00781
00782 if(slopeLeft>0)F[0] = F[0];
00783
00784 else if(slopeLeftStar>=0 && slopeLeft<=0) {
00785 F[0] = F[0] + slopeLeft*(Ustar[0] - vec[0]);
00786 }
00787
00788 else if(slopeM>=0 && slopeLeftStar<=0){
00789 F[0] = F[0] + slopeLeftStar*Ustarstar[0] - (slopeLeftStar - slopeLeft)*Ustar[0] - slopeLeft*vec[0];
00790 }
00791
00792 else if(slopeRightStar>=0 && slopeM<=0){
00793 F[0] = F[1] + slopeRightStar*Ustarstar[1] - (slopeRightStar - slopeRight)*Ustar[1] - slopeRight*vec[1];
00794 }
00795
00796 else if(slopeRight>=0 && slopeRightStar<=0){
00797 F[0] = F[1] + slopeRight*(Ustar[1] - vec[1]);
00798 }
00799
00800 else if(slopeRight<0)F[0] = F[1];
00801 }
00802
00803
00804
00805
00806
00807
00808 void RiemannSolver_HLL(VectorType vec[2], DataType rhoE[2], VectorType F[2], int position, DataType& slopeMax, int BFieldPosition) const
00809 {
00810 DataType slopeLeft, slopeRight;
00811 slope_HLL(vec,position,slopeLeft,slopeRight,slopeMax,BFieldPosition);
00812
00813
00814
00815 vec [0][1] = rhoE[0]; vec [1][1] = rhoE[1];
00816 vec [0][3] = vec[0][0]*vec[0][3]; vec [1][3] = vec[1][0]*vec[1][3];
00817 vec [0][4] = vec[0][0]*vec[0][4]; vec [1][4] = vec[1][0]*vec[1][4];
00818 vec [0][5] = vec[0][0]*vec[0][5]; vec [1][5] = vec[1][0]*vec[1][5];
00819
00820
00821
00822
00823
00824
00825
00826
00827 slopeLeft =std::min(slopeLeft,0.0);
00828 slopeRight=std::max(slopeRight,0.0);
00829
00830 F[0] = ( slopeRight*F[0]-slopeLeft*F[1] + slopeLeft*slopeRight*(vec[1]-vec[0]) )
00831 /(slopeRight-slopeLeft);
00832
00833 }
00834
00835
00836
00837
00838
00839 void computeMagneticAcousticFastWave (const VectorType vec[2], DataType& cfastLeft, DataType& cfastRight, int BFieldPosition) const
00840 {
00841
00842
00843 DataType four = DataType(4.), half = DataType(0.5);
00844 DataType a1l,a2l,a1r,a2r,
00845 rhol,prel,
00846 rhor,prer,
00847 Bxl,Byl,Bzl,Bl2,
00848 Bxr,Byr,Bzr,Br2;
00849
00850 rhol = vec[0][0];
00851 prel = vec[0][1];
00852 Bxl = vec[0][6];
00853 Byl = vec[0][7];
00854 Bzl = vec[0][8];
00855
00856 rhor = vec[1][0];
00857 prer = vec[1][1];
00858 Bxr = vec[1][6];
00859 Byr = vec[1][7];
00860 Bzr = vec[1][8];
00861
00862 Bl2 = Bxl*Bxl + Byl*Byl + Bzl*Bzl;
00863 Br2 = Bxr*Bxr + Byr*Byr + Bzr*Bzr;
00864
00865 a1l = (gamma*prel);
00866 a1r = (gamma*prer);
00867
00868 a2l = four * a1l * vec[0][BFieldPosition] * vec[0][BFieldPosition];
00869 a2r = four * a1r * vec[1][BFieldPosition] * vec[1][BFieldPosition];
00870
00871
00872
00873 cfastLeft = std::sqrt( half/rhol*( a1l + Bl2 + std::sqrt( (a1l + Bl2)*(a1l + Bl2) - a2l) ) );
00874 cfastRight = std::sqrt( half/rhor*( a1r + Br2 + std::sqrt( (a1r + Br2)*(a1r + Br2) - a2r) ) );
00875 }
00876
00877
00878
00879
00880
00881
00882
00883
00884 void computeStateUstarsHLLD(const VectorType vec[2], DataType rhoE[2], VectorType Ustar[2], VectorType Ustarstar[2], int position,
00885 int BFieldPosition, DataType& slopeLeft, DataType& slopeRight, DataType& slopeLeftStar,
00886 DataType& slopeRightStar, DataType& slopeM) const
00887 {
00888 DataType rhol,rhor,prel,prer,psil,psir,vxl,vxr,vyl,vyr,vzl,vzr,Bxr,Bxl,Byl,Byr,Bzl,Bzr,vl,vr,mf,Bl,Br,pTl,pTr,vBl,vBr;
00889 DataType rhols=0.,rhors=0.,vxls=0.,vxrs=0.,vyls=0.,vyrs=0.,vzls=0.,vzrs=0.,Bxls=0.,Bxrs=0.,Byls=0.,Byrs=0.,Bzls=0.,Bzrs=0.,pTs,vBls,vBrs,Els,Ers;
00890 DataType rholss=0.,rhorss=0.,vxlss=0.,vxrss=0.,vylss=0.,vyrss=0.,vzlss=0.,vzrss=0.;
00891 DataType Bxlss=0.,Bxrss=0.,Bylss=0.,Byrss=0.,Bzlss=0.,Bzrss=0.,vBlss=0.,vBrss=0.,Elss=0.,Erss=0.,Bsign;
00892 DataType half=DataType(0.5), auxl=0.,auxr=0.;
00893
00894
00895 rhol = vec[0][0];
00896 prel = vec[0][1];
00897 psil = vec[0][2];
00898 vxl = vec[0][3];
00899 vyl = vec[0][4];
00900 vzl = vec[0][5];
00901 Bxl = vec[0][6];
00902 Byl = vec[0][7];
00903 Bzl = vec[0][8];
00904
00905 rhor = vec[1][0];
00906 prer = vec[1][1];
00907 psir = vec[1][2];
00908 vxr = vec[1][3];
00909 vyr = vec[1][4];
00910 vzr = vec[1][5];
00911 Bxr = vec[1][6];
00912 Byr = vec[1][7];
00913 Bzr = vec[1][8];
00914
00915
00916 vl = vec[0][position];
00917 vr = vec[1][position];
00918
00919
00920 mf = (vec[0][BFieldPosition]+vec[1][BFieldPosition])/DataType(2.0);
00921
00922
00923 Bl = std::sqrt( Bxl*Bxl + Byl*Byl + Bzl*Bzl );
00924 Br = std::sqrt( Bxr*Bxr + Byr*Byr + Bzr*Bzr );
00925
00926 vBl = vxl*Bxl + vyl*Byl + vzl*Bzl;
00927 vBr = vxr*Bxr + vyr*Byr + vzr*Bzr;
00928
00929
00930
00931
00932
00933
00934
00935
00936
00937
00938
00939
00940
00941
00942
00943
00944
00945
00946
00947
00948 pTl = prel + half*Bl*Bl;
00949 pTr = prer + half*Br*Br;
00950
00951
00952 slopeM = ((slopeRight - vr)*rhor*vr - (slopeLeft - vl)*rhol*vl - pTr +pTl)/
00953 ((slopeRight - vr)*rhor - (slopeLeft - vl)*rhol);
00954
00955
00956
00957
00958
00959 rhols = rhol*(slopeLeft-vl)/(slopeLeft-slopeM);
00960 rhors = rhor*(slopeRight-vr)/(slopeRight-slopeM);
00961
00962 if(position == 3 && BFieldPosition == 6){
00963
00964
00965 vxls = slopeM;
00966 vxrs = slopeM;
00967
00968
00969 Bxls = mf;
00970 Bxrs = mf;
00971
00972
00973 auxl= (rhol*(slopeLeft - vl)*(slopeLeft - slopeM)-mf*mf);
00974 auxr= (rhor*(slopeRight - vr)*(slopeRight - slopeM)-mf*mf);
00975
00976 if( std::abs(auxl)<DataType(1.0E-16) ||
00977 ( ( ( std::abs(slopeM -vl) ) < DataType(1.0E-16) ) &&
00978 ( ( std::abs(Byl) + std::abs(Bzl) ) < DataType(1.0E-16) ) &&
00979 ( ( mf*mf ) > (gamma*prel) ) ) ) {
00980
00981
00982
00983
00984 vyls = vyl;
00985 vzls = vzl;
00986 Byls = Byl;
00987 Bzls = Bzl;
00988 }
00989 else{
00990
00991 vyls = vyl - mf*Byl*(slopeM-vl)/auxl;
00992
00993 vzls = vzl - mf*Bzl*(slopeM-vl)/auxl;
00994
00995 Byls = Byl*(rhol*(slopeLeft-vl)*(slopeLeft-vl) - mf*mf)/auxl;
00996
00997 Bzls = Bzl*(rhol*(slopeLeft-vl)*(slopeLeft-vl) - mf*mf)/auxl;
00998 }
00999
01000
01001
01002
01003 if( std::abs(auxr)<DataType(1.0E-16) ||
01004 ( ( ( std::abs(slopeM -vr) ) < DataType(1.0E-16) ) &&
01005 ( ( std::abs(Byr) + std::abs(Bzr) ) < DataType(1.0E-16) ) &&
01006 ( ( mf*mf ) > (gamma*prer) ) ) ) {
01007
01008 vyrs = vyr;
01009 vzrs = vzr;
01010 Byrs = Byr;
01011 Bzrs = Bzr;
01012 }
01013 else{
01014
01015 vyrs = vyr - mf*Byr*(slopeM-vr)/auxr;
01016
01017 vzrs = vzr - mf*Bzr*(slopeM-vr)/auxr;
01018
01019 Byrs = Byr*(rhor*(slopeRight-vr)*(slopeRight-vr) - mf*mf)/auxr;
01020
01021 Bzrs = Bzr*(rhor*(slopeRight-vr)*(slopeRight-vr) - mf*mf)/auxr;
01022 }
01023
01024 }
01025 else if(position == 4 && BFieldPosition == 7){
01026
01027
01028 vyls = slopeM;
01029 vyrs = slopeM;
01030
01031 Byls = mf;
01032 Byrs = mf;
01033
01034 if( std::abs(auxl)<DataType(1.0E-16) ||
01035 ( ( ( std::abs(slopeM -vl) ) < DataType(1.0E-16) ) &&
01036 ( ( std::abs(Bxl) + std::abs(Bzl) ) < DataType(1.0E-16) ) &&
01037 ( ( mf*mf ) > (gamma*prel) ) ) ) {
01038
01039 vxls = vxl;
01040 vzls = vzl;
01041 Bxls = Bxl;
01042 Bzls = Bzl;
01043 }
01044 else{
01045
01046 vxls = vxl - mf*Bxl*(slopeM-vl)/auxl;
01047
01048 vzls = vzl - mf*Bzl*(slopeM-vl)/auxl;
01049
01050 Bxls = Bxl*(rhol*(slopeLeft-vl)*(slopeLeft-vl) - mf*mf)/auxl;
01051
01052 Bzls = Bzl*(rhol*(slopeLeft-vl)*(slopeLeft-vl) - mf*mf)/auxl;
01053 }
01054
01055 if( std::abs(auxr)<DataType(1.0E-16) ||
01056 ( ( ( std::abs(slopeM -vr) ) < DataType(1.0E-16) ) &&
01057 ( ( std::abs(Bxr) + std::abs(Bzr) ) < DataType(1.0E-16) ) &&
01058 ( ( mf*mf ) > (gamma*prer) ) ) ) {
01059
01060
01061 vxrs = vxr;
01062 vzrs = vzr;
01063 Bxrs = Bxr;
01064 Bzrs = Bzr;
01065 }
01066 else{
01067
01068 vxrs = vxr - mf*Bxr*(slopeM-vr)/auxr;
01069
01070 Bxrs = Bxr*(rhor*(slopeRight-vr)*(slopeRight-vr) - mf*mf)/auxr;
01071
01072 vzrs = vzr - mf*Bzr*(slopeM-vr)/auxr;
01073
01074 Bzrs = Bzr*(rhor*(slopeRight-vr)*(slopeRight-vr) - mf*mf)/auxr;
01075 }
01076 }
01077
01078
01079 pTs = ((slopeRight - vr)*rhor*pTl - (slopeLeft - vl)*rhol*pTr + rhol*rhor*(slopeRight-vr)*(slopeLeft-vl)*(vr-vl))/
01080 ((slopeRight - vr)*rhor - (slopeLeft - vl)*rhol);
01081
01082
01083
01084
01085 vBls = vxls*Bxls + vyls*Byls + vzls*Bzls;
01086 vBrs = vxrs*Bxrs + vyrs*Byrs + vzrs*Bzrs;
01087
01088
01089 Els = ((slopeLeft-vl) *rhoE[0] - pTl*vl + pTs*slopeM + mf*(vBl - vBls))/(slopeLeft - slopeM);
01090 Ers = ((slopeRight-vr)*rhoE[1] - pTr*vr + pTs*slopeM + mf*(vBr - vBrs))/(slopeRight - slopeM);
01091
01092
01093
01094
01095 Ustar[0][0] = rhols;
01096 Ustar[0][1] = Els;
01097 Ustar[0][2] = psil;
01098 Ustar[0][3] = rhols*vxls;
01099 Ustar[0][4] = rhols*vyls;
01100 Ustar[0][5] = rhols*vzls;
01101 Ustar[0][6] = Bxls;
01102 Ustar[0][7] = Byls;
01103 Ustar[0][8] = Bzls;
01104
01105 Ustar[1][0] = rhors;
01106 Ustar[1][1] = Ers;
01107 Ustar[1][2] = psir;
01108 Ustar[1][3] = rhors*vxrs;
01109 Ustar[1][4] = rhors*vyrs;
01110 Ustar[1][5] = rhors*vzrs;
01111 Ustar[1][6] = Bxrs;
01112 Ustar[1][7] = Byrs;
01113 Ustar[1][8] = Bzrs;
01114
01115
01116
01117 if(mf > 0) Bsign = 1;
01118 else Bsign = -1;
01119
01120
01121
01122
01123 DataType rholsqrt = std::sqrt(rhols);
01124 DataType rhorsqrt = std::sqrt(rhors);
01125 DataType rhosqrtsum = rholsqrt+rhorsqrt;
01126 DataType rhopsqrt = std::sqrt(rhols*rhors);
01127
01128
01129 if((mf*mf/std::min((Bl*Bl),(Br*Br))) < DataType(1.0E-16))
01130 {
01131 Ustarstar[0] = Ustar[0];
01132 Ustarstar[1] = Ustar[1];
01133 }
01134 else{
01135
01136 rholss = rhols;
01137 rhorss = rhors;
01138 if(position==3 && BFieldPosition==6 ){
01139
01140
01141 vxlss = slopeM;
01142 vxrss = slopeM;
01143
01144 vylss = (rholsqrt*vyls + rhorsqrt*vyrs + (Byrs - Byls)*Bsign)/rhosqrtsum;
01145 vyrss = vylss;
01146
01147
01148
01149 Bxlss = mf;
01150 Bxrss = mf;
01151
01152 Bylss = (rholsqrt*Byrs + rhorsqrt*Byls + rhopsqrt*(vyrs - vyls)*Bsign)/rhosqrtsum;
01153 Byrss = Bylss;
01154 }
01155 else if(position==4 && BFieldPosition==7 ){
01156
01157
01158 vxlss = (rholsqrt*vxls + rhorsqrt*vxrs + (Bxrs - Bxls)*Bsign)/rhosqrtsum;
01159 vxrss = vxlss;
01160
01161 vylss = slopeM;
01162 vyrss = slopeM;
01163
01164
01165
01166 Bxlss = (rholsqrt*Bxrs + rhorsqrt*Bxls + rhopsqrt*(vxrs - vxls)*Bsign)/rhosqrtsum;
01167 Bxrss = Bxlss;
01168
01169 Bylss = mf;
01170 Byrss = mf;
01171 }
01172
01173
01174 vzlss = (rholsqrt*vzls + rhorsqrt*vzrs + (Bzrs - Bzls)*Bsign)/rhosqrtsum;
01175 vzrss = vzlss;
01176
01177 Bzlss = (rholsqrt*Bzrs + rhorsqrt*Bzls + rhopsqrt*(vzrs - vzls)*Bsign)/rhosqrtsum;
01178 Bzrss = Bzlss;
01179 }
01180
01181 vBlss = vxlss*Bxlss + vylss*Bylss + vzlss*Bzlss;
01182 vBrss = vxrss*Bxrss + vyrss*Byrss + vzrss*Bzrss;
01183
01184
01185 Elss = Els - rholsqrt*(vBls - vBlss)*Bsign;
01186 Erss = Ers + rhorsqrt*(vBrs - vBrss)*Bsign;
01187
01188
01189
01190 Ustarstar[0][0] = rholss;
01191 Ustarstar[0][1] = Elss;
01192 Ustarstar[0][2] = psil;
01193 Ustarstar[0][3] = rholss*vxlss;
01194 Ustarstar[0][4] = rholss*vylss;
01195 Ustarstar[0][5] = rholss*vzlss;
01196 Ustarstar[0][6] = Bxlss;
01197 Ustarstar[0][7] = Bylss;
01198 Ustarstar[0][8] = Bzlss;
01199
01200 Ustarstar[1][0] = rhorss;
01201 Ustarstar[1][1] = Erss;
01202 Ustarstar[1][2] = psir;
01203 Ustarstar[1][3] = rhorss*vxrss;
01204 Ustarstar[1][4] = rhorss*vyrss;
01205 Ustarstar[1][5] = rhorss*vzrss;
01206 Ustarstar[1][6] = Bxrss;
01207 Ustarstar[1][7] = Byrss;
01208 Ustarstar[1][8] = Bzrss;
01209
01210
01211 slopeLeftStar = slopeM - std::abs(mf)/rholsqrt;
01212 slopeRightStar = slopeM + std::abs(mf)/rhorsqrt;
01213 }
01214
01215
01216
01217
01218
01219
01220 void slope_HLL (const VectorType vec[2], int position, DataType& slopeLeft, DataType& slopeRight, DataType& slopeMax, int BFieldPosition) const
01221 {
01222 DataType vr,vl,cfastl,cfastr;
01223 computeMagneticAcousticFastWave(vec,cfastl,cfastr,BFieldPosition);
01224
01225 vl = vec[0][position];
01226 vr = vec[1][position];
01227
01228
01229 slopeLeft = std::min(std::min(vl - cfastl,vr - cfastr),DataType(0.));
01230 slopeRight = std::max(std::max(vl + cfastl,vr + cfastr),DataType(0.));
01231
01232 slopeMax = std::max(std::max(std::abs(slopeLeft),std::abs(slopeRight)),slopeMax);
01233 }
01234
01235
01236
01237
01238
01239
01240
01241 void slope_HLLD (const VectorType vec[2], int position, DataType& slopeLeft, DataType& slopeRight, DataType& slopeMax, int BFieldPosition) const
01242 {
01243 DataType vr,vl,cfastl,cfastr;
01244 computeMagneticAcousticFastWave(vec,cfastl,cfastr,BFieldPosition);
01245
01246 vl = vec[0][position];
01247 vr = vec[1][position];
01248
01249
01250 slopeLeft = std::min(std::min(vl,vr) - std::max(cfastl,cfastr),DataType(0.));
01251 slopeRight = std::max(std::max(vl,vr) + std::max(cfastl,cfastr),DataType(0.));
01252 slopeMax = std::max(std::max(std::abs(slopeLeft),std::abs(slopeRight)),slopeMax);
01253 }
01254
01255
01256
01257
01258
01259 void fluxCorrectionX(vec_grid_data_type &vecm, vec_grid_data_type &vec, vec_grid_data_type &flux, DataType ch) const
01260 {
01261 DataType psil,psir,bl,br;
01262 int Nx = vec.extents(0), Ny = vec.extents(1);
01263 VectorType *meshm = (VectorType *)vecm.databuffer();
01264 VectorType *mesh = (VectorType *)vec.databuffer();
01265 VectorType *F = (VectorType *)flux.databuffer();
01266
01267 for(int j=0;j<Ny;j++){
01268 for(int i=1;i<Nx;i++){
01269 register int pos = base::idx(i,j,Nx), posm = base::idx(i-1,j,Nx);
01270
01271 psil = meshm[posm](2);
01272 psir = mesh[pos](2);
01273 bl = meshm[posm](6);
01274 br = mesh[pos](6);
01275
01276 F[pos](6) += psil + DataType(0.5)*(psir-psil)-ch*DataType(0.5)*(br-bl);
01277 F[pos](2) = (bl + DataType(0.5)*(br-bl)-DataType(0.5)*(psir-psil)/ch)*ch*ch;
01278 }
01279 }
01280 }
01281
01282 void fluxCorrectionY(vec_grid_data_type &vecm, vec_grid_data_type &vec, vec_grid_data_type &flux, DataType ch) const
01283 {
01284 DataType psil,psir,bl,br;
01285 int Nx = vec.extents(0), Ny = vec.extents(1);
01286 VectorType *meshm = (VectorType *)vecm.databuffer();
01287 VectorType *mesh = (VectorType *)vec.databuffer();
01288 VectorType *F = (VectorType *)flux.databuffer();
01289
01290 for(int j=1;j<Ny;j++){
01291 for(int i=0;i<Nx;i++){
01292 register int pos = base::idx(i,j,Nx), posm = base::idx(i,j-1,Nx);
01293
01294 psil = meshm[posm](2);
01295 psir = mesh[pos](2);
01296 bl = meshm[posm](7);
01297 br = mesh[pos](7);
01298
01299 F[pos](7) += psil + DataType(0.5)*(psir-psil)-ch*DataType(0.5)*(br-bl);
01300 F[pos](2) = (bl + DataType(0.5)*(br-bl)-DataType(0.5)*(psir-psil)/ch)*ch*ch;
01301 }
01302 }
01303 }
01304
01305
01306
01307
01308 void conservationStep(vec_grid_data_type &vec, grid_data_type &rE, vec_grid_data_type &fluxx,
01309 vec_grid_data_type &fluxy, double dt, double dx, double dy) const
01310 {
01311
01312 primitive2conservative(vec,rE);
01313
01314
01315
01316 DataType dtx = dt/dx;
01317 DataType dty = dt/dy;
01318 int Nx = vec.extents(0), Ny = vec.extents(1);
01319 int mxs = base::NGhosts(), mys = base::NGhosts();
01320 int mxe = Nx-base::NGhosts()-1, mye = Ny-base::NGhosts()-1;
01321 VectorType *mesh = (VectorType *)vec.databuffer();
01322 VectorType *F = (VectorType *)fluxx.databuffer();
01323 VectorType *G = (VectorType *)fluxy.databuffer();
01324
01325 for(int j=mys;j<=mye;j++){
01326 for(int i=mxs;i<=mxe;i++)
01327 mesh[base::idx(i,j,Nx)] -= dtx*(F[base::idx(i+1,j,Nx)]-F[base::idx(i,j,Nx)]) +
01328 dty*(G[base::idx(i,j+1,Nx)]-G[base::idx(i,j,Nx)]);
01329 }
01330 }
01331
01332
01333
01334 void computeCorrection(vec_grid_data_type &vec, grid_data_type &rE, DataType ch, double dt, double dx, double dy) const
01335 {
01336 int Nx = vec.extents(0), Ny = vec.extents(1);
01337 int mxs = base::NGhosts(), mys = base::NGhosts();
01338 int mxe = Nx-base::NGhosts()-1, mye = Ny-base::NGhosts()-1;
01339 VectorType *mesh = (VectorType *)vec.databuffer();
01340 DataType *rhoE = (DataType *)rE.databuffer();
01341 DataType divb,dx2 = dx*2.0,dy2 = dy*2.0,rho,dpsidx,dpsidy,Bx,By,Bz;
01342
01343 for(int i=mxs;i<=mxe;i++){
01344 for(int j=mys;j<=mye;j++){
01345 register int pos = base::idx(i,j,Nx);
01346 register int pxm = base::idx(i-1,j,Nx), pxp = base::idx(i+1,j,Nx);
01347 register int pym = base::idx(i,j-1,Nx), pyp = base::idx(i,j+1,Nx);
01348
01349 dpsidx = (mesh[pxp](2)-mesh[pxm](2))/dx2;
01350 dpsidy = (mesh[pyp](2)-mesh[pym](2))/dy2;
01351 divb = (mesh[pxp](6)-mesh[pxm](6))/dx2 + (mesh[pyp](7)-mesh[pym](7))/dy2;
01352
01353 rho= mesh[pos](0);
01354 Bx = mesh[pos](6);
01355 By = mesh[pos](7);
01356 Bz = mesh[pos](8);
01357 mesh[pos](3) = mesh[pos](3)*rho - dt*Bx*divb;
01358 mesh[pos](4) = mesh[pos](4)*rho - dt*By*divb;
01359 mesh[pos](5) = mesh[pos](5)*rho - dt*Bz*divb;
01360 mesh[pos](2) = mesh[pos](2)*DataType(exp(-ch*dt/cr));
01361 mesh[pos](1) = rhoE[pos] - dt*(Bx*dpsidx + By*dpsidy);
01362 }
01363 }
01364 }
01365
01366 public:
01367
01368 void primitive2conservative(vec_grid_data_type &vec, grid_data_type &rE) const
01369 {
01370 DataType rho;
01371 int Nx = vec.extents(0), Ny = vec.extents(1);
01372 VectorType *mesh = (VectorType *)vec.databuffer();
01373 DataType *rhoE = (DataType *)rE.databuffer();
01374 for(int j=0;j<Ny;j++){
01375 for(int i=0;i<Nx;i++){
01376 register int pos = base::idx(i,j,Nx);
01377 rho = mesh[pos](0);
01378 mesh[pos](3) = mesh[pos](3)*rho;
01379 mesh[pos](4) = mesh[pos](4)*rho;
01380 mesh[pos](5) = mesh[pos](5)*rho;
01381 mesh[pos](1) = rhoE[pos];
01382 }
01383 }
01384 }
01385
01386
01387 void conservative2primitive(vec_grid_data_type &vec, grid_data_type &rE) const
01388 {
01389 DataType rho,vx,vy,vz,Bx,By,Bz;
01390 int Nx = vec.extents(0), Ny = vec.extents(1);
01391 VectorType *mesh = (VectorType *)vec.databuffer();
01392 DataType *rhoE = (DataType *)rE.databuffer();
01393
01394 for(int j=0;j<Ny;j++){
01395 for(int i=0;i<Nx;i++){
01396 register int pos = base::idx(i,j,Nx);
01397 rho = mesh[pos](0);
01398 rhoE[pos] = mesh[pos](1);
01399 vx = mesh[pos](3) = mesh[pos](3)/rho;
01400 vy = mesh[pos](4) = mesh[pos](4)/rho;
01401 vz = mesh[pos](5) = mesh[pos](5)/rho;
01402
01403 vx*=vx;
01404 vy*=vy;
01405 vz*=vz;
01406
01407 Bx = mesh[pos](6);
01408 By = mesh[pos](7);
01409 Bz = mesh[pos](8);
01410
01411 Bx*=Bx;
01412 By*=By;
01413 Bz*=Bz;
01414
01415 mesh[pos](1) = ( mesh[pos](1)-0.5*(vx+vy+vz)*rho - 0.5*(Bx+By+Bz))*(gamma-(double)1.0);
01416 }
01417 }
01418 }
01419
01420
01421
01422 void computeRhoE(vec_grid_data_type &vec, grid_data_type &rE) const
01423 {
01424 DataType rho, pre, vx, vy, vz, vx2, vy2, vz2, v2,
01425 Bx, By, Bz, Bx2, By2, Bz2, B2,
01426 half = DataType(0.5);
01427 int Nx = vec.extents(0), Ny = vec.extents(1);
01428 VectorType *mesh = (VectorType *)vec.databuffer();
01429 DataType *rhoE = (DataType *)rE.databuffer();
01430
01431 for(int j=0;j<Ny;j++)
01432 for(int i=0;i<Nx;i++){
01433 register int pos = base::idx(i,j,Nx);
01434
01435 rho = mesh[pos](0);
01436 pre = mesh[pos](1);
01437 vx = mesh[pos](3);
01438 vy = mesh[pos](4);
01439 vz = mesh[pos](5);
01440 Bx = mesh[pos](6);
01441 By = mesh[pos](7);
01442 Bz = mesh[pos](8);
01443
01444 Bx2 = Bx*Bx; By2 = By*By; Bz2 = Bz*Bz;
01445 B2 = half*(Bz2+Bx2+By2);
01446
01447 vx2 = vx*vx; vy2 = vy*vy; vz2 = vz*vz;
01448 v2 = half*(vz2+vx2+vy2);
01449
01450 rhoE[pos] = rho*v2 + pre/(gamma-1)+B2;
01451 }
01452 }
01453
01454 DataType Gamma() const { return gamma; }
01455
01456
01457
01458
01459
01460 inline void Reconstruct(DataType &r, DataType &r1, DataType &r2, DataType &rl, DataType &rr) const {
01461 DataType q1 = r - r1, q2 = r2 - r;
01462 DataType sl1 = SlopeLimiter(q1,q2)*q1, sl2 = SlopeLimiter(q2,q1)*q2;
01463 rl = r - 0.25*((1.+om)*sl1 + (1.-om)*sl2);
01464 rr = r + 0.25*((1.-om)*sl1 + (1.+om)*sl2);
01465 }
01466
01467 inline DataType SlopeLimiter(DataType &a, DataType &b) const {
01468 if (a==0. || Limiter<=0) return 1.;
01469 DataType r = b/a;
01470 switch(Limiter) {
01471 case 1:
01472
01473 return std::max(0.,std::min(1.,r));
01474 case 2:
01475
01476 return std::max(0.,std::max(std::min(1.,2.*r),std::min(2.,r)));
01477 case 3:
01478
01479 return (r+std::abs(r))/(1.+std::abs(r));
01480 case 4:
01481
01482 return std::max(0.,std::min((1.+r)*0.5,std::min(2.,2*r)));
01483 case 5:
01484
01485 return std::max(0.,r*(1.+r)/(1.+r*r));
01486 default:
01487 return 1.;
01488 }
01489 }
01490
01491 private:
01492 int SchemeNb, no_clean, Limiter;
01493 DataType gamma, epsilon, cr, om;
01494 double cfl_clean;
01495 };
01496
01497
01498
01499
01500
01501
01502
01503
01504
01505
01506
01507
01508
01509
01510
01511
01512
01513
01514
01515
01516
01517
01518
01519
01520
01521
01522
01523
01524
01525
01526
01527
01528
01529
01530
01531
01532
01533
01534
01535
01536
01537
01538
01539
01540
01541
01542
01543
01544
01545
01546
01547
01548
01549
01550
01551
01552
01553
01554
01555
01556
01557
01558
01559
01560
01561
01562
01563
01564
01565
01566
01567
01568
01569
01570
01571
01572
01573
01574
01575
01576
01577
01578
01579
01580
01581
01582
01583
01584
01585
01586
01587
01588
01589
01590
01591
01592
01593
01594
01595
01596
01597
01598
01599
01600
01601
01602
01603
01604
01605
01606
01607
01608
01609
01610
01611
01612
01613
01614
01615
01616
01617
01618
01619
01620
01621
01622
01623
01624
01625
01626
01627
01628
01629
01630
01631
01632
01633
01634
01635
01636
01637
01638
01639
01640
01641
01642
01643
01644
01645
01646
01647
01648
01649
01650
01651
01652
01653
01654 #endif