amroc/mhd/EGLM2D.h

Go to the documentation of this file.
00001 // -*- C++ -*-
00002 
00003 // Copyright (C) 2013 Oak Ridge National Laboratory
00004 // Ralf Deiterding, deiterdingr@ornl.gov
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     // Divergence-free correction constant
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.);     // x-direction and y-direction slopes
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       // Max slope between x-slope and y-slope
00108       DataType slopeMax = std::max(xSlopeMax/dx,ySlopeMax/dy);
00109       cfl_grid = dt*slopeMax;
00110       
00111       // Time evolution
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       // Divergence-free source terms
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   // Set just the one layer of ghost cells at the boundary 
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     // Wall boundary conditions involve a mirroring of all primitive values,
00336     // then invert normal velocity component
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         // Add boundary velocities if available
00344         if (naux>=2) {
00345           u += aux[naux*n]; 
00346           v += aux[naux*n+1];
00347         } 
00348         // Prescribe only normal velocity vector in ghost cells
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     // Extrapolate values on the boundary itself
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         // Extrapolation of velocity vector in normal direction
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         // Extrapolation of magnetic field vector in normal direction
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   // It computes the energy E, the fluxes Fx and Fy and the conservative variables: density, pression, psi, density*vx, density*vy, density*vz, Bx, By, Bz.
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         // psi = mesh[pos](2);
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   // It computes the energy E, the fluxes Fx and Fy and the conservative variables: density, pression, psi, density*vx, density*vy, density*vz, Bx, By, Bz.
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         // psi = mesh[pos](2);
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   // It computes the fluxes of the VAR conservative variables: density, rho*energy, psi, density*vx, density*vy, density*vz, Bx, By, Bz.
00761   // vec[2][VAR] = conservative variables (left and right).
00762   // F [2][VAR] = numerical flux (left and right) ff.
00763   // Output in F(0)[VAR].
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         // compute U* and U**
00771     computeStateUstarsHLLD (vec,rhoE,Ustar,Ustarstar,position,BFieldPosition,slopeLeft,slopeRight,slopeLeftStar,slopeRightStar,slopeM);
00772     
00773     // We are creating the conservative variables
00774     // vec[](0)-> rho, vec[](1)->rhoE, vec[][3,4,5]-> vx,vy,vz.
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         //Flux Function - Equation 66
00781         //F_L
00782         if(slopeLeft>0)F[0] = F[0];
00783         //F-star left //  FL=FLstar
00784         else if(slopeLeftStar>=0 && slopeLeft<=0) {
00785           F[0] = F[0] + slopeLeft*(Ustar[0] - vec[0]);
00786         }
00787         //F-star-star left
00788         else if(slopeM>=0 && slopeLeftStar<=0){
00789           F[0] = F[0] + slopeLeftStar*Ustarstar[0] - (slopeLeftStar - slopeLeft)*Ustar[0] - slopeLeft*vec[0];
00790         }
00791         //F-star-star right
00792         else if(slopeRightStar>=0 && slopeM<=0){
00793           F[0] = F[1] + slopeRightStar*Ustarstar[1] - (slopeRightStar - slopeRight)*Ustar[1] - slopeRight*vec[1];
00794         }
00795         //F-star right
00796         else if(slopeRight>=0 && slopeRightStar<=0){
00797           F[0] = F[1] + slopeRight*(Ustar[1] - vec[1]);
00798         }
00799         //F_R
00800         else if(slopeRight<0)F[0] = F[1];
00801   }
00802 
00803 //-------------------------------------------------------------------------------------------------------------------------
00804   // It computes the fluxes of the VAR conservative variables: density, rho*energy, psi, density*vx, density*vy, density*vz, Bx, By, Bz.
00805   // vec[2][VAR] = conservative variables (left and right).
00806   // F [2][VAR] = numerical flux (left and right) ff.
00807   // Output in F(0)[VAR].
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     // We are creating the conservative variables
00814     // vec[](0)-> rho, vec[](1)->rhoE, vec[][3,4,5]-> vx,vy,vz.
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     // eq. 10 and 11 could be unified if sl=min(sL,0) and sr=max(sr,0)
00821     /*if(slopeLeft>0)F[0] = F[0];  
00822     else if(slopeRight>=0 && slopeLeft<=0)
00823       F[0] = ( slopeRight*F[0]-slopeLeft*F[1] + slopeLeft*slopeRight*(vec[1]-vec[0]) )
00824             /(slopeRight-slopeLeft);
00825     else if(slopeRight<0) F[0] = F[1];
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   //-------------------------------------------------------------------------------- MD:20130419
00838   // It computes the magnetic-acoustic fast wave. Eq 3, MK:JCP:2005
00839   void computeMagneticAcousticFastWave (const VectorType vec[2], DataType& cfastLeft, DataType& cfastRight, int BFieldPosition) const
00840   {
00841     // BFieldPosition = 6 --> Bx
00842     // BFieldPosition = 7 --> By
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; // notation -> Bl2 = (B_Left)^2
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     //  the magnetic-acoustic fast wave. Eq 3:
00872     // c_f = \sqrt { \frac{\gammap + |\mathbb{B}|^2 + \sqrt{ (\gamma p +|\mathbb{B}|^2)^2 - 4 \gamma p B_x^2 or B_y^2}}{2 \rho} }
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   // It computes the slope. MK:2005  HLLD solver
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     //Variables
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     //v_x and v_y
00916     vl = vec[0][position];
00917     vr = vec[1][position];
00918   
00919     //|B_x| and |B_y|
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    /* //|B|
00930     if(position==3) // direction is x-axis
00931       {
00932         Bl = std::sqrt( mf*mf + Byl*Byl + Bzl*Bzl );
00933         Br = std::sqrt( mf*mf + Byr*Byr + Bzr*Bzr );
00934         //Inner product v.B
00935         vBl = vxl*mf + vyl*Byl + vzl*Bzl;
00936         vBr = vxr*mf + vyr*Byr + vzr*Bzr;
00937       }
00938     else{ //By->mf
00939       Bl = std::sqrt( Bxl*Bxl + mf*mf + Bzl*Bzl );
00940       Br = std::sqrt( Bxr*Bxr + mf*mf + Bzr*Bzr );
00941       //Inner product v.B
00942       vBl = vxl*Bxl + vyl*mf + vzl*Bzl;
00943       vBr = vxr*Bxr + vyr*mf + vzr*Bzr;
00944     }
00945     */
00946     
00947     //Total Pressure eq 2
00948     pTl = prel + half*Bl*Bl;
00949     pTr = prer + half*Br*Br;
00950     
00951     // S_M - Equation 38
00952     slopeM = ((slopeRight - vr)*rhor*vr - (slopeLeft - vl)*rhol*vl - pTr +pTl)/
00953              ((slopeRight - vr)*rhor    - (slopeLeft - vl)*rhol);
00954     
00955     //variables U-star
00956     
00957     //density - Equation 43
00958              
00959     rhols = rhol*(slopeLeft-vl)/(slopeLeft-slopeM);
00960     rhors = rhor*(slopeRight-vr)/(slopeRight-slopeM);
00961     
00962     if(position == 3 && BFieldPosition == 6){
00963       //velocities
00964       //Equation 39
00965       vxls = slopeM;
00966       vxrs = slopeM;
00967       
00968       //B_x
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         //if( std::abs(auxl)<DataType(1.0E-16)){  
00983 //      rhols=rhol;
00984         vyls = vyl;
00985         vzls = vzl;
00986         Byls = Byl;
00987         Bzls = Bzl;
00988       }
00989       else{
00990         //Equation 44
00991         vyls = vyl - mf*Byl*(slopeM-vl)/auxl;
00992         //Equation 46
00993         vzls = vzl - mf*Bzl*(slopeM-vl)/auxl;
00994         //Equation 45
00995         Byls = Byl*(rhol*(slopeLeft-vl)*(slopeLeft-vl) - mf*mf)/auxl;
00996         //Equation 47
00997         Bzls = Bzl*(rhol*(slopeLeft-vl)*(slopeLeft-vl) - mf*mf)/auxl;
00998       }
00999       
01000 //      if( ( std::abs(slopeM -vr) < DataType(1.0E-16) ) || 
01001 //        ( (std::abs(Byr) + std::abs(Bzr)) < DataType(1.0E-16) ) || 
01002 //        ( (mf*mf) >  (gamma*prer) ) ) {
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 //        rhors=rhor;
01008         vyrs = vyr;
01009         vzrs = vzr;
01010         Byrs = Byr;
01011         Bzrs = Bzr;
01012       }
01013       else{
01014         //Equation 44
01015         vyrs = vyr - mf*Byr*(slopeM-vr)/auxr;
01016         //Equation 46
01017         vzrs = vzr - mf*Bzr*(slopeM-vr)/auxr;
01018         //Equation 45;
01019         Byrs = Byr*(rhor*(slopeRight-vr)*(slopeRight-vr) - mf*mf)/auxr;
01020         //Equation 47
01021         Bzrs = Bzr*(rhor*(slopeRight-vr)*(slopeRight-vr) - mf*mf)/auxr;
01022       }
01023       
01024     }
01025     else if(position == 4 && BFieldPosition == 7){
01026       //velocities
01027       //Equation 39
01028       vyls = slopeM;
01029       vyrs = slopeM;
01030       //B_y
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) ) ) ) {  //       if( std::abs(auxl)<DataType(1.0E-16)){
01038 //      rhols=rhol;
01039         vxls = vxl;
01040         vzls = vzl;
01041         Bxls = Bxl;
01042         Bzls = Bzl;
01043       }
01044       else{
01045         //Equation 44
01046         vxls = vxl - mf*Bxl*(slopeM-vl)/auxl;
01047         //Equation 46
01048         vzls = vzl - mf*Bzl*(slopeM-vl)/auxl;
01049         //Equation 45
01050         Bxls = Bxl*(rhol*(slopeLeft-vl)*(slopeLeft-vl) - mf*mf)/auxl;
01051         //Equation 47
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 //      rhors=rhor;
01061         vxrs = vxr;
01062         vzrs = vzr;
01063         Bxrs = Bxr;
01064         Bzrs = Bzr;
01065       }
01066       else{
01067         //Equation 44
01068         vxrs = vxr - mf*Bxr*(slopeM-vr)/auxr;
01069         //Equation 45
01070         Bxrs = Bxr*(rhor*(slopeRight-vr)*(slopeRight-vr) - mf*mf)/auxr;
01071         //Equation 46
01072         vzrs = vzr - mf*Bzr*(slopeM-vr)/auxr;
01073         //Equation 47
01074         Bzrs = Bzr*(rhor*(slopeRight-vr)*(slopeRight-vr) - mf*mf)/auxr;
01075       }
01076     }
01077     
01078      //total pressure - Equation 41
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    // As at this point we have SM= SlopeM,  this eq 41 can be simplfied by  Eq 23, the total pressure must be continuos acrros the middle wave 
01082    // pTs = pTr + rhor*(slopeRight-vr)*(slopeM-vr);
01083 
01084     //inner product vs*Bs
01085     vBls = vxls*Bxls + vyls*Byls + vzls*Bzls;
01086     vBrs = vxrs*Bxrs + vyrs*Byrs + vzrs*Bzrs;
01087     
01088     //energy - Equation 48
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     //U-star variables
01094     //Left
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     //Right
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     //Sign function
01117     if(mf > 0)  Bsign =  1;
01118     else        Bsign = -1;
01119    
01120     
01121     //variables U-star-start
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     //if(  mf*mf/min(Bl*Bl,Br*Br) < epsilon){
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     //density - Equation 49
01136       rholss = rhols;
01137       rhorss = rhors;
01138       if(position==3 && BFieldPosition==6 ){
01139         //velocities
01140         //Equation 39
01141         vxlss = slopeM;
01142         vxrss = slopeM;
01143         //Equation 59
01144         vylss = (rholsqrt*vyls + rhorsqrt*vyrs + (Byrs - Byls)*Bsign)/rhosqrtsum;
01145         vyrss = vylss; // (rholsqrt*vyls + rhorsqrt*vyrs + (Byrs - Byls)*Bsign)/rhosqrtsum;
01146         
01147         //Magnetic Field components
01148         //B_x
01149         Bxlss = mf;
01150         Bxrss = mf;
01151         //Equation 61
01152         Bylss = (rholsqrt*Byrs + rhorsqrt*Byls + rhopsqrt*(vyrs - vyls)*Bsign)/rhosqrtsum;
01153         Byrss = Bylss; // (rholsqrt*Byrs + rhorsqrt*Byls + rhopsqrt*(vyrs - vyls)*Bsign)/rhosqrtsum;
01154       }
01155       else if(position==4 && BFieldPosition==7 ){
01156         //velocities
01157         //Equation 59
01158         vxlss = (rholsqrt*vxls + rhorsqrt*vxrs + (Bxrs - Bxls)*Bsign)/rhosqrtsum;
01159         vxrss =  vxlss; //(rholsqrt*vxls + rhorsqrt*vxrs + (Bxrs - Bxls)*Bsign)/rhosqrtsum;
01160         //Equation 39
01161         vylss = slopeM;
01162         vyrss = slopeM;
01163         
01164         //Magnetic Field components
01165         //Equation 61
01166         Bxlss = (rholsqrt*Bxrs + rhorsqrt*Bxls + rhopsqrt*(vxrs - vxls)*Bsign)/rhosqrtsum;
01167         Bxrss = Bxlss; // (rholsqrt*Bxrs + rhorsqrt*Bxls + rhopsqrt*(vxrs - vxls)*Bsign)/rhosqrtsum;
01168         //B_y
01169         Bylss = mf;
01170         Byrss = mf;
01171       }
01172 
01173       //Equation 60
01174       vzlss = (rholsqrt*vzls + rhorsqrt*vzrs + (Bzrs - Bzls)*Bsign)/rhosqrtsum;
01175       vzrss = vzlss; //(rholsqrt*vzls + rhorsqrt*vzrs + (Bzrs - Bzls)*Bsign)/rhosqrtsum;
01176       //Equation 62
01177       Bzlss = (rholsqrt*Bzrs + rhorsqrt*Bzls + rhopsqrt*(vzrs - vzls)*Bsign)/rhosqrtsum;
01178       Bzrss = Bzlss; //(rholsqrt*Bzrs + rhorsqrt*Bzls + rhopsqrt*(vzrs - vzls)*Bsign)/rhosqrtsum;
01179     }
01180     //inner product vss*Bss
01181     vBlss = vxlss*Bxlss + vylss*Bylss + vzlss*Bzlss;
01182     vBrss = vxrss*Bxrss + vyrss*Byrss + vzrss*Bzrss;
01183   
01184     //Enery - Equation 63
01185     Elss = Els - rholsqrt*(vBls - vBlss)*Bsign;
01186     Erss = Ers + rhorsqrt*(vBrs - vBrss)*Bsign;
01187   
01188     //U-star-star variables
01189     //left
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     //right
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     //Equation 51 - S-star left and S-star right
01211     slopeLeftStar  = slopeM - std::abs(mf)/rholsqrt;
01212     slopeRightStar = slopeM + std::abs(mf)/rhorsqrt;
01213   }
01214     
01215   
01216   //---------------------------------------------------------------------------
01217   // It computes the slope.
01218   // Velocity component in position = 3 for slope in x direction.
01219   // Velocity component in position = 4 for slope in y direction.
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     // Davis limiter - smalest and larger eingenvalues Eq 12 fromMK:JCP:2005 HLL
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   // It computes the slope.
01239   // Velocity component in position = 3 for slope in x direction.
01240   // Velocity component in position = 4 for slope in y direction.
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     //Eq. 67 from MK:JCP:2005 HLLD
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   //Conservation step.
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     // First step: primitive -> conservative form.
01312     primitive2conservative(vec,rE);
01313     
01314     // Second step: time evolution.
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   // It computes rhoE.
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;  // Eq after Eq 2, page 317 MK:2005
01451       }
01452   }
01453 
01454   DataType Gamma() const { return gamma; }
01455 
01456   // MUSCL reconstruction with slope-limiting
01457   // Linear reconstruction: om=0.d0, 2nd order spatial reconstuction: om!=0.d0 
01458   // 3rd order spatial reconstruction: om=1.d0/3.d0 
01459   // Higher order reconstructions are not conservative!
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       // Minmod
01473       return std::max(0.,std::min(1.,r));
01474     case 2:
01475       // Superbee
01476       return std::max(0.,std::max(std::min(1.,2.*r),std::min(2.,r)));
01477     case 3:
01478       // Van Leer
01479       return (r+std::abs(r))/(1.+std::abs(r));
01480     case 4:
01481       // Monotonized centered
01482       return std::max(0.,std::min((1.+r)*0.5,std::min(2.,2*r)));
01483     case 5:
01484       // Van Albada
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 /*OLD FUNCTIONS to be deleted
01498 
01499   //-------------------------------------------------------------------------------------------------------------------------
01500   // It computes the fluxes of the VAR conservative variables: density, rho*energy, psi, density*vx, density*vy, density*vz, Bx, By, Bz.
01501   // vec[2][VAR] = conservative variables (left and right).
01502   // F [2][VAR] = numerical flux (left and right) ff.
01503   // Output in F(0)[VAR].
01504   void RiemannSolver(VectorType vec[2], DataType rhoE[2], VectorType F[2], int position, DataType& slopeMax, int BFieldPosition) const
01505   {
01506     DataType slopeLeft, slopeRight,slopeLeftStar,slopeRightStar,slopeM;
01507     VectorType Ustar[2],Ustarstar[2];
01508     
01509     slope (vec,position,slopeLeft,slopeRight,slopeMax,BFieldPosition);
01510         // compute U* and U**
01511     computeStateUstarsHLLD (vec,rhoE,Ustar,Ustarstar,position,BFieldPosition,slopeLeft,slopeRight,slopeLeftStar,slopeRightStar,slopeM);
01512     
01513     // We are creating the conservative variables
01514     // vec[](0)-> rho, vec[](1)->rhoE, vec[][3,4,5]-> vx,vy,vz.
01515     vec [0][1] = rhoE[0];                       vec [1][1] = rhoE[1];
01516     vec [0][3] = vec[0][0]*vec[0][3];   vec [1][3] = vec[1][0]*vec[1][3];
01517     vec [0][4] = vec[0][0]*vec[0][4];   vec [1][4] = vec[1][0]*vec[1][4];
01518     vec [0][5] = vec[0][0]*vec[0][5];   vec [1][5] = vec[1][0]*vec[1][5];
01519     switch(SchemeNb)
01520       {
01521       case 1: // eq. 10 and 11 could be unified if sl=min(sL,0) and sr=max(sr,0)
01522         if(slopeLeft>0)F[0] = F[0];
01523         if(slopeRight>=0 && slopeLeft<=0)
01524           F[0] = ( slopeRight*F[0]-slopeLeft*F[1] + 
01525                    slopeLeft*slopeRight*(vec[1]-vec[0]) )
01526             /(slopeRight-slopeLeft);
01527         if(slopeRight<0)F[0] = F[1];
01528         break;
01529       
01530       case 2:     
01531         //Flux Function - Equation 66
01532         //F_L
01533         if(slopeLeft>0)F[0] = F[0];
01534         //F-star left //  FL=FLstar
01535         else if(slopeLeftStar>=0 && slopeLeft<=0) {
01536           F[0] = F[0] + slopeLeft*(Ustar[0] - vec[0]);
01537         }
01538         //F-star-star left
01539         else if(slopeM>=0 && slopeLeftStar<=0){
01540           F[0] = F[0] + slopeLeftStar*Ustarstar[0] - (slopeLeftStar - slopeLeft)*Ustar[0] - slopeLeft*vec[0];
01541         }
01542         //F-star-star right
01543         else if(slopeRightStar>=0 && slopeM<=0){
01544           F[0] = F[1] + slopeRightStar*Ustarstar[1] - (slopeRightStar - slopeRight)*Ustar[1] - slopeRight*vec[1];
01545         }
01546         //F-star right
01547         else if(slopeRight>=0 && slopeRightStar<=0){
01548           F[0] = F[1] + slopeRight*(Ustar[1] - vec[1]);
01549         }
01550         //F_R
01551         else if(slopeRight<0)F[0] = F[1];
01552         break;
01553       }
01554   }
01555 
01556   //---------------------------------------------------------------------------
01557   // It computes the slope.
01558   // Velocity component in position = 3 for slope in x direction.
01559   // Velocity component in position = 4 for slope in y direction.
01560   void slope (const VectorType vec[2], int position, DataType& slopeLeft, DataType& slopeRight, DataType& slopeMax, int BFieldPosition) const
01561   {
01562     DataType vr,vl,cfastl,cfastr;
01563     computeMagneticAcousticFastWave(vec,cfastl,cfastr,BFieldPosition);
01564     
01565     vl = vec[0][position];
01566     vr = vec[1][position];
01567     
01568     if(SchemeNb==1) {// Davis limiter - smalest and larger eingenvalues Eq 12 fromMK:JCP:2005 HLL
01569       slopeLeft  = std::min(std::min(vl - cfastl,vr - cfastr),DataType(0.));
01570       slopeRight = std::max(std::max(vl + cfastl,vr + cfastr),DataType(0.));
01571     }
01572     else if(SchemeNb==2) {//Eq. 67 from MK:JCP:2005 HLLD
01573       slopeLeft  = std::min(std::min(vl,vr) - std::max(cfastl,cfastr),DataType(0.)); 
01574       slopeRight = std::max(std::max(vl,vr) + std::max(cfastl,cfastr),DataType(0.));
01575     }
01576     
01577     slopeMax = std::max(std::max(std::abs(slopeLeft),std::abs(slopeRight)),slopeMax);
01578   }
01579 
01580   
01581   
01582   
01583   
01584   
01585   
01586   //--------------------------------------------------------------------------------
01587   // It computes the magnetic-acoustic fast wave
01588   
01589   void eigenvalue (const VectorType vec[2], DataType& eigenvalueleft, DataType& eigenvalueright, int BFieldPosition) const
01590   {
01591     // BFieldPosition = 6 --> Bx
01592     // BFieldPosition = 7 --> By
01593     DataType four = DataType(4.), half = DataType(0.5);
01594     DataType a1l,a2l,a3l,a1r,a2r,a3r,rhol,prel,//psil,
01595       Bxl,Byl,Bzl,Bl,
01596       //vxl,vyl,vzl,
01597       rhor,prer,//psir,
01598       Bxr,Byr,Bzr,Br;//,vxr,vyr,vzr;
01599     
01600     rhol = vec[0][0];
01601     prel = vec[0][1];
01602     //psil = vec(0)[2];
01603     //vxl  = vec(0)[3];
01604     //vyl  = vec(0)[4];
01605     //vzl  = vec(0)[5];
01606     Bxl  = std::abs(vec[0][6]);
01607     Byl  = std::abs(vec[0][7]);
01608     Bzl  = vec[0][8];
01609                 
01610     rhor = vec[1][0];
01611     prer = vec[1][1];
01612     //psir = vec(1)[2];
01613     //vxr  = vec(1)[3];
01614     //vyr  = vec(1)[4];
01615     //vzr  = vec(1)[5];
01616     Bxr  = std::abs(vec[1][6]);
01617     Byr  = std::abs(vec[1][7]);
01618     Bzr  = vec[1][8];
01619     
01620     a1l = std::sqrt( (gamma*prel)/rhol );
01621     a1r = std::sqrt( (gamma*prer)/rhor );
01622     
01623     if(SchemeNb==2){
01624       if(BFieldPosition==6){
01625         Bl = std::sqrt( Bxl*Bxl + Byl*Byl + Bzl*Bzl );
01626         Br = std::sqrt( Bxl*Bxl + Byr*Byr + Bzr*Bzr );
01627       }
01628       else{
01629         Bl = std::sqrt( Bxl*Bxl + Byl*Byl + Bzl*Bzl );
01630         Br = std::sqrt( Bxr*Bxr + Byl*Byl + Bzr*Bzr );
01631       }
01632       a2l = std::sqrt( (std::abs(vec[0][BFieldPosition])*std::abs(vec[0][BFieldPosition]))/rhol );
01633       a2r = std::sqrt( (std::abs(vec[1][BFieldPosition])*std::abs(vec[1][BFieldPosition]))/rhor );  //MD:20130419 here vec[0] before I chenged for 1. Ok this changes.
01634     }
01635     else{  
01636       Bl = std::sqrt( Bxl*Bxl + Byl*Byl + Bzl*Bzl );
01637       Br = std::sqrt( Bxr*Bxr + Byr*Byr + Bzr*Bzr );
01638       a2l = std::sqrt( (std::abs(vec[0][BFieldPosition])*std::abs(vec[0][BFieldPosition]))/rhol );
01639       a2r = std::sqrt( (std::abs(vec[1][BFieldPosition])*std::abs(vec[1][BFieldPosition]))/rhor );
01640     }
01641     
01642     a3l = std::sqrt( (Bl*Bl)/rhol );
01643     a3r = std::sqrt( (Br*Br)/rhor );
01644     
01645     a1l*=a1l;   a2l*=a2l;       a3l*=a3l;
01646     a1r*=a1r;   a2r*=a2r;       a3r*=a3r;
01647     eigenvalueleft  = std::sqrt( half*( a1l + a3l + std::sqrt( (a1l + a3l)*(a1l + a3l) - four*a1l*a2l) ) );
01648     eigenvalueright = std::sqrt( half*( a1r + a3r + std::sqrt( (a1r + a3r)*(a1r + a3r) - four*a1r*a2r) ) );
01649   }
01650 */
01651   
01652 //.........................................................................................
01653 
01654 #endif