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
00026
00027
00028
00029
00030
00031
00032
00033
00034
00035
00036 #include <float.h>
00037 #include <math.h>
00038 #include <stdio.h>
00039 #include <stdlib.h>
00040 #include <string.h>
00041 #include "defs.h"
00042 #include "athena.h"
00043 #include "globals.h"
00044 #include "prototypes.h"
00045
00046 #ifndef SHEARING_BOX
00047 #error : The HB3 problem requires shearing-box to be enabled.
00048 #endif
00049
00050
00051
00052
00053
00054
00055
00056
00057
00058
00059
00060
00061
00062
00063
00064 static double ran2(long int *idum);
00065 static Real UnstratifiedDisk(const Real x1, const Real x2, const Real x3);
00066 static Real expr_dV3(const GridS *pG, const int i, const int j, const int k);
00067 static Real hst_rho_Vx_dVy(const GridS *pG,const int i,const int j,const int k);
00068 static Real hst_dEk(const GridS *pG, const int i, const int j, const int k);
00069 static Real hst_E_total(const GridS *pG, const int i, const int j, const int k);
00070 #ifdef MHD
00071 static Real hst_Bx(const GridS *pG, const int i, const int j, const int k);
00072 static Real hst_By(const GridS *pG, const int i, const int j, const int k);
00073 static Real hst_Bz(const GridS *pG, const int i, const int j, const int k);
00074 static Real hst_BxBy(const GridS *pG, const int i, const int j, const int k);
00075 #endif
00076
00077
00078 static Real Lx,Lz;
00079
00080
00081
00082
00083
00084 void problem(DomainS *pDomain)
00085 {
00086 GridS *pGrid=(pDomain->Grid);
00087 int is = pGrid->is, ie = pGrid->ie;
00088 int js = pGrid->js, je = pGrid->je;
00089 int ks = pGrid->ks;
00090 int i,j,ipert,ifield;
00091 long int iseed = -1;
00092 Real x1, x2, x3, x1min, x1max, x2min, x2max;
00093 Real den = 1.0, pres = 1.0e-5, rd, rp, rvx;
00094 Real beta,B0,kx,kz,amp;
00095 double rval;
00096 static int frst=1;
00097
00098
00099 ShBoxCoord = xz;
00100
00101 if (pGrid->Nx[1] == 1){
00102 ath_error("[problem]: HB3 only works on a 2D grid\n");
00103 }
00104
00105 if (pGrid->Nx[2] > 1){
00106 ath_error("[problem]: HB3 does not work on 3D grid\n");
00107 }
00108
00109
00110 x1min = pDomain->RootMinX[0];
00111 x1max = pDomain->RootMaxX[0];
00112 Lx = x1max - x1min;
00113 kx = 2.0*PI/Lx;
00114
00115 x2min = pDomain->RootMinX[1];
00116 x2max = pDomain->RootMaxX[1];
00117 Lz = x2max - x2min;
00118 kz = 2.0*PI/Lz;
00119
00120
00121 Omega_0 = par_getd_def("problem","omega",1.0e-3);
00122 qshear = par_getd_def("problem","qshear",1.5);
00123 amp = par_getd("problem","amp");
00124 beta = par_getd("problem","beta");
00125 B0 = sqrt((double)(2.0*pres/beta));
00126 ifield = par_geti_def("problem","ifield", 1);
00127 ipert = par_geti_def("problem","ipert", 1);
00128
00129 for (j=js; j<=je; j++) {
00130 for (i=is; i<=ie; i++) {
00131 cc_pos(pGrid,i,j,ks,&x1,&x2,&x3);
00132
00133
00134
00135
00136
00137
00138
00139 if (ipert == 1) {
00140 rval = 1.0 + amp*(ran2(&iseed) - 0.5);
00141 #ifdef ADIABATIC
00142 rp = rval*pres;
00143 rd = pow(rval,1.0/Gamma)*den;
00144 #else
00145 rd = rval*den;
00146 #endif
00147 rvx = 0.0;
00148 }
00149 if (ipert == 2) {
00150 rp = pres;
00151 rd = den*(1.0 + 0.1*sin((double)kx*x1));
00152 #ifdef ADIABATIC
00153 rvx = amp*sqrt(Gamma*pres/den);
00154 #else
00155 rvx = amp*sqrt(pres/den);
00156 #endif
00157 }
00158 if (ipert == 3) {
00159 rval = 1.0 + amp*(ran2(&iseed) - 0.5);
00160 #ifdef ADIABATIC
00161 rp = rval*pres;
00162 rd = den;
00163 #else
00164 rd = rval*den;
00165 #endif
00166 rvx = 0.0;
00167 }
00168 if (ipert == 4) {
00169 rp = pres;
00170 rd = den;
00171 rvx = amp*sin((double)kz*x2);
00172 }
00173
00174
00175
00176 pGrid->U[ks][j][i].d = rd;
00177 pGrid->U[ks][j][i].M1 = rd*rvx;
00178 pGrid->U[ks][j][i].M2 = 0.0;
00179 #ifdef FARGO
00180 pGrid->U[ks][j][i].M3 = 0.0;
00181 #else
00182 pGrid->U[ks][j][i].M3 = -rd*qshear*Omega_0*x1;
00183 #endif
00184 #ifdef ADIABATIC
00185 pGrid->U[ks][j][i].E = rp/Gamma_1
00186 + 0.5*(SQR(pGrid->U[ks][j][i].M1) + SQR(pGrid->U[ks][j][i].M3))/rd;
00187 #endif
00188
00189
00190
00191
00192
00193 #ifdef MHD
00194 if (ifield == 1) {
00195 pGrid->U[ks][j][i].B1c = 0.0;
00196 pGrid->U[ks][j][i].B2c = B0*(sin((double)kx*x1));
00197 pGrid->U[ks][j][i].B3c = 0.0;
00198 pGrid->B1i[ks][j][i] = 0.0;
00199 pGrid->B2i[ks][j][i] = B0*(sin((double)kx*x1));
00200 pGrid->B3i[ks][j][i] = 0.0;
00201 if (i==ie) pGrid->B1i[ks][j][ie+1] = 0.0;
00202 if (j==je) pGrid->B2i[ks][je+1][i] = B0*(sin((double)kx*x1));
00203 }
00204 if (ifield == 2) {
00205 pGrid->U[ks][j][i].B1c = 0.0;
00206 pGrid->U[ks][j][i].B2c = B0;
00207 pGrid->U[ks][j][i].B3c = 0.0;
00208 pGrid->B1i[ks][j][i] = 0.0;
00209 pGrid->B2i[ks][j][i] = B0;
00210 pGrid->B3i[ks][j][i] = 0.0;
00211 if (i==ie) pGrid->B1i[ks][j][ie+1] = 0.0;
00212 if (j==je) pGrid->B2i[ks][je+1][i] = B0;
00213 }
00214 #ifdef ADIABATIC
00215 pGrid->U[ks][j][i].E += 0.5*(SQR(pGrid->U[ks][j][i].B1c)
00216 + SQR(pGrid->U[ks][j][i].B2c) + SQR(pGrid->U[ks][j][i].B3c));
00217 #endif
00218 #endif
00219 }
00220 }
00221
00222
00223 #ifdef RESISTIVITY
00224 eta_Ohm = par_getd_def("problem","eta_O",0.0);
00225 Q_Hall = par_getd_def("problem","Q_H",0.0);
00226 Q_AD = par_getd_def("problem","Q_A",0.0);
00227 #endif
00228
00229
00230
00231 ShearingBoxPot = UnstratifiedDisk;
00232
00233
00234
00235 if (frst == 1) {
00236 dump_history_enroll(hst_dEk, "<0.5rho(Vx^2+4dVy^2)>");
00237 dump_history_enroll(hst_rho_Vx_dVy, "<rho Vx dVy>");
00238 dump_history_enroll(hst_E_total, "<E + rho Phi>");
00239 #ifdef MHD
00240 dump_history_enroll(hst_Bx, "<Bx>");
00241 dump_history_enroll(hst_By, "<By>");
00242 dump_history_enroll(hst_Bz, "<Bz>");
00243 dump_history_enroll(hst_BxBy, "<-Bx By>");
00244 #endif
00245 frst = 0;
00246 }
00247
00248 return;
00249 }
00250
00251
00252
00253
00254
00255
00256
00257
00258
00259
00260
00261
00262 void problem_write_restart(MeshS *pM, FILE *fp)
00263 {
00264 return;
00265 }
00266
00267
00268
00269
00270
00271
00272 void problem_read_restart(MeshS *pM, FILE *fp)
00273 {
00274 Real x1min, x1max;
00275
00276 Omega_0 = par_getd_def("problem","omega",1.0e-3);
00277 qshear = par_getd_def("problem","qshear",1.5);
00278
00279 #ifdef RESISTIVITY
00280 eta_Ohm = par_getd_def("problem","eta_O",0.0);
00281 Q_Hall = par_getd_def("problem","Q_H",0.0);
00282 Q_AD = par_getd_def("problem","Q_A",0.0);
00283 #endif
00284
00285
00286 x1min = pM->RootMinX[0];
00287 x1max = pM->RootMaxX[0];
00288 Lx = x1max - x1min;
00289
00290 ShearingBoxPot = UnstratifiedDisk;
00291
00292 return;
00293 }
00294
00295
00296 ConsFun_t get_usr_expr(const char *expr)
00297 {
00298 if(strcmp(expr,"dVy")==0) return expr_dV3;
00299 return NULL;
00300 }
00301
00302 VOutFun_t get_usr_out_fun(const char *name){
00303 return NULL;
00304 }
00305
00306 void Userwork_in_loop(MeshS *pM)
00307 {
00308 }
00309
00310 void Userwork_after_loop(MeshS *pM)
00311 {
00312 }
00313
00314
00315
00316
00317
00318 #define IM1 2147483563
00319 #define IM2 2147483399
00320 #define AM (1.0/IM1)
00321 #define IMM1 (IM1-1)
00322 #define IA1 40014
00323 #define IA2 40692
00324 #define IQ1 53668
00325 #define IQ2 52774
00326 #define IR1 12211
00327 #define IR2 3791
00328 #define NTAB 32
00329 #define NDIV (1+IMM1/NTAB)
00330 #define RNMX (1.0-DBL_EPSILON)
00331
00332
00333
00334
00335
00336
00337
00338
00339
00340
00341
00342
00343
00344
00345 double ran2(long int *idum)
00346 {
00347 int j;
00348 long int k;
00349 static long int idum2=123456789;
00350 static long int iy=0;
00351 static long int iv[NTAB];
00352 double temp;
00353
00354 if (*idum <= 0) {
00355 if (-(*idum) < 1) *idum=1;
00356 else *idum = -(*idum);
00357 idum2=(*idum);
00358 for (j=NTAB+7;j>=0;j--) {
00359 k=(*idum)/IQ1;
00360 *idum=IA1*(*idum-k*IQ1)-k*IR1;
00361 if (*idum < 0) *idum += IM1;
00362 if (j < NTAB) iv[j] = *idum;
00363 }
00364 iy=iv[0];
00365 }
00366 k=(*idum)/IQ1;
00367 *idum=IA1*(*idum-k*IQ1)-k*IR1;
00368 if (*idum < 0) *idum += IM1;
00369 k=idum2/IQ2;
00370 idum2=IA2*(idum2-k*IQ2)-k*IR2;
00371 if (idum2 < 0) idum2 += IM2;
00372 j=(int)(iy/NDIV);
00373 iy=iv[j]-idum2;
00374 iv[j] = *idum;
00375 if (iy < 1) iy += IMM1;
00376 if ((temp=AM*iy) > RNMX) return RNMX;
00377 else return temp;
00378 }
00379
00380 #undef IM1
00381 #undef IM2
00382 #undef AM
00383 #undef IMM1
00384 #undef IA1
00385 #undef IA2
00386 #undef IQ1
00387 #undef IQ2
00388 #undef IR1
00389 #undef IR2
00390 #undef NTAB
00391 #undef NDIV
00392 #undef RNMX
00393
00394
00395
00396
00397
00398 static Real UnstratifiedDisk(const Real x1, const Real x2, const Real x3){
00399 Real phi=0.0;
00400 #ifndef FARGO
00401 phi -= qshear*SQR(Omega_0*x1);
00402 #endif
00403 return phi;
00404 }
00405
00406
00407
00408
00409
00410
00411 static Real expr_dV3(const GridS *pG, const int i, const int j, const int k)
00412 {
00413 Real x1,x2,x3;
00414 cc_pos(pG,i,j,k,&x1,&x2,&x3);
00415 #ifdef FARGO
00416 return pG->U[k][j][i].M3/pG->U[k][j][i].d;
00417 #else
00418 return (pG->U[k][j][i].M3/pG->U[k][j][i].d + qshear*Omega_0*x1);
00419 #endif
00420 }
00421
00422
00423
00424
00425
00426
00427 static Real hst_rho_Vx_dVy(const GridS *pG, const int i, const int j, const int k)
00428 {
00429 Real x1,x2,x3;
00430 cc_pos(pG,i,j,k,&x1,&x2,&x3);
00431 #ifdef FARGO
00432 return pG->U[k][j][i].M1*pG->U[k][j][i].M3/pG->U[k][j][i].d;
00433 #else
00434 return pG->U[k][j][i].M1*(pG->U[k][j][i].M3/pG->U[k][j][i].d + qshear*Omega_0*x1);
00435 #endif
00436 }
00437
00438
00439
00440
00441
00442
00443 static Real hst_dEk(const GridS *pG, const int i, const int j, const int k)
00444 {
00445 Real x1,x2,x3;
00446 Real dMy, dE;
00447
00448 #ifdef FARGO
00449 dMy = pG->U[k][j][i].M3;
00450 #else
00451 dMy = (pG->U[k][j][i].M3 + qshear*Omega_0*x1*pG->U[k][j][i].d);
00452 #endif
00453 dE = 0.5*(pG->U[k][j][i].M1*pG->U[k][j][i].M1 + 4.0*dMy*dMy)/pG->U[k][j][i].d;
00454
00455 return dE;
00456 }
00457
00458
00459
00460
00461
00462
00463
00464
00465 static Real hst_E_total(const GridS *pG, const int i, const int j, const int k)
00466 {
00467 #ifdef ADIABATIC
00468 Real x1,x2,x3,phi;
00469 cc_pos(pG,i,j,k,&x1,&x2,&x3);
00470 phi = UnstratifiedDisk(x1, x2, x3);
00471
00472 return pG->U[k][j][i].E + pG->U[k][j][i].d*phi;
00473 #else
00474 return 0.0;
00475 #endif
00476 }
00477
00478
00479
00480
00481
00482 #ifdef MHD
00483
00484
00485 static Real hst_Bx(const GridS *pG, const int i, const int j, const int k)
00486 {
00487 return pG->U[k][j][pG->is].B1c;
00488 }
00489
00490
00491
00492 static Real hst_By(const GridS *pG, const int i, const int j, const int k)
00493 {
00494 return pG->U[k][pG->js][i].B2c;
00495 }
00496
00497
00498
00499 static Real hst_Bz(const GridS *pG, const int i, const int j, const int k)
00500 {
00501 return pG->U[pG->ks][j][i].B3c;
00502 }
00503
00504
00505
00506
00507 static Real hst_BxBy(const GridS *pG, const int i, const int j, const int k)
00508 {
00509 return -pG->U[k][j][i].B1c*pG->U[k][j][i].B2c;
00510 }
00511
00512 #endif