00001 #include "copyright.h"
00002
00003
00004
00005
00006
00007
00008 #include <float.h>
00009 #include <math.h>
00010 #include <stdio.h>
00011 #include <stdlib.h>
00012 #include "defs.h"
00013 #include "athena.h"
00014 #include "globals.h"
00015 #include "prototypes.h"
00016
00017
00018
00019
00020
00021
00022
00023
00024
00025 static Real grav_pot(const Real x1, const Real x2, const Real x3);
00026 static void outflow_ix1(Grid *pGrid);
00027 static void outflow_ox1(Grid *pGrid);
00028 static void outflow_ox2(Grid *pGrid);
00029 static void outflow_ix3(Grid *pGrid);
00030 static void outflow_ox3(Grid *pGrid);
00031
00032
00033
00034
00035
00036 void problem(Grid *pGrid, Domain *pDomain)
00037 {
00038 int i=0,j=0,k=0;
00039 int is,ie,js,je,ks,ke;
00040 int ifield;
00041 Real x1,x2,x3,r2,rad_bubble=0.25;;
00042 #ifdef MHD
00043 Real b0;
00044 #endif
00045 int Nx1, Nx2, Nx3;
00046 int ixs, jxs, kxs;
00047
00048 is = pGrid->is; ie = pGrid->ie;
00049 js = pGrid->js; je = pGrid->je;
00050 ks = pGrid->ks; ke = pGrid->ke;
00051
00052 Nx1 = par_geti("grid","Nx1");
00053 Nx2 = par_geti("grid","Nx2");
00054 Nx3 = par_geti("grid","Nx3");
00055
00056
00057 #ifdef MHD
00058 b0 = par_getd("problem","b0");
00059 ifield = par_getd("problem","ifield");
00060 #endif
00061
00062
00063
00064
00065 for (k=ks; k<=ke; k++) {
00066 for (j=js; j<=je; j++) {
00067 for (i=is; i<=ie; i++) {
00068 cc_pos(pGrid,i,j,k,&x1,&x2,&x3);
00069 r2 = x1*x1 + x2*x2 + x3*x3;
00070 pGrid->U[k][j][i].d = pow((1.0 + r2),-0.75);
00071 pGrid->U[k][j][i].M1 = 0.0;
00072 pGrid->U[k][j][i].M2 = 0.0;
00073 pGrid->U[k][j][i].M3 = 0.0;
00074 pGrid->U[k][j][i].E = pow((1.0 + r2),-0.75)/(Gamma*Gamma_1);
00075
00076 r2 = x1*x1 + (x2-0.25)*(x2-0.25) + x3*x3;
00077 if (r2 < (rad_bubble*rad_bubble)) {
00078 pGrid->U[k][j][i].d = 0.01;
00079 }
00080 #ifdef MHD //choose ifield
00081 if (ifield==1){
00082 pGrid->B1i[k][j][i] = b0;
00083 pGrid->U[k][j][i].B1c = b0;
00084 pGrid->U[k][j][i].E += 0.5*b0*b0;
00085 pGrid->B1i[k][j][ie+1] = b0;
00086 }
00087 if (ifield==2){
00088 pGrid->B2i[k][j][i] = b0;
00089 pGrid->U[k][j][i].B2c = b0;
00090 pGrid->U[k][j][i].E += 0.5*b0*b0;
00091 pGrid->B2i[k][je+1][i] = b0;
00092 }
00093 if (ifield==3){
00094 pGrid->B3i[k][j][i] = b0;
00095 pGrid->U[k][j][i].B3c = b0;
00096 pGrid->U[k][j][i].E += 0.5*b0*b0;
00097 pGrid->B3i[ke+1][je][i] = b0;
00098 }
00099 #endif
00100 }
00101 }
00102 }
00103
00104
00105
00106 StaticGravPot = grav_pot;
00107 set_bvals_mhd_fun(left_x1, outflow_ix1);
00108 set_bvals_mhd_fun(right_x1, outflow_ox1);
00109 set_bvals_mhd_fun(right_x2, outflow_ox2);
00110 set_bvals_mhd_fun(left_x3, outflow_ix3);
00111 set_bvals_mhd_fun(right_x3, outflow_ox3);
00112
00113
00114
00115 #ifdef OHMIC
00116 eta_Ohm = par_getd("problem","eta");
00117 #endif
00118 #ifdef NAVIER_STOKES
00119 nu_V = par_getd("problem","nu");
00120 #endif
00121 #ifdef BRAGINSKII
00122 nu_V = par_getd("problem","nu");
00123 #endif
00124
00125 return;
00126 }
00127
00128
00129
00130
00131
00132
00133
00134
00135
00136
00137
00138
00139 void problem_write_restart(Grid *pG, Domain *pD, FILE *fp)
00140 {
00141 return;
00142 }
00143
00144
00145
00146
00147
00148
00149 void problem_read_restart(Grid *pG, Domain *pD, FILE *fp)
00150 {
00151 StaticGravPot = grav_pot;
00152 set_bvals_mhd_fun(left_x1, outflow_ix1);
00153 set_bvals_mhd_fun(right_x1, outflow_ox1);
00154 set_bvals_mhd_fun(right_x2, outflow_ox2);
00155 set_bvals_mhd_fun(left_x3, outflow_ix3);
00156 set_bvals_mhd_fun(right_x3, outflow_ox3);
00157
00158 return;
00159 }
00160
00161 Gasfun_t get_usr_expr(const char *expr)
00162 {
00163 return NULL;
00164 }
00165
00166 VGFunout_t get_usr_out_fun(const char *name){
00167 return NULL;
00168 }
00169
00170 void Userwork_in_loop(Grid *pGrid, Domain *pDomain)
00171 {
00172 }
00173
00174 void Userwork_after_loop(Grid *pGrid, Domain *pDomain)
00175 {
00176 }
00177
00178
00179
00180
00181
00182 static Real grav_pot(const Real x1, const Real x2, const Real x3)
00183 {
00184 Real r2;
00185 r2 = x1*x1 + x2*x2 + x3*x3;
00186 return (0.75/Gamma)*log(1.0+r2);
00187 }
00188
00189
00190
00191
00192
00193
00194 static void outflow_ix1(Grid *pGrid)
00195 {
00196 int is = pGrid->is;
00197 int js = pGrid->js, je = pGrid->je;
00198 int ks = pGrid->ks, ke = pGrid->ke;
00199 int i,j,k;
00200 Real r2,x1,x2,x3,d0,p0;
00201 #ifdef MHD
00202 int ibc_x1,ju,ku;
00203 #endif
00204
00205 for (k=ks; k<=ke; k++) {
00206 for (j=js; j<=je; j++) {
00207 for (i=1; i<=nghost; i++) {
00208 cc_pos(pGrid,is,j,k,&x1,&x2,&x3);
00209 r2 = x1*x1 + x2*x2 + x3*x3;
00210 d0 = pGrid->U[k][j][is].d/pow((1.0 + r2),-0.75);
00211
00212 p0 = pGrid->U[k][j][is].E - 0.5*(SQR(pGrid->U[k][j][is].M1)
00213 + SQR(pGrid->U[k][j][is].M2) + SQR(pGrid->U[k][j][is].M3))
00214 /pGrid->U[k][j][is].d;
00215 p0 /= pow((1.0 + r2),-0.75);
00216
00217 cc_pos(pGrid,(is-i),j,k,&x1,&x2,&x3);
00218 r2 = x1*x1 + x2*x2 + x3*x3;
00219 pGrid->U[k][j][is-i].d = d0*pow((1.0 + r2),-0.75);
00220 pGrid->U[k][j][is-i].M1 = 0.0;
00221 pGrid->U[k][j][is-i].M2 = 0.0;
00222 pGrid->U[k][j][is-i].M3 = 0.0;
00223 pGrid->U[k][j][is-i].E = p0*pow((1.0 + r2),-0.75);
00224 }
00225 }
00226 }
00227
00228 #ifdef MHD
00229 for (k=ks; k<=ke; k++) {
00230 for (j=js; j<=je; j++) {
00231 for (i=1; i<=nghost; i++) {
00232 pGrid->B1i[k][j][is-i] = pGrid->B1i[k][j][is];
00233 pGrid->U[k][j][is-i].B1c = pGrid->U[k][j][is].B1c;
00234 }
00235 }
00236 }
00237
00238 if (pGrid->Nx2 > 1) ju=je+1; else ju=je;
00239 for (k=ks; k<=ke; k++) {
00240 for (j=js; j<=ju; j++) {
00241 for (i=1; i<=nghost; i++) {
00242 pGrid->B2i[k][j][is-i] = pGrid->B2i[k][j][is];
00243 pGrid->U[k][j][is-i].B2c = pGrid->U[k][j][is].B2c;
00244 }
00245 }
00246 }
00247
00248 if (pGrid->Nx3 > 1) ku=ke+1; else ku=ke;
00249 for (k=ks; k<=ku; k++) {
00250 for (j=js; j<=je; j++) {
00251 for (i=1; i<=nghost; i++) {
00252 pGrid->B3i[k][j][is-i] = pGrid->B3i[k][j][is];
00253 pGrid->U[k][j][is-i].B3c = pGrid->U[k][j][is].B3c;
00254 }
00255 }
00256 }
00257 #endif
00258
00259 return;
00260 }
00261
00262
00263
00264
00265
00266
00267 static void outflow_ox1(Grid *pGrid)
00268 {
00269 int ie = pGrid->ie;
00270 int js = pGrid->js, je = pGrid->je;
00271 int ks = pGrid->ks, ke = pGrid->ke;
00272 int i,j,k;
00273 Real r2,x1,x2,x3,d0,p0;
00274 #ifdef MHD
00275 int obc_x1,ju,ku;
00276 #endif
00277
00278 for (k=ks; k<=ke; k++) {
00279 for (j=js; j<=je; j++) {
00280 for (i=1; i<=nghost; i++) {
00281 cc_pos(pGrid,ie,j,k,&x1,&x2,&x3);
00282 r2 = x1*x1 + x2*x2 + x3*x3;
00283 d0 = pGrid->U[k][j][ie].d/pow((1.0 + r2),-0.75);
00284
00285 p0 = pGrid->U[k][j][ie].E - 0.5*(SQR(pGrid->U[k][j][ie].M1)
00286 + SQR(pGrid->U[k][j][ie].M2) + SQR(pGrid->U[k][j][ie].M3))
00287 /pGrid->U[k][j][ie].d;
00288 p0 /= pow((1.0 + r2),-0.75);
00289
00290 cc_pos(pGrid,(ie+i),j,k,&x1,&x2,&x3);
00291 r2 = x1*x1 + x2*x2 + x3*x3;
00292 pGrid->U[k][j][ie+i].d = d0*pow((1.0 + r2),-0.75);
00293 pGrid->U[k][j][ie+i].M1 = 0.0;
00294 pGrid->U[k][j][ie+i].M2 = 0.0;
00295 pGrid->U[k][j][ie+i].M3 = 0.0;
00296 pGrid->U[k][j][ie+i].E = p0*pow((1.0 + r2),-0.75);
00297 }
00298 }
00299 }
00300
00301 #ifdef MHD
00302
00303 for (k=ks; k<=ke; k++) {
00304 for (j=js; j<=je; j++) {
00305 for (i=1; i<=nghost; i++) {
00306 if (i>1) pGrid->B1i[k][j][ie+i] = pGrid->B1i[k][j][ie+1];
00307 pGrid->U[k][j][ie+i].B1c = pGrid->U[k][j][ie].B1c;
00308 }
00309 }
00310 }
00311
00312 if (pGrid->Nx2 > 1) ju=je+1; else ju=je;
00313 for (k=ks; k<=ke; k++) {
00314 for (j=js; j<=ju; j++) {
00315 for (i=1; i<=nghost; i++) {
00316 pGrid->B2i[k][j][ie+i] = pGrid->B2i[k][j][ie];
00317 pGrid->U[k][j][ie+i].B2c = pGrid->U[k][j][ie].B2c;
00318 }
00319 }
00320 }
00321
00322 if (pGrid->Nx3 > 1) ku=ke+1; else ku=ke;
00323 for (k=ks; k<=ku; k++) {
00324 for (j=js; j<=je; j++) {
00325 for (i=1; i<=nghost; i++) {
00326 pGrid->B3i[k][j][ie+i] = pGrid->B3i[k][j][ie];
00327 pGrid->U[k][j][ie+i].B3c = pGrid->U[k][j][ie].B3c;
00328 }
00329 }
00330 }
00331 #endif
00332
00333 return;
00334 }
00335
00336
00337
00338
00339
00340
00341 static void outflow_ox2(Grid *pGrid)
00342 {
00343 int je = pGrid->je;
00344 int ks = pGrid->ks, ke = pGrid->ke;
00345 int i,j,k,il,iu;
00346 Real r2,x1,x2,x3,d0,p0;
00347 #ifdef MHD
00348 int obc_x2,ku;
00349 #endif
00350
00351 if (pGrid->Nx1 > 1){
00352 iu = pGrid->ie + nghost;
00353 il = pGrid->is - nghost;
00354 } else {
00355 iu = pGrid->ie;
00356 il = pGrid->is;
00357 }
00358
00359 for (k=ks; k<=ke; k++) {
00360 for (j=1; j<=nghost; j++) {
00361 for (i=il; i<=iu; i++) {
00362 cc_pos(pGrid,i,je,k,&x1,&x2,&x3);
00363 r2 = x1*x1 + x2*x2 + x3*x3;
00364 d0 = pGrid->U[k][je][i].d/pow((1.0 + r2),-0.75);
00365
00366 p0 = pGrid->U[k][je][i].E - 0.5*(SQR(pGrid->U[k][je][i].M1)
00367 + SQR(pGrid->U[k][je][i].M2) + SQR(pGrid->U[k][je][i].M3))
00368 /pGrid->U[k][je][i].d;
00369 p0 /= pow((1.0 + r2),-0.75);
00370
00371 cc_pos(pGrid,i,(je+j),k,&x1,&x2,&x3);
00372 r2 = x1*x1 + x2*x2 + x3*x3;
00373 pGrid->U[k][je+j][i].d = d0*pow((1.0 + r2),-0.75);
00374 pGrid->U[k][je+j][i].M1 = 0.0;
00375 pGrid->U[k][je+j][i].M2 = 0.0;
00376 pGrid->U[k][je+j][i].M3 = 0.0;
00377 pGrid->U[k][je+j][i].E = p0*pow((1.0 + r2),-0.75);
00378 }
00379 }
00380 }
00381
00382 #ifdef MHD
00383 for (k=ks; k<=ke; k++) {
00384 for (j=1; j<=nghost; j++) {
00385 for (i=il; i<=iu; i++) {
00386 pGrid->B1i[k][je+j][i] = pGrid->B1i[k][je][i];
00387 pGrid->U[k][je+j][i].B1c = pGrid->U[k][je][i].B1c;
00388 }
00389 }
00390 }
00391
00392
00393 for (k=ks; k<=ke; k++) {
00394 for (j=1; j<=nghost; j++) {
00395 for (i=il; i<=iu; i++) {
00396 if (j>1) pGrid->B2i[k][je+j][i] = pGrid->B2i[k][je+1][i];
00397 pGrid->U[k][je+j][i].B2c = pGrid->U[k][je][i].B2c;
00398 }
00399 }
00400 }
00401
00402 if (pGrid->Nx3 > 1) ku=ke+1; else ku=ke;
00403 for (k=ks; k<=ku; k++) {
00404 for (j=1; j<=nghost; j++) {
00405 for (i=il; i<=iu; i++) {
00406 pGrid->B3i[k][je+j][i] = pGrid->B3i[k][je][i];
00407 pGrid->U[k][je+j][i].B3c = pGrid->U[k][je][i].B3c;
00408 }
00409 }
00410 }
00411 #endif
00412
00413 return;
00414 }
00415
00416
00417
00418
00419
00420
00421 static void outflow_ix3(Grid *pGrid)
00422 {
00423 int ks = pGrid->ks;
00424 int i,j,k,il,iu,jl,ju;
00425 Real r2,x1,x2,x3,d0,p0;
00426
00427 if (pGrid->Nx1 > 1){
00428 iu = pGrid->ie + nghost;
00429 il = pGrid->is - nghost;
00430 } else {
00431 iu = pGrid->ie;
00432 il = pGrid->is;
00433 }
00434 if (pGrid->Nx2 > 1){
00435 ju = pGrid->je + nghost;
00436 jl = pGrid->js - nghost;
00437 } else {
00438 ju = pGrid->je;
00439 jl = pGrid->js;
00440 }
00441
00442 for (k=1; k<=nghost; k++) {
00443 for (j=jl; j<=ju; j++) {
00444 for (i=il; i<=iu; i++) {
00445 cc_pos(pGrid,i,j,ks,&x1,&x2,&x3);
00446 r2 = x1*x1 + x2*x2 + x3*x3;
00447 d0 = pGrid->U[ks][j][i].d/pow((1.0 + r2),-0.75);
00448
00449 p0 = pGrid->U[ks][j][i].E - 0.5*(SQR(pGrid->U[ks][j][i].M1)
00450 + SQR(pGrid->U[ks][j][i].M2) + SQR(pGrid->U[ks][j][i].M3))
00451 /pGrid->U[ks][j][i].d;
00452 p0 /= pow((1.0 + r2),-0.75);
00453
00454 cc_pos(pGrid,i,j,(ks-k),&x1,&x2,&x3);
00455 r2 = x1*x1 + x2*x2 + x3*x3;
00456 pGrid->U[ks-k][j][i].d = d0*pow((1.0 + r2),-0.75);
00457 pGrid->U[ks-k][j][i].M1 = 0.0;
00458 pGrid->U[ks-k][j][i].M2 = 0.0;
00459 pGrid->U[ks-k][j][i].M3 = 0.0;
00460 pGrid->U[ks-k][j][i].E = p0*pow((1.0 + r2),-0.75);
00461 }
00462 }
00463 }
00464
00465 #ifdef MHD
00466 for (k=1; k<=nghost; k++) {
00467 for (j=jl; j<=ju; j++) {
00468 for (i=il; i<=iu; i++) {
00469 pGrid->B1i[ks-k][j][i] = pGrid->B1i[ks][j][i];
00470 pGrid->U[ks-k][j][i].B1c = pGrid->U[ks][j][i].B1c;
00471 }
00472 }
00473 }
00474
00475 for (k=1; k<=nghost; k++) {
00476 for (j=jl; j<=ju; j++) {
00477 for (i=il; i<=iu; i++) {
00478 pGrid->B2i[ks-k][j][i] = pGrid->B2i[ks][j][i];
00479 pGrid->U[ks-k][j][i].B2c = pGrid->U[ks][j][i].B2c;
00480 }
00481 }
00482 }
00483
00484 for (k=1; k<=nghost; k++) {
00485 for (j=jl; j<=ju; j++) {
00486 for (i=il; i<=iu; i++) {
00487 pGrid->B3i[ks-k][j][i] = pGrid->B3i[ks][j][i];
00488 pGrid->U[ks-k][j][i].B3c = pGrid->U[ks][j][i].B3c;
00489 }
00490 }
00491 }
00492 #endif
00493
00494 return;
00495 }
00496
00497
00498
00499
00500
00501
00502 static void outflow_ox3(Grid *pGrid)
00503 {
00504 int ke = pGrid->ke;
00505 int i,j,k ,il,iu,jl,ju;
00506 Real r2,x1,x2,x3,d0,p0;
00507
00508 if (pGrid->Nx1 > 1){
00509 iu = pGrid->ie + nghost;
00510 il = pGrid->is - nghost;
00511 } else {
00512 iu = pGrid->ie;
00513 il = pGrid->is;
00514 }
00515 if (pGrid->Nx2 > 1){
00516 ju = pGrid->je + nghost;
00517 jl = pGrid->js - nghost;
00518 } else {
00519 ju = pGrid->je;
00520 jl = pGrid->js;
00521 }
00522
00523 for (k=1; k<=nghost; k++) {
00524 for (j=jl; j<=ju; j++) {
00525 for (i=il; i<=iu; i++) {
00526 cc_pos(pGrid,i,j,ke,&x1,&x2,&x3);
00527 r2 = x1*x1 + x2*x2 + x3*x3;
00528 d0 = pGrid->U[ke][j][i].d/pow((1.0 + r2),-0.75);
00529
00530 p0 = pGrid->U[ke][j][i].E - 0.5*(SQR(pGrid->U[ke][j][i].M1)
00531 + SQR(pGrid->U[ke][j][i].M2) + SQR(pGrid->U[ke][j][i].M3))
00532 /pGrid->U[ke][j][i].d;
00533 p0 /= pow((1.0 + r2),-0.75);
00534
00535 cc_pos(pGrid,i,j,(ke+k),&x1,&x2,&x3);
00536 r2 = x1*x1 + x2*x2 + x3*x3;
00537 pGrid->U[ke+k][j][i].d = d0*pow((1.0 + r2),-0.75);
00538 pGrid->U[ke+k][j][i].M1 = 0.0;
00539 pGrid->U[ke+k][j][i].M2 = 0.0;
00540 pGrid->U[ke+k][j][i].M3 = 0.0;
00541 pGrid->U[ke+k][j][i].E = p0*pow((1.0 + r2),-0.75);
00542 }
00543 }
00544 }
00545
00546 #ifdef MHD
00547 for (k=1; k<=nghost; k++) {
00548 for (j=jl; j<=ju; j++) {
00549 for (i=il; i<=iu; i++) {
00550 pGrid->B1i[ke+k][j][i] = pGrid->B1i[ke][j][i];
00551 pGrid->U[ke+k][j][i].B1c = pGrid->U[ke][j][i].B1c;
00552 }
00553 }
00554 }
00555
00556
00557 for (k=1; k<=nghost; k++) {
00558 for (j=jl; j<=ju; j++) {
00559 for (i=il; i<=iu; i++) {
00560 pGrid->B2i[ke+k][j][i] = pGrid->B2i[ke][j][i];
00561 pGrid->U[ke+k][j][i].B2c = pGrid->U[ke][j][i].B2c;
00562 }
00563 }
00564 }
00565
00566
00567 for (k=1; k<=nghost; k++) {
00568 for (j=jl; j<=ju; j++) {
00569 for (i=il; i<=iu; i++) {
00570 if (k>1) pGrid->B3i[ke+k][j][i] = pGrid->B3i[ke+1][j][i];
00571 pGrid->U[ke+k][j][i].B3c = pGrid->U[ke][j][i].B3c;
00572 }
00573 }
00574 }
00575 #endif
00576
00577 return;
00578 }