Athena: prob/cyladvect.c Source File

Go to the documentation of this file.
00001 #include "copyright.h"
00002 /*============================================================================*/
00003 /*! \file cyladvect.c
00004  *  \brief A simple density-pulse advection test in cylindrical coordinates with
00005  * no pressure or tension forces.
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  *============================================================================*/
00023 
00024 static Real omega0,rho0,bz0,vz0,Pgas0,amp,R0,phi0,z0,rad,alpha,x2min,x2max;
00025 static int iprob;
00026 static Real grav_pot(const Real x1, const Real x2, const Real x3);
00027 static Real grav_acc(const Real x1, const Real x2, const Real x3);
00028 Real d(const Real x1, const Real x2, const Real x3);
00029 Real M2(const Real x1, const Real x2, const Real x3);
00030 void cyladvect_ix1(GridS *pG);
00031 void cyladvect_ox1(GridS *pG);
00032 static ConsS ***RootSoln=NULL;
00033 
00034 /*=========================== PUBLIC FUNCTIONS ===============================*/
00035 /*----------------------------------------------------------------------------*/
00036 /* problem:  */
00037 void problem(DomainS *pDomain)
00038 {
00039   GridS *pG = pDomain->Grid;
00040   int i,j,k;
00041   int is,ie,il,iu,js,je,jl,ju,ks,ke,kl,ku;
00042   int nx1,nx2,nx3; 
00043   Real x1,x2,x3;
00044   Real Eint,Emag,Ekin;
00045   int mask=0;
00046 
00047   is = pG->is;  ie = pG->ie;  nx1 = ie-is+1;
00048   js = pG->js;  je = pG->je;  nx2 = je-js+1;
00049   ks = pG->ks;  ke = pG->ke;  nx3 = ke-ks+1;
00050 
00051   il = is-nghost*(nx1>1);  iu = ie+nghost*(nx1>1);  nx1 = iu-il+1;
00052   jl = js-nghost*(nx2>1);  ju = je+nghost*(nx2>1);  nx2 = ju-jl+1;
00053   kl = ks-nghost*(nx3>1);  ku = ke+nghost*(nx3>1);  nx3 = ku-kl+1;
00054 
00055 #ifndef CYLINDRICAL
00056   ath_error("[cyladvect]: This problem only works in cylindrical!\n");
00057 #endif
00058 
00059   if (nx1==1) {
00060     ath_error("[cyladvect]: This problem can only be run in 2D or 3D!\n");
00061   }
00062   else if (nx2==1 && nx3>1) {
00063     ath_error("[cyladvect]: Only (R,phi) can be used in 2D!\n");
00064   }
00065 
00066   /* ALLOCATE MEMORY FOR SOLUTION */
00067   if ((RootSoln = (ConsS***)calloc_3d_array(nx3,nx2,nx1,sizeof(ConsS))) == NULL)
00068     ath_error("[cyladvect]: Error allocating memory for solution\n");
00069 
00070   /* PARSE INPUT FILE */
00071   x2min  = par_getd("domain1", "x2min");
00072   x2max  = par_getd("domain1", "x2max");
00073   omega0 = par_getd("problem", "omega0");
00074   rho0   = par_getd("problem", "rho0");
00075   bz0    = par_getd("problem", "bz0");
00076   vz0    = par_getd("problem", "vz0");
00077   Pgas0  = par_getd("problem", "Pgas0");
00078   amp    = par_getd("problem", "amp");
00079   R0     = par_getd("problem", "R0");
00080   phi0   = par_getd("problem", "phi0");
00081   z0     = par_getd("problem", "z0");
00082   rad    = par_getd("problem", "rad");
00083   iprob  = par_geti("problem", "iprob");
00084 
00085   /* PERIOD IS ONE GRID-CROSSING */
00086   alpha = 2.0*PI/(x2max-x2min);
00087 
00088   /* SET DENSITY AND PHI-MOMENTUM (MUST USE VOLUME CENTER) */
00089   for (k=kl; k<=ku; k++) {
00090     for (j=jl; j<=ju; j++) {
00091       for (i=il; i<=iu; i++) {
00092         cc_pos(pG,i,j,k,&x1,&x2,&x3);
00093         memset(&(pG->U[k][j][i]),0.0,sizeof(ConsS));
00094 
00095         switch (iprob) {
00096           case 1:
00097                   mask = (int)((x1>=R0-rad) && (x1<=R0+rad) && (x2>=phi0-rad/R0) && (x2<=phi0+rad/R0));
00098                   pG->U[k][j][i].d  = rho0*(1.0 + amp*mask);
00099                   pG->U[k][j][i].M2 = pG->U[k][j][i].d*x1vc(pG,i)*omega0;
00100 #ifdef MHD
00101                   pG->U[k][j][i].B3c = bz0*mask;
00102                   pG->B3i[k][j][i]   = bz0*mask;
00103 #endif
00104                   break;
00105           case 2:
00106           case 3:
00107                   pG->U[k][j][i].d  = avg2d(d, pG,i,j,k);
00108                   pG->U[k][j][i].M2 = avg2d(M2,pG,i,j,k);
00109 #ifdef MHD
00110                   pG->U[k][j][i].B3c = bz0;
00111                   pG->B3i[k][j][i]   = bz0;
00112 #endif
00113                   break;
00114           default:
00115                   ath_error("[cyladvect]:  Not an accepted problem number\n");
00116         }
00117 
00118         pG->U[k][j][i].M3 = pG->U[k][j][i].d*vz0;
00119 
00120 #ifndef ISOTHERMAL
00121         /* INITIALIZE TOTAL ENERGY */
00122         Emag = 0.0;
00123 #ifdef MHD
00124         Emag = 0.5*SQR(pG->U[k][j][i].B3c);
00125 #endif
00126         Eint = (Pgas0-Emag)/Gamma_1;
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 
00129         pG->U[k][j][i].E = Eint + Emag + Ekin;
00130 #endif /* ISOTHERMAL */
00131 
00132         /* SAVE SOLUTION */
00133         RootSoln[k][j][i] = pG->U[k][j][i];
00134       }
00135     }
00136   }
00137 
00138   StaticGravPot = grav_pot;
00139   x1GravAcc = grav_acc;
00140   if (iprob==2) {
00141     bvals_mhd_fun(pDomain, left_x1,  cyladvect_ix1);
00142     bvals_mhd_fun(pDomain, right_x1, cyladvect_ox1);
00143   } else {
00144     bvals_mhd_fun(pDomain, left_x1,  do_nothing_bc);
00145     bvals_mhd_fun(pDomain, right_x1, do_nothing_bc);
00146   }
00147 
00148   return;
00149 }
00150 
00151 
00152 /*==============================================================================
00153  * PROBLEM USER FUNCTIONS:
00154  * problem_write_restart() - writes problem-specific user data to restart files
00155  * problem_read_restart()  - reads problem-specific user data from restart files
00156  * get_usr_expr()          - sets pointer to expression for special output data
00157  * get_usr_out_fun()       - returns a user defined output function pointer
00158  * Userwork_in_loop        - problem specific work IN     main loop
00159  * Userwork_after_loop     - problem specific work AFTER  main loop
00160  *----------------------------------------------------------------------------*/
00161 
00162 void problem_write_restart(MeshS *pM, FILE *fp)
00163 {
00164   return;
00165 }
00166 
00167 void problem_read_restart(MeshS *pM, FILE *fp)
00168 {
00169   return;
00170 }
00171 
00172 ConsFun_t get_usr_expr(const char *expr)
00173 {
00174   return NULL;
00175 }
00176 
00177 VOutFun_t get_usr_out_fun(const char *name){
00178   return NULL;
00179 }
00180 
00181 void Userwork_in_loop(MeshS *pM)
00182 {
00183 #ifdef MHD
00184   printf("Max divB = %1.10e\n", compute_div_b(pM->Domain[0][0].Grid));
00185 #endif
00186 }
00187 
00188 void Userwork_after_loop(MeshS *pM)
00189 {
00190   compute_l1_error("CylAdvect", pM, RootSoln, 0);
00191 }
00192 
00193 /*=========================== PRIVATE FUNCTIONS ==============================*/
00194 
00195 /*! \fn Real d(const Real x1, const Real x2, const Real x3) {
00196  *  \brief Density */
00197 Real d(const Real x1, const Real x2, const Real x3) {
00198   Real x,y,z,x0,y0,r;
00199 
00200   if (iprob==2) {  /* SINE-WAVE, HYDRO ONLY */
00201     return rho0*(1.0 + amp*sin(alpha*x2-phi0));
00202   }
00203   else if (iprob==3) {  /* GAUSSIAN PULSE, HYDRO ONLY */
00204     x  = x1*cos(x2);    y  = x1*sin(x2);    z  = x3;
00205     x0 = R0*cos(phi0);  y0 = R0*sin(phi0);
00206     r = sqrt(SQR(x-x0) + SQR(y-y0) + SQR(z-z0));
00207     return rho0*(1.0 + exp(-0.5*SQR(3.0*r/rad)));
00208   }
00209 
00210   return rho0;
00211 }
00212 
00213 /*! \fn Real M2(const Real x1, const Real x2, const Real x3) 
00214  *  \brief 2-component of momentum  */
00215 Real M2(const Real x1, const Real x2, const Real x3) {
00216   return d(x1,x2,x3)*omega0*x1;
00217 }
00218 
00219 /*! \fn static Real grav_pot(const Real x1, const Real x2, const Real x3)
00220  *  \brief Gravitational potential */
00221 static Real grav_pot(const Real x1, const Real x2, const Real x3) {
00222   return 0.5*SQR(x1*omega0);
00223 }
00224 
00225 /*! \fn static Real grav_acc(const Real x1, const Real x2, const Real x3) 
00226  *  \brief Gravitational acceleration */
00227 static Real grav_acc(const Real x1, const Real x2, const Real x3) {
00228   return x1*SQR(omega0);
00229 }
00230 
00231 /*----------------------------------------------------------------------------*/
00232 /*! \fn void cyladvect_ix1(GridS *pG)
00233  *  \brief Inner-R boundary conditions.  d, M2, B1, B1i, and P are all
00234  *   functions of R, phi, and t.
00235  */
00236 
00237 void cyladvect_ix1(GridS *pG)
00238 {
00239   int is = pG->is;
00240   int js = pG->js, je = pG->je;
00241   int ks = pG->ks, ke = pG->ke;
00242   int i,j,k;
00243   Real Eint,Emag,Ekin;
00244 
00245   phi0 = alpha*(omega0*pG->time + x2min);
00246   for (k=ks; k<=ke; k++) {
00247     for (j=js; j<=je; j++) {
00248       for (i=1; i<=nghost; i++) {
00249         pG->U[k][j][is-i].d   = avg2d(d,pG,is-i,j,k);
00250         pG->U[k][j][is-i].M2  = avg2d(M2,pG,is-i,j,k);
00251         pG->U[k][j][is-i].M3  = pG->U[k][j][is-i].d*vz0;
00252 #ifdef MHD
00253         pG->U[k][j][is-i].B3c = bz0;
00254         pG->B3c[k][j][is-i]   = bz0;
00255 #endif
00256 
00257 #ifndef ISOTHERMAL
00258         Emag = 0.0;
00259 #ifdef MHD
00260         Emag = 0.5*SQR(pG->U[k][j][is-i].B3c);
00261 #endif
00262         Eint = (Pgas0-Emag)/Gamma_1;
00263         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;
00264 
00265         pG->U[k][j][is-i].E = Eint + Emag + Ekin;
00266 #endif /* ISOTHERMAL */
00267       }
00268     }
00269   }
00270 
00271   return;
00272 }
00273 
00274 /*----------------------------------------------------------------------------*/
00275 /*! \fn void cyladvect_ox1(GridS *pG)
00276  *  \brief B_R = B_0/R boundary conditions, Outer x1 boundary
00277  */
00278 
00279 void cyladvect_ox1(GridS *pG)
00280 {
00281   int ie = pG->ie;
00282   int js = pG->js, je = pG->je;
00283   int ks = pG->ks, ke = pG->ke;
00284   int i,j,k;
00285   Real Eint,Emag,Ekin;
00286 
00287   phi0 = alpha*(omega0*pG->time + x2min);
00288   for (k=ks; k<=ke; k++) {
00289     for (j=js; j<=je; j++) {
00290       for (i=1; i<=nghost; i++) {
00291         pG->U[k][j][ie+i].d   = avg2d(d,pG,ie+i,j,k);
00292         pG->U[k][j][ie+i].M2  = avg2d(M2,pG,ie+i,j,k);
00293         pG->U[k][j][ie+i].M3  = pG->U[k][j][ie+i].d*vz0;
00294 #ifdef MHD
00295         pG->U[k][j][ie+i].B3c = bz0;
00296         pG->B3c[k][j][ie+i]   = bz0;
00297 #endif
00298 
00299 #ifndef ISOTHERMAL
00300         Emag = 0.0;
00301 #ifdef MHD
00302         Emag = 0.5*SQR(pG->U[k][j][ie+i].B3c);
00303 #endif
00304         Eint = (Pgas0-Emag)/Gamma_1;
00305         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;
00306 
00307         pG->U[k][j][ie+i].E = Eint + Emag + Ekin;
00308 #endif /* ISOTHERMAL */
00309       }
00310     }
00311   }
00312 
00313   return;
00314 }