Athena: prob/strat.c Source File

Go to the documentation of this file.
00001 #include "copyright.h"
00002 /*============================================================================*/
00003 /*! \file strat.c
00004  *  \brief Problem generator for stratified 3D shearing sheet.
00005  *
00006  * PURPOSE:  Problem generator for stratified 3D shearing sheet.  Based on the 
00007  *   initial conditions described in "Three-dimensional Magnetohydrodynamic
00008  *   Simulations of Vertically Stratified Accretion Disks" by Stone, Hawley,
00009  *   Gammie & Balbus.
00010  *
00011  * Several different field configurations and perturbations are possible:
00012  * -  ifield = 1 - Bz=B0sin(kx*x1) field with zero-net-flux [default] (kx input)
00013  * -  ifield = 2 - uniform Bz
00014  * -  ifield = 3 - B=(0,B0cos(kx*x1),B0sin(kx*x1))= zero-net flux w helicity
00015  * -  ifield = 4 - uniform By, but only for |z|<2
00016  *
00017  * - ipert = 1 - random perturbations to P and V [default, used by HGB]
00018  *
00019  * Code must be configured using --enable-shearing-box
00020  *
00021  * REFERENCE: Stone, J., Hawley, J. & Balbus, S. A., ApJ 463, 656-673 (1996)
00022  *            Hawley, J. F. & Balbus, S. A., ApJ 400, 595-609 (1992)          */
00023 /*============================================================================*/
00024 
00025 #include <float.h>
00026 #include <math.h>
00027 
00028 #include <stdlib.h>
00029 #include <string.h>
00030 #include "defs.h"
00031 #include "athena.h"
00032 #include "globals.h"
00033 #include "prototypes.h"
00034 
00035 /*==============================================================================
00036  * PRIVATE FUNCTION PROTOTYPES:
00037  * ran2()           - random number generator from NR
00038  * UnstratifiedDisk() - tidal potential in 3D shearing box
00039  * VertGrav()         - potential for vertical component of gravity
00040  * expr_*()         - computes new output variables
00041  * hst_*            - adds new history variables
00042  *============================================================================*/
00043 
00044 static double ran2(long int *idum);
00045 static Real UnstratifiedDisk(const Real x1, const Real x2, const Real x3);
00046 static Real VertGrav(const Real x1, const Real x2, const Real x3);
00047 static Real expr_dV2(const GridS *pG, const int i, const int j, const int k);
00048 static Real expr_beta(const GridS *pG, const int i, const int j, const int k);
00049 static Real expr_ME(const GridS *pG, const int i, const int j, const int k);
00050 static Real expr_KE(const GridS *pG, const int i, const int j, const int k);
00051 static Real hst_rho_Vx_dVy(const GridS *pG,const int i,const int j,const int k);
00052 static Real hst_rho_dVy2(const GridS *pG, const int i, const int j, const int k);
00053 #ifdef ADIABATIC
00054 static Real hst_E_total(const GridS *pG, const int i, const int j, const int k);
00055 #endif
00056 #ifdef MHD
00057 static Real hst_Bx(const GridS *pG, const int i, const int j, const int k);
00058 static Real hst_By(const GridS *pG, const int i, const int j, const int k);
00059 static Real hst_Bz(const GridS *pG, const int i, const int j, const int k);
00060 static Real hst_BxBy(const GridS *pG, const int i, const int j, const int k);
00061 #endif /* MHD */
00062 
00063 /* top and bottom of root Domain, shared with outputs, etc. */
00064 static Real ztop, zbtm;
00065 
00066 /*=========================== PUBLIC FUNCTIONS =================================
00067  *============================================================================*/
00068 /*----------------------------------------------------------------------------*/
00069 /* problem:  */
00070 
00071 void problem(DomainS *pDomain)
00072 {
00073   GridS *pGrid = (pDomain->Grid);
00074   FILE *fp;
00075   Real xFP[160],dFP[160],vxFP[160],vyFP[160];
00076   int is = pGrid->is, ie = pGrid->ie;
00077   int js = pGrid->js, je = pGrid->je;
00078   int ks = pGrid->ks, ke = pGrid->ke;
00079   int ixs,jxs,kxs,i,j,k,ipert,ifield;
00080   long int iseed = -1; /* Initialize on the first call to ran2 */
00081   Real x1,x2,x3,xmin,xmax,Lx,Ly,Lz;
00082   Real den=1.0, pres=1.0e-6, rd, rp, rvx, rvy, rvz;
00083   Real beta,B0,kx,amp;
00084   int nwx,nwy,nwz;  /* input number of waves per Lx,Ly,Lz [default=1] */
00085   double rval;
00086   static int frst=1;  /* flag so new history variables enrolled only once */
00087 
00088   if (pGrid->Nx[1] == 1){
00089     ath_error("[problem]: HGB only works on a 2D or 3D grid\n");
00090   }
00091 
00092 /* Read problem parameters.  Note Omega set to 10^{-3} by default */
00093 #ifdef ISOTHERMAL
00094   pres=den*Iso_csound2;
00095 #endif
00096   Omega_0 = par_getd_def("problem","omega",1.0e-3);
00097   qshear  = par_getd_def("problem","qshear",1.5);
00098   amp = par_getd("problem","amp");
00099   beta = par_getd("problem","beta");
00100   B0 = sqrt((double)(2.0*pres/beta));
00101   ifield = par_geti_def("problem","ifield", 1);
00102   ipert = par_geti_def("problem","ipert", 1);
00103 
00104 /* Ensure a different initial random seed for each process in an MPI calc. */
00105   ixs = pGrid->Disp[0];
00106   jxs = pGrid->Disp[1];
00107   kxs = pGrid->Disp[2];
00108   iseed = -1 - (ixs + pDomain->Nx[0]*(jxs + pDomain->Nx[1]*kxs));
00109 
00110 /* Initialize boxsize */
00111   ztop = pDomain->RootMaxX[2];
00112   zbtm = pDomain->RootMinX[2];
00113   Lx = pDomain->RootMaxX[0] - pDomain->RootMinX[0];
00114   Ly = pDomain->RootMaxX[1] - pDomain->RootMinX[1];
00115   Lz = pDomain->RootMaxX[2] - pDomain->RootMinX[2];
00116 
00117 /* initialize wavenumbers, given input number of waves per L */
00118   nwx = par_geti_def("problem","nwx",1);
00119   kx = (2.0*PI/Lx)*((double)nwx);
00120 
00121   for (k=ks; k<=ke; k++) {
00122   for (j=js; j<=je; j++) {
00123     for (i=is; i<=ie; i++) {
00124       cc_pos(pGrid,i,j,k,&x1,&x2,&x3);
00125 
00126 /* Initialize perturbations
00127  *  ipert = 1 - random perturbations to P and V [default, used by HGB]
00128  */
00129       if (ipert == 1) {
00130         rval = amp*(ran2(&iseed) - 0.5);
00131         rd = den*exp(-x3*x3);
00132 #ifdef ADIABATIC
00133         rp = pres*rd*(1.0 + 2.0*rval);
00134 #else
00135         rd *= (1.0 + 2.0*rval);
00136 #endif
00137 /* To conform to HGB, the perturbations to V/Cs are (1/5)amp/sqrt(Gamma)  */
00138         rval = amp*(ran2(&iseed) - 0.5);
00139         rvx = 0.4*rval*sqrt(pres/den);
00140 
00141         rval = amp*(ran2(&iseed) - 0.5);
00142         rvy = 0.4*rval*sqrt(pres/den);
00143 
00144         rval = amp*(ran2(&iseed) - 0.5);
00145         rvz = 0.4*rval*sqrt(pres/den);
00146       }
00147 
00148 /* Initialize d, M, and P.  For 3D shearing box M1=Vx, M2=Vy, M3=Vz
00149  * With FARGO do not initialize the background shear */ 
00150 
00151       pGrid->U[k][j][i].d  = rd;
00152       pGrid->U[k][j][i].M1 = rd*rvx;
00153       pGrid->U[k][j][i].M2 = rd*rvy;
00154 #ifndef FARGO
00155       pGrid->U[k][j][i].M2 -= rd*(qshear*Omega_0*x1);
00156 #endif
00157       pGrid->U[k][j][i].M3 = rd*rvz;
00158 #ifdef ADIABATIC
00159       pGrid->U[k][j][i].E = rp/Gamma_1
00160         + 0.5*(SQR(pGrid->U[k][j][i].M1) + SQR(pGrid->U[k][j][i].M2) 
00161              + SQR(pGrid->U[k][j][i].M3))/rd;
00162 #endif
00163 
00164 /* Initialize magnetic field.  For 3D shearing box B1=Bx, B2=By, B3=Bz
00165  *  ifield = 1 - Bz=B0 sin(x1) field with zero-net-flux [default]
00166  *  ifield = 2 - uniform Bz
00167  *  ifield = 3 - B=(0,B0cos(kx*x1),B0sin(kx*x1))= zero-net flux w helicity
00168  *  ifield = 4 - uniform By, but only for |z|<2
00169  */
00170 #ifdef MHD
00171       pGrid->B1i[k][j][i] = 0.0;
00172       pGrid->B2i[k][j][i] = 0.0;
00173       pGrid->B3i[k][j][i] = 0.0;
00174       if (i==ie) pGrid->B1i[k][j][ie+1] = 0.0;
00175       if (j==je) pGrid->B2i[k][je+1][i] = 0.0;
00176       if (k==ke) pGrid->B3i[ke+1][j][i] = 0.0;
00177 
00178       if (ifield == 1) {
00179         pGrid->B1i[k][j][i] = 0.0;
00180         pGrid->B2i[k][j][i] = 0.0;
00181         pGrid->B3i[k][j][i] = B0*(sin((double)kx*x1));
00182         if (i==ie) pGrid->B1i[k][j][ie+1] = 0.0;
00183         if (j==je) pGrid->B2i[k][je+1][i] = 0.0;
00184         if (k==ke) pGrid->B3i[ke+1][j][i] = B0*(sin((double)kx*x1));
00185       }
00186       if (ifield == 2) {
00187         pGrid->B1i[k][j][i] = 0.0;
00188         pGrid->B2i[k][j][i] = 0.0;
00189         pGrid->B3i[k][j][i] = B0;
00190         if (i==ie) pGrid->B1i[k][j][ie+1] = 0.0;
00191         if (j==je) pGrid->B2i[k][je+1][i] = 0.0;
00192         if (k==ke) pGrid->B3i[ke+1][j][i] = B0;
00193       }
00194       if (ifield == 3) {
00195         pGrid->B1i[k][j][i] = 0.0;
00196         pGrid->B2i[k][j][i] = B0*(cos((double)kx*x1));
00197         pGrid->B3i[k][j][i] = B0*(sin((double)kx*x1));
00198         if (i==ie) pGrid->B1i[k][j][ie+1] = 0.0;
00199         if (j==je) pGrid->B2i[k][je+1][i] = B0*(cos((double)kx*x1));
00200         if (k==ke) pGrid->B3i[ke+1][j][i] = B0*(sin((double)kx*x1));
00201       }
00202       if (ifield == 4 && fabs(x3) < 2.0) {
00203         pGrid->B1i[k][j][i] = 0.0;
00204         pGrid->B2i[k][j][i] = B0;
00205         pGrid->B3i[k][j][i] = 0.0;
00206         if (i==ie) pGrid->B1i[k][j][ie+1] = 0.0;
00207         if (j==je) pGrid->B2i[k][je+1][i] = B0;
00208         if (k==ke) pGrid->B3i[ke+1][j][i] = 0.0;
00209       }
00210 
00211 #ifdef ADIABATIC
00212       pGrid->U[k][j][i].E += 0.5*(SQR(pGrid->U[k][j][i].B1c)
00213          + SQR(pGrid->U[k][j][i].B2c) + SQR(pGrid->U[k][j][i].B3c));
00214 #endif
00215 #endif /* MHD */
00216     }
00217   }}
00218 #ifdef MHD
00219   for (k=ks; k<=ke; k++) {
00220     for (j=js; j<=je; j++) {
00221       for (i=is; i<=ie; i++) {
00222         pGrid->U[k][j][i].B1c = 0.5*(pGrid->B1i[k][j][i]+pGrid->B1i[k][j][i+1]);
00223         pGrid->U[k][j][i].B2c = 0.5*(pGrid->B2i[k][j][i]+pGrid->B2i[k][j+1][i]);
00224         pGrid->U[k][j][i].B3c = 0.5*(pGrid->B3i[k][j][i]+pGrid->B3i[k+1][j][i]);
00225       }
00226     }
00227   }
00228 #endif /* MHD */
00229 
00230 /* With viscosity and/or resistivity, read eta_Ohm and nu_V */
00231 #ifdef OHMIC
00232   eta_Ohm = par_getd("problem","eta");
00233 #endif
00234 #ifdef NAVIER_STOKES
00235   nu_V = par_getd("problem","nu");
00236 #endif
00237 
00238 /* enroll gravitational potential function */
00239 
00240   StaticGravPot = VertGrav;
00241   ShearingBoxPot = UnstratifiedDisk;
00242 
00243 /* enroll new history variables */
00244 
00245   if (frst == 1) {
00246     dump_history_enroll(hst_rho_Vx_dVy, "<rho Vx dVy>");
00247     dump_history_enroll(hst_rho_dVy2, "<rho dVy^2>");
00248 #ifdef ADIABATIC
00249     dump_history_enroll(hst_E_total, "<E + rho Phi>");
00250 #endif
00251 #ifdef MHD
00252     dump_history_enroll(hst_Bx, "<Bx>");
00253     dump_history_enroll(hst_By, "<By>");
00254     dump_history_enroll(hst_Bz, "<Bz>");
00255     dump_history_enroll(hst_BxBy, "<-Bx By>");
00256 #endif /* MHD */
00257     frst = 0;
00258   }
00259 
00260   return;
00261 }
00262 
00263 /*==============================================================================
00264  * PUBLIC PROBLEM USER FUNCTIONS:
00265  * problem_write_restart() - writes problem-specific user data to restart files
00266  * problem_read_restart()  - reads problem-specific user data from restart files
00267  * get_usr_expr()          - sets pointer to expression for special output data
00268  * get_usr_out_fun()       - returns a user defined output function pointer
00269  * get_usr_par_prop()      - returns a user defined particle selection function
00270  * Userwork_in_loop        - problem specific work IN     main loop
00271  * Userwork_after_loop     - problem specific work AFTER  main loop
00272  *----------------------------------------------------------------------------*/
00273 
00274 void problem_write_restart(MeshS *pM, FILE *fp)
00275 {
00276   return;
00277 }
00278 
00279 /*
00280  * 'problem_read_restart' must enroll gravity on restarts
00281  */
00282 
00283 void problem_read_restart(MeshS *pM, FILE *fp)
00284 {
00285   Omega_0 = par_getd_def("problem","omega",1.0e-3);
00286   qshear  = par_getd_def("problem","qshear",1.5);
00287 
00288 /* With viscosity and/or resistivity, read eta_Ohm and nu_V */
00289 #ifdef OHMIC
00290   eta_Ohm = par_getd("problem","eta");
00291 #endif
00292 #ifdef NAVIER_STOKES
00293   nu_V = par_getd("problem","nu");
00294 #endif
00295 
00296 /* enroll gravitational potential function */
00297 
00298   StaticGravPot = VertGrav;
00299   ShearingBoxPot = UnstratifiedDisk;
00300 
00301 /* enroll new history variables */
00302 
00303   dump_history_enroll(hst_rho_Vx_dVy, "<rho Vx dVy>");
00304   dump_history_enroll(hst_rho_dVy2, "<rho dVy^2>");
00305 #ifdef ADIABATIC
00306   dump_history_enroll(hst_E_total, "<E + rho Phi>");
00307 #endif
00308 #ifdef MHD
00309   dump_history_enroll(hst_Bx, "<Bx>");
00310   dump_history_enroll(hst_By, "<By>");
00311   dump_history_enroll(hst_Bz, "<Bz>");
00312   dump_history_enroll(hst_BxBy, "<-Bx By>");
00313 #endif /* MHD */
00314 
00315   return;
00316 }
00317 
00318 /* Get_user_expression computes dVy */
00319 ConsFun_t get_usr_expr(const char *expr)
00320 {
00321   if(strcmp(expr,"dVy")==0) return expr_dV2;
00322   else if(strcmp(expr,"beta")==0) return expr_beta;
00323   else if(strcmp(expr,"ME")==0) return expr_ME;
00324   else if(strcmp(expr,"KE")==0) return expr_KE;
00325   else if(strcmp(expr,"BxBy")==0) return hst_BxBy;
00326   return NULL;
00327 }
00328 
00329 VOutFun_t get_usr_out_fun(const char *name){
00330   return NULL;
00331 }
00332 
00333 void Userwork_in_loop(MeshS *pM)
00334 {
00335 }
00336 
00337 void Userwork_after_loop(MeshS *pM)
00338 {
00339 }
00340 
00341 /*------------------------------------------------------------------------------
00342  */
00343 
00344 #define IM1 2147483563
00345 #define IM2 2147483399
00346 #define AM (1.0/IM1)
00347 #define IMM1 (IM1-1)
00348 #define IA1 40014
00349 #define IA2 40692
00350 #define IQ1 53668
00351 #define IQ2 52774
00352 #define IR1 12211
00353 #define IR2 3791
00354 #define NTAB 32
00355 #define NDIV (1+IMM1/NTAB)
00356 #define RNMX (1.0-DBL_EPSILON)
00357 
00358 /*! \fn double ran2(long int *idum)
00359  *  \brief Extracted from the Numerical Recipes in C (version 2) code.  Modified
00360  *   to use doubles instead of floats. -- T. A. Gardiner -- Aug. 12, 2003 
00361  *
00362  * Long period (> 2 x 10^{18}) random number generator of L'Ecuyer
00363  * with Bays-Durham shuffle and added safeguards.  Returns a uniform
00364  * random deviate between 0.0 and 1.0 (exclusive of the endpoint
00365  * values).  Call with idum = a negative integer to initialize;
00366  * thereafter, do not alter idum between successive deviates in a
00367  * sequence.  RNMX should appriximate the largest floating point value
00368  * that is less than 1.
00369  */
00370 
00371 double ran2(long int *idum)
00372 {
00373   int j;
00374   long int k;
00375   static long int idum2=123456789;
00376   static long int iy=0;
00377   static long int iv[NTAB];
00378   double temp;
00379 
00380   if (*idum <= 0) { /* Initialize */
00381     if (-(*idum) < 1) *idum=1; /* Be sure to prevent idum = 0 */
00382     else *idum = -(*idum);
00383     idum2=(*idum);
00384     for (j=NTAB+7;j>=0;j--) { /* Load the shuffle table (after 8 warm-ups) */
00385       k=(*idum)/IQ1;
00386       *idum=IA1*(*idum-k*IQ1)-k*IR1;
00387       if (*idum < 0) *idum += IM1;
00388       if (j < NTAB) iv[j] = *idum;
00389     }
00390     iy=iv[0];
00391   }
00392   k=(*idum)/IQ1;                 /* Start here when not initializing */
00393   *idum=IA1*(*idum-k*IQ1)-k*IR1; /* Compute idum=(IA1*idum) % IM1 without */
00394   if (*idum < 0) *idum += IM1;   /* overflows by Schrage's method */
00395   k=idum2/IQ2;
00396   idum2=IA2*(idum2-k*IQ2)-k*IR2; /* Compute idum2=(IA2*idum) % IM2 likewise */
00397   if (idum2 < 0) idum2 += IM2;
00398   j=(int)(iy/NDIV);              /* Will be in the range 0...NTAB-1 */
00399   iy=iv[j]-idum2;                /* Here idum is shuffled, idum and idum2 */
00400   iv[j] = *idum;                 /* are combined to generate output */
00401   if (iy < 1) iy += IMM1;
00402   if ((temp=AM*iy) > RNMX) return RNMX; /* No endpoint values */
00403   else return temp;
00404 }
00405 
00406 #undef IM1
00407 #undef IM2
00408 #undef AM
00409 #undef IMM1
00410 #undef IA1
00411 #undef IA2
00412 #undef IQ1
00413 #undef IQ2
00414 #undef IR1
00415 #undef IR2
00416 #undef NTAB
00417 #undef NDIV
00418 #undef RNMX
00419 
00420 /*----------------------------------------------------------------------------*/
00421 /*! \fn static Real UnstratifiedDisk(const Real x1, const Real x2,const Real x3)
00422  *  \brief tidal potential in 3D shearing box */
00423 static Real UnstratifiedDisk(const Real x1, const Real x2, const Real x3)
00424 {
00425   Real phi=0.0;
00426 #ifndef FARGO
00427   phi -= qshear*Omega_0*Omega_0*x1*x1;
00428 #endif
00429   return phi;
00430 }
00431 
00432 /*----------------------------------------------------------------------------*/
00433 /*! \fn static Real VertGrav(const Real x1, const Real x2, const Real x3)
00434  *  \brief potential for vertical component of gravity */
00435 static Real VertGrav(const Real x1, const Real x2, const Real x3)
00436 {
00437   Real phi=0.0,z;
00438 
00439   /* Ensure vertical periodicity in ghost zones */
00440   if(x3 > ztop)
00441     z=x3-ztop+zbtm;
00442   else if (x3 < zbtm)
00443     z=x3-zbtm+ztop;
00444   else
00445     z=x3;
00446 
00447   phi += 0.5*Omega_0*Omega_0*
00448    (SQR(fabs(ztop)-sqrt(SQR(fabs(ztop)-fabs(z)) + 0.01)));
00449 
00450   return phi;
00451 }
00452 
00453 /*----------------------------------------------------------------------------*/
00454 /*! \fn static Real expr_dV2(const GridS *pG, const int i, const int j, 
00455  *                           const int k)
00456  *  \brief Computes delta(Vy) 
00457  */
00458 static Real expr_dV2(const GridS *pG, const int i, const int j, const int k)
00459 {
00460   Real x1,x2,x3;
00461   cc_pos(pG,i,j,k,&x1,&x2,&x3);
00462 #ifdef FARGO
00463   return (pG->U[k][j][i].M2/pG->U[k][j][i].d);
00464 #else
00465   return (pG->U[k][j][i].M2/pG->U[k][j][i].d + qshear*Omega_0*x1);
00466 #endif
00467 }
00468 
00469 /*----------------------------------------------------------------------------*/
00470 /*! \fn static Real expr_beta(const GridS *pG, const int i, const int j, 
00471  *                            const int k)
00472  *  \brief Computes beta=P/(B^2/8pi)  
00473  */
00474 static Real expr_beta(const GridS *pG, const int i, const int j, const int k)
00475 {
00476   Real x1,x2,x3,B2;
00477   cc_pos(pG,i,j,k,&x1,&x2,&x3);
00478 #ifdef MHD
00479   B2=pG->U[k][j][i].B1c*pG->U[k][j][i].B1c;
00480   B2+=pG->U[k][j][i].B2c*pG->U[k][j][i].B2c;
00481   B2+=pG->U[k][j][i].B3c*pG->U[k][j][i].B3c;
00482 
00483 #ifdef ISOTHERMAL
00484   return (2.0*Iso_csound2*pG->U[k][j][i].d/B2);
00485 #else
00486   return 0.0;
00487 #endif
00488 
00489 #else
00490   return 0.0;
00491 #endif /* MHD */
00492 }
00493 
00494 /*----------------------------------------------------------------------------*/
00495 /*! \fn static Real expr_ME(const GridS *pG, const int i, const int j, 
00496  *                          const int k)
00497  *  \brief  Computes B^2/8pi
00498  */
00499 static Real expr_ME(const GridS *pG, const int i, const int j, const int k)
00500 {
00501   Real x1,x2,x3,B2;
00502   cc_pos(pG,i,j,k,&x1,&x2,&x3);
00503 #ifdef MHD
00504   B2=pG->U[k][j][i].B1c*pG->U[k][j][i].B1c;
00505   B2+=pG->U[k][j][i].B2c*pG->U[k][j][i].B2c;
00506   B2+=pG->U[k][j][i].B3c*pG->U[k][j][i].B3c;
00507   return (B2/2.0);
00508 #else
00509   return NULL;
00510 #endif
00511 }
00512 /*----------------------------------------------------------------------------*/
00513 /*! \fn static Real expr_KE(const GridS *pG, const int i, const int j, 
00514  *                          const int k)
00515  *  \brief Computes dens*(Vx^2+Vy^2+Vz^2)/2
00516  */
00517 static Real expr_KE(const GridS *pG, const int i, const int j, const int k)
00518 {
00519   Real x1,x2,x3,Vy,Vx,Vz;
00520   cc_pos(pG,i,j,k,&x1,&x2,&x3);
00521 #ifdef FARGO
00522   Vy = (pG->U[k][j][i].M2/pG->U[k][j][i].d);
00523 #else
00524   Vy = (pG->U[k][j][i].M2/pG->U[k][j][i].d + qshear*Omega_0*x1);
00525 #endif
00526   Vx = pG->U[k][j][i].M1/pG->U[k][j][i].d;
00527   Vz = pG->U[k][j][i].M3/pG->U[k][j][i].d;
00528 
00529   return pG->U[k][j][i].d*(Vx*Vx + Vy*Vy + Vz*Vz)/2.0;
00530 
00531 }
00532 
00533 /*------------------------------------------------------------------------------
00534  * Hydro history variables:
00535  * hst_rho_Vx_dVy: Reynolds stress, added as history variable.
00536  * hst_rho_dVy2: KE in y-velocity fluctuations
00537  * hst_E_total: total energy (including tidal potential).
00538  */
00539 /*! \fn static Real hst_rho_Vx_dVy(const GridS *pG,const int i,const int j, 
00540  *                                const int k)
00541  *  \brief Reynolds stress, added as history variable. */
00542 static Real hst_rho_Vx_dVy(const GridS *pG,const int i,const int j, const int k)
00543 {
00544   Real x1,x2,x3;
00545   cc_pos(pG,i,j,k,&x1,&x2,&x3);
00546 #ifdef FARGO
00547   return pG->U[k][j][i].M1*(pG->U[k][j][i].M2/pG->U[k][j][i].d);
00548 #else
00549   return pG->U[k][j][i].M1*
00550     (pG->U[k][j][i].M2/pG->U[k][j][i].d + qshear*Omega_0*x1);
00551 #endif
00552 }
00553 
00554 /*! \fn static Real hst_rho_dVy2(const GridS *pG, const int i, const int j, 
00555  *                              const int k)
00556  *  \brief KE in y-velocity fluctuations */
00557 static Real hst_rho_dVy2(const GridS *pG, const int i, const int j, const int k)
00558 {
00559   Real x1,x2,x3,dVy;
00560   cc_pos(pG,i,j,k,&x1,&x2,&x3);
00561 #ifdef FARGO
00562   dVy = (pG->U[k][j][i].M2/pG->U[k][j][i].d);
00563 #else
00564   dVy = (pG->U[k][j][i].M2/pG->U[k][j][i].d + qshear*Omega_0*x1);
00565 #endif
00566   return pG->U[k][j][i].d*dVy*dVy;
00567 }
00568 
00569 #ifdef ADIABATIC
00570 /*! \fn static Real hst_E_total(const GridS *pG, const int i, const int j, 
00571  *                              const int k)
00572  *  \brief total energy (including tidal potential). */
00573 static Real hst_E_total(const GridS *pG, const int i, const int j, const int k)
00574 {
00575   Real x1,x2,x3,phi;
00576   cc_pos(pG,i,j,k,&x1,&x2,&x3);
00577   phi = UnstratifiedDisk(x1, x2, x3);
00578 
00579   return pG->U[k][j][i].E + pG->U[k][j][i].d*phi;
00580 }
00581 #endif /* ADIABATIC */
00582 
00583 /*------------------------------------------------------------------------------
00584  * MHD history variables
00585  * hst_Bx, etc.: Net flux, and Maxwell stress, added as history variables
00586  */
00587 
00588 #ifdef MHD
00589 /*! \fn static Real hst_Bx(const GridS *pG, const int i,const int j,const int k)
00590  *  \brief x-component of magnetic field */
00591 static Real hst_Bx(const GridS *pG, const int i, const int j, const int k)
00592 {
00593   return pG->U[k][j][i].B1c;
00594 }
00595 
00596 /*! \fn static Real hst_By(const GridS *pG,const int i,const int j,const int k) 
00597  *  \brief y-component of magnetic field */
00598 static Real hst_By(const GridS *pG, const int i, const int j, const int k)
00599 {
00600   return pG->U[k][j][i].B2c;
00601 }
00602 
00603 /*! \fn static Real hst_Bz(const GridS *pG, const int i, const int j, 
00604  *                        const int k)
00605  *  \brief z-component of magnetic field */
00606 static Real hst_Bz(const GridS *pG, const int i, const int j, const int k)
00607 {
00608   return pG->U[k][j][i].B3c;
00609 }
00610 
00611 /*! \fn static Real hst_BxBy(const GridS *pG, const int i, const int j, 
00612  *                           const int k)
00613  *  \brief Maxwell stress */
00614 static Real hst_BxBy(const GridS *pG, const int i, const int j, const int k)
00615 {
00616   return -pG->U[k][j][i].B1c*pG->U[k][j][i].B2c;
00617 }
00618 #endif /* MHD */