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
00019 static Real q, r0, rhomax, r_in, rho0, e0, dcut, beta, seed;
00020
00021
00022 static Real f, C, Kbar, n;
00023
00024
00025
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
00031 }
00032
00033
00034
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
00042
00043
00044
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
00053 return d;
00054 }
00055
00056
00057
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
00062 return x1*pG->dx1*pG->dx2*pG->dx3;
00063 }
00064
00065
00066
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
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
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
00130
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
00159 Ta = f*pow(x1, -1.0*q1) - C;
00160 T = (1.0/Ta) + 2.0;
00161
00162 if ( (x1 >= r_in) && (rad <= T) ) {
00163 rhoa = density(x1, x2, x3);
00164 IntE = pow(rhoa,Gamma)*Kbar/((n+1.0)*Gamma_1);
00165
00166 IntE = IntE*(1.0 - seed*sin(x2));
00167 if ((IntE >= e0) && (rhoa >= rho0)) {
00168
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
00174
00175 }
00176 }
00177 }
00178 }
00179 }
00180
00181
00182
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
00196 Ap[k][j][i] = SQR(rhoa-dcut);
00197 }
00198 }
00199 }
00200 }
00201
00202
00203
00204 #ifdef MHD
00205
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
00214 pG->B1i[k][j][i] = -(Ap[k+1][j][i] - Ap[k][j][i])/pG->dx3;
00215
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
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
00240
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
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
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
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
00312
00313
00314
00315
00316
00317
00318
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
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
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;
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 }