Athena: prob/ti.c Source File

Go to the documentation of this file.
00001 #include "copyright.h"
00002 /*============================================================================*/
00003 /*! \file ti.c
00004  *  \brief Problem generator for thermal instability.
00005  *
00006  *  PURPOSE: Problem generator for thermal instability.  Setup to use the 
00007  *   Koyama & Inutsuka cooling function, so must use cgs units for calculations
00008  *
00009  *  TI test w/ Isotropic Conduction - hydro
00010  *   The coefficient of conductivity kappa should be defined in cgs unit
00011  *
00012  *  TI test w/ Anisotropic Conduction - mhd
00013  *   The coefficient of Spitzer Coulombic conductivity should be given.
00014  *
00015  *  Various test cases are possible
00016  *- (iprob=1) : TI test w/ isotropic conduction in 1D - sinusoidal fluctuations
00017  *- (iprob=2) : TI test w/ isotropic conduction in 2D
00018  *                       with sinusoidal fluctuations along the diagonal        
00019  *- (iprob=3) : TI test w/ conduction in 2D with random pressure purtubation
00020  *              in case of anisotropic conduction : mhd - Diagonal field line
00021  *- (iprob=4) : TI test w/ conduction in 2D with random pressure purtubation
00022  *              in case of anisotropic conduction : mhd - Tangled field line  */
00023 /*============================================================================*/
00024 
00025 #include <float.h>
00026 #include <math.h>
00027 #include <stdio.h>
00028 #include <stdlib.h>
00029 #include "defs.h"
00030 #include "athena.h"
00031 #include "globals.h"
00032 #include "prototypes.h"
00033 
00034 /* These constants are in cgs */
00035 static const Real mbar = (1.27)*(1.6733e-24);
00036 static const Real kb = 1.380658e-16;
00037 static const Real pi = 3.14159265358979323846;
00038 
00039 /*==============================================================================
00040  * PRIVATE FUNCTION PROTOTYPES:
00041  * ran2()
00042  *============================================================================*/
00043 
00044 static double ran2(long int *idum);
00045 static Real logd(const Grid *pG, const int i, const int j, const int k);
00046 
00047 /*=========================== PUBLIC FUNCTIONS ===============================*/
00048 /*----------------------------------------------------------------------------*/
00049 /* problem:  */
00050 
00051 void problem(Grid *pGrid, Domain *pDomain)
00052 {
00053   int i=0,j=0,k=0,n=0,m=0,l=0;
00054   int is,ie,js,je,ks,ke,iprob,
00055       nkx=16,nky=8,nkz=8,nx1,nx2,nx3;
00056   Real P_k, n0, T0, kappa,chi, x1L, x2L, x3L, x12L, num,
00057                          drho, dp,
00058        krho, krho_i, krho_ij, angle, phi, phi2, phi3, ksi, tmp,
00059                          rho0, rho1, n_anal, r, r_ij, r_proj, 
00060        beta, b0, Btot=0.0,Bsum=0.0,
00061        ***az,***ax,***ay,amp,kx[16],ky[8],kz[8],xa,ya,za;
00062   long int iseed = -1, rseed;
00063 
00064   is = pGrid->is; ie = pGrid->ie;
00065   js = pGrid->js; je = pGrid->je;
00066   ks = pGrid->ks; ke = pGrid->ke;
00067 
00068 /* Read problem parameters */
00069   P_k   = par_getd("problem","P_k");
00070   n0    = par_getd("problem","n0");
00071   x1L           = par_getd("grid","x1max");
00072         x2L   = par_getd("grid","x2max");
00073   x3L   = par_getd("grid","x3max");
00074         iprob = par_geti("problem","iprob");
00075 
00076 #ifdef MHD
00077         beta  = par_getd("problem","beta");
00078 //  b0    = par_getd("problem","b0");
00079         b0 = sqrt(2.0*P_k*kb/beta);
00080 #endif
00081 #ifdef ISOTROPIC_CONDUCTION
00082   kappa = par_getd("problem","kappa");
00083   kappa_T = (mbar/kb)*kappa;
00084 #endif
00085 #ifdef ANISOTROPIC_CONDUCTION
00086   chi = par_getd("problem","chi");
00087   chi_C = (mbar/kb)*chi;
00088 #endif
00089 
00090 
00091 /* Set up Initial Temperature and density */
00092         T0 = P_k / n0;
00093   rho0 = n0 * mbar;
00094 
00095 /* For (iprob=1) -- TI w/ isotropic conduction - growth rate test in 1D
00096  *   2 with sinusoidal fluctuations in density, pressure and velocity  
00097  *    drho (initial fluctuation amplitude in density) should be given. 
00098  *    krho = 2pi/x1L * n */
00099 if(iprob==1) {
00100   drho  = par_getd("problem","drho");
00101   n_anal = par_getd("problem","n_anal"); // analytic growth rate - Field 1965
00102   num   = par_getd("problem","num"); // wavenumber n
00103 
00104   krho_i = 2.0 * pi * num / (float)(ie-is+1);
00105   krho = 2.0 * pi * num / x1L;
00106         rho1 = rho0 * drho;
00107   for (k=ks; k<=ke; k++) {
00108   for (j=js; j<=je; j++) {
00109     for (i=is; i<=ie; i++) {
00110       pGrid->U[k][j][i].d = rho0 + rho1*cos(krho_i*(float)(i-is));
00111       pGrid->U[k][j][i].M1 = pGrid->U[k][j][i].d * n_anal * rho1 * sin(krho_i*(float)(i-is)) / rho0 / krho * (-1.0);
00112       pGrid->U[k][j][i].M2 = 0.0;
00113       pGrid->U[k][j][i].M3 = 0.0;
00114       pGrid->U[k][j][i].E = (n0*kb*T0)/Gamma_1 - n_anal*n_anal/ (krho * krho * Gamma_1) * rho1*cos(krho_i*(float)(i-is))
00115              + 0.5*(SQR(pGrid->U[k][j][i].M1) + SQR(pGrid->U[k][j][i].M2)
00116              + SQR(pGrid->U[k][j][i].M3))/pGrid->U[k][j][i].d;
00117        }
00118      }
00119    } 
00120  }
00121 
00122 
00123 /* For (iprob=2) -- TI w/ isotropic conduction - growth rate test in 2D
00124  *    with sinusoidal fluctuations in density, pressure and velocity
00125  *    drho (initial fluctuation amplitude in density) should be given.
00126  *    krho = 2pi/x12L * n  - x12L (diagonal length of the 2D box)      */
00127 if(iprob==2) {
00128         angle = atan(x2L/x1L);
00129         x12L = pow(x1L*x1L + x2L*x2L, 0.5);     
00130         r_ij = pow(pow((float)(ie-is),2)+pow((float)(je-js),2),0.5);
00131         krho = 2.0 * pi * num / x12L;
00132         krho_ij = 2.0 * pi * num / r_ij;
00133         rho1 = rho0 * drho;
00134   for (k=ks; k<=ke; k++) {
00135   for (j=js; j<=je; j++) {
00136     for (i=is; i<=ie; i++) {
00137       r = pow(pow((float)(i-is),2)+pow((float)(j-js),2),0.5);
00138       if ((i-is) == 0) {
00139       phi = 0;
00140       }
00141       else{
00142       phi = atan((float)(j-js)/(float)(i-is));
00143       }
00144       r_proj = r * sin(phi + angle) / sin ( 2 * angle);
00145       pGrid->U[k][j][i].d = rho0 + rho1 * cos(krho_ij * r_proj);
00146       pGrid->U[k][j][i].M1 = pGrid->U[k][j][i].d * n_anal * rho1 * sin(krho_ij * r_proj) / rho0 / krho * cos(angle) * (-1.0);
00147       pGrid->U[k][j][i].M2 = pGrid->U[k][j][i].d * n_anal * rho1 * sin(krho_ij * r_proj) / rho0 / krho * sin(angle) * (-1.0);
00148       pGrid->U[k][j][i].M3 = 0.0;
00149       pGrid->U[k][j][i].E = (n0*kb*T0)/Gamma_1 
00150              - n_anal*n_anal/ (krho * krho * Gamma_1) * rho1 * cos(krho_ij * r_proj)
00151              + 0.5*(SQR(pGrid->U[k][j][i].M1) + SQR(pGrid->U[k][j][i].M2)
00152              + SQR(pGrid->U[k][j][i].M3))/pGrid->U[k][j][i].d;
00153        }
00154      }
00155    }
00156  }
00157 
00158 /* For (iprob=3) -- TI w/ (an)isotropic conduction - growth rate test in 2D
00159  *    with RANDOM fluctuations in pressure
00160  *    dp (initial fluctuation amplitude in pressure) should be given.
00161  *    iprob=3 - if MHD - diagonal streight field line / 
00162                 if HYDRO - isotropic conduction test with random perturbation in pressure  */
00163  if(iprob==3) {
00164   dp  = par_getd("problem","dp");
00165   angle = atan(x2L/x1L);
00166   for (k=ks; k<=ke; k++) {
00167   for (j=js; j<=je; j++) {
00168     for (i=is; i<=ie; i++) {
00169       pGrid->U[k][j][i].d = rho0; 
00170       pGrid->U[k][j][i].M1 = 0.0; 
00171       pGrid->U[k][j][i].M2 = 0.0;  
00172       pGrid->U[k][j][i].M3 = 0.0;
00173       pGrid->U[k][j][i].E = (n0*kb*T0)/Gamma_1 + (n0*kb*T0)/Gamma_1 * (ran2(&rseed) - 0.5) * dp
00174              + 0.5*(SQR(pGrid->U[k][j][i].M1) + SQR(pGrid->U[k][j][i].M2)
00175              + SQR(pGrid->U[k][j][i].M3))/pGrid->U[k][j][i].d;
00176 #ifdef MHD
00177       pGrid->B1i[k][j][i] = b0 * cos(angle);
00178       if (i == ie) pGrid->B1i[k][j][i+1] = b0 * cos(angle);
00179       pGrid->U[k][j][i].B1c = b0 * cos(angle);
00180 
00181       pGrid->B2i[k][j][i] = b0 * sin(angle);
00182       if (j==je) pGrid->B2i[k][j+1][i] = b0 * sin(angle);
00183       pGrid->U[k][j][i].B2c = b0 * sin(angle);
00184 
00185       pGrid->B3i[k][j][i] = 0.0;
00186       if (k==ke && pGrid->Nx3 > 1) pGrid->B3i[k+1][j][i] = 0.0;
00187       pGrid->U[k][j][i].B3c = 0.0;
00188       pGrid->U[k][j][i].E += 0.5*(SQR(pGrid->U[k][j][i].B1c)
00189                            + SQR(pGrid->U[k][j][i].B2c) + SQR(pGrid->U[k][j][i].B3c));
00190 #endif /* MHD */
00191       }
00192     }
00193    }
00194 }
00195 
00196 /* For (iprob=4) -- TI w/ (an)isotropic conduction - growth rate test in 2D
00197  *    with RANDOM fluctuations in pressure
00198  *    dp (initial fluctuation amplitude in pressure) should be given.
00199  *    iprob=4 - tangled field line */
00200 if(iprob==4){
00201 
00202   nx1 = (ie-is)+1 + 2*nghost;
00203   nx2 = (je-js)+1 + 2*nghost;
00204   if ((nx1 == 1) || (nx2 == 1)) {
00205     ath_error("[TI-tangledB]: This problem can only be run with Nx1>1\n");
00206   }
00207   if ((az = (Real***)calloc_3d_array(nx3,nx2, nx1, sizeof(Real))) == NULL) {
00208     ath_error("[TI-tangledB]: Error allocating memory for vector pot\n");
00209   }
00210 
00211 /* Vector Potential Generation*/
00212   for (n=0; n<=nkx-1; n++){
00213     kx[n] = (1.+(float)(n))*(2.0*pi)/2.0;
00214   }
00215         kx[nkx-1] = kx[0] * (-1.0);
00216   for (m=0; m<=nky-1; m++){
00217     ky[m] = (1.+(float)(m))*(2.0*pi);
00218   }
00219   ky[nky-1] = ky[0] * (-1.0);
00220 
00221   for (k=ks; k=ks; k++) {
00222   for (j=js; j<=je+1; j++) {
00223     for (i=is; i<=ie+1; i++) {
00224       az[k][j][i] = 0.0;
00225      }
00226    }
00227   }
00228   for (n=0; n<=nkx-1; n++){
00229     for (m=0; m<=nky-1; m++){
00230       phi = 2.0 * pi * ran2(&rseed);
00231 //                      printf("phi=%f",phi);
00232       amp = (1.0/sqrt( kx[n]*kx[n] + ky[m]*ky[m]));
00233      for (j=js; j<=je+1; j++) {
00234        ya=(float)(j-js)*(1.0/(float)(je-js+1));
00235 /*                       if(j == js | j==je+1){
00236                                 printf("ya=%f\n",ya);
00237        }*/
00238       for (i=is; i<=ie+1; i++) {
00239        xa=(float)(i-is)*(2.0/(float)(ie-is+1));
00240         az[ks][j][i] += amp*sin(kx[n]*xa + ky[m]*ya + phi);
00241        }
00242      }
00243     }
00244   }
00245 /*Vector Potential for 2D generation end*/
00246 
00247 #ifdef MHD
00248   dp  = par_getd("problem","dp");
00249   for (k=ks; k<=ke; k++) {
00250   for (j=js; j<=je; j++) {
00251     for (i=is; i<=ie; i++) {
00252       pGrid->U[k][j][i].d = rho0;
00253       pGrid->U[k][j][i].M1 = 0.0;
00254       pGrid->U[k][j][i].M2 = 0.0;
00255       pGrid->U[k][j][i].M3 = 0.0;
00256       pGrid->U[k][j][i].E = (n0*kb*T0)/Gamma_1 + (n0*kb*T0)/Gamma_1 * (ran2(&rseed) - 0.5) * dp
00257              + 0.5*(SQR(pGrid->U[k][j][i].M1) + SQR(pGrid->U[k][j][i].M2)
00258              + SQR(pGrid->U[k][j][i].M3))/pGrid->U[k][j][i].d;
00259       pGrid->B1i[k][j][i] = (az[ks][j+1][i] - az[ks][j][i]);
00260       pGrid->B2i[k][j][i] =-(az[ks][j][i+1] - az[ks][j][i]);
00261 
00262       pGrid->B3i[k][j][i] = 0.0;
00263       if (k==ke && pGrid->Nx3 > 1) pGrid->B3i[k+1][j][i] = 0.0;
00264 
00265       Bsum += (pGrid->B1i[k][j][i]*pGrid->B1i[k][j][i] + 
00266              pGrid->B2i[k][j][i]*pGrid->B2i[k][j][i] + pGrid->B3i[k][j][i]*pGrid->B3i[k][j][i]);
00267        }
00268      }
00269    }
00270         Btot = sqrt(Bsum/((float)(je-js+1) * (float)(ie-is+1)));
00271   printf("Btot=%f\n",Btot);
00272         Bsum = 0.0;
00273 
00274   for (k=ks; k<=ke; k++) {
00275   for (j=js; j<=je; j++) {
00276     for (i=is; i<=ie; i++) {
00277       pGrid->B1i[k][j][i] = pGrid->B1i[k][j][i] * b0 / Btot;
00278       if (i == ie) pGrid->B1i[k][j][i+1] = pGrid->B1i[k][j][is];
00279       pGrid->B2i[k][j][i] = pGrid->B2i[k][j][i] * b0 / Btot;
00280       if (j == je) pGrid->B2i[k][j+1][i] = pGrid->B2i[k][js][i];
00281       Bsum += (pGrid->B1i[k][j][i]*pGrid->B1i[k][j][i] +
00282              pGrid->B2i[k][j][i]*pGrid->B2i[k][j][i] + pGrid->B3i[k][j][i]*pGrid->B3i[k][j][i]);
00283        }
00284      }
00285    }
00286 
00287   for (k=ks; k<=ke; k++) {
00288   for (j=js; j<=je; j++) {
00289     for (i=is; i<=ie; i++) {
00290       pGrid->U[k][j][i].B1c = 0.5*(pGrid->B1i[k][j][i]+pGrid->B1i[k][j][i+1]);
00291       pGrid->U[k][j][i].B2c = 0.5*(pGrid->B2i[k][j][i]+pGrid->B2i[k][j+1][i]);
00292       pGrid->U[k][j][i].B3c = 0.0;
00293       pGrid->U[k][j][i].E += 0.5*(SQR(pGrid->U[k][j][i].B1c)
00294                            + SQR(pGrid->U[k][j][i].B2c) + SQR(pGrid->U[k][j][i].B3c));
00295        }
00296      }
00297    }
00298   Btot = Bsum/((float)(je-js+1) * (float)(ie-is+1));
00299   printf("B^2tot,b0^2=%e  %e \n",Btot,b0*b0);
00300 #endif //MHD
00301  }
00302 
00303 /* For (iprob=5) -- TI w/ (an)isotropic conduction - growth rate test in 3D
00304  *    with RANDOM fluctuations in pressure
00305  *    dp (initial fluctuation amplitude in pressure) should be given.
00306  *    iprob=5 - if MHD - diagonal streight field line /
00307                 if HYDRO - isotropic conduction test with random perturbation in pressure  */
00308  if(iprob==5) {
00309   dp  = par_getd("problem","dp");
00310   angle = atan(x2L/x1L);
00311   x12L = pow(x1L*x1L + x2L*x2L, 0.5);
00312         ksi = atan(x3L/x12L);
00313   for (k=ks; k<=ke; k++) {
00314   for (j=js; j<=je; j++) {
00315     for (i=is; i<=ie; i++) {
00316       pGrid->U[k][j][i].d = rho0;
00317       pGrid->U[k][j][i].M1 = 0.0;
00318       pGrid->U[k][j][i].M2 = 0.0;
00319       pGrid->U[k][j][i].M3 = 0.0;
00320       pGrid->U[k][j][i].E = (n0*kb*T0)/Gamma_1 + (n0*kb*T0)/Gamma_1 * (ran2(&rseed) - 0.5) * dp
00321              + 0.5*(SQR(pGrid->U[k][j][i].M1) + SQR(pGrid->U[k][j][i].M2)
00322              + SQR(pGrid->U[k][j][i].M3))/pGrid->U[k][j][i].d;
00323 #ifdef MHD
00324       pGrid->B1i[k][j][i] = b0 * cos(angle) * cos(ksi);
00325       if (i == ie) pGrid->B1i[k][j][i+1] = b0 * cos(angle) * cos(ksi);
00326       pGrid->U[k][j][i].B1c = b0 * cos(angle) * cos(ksi);
00327 
00328       pGrid->B2i[k][j][i] = b0 * sin(angle) * cos(ksi);
00329       if (j==je) pGrid->B2i[k][j+1][i] = b0 * sin(angle) * cos(ksi);
00330       pGrid->U[k][j][i].B2c = b0 * sin(angle) * cos(ksi);
00331 
00332       pGrid->B3i[k][j][i] = b0 * sin(ksi);
00333       if (k==ke && pGrid->Nx3 > 1) pGrid->B3i[k+1][j][i] = b0 * sin(ksi);
00334       pGrid->U[k][j][i].B3c = b0 * sin(ksi);
00335       pGrid->U[k][j][i].E += 0.5*(SQR(pGrid->U[k][j][i].B1c)
00336                            + SQR(pGrid->U[k][j][i].B2c) + SQR(pGrid->U[k][j][i].B3c));
00337 #endif /* MHD */
00338       }
00339     }
00340    }
00341 }
00342 
00343 
00344 /* For (iprob=6) -- TI w/ (an)isotropic conduction - growth rate test in *3D*
00345  *    with RANDOM fluctuations in pressure
00346  *    dp (initial fluctuation amplitude in pressure) should be given.
00347  *    iprob=6 - in mhd, tangled field line */
00348 if(iprob==6){
00349   nx1 = (ie-is)+1 + 2*nghost;
00350   nx2 = (je-js)+1 + 2*nghost;
00351   nx3 = (ke-ks)+1 + 2*nghost;
00352   if ((nx1 == 1) || (nx2 == 1) || (nx3 == 1)) {
00353     ath_error("[TI-tangledB]: This problem can only be run with Nx1>1\n");
00354   }
00355   if ((ay = (Real***)calloc_3d_array(nx3, nx2, nx1, sizeof(Real))) == NULL) {
00356     ath_error("[TI-tangledB]: Error allocating memory for vector pot\n");
00357   }
00358   if ((az = (Real***)calloc_3d_array(nx3, nx2, nx1, sizeof(Real))) == NULL) {
00359     ath_error("[TI-tangledB]: Error allocating memory for vector pot\n");
00360   }
00361   if ((ax = (Real***)calloc_3d_array(nx3, nx2, nx1, sizeof(Real))) == NULL) {
00362     ath_error("[TI-tangledB]: Error allocating memory for vector pot\n");
00363   }
00364 
00365 /* Vector Potential Generation*/
00366   for (n=0; n<=nkx-1; n++){
00367     kx[n] = (1.+(float)(n))*(2.0*pi)/2.0;
00368   }
00369   kx[nkx-1] = kx[0] * (-1.0);
00370   for (m=0; m<=nky-1; m++){
00371     ky[m] = (1.+(float)(m))*(2.0*pi);
00372   }
00373   ky[nky-1] = ky[0] * (-1.0);
00374   for (l=0; l<=nkz-1; l++){
00375     kz[l] = (1.+(float)(l))*(2.0*pi);
00376   }
00377   kz[nkz-1] = kz[0] * (-1.0);
00378 
00379   for (k=ks; k<=ke+1; k++) {
00380    for (j=js; j<=je+1; j++) {
00381     for (i=is; i<=ie+1; i++) {
00382       ax[k][j][i] = 0.0;
00383       ay[k][j][i] = 0.0;
00384       az[k][j][i] = 0.0;
00385     }
00386    }
00387   }
00388 
00389   for (n=0; n<=nkx-1; n++){
00390    for (m=0; m<=nky-1; m++){
00391     for (l=0; l<=nkz-1; l++){
00392       phi = 2.0 * pi * ran2(&rseed);
00393       phi2 = 2.0 * pi * ran2(&rseed);
00394       phi3 = 2.0 * pi * ran2(&rseed);
00395       amp = (1.0/sqrt(kx[n]*kx[n] + ky[m]*ky[m] + kz[l]*kz[l]));
00396      for (k=ks; k<=ke+1; k++) {
00397        za=(float)(k-ks)*(1.0/(float)(ke-ks+1));
00398      for (j=js; j<=je+1; j++) {
00399        ya=(float)(j-js)*(1.0/(float)(je-js+1));
00400       for (i=is; i<=ie+1; i++) {
00401        xa=(float)(i-is)*(2.0/(float)(ie-is+1));
00402        tmp = kx[n]*xa + ky[m]*ya + kz[l]*za;
00403         az[k][j][i] += amp*sin(tmp + phi);
00404         ay[k][j][i] += amp*sin(tmp + phi2);
00405         ax[k][j][i] += amp*sin(tmp + phi3);
00406        }
00407       }
00408      }
00409     }
00410    }
00411   }
00412 
00413   printf("za=%f\n",za);
00414   Bsum = 0.0;
00415 #ifdef MHD
00416   dp  = par_getd("problem","dp");
00417   for (k=ks; k<=ke; k++) {
00418   for (j=js; j<=je; j++) {
00419     for (i=is; i<=ie; i++) {
00420       pGrid->U[k][j][i].d = rho0;
00421       pGrid->U[k][j][i].M1 = 0.0;
00422       pGrid->U[k][j][i].M2 = 0.0;
00423       pGrid->U[k][j][i].M3 = 0.0;
00424       pGrid->U[k][j][i].E = (n0*kb*T0)/Gamma_1 + (n0*kb*T0)/Gamma_1 * (ran2(&rseed) - 0.5) * dp
00425              + 0.5*(SQR(pGrid->U[k][j][i].M1) + SQR(pGrid->U[k][j][i].M2)
00426              + SQR(pGrid->U[k][j][i].M3))/pGrid->U[k][j][i].d;
00427       pGrid->B1i[k][j][i] = (az[k][j+1][i] - az[k][j][i]) -
00428                             (ay[k+1][j][i] - ay[k][j][i]);
00429       pGrid->B2i[k][j][i] = (ax[k+1][j][i] - ax[k][j][i]) -
00430                             (az[k][j][i+1] - az[k][j][i]);
00431       pGrid->B3i[k][j][i] = (ay[k][j][i+1] - ay[k][j][i]) -
00432                             (ax[k][j+1][i] - ax[k][j][i]);
00433 
00434       Bsum += sqrt(pGrid->B1i[k][j][i]*pGrid->B1i[k][j][i] +
00435              pGrid->B2i[k][j][i]*pGrid->B2i[k][j][i] + pGrid->B3i[k][j][i]*pGrid->B3i[k][j][i]);
00436        }
00437      }
00438    }
00439   Btot = Bsum/((float)(je-js+1) * (float)(ie-is+1) * (float)(ke-ks+1));
00440   printf("Btot=%f\n",Btot);
00441   Bsum = 0.0;
00442   for (k=ks; k<=ke; k++) {
00443   for (j=js; j<=je; j++) {
00444     for (i=is; i<=ie; i++) {
00445       pGrid->B1i[k][j][i] = pGrid->B1i[k][j][i] * b0 / Btot;
00446       if (i == ie) pGrid->B1i[k][j][i+1] = pGrid->B1i[k][j][is];
00447       pGrid->B2i[k][j][i] = pGrid->B2i[k][j][i] * b0 / Btot;
00448       if (j == je) pGrid->B2i[k][j+1][i] = pGrid->B2i[k][js][i];
00449       pGrid->B3i[k][j][i] = pGrid->B3i[k][j][i] * b0 / Btot;
00450       if (k == ke) pGrid->B3i[k+1][j][i] = pGrid->B3i[ks][j][i];
00451 
00452       Bsum += sqrt(pGrid->B1i[k][j][i]*pGrid->B1i[k][j][i] +
00453              pGrid->B2i[k][j][i]*pGrid->B2i[k][j][i] + pGrid->B3i[k][j][i]*pGrid->B3i[k][j][i]);
00454        }
00455      }
00456    }
00457 
00458   for (k=ks; k<=ke; k++) {
00459   for (j=js; j<=je; j++) {
00460     for (i=is; i<=ie; i++) {
00461       pGrid->U[k][j][i].B1c = 0.5*(pGrid->B1i[k][j][i]+pGrid->B1i[k][j][i+1]);
00462       pGrid->U[k][j][i].B2c = 0.5*(pGrid->B2i[k][j][i]+pGrid->B2i[k][j+1][i]);
00463       pGrid->U[k][j][i].B3c = 0.5*(pGrid->B3i[k][j][i]+pGrid->B3i[k+1][j][i]);
00464       pGrid->U[k][j][i].E += 0.5*(SQR(pGrid->U[k][j][i].B1c)
00465                            + SQR(pGrid->U[k][j][i].B2c) + SQR(pGrid->U[k][j][i].B3c));
00466        }
00467      }
00468    }
00469   Btot = Bsum/((float)(je-js+1) * (float)(ie-is+1) * (float)(ke-ks+1));
00470   printf("Btot/b0=%e  %e\n",Btot,b0);
00471 #endif //MHD
00472  }
00473 //printf("Btot=%f\n",xa);
00474   CoolingFunc = KoyInut;
00475 }
00476 
00477 /*==============================================================================
00478  * PROBLEM USER FUNCTIONS:
00479  * problem_write_restart() - writes problem-specific user data to restart files
00480  * problem_read_restart()  - reads problem-specific user data from restart files
00481  * get_usr_expr()          - sets pointer to expression for special output data
00482  * get_usr_out_fun()       - returns a user defined output function pointer
00483  * Userwork_in_loop        - problem specific work IN     main loop
00484  * Userwork_after_loop     - problem specific work AFTER  main loop
00485  *----------------------------------------------------------------------------*/
00486 
00487 void problem_write_restart(Grid *pG, Domain *pD, FILE *fp)
00488 {
00489   return;
00490 }
00491 
00492 void problem_read_restart(Grid *pG, Domain *pD, FILE *fp)
00493 {
00494   return;
00495 }
00496 
00497 /*! \fn static Real logd(const Grid *pG, const int i, const int j, const int k)
00498  *  \brief Log10 of density */
00499 static Real logd(const Grid *pG, const int i, const int j, const int k)
00500 {
00501   return log10(pG->U[k][j][i].d);
00502 }
00503 
00504 Gasfun_t get_usr_expr(const char *expr)
00505 {
00506   if(strcmp(expr,"logd")==0) return logd;
00507   return NULL;
00508 }
00509 
00510 VGFunout_t get_usr_out_fun(const char *name){
00511   return NULL;
00512 }
00513 
00514 
00515 void Userwork_in_loop(Grid *pGrid, Domain *pDomain)
00516 {
00517   return;
00518 }
00519 
00520 void Userwork_after_loop(Grid *pGrid, Domain *pDomain)
00521 {
00522   return;
00523 }
00524 
00525 /*=========================== PRIVATE FUNCTIONS ==============================*/
00526 /*------------------------------------------------------------------------------
00527  */
00528 
00529 #define IM1 2147483563
00530 #define IM2 2147483399
00531 #define AM (1.0/IM1)
00532 #define IMM1 (IM1-1)
00533 #define IA1 40014
00534 #define IA2 40692
00535 #define IQ1 53668
00536 #define IQ2 52774
00537 #define IR1 12211
00538 #define IR2 3791
00539 #define NTAB 32
00540 #define NDIV (1+IMM1/NTAB)
00541 #define RNMX (1.0-DBL_EPSILON)
00542 
00543 /*! \fn double ran2(long int *idum){
00544  *  \brief Extracted from the Numerical Recipes in C (version 2) code.  Modified
00545  *   to use doubles instead of floats. -- T. A. Gardiner -- Aug. 12, 2003
00546  *
00547  * Long period (> 2 x 10^{18}) random number generator of L'Ecuyer
00548  * with Bays-Durham shuffle and added safeguards.  Returns a uniform
00549  * random deviate between 0.0 and 1.0 (exclusive of the endpoint
00550  * values).  Call with idum = a negative integer to initialize;
00551  * thereafter, do not alter idum between successive deviates in a
00552  * sequence.  RNMX should appriximate the largest floating point value
00553  * that is less than 1. 
00554 
00555  */
00556 double ran2(long int *idum){
00557   int j;
00558   long int k;
00559   static long int idum2=123456789;
00560   static long int iy=0;
00561   static long int iv[NTAB];
00562   double temp;
00563 
00564   if (*idum <= 0) { /* Initialize */
00565     if (-(*idum) < 1) *idum=1; /* Be sure to prevent idum = 0 */
00566     else *idum = -(*idum);
00567     idum2=(*idum);
00568     for (j=NTAB+7;j>=0;j--) { /* Load the shuffle table (after 8 warm-ups) */
00569       k=(*idum)/IQ1;
00570       *idum=IA1*(*idum-k*IQ1)-k*IR1;
00571       if (*idum < 0) *idum += IM1;
00572       if (j < NTAB) iv[j] = *idum;
00573     }
00574     iy=iv[0];
00575   }
00576   k=(*idum)/IQ1;                 /* Start here when not initializing */
00577   *idum=IA1*(*idum-k*IQ1)-k*IR1; /* Compute idum=(IA1*idum) % IM1 without */
00578   if (*idum < 0) *idum += IM1;   /* overflows by Schrage's method */
00579   k=idum2/IQ2;
00580   idum2=IA2*(idum2-k*IQ2)-k*IR2; /* Compute idum2=(IA2*idum) % IM2 likewise */
00581   if (idum2 < 0) idum2 += IM2;
00582   j=(int)(iy/NDIV);              /* Will be in the range 0...NTAB-1 */
00583   iy=iv[j]-idum2;                /* Here idum is shuffled, idum and idum2 */
00584   iv[j] = *idum;                 /* are combined to generate output */
00585   if (iy < 1) iy += IMM1;
00586   if ((temp=AM*iy) > RNMX) return RNMX; /* No endpoint values */
00587   else return temp;
00588 }
00589 
00590 #undef IM1
00591 #undef IM2
00592 #undef AM
00593 #undef IMM1
00594 #undef IA1
00595 #undef IA2
00596 #undef IQ1
00597 #undef IQ2
00598 #undef IR1
00599 #undef IR2
00600 #undef NTAB
00601 #undef NDIV
00602 #undef RNMX