Athena: prob/planet-disk.c Source File

Go to the documentation of this file.
00001 #include "copyright.h"
00002 /*============================================================================*/
00003 /*! \file planet-disk.c
00004  *  \brief Problem generator for planet embedded in a disk, using the
00005  *   shearing sheet approximation.
00006  *
00007  * Code must be configured using --enable-shearing-box */
00008 /*============================================================================*/
00009 #include <float.h>
00010 #include <math.h>
00011 #include <stdlib.h>
00012 #include <string.h>
00013 #include "defs.h"
00014 #include "athena.h"
00015 #include "globals.h"
00016 #include "prototypes.h"
00017 
00018 /*==============================================================================
00019  * PRIVATE FUNCTION PROTOTYPES:
00020  * PlanetPot()   - static gravitational potential of planet
00021  * UnstratifiedDisk() - tidal potential in shearing box
00022  * expr_dV2()    - computes delta(Vy)
00023  * hst_*         - new history variables
00024  *============================================================================*/
00025 
00026 
00027 void constant_iib(GridS *pGrid);
00028 void constant_oib(GridS *pGrid);
00029 
00030 static Real Mp,Rsoft,Xplanet,Yplanet,Zplanet;
00031 static Real ramp_time,insert_time;
00032 static Real PlanetPot(const Real x1, const Real x2, const Real x3);
00033 static Real UnstratifiedDisk(const Real x1, const Real x2, const Real x3);
00034 static Real expr_dV2(const GridS *pG, const int i, const int j, const int k);
00035 static Real hst_rho_Vx_dVy(const GridS *pG,const int i,const int j,const int k);
00036 static Real hst_rho_dVy2(const GridS *pG,const int i, const int j, const int k);
00037 #ifdef ADIABATIC
00038 static Real hst_E_total(const GridS *pG, const int i, const int j, const int k);
00039 #endif
00040 
00041 /*=========================== PUBLIC FUNCTIONS ===============================*/
00042 /* problem:  */
00043 
00044 void problem(DomainS *pDomain)
00045 {
00046   GridS *pGrid = pDomain->Grid;
00047   int is = pGrid->is, ie = pGrid->ie;
00048   int js = pGrid->js, je = pGrid->je;
00049   int ks = pGrid->ks, ke = pGrid->ke;
00050   int i,j,k,BCFlag;
00051   Real x1,x2,x3;
00052   Real den = 1.0, pres = 1.0e-6;
00053   static int frst=1;  /* flag so new history variables enrolled only once */
00054 
00055 #ifdef SHEARING_BOX
00056 /* specify xy (r-phi) plane */
00057   ShBoxCoord = xy;
00058 #endif
00059 
00060 /* Read problem parameters.  Note Omega_0 set to 10^{-3} by default */
00061 #ifdef SHEARING_BOX
00062   Omega_0 = par_getd_def("problem","omega",1.0e-3);
00063   qshear  = par_getd_def("problem","qshear",1.5);
00064 #endif
00065   Mp      = par_getd_def("problem","Mplanet",0.0);
00066   Xplanet = par_getd_def("problem","Xplanet",0.0);
00067   Yplanet = par_getd_def("problem","Yplanet",0.0);
00068   Zplanet = par_getd_def("problem","Zplanet",0.0);
00069   Rsoft   = par_getd_def("problem","Rsoft",0.1);
00070   ramp_time = 0.0;
00071   insert_time = par_getd_def("problem","insert_time",0.0);
00072 
00073 /* Compute field strength based on beta.  */
00074 #ifdef ISOTHERMAL
00075   pres = Iso_csound2;
00076 #endif
00077 
00078   for (k=ks; k<=ke; k++) {
00079   for (j=js; j<=je; j++) {
00080     for (i=is; i<=ie; i++) {
00081       cc_pos(pGrid,i,j,k,&x1,&x2,&x3);
00082 
00083 /* Initialize d, M, and P.  With FARGO do not initialize the background shear */
00084 
00085       pGrid->U[k][j][i].d  = den;
00086       pGrid->U[k][j][i].M1 = 0.0;
00087       pGrid->U[k][j][i].M2 = 0.0;
00088 #ifdef SHEARING_BOX
00089 #ifndef FARGO
00090       pGrid->U[k][j][i].M2 -= den*(qshear*Omega_0*x1);
00091 #endif
00092 #endif
00093       pGrid->U[k][j][i].M3 = 0.0;
00094 #ifdef ADIABATIC
00095       pGrid->U[k][j][i].E = pres/Gamma_1
00096         + 0.5*(SQR(pGrid->U[k][j][i].M1) + SQR(pGrid->U[k][j][i].M2) 
00097              + SQR(pGrid->U[k][j][i].M3))/den;
00098 #endif
00099 
00100     }
00101   }}
00102 
00103 /* enroll gravitational potential of planet & shearing-box potential fns */
00104 
00105   StaticGravPot = PlanetPot;
00106   ShearingBoxPot = UnstratifiedDisk;
00107 
00108 /* enroll new history variables, only once  */
00109 
00110   if (frst == 1) {
00111     dump_history_enroll(hst_rho_Vx_dVy, "<rho Vx dVy>");
00112     dump_history_enroll(hst_rho_dVy2, "<rho dVy^2>");
00113 #ifdef ADIABATIC
00114     dump_history_enroll(hst_E_total, "<E + rho Phi>");
00115 #endif
00116     frst = 0;
00117   }
00118 
00119 /* With viscosity and/or resistivity, read eta_Ohm and nu_V */
00120 #ifdef NAVIER_STOKES
00121   nu_V = par_getd("problem","nu");
00122 #endif
00123 
00124 /* Enroll outflow BCs if perdiodic BCs NOT selected.  This assumes the root
00125  * level grid is specified by the <domain1> block in the input file */
00126 
00127   BCFlag = par_geti_def("domain1","bc_ix1",0);
00128   if (BCFlag != 4) {
00129     if (pDomain->Disp[0] == 0) bvals_mhd_fun(pDomain, left_x1,  constant_iib);
00130   }
00131   BCFlag = par_geti_def("domain1","bc_ox1",0);
00132   if (BCFlag != 4) {
00133     if (pDomain->MaxX[0] == pDomain->RootMaxX[0])
00134       bvals_mhd_fun(pDomain, right_x1, constant_oib);
00135   }
00136 
00137   return;
00138 }
00139 
00140 /*==============================================================================
00141  * PUBLIC PROBLEM USER FUNCTIONS:
00142  * problem_write_restart() - writes problem-specific user data to restart files
00143  * problem_read_restart()  - reads problem-specific user data from restart files
00144  * get_usr_expr()          - sets pointer to expression for special output data
00145  * get_usr_out_fun()       - returns a user defined output function pointer
00146  * get_usr_par_prop()      - returns a user defined particle selection function
00147  * Userwork_in_loop        - problem specific work IN     main loop
00148  * Userwork_after_loop     - problem specific work AFTER  main loop
00149  *----------------------------------------------------------------------------*/
00150 
00151 void problem_write_restart(MeshS *pM, FILE *fp)
00152 {
00153   return;
00154 }
00155 
00156 /* 'problem_read_restart' must enroll gravity on restarts */
00157 
00158 void problem_read_restart(MeshS *pM, FILE *fp)
00159 {
00160   int nl,nd,BCFlag_ix1,BCFlag_ox1;
00161 /* Read Omega, and with viscosity and/or resistivity, read eta_Ohm and nu_V */
00162 
00163 #ifdef SHEARING_BOX
00164   Omega_0 = par_getd_def("problem","omega",1.0e-3);
00165   qshear  = par_getd_def("problem","qshear",1.5);
00166 #endif
00167   Mp      = par_getd_def("problem","Mplanet",0.0);
00168   Xplanet = par_getd_def("problem","Xplanet",0.0);
00169   Yplanet = par_getd_def("problem","Yplanet",0.0);
00170   Zplanet = par_getd_def("problem","Zplanet",0.0);
00171   Rsoft   = par_getd_def("problem","Rsoft",0.1);
00172   ramp_time = 0.0;
00173   insert_time = par_getd_def("problem","insert_time",0.0);
00174 #ifdef NAVIER_STOKES
00175   nu_V = par_getd("problem","nu");
00176 #endif
00177 
00178 /* enroll gravitational potential of planet & shearing-box potential fns */
00179 
00180   StaticGravPot = PlanetPot;
00181   ShearingBoxPot = UnstratifiedDisk;
00182 
00183 /* enroll new history variables */
00184 
00185   dump_history_enroll(hst_rho_Vx_dVy, "<rho Vx dVy>");
00186   dump_history_enroll(hst_rho_dVy2, "<rho dVy^2>");
00187 #ifdef ADIABATIC
00188   dump_history_enroll(hst_E_total, "<E + rho Phi>");
00189 #endif
00190 
00191   BCFlag_ix1 = par_geti_def("domain1","bc_ix1",0);
00192   BCFlag_ox1 = par_geti_def("domain1","bc_ox1",0);
00193   for (nl=0; nl<(pM->NLevels); nl++){
00194     for (nd=0; nd<(pM->DomainsPerLevel[nl]); nd++){
00195       if (pM->Domain[nl][nd].Disp[0] == 0 && BCFlag_ix1 != 4) 
00196         bvals_mhd_fun(&(pM->Domain[nl][nd]), left_x1,  constant_iib);
00197       if (pM->Domain[nl][nd].MaxX[0] == pM->Domain[nl][nd].RootMaxX[0] 
00198           && BCFlag_ox1 != 4)
00199         bvals_mhd_fun(&(pM->Domain[nl][nd]), right_x1, constant_oib);
00200     }
00201   }
00202 
00203   return;
00204 }
00205 
00206 /* Get_user_expression computes dVy */
00207 ConsFun_t get_usr_expr(const char *expr)
00208 {
00209   if(strcmp(expr,"dVy")==0) return expr_dV2;
00210   return NULL;
00211 }
00212 
00213 VOutFun_t get_usr_out_fun(const char *name){
00214   return NULL;
00215 }
00216 
00217 void Userwork_in_loop(MeshS *pM)
00218 {
00219   ramp_time = pM->time;
00220 }
00221 
00222 void Userwork_after_loop(MeshS *pM)
00223 {
00224 }
00225 
00226 /*------------------------------------------------------------------------------
00227  * PlanetPot:
00228  */
00229 /*! \fn static Real PlanetPot(const Real x1, const Real x2, const Real x3)
00230  *  \brief static gravitational potential of planet */
00231 static Real PlanetPot(const Real x1, const Real x2, const Real x3)
00232 {
00233   Real rad,phi=0.0;
00234   rad = sqrt(SQR(x1-Xplanet) + SQR(x2-Yplanet) + SQR(x3-Zplanet));
00235   phi = -1.0*MIN(1.0,(ramp_time/(insert_time+0.0001)))*Mp/(rad+Rsoft);
00236   return phi;
00237 }
00238 
00239 /*------------------------------------------------------------------------------
00240  *! \fn static Real UnstratifiedDisk(const Real x1, const Real x2,const Real x3)
00241  *  \brief tidal potential in shearing box
00242  */
00243 
00244 static Real UnstratifiedDisk(const Real x1, const Real x2, const Real x3)
00245 {
00246   Real phi=0.0;
00247 #ifdef SHEARING_BOX
00248 #ifndef FARGO
00249   phi -= qshear*Omega_0*Omega_0*x1*x1;
00250 #endif
00251 #endif
00252   return phi;
00253 }
00254 
00255 /*----------------------------------------------------------------------------*/
00256 /*! \fn static Real expr_dV2(const GridS *pG, const int i, const int j, 
00257  *                           const int k)
00258  *  \brief Computes delta(Vy) 
00259  */
00260 
00261 static Real expr_dV2(const GridS *pG, const int i, const int j, const int k)
00262 {
00263   Real x1,x2,x3;
00264   cc_pos(pG,i,j,k,&x1,&x2,&x3);
00265 #ifdef SHEARING_BOX
00266 #ifdef FARGO
00267   return (pG->U[k][j][i].M2/pG->U[k][j][i].d);
00268 #else
00269   return (pG->U[k][j][i].M2/pG->U[k][j][i].d + qshear*Omega_0*x1);
00270 #endif
00271 #endif
00272 }
00273 
00274 /*---------------------------------------------------------------------------*/
00275 /* Hydro history variables:
00276  * hst_rho_Vx_dVy: Reynolds stress, added as history variable.
00277  * hst_rho_dVy2: KE in y-velocity fluctuations
00278  * hst_E_total: total energy (including tidal potential).
00279  */
00280 
00281 /*! \fn static Real hst_rho_Vx_dVy(const GridS *pG,const int i,const int j, 
00282  *                                 const int k)
00283  *  \brief Reynolds stress, added as history variable.*/
00284 static Real hst_rho_Vx_dVy(const GridS *pG,const int i,const int j, const int k)
00285 {
00286   Real x1,x2,x3;
00287   cc_pos(pG,i,j,k,&x1,&x2,&x3);
00288 #ifdef SHEARING_BOX
00289 #ifdef FARGO
00290   return pG->U[k][j][i].M1*(pG->U[k][j][i].M2/pG->U[k][j][i].d);
00291 #else
00292   return pG->U[k][j][i].M1*
00293     (pG->U[k][j][i].M2/pG->U[k][j][i].d + qshear*Omega_0*x1);
00294 #endif
00295 #endif
00296 }
00297 
00298 /*! \fn static Real hst_rho_dVy2(const GridS *pG, const int i, const int j,  
00299  *                              const int k)
00300  *  \brief KE in y-velocity fluctuations */
00301 static Real hst_rho_dVy2(const GridS *pG, const int i, const int j, const int k)
00302 {
00303   Real x1,x2,x3,dVy;
00304   cc_pos(pG,i,j,k,&x1,&x2,&x3);
00305 #ifdef SHEARING_BOX
00306 #ifdef FARGO
00307   dVy = (pG->U[k][j][i].M2/pG->U[k][j][i].d);
00308 #else
00309   dVy = (pG->U[k][j][i].M2/pG->U[k][j][i].d + qshear*Omega_0*x1);
00310 #endif
00311 #endif
00312   return pG->U[k][j][i].d*dVy*dVy;
00313 }
00314 
00315 #ifdef ADIABATIC
00316 /*! \fn static Real hst_E_total(const GridS *pG, const int i, const int j, 
00317  *                              const int k)
00318  *  \brief total energy (including tidal potential). */
00319 static Real hst_E_total(const GridS *pG, const int i, const int j, const int k)
00320 {
00321   Real x1,x2,x3,phi;
00322   cc_pos(pG,i,j,k,&x1,&x2,&x3);
00323   phi = UnstratifiedDisk(x1, x2, x3);
00324 
00325   return pG->U[k][j][i].E + pG->U[k][j][i].d*phi;
00326 }
00327 #endif /* ADIABATIC */
00328 
00329 /*---------------------------------------------------------------------------*/
00330 /*! \fn void constant_iib(GridS *pGrid)
00331  *  \brief Sets boundary condition on left X boundary (iib) to constant
00332  * state (initial values).
00333  */
00334 
00335 void constant_iib(GridS *pGrid)
00336 {
00337   int is = pGrid->is;
00338   int js = pGrid->js, je = pGrid->je;
00339   int ks = pGrid->ks, ke = pGrid->ke;
00340   int i,j,k;
00341   Real x1,x2,x3;
00342   Real den = 1.0, pres = 1.0e-6;
00343 
00344   for (k=ks; k<=ke; k++) {
00345     for (j=js; j<=je; j++) {
00346       for (i=1; i<=nghost; i++) {
00347       cc_pos(pGrid,(is-i),j,k,&x1,&x2,&x3);
00348 
00349 /* Initialize d, M, and P.  With FARGO do not initialize the background shear */
00350 
00351       pGrid->U[k][j][is-i].d  = den;
00352       pGrid->U[k][j][is-i].M1 = 0.0;
00353       pGrid->U[k][j][is-i].M2 = 0.0;
00354 #ifdef SHEARING_BOX
00355 #ifndef FARGO
00356       pGrid->U[k][j][is-i].M2 -= den*(qshear*Omega_0*x1);
00357 #endif
00358 #endif
00359       pGrid->U[k][j][is-i].M3 = 0.0;
00360 #ifdef ADIABATIC
00361       pGrid->U[k][j][is-i].E = pres/Gamma_1
00362         + 0.5*(SQR(pGrid->U[k][j][is-i].M1) + SQR(pGrid->U[k][j][is-i].M2)
00363              + SQR(pGrid->U[k][j][is-i].M3))/den;
00364 #endif
00365       }
00366     }
00367   } 
00368 }
00369 
00370 /*---------------------------------------------------------------------------*/
00371 /*! \fn void constant_oib(GridS *pGrid)
00372  *  \brief  Sets boundary condition on right X boundary (oib) to constant
00373  * state (initial values).
00374  */
00375 
00376 void constant_oib(GridS *pGrid)
00377 {
00378   int is = pGrid->is, ie = pGrid->ie;
00379   int js = pGrid->js, je = pGrid->je;
00380   int ks = pGrid->ks, ke = pGrid->ke;
00381   int i,j,k;
00382   Real x1,x2,x3;
00383   Real den = 1.0, pres = 1.0e-6;
00384 
00385   for (k=ks; k<=ke; k++) {
00386     for (j=js; j<=je; j++) {
00387       for (i=1; i<=nghost; i++) {
00388       cc_pos(pGrid,(ie+i),j,k,&x1,&x2,&x3);
00389 
00390 /* Initialize d, M, and P.  With FARGO do not initialize the background shear */
00391 
00392       pGrid->U[k][j][ie+i].d  = den;
00393       pGrid->U[k][j][ie+i].M1 = 0.0;
00394       pGrid->U[k][j][ie+i].M2 = 0.0;
00395 #ifdef SHEARING_BOX
00396 #ifndef FARGO
00397       pGrid->U[k][j][ie+i].M2 -= den*(qshear*Omega_0*x1);
00398 #endif
00399 #endif
00400       pGrid->U[k][j][ie+i].M3 = 0.0;
00401 #ifdef ADIABATIC
00402       pGrid->U[k][j][ie+i].E = pres/Gamma_1
00403         + 0.5*(SQR(pGrid->U[k][j][ie+i].M1) + SQR(pGrid->U[k][j][ie+i].M2)
00404              + SQR(pGrid->U[k][j][ie+i].M3))/den;
00405 #endif
00406       }
00407     }
00408   }
00409 }
00410