00001 #include "copyright.h"
00002
00003
00004
00005
00006
00007
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
00020
00021
00022
00023
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
00042
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;
00054
00055 #ifdef SHEARING_BOX
00056
00057 ShBoxCoord = xy;
00058 #endif
00059
00060
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
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
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
00104
00105 StaticGravPot = PlanetPot;
00106 ShearingBoxPot = UnstratifiedDisk;
00107
00108
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
00120 #ifdef NAVIER_STOKES
00121 nu_V = par_getd("problem","nu");
00122 #endif
00123
00124
00125
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
00142
00143
00144
00145
00146
00147
00148
00149
00150
00151 void problem_write_restart(MeshS *pM, FILE *fp)
00152 {
00153 return;
00154 }
00155
00156
00157
00158 void problem_read_restart(MeshS *pM, FILE *fp)
00159 {
00160 int nl,nd,BCFlag_ix1,BCFlag_ox1;
00161
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
00179
00180 StaticGravPot = PlanetPot;
00181 ShearingBoxPot = UnstratifiedDisk;
00182
00183
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
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
00228
00229
00230
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
00241
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
00257
00258
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
00276
00277
00278
00279
00280
00281
00282
00283
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
00299
00300
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
00317
00318
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
00328
00329
00330
00331
00332
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
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
00372
00373
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
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