Athena: prob/cylnewtmri.c Source File

Go to the documentation of this file.
00001 
00002 #include "copyright.h"
00003 /*============================================================================*/
00004 /*! \file cylnewtmri.c
00005  * \brief Test of the evolution/saturation of the MRI in an unstratified/uniform
00006  * Newtonian disk
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         // Define Bp so that the critical M is constant in radius
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 /*=========================== PUBLIC FUNCTIONS ===============================*/
00124 /*----------------------------------------------------------------------------*/
00125 /* problem:  */
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 /* MHD */
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         // Enroll history functions
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  * PROBLEM USER FUNCTIONS:
00233  * problem_write_restart() - writes problem-specific user data to restart files
00234  * problem_read_restart()  - reads problem-specific user data from restart files
00235  * get_usr_expr()          - sets pointer to expression for special output data
00236  * get_usr_out_fun()       - returns a user defined output function pointer
00237  * Userwork_in_loop        - problem specific work IN     main loop
00238  * Userwork_after_loop     - problem specific work AFTER  main loop
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         // Note, this assumes Fargo to calculate perturbation Vp
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         // Enroll history functions
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 // Private functions
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   /* Calculate total gas/magnetic pressure on local tile */
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   //calculate and use scaling factor so that the correct beta is ensured
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 /* MHD */
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; /* j-upper, k-upper */
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         // Calculate angular momentum in rotating frame
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                                 // Enforce diode
00601                                 pGrid->U[k][j][is-i].M1 = MIN(pGrid->U[k][j][is-i].M1,0.0);
00602                                 // Calculate Keplerian velocity
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 /* B1i is not set at i=is-nghost */
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 /* MHD */
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; /* j-upper, k-upper */
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         // Calculate angular momentum in rotating frame
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                                 // Enforce diode
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 /* i=ie+1 is not a boundary condition for the interface field B1i */
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 /* MHD */
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         /* Upper and Lower bounds on i,j,k for data dump */
00787         int i, j, k, il, iu, jl, ju, kl, ku;
00788         int big_end = ath_big_endian();
00789         int ndata0;
00790         float *data; /* points to 3*ndata0 allocated floats */
00791         double x1m, x2m, x3m, x1M, x1, x2, x3, dR;
00792 
00793         /* Loop over all Domains in Mesh, and output Grid data */
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         // Exclude tiles not in the subvolume
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                 // Get first index in range
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         /* construct filename, open file */
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         /* Allocate memory for temporary array of floats */
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         /* There are five basic parts to the VTK "legacy" file format.  */
00844         /*  1. Write file version and identifier */
00845 
00846         fprintf(pfile, "# vtk DataFile Version 2.0\n");
00847 
00848         /*  2. Header */
00849 
00850         fprintf(pfile, "Subvolume variables at time= %e, level= %i, domain= %i\n",
00851                         pGrid->time, 0, 0);
00852 
00853         /*  3. File format */
00854 
00855         fprintf(pfile, "BINARY\n");
00856 
00857         /*  4. Dataset structure */
00858 
00859         /* Set the Grid origin */
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         /*  5. Data  */
00876 
00877         fprintf(pfile, "CELL_DATA %d \n", (iu - il + 1) * (ju - jl + 1) * (ku - kl
00878                         + 1));
00879 
00880         /* Write density */
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         /* Write momentum or velocity */
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         /* Write cell centered B */
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         /* close file and free memory */
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 /* Add a white space to the format, setup format for integer zone columns */
00949   if(pOut->dat_fmt == NULL){
00950      sprintf(fmt," %%12.8e"); /* Use a default format */
00951   }
00952   else{
00953     sprintf(fmt," %s",pOut->dat_fmt);
00954   }
00955 
00956 /* compute 1D array of data */
00957   data = OutData1(pGrid,pOut,&nx1);
00958   if (data == NULL) return;  /* slice not in range of Grid */
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 /* open filename */
00970   if (pFile == NULL) {
00971     ath_error("[output_tab]: Unable to open tab file %s\n",fname);
00972   }
00973 
00974 /* write data */
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 /* Compute and store global min/max, for output at end of run */
00988   pOut->gmin = MIN(dmin,pOut->gmin);
00989   pOut->gmax = MAX(dmax,pOut->gmax);
00990 
00991   fclose(pFile);
00992   free_1d_array(data); /* Free the memory we malloc'd */
00993 }