00001 #include "copyright.h"
00002
00003
00004
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 static Real br, bphi, omega, rho, pgas;
00019 static int iprob;
00020
00021
00022
00023
00024 static Real grav_pot(const Real x1, const Real x2, const Real x3) {
00025 switch(iprob) {
00026 case 0:
00027 return 0.0;
00028 break;
00029 case 1:
00030 return 0.5*SQR(x1*omega);
00031 break;
00032 case 2:
00033 return -1.0/(rho*x1) + 0.5*SQR(x1*omega);
00034 break;
00035 case 3:
00036 return x1/rho + 0.5*SQR(x1*omega);
00037 break;
00038 case 4:
00039 return 0.0;
00040 break;
00041 default:
00042 ath_error("[cylunif]: Invalid problem number!\n");
00043 }
00044 }
00045
00046
00047
00048 static Real grav_acc(const Real x1, const Real x2, const Real x3) {
00049 switch(iprob) {
00050 case 0:
00051 return 0.0;
00052 break;
00053 case 1:
00054 return x1*SQR(omega);
00055 break;
00056 case 2:
00057 return 1.0/(rho*SQR(x1)) + x1*SQR(omega);
00058 break;
00059 case 3:
00060 return 1.0/rho + x1*SQR(omega);
00061 break;
00062 case 4:
00063 return 0.0;
00064 break;
00065 default:
00066 ath_error("[cylunif]: Invalid problem number!\n");
00067 }
00068 }
00069
00070 static Gas ***Soln=NULL;
00071
00072
00073
00074
00075
00076
00077 void problem(Grid *pG, Domain *pDomain)
00078 {
00079 int i,j,k;
00080 int is,ie,il,iu,js,je,jl,ju,ks,ke,kl,ku;
00081 int nx1, nx2, nx3;
00082 Real x1, x2, x3, y1, y2, y3, r, noise, r1, r2;
00083 Real x1min, x1max, x2min, x2max, pb;
00084 Real ptotal;
00085
00086 is = pG->is; ie = pG->ie; nx1 = ie-is+1;
00087 js = pG->js; je = pG->je; nx2 = je-js+1;
00088 ks = pG->ks; ke = pG->ke; nx3 = ke-ks+1;
00089
00090 il = is-nghost*(nx1>1); iu = ie+nghost*(nx1>1); nx1 = iu-il+1;
00091 jl = js-nghost*(nx2>1); ju = je+nghost*(nx2>1); nx2 = ju-jl+1;
00092 kl = ks-nghost*(nx3>1); ku = ke+nghost*(nx3>1); nx3 = ku-kl+1;
00093
00094 #ifndef CYLINDRICAL
00095 ath_error("[cylunif]: This problem only works in cylindrical!\n");
00096 #endif
00097
00098
00099 if ((Soln = (Gas***)calloc_3d_array(nx3,nx2,nx1,sizeof(Gas))) == NULL)
00100 ath_error("[cylunif]: Error allocating memory for solution\n");
00101
00102
00103 omega = par_getd("problem", "omega");
00104 br = par_getd("problem", "br");
00105 bphi = par_getd("problem", "bphi");
00106 rho = par_getd("problem", "rho");
00107 pgas = par_getd("problem", "pgas");
00108 iprob = par_geti("problem","iprob");
00109
00110
00111
00112
00113
00114 for (k=kl; k<=ku; k++) {
00115 for (j=jl; j<=ju; j++) {
00116 for (i=il; i<=iu; i++) {
00117 cc_pos(pG,i,j,k,&x1,&x2,&x3);
00118 vc_pos(pG,i,j,k,&y1,&y2,&y3);
00119 memset(&(pG->U[k][j][i]),0.0,sizeof(Gas));
00120
00121
00122 switch(iprob) {
00123 case 0:
00124 pG->U[k][j][i].d = rho;
00125 pG->U[k][j][i].E = pgas/Gamma_1;
00126 break;
00127 case 1:
00128 pG->U[k][j][i].d = rho;
00129 pG->U[k][j][i].M2 = pG->U[k][j][i].d*y1*omega;
00130 pG->U[k][j][i].E = pgas/Gamma_1;
00131 break;
00132 case 2:
00133 pG->U[k][j][i].d = rho;
00134 pG->U[k][j][i].M2 = pG->U[k][j][i].d*y1*omega;
00135 pG->U[k][j][i].E = (pgas/x1)/Gamma_1 + 0.5*SQR(pG->U[k][j][i].M2)/pG->U[k][j][i].d;
00136 break;
00137 case 3:
00138 pG->U[k][j][i].d = rho;
00139 pG->U[k][j][i].M2 = pG->U[k][j][i].d*y1*omega;
00140
00141 pG->U[k][j][i].E = (pgas+5.0-y1)/Gamma_1 + 0.5*SQR(pG->U[k][j][i].M2)/pG->U[k][j][i].d;
00142 #ifdef MHD
00143 pG->U[k][j][i].B1c = br/x1;
00144 pG->B1i[k][j][i] = br/(x1-0.5*pG->dx1);
00145 pG->U[k][j][i].E += 0.5*SQR(pG->U[k][j][i].B1c);
00146 #endif
00147 break;
00148 case 4:
00149 pG->U[k][j][i].d = rho*(SQR(x1)+0.25*SQR(pG->dx1));
00150 pG->U[k][j][i].E = pgas/Gamma_1;
00151 break;
00152 default:
00153 ath_error("[cylunif]: Invalid problem number!\n");
00154 }
00155
00156
00157 Soln[k][j][i] = pG->U[k][j][i];
00158 }
00159 }
00160 }
00161
00162
00163
00164 StaticGravPot = grav_pot;
00165 x1GravAcc = grav_acc;
00166 set_bvals_mhd_fun(left_x1,do_nothing_bc);
00167 set_bvals_mhd_fun(right_x1,do_nothing_bc);
00168
00169 return;
00170 }
00171
00172
00173
00174
00175
00176
00177
00178
00179
00180
00181
00182
00183
00184
00185 void problem_write_restart(Grid *pG, Domain *pD, FILE *fp)
00186 {
00187 return;
00188 }
00189
00190 void problem_read_restart(Grid *pG, Domain *pD, FILE *fp)
00191 {
00192 return;
00193 }
00194
00195 Gasfun_t get_usr_expr(const char *expr)
00196 {
00197 return NULL;
00198 }
00199
00200 VGFunout_t get_usr_out_fun(const char *name){
00201 return NULL;
00202 }
00203
00204 void Userwork_in_loop(Grid *pGrid, Domain *pDomain)
00205 {
00206 }
00207
00208 void Userwork_after_loop(Grid *pGrid, Domain *pDomain)
00209 {
00210 compute_l1_error("CylUnif", pGrid, pDomain, Soln, 0);
00211 }