00001
00002 #include "copyright.h"
00003
00004
00005
00006
00007
00008
00009
00010 #include <math.h>
00011 #include <stdio.h>
00012 #include <stdlib.h>
00013 #include <string.h>
00014 #include "defs.h"
00015 #include "athena.h"
00016 #include "globals.h"
00017 #include "prototypes.h"
00018
00019 #define SEED 31337
00020
00021 static Real rho0,Amp,Beta, R0, rhomin, Field, Hbc, Mbc, Mc;
00022 Real Mollifier(Real Center, Real Scl, Real x);
00023 void ScaleToBeta(GridS *pG, Real beta);
00024 void disk_ir(GridS *pG);
00025 void disk_or(GridS *pG);
00026 static Real ChiMag(Real R);
00027 static Real ChiSub(Real R);
00028 static Real Mrp(const GridS *pG, const int i, const int j, const int k);
00029 static Real Trp(const GridS *pG, const int i, const int j, const int k);
00030 static Real Pb(const GridS *pG, const int i, const int j, const int k);
00031 static Real Vaz(const GridS *pG, const int i, const int j, const int k);
00032 static Real Br(const GridS *pG, const int i, const int j, const int k);
00033 static Real Bp(const GridS *pG, const int i, const int j, const int k);
00034 static Real Bz(const GridS *pG, const int i, const int j, const int k);
00035 static Real Vr(const GridS *pG, const int i, const int j, const int k);
00036 static Real Vp(const GridS *pG, const int i, const int j, const int k);
00037 static Real Vz(const GridS *pG, const int i, const int j, const int k);
00038 static Real MdotR1(const GridS *pG, const int i, const int j, const int k);
00039 static Real MdotR2(const GridS *pG, const int i, const int j, const int k);
00040 static Real MdotR3(const GridS *pG, const int i, const int j, const int k);
00041 static Real MdotR4(const GridS *pG, const int i, const int j, const int k);
00042 static Real Msub(const GridS *pG, const int i, const int j, const int k);
00043 static Real Mrpsub(const GridS *pG, const int i, const int j, const int k);
00044 static Real Bzsub(const GridS *pG, const int i, const int j, const int k);
00045 static Real Bpsub(const GridS *pG, const int i, const int j, const int k);
00046 static Real Pbsub(const GridS *pG, const int i, const int j, const int k);
00047 static void dump_vtksub(MeshS *pM, OutputS *pOut);
00048 void out_ktab(MeshS *pM, OutputS *pOut);
00049
00050 static Real grav_pot(const Real x1, const Real x2, const Real x3) {
00051
00052 return (-1.0/x1);
00053 }
00054
00055 static Real grav_acc(const Real x1, const Real x2, const Real x3) {
00056 Real g;
00057 g = (1.0/SQR(x1));
00058
00059 return g;
00060
00061 }
00062
00063 static Real Omega(const Real R) {
00064 Real Arg;
00065 Arg = pow(R,1.5);
00066 return (1.0/Arg);
00067 }
00068 static Real Shear(const Real R) {
00069 return 1.5;
00070 }
00071
00072 static Real vphi(const Real x1, const Real x2, const Real x3) {
00073 return x1*Omega(x1);
00074 }
00075
00076
00077 static Real BzZero(const Real R, const Real phi, const Real z) {
00078 Real Arg, bz, H0, n, Rs, I;
00079
00080 H0 = sqrt(2.0)*Iso_csound/Omega(R0);
00081 Arg = 2.0*PI*(R-R0)/H0;
00082 n = floor( (R-1.2)/H0 );
00083 Rs = 1.3 + n*H0;
00084 I = ChiMag(R);
00085 bz = I*sin(Arg)*(Rs/R)*Omega(Rs);
00086 return bz;
00087 }
00088
00089 static Real BzNet(const Real R, const Real phi, const Real z) {
00090 Real bz, I;
00091
00092 I = ChiMag(R);
00093 bz = I*Omega(R);
00094 return bz;
00095 }
00096
00097 static Real BpNet(const Real R, const Real phi, const Real z) {
00098 Real Bp, I;
00099 I = ChiMag(R);
00100 Bp = sqrt(rho0/R)*I/Mc;
00101
00102 return Bp;
00103 }
00104
00105 static Real ChiMag(Real R) {
00106 Real I;
00107 if ((R >= 1.2) && (R <= 3.8)) {
00108 I = 1.0;
00109 } else {
00110 I = 0.0;
00111 }
00112 return I;
00113 }
00114 static Real ChiSub(Real R) {
00115 Real I;
00116 if ((R >= 1.6) && (R <= 2.4)) {
00117 I = 1.0;
00118 } else {
00119 I = 0.0;
00120 }
00121 return I;
00122 }
00123
00124
00125
00126
00127 void problem(DomainS *pDomain)
00128 {
00129 int i,j,k;
00130 int is,ie,il,iu,js,je,jl,ju,ks,ke,kl,ku;
00131 int nx1,nx2,nx3,myid=0;
00132 Real x1,x2,x3,y1, r;
00133 GridS *pG = pDomain->Grid;
00134
00135 is = pG->is; ie = pG->ie; nx1 = ie-is+1;
00136 js = pG->js; je = pG->je; nx2 = je-js+1;
00137 ks = pG->ks; ke = pG->ke; nx3 = ke-ks+1;
00138
00139 il = is-nghost*(nx1>1); iu = ie+nghost*(nx1>1); nx1 = iu-il+1;
00140 jl = js-nghost*(nx2>1); ju = je+nghost*(nx2>1); nx2 = ju-jl+1;
00141 kl = ks-nghost*(nx3>1); ku = ke+nghost*(nx3>1); nx3 = ku-kl+1;
00142
00143 rho0 = par_getd_def("problem", "rho0", 100.0);
00144 Amp = par_getd_def("problem", "Amp", 1.0e-2);
00145 Beta = par_getd_def("problem", "Beta",100.0);
00146 R0 = par_getd_def("problem", "R0",2.0);
00147 rhomin = par_getd_def("problem","rhomin",1.0e-3);
00148 Field = par_getd_def("problem","Field",0);
00149 Hbc = par_getd_def("problem","Hbc",1);
00150 Mbc = par_getd_def("problem","Mbc",1);
00151 Mc = par_getd_def("problem","Mc",20.0);
00152
00153 srand(SEED + myID_Comm_world);
00154 for (k=kl; k<=ku; k++) {
00155 for (j=jl; j<=ju; j++) {
00156 for (i=il; i<=iu; i++) {
00157 cc_pos(pG,i,j,k,&x1,&x2,&x3);
00158 x1 = x1vc(pG,i);
00159 r = 2.0*rand()/((double)RAND_MAX) - 1.0;
00160 pG->U[k][j][i].d = rho0;
00161 #ifdef FARGO
00162 pG->U[k][j][i].M2 = 0.0;
00163 #else
00164 pG->U[k][j][i].M2 = pG->U[k][j][i].d*avg1d(vphi,pG,i,j,k);
00165 #endif
00166 pG->U[k][j][i].M2 += pG->U[k][j][i].d*Amp*r*Iso_csound*ChiMag(x1);
00167 r = 2.0*rand()/((double)RAND_MAX)-1.0;
00168 pG->U[k][j][i].M1 = 0.0;
00169 pG->U[k][j][i].M3 = pG->U[k][j][i].d*Amp*r*Iso_csound*ChiMag(x1);
00170 #ifdef MHD
00171 pG->U[k][j][i].B1c = 0.0;
00172
00173 if (Field == 2) {
00174 pG->U[k][j][i].B2c = BpNet(x1,x2,x3);
00175 } else {
00176 pG->U[k][j][i].B2c = 0.0;
00177 }
00178 if (Field == 0) {
00179 pG->U[k][j][i].B3c = BzZero(x1,x2,x3);
00180 } else if (Field == 1) {
00181 pG->U[k][j][i].B3c = BzNet(x1,x2,x3);
00182 } else {
00183 pG->U[k][j][i].B3c = 0.0;
00184 }
00185
00186 pG->B1i[k][j][i] = 0.0;
00187 pG->B2i[k][j][i] = pG->U[k][j][i].B2c;
00188 pG->B3i[k][j][i] = pG->U[k][j][i].B3c;
00189
00190 #endif
00191
00192 }
00193 }
00194 }
00195 #ifdef MHD
00196 if (Field != 2) {
00197 ScaleToBeta(pG,Beta);
00198 }
00199 #endif
00200
00201 StaticGravPot = grav_pot;
00202 x1GravAcc = grav_acc;
00203 bvals_mhd_fun(pDomain,left_x1,disk_ir);
00204 bvals_mhd_fun(pDomain,right_x1,disk_or);
00205 #ifdef FARGO
00206 OrbitalProfile = Omega;
00207 ShearProfile = Shear;
00208 #endif
00209
00210
00211 #ifdef MHD
00212 dump_history_enroll(Br, "<Br>");
00213 dump_history_enroll(Bp, "<Bp>");
00214 dump_history_enroll(Bz, "<Bz>");
00215 dump_history_enroll(Mrp, "<Mrp>");
00216 dump_history_enroll(Trp, "<Trp>");
00217 dump_history_enroll(MdotR1, "<MdotR1>");
00218 dump_history_enroll(MdotR2, "<MdotR2>");
00219 dump_history_enroll(MdotR3, "<MdotR3>");
00220 dump_history_enroll(MdotR4, "<MdotR4>");
00221 dump_history_enroll(Msub, "<Msub>");
00222 dump_history_enroll(Mrpsub, "<Mrpsub>");
00223 dump_history_enroll(Bpsub, "<Bpsub>");
00224 dump_history_enroll(Bzsub, "<Bzsub>");
00225 dump_history_enroll(Pbsub, "<Pbsub>");
00226
00227 #endif
00228 return;
00229
00230 }
00231
00232
00233
00234
00235
00236
00237
00238
00239
00240 static Real MdotR1(const GridS *pG, const int i, const int j, const int k) {
00241 Real dR=pG->dx1, R0=1.0, R, phi, z, Md;
00242
00243 cc_pos(pG,i,j,k,&R,&phi,&z);
00244 if ((R0 > (R-dR)) && (R0 < (R+dR))) {
00245 Md = -1.0*pG->U[k][j][i].M1/dR;
00246 } else {
00247 Md = 0.0;
00248 }
00249 return Md;
00250
00251 }
00252
00253 static Real MdotR2(const GridS *pG, const int i, const int j, const int k) {
00254 Real dR=pG->dx1, R0=2.0, R, phi, z, Md;
00255
00256 cc_pos(pG,i,j,k,&R,&phi,&z);
00257 if ((R0 > (R-dR)) && (R0 < (R+dR))) {
00258 Md = -1.0*pG->U[k][j][i].M1/dR;
00259 } else {
00260 Md = 0.0;
00261 }
00262 return Md;
00263 }
00264
00265 static Real MdotR3(const GridS *pG, const int i, const int j, const int k) {
00266 Real dR=pG->dx1, R0=3.0, R, phi, z, Md;
00267
00268 cc_pos(pG,i,j,k,&R,&phi,&z);
00269 if ((R0 > (R-dR)) && (R0 < (R+dR))) {
00270 Md = -1.0*pG->U[k][j][i].M1/dR;
00271 } else {
00272 Md = 0.0;
00273 }
00274 return Md;
00275 }
00276 static Real MdotR4(const GridS *pG, const int i, const int j, const int k) {
00277 Real dR=pG->dx1, R0=4.0, R, phi, z, Md;
00278
00279 cc_pos(pG,i,j,k,&R,&phi,&z);
00280 if ((R0 > (R-dR)) && (R0 < (R+dR))) {
00281 Md = -1.0*pG->U[k][j][i].M1/dR;
00282 } else {
00283 Md = 0.0;
00284 }
00285 return Md;
00286 }
00287 static Real Pb(const GridS *pG, const int i, const int j, const int k)
00288 {
00289 #ifdef MHD
00290 return ( 0.5*(SQR(pG->U[k][j][i].B1c) + SQR(pG->U[k][j][i].B2c) + SQR(pG->U[k][j][i].B3c)) );
00291 #else
00292 return 0.0;
00293 #endif
00294 }
00295
00296 static Real Mrp(const GridS *pG, const int i, const int j, const int k)
00297 {
00298 #ifdef MHD
00299 return ( -1.0*pG->U[k][j][i].B1c*pG->U[k][j][i].B2c );
00300 #else
00301 return 0.0;
00302 #endif
00303 }
00304
00305 static Real Trp(const GridS *pG, const int i, const int j, const int k)
00306 {
00307
00308 return ((pG->U[k][j][i].M1*pG->U[k][j][i].M2)/(pG->U[k][j][i].d));
00309 }
00310 static Real Vaz(const GridS *pG, const int i, const int j, const int k)
00311 {
00312 #ifdef MHD
00313 return (fabs(pG->U[k][j][i].B3c)/sqrt(pG->U[k][j][i].d));
00314 #else
00315 return 0.0;
00316 #endif
00317 }
00318 static Real Br(const GridS *pG, const int i, const int j, const int k) {
00319 #ifdef MHD
00320 return pG->U[k][j][i].B1c;
00321 #else
00322 return 0.0;
00323 #endif
00324
00325 }
00326 static Real Bp(const GridS *pG, const int i, const int j, const int k) {
00327 #ifdef MHD
00328 return pG->U[k][j][i].B2c;
00329 #else
00330 return 0.0;
00331 #endif
00332
00333 }
00334 static Real Bz(const GridS *pG, const int i, const int j, const int k) {
00335 #ifdef MHD
00336 return pG->U[k][j][i].B3c;
00337 #else
00338 return 0.0;
00339 #endif
00340
00341 }
00342 static Real Vr(const GridS *pG, const int i, const int j, const int k) {
00343 return (pG->U[k][j][i].M1/pG->U[k][j][i].d);
00344
00345 }
00346 static Real Vp(const GridS *pG, const int i, const int j, const int k) {
00347 return (pG->U[k][j][i].M2/pG->U[k][j][i].d);
00348
00349 }
00350 static Real Vz(const GridS *pG, const int i, const int j, const int k) {
00351 return (pG->U[k][j][i].M3/pG->U[k][j][i].d);
00352
00353 }
00354
00355 static Real Msub(const GridS *pG, const int i, const int j, const int k) {
00356 Real I, R;
00357
00358 R = x1vc(pG,i);
00359 I = ChiSub(R);
00360 return (pG->U[k][j][i].d*I);
00361 }
00362
00363 static Real Mrpsub(const GridS *pG, const int i, const int j, const int k) {
00364 Real I, R;
00365
00366 R = x1vc(pG,i);
00367 I = ChiSub(R);
00368 #ifdef MHD
00369 return (-1.0*pG->U[k][j][i].B1c*pG->U[k][j][i].B2c*I);
00370 #else
00371 return 0.0;
00372 #endif
00373
00374 }
00375
00376 static Real Bpsub(const GridS *pG, const int i, const int j, const int k) {
00377 Real I, R;
00378
00379 R = x1vc(pG,i);
00380 I = ChiSub(R);
00381 #ifdef MHD
00382 return (pG->U[k][j][i].B2c*I);
00383 #else
00384 return 0.0;
00385 #endif
00386
00387 }
00388
00389 static Real Bzsub(const GridS *pG, const int i, const int j, const int k) {
00390 Real I, R;
00391
00392 R = x1vc(pG,i);
00393 I = ChiSub(R);
00394 #ifdef MHD
00395 return (pG->U[k][j][i].B3c*I);
00396 #else
00397 return 0.0;
00398 #endif
00399
00400 }
00401
00402 static Real Pbsub(const GridS *pG, const int i, const int j, const int k) {
00403 Real I, R;
00404
00405 R = x1vc(pG,i);
00406 I = ChiSub(R);
00407 #ifdef MHD
00408 return (0.5*(SQR(pG->U[k][j][i].B1c) + SQR(pG->U[k][j][i].B2c) + SQR(pG->U[k][j][i].B3c))*I);
00409 #else
00410 return 0.0;
00411 #endif
00412
00413 }
00414 void problem_write_restart(MeshS *pM, FILE *fp)
00415 {
00416 return;
00417 }
00418
00419 void problem_read_restart(MeshS *pM, FILE *fp)
00420 {
00421
00422 R0 = par_getd_def("problem", "R0",2.0);
00423 Hbc = par_getd_def("problem","Hbc",1);
00424 Mbc = par_getd_def("problem","Mbc",1);
00425
00426 StaticGravPot = grav_pot;
00427 x1GravAcc = grav_acc;
00428 bvals_mhd_fun(&(pM->Domain[0][0]),left_x1,disk_ir);
00429 bvals_mhd_fun(&(pM->Domain[0][0]),right_x1,disk_or);
00430 #ifdef FARGO
00431 OrbitalProfile = Omega;
00432 ShearProfile = Shear;
00433 #endif
00434
00435
00436 #ifdef MHD
00437 dump_history_enroll(Br, "<Br>");
00438 dump_history_enroll(Bp, "<Bp>");
00439 dump_history_enroll(Bz, "<Bz>");
00440 dump_history_enroll(Mrp, "<Mrp>");
00441 dump_history_enroll(Trp, "<Trp>");
00442 dump_history_enroll(MdotR1, "<MdotR1>");
00443 dump_history_enroll(MdotR2, "<MdotR2>");
00444 dump_history_enroll(MdotR3, "<MdotR3>");
00445 dump_history_enroll(MdotR4, "<MdotR4>");
00446 dump_history_enroll(Msub, "<Msub>");
00447 dump_history_enroll(Mrpsub, "<Mrpsub>");
00448 dump_history_enroll(Bpsub, "<Bpsub>");
00449 dump_history_enroll(Bzsub, "<Bzsub>");
00450 dump_history_enroll(Pbsub, "<Pbsub>");
00451
00452
00453 #endif
00454
00455 return;
00456 }
00457
00458 ConsFun_t get_usr_expr(const char *expr)
00459 {
00460 if (strcmp(expr,"Pb")==0) return Pb;
00461 else if(strcmp(expr,"Mrp")==0) return Mrp;
00462 else if(strcmp(expr,"Trp")==0) return Trp;
00463 else if(strcmp(expr,"Vaz")==0) return Vaz;
00464 else if(strcmp(expr,"Br")==0) return Br;
00465 else if(strcmp(expr,"Bp")==0) return Bp;
00466 else if(strcmp(expr,"Bz")==0) return Bz;
00467 else if(strcmp(expr,"Vr")==0) return Vr;
00468 else if(strcmp(expr,"Vp")==0) return Vp;
00469 else if(strcmp(expr,"Vz")==0) return Vz;
00470 else if(strcmp(expr,"Msub")==0) return Msub;
00471 else if(strcmp(expr,"Mrpsub")==0) return Mrpsub;
00472 else if(strcmp(expr,"Bpsub")==0) return Bpsub;
00473 else if(strcmp(expr,"Bzsub")==0) return Bzsub;
00474 else if(strcmp(expr,"Pbsub")==0) return Pbsub;
00475
00476 return NULL;
00477 }
00478
00479 VOutFun_t get_usr_out_fun(const char *name){
00480 if(strcmp(name,"vtksub")==0) {
00481 return dump_vtksub;
00482 } else if(strcmp(name,"ktab")==0) {
00483 return out_ktab;
00484 }
00485 return NULL;
00486 }
00487
00488 void Userwork_in_loop(MeshS *pM)
00489 {
00490 }
00491
00492 void Userwork_after_loop(MeshS *pM)
00493 {
00494 }
00495
00496
00497 Real Mollifier(Real Center, Real Scl, Real z) {
00498 Real Arg, a, M;
00499
00500 a = fabs( (z-Center)/Scl );
00501
00502 if (a + TINY_NUMBER < 1.0) {
00503 Arg = 1.0 - SQR(a);
00504 M = exp(-1.0/Arg);
00505 } else {
00506 M = 0.0;
00507 }
00508 return M;
00509 }
00510
00511 void ScaleToBeta(GridS *pG, Real beta) {
00512 #ifdef MHD
00513 int i,j,k;
00514 int is,ie,js,je,ks,ke,nx1,nx2,nx3;
00515 int il,iu,jl,ju,kl,ku;
00516 Real Pgas = 0.0, Pb = 0.0, TotPgas = 0.0, TotPb = 0.0, scl = 0.0, CellVol;
00517
00518 is = pG->is; ie = pG->ie;
00519 js = pG->js; je = pG->je;
00520 ks = pG->ks; ke = pG->ke;
00521 nx1 = (ie-is)+1 + 2*nghost;
00522 nx2 = (je-js)+1 + 2*nghost;
00523 nx3 = (ke-ks)+1 + 2*nghost;
00524 il = is - nghost*(ie > is);
00525 jl = js - nghost*(je > js);
00526 kl = ks - nghost*(ke > ks);
00527 iu = ie + nghost*(ie > is);
00528 ju = je + nghost*(je > js);
00529 ku = ke + nghost*(ke > ks);
00530
00531
00532 for (k=ks; k<=ke; k++) {
00533 for (j=js; j<=je; j++) {
00534 for (i=is; i<=ie; i++) {
00535 CellVol = (pG->r[i])*(pG->dx1)*(pG->dx2)*(pG->dx3);
00536 #ifndef ISOTHERMAL
00537 Pgas += (Gamma-1)*pG->U[k][j][i].E*CellVol;
00538 #else
00539 Pgas += Iso_csound2*pG->U[k][j][i].d*CellVol;
00540 #endif
00541 Pb += 0.5*(SQR(pG->U[k][j][i].B1c) + SQR(pG->U[k][j][i].B2c)
00542 + SQR(pG->U[k][j][i].B3c))*CellVol;
00543 }
00544 }
00545 }
00546
00547 #ifdef MPI_PARALLEL
00548 MPI_Reduce(&Pgas, &TotPgas, 1, MPI_DOUBLE, MPI_SUM, 0, MPI_COMM_WORLD);
00549 MPI_Reduce(&Pb, &TotPb, 1, MPI_DOUBLE, MPI_SUM, 0, MPI_COMM_WORLD);
00550 MPI_Bcast(&TotPgas, 1, MPI_DOUBLE, 0, MPI_COMM_WORLD);
00551 MPI_Bcast(&TotPb, 1, MPI_DOUBLE, 0, MPI_COMM_WORLD);
00552 #else
00553 TotPgas = Pgas;
00554 TotPb = Pb;
00555 #endif //PARALLEL
00556
00557 scl = sqrt(TotPgas/(TotPb*beta));
00558 for (k=kl; k<=ku; k++) {
00559 for (j=jl; j<=ju; j++) {
00560 for (i=il; i<=iu; i++) {
00561 pG->U[k][j][i].B1c *= scl;
00562 pG->U[k][j][i].B2c *= scl;
00563 pG->U[k][j][i].B3c *= scl;
00564 pG->B1i[k][j][i] *= scl;
00565 pG->B2i[k][j][i] *= scl;
00566 pG->B3i[k][j][i] *= scl;
00567 }
00568 }
00569 }
00570 #endif
00571 }
00572
00573 void disk_ir(GridS *pGrid) {
00574 int is = pGrid->is;
00575 int js = pGrid->js, je = pGrid->je;
00576 int ks = pGrid->ks, ke = pGrid->ke;
00577 int i,j,k;
00578 #ifdef MHD
00579 int ju, ku;
00580 #endif
00581 Real Vkep,R,p,z, RBr, Lper;
00582
00583 for (k=ks; k<=ke; k++) {
00584 for (j=js; j<=je; j++) {
00585 cc_pos(pGrid,is,j,k,&R,&p,&z);
00586 #ifdef MHD
00587 RBr = (R-0.5*pGrid->dx1)*pGrid->B1i[k][j][is];
00588 #endif
00589
00590 #ifdef FARGO
00591 Lper = R*pGrid->U[k][j][is].M2;
00592 #else
00593 Vkep = sqrt(R*(*x1GravAcc)(R,p,z));
00594 Lper = R*pGrid->U[k][j][is].M2 - R*pGrid->U[k][j][is].d*Vkep;
00595 #endif
00596 Lper = MIN(Lper,0.0);
00597
00598 for (i=1; i<=nghost; i++) {
00599 pGrid->U[k][j][is-i] = pGrid->U[k][j][is];
00600
00601 pGrid->U[k][j][is-i].M1 = MIN(pGrid->U[k][j][is-i].M1,0.0);
00602
00603 cc_pos(pGrid,is-i,j,k,&R,&p,&z);
00604 Vkep = sqrt(R*(*x1GravAcc)(R,p,z));
00605 #ifdef FARGO
00606 if (Hbc == 1) {
00607 pGrid->U[k][j][is-i].M2 = 0.0;
00608 } else {
00609 pGrid->U[k][j][is-i].M2 = Lper/R;
00610 }
00611 #else
00612 if (Hbc == 1) {
00613 pGrid->U[k][j][is-i].M2 = pGrid->U[k][j][is].d*Vkep;
00614 } else {
00615 pGrid->U[k][j][is-i].M2 = Lper/R + pGrid->U[k][j][is-i].d*Vkep;
00616 }
00617 #endif
00618 #ifdef MHD
00619 if (Mbc == 0) {
00620 pGrid->U[k][j][is-i].B1c = RBr/R;
00621 pGrid->U[k][j][is-i].B2c = 0.0;
00622 pGrid->U[k][j][is-i].B3c = 0.0;
00623 }
00624 #endif
00625 }
00626 }
00627 }
00628
00629 #ifdef MHD
00630
00631 for (k=ks; k<=ke; k++) {
00632 for (j=js; j<=je; j++) {
00633 cc_pos(pGrid,is,j,k,&R,&p,&z);
00634 RBr = (R-0.5*pGrid->dx1)*pGrid->B1i[k][j][is];
00635 for (i=1; i<=nghost-1; i++) {
00636 cc_pos(pGrid,is,j,k,&R,&p,&z);
00637 if (Mbc == 0) {
00638 pGrid->B1i[k][j][is-i] = RBr/(R-0.5*pGrid->dx1);
00639 } else {
00640 pGrid->B1i[k][j][is-i] = pGrid->B1i[k][j][is];
00641 }
00642 }
00643 }
00644 }
00645
00646 if (pGrid->Nx[1] > 1) ju=je+1; else ju=je;
00647 for (k=ks; k<=ke; k++) {
00648 for (j=js; j<=ju; j++) {
00649 for (i=1; i<=nghost; i++) {
00650 if (Mbc == 1) {
00651 pGrid->B2i[k][j][is-i] = pGrid->B2i[k][j][is];
00652 } else {
00653 pGrid->B2i[k][j][is-i] = 0.0;
00654 }
00655 }
00656 }
00657 }
00658
00659 if (pGrid->Nx[2] > 1) ku=ke+1; else ku=ke;
00660 for (k=ks; k<=ku; k++) {
00661 for (j=js; j<=je; j++) {
00662 for (i=1; i<=nghost; i++) {
00663 if (Mbc == 1) {
00664 pGrid->B3i[k][j][is-i] = pGrid->B3i[k][j][is];
00665 } else {
00666 pGrid->B3i[k][j][is-i] = 0.0;
00667 }
00668
00669 }
00670 }
00671 }
00672 #endif
00673
00674 return;
00675
00676
00677 }
00678
00679 void disk_or(GridS *pGrid) {
00680
00681 int ie = pGrid->ie;
00682 int js = pGrid->js, je = pGrid->je;
00683 int ks = pGrid->ks, ke = pGrid->ke;
00684 int i,j,k;
00685 #ifdef MHD
00686 int ju, ku;
00687 #endif
00688 Real Vkep, R,p,z, RBr, Lper;
00689
00690 for (k=ks; k<=ke; k++) {
00691 for (j=js; j<=je; j++) {
00692 #ifdef MHD
00693 RBr = pGrid->ri[ie+1]*pGrid->B1i[k][j][ie+1];
00694 #endif
00695
00696 #ifdef FARGO
00697 Lper = pGrid->r[ie]*pGrid->U[k][j][ie].M2;
00698 #else
00699 cc_pos(pGrid,ie,j,k,&R,&p,&z);
00700 Vkep = sqrt(R*(*x1GravAcc)(R,p,z));
00701 Lper = R*pGrid->U[k][j][ie].M2 - x1p*pGrid->U[k][j][ie]*Vkep;
00702 #endif
00703 for (i=1; i<=nghost; i++) {
00704 pGrid->U[k][j][ie+i] = pGrid->U[k][j][ie];
00705
00706 pGrid->U[k][j][ie+i].M1 = MAX(pGrid->U[k][j][ie+i].M1,0.0);
00707 cc_pos(pGrid,ie+i,j,k,&R,&p,&z);
00708 Vkep = sqrt(R*(*x1GravAcc)(R,p,z));
00709 #ifdef FARGO
00710 if (Hbc == 1) {
00711 pGrid->U[k][j][ie+i].M2 = 0.0;
00712 } else {
00713 pGrid->U[k][j][ie+i].M2 = Lper/R;
00714 }
00715 #else
00716 if (Hbc == 1) {
00717 pGrid->U[k][j][ie+i].M2 = pGrid->U[k][j][ie].d*Vkep;
00718 } else {
00719 pGrid->U[k][j][ie+i].M2 = Lper/R + pGrid->U[k][j][ie+i].d*Vkep;
00720 }
00721 #endif
00722 #ifdef MHD
00723 if (Mbc == 0) {
00724 pGrid->U[k][j][ie+i].B1c = RBr/R;
00725 pGrid->U[k][j][ie+i].B2c = 0.0;
00726 pGrid->U[k][j][ie+i].B3c = 0.0;
00727 }
00728 #endif
00729
00730 }
00731 }
00732 }
00733
00734 #ifdef MHD
00735
00736 for (k=ks; k<=ke; k++) {
00737 for (j=js; j<=je; j++) {
00738 RBr = pGrid->ri[ie+1]*pGrid->B1i[k][j][ie+1];
00739 for (i=2; i<=nghost; i++) {
00740 cc_pos(pGrid,ie+i,j,k,&R,&p,&z);
00741 if (Mbc == 0) {
00742 pGrid->B1i[k][j][ie+i] = RBr/R;
00743 } else {
00744 pGrid->B1i[k][j][ie+i] = pGrid->B1i[k][j][ie];
00745 }
00746 }
00747 }
00748 }
00749
00750 if (pGrid->Nx[1] > 1) ju=je+1; else ju=je;
00751 for (k=ks; k<=ke; k++) {
00752 for (j=js; j<=ju; j++) {
00753 for (i=1; i<=nghost; i++) {
00754 if (Mbc == 0) {
00755 pGrid->B2i[k][j][ie+i] = 0.0;
00756 } else {
00757 pGrid->B2i[k][j][ie+i] = pGrid->B2i[k][j][ie];
00758 }
00759 }
00760 }
00761 }
00762
00763 if (pGrid->Nx[2] > 1) ku=ke+1; else ku=ke;
00764 for (k=ks; k<=ku; k++) {
00765 for (j=js; j<=je; j++) {
00766 for (i=1; i<=nghost; i++) {
00767 if (Mbc == 0) {
00768 pGrid->B3i[k][j][ie+i] = 0.0;
00769 } else {
00770 pGrid->B3i[k][j][ie+i] = pGrid->B3i[k][j][ie];
00771 }
00772 }
00773 }
00774 }
00775 #endif
00776
00777 return;
00778
00779 }
00780
00781 static void dump_vtksub(MeshS *pM, OutputS *pOut) {
00782 GridS *pGrid;
00783 FILE *pfile;
00784 char *fname, *plev = NULL, *pdom = NULL;
00785 char levstr[8], domstr[8];
00786
00787 int i, j, k, il, iu, jl, ju, kl, ku;
00788 int big_end = ath_big_endian();
00789 int ndata0;
00790 float *data;
00791 double x1m, x2m, x3m, x1M, x1, x2, x3, dR;
00792
00793
00794
00795 pGrid = pM->Domain[0][0].Grid;
00796
00797 il = pGrid->is, iu = pGrid->ie;
00798 jl = pGrid->js, ju = pGrid->je;
00799 kl = pGrid->ks, ku = pGrid->ke;
00800 x1m = pGrid->MinX[0];
00801 x2m = pGrid->MinX[1];
00802 x3m = pGrid->MinX[2];
00803
00804 x1M = pGrid->MaxX[0];
00805 dR = pGrid->dx1;
00806
00807 if ((x1M <= 1.6) || (x1m >= 2.4)) {
00808 return;
00809 }
00810 il = -1;
00811 for (i = pGrid->is; i <= pGrid->ie; i++) {
00812 cc_pos(pGrid, i, jl, kl, &x1, &x2, &x3);
00813
00814 if ((x1 + 0.5 * dR >= 1.6) && (il < 0)) {
00815 il = i;
00816 x1m = x1 - 0.5 * dR;
00817 }
00818 if (x1 - 0.5 * dR <= 2.4) {
00819 iu = i;
00820 }
00821 }
00822
00823 ndata0 = iu - il + 1;
00824
00825
00826 if ((fname = ath_fname(plev, pM->outfilename, plev, pdom, num_digit,
00827 pOut->num, "sub", "vtk")) == NULL) {
00828 ath_error("[dump_vtk]: Error constructing filename\n");
00829 }
00830
00831 if ((pfile = fopen(fname, "w")) == NULL) {
00832 ath_error("[dump_vtk]: Unable to open vtk dump file\n");
00833 return;
00834 }
00835
00836
00837
00838 if ((data = (float *) malloc(3 * ndata0 * sizeof(float))) == NULL) {
00839 ath_error("[dump_vtk]: malloc failed for temporary array\n");
00840 return;
00841 }
00842
00843
00844
00845
00846 fprintf(pfile, "# vtk DataFile Version 2.0\n");
00847
00848
00849
00850 fprintf(pfile, "Subvolume variables at time= %e, level= %i, domain= %i\n",
00851 pGrid->time, 0, 0);
00852
00853
00854
00855 fprintf(pfile, "BINARY\n");
00856
00857
00858
00859
00860
00861 fprintf(pfile, "DATASET STRUCTURED_POINTS\n");
00862 if (pGrid->Nx[1] == 1) {
00863 fprintf(pfile, "DIMENSIONS %d %d %d\n", iu - il + 2, 1, 1);
00864 } else {
00865 if (pGrid->Nx[2] == 1) {
00866 fprintf(pfile, "DIMENSIONS %d %d %d\n", iu - il + 2, ju - jl + 2, 1);
00867 } else {
00868 fprintf(pfile, "DIMENSIONS %d %d %d\n", iu - il + 2, ju - jl + 2, ku - kl
00869 + 2);
00870 }
00871 }
00872 fprintf(pfile, "ORIGIN %e %e %e \n", x1m, x2m, x3m);
00873 fprintf(pfile, "SPACING %e %e %e \n", pGrid->dx1, pGrid->dx2, pGrid->dx3);
00874
00875
00876
00877 fprintf(pfile, "CELL_DATA %d \n", (iu - il + 1) * (ju - jl + 1) * (ku - kl
00878 + 1));
00879
00880
00881
00882 fprintf(pfile, "SCALARS density float\n");
00883 fprintf(pfile, "LOOKUP_TABLE default\n");
00884 for (k = kl; k <= ku; k++) {
00885 for (j = jl; j <= ju; j++) {
00886 for (i = il; i <= iu; i++) {
00887 data[i - il] = (float) pGrid->U[k][j][i].d;
00888 }
00889 if (!big_end)
00890 ath_bswap(data, sizeof(float), iu - il + 1);
00891 fwrite(data, sizeof(float), (size_t) ndata0, pfile);
00892 }
00893 }
00894
00895
00896
00897 fprintf(pfile, "\nVECTORS momentum float\n");
00898 for (k = kl; k <= ku; k++) {
00899 for (j = jl; j <= ju; j++) {
00900 for (i = il; i <= iu; i++) {
00901 data[3 * (i - il)] = (float) pGrid->U[k][j][i].M1;
00902 data[3 * (i - il) + 1] = (float) pGrid->U[k][j][i].M2;
00903 data[3 * (i - il) + 2] = (float) pGrid->U[k][j][i].M3;
00904 }
00905 if (!big_end)
00906 ath_bswap(data, sizeof(float), 3 * (iu - il + 1));
00907 fwrite(data, sizeof(float), (size_t) (3 * ndata0), pfile);
00908 }
00909 }
00910
00911
00912
00913 #ifdef MHD
00914 fprintf(pfile, "\nVECTORS cell_centered_B float\n");
00915 for (k = kl; k <= ku; k++) {
00916 for (j = jl; j <= ju; j++) {
00917 for (i = il; i <= iu; i++) {
00918 data[3 * (i - il)] = (float) pGrid->U[k][j][i].B1c;
00919 data[3 * (i - il) + 1] = (float) pGrid->U[k][j][i].B2c;
00920 data[3 * (i - il) + 2] = (float) pGrid->U[k][j][i].B3c;
00921 }
00922 if (!big_end)
00923 ath_bswap(data, sizeof(float), 3 * (iu - il + 1));
00924 fwrite(data, sizeof(float), (size_t) (3 * ndata0), pfile);
00925 }
00926 }
00927
00928 #endif
00929
00930
00931
00932 fclose(pfile);
00933 free(data);
00934
00935 return;
00936 }
00937
00938 void out_ktab(MeshS *pM, OutputS *pOut)
00939 {
00940 GridS *pGrid=pM->Domain[0][0].Grid;
00941 int i,nx1;
00942 FILE *pFile;
00943 char fmt[80],*fname,*plev=NULL,*pdom=NULL;
00944 char levstr[8],domstr[8];
00945 Real *data=NULL;
00946 Real dmin, dmax, xworld;
00947
00948
00949 if(pOut->dat_fmt == NULL){
00950 sprintf(fmt," %%12.8e");
00951 }
00952 else{
00953 sprintf(fmt," %s",pOut->dat_fmt);
00954 }
00955
00956
00957 data = OutData1(pGrid,pOut,&nx1);
00958 if (data == NULL) return;
00959
00960 minmax1(data,nx1,&dmin,&dmax);
00961
00962
00963 if((fname = ath_fname(plev,pM->outfilename,0,0,0,0,
00964 pOut->id,"tab")) == NULL){
00965 ath_error("[output_tab]: Error constructing filename\n");
00966 }
00967
00968 pFile = fopen(fname,"a");
00969
00970 if (pFile == NULL) {
00971 ath_error("[output_tab]: Unable to open tab file %s\n",fname);
00972 }
00973
00974
00975
00976 if (pOut->num == 0) {
00977 for (i=0; i<nx1; i++) {
00978 fprintf(pFile,"%12d\t",pGrid->Disp[0]+i);
00979 }
00980 fprintf(pFile,"\n");
00981 }
00982 for (i=0; i<nx1; i++) {
00983 fprintf(pFile,fmt,data[i]);
00984 fprintf(pFile,"\t");
00985 }
00986 fprintf(pFile,"\n");
00987
00988 pOut->gmin = MIN(dmin,pOut->gmin);
00989 pOut->gmax = MAX(dmax,pOut->gmax);
00990
00991 fclose(pFile);
00992 free_1d_array(data);
00993 }