Athena: prob/rt.c Source File

Go to the documentation of this file.
00001 #include "copyright.h"
00002 /*============================================================================*/
00003 /*! \file rt.c
00004  *  \brief Problem generator for RT instabilty.
00005  *
00006  * PURPOSE: Problem generator for RT instabilty.  Gravitational pot. is
00007  *   hardwired to be 0.1z. Density difference is hardwired to be 2.0 in 2D, and
00008  *   is set by the input parameter <problem>/rhoh in 3D (default value is 3.0).
00009  *   This reproduces 2D results of Liska & Wendroff, 3D results of
00010  *   Dimonte et al.
00011  * 
00012  * FOR 2D HYDRO:
00013  * Problem domain should be -1/6 < x < 1/6; -0.5 < y < 0.5 with gamma=1.4 to
00014  * match Liska & Wendroff. Interface is at y=0; perturbation added to Vy
00015  * Gravity acts in the y-direction.  Special reflecting boundary conditions
00016  *   added in x2 to improve hydrostatic eqm (prevents launching of weak waves)
00017  * Atwood number A = (d2-d1)/(d2+d1) = 1/3
00018  *
00019  * FOR 3D:
00020  * Problem domain should be -.05 < x < .05; -.05 < y < .05, -.1 < z < .1
00021  * Use gamma=5/3 to match Dimonte et al.
00022  * Interface is at z=0; perturbation added to Vz
00023  * Gravity acts in the z-direction.  Special reflecting boundary conditions
00024  *   added in x3 to improve hydrostatic eqm (prevents launching of weak waves)
00025  * Atwood number A = (d2-d1)/(d2+d1) = 1/2
00026  *
00027  * PRIVATE FUNCTION PROTOTYPES:
00028  * - ran2() - random number generator from NR
00029  * - reflect_ix2() - sets BCs on L-x2 (left edge) of grid used in 2D
00030  * - reflect_ox2() - sets BCs on R-x2 (right edge) of grid used in 2D
00031  * - reflect_ix3() - sets BCs on L-x3 (left edge) of grid used in 3D
00032  * - reflect_ox3() - sets BCs on R-x3 (right edge) of grid used in 3D
00033  * - grav_pot2() - gravitational potential for 2D problem (accn in Y)
00034  * - grav_pot3() - gravitational potential for 3D problem (accn in Z)
00035  *
00036  * REFERENCE: R. Liska & B. Wendroff, SIAM J. Sci. Comput., 25, 995 (2003)    */
00037 /*============================================================================*/
00038 
00039 #include <float.h>
00040 #include <math.h>
00041 #include <stdio.h>
00042 #include <stdlib.h>
00043 #include "defs.h"
00044 #include "athena.h"
00045 #include "globals.h"
00046 #include "prototypes.h"
00047 
00048 /*==============================================================================
00049  * PRIVATE FUNCTION PROTOTYPES:
00050  * ran2() - random number generator from NR
00051  * reflect_ix2() - sets BCs on L-x2 (left edge) of grid used in 2D
00052  * reflect_ox2() - sets BCs on R-x2 (right edge) of grid used in 2D
00053  * reflect_ix3() - sets BCs on L-x3 (left edge) of grid used in 3D
00054  * reflect_ox3() - sets BCs on R-x3 (right edge) of grid used in 3D
00055  * grav_pot2() - gravitational potential for 2D problem (accn in Y)
00056  * grav_pot3() - gravitational potential for 3D problem (accn in Z)
00057  *============================================================================*/
00058 
00059 static double ran2(long int *idum);
00060 static void reflect_ix2(GridS *pGrid);
00061 static void reflect_ox2(GridS *pGrid);
00062 static void reflect_ix3(GridS *pGrid);
00063 static void reflect_ox3(GridS *pGrid);
00064 static Real grav_pot2(const Real x1, const Real x2, const Real x3);
00065 static Real grav_pot3(const Real x1, const Real x2, const Real x3);
00066 
00067 /*=========================== PUBLIC FUNCTIONS ===============================*/
00068 /*----------------------------------------------------------------------------*/
00069 /* problem:  */
00070 
00071 void problem(DomainS *pDomain)
00072 {
00073   GridS *pGrid = pDomain->Grid;
00074   int i=0,j=0,k=0;
00075   int is,ie,js,je,ks,ke,iprob;
00076   long int iseed = -1;
00077   Real amp,x1,x2,x3,lx,ly,lz,rhoh,L_rot,fact;
00078 #ifdef MHD
00079   Real b0,angle;
00080 #endif
00081   int ixs, jxs, kxs;
00082 
00083   is = pGrid->is;  ie = pGrid->ie;
00084   js = pGrid->js;  je = pGrid->je;
00085   ks = pGrid->ks;  ke = pGrid->ke;
00086 
00087   lx = pDomain->RootMaxX[0] - pDomain->RootMinX[0];
00088   ly = pDomain->RootMaxX[1] - pDomain->RootMinX[1];
00089   lz = pDomain->RootMaxX[2] - pDomain->RootMinX[2];
00090 
00091 /* Ensure a different initial random seed for each process in an MPI calc. */
00092   ixs = pGrid->Disp[0];
00093   jxs = pGrid->Disp[1];
00094   kxs = pGrid->Disp[2];
00095   iseed = -1 - (ixs + pDomain->Nx[0]*(jxs + pDomain->Nx[1]*kxs));
00096 
00097 /* Read perturbation amplitude, problem switch, background density */
00098   amp = par_getd("problem","amp");
00099   iprob = par_geti("problem","iprob");
00100   rhoh  = par_getd_def("problem","rhoh",3.0);
00101 /* Distance over which field is rotated */
00102   L_rot  = par_getd_def("problem","L_rot",0.0);
00103 
00104 /* Read magnetic field strength, angle [should be in degrees, 0 is along +ve
00105  * X-axis (no rotation)] */
00106 #ifdef MHD
00107   b0 = par_getd("problem","b0");
00108   angle = par_getd("problem","angle");
00109   angle = (angle/180.)*PI;
00110 #endif
00111 
00112 /* 2D PROBLEM --------------------------------------------------------------- */
00113 /* Initialize two fluids with interface at y=0.0.  Pressure scaled to give a
00114  * sound speed of 1 at the interface in the light (lower, d=1) fluid 
00115  * Perturb V2 using single (iprob=1) or multiple (iprob=2) mode 
00116  */
00117 
00118   if (pGrid->Nx[2] == 1) {
00119   for (k=ks; k<=ke; k++) {
00120     for (j=js; j<=je; j++) {
00121       for (i=is; i<=ie; i++) {
00122         cc_pos(pGrid,i,j,k,&x1,&x2,&x3);
00123         pGrid->U[k][j][i].d = 1.0;
00124         pGrid->U[k][j][i].E = (1.0/Gamma - 0.1*x2)/Gamma_1;
00125         pGrid->U[k][j][i].M1 = 0.0;
00126         if (iprob == 1) {
00127           pGrid->U[k][j][i].M2 = amp/4.0*
00128             (1.0+cos(2.0*PI*x1/lx))*(1.0+cos(2.0*PI*x2/ly));
00129         }
00130         else {
00131           pGrid->U[k][j][i].M2 = amp*(ran2(&iseed) - 0.5)*
00132             (1.0+cos(2.0*PI*x2/ly));
00133         }
00134         pGrid->U[k][j][i].M3 = 0.0;
00135         if (x2 > 0.0) {
00136           pGrid->U[k][j][i].d = 2.0;
00137           pGrid->U[k][j][i].M2 *= 2.0;
00138           pGrid->U[k][j][i].E = (1.0/Gamma - 0.2*x2)/Gamma_1;
00139         }
00140         pGrid->U[k][j][i].E+=0.5*SQR(pGrid->U[k][j][i].M2)/pGrid->U[k][j][i].d;
00141 #ifdef MHD
00142         pGrid->B1i[k][j][i] = b0;
00143         pGrid->U[k][j][i].B1c = b0;
00144         pGrid->U[k][j][i].E += 0.5*b0*b0;
00145 #endif
00146       }
00147 #ifdef MHD
00148     pGrid->B1i[k][j][ie+1] = b0;
00149 #endif
00150     }
00151   }
00152 
00153 /* Enroll gravitational potential to give acceleration in y-direction for 2D
00154  * Use special boundary condition routines.  In 2D, gravity is in the
00155  * y-direction, so special boundary conditions needed for x2
00156 */
00157 
00158   StaticGravPot = grav_pot2;
00159   if (pDomain->Disp[1] == 0) bvals_mhd_fun(pDomain, left_x2,  reflect_ix2);
00160   if (pDomain->MaxX[1] == pDomain->RootMaxX[1])
00161     bvals_mhd_fun(pDomain, right_x2, reflect_ox2);
00162 
00163   } /* end of 2D initialization  */
00164 
00165 /* 3D PROBLEM ----------------------------------------------------------------*/
00166 /* Initialize two fluids with interface at z=0.0
00167  * Pressure scaled to give a sound speed of 1 at the interface
00168  * in the light (lower, d=1) fluid
00169  * iprob = 1 -- Perturb V3 using single mode
00170  * iprob = 2 -- Perturb V3 using multiple mode
00171  * iprob = 3 -- B in light fluid only, with multimode perturbation
00172  * iprob = 4 -- B rotated by "angle" at interface, multimode perturbation
00173  */
00174 
00175   if (pGrid->Nx[2] > 1) {
00176   for (k=ks; k<=ke; k++) {
00177     for (j=js; j<=je; j++) {
00178       for (i=is; i<=ie; i++) {
00179         cc_pos(pGrid,i,j,k,&x1,&x2,&x3);
00180         pGrid->U[k][j][i].d = 1.0;
00181         pGrid->U[k][j][i].E = (1.0/Gamma - 0.1*x3)/Gamma_1;
00182         pGrid->U[k][j][i].M1 = 0.0;
00183         pGrid->U[k][j][i].M2 = 0.0;
00184         if (iprob == 1) {
00185           pGrid->U[k][j][i].M3 = amp/8.0*(1.0+cos(2.0*PI*x1/lx))*
00186             (1.0+cos(2.0*PI*x2/ly))*(1.0+cos(2.0*PI*x3/lz));
00187         }
00188         else {
00189           pGrid->U[k][j][i].M3 = amp*(ran2(&iseed) - 0.5)*
00190             (1.0+cos(2.0*PI*x3/lz));
00191         }
00192         if (x3 > 0.0) {
00193           pGrid->U[k][j][i].d = rhoh;
00194           pGrid->U[k][j][i].M3 *= rhoh;
00195           pGrid->U[k][j][i].E = (1.0/Gamma - 0.1*rhoh*x3)/Gamma_1;
00196         }
00197         pGrid->U[k][j][i].E+=0.5*SQR(pGrid->U[k][j][i].M3)/pGrid->U[k][j][i].d;
00198 #ifdef MHD
00199         switch(iprob){
00200         case 3: /* B only in light fluid, do not add B^2 to E, total P const */
00201           if (x3 <= 0.0) {
00202             pGrid->B1i[k][j][i] = b0;
00203             if (i == ie) pGrid->B1i[k][j][ie+1] = b0;
00204             pGrid->U[k][j][i].B1c = b0;
00205           }
00206           break;
00207         case 4: /* discontinuous rotation of B by angle at interface */
00208           if (x3 <= 0.0) {
00209             pGrid->B1i[k][j][i] = b0;
00210             if (i == ie) pGrid->B1i[k][j][ie+1] = b0;
00211             pGrid->U[k][j][i].B1c = b0;
00212             pGrid->U[k][j][i].E += 0.5*b0*b0;
00213           }
00214           else {
00215             pGrid->B1i[k][j][i] = b0*cos(angle);
00216             pGrid->B2i[k][j][i] = b0*sin(angle);
00217             if (i == ie) pGrid->B1i[k][j][ie+1] = b0*cos(angle);
00218             if (j == je) pGrid->B2i[k][je+1][i] = b0*sin(angle);
00219             pGrid->U[k][j][i].B1c = b0*cos(angle);
00220             pGrid->U[k][j][i].B2c = b0*sin(angle);
00221             pGrid->U[k][j][i].E += 0.5*b0*b0;
00222           }
00223           break;
00224         case 5: /* rotation of B by angle over distance L_rot at interface */
00225           if (x3 <= (-L_rot/2.0)) {
00226             pGrid->B1i[k][j][i] = b0;
00227             if (i == ie) pGrid->B1i[k][j][ie+1] = b0;
00228             pGrid->U[k][j][i].B1c = b0;
00229             pGrid->U[k][j][i].E += 0.5*b0*b0;
00230           }
00231           else if (x3 >= (L_rot/2.0)) {
00232             pGrid->B1i[k][j][i] = b0*cos(angle);
00233             pGrid->B2i[k][j][i] = b0*sin(angle);
00234             if (i == ie) pGrid->B1i[k][j][ie+1] = b0*cos(angle);
00235             if (j == je) pGrid->B2i[k][je+1][i] = b0*sin(angle);
00236             pGrid->U[k][j][i].B1c = b0*cos(angle);
00237             pGrid->U[k][j][i].B2c = b0*sin(angle);
00238             pGrid->U[k][j][i].E += 0.5*b0*b0;
00239           }
00240           else {
00241             fact = ((L_rot/2.0)+x3)/L_rot;
00242             pGrid->B1i[k][j][i] = b0*cos(fact*angle);
00243             pGrid->B2i[k][j][i] = b0*sin(fact*angle);
00244             if (i == ie) pGrid->B1i[k][j][ie+1] = b0*cos(fact*angle);
00245             if (j == je) pGrid->B2i[k][je+1][i] = b0*sin(fact*angle);
00246             pGrid->U[k][j][i].B1c = b0*cos(fact*angle);
00247             pGrid->U[k][j][i].B2c = b0*sin(fact*angle);
00248             pGrid->U[k][j][i].E += 0.5*b0*b0;
00249           }
00250 
00251           break;
00252         default:
00253           pGrid->B1i[k][j][i] = b0;
00254           if (i == ie) pGrid->B1i[k][j][ie+1] = b0;
00255           pGrid->U[k][j][i].B1c = b0;
00256           pGrid->U[k][j][i].E += 0.5*b0*b0;
00257         }
00258 #endif
00259       }
00260     }
00261   }
00262 
00263 /* Enroll gravitational potential to give accn in z-direction for 3D
00264  * Use special boundary condition routines.  In 3D, gravity is in the
00265  * z-direction, so special boundary conditions needed for x3
00266  */
00267 
00268   StaticGravPot = grav_pot3;
00269 
00270   if (pDomain->Disp[2] == 0) bvals_mhd_fun(pDomain, left_x3,  reflect_ix3);
00271   if (pDomain->MaxX[2] == pDomain->RootMaxX[2])
00272     bvals_mhd_fun(pDomain, right_x3, reflect_ox3);
00273 
00274   } /* end of 3D initialization */
00275 
00276   return;
00277 }
00278 
00279 /*==============================================================================
00280  * PROBLEM USER FUNCTIONS:
00281  * problem_write_restart() - writes problem-specific user data to restart files
00282  * problem_read_restart()  - reads problem-specific user data from restart files
00283  * get_usr_expr()          - sets pointer to expression for special output data
00284  * get_usr_out_fun()       - returns a user defined output function pointer
00285  * get_usr_par_prop()      - returns a user defined particle selection function
00286  * Userwork_in_loop        - problem specific work IN     main loop
00287  * Userwork_after_loop     - problem specific work AFTER  main loop
00288  *----------------------------------------------------------------------------*/
00289 
00290 void problem_write_restart(MeshS *pM, FILE *fp)
00291 {
00292   return;
00293 }
00294 
00295 /*
00296  * 'problem_read_restart' must enroll special boundary value functions,
00297  *    and initialize gravity on restarts
00298  */
00299 
00300 void problem_read_restart(MeshS *pM, FILE *fp)
00301 {
00302   int nl,nd;
00303 
00304   if (pM->Nx[2] == 1) {
00305     StaticGravPot = grav_pot2;
00306     for (nl=0; nl<(pM->NLevels); nl++){
00307       for (nd=0; nd<(pM->DomainsPerLevel[nl]); nd++){
00308         bvals_mhd_fun(&(pM->Domain[nl][nd]), left_x2,  reflect_ix2);
00309         bvals_mhd_fun(&(pM->Domain[nl][nd]), right_x2, reflect_ox2);
00310       }
00311     }
00312   }
00313  
00314   if (pM->Nx[2] > 1) {
00315     StaticGravPot = grav_pot3;
00316     for (nl=0; nl<(pM->NLevels); nl++){
00317       for (nd=0; nd<(pM->DomainsPerLevel[nl]); nd++){
00318         bvals_mhd_fun(&(pM->Domain[nl][nd]), left_x3,  reflect_ix3);
00319         bvals_mhd_fun(&(pM->Domain[nl][nd]), right_x3, reflect_ox3);
00320       }
00321     }
00322   }
00323 
00324   return;
00325 }
00326 
00327 ConsFun_t get_usr_expr(const char *expr)
00328 {
00329   return NULL;
00330 }
00331 
00332 VOutFun_t get_usr_out_fun(const char *name){
00333   return NULL;
00334 }
00335 
00336 void Userwork_in_loop(MeshS *pM)
00337 {
00338 }
00339 
00340 void Userwork_after_loop(MeshS *pM)
00341 {
00342 }
00343 
00344 /*=========================== PRIVATE FUNCTIONS ==============================*/
00345 
00346 /*----------------------------------------------------------------------------*/
00347 
00348 #define IM1 2147483563
00349 #define IM2 2147483399
00350 #define AM (1.0/IM1)
00351 #define IMM1 (IM1-1)
00352 #define IA1 40014
00353 #define IA2 40692
00354 #define IQ1 53668
00355 #define IQ2 52774
00356 #define IR1 12211
00357 #define IR2 3791
00358 #define NTAB 32
00359 #define NDIV (1+IMM1/NTAB)
00360 #define RNMX (1.0-DBL_EPSILON)
00361 
00362 /*! \fn double ran2(long int *idum)
00363  *  \brief  Extracted from the Numerical Recipes in C (version 2) code.  
00364  *   Modified to use doubles instead of floats. - T. A. Gardiner - Aug. 12, 2003
00365  *   
00366  *
00367  * Long period (> 2 x 10^{18}) random number generator of L'Ecuyer
00368  * with Bays-Durham shuffle and added safeguards.  Returns a uniform
00369  * random deviate between 0.0 and 1.0 (exclusive of the endpoint
00370  * values).  Call with idum = a negative integer to initialize;
00371  * thereafter, do not alter idum between successive deviates in a
00372  * sequence.  RNMX should appriximate the largest floating point value
00373  * that is less than 1. 
00374  */
00375 
00376 double ran2(long int *idum)
00377 {
00378   int j;
00379   long int k;
00380   static long int idum2=123456789;
00381   static long int iy=0;
00382   static long int iv[NTAB];
00383   double temp;
00384 
00385   if (*idum <= 0) { /* Initialize */
00386     if (-(*idum) < 1) *idum=1; /* Be sure to prevent idum = 0 */
00387     else *idum = -(*idum);
00388     idum2=(*idum);
00389     for (j=NTAB+7;j>=0;j--) { /* Load the shuffle table (after 8 warm-ups) */
00390       k=(*idum)/IQ1;
00391       *idum=IA1*(*idum-k*IQ1)-k*IR1;
00392       if (*idum < 0) *idum += IM1;
00393       if (j < NTAB) iv[j] = *idum;
00394     }
00395     iy=iv[0];
00396   }
00397   k=(*idum)/IQ1;                 /* Start here when not initializing */
00398   *idum=IA1*(*idum-k*IQ1)-k*IR1; /* Compute idum=(IA1*idum) % IM1 without */
00399   if (*idum < 0) *idum += IM1;   /* overflows by Schrage's method */
00400   k=idum2/IQ2;
00401   idum2=IA2*(idum2-k*IQ2)-k*IR2; /* Compute idum2=(IA2*idum) % IM2 likewise */
00402   if (idum2 < 0) idum2 += IM2;
00403   j=(int)(iy/NDIV);              /* Will be in the range 0...NTAB-1 */
00404   iy=iv[j]-idum2;                /* Here idum is shuffled, idum and idum2 */
00405   iv[j] = *idum;                 /* are combined to generate output */
00406   if (iy < 1) iy += IMM1;
00407   if ((temp=AM*iy) > RNMX) return RNMX; /* No endpoint values */
00408   else return temp;
00409 }
00410 
00411 #undef IM1
00412 #undef IM2
00413 #undef AM
00414 #undef IMM1
00415 #undef IA1
00416 #undef IA2
00417 #undef IQ1
00418 #undef IQ2
00419 #undef IR1
00420 #undef IR2
00421 #undef NTAB
00422 #undef NDIV
00423 #undef RNMX
00424 
00425 /*----------------------------------------------------------------------------*/
00426 /*! \fn static void reflect_ix2(GridS *pGrid)
00427  *  \brief Special reflecting boundary functions in x2 for 2D sims
00428  */
00429 
00430 static void reflect_ix2(GridS *pGrid)
00431 {
00432   int js = pGrid->js;
00433   int ks = pGrid->ks, ke = pGrid->ke;
00434   int i,j,k,il,iu,ku; /* i-lower/upper;  k-upper */
00435 
00436   iu = pGrid->ie + nghost;
00437   il = pGrid->is - nghost;
00438 
00439   for (k=ks; k<=ke; k++) {
00440     for (j=1; j<=nghost; j++) {
00441       for (i=il; i<=iu; i++) {
00442         pGrid->U[k][js-j][i]    =  pGrid->U[k][js+(j-1)][i];
00443         pGrid->U[k][js-j][i].M2 = -pGrid->U[k][js-j][i].M2; /* reflect 2-mom. */
00444         pGrid->U[k][js-j][i].E +=  
00445           pGrid->U[k][js+(j-1)][i].d*0.1*(2*j-1)*pGrid->dx2/Gamma_1;
00446       }
00447     }
00448   }
00449 
00450 #ifdef MHD
00451   for (k=ks; k<=ke; k++) {
00452     for (j=1; j<=nghost; j++) {
00453       for (i=il; i<=iu; i++) {
00454         pGrid->B1i[k][js-j][i] = pGrid->B1i[k][js+(j-1)][i];
00455       }
00456     }
00457   }
00458 
00459   for (k=ks; k<=ke; k++) {
00460     for (j=1; j<=nghost; j++) {
00461       for (i=il; i<=iu; i++) {
00462         pGrid->B2i[k][js-j][i] = pGrid->B2i[k][js+(j-1)][i];
00463       }
00464     }
00465   }
00466 
00467   if (pGrid->Nx[2] > 1) ku=ke+1; else ku=ke;
00468   for (k=ks; k<=ku; k++) {
00469     for (j=1; j<=nghost; j++) {
00470       for (i=il; i<=iu; i++) {
00471         pGrid->B3i[k][js-j][i] = pGrid->B3i[k][js+(j-1)][i];
00472       }
00473     }
00474   }
00475 #endif
00476 
00477   return;
00478 }
00479 
00480 /*----------------------------------------------------------------------------*/
00481 /*! \fn static void reflect_ox2(GridS *pGrid)
00482  *  \brief Special reflecting boundary functions in x2 for 2D sims
00483  */
00484 
00485 static void reflect_ox2(GridS *pGrid)
00486 {
00487   int je = pGrid->je;
00488   int ks = pGrid->ks, ke = pGrid->ke, ku;
00489   int i,j,k,il,iu,jl,ju; /* i/j-lower/upper */
00490 
00491   iu = pGrid->ie + nghost;
00492   il = pGrid->is - nghost;
00493 
00494   if (pGrid->Nx[1] > 1){
00495     ju = pGrid->je + nghost;
00496     jl = pGrid->js - nghost;
00497   } else {
00498     ju = pGrid->je;
00499     jl = pGrid->js;
00500   }
00501 
00502   for (k=ks; k<=ke; k++) {
00503     for (j=1; j<=nghost; j++) {
00504       for (i=il; i<=iu; i++) {
00505         pGrid->U[k][je+j][i]    =  pGrid->U[k][je-(j-1)][i];
00506         pGrid->U[k][je+j][i].M2 = -pGrid->U[k][je+j][i].M2; /* reflect 2-mom. */
00507         pGrid->U[k][je+j][i].E -=
00508           pGrid->U[k][je-(j-1)][i].d*0.1*(2*j-1)*pGrid->dx2/Gamma_1;
00509       }
00510     }
00511   }
00512 
00513 #ifdef MHD
00514   for (k=ks; k<=ke; k++) {
00515     for (j=1; j<=nghost; j++) {
00516       for (i=il; i<=iu; i++) {
00517         pGrid->B1i[k][je+j][i] = pGrid->B1i[k][je-(j-1)][i];
00518       }
00519     }
00520   }
00521 
00522 /* j=je+1 is not set for the interface field B2i */
00523   for (k=ks; k<=ke; k++) {
00524     for (j=2; j<=nghost; j++) {
00525       for (i=il; i<=iu; i++) {
00526         pGrid->B2i[k][je+j][i] = pGrid->B2i[k][je-(j-2)][i];
00527       }
00528     }
00529   }
00530 
00531   if (pGrid->Nx[2] > 1) ku=ke+1; else ku=ke;
00532   for (k=ks; k<=ku; k++) {
00533     for (j=1; j<=nghost; j++) {
00534       for (i=il; i<=iu; i++) {
00535         pGrid->B3i[k][je+j][i] = pGrid->B3i[k][je-(j-1)][i];
00536       }
00537     }
00538   }
00539 #endif
00540 
00541   return;
00542 }
00543 
00544 /*----------------------------------------------------------------------------*/
00545 /*! \fn static void reflect_ix3(GridS *pGrid)
00546  *  \brief Special reflecting boundary functions in x3 for 2D sims
00547  */
00548 
00549 static void reflect_ix3(GridS *pGrid)
00550 {
00551   int ks = pGrid->ks;
00552   int i,j,k,il,iu,jl,ju; /* i-lower/upper;  j-lower/upper */
00553 
00554   iu = pGrid->ie + nghost;
00555   il = pGrid->is - nghost;
00556   if (pGrid->Nx[1] > 1){
00557     ju = pGrid->je + nghost;
00558     jl = pGrid->js - nghost;
00559   } else {
00560     ju = pGrid->je;
00561     jl = pGrid->js;
00562   }
00563 
00564   for (k=1; k<=nghost; k++) {
00565     for (j=jl; j<=ju; j++) {
00566       for (i=il; i<=iu; i++) {
00567         pGrid->U[ks-k][j][i]    =  pGrid->U[ks+(k-1)][j][i];
00568         pGrid->U[ks-k][j][i].M3 = -pGrid->U[ks-k][j][i].M3; /* reflect 3-mom. */
00569         pGrid->U[ks-k][j][i].E +=
00570           pGrid->U[ks+(k-1)][j][i].d*0.1*(2*k-1)*pGrid->dx3/Gamma_1;
00571       }
00572     }
00573   }
00574 
00575 #ifdef MHD
00576   for (k=1; k<=nghost; k++) {
00577     for (j=jl; j<=ju; j++) {
00578       for (i=il; i<=iu; i++) {
00579         pGrid->B1i[ks-k][j][i] = pGrid->B1i[ks+(k-1)][j][i];
00580       }
00581     }
00582   }
00583   for (k=1; k<=nghost; k++) {
00584     for (j=jl; j<=ju; j++) {
00585       for (i=il; i<=iu; i++) {
00586         pGrid->B2i[ks-k][j][i] = pGrid->B2i[ks+(k-1)][j][i];
00587       }
00588     }
00589   }
00590 
00591   for (k=1; k<=nghost; k++) {
00592     for (j=jl; j<=ju; j++) {
00593       for (i=il; i<=iu; i++) {
00594         pGrid->B3i[ks-k][j][i] = pGrid->B3i[ks+(k-1)][j][i];
00595       }
00596     }
00597   }
00598 #endif
00599 
00600   return;
00601 }
00602 
00603 /*----------------------------------------------------------------------------*/
00604 /*! \fn static void reflect_ox3(GridS *pGrid)
00605  *  \brief Special reflecting boundary functions in x3 for 3D sims
00606  */
00607 
00608 static void reflect_ox3(GridS *pGrid)
00609 {
00610   int ke = pGrid->ke;
00611   int i,j,k ,il,iu,jl,ju; /* i-lower/upper;  j-lower/upper */
00612 
00613   iu = pGrid->ie + nghost;
00614   il = pGrid->is - nghost;
00615   if (pGrid->Nx[1] > 1){
00616     ju = pGrid->je + nghost;
00617     jl = pGrid->js - nghost;
00618   } else {
00619     ju = pGrid->je;
00620     jl = pGrid->js;
00621   }
00622   for (k=1; k<=nghost; k++) {
00623     for (j=jl; j<=ju; j++) {
00624       for (i=il; i<=iu; i++) {
00625         pGrid->U[ke+k][j][i]    =  pGrid->U[ke-(k-1)][j][i];
00626         pGrid->U[ke+k][j][i].M3 = -pGrid->U[ke+k][j][i].M3; /* reflect 3-mom. */
00627         pGrid->U[ke+k][j][i].E -=
00628           pGrid->U[ke-(k-1)][j][i].d*0.1*(2*k-1)*pGrid->dx3/Gamma_1;
00629       }
00630     }
00631   }
00632 
00633 #ifdef MHD
00634   for (k=1; k<=nghost; k++) {
00635     for (j=jl; j<=ju; j++) {
00636       for (i=il; i<=iu; i++) {
00637         pGrid->B1i[ke+k][j][i] = pGrid->B1i[ke-(k-1)][j][i];
00638       }
00639     }
00640   }
00641 
00642   for (k=1; k<=nghost; k++) {
00643     for (j=jl; j<=ju; j++) {
00644       for (i=il; i<=iu; i++) {
00645         pGrid->B2i[ke+k][j][i] = pGrid->B2i[ke-(k-1)][j][i];
00646       }
00647     }
00648   }
00649 
00650 /* Note that k=ke+1 is not a boundary condition for the interface field B3i */
00651   for (k=2; k<=nghost; k++) {
00652     for (j=jl; j<=ju; j++) {
00653       for (i=il; i<=iu; i++) {
00654         pGrid->B3i[ke+k][j][i] = pGrid->B3i[ke-(k-1)][j][i];
00655       }
00656     }
00657   }
00658 #endif
00659 
00660   return;
00661 }
00662 
00663 /*----------------------------------------------------------------------------*/
00664 /*! \fn static Real grav_pot2(const Real x1, const Real x2, const Real x3)
00665  *  \brief Gravitational potential; g = 0.1
00666  */
00667 
00668 static Real grav_pot2(const Real x1, const Real x2, const Real x3)
00669 {
00670   return 0.1*x2;
00671 }
00672 /*! \fn static Real grav_pot3(const Real x1, const Real x2, const Real x3)
00673  *  \brief Gravitational potential; g = 0.1
00674  */
00675 static Real grav_pot3(const Real x1, const Real x2, const Real x3)
00676 {
00677   return 0.1*x3;
00678 }