Athena: prob/hgb.c Source File

Go to the documentation of this file.
00001 #include "copyright.h"
00002 /*============================================================================*/
00003 /*! \file hgb.c
00004  *  \brief Problem generator for 3D shearing sheet.
00005  *
00006  * PURPOSE:  Problem generator for 3D shearing sheet.  Based on the initial
00007  *   conditions described in "Local Three-dimensional Magnetohydrodynamic
00008  *   Simulations of Accretion Disks" by Hawley, Gammie & Balbus, or HGB.
00009  *
00010  * Several different field configurations and perturbations are possible:
00011  *
00012  *- ifield = 0 - uses field set by choice of ipert flag
00013  *- ifield = 1 - Bz=B0sin(kx*x1) field with zero-net-flux [default] (kx input)
00014  *- ifield = 2 - uniform Bz
00015  *- ifield = 3 - B=(0,B0cos(kx*x1),B0sin(kx*x1))= zero-net flux w helicity
00016  *- ifield = 4 - B=(0,B0/sqrt(2),B0/sqrt(2))= net toroidal+vertical field
00017  *- ifield = 5 - uniform By
00018  *
00019  *- ipert = 1 - random perturbations to P and V [default, used by HGB]
00020  *- ipert = 2 - uniform Vx=amp (epicyclic wave test)
00021  *- ipert = 3 - J&G vortical shwave (hydro test)
00022  *- ipert = 4 - nonlinear density wave test of Fromang & Papaloizou
00023  *- ipert = 5 - 2nd MHD shwave test of JGG (2008) -- their figure 9
00024  *- ipert = 6 - 3rd MHD shwave test of JGG (2008) -- their figure 11
00025  *- ipert = 7 - nonlinear shearing wave test of Heinemann & Papaloizou (2008)
00026  *
00027  * To run simulations of stratified disks (including vertical gravity), use the
00028  * strat.c problem generator.
00029  *
00030  * Code must be configured using --enable-shearing-box
00031  *
00032  * REFERENCE: Hawley, J. F. & Balbus, S. A., ApJ 400, 595-609 (1992).
00033  *            Johnson, Guan, & Gammie, ApJSupp, (2008)                        */
00034 /*============================================================================*/
00035 
00036 #include <float.h>
00037 #include <math.h>
00038 
00039 #include <stdlib.h>
00040 #include <string.h>
00041 #include "defs.h"
00042 #include "athena.h"
00043 #include "globals.h"
00044 #include "prototypes.h"
00045 
00046 Real Lx,Ly,Lz; /* root grid size, global to share with output functions */
00047 
00048 /*==============================================================================
00049  * PRIVATE FUNCTION PROTOTYPES:
00050  * ran2()          - random number generator from NR
00051  * UnstratifiedDisk() - tidal potential in 3D shearing box
00052  * expr_dV2()       - computes delta(Vy)
00053  * hst_*            - new history variables
00054  *============================================================================*/
00055 
00056 static double ran2(long int *idum);
00057 static Real UnstratifiedDisk(const Real x1, const Real x2, const Real x3);
00058 static Real expr_dV2(const GridS *pG, const int i, const int j, const int k);
00059 static Real expr_Jsq(const GridS *pG, const int i, const int j, const int k);
00060 static Real hst_rho_Vx_dVy(const GridS *pG,const int i,const int j,const int k);
00061 static Real hst_rho_dVy2(const GridS *pG,const int i, const int j, const int k);
00062 #ifdef ADIABATIC
00063 static Real hst_E_total(const GridS *pG, const int i, const int j, const int k);
00064 #endif
00065 #ifdef MHD
00066 static Real hst_Bx(const GridS *pG, const int i, const int j, const int k);
00067 static Real hst_By(const GridS *pG, const int i, const int j, const int k);
00068 static Real hst_Bz(const GridS *pG, const int i, const int j, const int k);
00069 static Real hst_BxBy(const GridS *pG, const int i, const int j, const int k);
00070 static Real hst_dEw2(const GridS *pG, const int i, const int j, const int k);
00071 static Real hst_dBy(const GridS *pG, const int i, const int j, const int k);
00072 #endif /* MHD */
00073 
00074 /*=========================== PUBLIC FUNCTIONS =================================
00075  *============================================================================*/
00076 /*----------------------------------------------------------------------------*/
00077 /* problem:  */
00078 
00079 void problem(DomainS *pDomain)
00080 {
00081   GridS *pGrid = pDomain->Grid;
00082   FILE *fp;
00083   Real xFP[160],dFP[160],vxFP[160],vyFP[160];
00084   int is = pGrid->is, ie = pGrid->ie;
00085   int js = pGrid->js, je = pGrid->je;
00086   int ks = pGrid->ks, ke = pGrid->ke;
00087   int ixs,jxs,kxs,i,j,k,ipert,ifield;
00088   long int iseed = -1; /* Initialize on the first call to ran2 */
00089   Real x1,x2,x3,xmin,xmax;
00090   Real den = 1.0, pres = 1.0, rd, rp, rvx, rvy, rvz, rbx, rby, rbz;
00091   Real beta=1.0,B0,kx,ky,kz,amp;
00092   int nwx,nwy,nwz;  /* input number of waves per Lx,Ly,Lz [default=1] */
00093   double rval;
00094   static int frst=1;  /* flag so new history variables enrolled only once */
00095 
00096   if (pGrid->Nx[1] == 1){
00097     ath_error("[problem]: HGB only works on a 2D or 3D grid\n");
00098   }
00099 
00100 /* Read problem parameters.  Note Omega_0 set to 10^{-3} by default */
00101   Omega_0 = par_getd_def("problem","Omega",1.0e-3);
00102   qshear  = par_getd_def("problem","qshear",1.5);
00103   amp = par_getd("problem","amp");
00104 #ifdef MHD
00105   beta = par_getd("problem","beta");
00106   ifield = par_geti_def("problem","ifield", 1);
00107 #endif
00108   ipert = par_geti_def("problem","ipert", 1);
00109 
00110 /* Compute field strength based on beta.  */
00111 #ifdef ISOTHERMAL
00112   pres = Iso_csound2;
00113 #else
00114   pres = par_getd("problem","pres");
00115 #endif
00116   B0 = sqrt((double)(2.0*pres/beta));
00117 
00118 /* Ensure a different initial random seed for each process in an MPI calc. */
00119   ixs = pGrid->Disp[0];
00120   jxs = pGrid->Disp[1];
00121   kxs = pGrid->Disp[2];
00122   iseed = -1 - (ixs + pDomain->Nx[0]*(jxs + pDomain->Nx[1]*kxs));
00123 
00124 /* Initialize boxsize */
00125   Lx = pDomain->RootMaxX[0] - pDomain->RootMinX[0];
00126   Ly = pDomain->RootMaxX[1] - pDomain->RootMinX[1];
00127   Lz = pDomain->RootMaxX[2] - pDomain->RootMinX[2];
00128 
00129 /* initialize wavenumbers, given input number of waves per L */
00130   nwx = par_geti_def("problem","nwx",1);
00131   nwy = par_geti_def("problem","nwy",1);
00132   nwz = par_geti_def("problem","nwz",1);
00133   kx = (2.0*PI/Lx)*((double)nwx);  /* nxw should be -ve for leading wave */
00134   ky = (2.0*PI/Ly)*((double)nwy);
00135   kz = (2.0*PI/Lz)*((double)nwz);
00136 
00137 /* For PF density wave test, read data from file */
00138 
00139   if (ipert == 4) {
00140     if (pGrid->Nx[0] == 160) {
00141       if((fp = fopen("Data-160-FPwave.dat","r")) == NULL)
00142          ath_error("Error opening Data-160-FPwave.dat\n");
00143       for (i=0; i<160; i++) {
00144         fscanf(fp,"%lf %lf %lf %lf",&xFP[i],&dFP[i],&vxFP[i],&vyFP[i]);
00145       }
00146     }
00147 
00148     if (pGrid->Nx[0] == 40) {
00149       if((fp = fopen("Data-40-FPwave.dat","r")) == NULL)
00150          ath_error("Error opening Data-40-FPwave.dat\n");
00151       for (i=0; i<40; i++) {
00152         fscanf(fp,"%lf %lf %lf %lf",&xFP[i],&dFP[i],&vxFP[i],&vyFP[i]);
00153       }
00154     }
00155 
00156     xmin = pDomain->RootMinX[0];
00157     if (xmin != -4.7965) ath_error("[hgb]: iprob=4 requires xmin=-4.7965\n");
00158     xmax = pDomain->RootMaxX[0];
00159     if (xmax != 4.7965) ath_error("[hgb]: iprob=4 requires xmax=4.7965\n");
00160   }
00161 
00162 /* Rescale amp to sound speed for ipert 2,3 */
00163 #ifdef ADIABATIC
00164   if (ipert == 2 || ipert == 3) amp *= sqrt(Gamma*pres/den);
00165 #else
00166   if (ipert == 2 || ipert == 3) amp *= Iso_csound;
00167 #endif
00168 
00169   for (k=ks; k<=ke; k++) {
00170   for (j=js; j<=je; j++) {
00171     for (i=is; i<=ie; i++) {
00172       cc_pos(pGrid,i,j,k,&x1,&x2,&x3);
00173 
00174 /* Initialize perturbations
00175  *  ipert = 1 - random perturbations to P and V [default, used by HGB]
00176  *  ipert = 2 - uniform Vx=amp (epicyclic wave test)
00177  *  ipert = 3 - vortical shwave (hydro test)
00178  *  ipert = 4 - Fromang & Papaloizou nonlinear density wave (hydro test)
00179  *  ipert = 5 & 6 - JGG MHD shwave tests
00180  *  ipert = 7 - Heinemann & Papaloizou (2008) nonlinear shwave (hydro test)
00181  */
00182       if (ipert == 1) {
00183         rval = amp*(ran2(&iseed) - 0.5);
00184 #ifdef ADIABATIC
00185         rp = pres*(1.0 + 2.0*rval);
00186         rd = den;
00187 #else
00188         rd = den*(1.0 + 2.0*rval);
00189 #endif
00190 /* To conform to HGB, the perturbations to V/Cs are (1/5)amp/sqrt(Gamma)  */
00191         rval = amp*(ran2(&iseed) - 0.5);
00192         rvx = 0.4*rval*sqrt(pres/den);
00193 
00194         rval = amp*(ran2(&iseed) - 0.5);
00195         rvy = 0.4*rval*sqrt(pres/den);
00196 
00197         rval = amp*(ran2(&iseed) - 0.5);
00198         rvz = 0.4*rval*sqrt(pres/den);
00199       }
00200       if (ipert == 2) {
00201         rp = pres;
00202         rd = den;
00203         rvx = amp;
00204         rvy = 0.0;
00205         rvz = 0.0;
00206       }
00207       if (ipert == 3) {
00208         rp = pres;
00209         rd = den;
00210         rvx = amp*sin((double)(kx*x1 + ky*x2));
00211         rvy = -amp*(kx/ky)*sin((double)(kx*x1 + ky*x2));
00212         rvz = 0.0;
00213       }
00214       if (ipert == 4) {
00215         rd = dFP[i-is];
00216         rvx = vxFP[i-is];
00217         rvy = vyFP[i-is] + qshear*Omega_0*x1; /*subtract mean vy*/
00218         rvz = 0.0;
00219       }
00220 /* Note the initial conditions in JGG for this test are incorrect. */
00221       if (ipert == 5) {
00222         ifield = 0;
00223         rd = den + 8.9525e-10*cos((double)(kx*x1 + ky*x2 + kz*x3 - PI/4.));
00224         rvx = 8.16589e-8*cos((double)(kx*x1 + ky*x2 + kz*x3 + PI/4.));
00225         rvy = 8.70641e-8*cos((double)(kx*x1 + ky*x2 + kz*x3 + PI/4.));
00226         rvz = 0.762537e-8*cos((double)(kx*x1 + ky*x2 + kz*x3 + PI/4.));
00227         rbx = -1.08076e-7;
00228         rbx *= cos((double)(kx*(x1-0.5*pGrid->dx1) + ky*x2 + kz*x3 - PI/4.));
00229         rby = 1.04172e-7;
00230         rby *= cos((double)(kx*x1 + ky*(x2-0.5*pGrid->dx2) + kz*x3 - PI/4.));
00231         rbz = -0.320324e-7;
00232         rbz *= cos((double)(kx*x1 + ky*x2 + kz*(x3-0.5*pGrid->dx3) - PI/4.));;
00233         rbz += (sqrt(15.0)/16.0)*(Omega_0/kz);
00234       } 
00235       if (ipert == 6) {
00236         ifield = 0;
00237         rd = den + 5.48082e-6*cos((double)(kx*x1 + ky*x2 + kz*x3));
00238         rvx = -4.5856e-6*cos((double)(kx*x1 + ky*x2 + kz*x3));
00239         rvy = 2.29279e-6*cos((double)(kx*x1 + ky*x2 + kz*x3));
00240         rvz = 2.29279e-6*cos((double)(kx*x1 + ky*x2 + kz*x3));
00241         rbx = 5.48082e-7;
00242         rbx *= cos((double)(kx*(x1-0.5*pGrid->dx1) + ky*x2 + kz*x3));
00243         rbx += (0.1);
00244         rby = 1.0962e-6;
00245         rby *= cos((double)(kx*x1 + ky*(x2-0.5*pGrid->dx2) + kz*x3));
00246         rby += (0.2);
00247         rbz = 0.0;
00248       }
00249       if (ipert == 7) {
00250 #ifdef ISOTHERMAL
00251         double kappa2 = 2.0*(2.0 - qshear)*Omega_0*Omega_0;
00252         double aa = (kx*kx + ky*ky)*Iso_csound*Iso_csound + kappa2;
00253         double bb = 2.0*qshear*Omega_0*ky*Iso_csound;
00254         double denom = aa*aa + bb*bb;
00255         double rd_hat =            (ky*Iso_csound*bb -2.0*Omega_0*aa)*amp/denom;
00256         double px_hat =-Iso_csound*(ky*Iso_csound*aa +2.0*Omega_0*bb)*amp/denom;
00257         double py_hat = (amp + ky*px_hat + (2.0-qshear)*Omega_0*rd_hat)/kx;
00258         rd  = 1.0 + rd_hat*cos((double)(kx*x1 + ky*x2));
00259         rvx =       px_hat*sin((double)(kx*x1 + ky*x2))/rd;
00260         rvy =       py_hat*sin((double)(kx*x1 + ky*x2))/rd;
00261 #endif
00262         rvz = 0.0;
00263       }
00264 
00265 /* Initialize d, M, and P.  For 3D shearing box M1=Vx, M2=Vy, M3=Vz
00266  * With FARGO do not initialize the background shear */ 
00267 
00268       pGrid->U[k][j][i].d  = rd;
00269       pGrid->U[k][j][i].M1 = rd*rvx;
00270       pGrid->U[k][j][i].M2 = rd*rvy;
00271 #ifndef FARGO
00272       pGrid->U[k][j][i].M2 -= rd*(qshear*Omega_0*x1);
00273 #endif
00274       pGrid->U[k][j][i].M3 = rd*rvz;
00275 #ifdef ADIABATIC
00276       pGrid->U[k][j][i].E = rp/Gamma_1
00277         + 0.5*(SQR(pGrid->U[k][j][i].M1) + SQR(pGrid->U[k][j][i].M2) 
00278              + SQR(pGrid->U[k][j][i].M3))/rd;
00279 #endif
00280 
00281 /* Initialize magnetic field.  For 3D shearing box B1=Bx, B2=By, B3=Bz
00282  *  ifield = 0 - 
00283  *  ifield = 1 - Bz=B0 sin(x1) field with zero-net-flux [default]
00284  *  ifield = 2 - uniform Bz
00285  *  ifield = 3 - B=(0,B0cos(kx*x1),B0sin(kx*x1))= zero-net flux w helicity
00286  *  ifield = 4 - B=(0,B0/sqrt(2),B0/sqrt(2))= net toroidal+vertical field
00287  */
00288 #ifdef MHD
00289       if (ifield == 0) {
00290         pGrid->B1i[k][j][i] = rbx;
00291         pGrid->B2i[k][j][i] = rby;
00292         pGrid->B3i[k][j][i] = rbz;
00293         if (i==ie) pGrid->B1i[k][j][ie+1] =  pGrid->B1i[k][j][is];
00294         if (j==je) pGrid->B2i[k][je+1][i] =  pGrid->B2i[k][js][i];
00295         if (k==ke) pGrid->B3i[ke+1][j][i] =  pGrid->B3i[ks][j][i];
00296       }
00297       if (ifield == 1) {
00298         pGrid->B1i[k][j][i] = 0.0;
00299         pGrid->B2i[k][j][i] = 0.0;
00300         pGrid->B3i[k][j][i] = B0*(sin((double)kx*x1));
00301         if (i==ie) pGrid->B1i[k][j][ie+1] = 0.0;
00302         if (j==je) pGrid->B2i[k][je+1][i] = 0.0;
00303         if (k==ke) pGrid->B3i[ke+1][j][i] = B0*(sin((double)kx*x1));
00304       }
00305       if (ifield == 2) {
00306         pGrid->B1i[k][j][i] = 0.0;
00307         pGrid->B2i[k][j][i] = 0.0;
00308         pGrid->B3i[k][j][i] = B0;
00309         if (i==ie) pGrid->B1i[k][j][ie+1] = 0.0;
00310         if (j==je) pGrid->B2i[k][je+1][i] = 0.0;
00311         if (k==ke) pGrid->B3i[ke+1][j][i] = B0;
00312       }
00313       if (ifield == 3) {
00314         pGrid->B1i[k][j][i] = 0.0;
00315         pGrid->B2i[k][j][i] = B0*(cos((double)kx*x1));
00316         pGrid->B3i[k][j][i] = B0*(sin((double)kx*x1));
00317         if (i==ie) pGrid->B1i[k][j][ie+1] = 0.0;
00318         if (j==je) pGrid->B2i[k][je+1][i] = B0*(cos((double)kx*x1));
00319         if (k==ke) pGrid->B3i[ke+1][j][i] = B0*(sin((double)kx*x1));
00320       }
00321       if (ifield == 4) {
00322         pGrid->B1i[k][j][i] = 0.0;
00323         pGrid->B2i[k][j][i] = B0/sqrt(2);
00324         pGrid->B3i[k][j][i] = B0/sqrt(2);
00325         if (i==ie) pGrid->B1i[k][j][ie+1] = 0.0;
00326         if (j==je) pGrid->B2i[k][je+1][i] = B0/sqrt(2);
00327         if (k==ke) pGrid->B3i[ke+1][j][i] = B0/sqrt(2);
00328       }
00329       if (ifield == 5) {
00330         pGrid->B1i[k][j][i] = 0.0;
00331         pGrid->B2i[k][j][i] = B0;
00332         pGrid->B3i[k][j][i] = 0.0;
00333         if (i==ie) pGrid->B1i[k][j][ie+1] = 0.0;
00334         if (j==je) pGrid->B2i[k][je+1][i] = B0;
00335         if (k==ke) pGrid->B3i[ke+1][j][i] = 0.0;
00336       }
00337 #endif /* MHD */
00338     }
00339   }}
00340 #ifdef MHD
00341   for (k=ks; k<=ke; k++) {
00342     for (j=js; j<=je; j++) {
00343       for (i=is; i<=ie; i++) {
00344         pGrid->U[k][j][i].B1c = 0.5*(pGrid->B1i[k][j][i]+pGrid->B1i[k][j][i+1]);
00345         pGrid->U[k][j][i].B2c = 0.5*(pGrid->B2i[k][j][i]+pGrid->B2i[k][j+1][i]);
00346         pGrid->U[k][j][i].B3c = 0.5*(pGrid->B3i[k][j][i]+pGrid->B3i[k+1][j][i]);
00347 #ifdef ADIABATIC
00348       pGrid->U[k][j][i].E += 0.5*(SQR(pGrid->U[k][j][i].B1c)
00349          + SQR(pGrid->U[k][j][i].B2c) + SQR(pGrid->U[k][j][i].B3c));
00350 #endif
00351       }
00352     }
00353   }
00354 #endif /* MHD */
00355 
00356 /* enroll gravitational potential function */
00357 
00358   ShearingBoxPot = UnstratifiedDisk;
00359 
00360 /* enroll new history variables, only once  */
00361 
00362   if (frst == 1) {
00363     dump_history_enroll(hst_rho_Vx_dVy, "<rho Vx dVy>");
00364     dump_history_enroll(hst_rho_dVy2, "<rho dVy^2>");
00365 #ifdef ADIABATIC
00366     dump_history_enroll(hst_E_total, "<E + rho Phi>");
00367 #endif
00368 #ifdef MHD
00369     dump_history_enroll(hst_Bx, "<Bx>");
00370     dump_history_enroll(hst_By, "<By>");
00371     dump_history_enroll(hst_Bz, "<Bz>");
00372     dump_history_enroll(hst_BxBy, "<-Bx By>");
00373     if (ipert == 5) dump_history_enroll(hst_dEw2, "<dEw2>");
00374     if (ipert == 6) dump_history_enroll(hst_dBy, "<dBy>");
00375 #endif /* MHD */
00376     frst = 0;
00377   }
00378 
00379 /* With viscosity and/or resistivity, read eta_Ohm and nu */
00380 #ifdef RESISTIVITY
00381   eta_Ohm = par_getd_def("problem","eta_O",0.0);
00382   Q_Hall  = par_getd_def("problem","Q_H",0.0);
00383   Q_AD    = par_getd_def("problem","Q_A",0.0);
00384   d_ind   = par_getd_def("problem","d_ind",0.0);
00385 #endif
00386 #ifdef VISCOSITY
00387   nu_iso = par_getd_def("problem","nu_iso",0.0);
00388   nu_aniso = par_getd_def("problem","nu_aniso",0.0);
00389 #endif
00390 
00391   return;
00392 }
00393 
00394 /*==============================================================================
00395  * PUBLIC PROBLEM USER FUNCTIONS:
00396  * problem_write_restart() - writes problem-specific user data to restart files
00397  * problem_read_restart()  - reads problem-specific user data from restart files
00398  * get_usr_expr()          - sets pointer to expression for special output data
00399  * get_usr_out_fun()       - returns a user defined output function pointer
00400  * get_usr_par_prop()      - returns a user defined particle selection function
00401  * Userwork_in_loop        - problem specific work IN     main loop
00402  * Userwork_after_loop     - problem specific work AFTER  main loop
00403  *----------------------------------------------------------------------------*/
00404 
00405 void problem_write_restart(MeshS *pM, FILE *fp)
00406 {
00407   return;
00408 }
00409 
00410 /*
00411  * 'problem_read_restart' must enroll gravity on restarts
00412  */
00413 
00414 void problem_read_restart(MeshS *pM, FILE *fp)
00415 {
00416 /* Read Omega, and with viscosity and/or resistivity, read eta_Ohm and nu */
00417 
00418   Omega_0 = par_getd_def("problem","Omega",1.0e-3);
00419   qshear  = par_getd_def("problem","qshear",1.5);
00420 
00421 #ifdef RESISTIVITY
00422   eta_Ohm = par_getd_def("problem","eta_O",0.0);
00423   Q_Hall  = par_getd_def("problem","Q_H",0.0);
00424   Q_AD    = par_getd_def("problem","Q_A",0.0);
00425   d_ind   = par_getd_def("problem","d_ind",0.0);
00426 #endif
00427 
00428 #ifdef VISCOSITY
00429   nu_iso = par_getd_def("problem","nu_iso",0.0);
00430   nu_aniso = par_getd_def("problem","nu_aniso",0.0);
00431 #endif
00432 
00433 /* enroll gravitational potential function */
00434 
00435   ShearingBoxPot = UnstratifiedDisk;
00436 
00437 /* enroll new history variables */
00438 
00439   dump_history_enroll(hst_rho_Vx_dVy, "<rho Vx dVy>");
00440   dump_history_enroll(hst_rho_dVy2, "<rho dVy^2>");
00441 #ifdef ADIABATIC
00442   dump_history_enroll(hst_E_total, "<E + rho Phi>");
00443 #endif
00444 #ifdef MHD
00445   dump_history_enroll(hst_Bx, "<Bx>");
00446   dump_history_enroll(hst_By, "<By>");
00447   dump_history_enroll(hst_Bz, "<Bz>");
00448   dump_history_enroll(hst_BxBy, "<-Bx By>");
00449 #endif /* MHD */
00450 
00451   return;
00452 }
00453 
00454 /* Get_user_expression computes dVy */
00455 ConsFun_t get_usr_expr(const char *expr)
00456 {
00457   if (strcmp(expr,"dVy")==0) return expr_dV2;
00458   if (strcmp(expr,"Jsq")==0) return expr_Jsq;
00459   return NULL;
00460 }
00461 
00462 VOutFun_t get_usr_out_fun(const char *name){
00463   return NULL;
00464 }
00465 
00466 #ifdef RESISTIVITY
00467 void get_eta_user(GridS *pG, int i, int j, int k,
00468                              Real *eta_O, Real *eta_H, Real *eta_A)
00469 {
00470   *eta_O = 0.0;
00471   *eta_H = 0.0;
00472   *eta_A = 0.0;
00473 
00474   return;
00475 }
00476 #endif
00477 
00478 void Userwork_in_loop(MeshS *pM)
00479 {
00480 }
00481 
00482 void Userwork_after_loop(MeshS *pM)
00483 {
00484 }
00485 
00486 /*------------------------------------------------------------------------------
00487  * ran2: extracted from the Numerical Recipes in C (version 2) code.  Modified
00488  *   to use doubles instead of floats. -- T. A. Gardiner -- Aug. 12, 2003
00489  */
00490 
00491 #define IM1 2147483563
00492 #define IM2 2147483399
00493 #define AM (1.0/IM1)
00494 #define IMM1 (IM1-1)
00495 #define IA1 40014
00496 #define IA2 40692
00497 #define IQ1 53668
00498 #define IQ2 52774
00499 #define IR1 12211
00500 #define IR2 3791
00501 #define NTAB 32
00502 #define NDIV (1+IMM1/NTAB)
00503 #define RNMX (1.0-DBL_EPSILON)
00504 
00505 /* Long period (> 2 x 10^{18}) random number generator of L'Ecuyer
00506  * with Bays-Durham shuffle and added safeguards.  Returns a uniform
00507  * random deviate between 0.0 and 1.0 (exclusive of the endpoint
00508  * values).  Call with idum = a negative integer to initialize;
00509  * thereafter, do not alter idum between successive deviates in a
00510  * sequence.  RNMX should appriximate the largest floating point value
00511  * that is less than 1.
00512  */
00513 
00514 double ran2(long int *idum)
00515 {
00516   int j;
00517   long int k;
00518   static long int idum2=123456789;
00519   static long int iy=0;
00520   static long int iv[NTAB];
00521   double temp;
00522 
00523   if (*idum <= 0) { /* Initialize */
00524     if (-(*idum) < 1) *idum=1; /* Be sure to prevent idum = 0 */
00525     else *idum = -(*idum);
00526     idum2=(*idum);
00527     for (j=NTAB+7;j>=0;j--) { /* Load the shuffle table (after 8 warm-ups) */
00528       k=(*idum)/IQ1;
00529       *idum=IA1*(*idum-k*IQ1)-k*IR1;
00530       if (*idum < 0) *idum += IM1;
00531       if (j < NTAB) iv[j] = *idum;
00532     }
00533     iy=iv[0];
00534   }
00535   k=(*idum)/IQ1;                 /* Start here when not initializing */
00536   *idum=IA1*(*idum-k*IQ1)-k*IR1; /* Compute idum=(IA1*idum) % IM1 without */
00537   if (*idum < 0) *idum += IM1;   /* overflows by Schrage's method */
00538   k=idum2/IQ2;
00539   idum2=IA2*(idum2-k*IQ2)-k*IR2; /* Compute idum2=(IA2*idum) % IM2 likewise */
00540   if (idum2 < 0) idum2 += IM2;
00541   j=(int)(iy/NDIV);              /* Will be in the range 0...NTAB-1 */
00542   iy=iv[j]-idum2;                /* Here idum is shuffled, idum and idum2 */
00543   iv[j] = *idum;                 /* are combined to generate output */
00544   if (iy < 1) iy += IMM1;
00545   if ((temp=AM*iy) > RNMX) return RNMX; /* No endpoint values */
00546   else return temp;
00547 }
00548 
00549 #undef IM1
00550 #undef IM2
00551 #undef AM
00552 #undef IMM1
00553 #undef IA1
00554 #undef IA2
00555 #undef IQ1
00556 #undef IQ2
00557 #undef IR1
00558 #undef IR2
00559 #undef NTAB
00560 #undef NDIV
00561 #undef RNMX
00562 
00563 /*----------------------------------------------------------------------------*/
00564 /*! \fn static Real UnstratifiedDisk(const Real x1, const Real x2,const Real x3)
00565  *  \brief tidal potential in 3D shearing box
00566  */
00567 static Real UnstratifiedDisk(const Real x1, const Real x2, const Real x3)
00568 {
00569   Real phi=0.0;
00570 #ifndef FARGO
00571   phi -= qshear*Omega_0*Omega_0*x1*x1;
00572 #endif
00573   return phi;
00574 }
00575 
00576 /*----------------------------------------------------------------------------*/
00577 /*! \fn static Real expr_dV2(const GridS *pG, const int i, const int j, 
00578  *                           const int k)
00579  *  \brief Computes delta(Vy) 
00580  */
00581 static Real expr_dV2(const GridS *pG, const int i, const int j, const int k)
00582 {
00583   Real x1,x2,x3;
00584   cc_pos(pG,i,j,k,&x1,&x2,&x3);
00585 #ifdef FARGO
00586   return (pG->U[k][j][i].M2/pG->U[k][j][i].d);
00587 #else
00588   return (pG->U[k][j][i].M2/pG->U[k][j][i].d + qshear*Omega_0*x1);
00589 #endif
00590 }
00591 
00592 /*----------------------------------------------------------------------------*/
00593 /*! \fn static Real expr_Jsq(const GridS *pG, const int i, const int j, 
00594  *                           const int k)
00595  *  \brief Computes current density square
00596  */
00597 static Real expr_Jsq(const GridS *pG, const int i, const int j, const int k)
00598 {
00599   Real J1,J2,J3;
00600 
00601   J1 = (pG->B3i[k][j][i] - pG->B3i[k  ][j-1][i  ])/pG->dx2 -
00602        (pG->B2i[k][j][i] - pG->B2i[k-1][j  ][i  ])/pG->dx3;
00603   J2 = (pG->B1i[k][j][i] - pG->B1i[k-1][j  ][i  ])/pG->dx3 -
00604        (pG->B3i[k][j][i] - pG->B3i[k  ][j  ][i-1])/pG->dx1;
00605   J3 = (pG->B2i[k][j][i] - pG->B2i[k  ][j  ][i-1])/pG->dx1 -
00606        (pG->B1i[k][j][i] - pG->B1i[k  ][j-1][i  ])/pG->dx2;
00607 
00608   return SQR(J1)+SQR(J2)+SQR(J3);
00609 }
00610 
00611 /*------------------------------------------------------------------------------
00612  * Hydro history variables:
00613  * hst_rho_Vx_dVy: Reynolds stress, added as history variable.
00614  * hst_rho_dVy2: KE in y-velocity fluctuations
00615  * hst_E_total: total energy (including tidal potential).
00616  */
00617 
00618 /*! \fn static Real hst_rho_Vx_dVy(const GridS *pG,const int i,const int j, 
00619  *                                 const int k)
00620  *  \brief Reynolds stress, added as history variable.*/
00621 static Real hst_rho_Vx_dVy(const GridS *pG,const int i,const int j, const int k)
00622 {
00623   Real x1,x2,x3;
00624   cc_pos(pG,i,j,k,&x1,&x2,&x3);
00625 #ifdef FARGO
00626   return pG->U[k][j][i].M1*(pG->U[k][j][i].M2/pG->U[k][j][i].d);
00627 #else
00628   return pG->U[k][j][i].M1*
00629     (pG->U[k][j][i].M2/pG->U[k][j][i].d + qshear*Omega_0*x1);
00630 #endif
00631 }
00632 
00633 /*! \fn static Real hst_rho_dVy2(const GridS *pG, const int i, const int j, 
00634  *                               const int k)
00635  *  \brief KE in y-velocity fluctuations */
00636 static Real hst_rho_dVy2(const GridS *pG, const int i, const int j, const int k)
00637 {
00638   Real x1,x2,x3,dVy;
00639   cc_pos(pG,i,j,k,&x1,&x2,&x3);
00640 #ifdef FARGO
00641   dVy = (pG->U[k][j][i].M2/pG->U[k][j][i].d);
00642 #else
00643   dVy = (pG->U[k][j][i].M2/pG->U[k][j][i].d + qshear*Omega_0*x1);
00644 #endif
00645   return pG->U[k][j][i].d*dVy*dVy;
00646 }
00647 
00648 #ifdef ADIABATIC
00649 /*! \fn static Real hst_E_total(const GridS *pG, const int i, const int j, 
00650  *                              const int k)
00651  *  \brief total energy (including tidal potential). */
00652 static Real hst_E_total(const GridS *pG, const int i, const int j, const int k)
00653 {
00654   Real x1,x2,x3,phi;
00655   cc_pos(pG,i,j,k,&x1,&x2,&x3);
00656   phi = UnstratifiedDisk(x1, x2, x3);
00657 
00658   return pG->U[k][j][i].E + pG->U[k][j][i].d*phi;
00659 }
00660 #endif /* ADIABATIC */
00661 
00662 /*------------------------------------------------------------------------------
00663  * MHD history variables
00664  * hst_Bx, etc.: Net flux, and Maxwell stress, added as history variables
00665  */
00666 
00667 #ifdef MHD
00668 /*! \fn static Real hst_Bx(const GridS *pG, const int i,const int j,const int k)
00669  *  \brief x-component of B-field */
00670 static Real hst_Bx(const GridS *pG, const int i, const int j, const int k)
00671 {
00672   return pG->U[k][j][i].B1c;
00673 }
00674 
00675 /*! \fn static Real hst_By(const GridS *pG, const int i,const int j,const int k)
00676  *  \brief y-component of B-field */
00677 static Real hst_By(const GridS *pG, const int i, const int j, const int k)
00678 {
00679   return pG->U[k][j][i].B2c;
00680 }
00681 
00682 /*! \fn static Real hst_Bz(const GridS *pG, const int i,const int j,const int k)
00683  *  \brief z-component of B-field */
00684 static Real hst_Bz(const GridS *pG, const int i, const int j, const int k)
00685 {
00686   return pG->U[k][j][i].B3c;
00687 }
00688 
00689 /*! \fn static Real hst_BxBy(const GridS *pG, const int i, const int j, 
00690  *                           const int k)
00691  *  \brief Maxwell stress */
00692 static Real hst_BxBy(const GridS *pG, const int i, const int j, const int k)
00693 {
00694   return -pG->U[k][j][i].B1c*pG->U[k][j][i].B2c;
00695 }
00696 
00697 static Real hst_dEw2(const GridS *pG, const int i, const int j, const int k)
00698 {
00699   Real x1,x2,x3,dVx,dVy,dVz,dBz;
00700   cc_pos(pG,i,j,k,&x1,&x2,&x3);
00701   dBz = pG->U[k][j][i].B3c-(sqrt(15.0/16.0))/(2.0*PI)/sqrt(4.*PI);
00702   dVx = pG->U[k][j][i].M1/pG->U[k][j][i].d;
00703   dVy = pG->U[k][j][i].M2/pG->U[k][j][i].d + qshear*Omega_0*x1;
00704   dVz = pG->U[k][j][i].M3/pG->U[k][j][i].d;
00705   
00706 /*  return (dVx*dVx + dVy*dVy + dVz*dVz + pG->U[k][j][i].B1c*pG->U[k][j][i].B1c
00707     + pG->U[k][j][i].B2c*pG->U[k][j][i].B2c + dBz*dBz); */
00708   return (pG->U[k][j][i].B1c*pG->U[k][j][i].B1c
00709     + pG->U[k][j][i].B2c*pG->U[k][j][i].B2c + dBz*dBz); 
00710 }
00711 
00712 static Real hst_dBy(const GridS *pG, const int i, const int j, const int k)
00713 {
00714   double fkx, fky, fkz; /* Fourier kx, ky */
00715   double dBy;
00716   Real x1,x2,x3;
00717 
00718 /* Lx,Ly, and Lz are globals */
00719 
00720   fky = 2.0*PI/Ly;
00721   fkx = -4.0*PI/Lx + qshear*Omega_0*fky*pG->time;
00722   fkz = 2.0*PI/Lz;
00723 
00724 /* compute real part of Fourier mode, for comparison to JGG fig 11 */
00725   cc_pos(pG,i,j,k,&x1,&x2,&x3);
00726   dBy = 2.0*(pG->U[k][j][i].B2c - (0.2-0.15*Omega_0*pG->time));
00727   dBy *= cos(fkx*x1 + fky*x2 + fkz*x3);
00728 
00729   return dBy;
00730 }
00731 
00732 #endif /* MHD */
00733