Athena: prob/cylbr.c Source File

Go to the documentation of this file.
00001 #include "copyright.h"
00002 /*============================================================================*/
00003 /*! \file cylbr.c
00004  *  \brief A dynamic test of force balance using a B_R-only, time-dependent,
00005  * non-axisymmetric magnetic field.
00006  */
00007 /*============================================================================*/
00008 
00009 #include <math.h>
00010 #include <stdio.h>
00011 #include <stdlib.h>
00012 #include <string.h>
00013 #include "defs.h"
00014 #include "athena.h"
00015 #include "globals.h"
00016 #include "prototypes.h"
00017 
00018 /*==============================================================================
00019  * PRIVATE FUNCTION PROTOTYPES:
00020  * grav_pot() - gravitational potential
00021  * grav_acc() - gravitational acceleration
00022  * d()        - density
00023  * M2()       - phi-momentum
00024  * B1()       - magnetic field
00025  * B1i()      - R-interface magnetic field
00026  * Pgas()     - gas pressure
00027  * cylbr_ix1  - inner-R boundary condition
00028  * cylbr_ox1  - outer-R boundary condition
00029  *============================================================================*/
00030 
00031 static Real br0,omega0,vz0,rho0,pgas0,a,x2min,x2max,phi0,x1save,x3save;
00032 Real d(const Real x1, const Real x2, const Real x3);
00033 Real M2(const Real x1, const Real x2, const Real x3);
00034 Real B1(const Real x1, const Real x2, const Real x3);
00035 Real B1i(const Real x2);
00036 Real Pgas(const Real x1, const Real x2, const Real x3);
00037 static Real grav_pot(const Real x1, const Real x2, const Real x3);
00038 static Real grav_acc(const Real x1, const Real x2, const Real x3);
00039 void cylbr_ix1(GridS *pG);
00040 void cylbr_ox1(GridS *pG);
00041 static ConsS ***RootSoln=NULL;
00042 
00043 /*=========================== PUBLIC FUNCTIONS ===============================*/
00044 /*----------------------------------------------------------------------------*/
00045 /* problem:  */
00046 void problem(DomainS *pDomain)
00047 {
00048   GridS *pG = pDomain->Grid;
00049   int i,j,k;
00050   int is,ie,il,iu,js,je,jl,ju,ks,ke,kl,ku;
00051   int nx1,nx2,nx3;
00052   Real x1,x2,x3,x1a,x2a,x2b;
00053   Real Eint,Ekin,Emag;
00054 
00055   is = pG->is;  ie = pG->ie;  nx1 = ie-is+1;
00056   js = pG->js;  je = pG->je;  nx2 = je-js+1;
00057   ks = pG->ks;  ke = pG->ke;  nx3 = ke-ks+1;
00058 
00059   il = is-nghost*(nx1>1);  iu = ie+nghost*(nx1>1);  nx1 = iu-il+1;
00060   jl = js-nghost*(nx2>1);  ju = je+nghost*(nx2>1);  nx2 = ju-jl+1;
00061   kl = ks-nghost*(nx3>1);  ku = ke+nghost*(nx3>1);  nx3 = ku-kl+1;
00062 
00063 #ifndef CYLINDRICAL
00064   ath_error("[cylbr]: This problem only works in cylindrical!\n");
00065 #endif
00066 
00067   if (nx1==1) {
00068     ath_error("[cylbr]: This problem can only be run in 2D or 3D!\n");
00069   }
00070   else if (nx2==1 && nx3>1) {
00071     ath_error("[cylbr]: Only (R,phi) can be used in 2D!\n");
00072   }
00073 
00074   x2min = par_getd("grid","x2min");
00075   x2max = par_getd("grid","x2max");
00076   a = 2.0*PI/(x2max-x2min);
00077   phi0 = a*(omega0*pG->time + x2min);
00078 
00079   omega0 = par_getd("problem", "omega0");
00080   vz0    = par_getd("problem", "vz0");
00081   br0    = par_getd("problem", "br0");
00082   rho0   = par_getd("problem", "rho0");
00083   pgas0  = par_getd("problem", "pgas0");
00084 
00085   /* ALLOCATE MEMORY FOR SOLUTION */
00086   if ((RootSoln = (ConsS***)calloc_3d_array(nx3,nx2,nx1,sizeof(ConsS))) == NULL)
00087     ath_error("[cylbr]: Error allocating memory for solution\n");
00088 
00089   for (k=kl; k<=ku; k++) {
00090     for (j=jl; j<=ju; j++) {
00091       for (i=il; i<=iu; i++) {
00092         x1a = x1 - 0.5*pG->dx1;
00093         x2a = x2 - 0.5*pG->dx2;
00094         x2b = x2 + 0.5*pG->dx2;
00095         x1save = x1a;
00096         x3save = x3;
00097 
00098         pG->U[k][j][i].d   = avg2d(d,pG,i,j,k);
00099         pG->U[k][j][i].M2  = avg2d(M2,pG,i,j,k);
00100         pG->U[k][j][i].M3  = pG->U[k][j][i].d*vz0;
00101         pG->U[k][j][i].B1c = avg2d(B1,pG,i,j,k);
00102         pG->B1i[k][j][i]   = qsimp(B1i,x2a,x2b)/pG->dx2;
00103 
00104         Eint = avg2d(Pgas,pG,i,j,k)/Gamma_1;
00105         Emag = 0.5*(SQR(pG->U[k][j][i].B1c) + SQR(pG->U[k][j][i].B2c) + SQR(pG->U[k][j][i].B3c));
00106         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; 
00107         pG->U[k][j][i].E = Eint + Emag + Ekin;
00108 
00109         /* SAVE SOLUTION */
00110         RootSoln[k][j][i] = pG->U[k][j][i];
00111       }
00112     }
00113   }
00114 
00115   StaticGravPot = grav_pot;
00116   x1GravAcc = grav_acc;
00117   bvals_mhd_fun(pDomain,left_x1,cylbr_ix1);
00118   bvals_mhd_fun(pDomain,right_x1,cylbr_ox1);
00119 
00120   return;
00121 }
00122 
00123 /*==============================================================================
00124  * PROBLEM USER FUNCTIONS:
00125  * problem_write_restart() - writes problem-specific user data to restart files
00126  * problem_read_restart()  - reads problem-specific user data from restart files
00127  * get_usr_expr()          - sets pointer to expression for special output data
00128  * get_usr_out_fun()       - returns a user defined output function pointer
00129  * Userwork_in_loop        - problem specific work IN     main loop
00130  * Userwork_after_loop     - problem specific work AFTER  main loop
00131  *----------------------------------------------------------------------------*/
00132 
00133 void problem_write_restart(MeshS *pM, FILE *fp)
00134 {
00135   return;
00136 }
00137 
00138 void problem_read_restart(MeshS *pM, FILE *fp)
00139 {
00140   return;
00141 }
00142 
00143 ConsFun_t get_usr_expr(const char *expr)
00144 {
00145   return NULL;
00146 }
00147 
00148 VOutFun_t get_usr_out_fun(const char *name){
00149   return NULL;
00150 }
00151 
00152 void Userwork_in_loop(MeshS *pM)
00153 {
00154   printf("Max divB = %1.10e\n", compute_div_b(pM->Domain[0][0].Grid));
00155 }
00156 
00157 void Userwork_after_loop(MeshS *pM)
00158 {
00159   compute_l1_error("CylBR", pM, RootSoln, 1);
00160 }
00161 
00162 /*=========================== PRIVATE FUNCTIONS ==============================*/
00163 
00164 /*! \fn static Real grav_pot(const Real x1, const Real x2, const Real x3)
00165  *  \brief Gravitational potential */
00166 static Real grav_pot(const Real x1, const Real x2, const Real x3) {
00167   return 0.5*SQR(x1*omega0) - 0.5*SQR(br0)/SQR(x1);
00168 }
00169 
00170 /*! \fn static Real grav_acc(const Real x1, const Real x2, const Real x3)
00171  *  \brief Gravitational acceleration. */
00172 static Real grav_acc(const Real x1, const Real x2, const Real x3) {
00173   return x1*SQR(omega0) + SQR(br0)/pow(x1,3);
00174 }
00175 
00176 /*! \fn Real d(const Real x1, const Real x2, const Real x3) 
00177  *  \brief Density */
00178 Real d(const Real x1, const Real x2, const Real x3) {
00179   return rho0 + SQR(sin(a*x2-phi0));
00180 }
00181 
00182 /*! \fn Real M2(const Real x1, const Real x2, const Real x3) 
00183  *  \brief 2-component of momentum */
00184 Real M2(const Real x1, const Real x2, const Real x3) {
00185   return (rho0 + SQR(sin(a*x2-phi0)))*omega0*x1;
00186 }
00187 
00188 /*! \fn Real B1(const Real x1, const Real x2, const Real x3)
00189  *  \brief 1-component of B field */
00190 Real B1(const Real x1, const Real x2, const Real x3) {
00191   return br0*(cos(a*x2-phi0))/x1;
00192 }
00193 
00194 /*! \fn Real B1i(const Real x2) 
00195  *  \brief 1-component of B field at interface */
00196 Real B1i(const Real x2) {
00197   return B1(x1save,x2,x3save);
00198 }
00199 
00200 /*! \fn Real Pgas(const Real x1, const Real x2, const Real x3) 
00201  *  \brief gas pressure */
00202 Real Pgas(const Real x1, const Real x2, const Real x3) {
00203   return pgas0 + 0.5*SQR(br0)*(rho0 + SQR(sin(a*x2-phi0)))/SQR(x1);
00204 }
00205 
00206 /*----------------------------------------------------------------------------*/
00207 /*! \fn void cylbr_ix1(GridS *pG)
00208  *  \brief  Inner-R boundary conditions.  d, M2, B1, B1i, and P are all
00209  *   functions of R, phi, and t.
00210  */
00211 
00212 void cylbr_ix1(GridS *pG)
00213 {
00214   int is = pG->is;
00215   int js = pG->js, je = pG->je;
00216   int ks = pG->ks, ke = pG->ke;
00217   int i,j,k;
00218   Real x1,x2,x3,x1a,x2a,x2b;
00219   Real Eint,Emag,Ekin;
00220 
00221   phi0 = a*(omega0*pG->time + x2min);
00222   for (k=ks; k<=ke; k++) {
00223     for (j=js; j<=je; j++) {
00224       for (i=1; i<=nghost; i++) {
00225         cc_pos(pG,is-i,j,k,&x1,&x2,&x3);
00226         x1a = x1 - 0.5*pG->dx1;
00227         x2a = x2 - 0.5*pG->dx2;
00228         x2b = x2 + 0.5*pG->dx2;
00229         x1save = x1a;
00230         x3save = x3;
00231 
00232         pG->U[k][j][is-i].d   = avg2d(d,pG,is-i,j,k);
00233         pG->U[k][j][is-i].M2  = avg2d(M2,pG,is-i,j,k);
00234         pG->U[k][j][is-i].M3  = pG->U[k][j][is-i].d*vz0;
00235         pG->U[k][j][is-i].B1c = avg2d(B1,pG,is-i,j,k);
00236         pG->B1i[k][j][is-i]   = qsimp(B1i,x2a,x2b)/pG->dx2;
00237 
00238         Eint = avg2d(Pgas,pG,is-i,j,k)/Gamma_1;
00239         Emag = 0.5*(SQR(pG->U[k][j][is-i].B1c) + SQR(pG->U[k][j][is-i].B2c) + SQR(pG->U[k][j][is-i].B3c));
00240         Ekin = 0.5*(SQR(pG->U[k][j][is-i].M1) + SQR(pG->U[k][j][is-i].M2) + SQR(pG->U[k][j][is-i].M3))/pG->U[k][j][is-i].d; 
00241         pG->U[k][j][is-i].E = Eint + Emag + Ekin;
00242       }
00243     }
00244   }
00245 
00246   return;
00247 }
00248 
00249 /*----------------------------------------------------------------------------*/
00250 /*! \fn void cylbr_ox1(GridS *pG)
00251  *  \brief B_R = B_0/R boundary conditions, Outer x1 boundary
00252  */
00253 void cylbr_ox1(GridS *pG)
00254 {
00255   int ie = pG->ie;
00256   int js = pG->js, je = pG->je;
00257   int ks = pG->ks, ke = pG->ke;
00258   int i,j,k;
00259   Real x1,x2,x3,x1a,x2a,x2b;
00260   Real Eint,Emag,Ekin;
00261 
00262   phi0 = a*(omega0*pG->time + x2min);
00263   for (k=ks; k<=ke; k++) {
00264     for (j=js; j<=je; j++) {
00265       for (i=1; i<=nghost; i++) {
00266         cc_pos(pG,ie+i,j,k,&x1,&x2,&x3);
00267         x1a = x1 - 0.5*pG->dx1;
00268         x2a = x2 - 0.5*pG->dx2;  
00269         x2b = x2 + 0.5*pG->dx2;
00270         x1save = x1a;  
00271         x3save = x3;
00272 
00273         pG->U[k][j][ie+i].d   = avg2d(d,pG,ie+i,j,k);
00274         pG->U[k][j][ie+i].M2  = avg2d(M2,pG,ie+i,j,k);
00275         pG->U[k][j][ie+i].M3  = pG->U[k][j][ie+i].d*vz0;
00276         pG->U[k][j][ie+i].B1c = avg2d(B1,pG,ie+i,j,k);
00277         pG->B1i[k][j][ie+i]   = qsimp(B1i,x2a,x2b)/pG->dx2;
00278 
00279         Eint = avg2d(Pgas,pG,ie+i,j,k)/Gamma_1;
00280         Emag = 0.5*(SQR(pG->U[k][j][ie+i].B1c) + SQR(pG->U[k][j][ie+i].B2c) + SQR(pG->U[k][j][ie+i].B3c));
00281         Ekin = 0.5*(SQR(pG->U[k][j][ie+i].M1) + SQR(pG->U[k][j][ie+i].M2) + SQR(pG->U[k][j][ie+i].M3))/pG->U[k][j][ie+i].d; 
00282         pG->U[k][j][ie+i].E = Eint + Emag + Ekin;
00283       }
00284     }
00285   }
00286 
00287   return;
00288 }