Athena: prob/firehose.c Source File

Go to the documentation of this file.
00001 #include "copyright.h"
00002 /*============================================================================*/
00003 /*! \file firehose.c
00004  *  \brief Problem generator for firehose test of Braginskii viscosity.
00005  *
00006  * PURPOSE: Problem generator for firehose test of Braginskii viscosity.
00007  * - iprob=1: vortex test suggested by Steve Cowley, developed by Greg Hammett
00008  * - iprob=2: shear test from Appendix C.3.4 in Prateek Sharma's thesis.
00009  *
00010  * PRIVATE FUNCTION PROTOTYPES:
00011  * - ran2()    - random number generator from NR
00012  * - pbc_ix1() - sets BCs on L-x1 (left edge) of grid used in shear test
00013  * - pbc_ox1() - sets BCs on R-x1 (right edge) of grid used in shear test */
00014 /*============================================================================*/
00015 
00016 #include <float.h>
00017 #include <math.h>
00018 #include <stdio.h>
00019 #include <stdlib.h>
00020 #include "defs.h"
00021 #include "athena.h"
00022 #include "globals.h"
00023 #include "prototypes.h"
00024 
00025 /*==============================================================================
00026  * PRIVATE FUNCTION PROTOTYPES:
00027  * ran2()    - random number generator from NR
00028  * pbc_ix1() - sets BCs on L-x1 (left edge) of grid used in shear test
00029  * pbc_ox1() - sets BCs on R-x1 (right edge) of grid used in shear test
00030  *============================================================================*/
00031 
00032 static double ran2(long int *idum);
00033 static void pbc_ix1(Grid *pGrid);
00034 static void pbc_ox1(Grid *pGrid);
00035 
00036 /*----------------------------------------------------------------------------*/
00037 /* problem:   */
00038 
00039 void problem(Grid *pGrid, Domain *pDomain)
00040 {
00041   int i,is,ie,j,js,je,ks,nx1,nx2,iprob;
00042   long int iseed = -1; /* Initialize on the first call to ran2 */
00043   Real d0,p0,B0,v0,x1,x2,x3,x1min,x1max,Lx,k0,press,amp,kp;
00044 
00045   is = pGrid->is; ie = pGrid->ie;
00046   js = pGrid->js; je = pGrid->je;
00047   ks = pGrid->ks;
00048   nx1 = (ie-is)+1;
00049   nx2 = (je-js)+1;
00050   if ((nx1 == 1) || (nx2 == 1)) {
00051     ath_error("[firehose]: This problem can only be run in 2D\n");
00052   }
00053   d0 = 1.0;
00054   p0 = 1.0;
00055   B0 = par_getd("problem","B0");
00056   v0 = par_getd("problem","v0");
00057 #ifdef BRAGINSKII
00058   nu_V = par_getd("problem","nu");
00059 #endif
00060   amp = par_getd("problem","amp");
00061   kp =  par_getd("problem","kp");
00062   iprob =  par_getd("problem","iprob");
00063 
00064 /* Initialize wavenumber */
00065   x1min = par_getd("grid","x1min");
00066   x1max = par_getd("grid","x1max");
00067   Lx = x1max - x1min;
00068   k0 = 2.0*PI/Lx;
00069 
00070 /* iprob=1: Cowley/Hammett vortex test ---------------------------------------*/
00071 /* Initialize density, momentum, face-centered fields */
00072 
00073   if (iprob == 1) {
00074   for (j=js; j<=je; j++) {
00075     for (i=is; i<=ie; i++) {
00076 /* Calculate the cell center positions */
00077       cc_pos(pGrid,i,j,ks,&x1,&x2,&x3);
00078 
00079       pGrid->U[ks][j][i].d = d0;
00080       pGrid->U[ks][j][i].M1 = -d0*v0*cos(k0*x2)*sin(k0*x1);
00081       pGrid->U[ks][j][i].M2 = d0*v0*sin(k0*x2)*cos(k0*x1)+amp*cos(kp*x1);
00082       pGrid->U[ks][j][i].M3 = 0.0;
00083 #ifdef MHD
00084       pGrid->B1i[ks][j][i] = B0;
00085       pGrid->B2i[ks][j][i] = 0.0;
00086       pGrid->B3i[ks][j][i] = 0.0;
00087       pGrid->U[ks][j][i].B1c = B0;
00088       pGrid->U[ks][j][i].B2c = 0.0;
00089       pGrid->U[ks][j][i].B3c = 0.0;
00090 #endif /* MHD */
00091       press = p0 + 0.5*d0*v0*v0*(cos(k0*x1)*cos(k0*x1) + cos(k0*x2)*cos(k0*x2));
00092       pGrid->U[ks][j][i].E = press/Gamma_1
00093 #ifdef MHD
00094         + 0.5*(SQR(pGrid->U[ks][j][i].B1c) + SQR(pGrid->U[ks][j][i].B2c)
00095              + SQR(pGrid->U[ks][j][i].B3c))
00096 #endif /* MHD */
00097         + 0.5*(SQR(pGrid->U[ks][j][i].M1) + SQR(pGrid->U[ks][j][i].M2)
00098               + SQR(pGrid->U[ks][j][i].M3))/pGrid->U[ks][j][i].d;
00099     }
00100   }
00101 #ifdef MHD
00102 /* boundary conditions on interface B */
00103   for (j=js; j<=je; j++) {
00104     pGrid->B1i[ks][j][ie+1] = B0;
00105   }
00106   for (i=is; i<=ie; i++) {
00107     pGrid->B2i[ks][je+1][i] = 0.0;
00108   }
00109 #endif /* MHD */
00110   }
00111 
00112 /* iprob=2: Sharma shear test ----------------------------------------------- */
00113 /* Initialize density, momentum, face-centered fields */
00114 
00115   if (iprob == 2) {
00116   for (j=js; j<=je; j++) {
00117     for (i=is; i<=ie; i++) {
00118 /* Calculate the cell center positions */
00119       cc_pos(pGrid,i,j,ks,&x1,&x2,&x3);
00120 
00121       pGrid->U[ks][j][i].d = d0;
00122       pGrid->U[ks][j][i].M1 = d0*amp*(ran2(&iseed) - 0.5);
00123       pGrid->U[ks][j][i].M2 = d0*(amp*(ran2(&iseed) - 0.5) - 0.015*x1);
00124       pGrid->U[ks][j][i].M3 = 0.0;
00125 #ifdef MHD
00126       pGrid->B1i[ks][j][i] = B0;
00127       pGrid->B2i[ks][j][i] = B0;
00128       pGrid->B3i[ks][j][i] = 0.0;
00129       pGrid->U[ks][j][i].B1c = B0;
00130       pGrid->U[ks][j][i].B2c = B0;
00131       pGrid->U[ks][j][i].B3c = 0.0;
00132 #endif /* MHD */
00133       pGrid->U[ks][j][i].E = 0.1/Gamma_1
00134 #ifdef MHD
00135         + 0.5*(SQR(pGrid->U[ks][j][i].B1c) + SQR(pGrid->U[ks][j][i].B2c)
00136              + SQR(pGrid->U[ks][j][i].B3c))
00137 #endif /* MHD */
00138         + 0.5*(SQR(pGrid->U[ks][j][i].M1) + SQR(pGrid->U[ks][j][i].M2)
00139               + SQR(pGrid->U[ks][j][i].M3))/pGrid->U[ks][j][i].d;
00140     }
00141   }
00142 #ifdef MHD
00143 /* boundary conditions on interface B */
00144   for (j=js; j<=je; j++) {
00145     pGrid->B1i[ks][j][ie+1] = B0;
00146   }
00147   for (i=is; i<=ie; i++) {
00148     pGrid->B2i[ks][je+1][i] = B0;
00149   }
00150 #endif /* MHD */
00151 /* Enroll special BCs for shearing sheet */
00152   set_bvals_mhd_fun(left_x1,  pbc_ix1);
00153   set_bvals_mhd_fun(right_x1, pbc_ox1);
00154   }
00155 
00156   return;
00157 }
00158 
00159 /*==============================================================================
00160  * PROBLEM USER FUNCTIONS:
00161  * problem_write_restart() - writes problem-specific user data to restart files
00162  * problem_read_restart()  - reads problem-specific user data from restart files
00163  * get_usr_expr()          - sets pointer to expression for special output data
00164  * get_usr_out_fun()       - returns a user defined output function pointer
00165  * get_usr_par_prop()      - returns a user defined particle selection function
00166  * Userwork_in_loop        - problem specific work IN     main loop
00167  * Userwork_after_loop     - problem specific work AFTER  main loop
00168  *----------------------------------------------------------------------------*/
00169 
00170 void problem_write_restart(Grid *pG, Domain *pD, FILE *fp)
00171 {
00172   return;
00173 }
00174 
00175 void problem_read_restart(Grid *pG, Domain *pD, FILE *fp)
00176 {
00177   return;
00178 }
00179 
00180 Gasfun_t get_usr_expr(const char *expr)
00181 {
00182   return NULL;
00183 }
00184 
00185 VGFunout_t get_usr_out_fun(const char *name){
00186   return NULL;
00187 }
00188 
00189 #ifdef PARTICLES
00190 PropFun_t get_usr_par_prop(const char *name)
00191 {
00192   return NULL;
00193 }
00194 
00195 void gasvshift(const Real x1, const Real x2, const Real x3,
00196                                     Real *u1, Real *u2, Real *u3)
00197 {
00198   return;
00199 }
00200 
00201 void Userforce_particle(Vector *ft, const Real x1, const Real x2, const Real x3,
00202                                           Real *w1, Real *w2, Real *w3)
00203 {
00204   return;
00205 }
00206 #endif
00207 
00208 void Userwork_in_loop(Grid *pGrid, Domain *pDomain)
00209 {
00210 }
00211 
00212 void Userwork_after_loop(Grid *pGrid, Domain *pDomain)
00213 {
00214 }
00215 
00216 /*----------------------------------------------------------------------------*/
00217 
00218 #define IM1 2147483563
00219 #define IM2 2147483399
00220 #define AM (1.0/IM1)
00221 #define IMM1 (IM1-1)
00222 #define IA1 40014
00223 #define IA2 40692
00224 #define IQ1 53668
00225 #define IQ2 52774
00226 #define IR1 12211
00227 #define IR2 3791
00228 #define NTAB 32
00229 #define NDIV (1+IMM1/NTAB)
00230 #define RNMX (1.0-DBL_EPSILON)
00231 
00232 /*! \fn double ran2(long int *idum)
00233  *  \brief Extracted from the Numerical Recipes in C (version 2) code.  Modified
00234  *   to use doubles instead of floats. -- T. A. Gardiner -- Aug. 12, 2003
00235  * 
00236  * Long period (> 2 x 10^{18}) random number generator of L'Ecuyer
00237  * with Bays-Durham shuffle and added safeguards.  Returns a uniform
00238  * random deviate between 0.0 and 1.0 (exclusive of the endpoint
00239  * values).  Call with idum = a negative integer to initialize;
00240  * thereafter, do not alter idum between successive deviates in a
00241  * sequence.  RNMX should appriximate the largest floating point value
00242  * that is less than 1.
00243  */
00244 
00245 double ran2(long int *idum)
00246 {
00247   int j;
00248   long int k;
00249   static long int idum2=123456789;
00250   static long int iy=0;
00251   static long int iv[NTAB];
00252   double temp;
00253 
00254   if (*idum <= 0) { /* Initialize */
00255     if (-(*idum) < 1) *idum=1; /* Be sure to prevent idum = 0 */
00256     else *idum = -(*idum);
00257     idum2=(*idum);
00258     for (j=NTAB+7;j>=0;j--) { /* Load the shuffle table (after 8 warm-ups) */
00259       k=(*idum)/IQ1;
00260       *idum=IA1*(*idum-k*IQ1)-k*IR1;
00261       if (*idum < 0) *idum += IM1;
00262       if (j < NTAB) iv[j] = *idum;
00263     }
00264     iy=iv[0];
00265   }
00266   k=(*idum)/IQ1;                 /* Start here when not initializing */
00267   *idum=IA1*(*idum-k*IQ1)-k*IR1; /* Compute idum=(IA1*idum) % IM1 without */
00268   if (*idum < 0) *idum += IM1;   /* overflows by Schrage's method */
00269   k=idum2/IQ2;
00270   idum2=IA2*(idum2-k*IQ2)-k*IR2; /* Compute idum2=(IA2*idum) % IM2 likewise */
00271   if (idum2 < 0) idum2 += IM2;
00272   j=(int)(iy/NDIV);              /* Will be in the range 0...NTAB-1 */
00273   iy=iv[j]-idum2;                /* Here idum is shuffled, idum and idum2 */
00274   iv[j] = *idum;                 /* are combined to generate output */
00275   if (iy < 1) iy += IMM1;
00276   if ((temp=AM*iy) > RNMX) return RNMX; /* No endpoint values */
00277   else return temp;
00278 }
00279 
00280 #undef IM1
00281 #undef IM2
00282 #undef AM
00283 #undef IMM1
00284 #undef IA1
00285 #undef IA2
00286 #undef IQ1
00287 #undef IQ2
00288 #undef IR1
00289 #undef IR2
00290 #undef NTAB
00291 #undef NDIV
00292 #undef RNMX
00293 
00294 /*----------------------------------------------------------------------------*/
00295 /*! \fn static void pbc_ix1(Grid *pGrid)
00296  *  \brief Special PERIODIC boundary conditions, Inner x1 boundary
00297  */
00298 
00299 static void pbc_ix1(Grid *pGrid)
00300 {
00301   int is = pGrid->is, ie = pGrid->ie;
00302   int js = pGrid->js, je = pGrid->je;
00303   int ks = pGrid->ks, ke = pGrid->ke;
00304   int i,j,k;
00305 #ifdef MHD
00306   int ju, ku; /* j-upper, k-upper */
00307 #endif
00308   Real xmin,xmax,Lx;
00309   xmin = par_getd("grid","x1min");
00310   xmax = par_getd("grid","x1max");
00311   Lx = xmax - xmin;
00312 
00313   for (k=ks; k<=ke; k++) {
00314     for (j=js; j<=je; j++) {
00315       for (i=1; i<=nghost; i++) {
00316         pGrid->U[k][j][is-i] = pGrid->U[k][j][ie-(i-1)];
00317         pGrid->U[k][j][is-i].M2 =
00318           pGrid->U[k][j][is-i].M2 + pGrid->U[k][j][is-i].d*0.015*Lx;
00319 #ifdef SPECIAL_RELATIVITY
00320         pGrid->W[k][j][is-i] = pGrid->W[k][j][ie-(i-1)];
00321 #endif
00322       }
00323     }
00324   }
00325 
00326 #ifdef MHD
00327   for (k=ks; k<=ke; k++) {
00328     for (j=js; j<=je; j++) {
00329       for (i=1; i<=nghost; i++) {
00330         pGrid->B1i[k][j][is-i] = pGrid->B1i[k][j][ie-(i-1)];
00331       }
00332     }
00333   }
00334 
00335   if (pGrid->Nx2 > 1) ju=je+1; else ju=je;
00336   for (k=ks; k<=ke; k++) {
00337     for (j=js; j<=ju; j++) {
00338       for (i=1; i<=nghost; i++) {
00339         pGrid->B2i[k][j][is-i] = pGrid->B2i[k][j][ie-(i-1)];
00340       }
00341     }
00342   }
00343 
00344   if (pGrid->Nx3 > 1) ku=ke+1; else ku=ke;
00345   for (k=ks; k<=ku; k++) {
00346     for (j=js; j<=je; j++) {
00347       for (i=1; i<=nghost; i++) {
00348         pGrid->B3i[k][j][is-i] = pGrid->B3i[k][j][ie-(i-1)];
00349       }
00350     }
00351   }
00352 #endif /* MHD */
00353 
00354   return;
00355 }
00356 
00357 /*----------------------------------------------------------------------------*/
00358 /*! \fn static void pbc_ox1(Grid *pGrid)
00359  *  \brief Special PERIODIC boundary conditions, Outer x1 boundary
00360  */
00361 
00362 static void pbc_ox1(Grid *pGrid)
00363 {
00364   int is = pGrid->is, ie = pGrid->ie;
00365   int js = pGrid->js, je = pGrid->je;
00366   int ks = pGrid->ks, ke = pGrid->ke;
00367   int i,j,k;
00368 #ifdef MHD
00369   int ju, ku; /* j-upper, k-upper */
00370 #endif
00371   Real xmin,xmax,Lx;
00372   xmin = par_getd("grid","x1min");
00373   xmax = par_getd("grid","x1max");
00374   Lx = xmax - xmin;
00375 
00376   for (k=ks; k<=ke; k++) {
00377     for (j=js; j<=je; j++) {
00378       for (i=1; i<=nghost; i++) {
00379         pGrid->U[k][j][ie+i] = pGrid->U[k][j][is+(i-1)];
00380         pGrid->U[k][j][ie+i].M2 =
00381           pGrid->U[k][j][ie+i].M2 - pGrid->U[k][j][ie+i].d*0.015*Lx;
00382 #ifdef SPECIAL_RELATIVITY
00383         pGrid->W[k][j][ie+i] = pGrid->W[k][j][is+(i-1)];
00384 #endif
00385       }
00386     }
00387   }
00388 
00389 #ifdef MHD
00390 /* Note that i=ie+1 is not a boundary condition for the interface field B1i */
00391   for (k=ks; k<=ke; k++) {
00392     for (j=js; j<=je; j++) {
00393       for (i=2; i<=nghost; i++) {
00394         pGrid->B1i[k][j][ie+i] = pGrid->B1i[k][j][is+(i-1)];
00395       }
00396     }
00397   }
00398 
00399   if (pGrid->Nx2 > 1) ju=je+1; else ju=je;
00400   for (k=ks; k<=ke; k++) {
00401     for (j=js; j<=ju; j++) {
00402       for (i=1; i<=nghost; i++) {
00403         pGrid->B2i[k][j][ie+i] = pGrid->B2i[k][j][is+(i-1)];
00404       }
00405     }
00406   }
00407 
00408   if (pGrid->Nx3 > 1) ku=ke+1; else ku=ke;
00409   for (k=ks; k<=ku; k++) {
00410     for (j=js; j<=je; j++) {
00411       for (i=1; i<=nghost; i++) {
00412         pGrid->B3i[k][j][ie+i] = pGrid->B3i[k][j][is+(i-1)];
00413       }
00414     }
00415   }
00416 #endif /* MHD */
00417 
00418   return;
00419 }