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