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
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 Real Lx,Ly,Lz;
00047
00048
00049
00050
00051
00052
00053
00054
00055
00056 static double ran2(long int *idum);
00057 static Real UnstratifiedDisk(const Real x1, const Real x2, const Real x3);
00058 static Real expr_dV2(const GridS *pG, const int i, const int j, const int k);
00059 static Real expr_Jsq(const GridS *pG, const int i, const int j, const int k);
00060 static Real hst_rho_Vx_dVy(const GridS *pG,const int i,const int j,const int k);
00061 static Real hst_rho_dVy2(const GridS *pG,const int i, const int j, const int k);
00062 #ifdef ADIABATIC
00063 static Real hst_E_total(const GridS *pG, const int i, const int j, const int k);
00064 #endif
00065 #ifdef MHD
00066 static Real hst_Bx(const GridS *pG, const int i, const int j, const int k);
00067 static Real hst_By(const GridS *pG, const int i, const int j, const int k);
00068 static Real hst_Bz(const GridS *pG, const int i, const int j, const int k);
00069 static Real hst_BxBy(const GridS *pG, const int i, const int j, const int k);
00070 static Real hst_dEw2(const GridS *pG, const int i, const int j, const int k);
00071 static Real hst_dBy(const GridS *pG, const int i, const int j, const int k);
00072 #endif
00073
00074
00075
00076
00077
00078
00079 void problem(DomainS *pDomain)
00080 {
00081 GridS *pGrid = pDomain->Grid;
00082 FILE *fp;
00083 Real xFP[160],dFP[160],vxFP[160],vyFP[160];
00084 int is = pGrid->is, ie = pGrid->ie;
00085 int js = pGrid->js, je = pGrid->je;
00086 int ks = pGrid->ks, ke = pGrid->ke;
00087 int ixs,jxs,kxs,i,j,k,ipert,ifield;
00088 long int iseed = -1;
00089 Real x1,x2,x3,xmin,xmax;
00090 Real den = 1.0, pres = 1.0, rd, rp, rvx, rvy, rvz, rbx, rby, rbz;
00091 Real beta=1.0,B0,kx,ky,kz,amp;
00092 int nwx,nwy,nwz;
00093 double rval;
00094 static int frst=1;
00095
00096 if (pGrid->Nx[1] == 1){
00097 ath_error("[problem]: HGB only works on a 2D or 3D grid\n");
00098 }
00099
00100
00101 Omega_0 = par_getd_def("problem","Omega",1.0e-3);
00102 qshear = par_getd_def("problem","qshear",1.5);
00103 amp = par_getd("problem","amp");
00104 #ifdef MHD
00105 beta = par_getd("problem","beta");
00106 ifield = par_geti_def("problem","ifield", 1);
00107 #endif
00108 ipert = par_geti_def("problem","ipert", 1);
00109
00110
00111 #ifdef ISOTHERMAL
00112 pres = Iso_csound2;
00113 #else
00114 pres = par_getd("problem","pres");
00115 #endif
00116 B0 = sqrt((double)(2.0*pres/beta));
00117
00118
00119 ixs = pGrid->Disp[0];
00120 jxs = pGrid->Disp[1];
00121 kxs = pGrid->Disp[2];
00122 iseed = -1 - (ixs + pDomain->Nx[0]*(jxs + pDomain->Nx[1]*kxs));
00123
00124
00125 Lx = pDomain->RootMaxX[0] - pDomain->RootMinX[0];
00126 Ly = pDomain->RootMaxX[1] - pDomain->RootMinX[1];
00127 Lz = pDomain->RootMaxX[2] - pDomain->RootMinX[2];
00128
00129
00130 nwx = par_geti_def("problem","nwx",1);
00131 nwy = par_geti_def("problem","nwy",1);
00132 nwz = par_geti_def("problem","nwz",1);
00133 kx = (2.0*PI/Lx)*((double)nwx);
00134 ky = (2.0*PI/Ly)*((double)nwy);
00135 kz = (2.0*PI/Lz)*((double)nwz);
00136
00137
00138
00139 if (ipert == 4) {
00140 if (pGrid->Nx[0] == 160) {
00141 if((fp = fopen("Data-160-FPwave.dat","r")) == NULL)
00142 ath_error("Error opening Data-160-FPwave.dat\n");
00143 for (i=0; i<160; i++) {
00144 fscanf(fp,"%lf %lf %lf %lf",&xFP[i],&dFP[i],&vxFP[i],&vyFP[i]);
00145 }
00146 }
00147
00148 if (pGrid->Nx[0] == 40) {
00149 if((fp = fopen("Data-40-FPwave.dat","r")) == NULL)
00150 ath_error("Error opening Data-40-FPwave.dat\n");
00151 for (i=0; i<40; i++) {
00152 fscanf(fp,"%lf %lf %lf %lf",&xFP[i],&dFP[i],&vxFP[i],&vyFP[i]);
00153 }
00154 }
00155
00156 xmin = pDomain->RootMinX[0];
00157 if (xmin != -4.7965) ath_error("[hgb]: iprob=4 requires xmin=-4.7965\n");
00158 xmax = pDomain->RootMaxX[0];
00159 if (xmax != 4.7965) ath_error("[hgb]: iprob=4 requires xmax=4.7965\n");
00160 }
00161
00162
00163 #ifdef ADIABATIC
00164 if (ipert == 2 || ipert == 3) amp *= sqrt(Gamma*pres/den);
00165 #else
00166 if (ipert == 2 || ipert == 3) amp *= Iso_csound;
00167 #endif
00168
00169 for (k=ks; k<=ke; k++) {
00170 for (j=js; j<=je; j++) {
00171 for (i=is; i<=ie; i++) {
00172 cc_pos(pGrid,i,j,k,&x1,&x2,&x3);
00173
00174
00175
00176
00177
00178
00179
00180
00181
00182 if (ipert == 1) {
00183 rval = amp*(ran2(&iseed) - 0.5);
00184 #ifdef ADIABATIC
00185 rp = pres*(1.0 + 2.0*rval);
00186 rd = den;
00187 #else
00188 rd = den*(1.0 + 2.0*rval);
00189 #endif
00190
00191 rval = amp*(ran2(&iseed) - 0.5);
00192 rvx = 0.4*rval*sqrt(pres/den);
00193
00194 rval = amp*(ran2(&iseed) - 0.5);
00195 rvy = 0.4*rval*sqrt(pres/den);
00196
00197 rval = amp*(ran2(&iseed) - 0.5);
00198 rvz = 0.4*rval*sqrt(pres/den);
00199 }
00200 if (ipert == 2) {
00201 rp = pres;
00202 rd = den;
00203 rvx = amp;
00204 rvy = 0.0;
00205 rvz = 0.0;
00206 }
00207 if (ipert == 3) {
00208 rp = pres;
00209 rd = den;
00210 rvx = amp*sin((double)(kx*x1 + ky*x2));
00211 rvy = -amp*(kx/ky)*sin((double)(kx*x1 + ky*x2));
00212 rvz = 0.0;
00213 }
00214 if (ipert == 4) {
00215 rd = dFP[i-is];
00216 rvx = vxFP[i-is];
00217 rvy = vyFP[i-is] + qshear*Omega_0*x1;
00218 rvz = 0.0;
00219 }
00220
00221 if (ipert == 5) {
00222 ifield = 0;
00223 rd = den + 8.9525e-10*cos((double)(kx*x1 + ky*x2 + kz*x3 - PI/4.));
00224 rvx = 8.16589e-8*cos((double)(kx*x1 + ky*x2 + kz*x3 + PI/4.));
00225 rvy = 8.70641e-8*cos((double)(kx*x1 + ky*x2 + kz*x3 + PI/4.));
00226 rvz = 0.762537e-8*cos((double)(kx*x1 + ky*x2 + kz*x3 + PI/4.));
00227 rbx = -1.08076e-7;
00228 rbx *= cos((double)(kx*(x1-0.5*pGrid->dx1) + ky*x2 + kz*x3 - PI/4.));
00229 rby = 1.04172e-7;
00230 rby *= cos((double)(kx*x1 + ky*(x2-0.5*pGrid->dx2) + kz*x3 - PI/4.));
00231 rbz = -0.320324e-7;
00232 rbz *= cos((double)(kx*x1 + ky*x2 + kz*(x3-0.5*pGrid->dx3) - PI/4.));;
00233 rbz += (sqrt(15.0)/16.0)*(Omega_0/kz);
00234 }
00235 if (ipert == 6) {
00236 ifield = 0;
00237 rd = den + 5.48082e-6*cos((double)(kx*x1 + ky*x2 + kz*x3));
00238 rvx = -4.5856e-6*cos((double)(kx*x1 + ky*x2 + kz*x3));
00239 rvy = 2.29279e-6*cos((double)(kx*x1 + ky*x2 + kz*x3));
00240 rvz = 2.29279e-6*cos((double)(kx*x1 + ky*x2 + kz*x3));
00241 rbx = 5.48082e-7;
00242 rbx *= cos((double)(kx*(x1-0.5*pGrid->dx1) + ky*x2 + kz*x3));
00243 rbx += (0.1);
00244 rby = 1.0962e-6;
00245 rby *= cos((double)(kx*x1 + ky*(x2-0.5*pGrid->dx2) + kz*x3));
00246 rby += (0.2);
00247 rbz = 0.0;
00248 }
00249 if (ipert == 7) {
00250 #ifdef ISOTHERMAL
00251 double kappa2 = 2.0*(2.0 - qshear)*Omega_0*Omega_0;
00252 double aa = (kx*kx + ky*ky)*Iso_csound*Iso_csound + kappa2;
00253 double bb = 2.0*qshear*Omega_0*ky*Iso_csound;
00254 double denom = aa*aa + bb*bb;
00255 double rd_hat = (ky*Iso_csound*bb -2.0*Omega_0*aa)*amp/denom;
00256 double px_hat =-Iso_csound*(ky*Iso_csound*aa +2.0*Omega_0*bb)*amp/denom;
00257 double py_hat = (amp + ky*px_hat + (2.0-qshear)*Omega_0*rd_hat)/kx;
00258 rd = 1.0 + rd_hat*cos((double)(kx*x1 + ky*x2));
00259 rvx = px_hat*sin((double)(kx*x1 + ky*x2))/rd;
00260 rvy = py_hat*sin((double)(kx*x1 + ky*x2))/rd;
00261 #endif
00262 rvz = 0.0;
00263 }
00264
00265
00266
00267
00268 pGrid->U[k][j][i].d = rd;
00269 pGrid->U[k][j][i].M1 = rd*rvx;
00270 pGrid->U[k][j][i].M2 = rd*rvy;
00271 #ifndef FARGO
00272 pGrid->U[k][j][i].M2 -= rd*(qshear*Omega_0*x1);
00273 #endif
00274 pGrid->U[k][j][i].M3 = rd*rvz;
00275 #ifdef ADIABATIC
00276 pGrid->U[k][j][i].E = rp/Gamma_1
00277 + 0.5*(SQR(pGrid->U[k][j][i].M1) + SQR(pGrid->U[k][j][i].M2)
00278 + SQR(pGrid->U[k][j][i].M3))/rd;
00279 #endif
00280
00281
00282
00283
00284
00285
00286
00287
00288 #ifdef MHD
00289 if (ifield == 0) {
00290 pGrid->B1i[k][j][i] = rbx;
00291 pGrid->B2i[k][j][i] = rby;
00292 pGrid->B3i[k][j][i] = rbz;
00293 if (i==ie) pGrid->B1i[k][j][ie+1] = pGrid->B1i[k][j][is];
00294 if (j==je) pGrid->B2i[k][je+1][i] = pGrid->B2i[k][js][i];
00295 if (k==ke) pGrid->B3i[ke+1][j][i] = pGrid->B3i[ks][j][i];
00296 }
00297 if (ifield == 1) {
00298 pGrid->B1i[k][j][i] = 0.0;
00299 pGrid->B2i[k][j][i] = 0.0;
00300 pGrid->B3i[k][j][i] = B0*(sin((double)kx*x1));
00301 if (i==ie) pGrid->B1i[k][j][ie+1] = 0.0;
00302 if (j==je) pGrid->B2i[k][je+1][i] = 0.0;
00303 if (k==ke) pGrid->B3i[ke+1][j][i] = B0*(sin((double)kx*x1));
00304 }
00305 if (ifield == 2) {
00306 pGrid->B1i[k][j][i] = 0.0;
00307 pGrid->B2i[k][j][i] = 0.0;
00308 pGrid->B3i[k][j][i] = B0;
00309 if (i==ie) pGrid->B1i[k][j][ie+1] = 0.0;
00310 if (j==je) pGrid->B2i[k][je+1][i] = 0.0;
00311 if (k==ke) pGrid->B3i[ke+1][j][i] = B0;
00312 }
00313 if (ifield == 3) {
00314 pGrid->B1i[k][j][i] = 0.0;
00315 pGrid->B2i[k][j][i] = B0*(cos((double)kx*x1));
00316 pGrid->B3i[k][j][i] = B0*(sin((double)kx*x1));
00317 if (i==ie) pGrid->B1i[k][j][ie+1] = 0.0;
00318 if (j==je) pGrid->B2i[k][je+1][i] = B0*(cos((double)kx*x1));
00319 if (k==ke) pGrid->B3i[ke+1][j][i] = B0*(sin((double)kx*x1));
00320 }
00321 if (ifield == 4) {
00322 pGrid->B1i[k][j][i] = 0.0;
00323 pGrid->B2i[k][j][i] = B0/sqrt(2);
00324 pGrid->B3i[k][j][i] = B0/sqrt(2);
00325 if (i==ie) pGrid->B1i[k][j][ie+1] = 0.0;
00326 if (j==je) pGrid->B2i[k][je+1][i] = B0/sqrt(2);
00327 if (k==ke) pGrid->B3i[ke+1][j][i] = B0/sqrt(2);
00328 }
00329 if (ifield == 5) {
00330 pGrid->B1i[k][j][i] = 0.0;
00331 pGrid->B2i[k][j][i] = B0;
00332 pGrid->B3i[k][j][i] = 0.0;
00333 if (i==ie) pGrid->B1i[k][j][ie+1] = 0.0;
00334 if (j==je) pGrid->B2i[k][je+1][i] = B0;
00335 if (k==ke) pGrid->B3i[ke+1][j][i] = 0.0;
00336 }
00337 #endif
00338 }
00339 }}
00340 #ifdef MHD
00341 for (k=ks; k<=ke; k++) {
00342 for (j=js; j<=je; j++) {
00343 for (i=is; i<=ie; i++) {
00344 pGrid->U[k][j][i].B1c = 0.5*(pGrid->B1i[k][j][i]+pGrid->B1i[k][j][i+1]);
00345 pGrid->U[k][j][i].B2c = 0.5*(pGrid->B2i[k][j][i]+pGrid->B2i[k][j+1][i]);
00346 pGrid->U[k][j][i].B3c = 0.5*(pGrid->B3i[k][j][i]+pGrid->B3i[k+1][j][i]);
00347 #ifdef ADIABATIC
00348 pGrid->U[k][j][i].E += 0.5*(SQR(pGrid->U[k][j][i].B1c)
00349 + SQR(pGrid->U[k][j][i].B2c) + SQR(pGrid->U[k][j][i].B3c));
00350 #endif
00351 }
00352 }
00353 }
00354 #endif
00355
00356
00357
00358 ShearingBoxPot = UnstratifiedDisk;
00359
00360
00361
00362 if (frst == 1) {
00363 dump_history_enroll(hst_rho_Vx_dVy, "<rho Vx dVy>");
00364 dump_history_enroll(hst_rho_dVy2, "<rho dVy^2>");
00365 #ifdef ADIABATIC
00366 dump_history_enroll(hst_E_total, "<E + rho Phi>");
00367 #endif
00368 #ifdef MHD
00369 dump_history_enroll(hst_Bx, "<Bx>");
00370 dump_history_enroll(hst_By, "<By>");
00371 dump_history_enroll(hst_Bz, "<Bz>");
00372 dump_history_enroll(hst_BxBy, "<-Bx By>");
00373 if (ipert == 5) dump_history_enroll(hst_dEw2, "<dEw2>");
00374 if (ipert == 6) dump_history_enroll(hst_dBy, "<dBy>");
00375 #endif
00376 frst = 0;
00377 }
00378
00379
00380 #ifdef RESISTIVITY
00381 eta_Ohm = par_getd_def("problem","eta_O",0.0);
00382 Q_Hall = par_getd_def("problem","Q_H",0.0);
00383 Q_AD = par_getd_def("problem","Q_A",0.0);
00384 d_ind = par_getd_def("problem","d_ind",0.0);
00385 #endif
00386 #ifdef VISCOSITY
00387 nu_iso = par_getd_def("problem","nu_iso",0.0);
00388 nu_aniso = par_getd_def("problem","nu_aniso",0.0);
00389 #endif
00390
00391 return;
00392 }
00393
00394
00395
00396
00397
00398
00399
00400
00401
00402
00403
00404
00405 void problem_write_restart(MeshS *pM, FILE *fp)
00406 {
00407 return;
00408 }
00409
00410
00411
00412
00413
00414 void problem_read_restart(MeshS *pM, FILE *fp)
00415 {
00416
00417
00418 Omega_0 = par_getd_def("problem","Omega",1.0e-3);
00419 qshear = par_getd_def("problem","qshear",1.5);
00420
00421 #ifdef RESISTIVITY
00422 eta_Ohm = par_getd_def("problem","eta_O",0.0);
00423 Q_Hall = par_getd_def("problem","Q_H",0.0);
00424 Q_AD = par_getd_def("problem","Q_A",0.0);
00425 d_ind = par_getd_def("problem","d_ind",0.0);
00426 #endif
00427
00428 #ifdef VISCOSITY
00429 nu_iso = par_getd_def("problem","nu_iso",0.0);
00430 nu_aniso = par_getd_def("problem","nu_aniso",0.0);
00431 #endif
00432
00433
00434
00435 ShearingBoxPot = UnstratifiedDisk;
00436
00437
00438
00439 dump_history_enroll(hst_rho_Vx_dVy, "<rho Vx dVy>");
00440 dump_history_enroll(hst_rho_dVy2, "<rho dVy^2>");
00441 #ifdef ADIABATIC
00442 dump_history_enroll(hst_E_total, "<E + rho Phi>");
00443 #endif
00444 #ifdef MHD
00445 dump_history_enroll(hst_Bx, "<Bx>");
00446 dump_history_enroll(hst_By, "<By>");
00447 dump_history_enroll(hst_Bz, "<Bz>");
00448 dump_history_enroll(hst_BxBy, "<-Bx By>");
00449 #endif
00450
00451 return;
00452 }
00453
00454
00455 ConsFun_t get_usr_expr(const char *expr)
00456 {
00457 if (strcmp(expr,"dVy")==0) return expr_dV2;
00458 if (strcmp(expr,"Jsq")==0) return expr_Jsq;
00459 return NULL;
00460 }
00461
00462 VOutFun_t get_usr_out_fun(const char *name){
00463 return NULL;
00464 }
00465
00466 #ifdef RESISTIVITY
00467 void get_eta_user(GridS *pG, int i, int j, int k,
00468 Real *eta_O, Real *eta_H, Real *eta_A)
00469 {
00470 *eta_O = 0.0;
00471 *eta_H = 0.0;
00472 *eta_A = 0.0;
00473
00474 return;
00475 }
00476 #endif
00477
00478 void Userwork_in_loop(MeshS *pM)
00479 {
00480 }
00481
00482 void Userwork_after_loop(MeshS *pM)
00483 {
00484 }
00485
00486
00487
00488
00489
00490
00491 #define IM1 2147483563
00492 #define IM2 2147483399
00493 #define AM (1.0/IM1)
00494 #define IMM1 (IM1-1)
00495 #define IA1 40014
00496 #define IA2 40692
00497 #define IQ1 53668
00498 #define IQ2 52774
00499 #define IR1 12211
00500 #define IR2 3791
00501 #define NTAB 32
00502 #define NDIV (1+IMM1/NTAB)
00503 #define RNMX (1.0-DBL_EPSILON)
00504
00505
00506
00507
00508
00509
00510
00511
00512
00513
00514 double ran2(long int *idum)
00515 {
00516 int j;
00517 long int k;
00518 static long int idum2=123456789;
00519 static long int iy=0;
00520 static long int iv[NTAB];
00521 double temp;
00522
00523 if (*idum <= 0) {
00524 if (-(*idum) < 1) *idum=1;
00525 else *idum = -(*idum);
00526 idum2=(*idum);
00527 for (j=NTAB+7;j>=0;j--) {
00528 k=(*idum)/IQ1;
00529 *idum=IA1*(*idum-k*IQ1)-k*IR1;
00530 if (*idum < 0) *idum += IM1;
00531 if (j < NTAB) iv[j] = *idum;
00532 }
00533 iy=iv[0];
00534 }
00535 k=(*idum)/IQ1;
00536 *idum=IA1*(*idum-k*IQ1)-k*IR1;
00537 if (*idum < 0) *idum += IM1;
00538 k=idum2/IQ2;
00539 idum2=IA2*(idum2-k*IQ2)-k*IR2;
00540 if (idum2 < 0) idum2 += IM2;
00541 j=(int)(iy/NDIV);
00542 iy=iv[j]-idum2;
00543 iv[j] = *idum;
00544 if (iy < 1) iy += IMM1;
00545 if ((temp=AM*iy) > RNMX) return RNMX;
00546 else return temp;
00547 }
00548
00549 #undef IM1
00550 #undef IM2
00551 #undef AM
00552 #undef IMM1
00553 #undef IA1
00554 #undef IA2
00555 #undef IQ1
00556 #undef IQ2
00557 #undef IR1
00558 #undef IR2
00559 #undef NTAB
00560 #undef NDIV
00561 #undef RNMX
00562
00563
00564
00565
00566
00567 static Real UnstratifiedDisk(const Real x1, const Real x2, const Real x3)
00568 {
00569 Real phi=0.0;
00570 #ifndef FARGO
00571 phi -= qshear*Omega_0*Omega_0*x1*x1;
00572 #endif
00573 return phi;
00574 }
00575
00576
00577
00578
00579
00580
00581 static Real expr_dV2(const GridS *pG, const int i, const int j, const int k)
00582 {
00583 Real x1,x2,x3;
00584 cc_pos(pG,i,j,k,&x1,&x2,&x3);
00585 #ifdef FARGO
00586 return (pG->U[k][j][i].M2/pG->U[k][j][i].d);
00587 #else
00588 return (pG->U[k][j][i].M2/pG->U[k][j][i].d + qshear*Omega_0*x1);
00589 #endif
00590 }
00591
00592
00593
00594
00595
00596
00597 static Real expr_Jsq(const GridS *pG, const int i, const int j, const int k)
00598 {
00599 Real J1,J2,J3;
00600
00601 J1 = (pG->B3i[k][j][i] - pG->B3i[k ][j-1][i ])/pG->dx2 -
00602 (pG->B2i[k][j][i] - pG->B2i[k-1][j ][i ])/pG->dx3;
00603 J2 = (pG->B1i[k][j][i] - pG->B1i[k-1][j ][i ])/pG->dx3 -
00604 (pG->B3i[k][j][i] - pG->B3i[k ][j ][i-1])/pG->dx1;
00605 J3 = (pG->B2i[k][j][i] - pG->B2i[k ][j ][i-1])/pG->dx1 -
00606 (pG->B1i[k][j][i] - pG->B1i[k ][j-1][i ])/pG->dx2;
00607
00608 return SQR(J1)+SQR(J2)+SQR(J3);
00609 }
00610
00611
00612
00613
00614
00615
00616
00617
00618
00619
00620
00621 static Real hst_rho_Vx_dVy(const GridS *pG,const int i,const int j, const int k)
00622 {
00623 Real x1,x2,x3;
00624 cc_pos(pG,i,j,k,&x1,&x2,&x3);
00625 #ifdef FARGO
00626 return pG->U[k][j][i].M1*(pG->U[k][j][i].M2/pG->U[k][j][i].d);
00627 #else
00628 return pG->U[k][j][i].M1*
00629 (pG->U[k][j][i].M2/pG->U[k][j][i].d + qshear*Omega_0*x1);
00630 #endif
00631 }
00632
00633
00634
00635
00636 static Real hst_rho_dVy2(const GridS *pG, const int i, const int j, const int k)
00637 {
00638 Real x1,x2,x3,dVy;
00639 cc_pos(pG,i,j,k,&x1,&x2,&x3);
00640 #ifdef FARGO
00641 dVy = (pG->U[k][j][i].M2/pG->U[k][j][i].d);
00642 #else
00643 dVy = (pG->U[k][j][i].M2/pG->U[k][j][i].d + qshear*Omega_0*x1);
00644 #endif
00645 return pG->U[k][j][i].d*dVy*dVy;
00646 }
00647
00648 #ifdef ADIABATIC
00649
00650
00651
00652 static Real hst_E_total(const GridS *pG, const int i, const int j, const int k)
00653 {
00654 Real x1,x2,x3,phi;
00655 cc_pos(pG,i,j,k,&x1,&x2,&x3);
00656 phi = UnstratifiedDisk(x1, x2, x3);
00657
00658 return pG->U[k][j][i].E + pG->U[k][j][i].d*phi;
00659 }
00660 #endif
00661
00662
00663
00664
00665
00666
00667 #ifdef MHD
00668
00669
00670 static Real hst_Bx(const GridS *pG, const int i, const int j, const int k)
00671 {
00672 return pG->U[k][j][i].B1c;
00673 }
00674
00675
00676
00677 static Real hst_By(const GridS *pG, const int i, const int j, const int k)
00678 {
00679 return pG->U[k][j][i].B2c;
00680 }
00681
00682
00683
00684 static Real hst_Bz(const GridS *pG, const int i, const int j, const int k)
00685 {
00686 return pG->U[k][j][i].B3c;
00687 }
00688
00689
00690
00691
00692 static Real hst_BxBy(const GridS *pG, const int i, const int j, const int k)
00693 {
00694 return -pG->U[k][j][i].B1c*pG->U[k][j][i].B2c;
00695 }
00696
00697 static Real hst_dEw2(const GridS *pG, const int i, const int j, const int k)
00698 {
00699 Real x1,x2,x3,dVx,dVy,dVz,dBz;
00700 cc_pos(pG,i,j,k,&x1,&x2,&x3);
00701 dBz = pG->U[k][j][i].B3c-(sqrt(15.0/16.0))/(2.0*PI)/sqrt(4.*PI);
00702 dVx = pG->U[k][j][i].M1/pG->U[k][j][i].d;
00703 dVy = pG->U[k][j][i].M2/pG->U[k][j][i].d + qshear*Omega_0*x1;
00704 dVz = pG->U[k][j][i].M3/pG->U[k][j][i].d;
00705
00706
00707
00708 return (pG->U[k][j][i].B1c*pG->U[k][j][i].B1c
00709 + pG->U[k][j][i].B2c*pG->U[k][j][i].B2c + dBz*dBz);
00710 }
00711
00712 static Real hst_dBy(const GridS *pG, const int i, const int j, const int k)
00713 {
00714 double fkx, fky, fkz;
00715 double dBy;
00716 Real x1,x2,x3;
00717
00718
00719
00720 fky = 2.0*PI/Ly;
00721 fkx = -4.0*PI/Lx + qshear*Omega_0*fky*pG->time;
00722 fkz = 2.0*PI/Lz;
00723
00724
00725 cc_pos(pG,i,j,k,&x1,&x2,&x3);
00726 dBy = 2.0*(pG->U[k][j][i].B2c - (0.2-0.15*Omega_0*pG->time));
00727 dBy *= cos(fkx*x1 + fky*x2 + fkz*x3);
00728
00729 return dBy;
00730 }
00731
00732 #endif
00733