Athena: prob/cylblast.c Source File

Go to the documentation of this file.
00001 #include "copyright.h"
00002 /*============================================================================*/
00003 /*! \file cylblast.c
00004  *  \brief Problem generator for blast wave in cylindrical coords.
00005  *
00006  * PURPOSE: Problem generator for blast wave in cylindrical coords.  Can only
00007  *   be run in 2D or 3D.  Input parameters are:
00008  *   -  problem/radius = radius of field initial overpressured region
00009  *   -  problem/pamb   = ambient pressure
00010  *   -  problem/prat   = ratio of interior to ambient pressure
00011  *   -  problem/b0     = initial azimuthal magnetic field (units sqrt(Pamb))
00012  *   -  problem/rho0   = background density
00013  *   -  problem/omega0 = initial azimuthal flow angular velocity
00014  *
00015  * REFERENCE: P. Londrillo & L. Del Zanna, "High-order upwind schemes for
00016  *   multidimensional MHD", ApJ, 530, 508 (2000), and references therein. */
00017 /*============================================================================*/
00018 
00019 #include <math.h>
00020 #include <stdio.h>
00021 #include <stdlib.h>
00022 #include <string.h>
00023 #include "defs.h"
00024 #include "athena.h"
00025 #include "globals.h"
00026 #include "prototypes.h"
00027 
00028 /*==============================================================================
00029  * PRIVATE FUNCTION PROTOTYPES:
00030  * grav_pot() - gravitational potential
00031  * grav_acc() - gravitational acceleration
00032  * M2()       - phi-momentum
00033  *============================================================================*/
00034 
00035 static Real omega0,rho0;
00036 static Real grav_pot(const Real x1, const Real x2, const Real x3);
00037 static Real grav_acc(const Real x1, const Real x2, const Real x3);
00038 Real M2(const Real x1, const Real x2, const Real x3);
00039 
00040 /*=========================== PUBLIC FUNCTIONS ===============================*/
00041 /*----------------------------------------------------------------------------*/
00042 /* problem:   */
00043 
00044 void problem(DomainS *pDomain)
00045 {
00046   GridS *pG = pDomain->Grid;
00047   int i,j,k;
00048   int is,ie,js,je,ks,ke,nx1,nx2,nx3;
00049   int il,iu,jl,ju,kl,ku;
00050   Real r0,phi0,x0,y0,z0,angle,radius,prat,pamb,b0;
00051   Real x1,x2,x3,x2i;
00052   Real x,y,z,Eint,Emag,Ekin;
00053 
00054   is = pG->is;  ie = pG->ie;  nx1 = ie-is+1;
00055   js = pG->js;  je = pG->je;  nx2 = je-js+1;
00056   ks = pG->ks;  ke = pG->ke;  nx3 = ke-ks+1;
00057 
00058   il = is-nghost*(nx1>1);  iu = ie+nghost*(nx1>1);  nx1 = iu-il+1;
00059   jl = js-nghost*(nx2>1);  ju = je+nghost*(nx2>1);  nx2 = ju-jl+1;
00060   kl = ks-nghost*(nx3>1);  ku = ke+nghost*(nx3>1);  nx3 = ku-kl+1;
00061 
00062 #ifndef CYLINDRICAL
00063   ath_error("[cylblast]: This problem only works in cylindrical!\n");
00064 #endif
00065 
00066   if (nx1==1) {
00067     ath_error("[cylblast]: This problem can only be run in 2D or 3D!\n");
00068   }
00069   else if (nx2==1 && nx3>1) {
00070     ath_error("[cylblast]: Only (R,phi) can be used in 2D!\n");
00071   }
00072 
00073   /* READ IN INITIAL CONDITIONS */
00074   radius = par_getd("problem","radius");
00075   pamb   = par_getd("problem","pamb");
00076   prat   = par_getd("problem","prat");
00077   rho0   = par_getd("problem","rho0");
00078   omega0 = par_getd("problem","omega0");
00079   b0     = par_getd("problem","b0");
00080 
00081   /* PLACEMENT OF CENTER OF BLAST */
00082   r0   = par_getd("problem","r0");
00083   phi0 = par_getd("problem","phi0");
00084   z0   = par_getd("problem","z0");
00085 
00086   /* ORIENTATION OF FIELD W.R.T. POS. X-AXIS */
00087   angle = (PI/180.0)*par_getd("problem","angle");
00088 
00089   x0 = r0*cos(phi0);
00090   y0 = r0*sin(phi0);
00091 
00092   /* SET UP UNIFORM AMBIENT MEDIUM WITH CIRCULAR OVER-PRESSURED REGION */
00093   for (k=kl; k<=ku; k++) {
00094     for (j=jl; j<=ju; j++) {
00095       for (i=il; i<=iu; i++) {
00096         cc_pos(pG,i,j,k,&x1,&x2,&x3);
00097         x2i = x2 - 0.5*pG->dx2;
00098 
00099         pG->U[k][j][i].d  = rho0;
00100         pG->U[k][j][i].M1 = 0.0;
00101         pG->U[k][j][i].M2 = pG->U[k][j][i].d*x1*omega0;
00102         pG->U[k][j][i].M3 = 0.0;
00103 #ifdef MHD
00104         /* SET UP A PLANAR MAGNETIC FIELD IN THE X-Y (R-PHI) PLANE */
00105         pG->B1i[k][j][i]   = b0*(cos(angle)*cos(x2)+sin(angle)*sin(x2)); 
00106         pG->B2i[k][j][i]   = b0*(-cos(angle)*sin(x2i)+sin(angle)*cos(x2i));
00107         pG->B3i[k][j][i]   = 0.0;
00108         pG->U[k][j][i].B1c = b0*(cos(angle)*cos(x2)+sin(angle)*sin(x2));
00109         pG->U[k][j][i].B2c = b0*(-cos(angle)*sin(x2)+sin(angle)*cos(x2));
00110         pG->U[k][j][i].B3c = 0.0;
00111 #endif
00112         /* CARTESIAN POSITION OF CELL CENTER */
00113         x = x1*cos(x2);
00114         y = x1*sin(x2);
00115         z = x3;
00116 
00117         /* INITIALIZE TOTAL ENERGY */
00118 #ifndef ISOTHERMAL
00119         Eint = pamb/Gamma_1;
00120         if (SQR(x-x0) + SQR(y-y0) + SQR(z-z0) < SQR(radius)) {
00121           Eint *= prat;
00122         }
00123         Emag = 0.0;
00124 #ifdef MHD
00125         Emag = 0.5*(SQR(pG->U[k][j][i].B1c) + SQR(pG->U[k][j][i].B2c) + SQR(pG->U[k][j][i].B3c));
00126 #endif
00127         Ekin = 0.5*(SQR(pG->U[k][j][i].M1) + SQR(pG->U[k][j][i].M2) + SQR(pG->U[k][j][i].M3))/pG->U[k][j][i].d;
00128         pG->U[k][j][i].E = Eint + Emag + Ekin;
00129 #endif /* ISOTHERMAL */
00130       }
00131     }
00132   }
00133 
00134   /* Enroll the gravitational function and radial BC */
00135   StaticGravPot = grav_pot;
00136   x1GravAcc = grav_acc;
00137   bvals_mhd_fun(pDomain,left_x1, do_nothing_bc);
00138   bvals_mhd_fun(pDomain,right_x1,do_nothing_bc);
00139   bvals_mhd_fun(pDomain,left_x2, do_nothing_bc);
00140   bvals_mhd_fun(pDomain,right_x2,do_nothing_bc);
00141 
00142   return;
00143 }
00144 
00145 /*==============================================================================
00146  * PROBLEM USER FUNCTIONS:
00147  * problem_write_restart() - writes problem-specific user data to restart files
00148  * problem_read_restart()  - reads problem-specific user data from restart files
00149  * get_usr_expr()          - sets pointer to expression for special output data
00150  * get_usr_out_fun()       - returns a user defined output function pointer
00151  * Userwork_in_loop        - problem specific work IN     main loop
00152  * Userwork_after_loop     - problem specific work AFTER  main loop
00153  *----------------------------------------------------------------------------*/
00154 
00155 void problem_write_restart(MeshS *pM, FILE *fp)
00156 {
00157   return;
00158 }
00159 
00160 void problem_read_restart(MeshS *pM, FILE *fp)
00161 {
00162   return;
00163 }
00164 
00165 ConsFun_t get_usr_expr(const char *expr)
00166 {
00167   return NULL;
00168 }
00169 
00170 VOutFun_t get_usr_out_fun(const char *name){
00171   return NULL;
00172 }
00173 
00174 void Userwork_in_loop(MeshS *pM)
00175 {
00176 }
00177 
00178 void Userwork_after_loop(MeshS *pM)
00179 {
00180 }
00181 
00182 /*=========================== PRIVATE FUNCTIONS ==============================*/
00183 
00184 /*! \fn static Real grav_pot(const Real x1, const Real x2, const Real x3) {
00185  *  \brief  Gravitational potential*/
00186 static Real grav_pot(const Real x1, const Real x2, const Real x3) {
00187   return 0.5*SQR(x1*omega0);
00188 }
00189 
00190 /*! \fn static Real grav_acc(const Real x1, const Real x2, const Real x3) {
00191  *  \brief Gravitational acceleration  */
00192 static Real grav_acc(const Real x1, const Real x2, const Real x3) {
00193   return x1*SQR(omega0);
00194 }
00195 
00196 /*! \fn Real M2(const Real x1, const Real x2, const Real x3) 
00197  *  \brief 2-component of momentum */
00198 Real M2(const Real x1, const Real x2, const Real x3) {
00199   return rho0*omega0*x1;
00200 }