00001 #include "copyright.h"
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021
00022
00023
00024
00025 #include <float.h>
00026 #include <math.h>
00027
00028 #include <stdlib.h>
00029 #include <string.h>
00030 #include "defs.h"
00031 #include "athena.h"
00032 #include "globals.h"
00033 #include "prototypes.h"
00034
00035
00036
00037
00038
00039
00040
00041
00042
00043
00044 static double ran2(long int *idum);
00045 static Real UnstratifiedDisk(const Real x1, const Real x2, const Real x3);
00046 static Real VertGrav(const Real x1, const Real x2, const Real x3);
00047 static Real expr_dV2(const GridS *pG, const int i, const int j, const int k);
00048 static Real expr_beta(const GridS *pG, const int i, const int j, const int k);
00049 static Real expr_ME(const GridS *pG, const int i, const int j, const int k);
00050 static Real expr_KE(const GridS *pG, const int i, const int j, const int k);
00051 static Real hst_rho_Vx_dVy(const GridS *pG,const int i,const int j,const int k);
00052 static Real hst_rho_dVy2(const GridS *pG, const int i, const int j, const int k);
00053 #ifdef ADIABATIC
00054 static Real hst_E_total(const GridS *pG, const int i, const int j, const int k);
00055 #endif
00056 #ifdef MHD
00057 static Real hst_Bx(const GridS *pG, const int i, const int j, const int k);
00058 static Real hst_By(const GridS *pG, const int i, const int j, const int k);
00059 static Real hst_Bz(const GridS *pG, const int i, const int j, const int k);
00060 static Real hst_BxBy(const GridS *pG, const int i, const int j, const int k);
00061 #endif
00062
00063
00064 static Real ztop, zbtm;
00065
00066
00067
00068
00069
00070
00071 void problem(DomainS *pDomain)
00072 {
00073 GridS *pGrid = (pDomain->Grid);
00074 FILE *fp;
00075 Real xFP[160],dFP[160],vxFP[160],vyFP[160];
00076 int is = pGrid->is, ie = pGrid->ie;
00077 int js = pGrid->js, je = pGrid->je;
00078 int ks = pGrid->ks, ke = pGrid->ke;
00079 int ixs,jxs,kxs,i,j,k,ipert,ifield;
00080 long int iseed = -1;
00081 Real x1,x2,x3,xmin,xmax,Lx,Ly,Lz;
00082 Real den=1.0, pres=1.0e-6, rd, rp, rvx, rvy, rvz;
00083 Real beta,B0,kx,amp;
00084 int nwx,nwy,nwz;
00085 double rval;
00086 static int frst=1;
00087
00088 if (pGrid->Nx[1] == 1){
00089 ath_error("[problem]: HGB only works on a 2D or 3D grid\n");
00090 }
00091
00092
00093 #ifdef ISOTHERMAL
00094 pres=den*Iso_csound2;
00095 #endif
00096 Omega_0 = par_getd_def("problem","omega",1.0e-3);
00097 qshear = par_getd_def("problem","qshear",1.5);
00098 amp = par_getd("problem","amp");
00099 beta = par_getd("problem","beta");
00100 B0 = sqrt((double)(2.0*pres/beta));
00101 ifield = par_geti_def("problem","ifield", 1);
00102 ipert = par_geti_def("problem","ipert", 1);
00103
00104
00105 ixs = pGrid->Disp[0];
00106 jxs = pGrid->Disp[1];
00107 kxs = pGrid->Disp[2];
00108 iseed = -1 - (ixs + pDomain->Nx[0]*(jxs + pDomain->Nx[1]*kxs));
00109
00110
00111 ztop = pDomain->RootMaxX[2];
00112 zbtm = pDomain->RootMinX[2];
00113 Lx = pDomain->RootMaxX[0] - pDomain->RootMinX[0];
00114 Ly = pDomain->RootMaxX[1] - pDomain->RootMinX[1];
00115 Lz = pDomain->RootMaxX[2] - pDomain->RootMinX[2];
00116
00117
00118 nwx = par_geti_def("problem","nwx",1);
00119 kx = (2.0*PI/Lx)*((double)nwx);
00120
00121 for (k=ks; k<=ke; k++) {
00122 for (j=js; j<=je; j++) {
00123 for (i=is; i<=ie; i++) {
00124 cc_pos(pGrid,i,j,k,&x1,&x2,&x3);
00125
00126
00127
00128
00129 if (ipert == 1) {
00130 rval = amp*(ran2(&iseed) - 0.5);
00131 rd = den*exp(-x3*x3);
00132 #ifdef ADIABATIC
00133 rp = pres*rd*(1.0 + 2.0*rval);
00134 #else
00135 rd *= (1.0 + 2.0*rval);
00136 #endif
00137
00138 rval = amp*(ran2(&iseed) - 0.5);
00139 rvx = 0.4*rval*sqrt(pres/den);
00140
00141 rval = amp*(ran2(&iseed) - 0.5);
00142 rvy = 0.4*rval*sqrt(pres/den);
00143
00144 rval = amp*(ran2(&iseed) - 0.5);
00145 rvz = 0.4*rval*sqrt(pres/den);
00146 }
00147
00148
00149
00150
00151 pGrid->U[k][j][i].d = rd;
00152 pGrid->U[k][j][i].M1 = rd*rvx;
00153 pGrid->U[k][j][i].M2 = rd*rvy;
00154 #ifndef FARGO
00155 pGrid->U[k][j][i].M2 -= rd*(qshear*Omega_0*x1);
00156 #endif
00157 pGrid->U[k][j][i].M3 = rd*rvz;
00158 #ifdef ADIABATIC
00159 pGrid->U[k][j][i].E = rp/Gamma_1
00160 + 0.5*(SQR(pGrid->U[k][j][i].M1) + SQR(pGrid->U[k][j][i].M2)
00161 + SQR(pGrid->U[k][j][i].M3))/rd;
00162 #endif
00163
00164
00165
00166
00167
00168
00169
00170 #ifdef MHD
00171 pGrid->B1i[k][j][i] = 0.0;
00172 pGrid->B2i[k][j][i] = 0.0;
00173 pGrid->B3i[k][j][i] = 0.0;
00174 if (i==ie) pGrid->B1i[k][j][ie+1] = 0.0;
00175 if (j==je) pGrid->B2i[k][je+1][i] = 0.0;
00176 if (k==ke) pGrid->B3i[ke+1][j][i] = 0.0;
00177
00178 if (ifield == 1) {
00179 pGrid->B1i[k][j][i] = 0.0;
00180 pGrid->B2i[k][j][i] = 0.0;
00181 pGrid->B3i[k][j][i] = B0*(sin((double)kx*x1));
00182 if (i==ie) pGrid->B1i[k][j][ie+1] = 0.0;
00183 if (j==je) pGrid->B2i[k][je+1][i] = 0.0;
00184 if (k==ke) pGrid->B3i[ke+1][j][i] = B0*(sin((double)kx*x1));
00185 }
00186 if (ifield == 2) {
00187 pGrid->B1i[k][j][i] = 0.0;
00188 pGrid->B2i[k][j][i] = 0.0;
00189 pGrid->B3i[k][j][i] = B0;
00190 if (i==ie) pGrid->B1i[k][j][ie+1] = 0.0;
00191 if (j==je) pGrid->B2i[k][je+1][i] = 0.0;
00192 if (k==ke) pGrid->B3i[ke+1][j][i] = B0;
00193 }
00194 if (ifield == 3) {
00195 pGrid->B1i[k][j][i] = 0.0;
00196 pGrid->B2i[k][j][i] = B0*(cos((double)kx*x1));
00197 pGrid->B3i[k][j][i] = B0*(sin((double)kx*x1));
00198 if (i==ie) pGrid->B1i[k][j][ie+1] = 0.0;
00199 if (j==je) pGrid->B2i[k][je+1][i] = B0*(cos((double)kx*x1));
00200 if (k==ke) pGrid->B3i[ke+1][j][i] = B0*(sin((double)kx*x1));
00201 }
00202 if (ifield == 4 && fabs(x3) < 2.0) {
00203 pGrid->B1i[k][j][i] = 0.0;
00204 pGrid->B2i[k][j][i] = B0;
00205 pGrid->B3i[k][j][i] = 0.0;
00206 if (i==ie) pGrid->B1i[k][j][ie+1] = 0.0;
00207 if (j==je) pGrid->B2i[k][je+1][i] = B0;
00208 if (k==ke) pGrid->B3i[ke+1][j][i] = 0.0;
00209 }
00210
00211 #ifdef ADIABATIC
00212 pGrid->U[k][j][i].E += 0.5*(SQR(pGrid->U[k][j][i].B1c)
00213 + SQR(pGrid->U[k][j][i].B2c) + SQR(pGrid->U[k][j][i].B3c));
00214 #endif
00215 #endif
00216 }
00217 }}
00218 #ifdef MHD
00219 for (k=ks; k<=ke; k++) {
00220 for (j=js; j<=je; j++) {
00221 for (i=is; i<=ie; i++) {
00222 pGrid->U[k][j][i].B1c = 0.5*(pGrid->B1i[k][j][i]+pGrid->B1i[k][j][i+1]);
00223 pGrid->U[k][j][i].B2c = 0.5*(pGrid->B2i[k][j][i]+pGrid->B2i[k][j+1][i]);
00224 pGrid->U[k][j][i].B3c = 0.5*(pGrid->B3i[k][j][i]+pGrid->B3i[k+1][j][i]);
00225 }
00226 }
00227 }
00228 #endif
00229
00230
00231 #ifdef OHMIC
00232 eta_Ohm = par_getd("problem","eta");
00233 #endif
00234 #ifdef NAVIER_STOKES
00235 nu_V = par_getd("problem","nu");
00236 #endif
00237
00238
00239
00240 StaticGravPot = VertGrav;
00241 ShearingBoxPot = UnstratifiedDisk;
00242
00243
00244
00245 if (frst == 1) {
00246 dump_history_enroll(hst_rho_Vx_dVy, "<rho Vx dVy>");
00247 dump_history_enroll(hst_rho_dVy2, "<rho dVy^2>");
00248 #ifdef ADIABATIC
00249 dump_history_enroll(hst_E_total, "<E + rho Phi>");
00250 #endif
00251 #ifdef MHD
00252 dump_history_enroll(hst_Bx, "<Bx>");
00253 dump_history_enroll(hst_By, "<By>");
00254 dump_history_enroll(hst_Bz, "<Bz>");
00255 dump_history_enroll(hst_BxBy, "<-Bx By>");
00256 #endif
00257 frst = 0;
00258 }
00259
00260 return;
00261 }
00262
00263
00264
00265
00266
00267
00268
00269
00270
00271
00272
00273
00274 void problem_write_restart(MeshS *pM, FILE *fp)
00275 {
00276 return;
00277 }
00278
00279
00280
00281
00282
00283 void problem_read_restart(MeshS *pM, FILE *fp)
00284 {
00285 Omega_0 = par_getd_def("problem","omega",1.0e-3);
00286 qshear = par_getd_def("problem","qshear",1.5);
00287
00288
00289 #ifdef OHMIC
00290 eta_Ohm = par_getd("problem","eta");
00291 #endif
00292 #ifdef NAVIER_STOKES
00293 nu_V = par_getd("problem","nu");
00294 #endif
00295
00296
00297
00298 StaticGravPot = VertGrav;
00299 ShearingBoxPot = UnstratifiedDisk;
00300
00301
00302
00303 dump_history_enroll(hst_rho_Vx_dVy, "<rho Vx dVy>");
00304 dump_history_enroll(hst_rho_dVy2, "<rho dVy^2>");
00305 #ifdef ADIABATIC
00306 dump_history_enroll(hst_E_total, "<E + rho Phi>");
00307 #endif
00308 #ifdef MHD
00309 dump_history_enroll(hst_Bx, "<Bx>");
00310 dump_history_enroll(hst_By, "<By>");
00311 dump_history_enroll(hst_Bz, "<Bz>");
00312 dump_history_enroll(hst_BxBy, "<-Bx By>");
00313 #endif
00314
00315 return;
00316 }
00317
00318
00319 ConsFun_t get_usr_expr(const char *expr)
00320 {
00321 if(strcmp(expr,"dVy")==0) return expr_dV2;
00322 else if(strcmp(expr,"beta")==0) return expr_beta;
00323 else if(strcmp(expr,"ME")==0) return expr_ME;
00324 else if(strcmp(expr,"KE")==0) return expr_KE;
00325 else if(strcmp(expr,"BxBy")==0) return hst_BxBy;
00326 return NULL;
00327 }
00328
00329 VOutFun_t get_usr_out_fun(const char *name){
00330 return NULL;
00331 }
00332
00333 void Userwork_in_loop(MeshS *pM)
00334 {
00335 }
00336
00337 void Userwork_after_loop(MeshS *pM)
00338 {
00339 }
00340
00341
00342
00343
00344 #define IM1 2147483563
00345 #define IM2 2147483399
00346 #define AM (1.0/IM1)
00347 #define IMM1 (IM1-1)
00348 #define IA1 40014
00349 #define IA2 40692
00350 #define IQ1 53668
00351 #define IQ2 52774
00352 #define IR1 12211
00353 #define IR2 3791
00354 #define NTAB 32
00355 #define NDIV (1+IMM1/NTAB)
00356 #define RNMX (1.0-DBL_EPSILON)
00357
00358
00359
00360
00361
00362
00363
00364
00365
00366
00367
00368
00369
00370
00371 double ran2(long int *idum)
00372 {
00373 int j;
00374 long int k;
00375 static long int idum2=123456789;
00376 static long int iy=0;
00377 static long int iv[NTAB];
00378 double temp;
00379
00380 if (*idum <= 0) {
00381 if (-(*idum) < 1) *idum=1;
00382 else *idum = -(*idum);
00383 idum2=(*idum);
00384 for (j=NTAB+7;j>=0;j--) {
00385 k=(*idum)/IQ1;
00386 *idum=IA1*(*idum-k*IQ1)-k*IR1;
00387 if (*idum < 0) *idum += IM1;
00388 if (j < NTAB) iv[j] = *idum;
00389 }
00390 iy=iv[0];
00391 }
00392 k=(*idum)/IQ1;
00393 *idum=IA1*(*idum-k*IQ1)-k*IR1;
00394 if (*idum < 0) *idum += IM1;
00395 k=idum2/IQ2;
00396 idum2=IA2*(idum2-k*IQ2)-k*IR2;
00397 if (idum2 < 0) idum2 += IM2;
00398 j=(int)(iy/NDIV);
00399 iy=iv[j]-idum2;
00400 iv[j] = *idum;
00401 if (iy < 1) iy += IMM1;
00402 if ((temp=AM*iy) > RNMX) return RNMX;
00403 else return temp;
00404 }
00405
00406 #undef IM1
00407 #undef IM2
00408 #undef AM
00409 #undef IMM1
00410 #undef IA1
00411 #undef IA2
00412 #undef IQ1
00413 #undef IQ2
00414 #undef IR1
00415 #undef IR2
00416 #undef NTAB
00417 #undef NDIV
00418 #undef RNMX
00419
00420
00421
00422
00423 static Real UnstratifiedDisk(const Real x1, const Real x2, const Real x3)
00424 {
00425 Real phi=0.0;
00426 #ifndef FARGO
00427 phi -= qshear*Omega_0*Omega_0*x1*x1;
00428 #endif
00429 return phi;
00430 }
00431
00432
00433
00434
00435 static Real VertGrav(const Real x1, const Real x2, const Real x3)
00436 {
00437 Real phi=0.0,z;
00438
00439
00440 if(x3 > ztop)
00441 z=x3-ztop+zbtm;
00442 else if (x3 < zbtm)
00443 z=x3-zbtm+ztop;
00444 else
00445 z=x3;
00446
00447 phi += 0.5*Omega_0*Omega_0*
00448 (SQR(fabs(ztop)-sqrt(SQR(fabs(ztop)-fabs(z)) + 0.01)));
00449
00450 return phi;
00451 }
00452
00453
00454
00455
00456
00457
00458 static Real expr_dV2(const GridS *pG, const int i, const int j, const int k)
00459 {
00460 Real x1,x2,x3;
00461 cc_pos(pG,i,j,k,&x1,&x2,&x3);
00462 #ifdef FARGO
00463 return (pG->U[k][j][i].M2/pG->U[k][j][i].d);
00464 #else
00465 return (pG->U[k][j][i].M2/pG->U[k][j][i].d + qshear*Omega_0*x1);
00466 #endif
00467 }
00468
00469
00470
00471
00472
00473
00474 static Real expr_beta(const GridS *pG, const int i, const int j, const int k)
00475 {
00476 Real x1,x2,x3,B2;
00477 cc_pos(pG,i,j,k,&x1,&x2,&x3);
00478 #ifdef MHD
00479 B2=pG->U[k][j][i].B1c*pG->U[k][j][i].B1c;
00480 B2+=pG->U[k][j][i].B2c*pG->U[k][j][i].B2c;
00481 B2+=pG->U[k][j][i].B3c*pG->U[k][j][i].B3c;
00482
00483 #ifdef ISOTHERMAL
00484 return (2.0*Iso_csound2*pG->U[k][j][i].d/B2);
00485 #else
00486 return 0.0;
00487 #endif
00488
00489 #else
00490 return 0.0;
00491 #endif
00492 }
00493
00494
00495
00496
00497
00498
00499 static Real expr_ME(const GridS *pG, const int i, const int j, const int k)
00500 {
00501 Real x1,x2,x3,B2;
00502 cc_pos(pG,i,j,k,&x1,&x2,&x3);
00503 #ifdef MHD
00504 B2=pG->U[k][j][i].B1c*pG->U[k][j][i].B1c;
00505 B2+=pG->U[k][j][i].B2c*pG->U[k][j][i].B2c;
00506 B2+=pG->U[k][j][i].B3c*pG->U[k][j][i].B3c;
00507 return (B2/2.0);
00508 #else
00509 return NULL;
00510 #endif
00511 }
00512
00513
00514
00515
00516
00517 static Real expr_KE(const GridS *pG, const int i, const int j, const int k)
00518 {
00519 Real x1,x2,x3,Vy,Vx,Vz;
00520 cc_pos(pG,i,j,k,&x1,&x2,&x3);
00521 #ifdef FARGO
00522 Vy = (pG->U[k][j][i].M2/pG->U[k][j][i].d);
00523 #else
00524 Vy = (pG->U[k][j][i].M2/pG->U[k][j][i].d + qshear*Omega_0*x1);
00525 #endif
00526 Vx = pG->U[k][j][i].M1/pG->U[k][j][i].d;
00527 Vz = pG->U[k][j][i].M3/pG->U[k][j][i].d;
00528
00529 return pG->U[k][j][i].d*(Vx*Vx + Vy*Vy + Vz*Vz)/2.0;
00530
00531 }
00532
00533
00534
00535
00536
00537
00538
00539
00540
00541
00542 static Real hst_rho_Vx_dVy(const GridS *pG,const int i,const int j, const int k)
00543 {
00544 Real x1,x2,x3;
00545 cc_pos(pG,i,j,k,&x1,&x2,&x3);
00546 #ifdef FARGO
00547 return pG->U[k][j][i].M1*(pG->U[k][j][i].M2/pG->U[k][j][i].d);
00548 #else
00549 return pG->U[k][j][i].M1*
00550 (pG->U[k][j][i].M2/pG->U[k][j][i].d + qshear*Omega_0*x1);
00551 #endif
00552 }
00553
00554
00555
00556
00557 static Real hst_rho_dVy2(const GridS *pG, const int i, const int j, const int k)
00558 {
00559 Real x1,x2,x3,dVy;
00560 cc_pos(pG,i,j,k,&x1,&x2,&x3);
00561 #ifdef FARGO
00562 dVy = (pG->U[k][j][i].M2/pG->U[k][j][i].d);
00563 #else
00564 dVy = (pG->U[k][j][i].M2/pG->U[k][j][i].d + qshear*Omega_0*x1);
00565 #endif
00566 return pG->U[k][j][i].d*dVy*dVy;
00567 }
00568
00569 #ifdef ADIABATIC
00570
00571
00572
00573 static Real hst_E_total(const GridS *pG, const int i, const int j, const int k)
00574 {
00575 Real x1,x2,x3,phi;
00576 cc_pos(pG,i,j,k,&x1,&x2,&x3);
00577 phi = UnstratifiedDisk(x1, x2, x3);
00578
00579 return pG->U[k][j][i].E + pG->U[k][j][i].d*phi;
00580 }
00581 #endif
00582
00583
00584
00585
00586
00587
00588 #ifdef MHD
00589
00590
00591 static Real hst_Bx(const GridS *pG, const int i, const int j, const int k)
00592 {
00593 return pG->U[k][j][i].B1c;
00594 }
00595
00596
00597
00598 static Real hst_By(const GridS *pG, const int i, const int j, const int k)
00599 {
00600 return pG->U[k][j][i].B2c;
00601 }
00602
00603
00604
00605
00606 static Real hst_Bz(const GridS *pG, const int i, const int j, const int k)
00607 {
00608 return pG->U[k][j][i].B3c;
00609 }
00610
00611
00612
00613
00614 static Real hst_BxBy(const GridS *pG, const int i, const int j, const int k)
00615 {
00616 return -pG->U[k][j][i].B1c*pG->U[k][j][i].B2c;
00617 }
00618 #endif