Athena: prob/hb3.c Source File

Go to the documentation of this file.
00001 #include "copyright.h"
00002 /*============================================================================*/
00003 /*! \file hb3.c
00004  *  \brief Problem generator for 2D MRI simulations using the shearing sheet
00005  *   based on "A powerful local shear instability in weakly magnetized disks.
00006  *
00007  * PURPOSE: Problem generator for 2D MRI simulations using the shearing sheet
00008  *   based on "A powerful local shear instability in weakly magnetized disks.
00009  *   III - Long-term evolution in a shearing sheet" by Hawley & Balbus.  This
00010  *   is the third of the HB papers on the MRI, thus hb3.
00011  *
00012  * Several different perturbations and field configurations are possible:
00013  * - ipert = 1 - isentropic perturbations to P & d [default]
00014  * - ipert = 2 - uniform Vx=amp, sinusoidal density
00015  * - ipert = 3 - random perturbations to P [used by HB]
00016  * - ipert = 4 - sinusoidal perturbation to Vx in z
00017  *
00018  * - ifield = 1 - Bz=B0 sin(x1) field with zero-net-flux [default]
00019  * - ifield = 2 - uniform Bz
00020  *
00021  * PRIVATE FUNCTION PROTOTYPES:
00022  * - ran2() - random number generator from NR
00023  * - UnstratifiedDisk() - tidal potential in 2D shearing box
00024  * - expr_dV3() - computes delta(Vy)
00025  * - hst_rho_Vx_dVy () - new history variable
00026  * - hst_E_total() - new history variable
00027  * - hst_dEk() - new history variable
00028  * - hst_Bx()  - new history variable
00029  * - hst_By()  - new history variable
00030  * - hst_Bz()  - new history variable
00031  * - hst_BxBy() - new history variable
00032  *
00033  * REFERENCE: Hawley, J. F. & Balbus, S. A., ApJ 400, 595-609 (1992).*/
00034 /*============================================================================*/
00035 
00036 #include <float.h>
00037 #include <math.h>
00038 #include <stdio.h>
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 #ifndef SHEARING_BOX
00047 #error : The HB3 problem requires shearing-box to be enabled.
00048 #endif
00049 
00050 /*==============================================================================
00051  * PRIVATE FUNCTION PROTOTYPES:
00052  * ran2() - random number generator from NR
00053  * UnstratifiedDisk() - tidal potential in 2D shearing box
00054  * expr_dV3() - computes delta(Vy)
00055  * hst_rho_Vx_dVy () - new history variable
00056  * hst_E_total() - new history variable
00057  * hst_dEk() - new history variable
00058  * hst_Bx()  - new history variable
00059  * hst_By()  - new history variable
00060  * hst_Bz()  - new history variable
00061  * hst_BxBy() - new history variable
00062  *============================================================================*/
00063 
00064 static double ran2(long int *idum);
00065 static Real UnstratifiedDisk(const Real x1, const Real x2, const Real x3);
00066 static Real expr_dV3(const GridS *pG, const int i, const int j, const int k);
00067 static Real hst_rho_Vx_dVy(const GridS *pG,const int i,const int j,const int k);
00068 static Real hst_dEk(const GridS *pG, const int i, const int j, const int k);
00069 static Real hst_E_total(const GridS *pG, const int i, const int j, const int k);
00070 #ifdef MHD
00071 static Real hst_Bx(const GridS *pG, const int i, const int j, const int k);
00072 static Real hst_By(const GridS *pG, const int i, const int j, const int k);
00073 static Real hst_Bz(const GridS *pG, const int i, const int j, const int k);
00074 static Real hst_BxBy(const GridS *pG, const int i, const int j, const int k);
00075 #endif
00076 
00077 /* boxsize, made a global variable so can be accessed by bval, etc. routines */
00078 static Real Lx,Lz;
00079 
00080 /*=========================== PUBLIC FUNCTIONS ===============================*/
00081 /*----------------------------------------------------------------------------*/
00082 /* problem:  */
00083 
00084 void problem(DomainS *pDomain)
00085 {
00086   GridS *pGrid=(pDomain->Grid);
00087   int is = pGrid->is, ie = pGrid->ie;
00088   int js = pGrid->js, je = pGrid->je;
00089   int ks = pGrid->ks;
00090   int i,j,ipert,ifield;
00091   long int iseed = -1; /* Initialize on the first call to ran2 */
00092   Real x1, x2, x3, x1min, x1max, x2min, x2max;
00093   Real den = 1.0, pres = 1.0e-5, rd, rp, rvx;
00094   Real beta,B0,kx,kz,amp;
00095   double rval;
00096   static int frst=1;  /* flag so new history variables enrolled only once */
00097 
00098 /* specify xz (r-z) plane */
00099   ShBoxCoord = xz;
00100 
00101   if (pGrid->Nx[1] == 1){
00102     ath_error("[problem]: HB3 only works on a 2D grid\n");
00103   }
00104 
00105   if (pGrid->Nx[2] > 1){
00106     ath_error("[problem]: HB3 does not work on 3D grid\n");
00107   }
00108 
00109 /* Initialize boxsize */
00110   x1min = pDomain->RootMinX[0];
00111   x1max = pDomain->RootMaxX[0];
00112   Lx = x1max - x1min;
00113   kx = 2.0*PI/Lx;
00114 
00115   x2min = pDomain->RootMinX[1];
00116   x2max = pDomain->RootMaxX[1];
00117   Lz = x2max - x2min;
00118   kz = 2.0*PI/Lz;
00119 
00120 /* Read problem parameters */
00121   Omega_0 = par_getd_def("problem","omega",1.0e-3);
00122   qshear  = par_getd_def("problem","qshear",1.5);
00123   amp = par_getd("problem","amp");
00124   beta = par_getd("problem","beta");
00125   B0 = sqrt((double)(2.0*pres/beta));
00126   ifield = par_geti_def("problem","ifield", 1);
00127   ipert = par_geti_def("problem","ipert", 1);
00128 
00129   for (j=js; j<=je; j++) {
00130     for (i=is; i<=ie; i++) {
00131       cc_pos(pGrid,i,j,ks,&x1,&x2,&x3);
00132 
00133 /* Initialize perturbations
00134  *  ipert = 1 - isentropic perturbations to P & d [default]
00135  *  ipert = 2 - uniform Vx=amp, sinusoidal density
00136  *  ipert = 3 - random perturbations to P [used by HB]
00137  *  ipert = 4 - sinusoidal perturbation to Vx in z
00138  */
00139       if (ipert == 1) {
00140         rval = 1.0 + amp*(ran2(&iseed) - 0.5);
00141 #ifdef ADIABATIC
00142         rp = rval*pres;
00143         rd = pow(rval,1.0/Gamma)*den;
00144 #else
00145         rd = rval*den;
00146 #endif
00147         rvx = 0.0;
00148       }
00149       if (ipert == 2) {
00150         rp = pres;
00151         rd = den*(1.0 + 0.1*sin((double)kx*x1));
00152 #ifdef ADIABATIC
00153         rvx = amp*sqrt(Gamma*pres/den);
00154 #else
00155         rvx = amp*sqrt(pres/den);
00156 #endif
00157       }
00158       if (ipert == 3) {
00159         rval = 1.0 + amp*(ran2(&iseed) - 0.5);
00160 #ifdef ADIABATIC
00161         rp = rval*pres;
00162         rd = den;
00163 #else
00164         rd = rval*den;
00165 #endif
00166         rvx = 0.0;
00167       }
00168       if (ipert == 4) {
00169         rp = pres;
00170         rd = den;
00171         rvx = amp*sin((double)kz*x2);
00172       }
00173 
00174 /* Initialize d, M, and P.  For 2D shearing box M1=Vx, M2=Vz, M3=Vy */ 
00175 
00176       pGrid->U[ks][j][i].d  = rd;
00177       pGrid->U[ks][j][i].M1 = rd*rvx;
00178       pGrid->U[ks][j][i].M2 = 0.0;
00179 #ifdef FARGO
00180       pGrid->U[ks][j][i].M3 = 0.0;
00181 #else
00182       pGrid->U[ks][j][i].M3 = -rd*qshear*Omega_0*x1;
00183 #endif
00184 #ifdef ADIABATIC
00185       pGrid->U[ks][j][i].E = rp/Gamma_1
00186         + 0.5*(SQR(pGrid->U[ks][j][i].M1) + SQR(pGrid->U[ks][j][i].M3))/rd;
00187 #endif
00188 
00189 /* Initialize magnetic field.  For 2D shearing box B1=Bx, B2=Bz, B3=By
00190  *  ifield = 1 - Bz=B0 sin(x1) field with zero-net-flux [default]
00191  *  ifield = 2 - uniform Bz
00192  */
00193 #ifdef MHD
00194       if (ifield == 1) {
00195         pGrid->U[ks][j][i].B1c = 0.0;
00196         pGrid->U[ks][j][i].B2c = B0*(sin((double)kx*x1));
00197         pGrid->U[ks][j][i].B3c = 0.0;
00198         pGrid->B1i[ks][j][i] = 0.0;
00199         pGrid->B2i[ks][j][i] = B0*(sin((double)kx*x1));
00200         pGrid->B3i[ks][j][i] = 0.0;
00201         if (i==ie) pGrid->B1i[ks][j][ie+1] = 0.0;
00202         if (j==je) pGrid->B2i[ks][je+1][i] = B0*(sin((double)kx*x1));
00203       }
00204       if (ifield == 2) {
00205         pGrid->U[ks][j][i].B1c = 0.0;
00206         pGrid->U[ks][j][i].B2c = B0;
00207         pGrid->U[ks][j][i].B3c = 0.0;
00208         pGrid->B1i[ks][j][i] = 0.0;
00209         pGrid->B2i[ks][j][i] = B0;
00210         pGrid->B3i[ks][j][i] = 0.0;
00211         if (i==ie) pGrid->B1i[ks][j][ie+1] = 0.0;
00212         if (j==je) pGrid->B2i[ks][je+1][i] = B0;
00213       }
00214 #ifdef ADIABATIC
00215       pGrid->U[ks][j][i].E += 0.5*(SQR(pGrid->U[ks][j][i].B1c)
00216          + SQR(pGrid->U[ks][j][i].B2c) + SQR(pGrid->U[ks][j][i].B3c));
00217 #endif
00218 #endif /* MHD */
00219     }
00220   }
00221 
00222 /* With viscosity and/or resistivity, read eta_Ohm and nu_V */
00223 #ifdef RESISTIVITY
00224   eta_Ohm = par_getd_def("problem","eta_O",0.0);
00225   Q_Hall  = par_getd_def("problem","Q_H",0.0);
00226   Q_AD    = par_getd_def("problem","Q_A",0.0);
00227 #endif
00228 
00229 /* enroll gravitational potential function, shearing sheet BC functions */
00230 
00231   ShearingBoxPot = UnstratifiedDisk;
00232 
00233 /* enroll new history variables */
00234 
00235   if (frst == 1) {
00236     dump_history_enroll(hst_dEk, "<0.5rho(Vx^2+4dVy^2)>");
00237     dump_history_enroll(hst_rho_Vx_dVy, "<rho Vx dVy>");
00238     dump_history_enroll(hst_E_total, "<E + rho Phi>");
00239 #ifdef MHD
00240     dump_history_enroll(hst_Bx, "<Bx>");
00241     dump_history_enroll(hst_By, "<By>");
00242     dump_history_enroll(hst_Bz, "<Bz>");
00243     dump_history_enroll(hst_BxBy, "<-Bx By>");
00244 #endif /* MHD */
00245     frst = 0;
00246   }
00247 
00248   return;
00249 }
00250 
00251 /*==============================================================================
00252  * PROBLEM USER FUNCTIONS:
00253  * problem_write_restart() - writes problem-specific user data to restart files
00254  * problem_read_restart()  - reads problem-specific user data from restart files
00255  * get_usr_expr()          - sets pointer to expression for special output data
00256  * get_usr_out_fun()       - returns a user defined output function pointer
00257  * get_usr_par_prop()      - returns a user defined particle selection function
00258  * Userwork_in_loop        - problem specific work IN     main loop
00259  * Userwork_after_loop     - problem specific work AFTER  main loop
00260  *----------------------------------------------------------------------------*/
00261 
00262 void problem_write_restart(MeshS *pM, FILE *fp)
00263 {
00264   return;
00265 }
00266 
00267 /*
00268  * 'problem_read_restart' must enroll special boundary value functions,
00269  *    and initialize gravity on restarts
00270  */
00271 
00272 void problem_read_restart(MeshS *pM, FILE *fp)
00273 {
00274   Real x1min, x1max;
00275 
00276   Omega_0 = par_getd_def("problem","omega",1.0e-3);
00277   qshear  = par_getd_def("problem","qshear",1.5);
00278 
00279 #ifdef RESISTIVITY
00280   eta_Ohm = par_getd_def("problem","eta_O",0.0);
00281   Q_Hall  = par_getd_def("problem","Q_H",0.0);
00282   Q_AD    = par_getd_def("problem","Q_A",0.0);
00283 #endif
00284 
00285 /* Must recompute global variable Lx needed by BC routines */
00286   x1min = pM->RootMinX[0];
00287   x1max = pM->RootMaxX[0];
00288   Lx = x1max - x1min;
00289 
00290   ShearingBoxPot = UnstratifiedDisk;
00291 
00292   return;
00293 }
00294 
00295 /* Get_user_expression computes dVy */
00296 ConsFun_t get_usr_expr(const char *expr)
00297 {
00298   if(strcmp(expr,"dVy")==0) return expr_dV3;
00299   return NULL;
00300 }
00301 
00302 VOutFun_t get_usr_out_fun(const char *name){
00303   return NULL;
00304 }
00305 
00306 void Userwork_in_loop(MeshS *pM)
00307 {
00308 }
00309 
00310 void Userwork_after_loop(MeshS *pM)
00311 {
00312 }
00313 
00314 /*=========================== PRIVATE FUNCTIONS ==============================*/
00315 
00316 /*----------------------------------------------------------------------------*/
00317 
00318 #define IM1 2147483563
00319 #define IM2 2147483399
00320 #define AM (1.0/IM1)
00321 #define IMM1 (IM1-1)
00322 #define IA1 40014
00323 #define IA2 40692
00324 #define IQ1 53668
00325 #define IQ2 52774
00326 #define IR1 12211
00327 #define IR2 3791
00328 #define NTAB 32
00329 #define NDIV (1+IMM1/NTAB)
00330 #define RNMX (1.0-DBL_EPSILON)
00331 
00332 /*! \fn double ran2(long int *idum)
00333  *  \brief Extracted from the Numerical Recipes in C (version 2) code.  Modified
00334  *   to use doubles instead of floats. -- T. A. Gardiner -- Aug. 12, 2003
00335  * 
00336  * Long period (> 2 x 10^{18}) random number generator of L'Ecuyer
00337  * with Bays-Durham shuffle and added safeguards.  Returns a uniform
00338  * random deviate between 0.0 and 1.0 (exclusive of the endpoint
00339  * values).  Call with idum = a negative integer to initialize;
00340  * thereafter, do not alter idum between successive deviates in a
00341  * sequence.  RNMX should appriximate the largest floating point value
00342  * that is less than 1.
00343  */
00344 
00345 double ran2(long int *idum)
00346 {
00347   int j;
00348   long int k;
00349   static long int idum2=123456789;
00350   static long int iy=0;
00351   static long int iv[NTAB];
00352   double temp;
00353 
00354   if (*idum <= 0) { /* Initialize */
00355     if (-(*idum) < 1) *idum=1; /* Be sure to prevent idum = 0 */
00356     else *idum = -(*idum);
00357     idum2=(*idum);
00358     for (j=NTAB+7;j>=0;j--) { /* Load the shuffle table (after 8 warm-ups) */
00359       k=(*idum)/IQ1;
00360       *idum=IA1*(*idum-k*IQ1)-k*IR1;
00361       if (*idum < 0) *idum += IM1;
00362       if (j < NTAB) iv[j] = *idum;
00363     }
00364     iy=iv[0];
00365   }
00366   k=(*idum)/IQ1;                 /* Start here when not initializing */
00367   *idum=IA1*(*idum-k*IQ1)-k*IR1; /* Compute idum=(IA1*idum) % IM1 without */
00368   if (*idum < 0) *idum += IM1;   /* overflows by Schrage's method */
00369   k=idum2/IQ2;
00370   idum2=IA2*(idum2-k*IQ2)-k*IR2; /* Compute idum2=(IA2*idum) % IM2 likewise */
00371   if (idum2 < 0) idum2 += IM2;
00372   j=(int)(iy/NDIV);              /* Will be in the range 0...NTAB-1 */
00373   iy=iv[j]-idum2;                /* Here idum is shuffled, idum and idum2 */
00374   iv[j] = *idum;                 /* are combined to generate output */
00375   if (iy < 1) iy += IMM1;
00376   if ((temp=AM*iy) > RNMX) return RNMX; /* No endpoint values */
00377   else return temp;
00378 }
00379 
00380 #undef IM1
00381 #undef IM2
00382 #undef AM
00383 #undef IMM1
00384 #undef IA1
00385 #undef IA2
00386 #undef IQ1
00387 #undef IQ2
00388 #undef IR1
00389 #undef IR2
00390 #undef NTAB
00391 #undef NDIV
00392 #undef RNMX
00393 
00394 /*----------------------------------------------------------------------------*/
00395 /*! \fn static Real UnstratifiedDisk(const Real x1, const Real x2,const Real x3)
00396  *  \brief  ShearingBoxPot 
00397  */
00398 static Real UnstratifiedDisk(const Real x1, const Real x2, const Real x3){
00399   Real phi=0.0;
00400 #ifndef FARGO
00401   phi -= qshear*SQR(Omega_0*x1);
00402 #endif
00403   return phi;
00404 }
00405 
00406 /*----------------------------------------------------------------------------*/
00407 /*! \fn static Real expr_dV3(const GridS *pG, const int i, const int j, 
00408  *                           const int k)
00409  *  \brief Computes delta(Vy) 
00410  */
00411 static Real expr_dV3(const GridS *pG, const int i, const int j, const int k)
00412 {
00413   Real x1,x2,x3;
00414   cc_pos(pG,i,j,k,&x1,&x2,&x3);
00415 #ifdef FARGO
00416   return pG->U[k][j][i].M3/pG->U[k][j][i].d;
00417 #else
00418   return (pG->U[k][j][i].M3/pG->U[k][j][i].d + qshear*Omega_0*x1);
00419 #endif
00420 }
00421 
00422 /*----------------------------------------------------------------------------*/
00423 /*! \fn static Real hst_rho_Vx_dVy(const GridS *pG, const int i, const int j, 
00424  *                                 const int k)
00425  *  \brief Reynolds stress, added as history variable.
00426  */
00427 static Real hst_rho_Vx_dVy(const GridS *pG, const int i, const int j, const int k)
00428 {
00429   Real x1,x2,x3;
00430   cc_pos(pG,i,j,k,&x1,&x2,&x3);
00431 #ifdef FARGO
00432   return pG->U[k][j][i].M1*pG->U[k][j][i].M3/pG->U[k][j][i].d;
00433 #else
00434   return pG->U[k][j][i].M1*(pG->U[k][j][i].M3/pG->U[k][j][i].d + qshear*Omega_0*x1);
00435 #endif
00436 }
00437 
00438 /*----------------------------------------------------------------------------*/
00439 /*! \fn static Real hst_dEk(const GridS *pG, const int i, const int j, 
00440  *                          const int k)
00441  *  \brief computes 0.5*(Vx^2 + 4(\delta Vy)^2), which for epicyclic motion
00442  *   is a constant, added as history variable */
00443 static Real hst_dEk(const GridS *pG, const int i, const int j, const int k)
00444 {
00445   Real x1,x2,x3;
00446   Real dMy, dE;
00447 
00448 #ifdef FARGO
00449   dMy = pG->U[k][j][i].M3;
00450 #else
00451   dMy = (pG->U[k][j][i].M3 + qshear*Omega_0*x1*pG->U[k][j][i].d);
00452 #endif
00453   dE = 0.5*(pG->U[k][j][i].M1*pG->U[k][j][i].M1 + 4.0*dMy*dMy)/pG->U[k][j][i].d;
00454 
00455   return dE;
00456 }
00457 
00458 /*------------------------------------------------------------------------------
00459  * hst_E_total: total energy (including tidal potential).
00460  */
00461 
00462 /*! \fn static Real hst_E_total(const GridS *pG, const int i, const int j, 
00463  *                              const int k)
00464  *  \brief Total energy, including tidal potential */
00465 static Real hst_E_total(const GridS *pG, const int i, const int j, const int k)
00466 {
00467 #ifdef ADIABATIC
00468   Real x1,x2,x3,phi;
00469   cc_pos(pG,i,j,k,&x1,&x2,&x3);
00470   phi = UnstratifiedDisk(x1, x2, x3);
00471 
00472   return pG->U[k][j][i].E + pG->U[k][j][i].d*phi;
00473 #else
00474   return 0.0;
00475 #endif
00476 }
00477 
00478 /*------------------------------------------------------------------------------
00479  * hst_Bx, etc.: Net flux, and Maxwell stress, added as history variables
00480  */
00481 
00482 #ifdef MHD
00483 /*! \fn static Real hst_Bx(const GridS *pG, const int i,const int j,const int k)
00484  *  \brief x-component of magnetic field */
00485 static Real hst_Bx(const GridS *pG, const int i, const int j, const int k)
00486 {
00487   return pG->U[k][j][pG->is].B1c;
00488 }
00489 
00490 /*! \fn static Real hst_By(const GridS *pG, const int i,const int j,const int k)
00491  *  \brief y-component of the magnetic field */
00492 static Real hst_By(const GridS *pG, const int i, const int j, const int k)
00493 {
00494   return pG->U[k][pG->js][i].B2c;
00495 }
00496 
00497 /*! \fn static Real hst_Bz(const GridS *pG, const int i,const int j,const int k)
00498  *  \brief z-component of the magnetic field */
00499 static Real hst_Bz(const GridS *pG, const int i, const int j, const int k)
00500 {
00501   return pG->U[pG->ks][j][i].B3c;
00502 }
00503 
00504 /*! \fn static Real hst_BxBy(const GridS *pG, const int i, const int j, 
00505  *                           const int k)
00506  *  \brief Maxwell stress */
00507 static Real hst_BxBy(const GridS *pG, const int i, const int j, const int k)
00508 {
00509   return -pG->U[k][j][i].B1c*pG->U[k][j][i].B2c;
00510 }
00511 
00512 #endif /* MHD */