Athena: prob/hkdisk.c Source File

Go to the documentation of this file.
00001 #include "copyright.h"
00002 /*============================================================================*/
00003 /*! \file hkdisk.c
00004  *  \brief Problem generator for Hawley Krolik disk (Specifically, GT4)
00005  */
00006 /*============================================================================*/
00007 
00008 #include <math.h>
00009 #include <stdio.h>
00010 #include <stdlib.h>
00011 #include <string.h>
00012 #include "defs.h"
00013 #include "athena.h"
00014 #include "globals.h"
00015 #include "prototypes.h"
00016 #include "cyl.h"
00017 
00018 // Input parameters
00019 static Real q, r0, rhomax, r_in, rho0, e0, dcut, beta, seed;
00020 
00021 // Derived quantities
00022 static Real f, C, Kbar, n;
00023 
00024 /*! \fn static Real grav_pot(const Real x1, const Real x2, const Real x3) 
00025  *  \brief Gravitational potential */
00026 static Real grav_pot(const Real x1, const Real x2, const Real x3) {
00027   Real rad;
00028   rad = sqrt( SQR(x1) + SQR(x3) );
00029   return -1.0/(rad-2.0);
00030 // return 0.0;
00031 }
00032 
00033 /*! \fn static Real grav_acc(const Real x1, const Real x2, const Real x3)
00034  *  \brief Gravitational acceleration */
00035 static Real grav_acc(const Real x1, const Real x2, const Real x3) {
00036   Real rad;
00037   rad = sqrt( SQR(x1) + SQR(x3) );
00038   return (1.0/(SQR(rad-2.0)))*(x1/rad);
00039 }
00040  
00041 //Private functions (x1,x2,x3) = (R,p,z)
00042 
00043 /*! \fn Real density(Real x1, Real x2, Real x3) 
00044  *  \brief Calculates the density at x1, x2, x3*/
00045 Real density(Real x1, Real x2, Real x3) {
00046   Real rad, temp, d;
00047   rad = sqrt( SQR(x1) + SQR(x3));
00048   temp = C + (1.0/(rad-2.0)) - f*pow(x1,-2.0*(q-1.0));
00049   temp = MAX(temp,0.0);
00050   d = pow(temp/Kbar,n)*(x1>=r_in);
00051   d = MAX(d,rho0);
00052 //   d = x1 < r_in ? rho0 : d;
00053   return d;
00054 }
00055 
00056 /*! \fn Real Volume(Grid *pG, int i, int j, int k) 
00057  *  \brief Calculates the volume of cell (i,j,k) */
00058 Real Volume(Grid *pG, int i, int j, int k) {
00059   Real x1,x2,x3;
00060   cc_pos(pG,i,j,k,&x1,&x2,&x3);
00061   // Volume_ijk = R_i * dR * dPhi * dZ
00062   return x1*pG->dx1*pG->dx2*pG->dx3;
00063 }
00064 
00065 /*----------------------------------------------------------------------------*/
00066 /* problem:   */
00067 #define VP(R) ((sqrt(r0)/(r0-2))*pow(R/r0,-q+1))
00068 void problem(Grid *pG, Domain *pDomain)
00069 {
00070   int i,j,k;
00071   int is,ie,js,je,ks,ke,nx1,nx2,nx3;
00072   int il,iu,jl,ju,kl,ku;
00073   Real rad, IntE, KinE, MagE;
00074   Real x1,x2,x3, x1i, x2i, x3i, T, Ta, rhoa,q1,r02;
00075   Real Pgas, Pb, TotPgas, TotPb;
00076   Real scl, ***Ap, ***Bri, ***Bzi;
00077 Real divB=0.0, maxdivB=0.0;
00078 
00079   is = pG->is; ie = pG->ie;
00080   js = pG->js; je = pG->je;
00081   ks = pG->ks; ke = pG->ke;
00082   nx1 = (ie-is)+1 + 2*nghost;
00083   nx2 = (je-js)+1 + 2*nghost;
00084   nx3 = (ke-ks)+1 + 2*nghost;
00085   il = is - nghost*(ie > is);
00086   jl = js - nghost*(je > js);
00087   kl = ks - nghost*(ke > ks);
00088   iu = ie + nghost*(ie > is);
00089   ju = je + nghost*(je > js);
00090   ku = ke + nghost*(ke > ks);
00091 
00092 
00093   if ((Ap = (Real***)calloc_3d_array(nx3+1, nx2+1, nx1+1, sizeof(Real))) == NULL) {
00094     ath_error("[HK-Disk]: Error allocating memory for vector pot\n");
00095   }
00096 
00097   if ((Bri = (Real***)calloc_3d_array(nx3+1, nx2+1, nx1+1, sizeof(Real))) == NULL) {
00098     ath_error("[HK-Disk]: Error allocating memory for Bri\n");
00099   }
00100 
00101   if ((Bzi = (Real***)calloc_3d_array(nx3+1, nx2+1, nx1+1, sizeof(Real))) == NULL) {
00102     ath_error("[HK-Disk]: Error allocating memory for Bzi\n");
00103   }
00104 
00105   /* Read initial conditions */
00106   q      = par_getd("problem","q");
00107   r0     = par_getd("problem","r0");
00108   rhomax = par_getd("problem","rhomax");
00109   r_in   = par_getd("problem","r_in");
00110   rho0   = par_getd("problem","rho0");
00111   e0     = par_getd("problem","e0");
00112         seed   = par_getd("problem","seed");
00113 
00114 #ifdef MHD
00115   dcut = par_getd("problem","dcut");
00116   beta = par_getd("problem","beta");
00117 #endif
00118 
00119 
00120   // Set up physical parameters
00121   n = 1.0/Gamma_1;
00122   q1 = 2.0*(q-1.0);
00123   r02 = (r0-2.0);
00124 
00125   f = pow(r0, 2.0*q-1.0)/(q1*pow(r02,2.0));
00126   C = f*pow(r_in,-1.0*q1) - 1.0/(r_in-2.0);
00127   Kbar = (C + (1.0/r02) - f*pow(r0,-1.0*q1))/pow(rhomax, 1.0/n);
00128 
00129   // Initialize data structures, set density, energy, velocity and Ap
00130   // Loop over all cells (including ghosts)
00131   for (k=kl; k<=ku; k++) {
00132     for (j=jl; j<=ju; j++) {
00133       for (i=il; i<=iu; i++) {
00134         cc_pos(pG,i,j,k,&x1,&x2,&x3);
00135         rad = sqrt(SQR(x1) + SQR(x3));
00136 
00137         x1i = x1 - 0.5*pG->dx1;
00138         x2i = x2 - 0.5*pG->dx2;
00139         x3i = x3 - 0.5*pG->dx3;
00140 
00141         pG->U[k][j][i].d = rho0;
00142         pG->U[k][j][i].M1 = 0.0;
00143         pG->U[k][j][i].M2 = VP(x1)*rho0;
00144         pG->U[k][j][i].M3 = 0.0;
00145         pG->U[k][j][i].E   = e0;
00146         Ap[k][j][i] = 0.0;
00147 
00148 #ifdef MHD
00149         pG->B1i[k][j][i] = 0.0;
00150         pG->B2i[k][j][i] = 0.0;
00151         pG->B3i[k][j][i] = 0.0;
00152         pG->U[k][j][i].B1c = 0.0;
00153         pG->U[k][j][i].B2c = 0.0;
00154         pG->U[k][j][i].B3c = 0.0;
00155 
00156 #endif
00157 
00158         // Set up torus
00159         Ta = f*pow(x1, -1.0*q1) - C;
00160         T = (1.0/Ta) + 2.0;
00161 
00162         if ( (x1 >= r_in) && (rad <= T) ) { //Checks to see if cell is in torus
00163           rhoa = density(x1, x2, x3);
00164           IntE = pow(rhoa,Gamma)*Kbar/((n+1.0)*Gamma_1);
00165                                         // Add pressure fluctuations to seed instability
00166                                         IntE = IntE*(1.0 - seed*sin(x2));
00167           if ((IntE >= e0) && (rhoa >= rho0)) {
00168                                                 //If the values are above cutoff, set up the cell
00169             pG->U[k][j][i].d = rhoa;
00170             pG->U[k][j][i].M2 = VP(x1)*pG->U[k][j][i].d;
00171             pG->U[k][j][i].E = IntE;
00172                                                 
00173             //Note, at this point, E holds only the internal energy.  This must
00174             //be fixed later
00175           }
00176         }
00177       }
00178     }
00179   }
00180 
00181 
00182   // Calculate density at corner and set up Ap if appropriate
00183   for (k=kl; k<=ku+1; k++) {
00184     for (j=jl; j<=ju+1; j++) {
00185       for (i=il; i<=iu+1; i++) {
00186         cc_pos(pG,i,j,k,&x1,&x2,&x3);
00187         rad = sqrt(SQR(x1) + SQR(x3));
00188         x1i = x1 - 0.5*pG->dx1;
00189         x2i = x2 - 0.5*pG->dx2;
00190         x3i = x3 - 0.5*pG->dx3;
00191 
00192         rhoa = density(x1i,x2i,x3i);
00193 
00194         if (rhoa >= dcut) {
00195           // Ap = (density-cutoff)^2
00196           Ap[k][j][i] = SQR(rhoa-dcut);
00197         }
00198       }
00199     }
00200   }
00201 
00202 
00203 
00204 #ifdef MHD
00205   // Set up interface magnetic fields by using Ap
00206   for (k=kl; k<=ku; k++) {
00207     for (j=jl; j<=ju; j++) {
00208       for (i=il; i<=iu; i++) {
00209         cc_pos(pG,i,j,k,&x1,&x2,&x3);
00210         x1i = x1 - 0.5*pG->dx1;
00211         x2i = x2 - 0.5*pG->dx2;
00212         x3i = x3 - 0.5*pG->dx3;
00213         // Br = -dAp/dz
00214         pG->B1i[k][j][i] = -(Ap[k+1][j][i] - Ap[k][j][i])/pG->dx3;
00215         // Bz = (1/R)*d(R*Ap)/dr
00216         pG->B3i[k][j][i] = (Ap[k][j][i+1]*(x1i+pG->dx1) - Ap[k][j][i]*x1i)/(x1*pG->dx1);
00217       }
00218     }
00219   }
00220 
00221   // Calculate cell centered fields
00222   for (k=kl; k<=ku; k++) {
00223     for (j=jl; j<=ju; j++) {
00224       for (i=il; i<=iu; i++) {
00225         cc_pos(pG,i,j,k,&x1,&x2,&x3);
00226         if (i==iu)
00227           pG->U[k][j][i].B1c = pG->B1i[k][j][i];
00228         else
00229           pG->U[k][j][i].B1c = 0.5*((x1-0.5*pG->dx1)*pG->B1i[k][j][i] + (x1+0.5*pG->dx1)*pG->B1i[k][j][i+1])/x1;
00230         if (k==ku)
00231           pG->U[k][j][i].B3c = pG->B3i[k][j][i];
00232         else
00233           pG->U[k][j][i].B3c = 0.5*(pG->B3i[k+1][j][i] + pG->B3i[k][j][i]);
00234       }
00235     }
00236   }
00237 #endif //MHD
00238 #ifdef MHD
00239   // Calculate scaling factor to satisfy beta, specifically Pgas and Pb per tile
00240   // Don't loop over ghosts
00241   Pgas = 0.0;
00242   Pb   = 0.0;
00243 
00244   for (k=ks; k<=ke; k++) {
00245     for (j=js; j<=je; j++) {
00246       for (i=is; i<=ie; i++) {
00247         Pgas += (Gamma-1)*pG->U[k][j][i].E*Volume(pG,i,j,k);
00248         Pb += 0.5*(SQR(pG->U[k][j][i].B1c) + SQR(pG->U[k][j][i].B2c)
00249               + SQR(pG->U[k][j][i].B3c))*Volume(pG,i,j,k);
00250       }
00251     }
00252   }
00253 #endif //MHD
00254 #ifdef MPI_PARALLEL
00255   MPI_Reduce(&Pgas, &TotPgas, 1, MPI_DOUBLE, MPI_SUM, 0, MPI_COMM_WORLD);
00256   MPI_Reduce(&Pb, &TotPb, 1, MPI_DOUBLE, MPI_SUM, 0, MPI_COMM_WORLD);
00257   if (pG->my_id == 0) {
00258     printf("Total gas pressure = %f\n", TotPgas);
00259     printf("Total magnetic pressure = %f\n", TotPb);
00260   }
00261   MPI_Bcast(&TotPgas, 1, MPI_DOUBLE, 0, MPI_COMM_WORLD);
00262   MPI_Bcast(&TotPb, 1, MPI_DOUBLE, 0, MPI_COMM_WORLD);
00263 #else
00264   TotPgas = Pgas;
00265   TotPb = Pb;
00266 #endif //PARALLEL
00267 
00268 #ifdef MHD
00269   //calculate and use scaling factor so that the correct beta is ensured
00270   scl = sqrt(TotPgas/(TotPb*beta));
00271   printf("Using magnetic scaling factor %f\n", scl);
00272   for (k=kl; k<=ku; k++) {
00273     for (j=jl; j<=ju; j++) {
00274       for (i=il; i<=iu; i++) {
00275         pG->U[k][j][i].B1c *= scl;
00276         pG->U[k][j][i].B3c *= scl;
00277         pG->B1i[k][j][i]   *= scl;
00278         pG->B3i[k][j][i]   *= scl;
00279       }
00280     }
00281   }
00282 #endif
00283 
00284   // Fix E to hold the total energy (Kin + Int + Mag) (Loop over ghosts)
00285   for (k=kl; k<=ku; k++) {
00286     for (j=jl; j<=ju; j++) {
00287       for (i=il; i<=iu; i++) {
00288         KinE = 0.5*(SQR(pG->U[k][j][i].M1) + SQR(pG->U[k][j][i].M2) 
00289                 + SQR(pG->U[k][j][i].M3))/(pG->U[k][j][i].d);
00290         MagE = 0.0;
00291 #ifdef MHD
00292         MagE = 0.5*(SQR(pG->U[k][j][i].B1c) + SQR(pG->U[k][j][i].B2c)
00293                 + SQR(pG->U[k][j][i].B3c));
00294 #endif
00295         pG->U[k][j][i].E += KinE + MagE;
00296       }
00297     }
00298   }
00299   /* Enroll the gravitational function and radial BC */
00300   StaticGravPot = grav_pot;
00301   x1GravAcc = grav_acc;
00302   set_bvals_fun(left_x1,disk_ir_bc);
00303   set_bvals_fun(right_x1,disk_or_bc);
00304   set_bvals_fun(left_x3,diode_outflow_ix3);
00305   set_bvals_fun(right_x3,diode_outflow_ox3);
00306 
00307   return;
00308 }
00309 
00310 /*==============================================================================
00311  * PROBLEM USER FUNCTIONS:
00312  * problem_write_restart() - writes problem-specific user data to restart files
00313  * problem_read_restart()  - reads problem-specific user data from restart files
00314  * get_usr_expr()          - sets pointer to expression for special output data
00315  * Userwork_in_loop        - problem specific work IN     main loop
00316  * Userwork_after_loop     - problem specific work AFTER  main loop
00317  * current() - computes x3-component of current
00318  * Bp2()     - computes magnetic pressure (Bx2 + By2)
00319  *----------------------------------------------------------------------------*/
00320 
00321 void problem_write_restart(Grid *pG, Domain *pD, FILE *fp)
00322 {
00323   return;
00324 }
00325 
00326 void problem_read_restart(Grid *pG, Domain *pD, FILE *fp)
00327 {
00328   q = par_getd("problem","q");
00329   r0 = par_getd("problem","r0");
00330   rhomax = par_getd("problem","rhomax");
00331   r_in = par_getd("problem","r_in");
00332   rho0 = par_getd("problem","rho0");
00333   e0 = par_getd("problem","e0");
00334         seed = par_getd("problem","seed");
00335 
00336 #ifdef MHD
00337   dcut = par_getd("problem","dcut");
00338   beta = par_getd("problem","beta");
00339 #endif
00340 
00341   /* Enroll the gravitational function and radial BC */
00342   StaticGravPot = grav_pot;
00343   x1GravAcc = grav_acc;
00344   set_bvals_fun(left_x1,disk_ir_bc);
00345   set_bvals_fun(right_x1,disk_or_bc);
00346   set_bvals_fun(left_x3,diode_outflow_ix3);
00347   set_bvals_fun(right_x3,diode_outflow_ox3);
00348 
00349   return;
00350 }
00351 
00352 
00353 Gasfun_t get_usr_expr(const char *expr)
00354 {
00355   return NULL;
00356 }
00357 
00358 void Userwork_in_loop(Grid *pG, Domain *pDomain)
00359 {
00360   int i,j,k,is,ie,js,je,ks,ke, nx1, nx2, nx3, il, iu, jl,ju,kl,ku;
00361         int prote, protd;
00362   Real IntE, KinE, MagE=0.0, x1, x2, x3, DivB;
00363   static Real TotMass=0.0;
00364 
00365   is = pG->is; ie = pG->ie;
00366   js = pG->js; je = pG->je;
00367   ks = pG->ks; ke = pG->ke;
00368   nx1 = (ie-is)+1 + 2*nghost;
00369   nx2 = (je-js)+1 + 2*nghost;
00370   nx3 = (ke-ks)+1 + 2*nghost;
00371   il = is - nghost*(ie > is);
00372   jl = js - nghost*(je > js);
00373   kl = ks - nghost*(ke > ks);
00374   iu = ie + nghost*(ie > is);
00375   ju = je + nghost*(je > js);
00376   ku = ke + nghost*(ke > ks);
00377 
00378         // Verify divB
00379   protd = 0;
00380   prote = 0;
00381 #ifdef MHD
00382   DivB = compute_div_b(pG);
00383 #endif 
00384 
00385 
00386   for (k=kl; k<=ku; k++) {
00387     for (j=jl; j<=ju; j++) {
00388       for (i=il; i<=iu; i++) {
00389         cc_pos(pG,i,j,k,&x1,&x2,&x3);
00390         if (isnan(pG->U[k][j][i].d) || isnan(pG->U[k][j][i].E)) {
00391           printf("At pos (%f,%f,%f), Den = %f, E = %f\n", x1,x2,x3,pG->U[k][j][i].d, pG->U[k][j][i].E);
00392         }
00393 
00394         KinE = 0.5*(SQR(pG->U[k][j][i].M1) + SQR(pG->U[k][j][i].M2) 
00395                 + SQR(pG->U[k][j][i].M3))/(pG->U[k][j][i].d);
00396 #ifdef MHD
00397         MagE = 0.5*(SQR(pG->U[k][j][i].B1c) + SQR(pG->U[k][j][i].B2c)
00398                 + SQR(pG->U[k][j][i].B3c));
00399 #endif
00400         IntE = pG->U[k][j][i].E - KinE - MagE;
00401 
00402         if (pG->U[k][j][i].d < rho0) {
00403 
00404           pG->U[k][j][i].d = rho0;
00405 
00406           KinE = 0.5*(SQR(pG->U[k][j][i].M1) + SQR(pG->U[k][j][i].M2)
00407                   + SQR(pG->U[k][j][i].M3))/(pG->U[k][j][i].d);
00408 
00409           IntE = e0;  // Set this as well to keep c_s reasonable
00410 
00411           pG->U[k][j][i].E = IntE + KinE + MagE;
00412                                         protd++;
00413                                         prote++;
00414         }
00415 
00416 
00417 
00418         else if (IntE < e0) {
00419           pG->U[k][j][i].E = e0 + KinE + MagE;
00420                                         prote++;
00421         }
00422         IntE = pG->U[k][j][i].E - KinE - MagE;
00423 
00424       }
00425     }
00426   }
00427 
00428 #ifdef MPI_PARALLEL
00429         if (pG->my_id == 0) {
00430         printf("\tDivergence @ Orbit %2.3f = %e\n",(pG->time)/61.6, DivB);
00431         if ((protd+prote) > 0) {
00432         printf("\tProtection enforced (D=%d,E=%d), Cumulative Mass = %2.5f\n", protd, prote,TotMass);
00433       }
00434   }
00435 #else
00436   printf("\tDivergence @ Orbit %2.3f = %e\n",(pG->time)/61.6, DivB);
00437   if ((protd+prote) > 0) {
00438         printf("\tProtection enforced (D=%d,E=%d), Cumulative Mass = %2.5f\n", protd, prote,TotMass);
00439   }
00440 #endif //PARALLEL
00441         
00442 
00443 }
00444 
00445 void Userwork_after_loop(Grid *pG, Domain *pDomain)
00446 {
00447 }