00001 #include "copyright.h"
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021
00022
00023
00024
00025
00026
00027
00028
00029
00030
00031
00032
00033
00034
00035
00036
00037
00038
00039 #include <float.h>
00040 #include <math.h>
00041 #include <stdio.h>
00042 #include <stdlib.h>
00043 #include "defs.h"
00044 #include "athena.h"
00045 #include "globals.h"
00046 #include "prototypes.h"
00047
00048
00049
00050
00051
00052
00053
00054
00055
00056
00057
00058
00059 static double ran2(long int *idum);
00060 static void reflect_ix2(GridS *pGrid);
00061 static void reflect_ox2(GridS *pGrid);
00062 static void reflect_ix3(GridS *pGrid);
00063 static void reflect_ox3(GridS *pGrid);
00064 static Real grav_pot2(const Real x1, const Real x2, const Real x3);
00065 static Real grav_pot3(const Real x1, const Real x2, const Real x3);
00066
00067
00068
00069
00070
00071 void problem(DomainS *pDomain)
00072 {
00073 GridS *pGrid = pDomain->Grid;
00074 int i=0,j=0,k=0;
00075 int is,ie,js,je,ks,ke,iprob;
00076 long int iseed = -1;
00077 Real amp,x1,x2,x3,lx,ly,lz,rhoh,L_rot,fact;
00078 #ifdef MHD
00079 Real b0,angle;
00080 #endif
00081 int ixs, jxs, kxs;
00082
00083 is = pGrid->is; ie = pGrid->ie;
00084 js = pGrid->js; je = pGrid->je;
00085 ks = pGrid->ks; ke = pGrid->ke;
00086
00087 lx = pDomain->RootMaxX[0] - pDomain->RootMinX[0];
00088 ly = pDomain->RootMaxX[1] - pDomain->RootMinX[1];
00089 lz = pDomain->RootMaxX[2] - pDomain->RootMinX[2];
00090
00091
00092 ixs = pGrid->Disp[0];
00093 jxs = pGrid->Disp[1];
00094 kxs = pGrid->Disp[2];
00095 iseed = -1 - (ixs + pDomain->Nx[0]*(jxs + pDomain->Nx[1]*kxs));
00096
00097
00098 amp = par_getd("problem","amp");
00099 iprob = par_geti("problem","iprob");
00100 rhoh = par_getd_def("problem","rhoh",3.0);
00101
00102 L_rot = par_getd_def("problem","L_rot",0.0);
00103
00104
00105
00106 #ifdef MHD
00107 b0 = par_getd("problem","b0");
00108 angle = par_getd("problem","angle");
00109 angle = (angle/180.)*PI;
00110 #endif
00111
00112
00113
00114
00115
00116
00117
00118 if (pGrid->Nx[2] == 1) {
00119 for (k=ks; k<=ke; k++) {
00120 for (j=js; j<=je; j++) {
00121 for (i=is; i<=ie; i++) {
00122 cc_pos(pGrid,i,j,k,&x1,&x2,&x3);
00123 pGrid->U[k][j][i].d = 1.0;
00124 pGrid->U[k][j][i].E = (1.0/Gamma - 0.1*x2)/Gamma_1;
00125 pGrid->U[k][j][i].M1 = 0.0;
00126 if (iprob == 1) {
00127 pGrid->U[k][j][i].M2 = amp/4.0*
00128 (1.0+cos(2.0*PI*x1/lx))*(1.0+cos(2.0*PI*x2/ly));
00129 }
00130 else {
00131 pGrid->U[k][j][i].M2 = amp*(ran2(&iseed) - 0.5)*
00132 (1.0+cos(2.0*PI*x2/ly));
00133 }
00134 pGrid->U[k][j][i].M3 = 0.0;
00135 if (x2 > 0.0) {
00136 pGrid->U[k][j][i].d = 2.0;
00137 pGrid->U[k][j][i].M2 *= 2.0;
00138 pGrid->U[k][j][i].E = (1.0/Gamma - 0.2*x2)/Gamma_1;
00139 }
00140 pGrid->U[k][j][i].E+=0.5*SQR(pGrid->U[k][j][i].M2)/pGrid->U[k][j][i].d;
00141 #ifdef MHD
00142 pGrid->B1i[k][j][i] = b0;
00143 pGrid->U[k][j][i].B1c = b0;
00144 pGrid->U[k][j][i].E += 0.5*b0*b0;
00145 #endif
00146 }
00147 #ifdef MHD
00148 pGrid->B1i[k][j][ie+1] = b0;
00149 #endif
00150 }
00151 }
00152
00153
00154
00155
00156
00157
00158 StaticGravPot = grav_pot2;
00159 if (pDomain->Disp[1] == 0) bvals_mhd_fun(pDomain, left_x2, reflect_ix2);
00160 if (pDomain->MaxX[1] == pDomain->RootMaxX[1])
00161 bvals_mhd_fun(pDomain, right_x2, reflect_ox2);
00162
00163 }
00164
00165
00166
00167
00168
00169
00170
00171
00172
00173
00174
00175 if (pGrid->Nx[2] > 1) {
00176 for (k=ks; k<=ke; k++) {
00177 for (j=js; j<=je; j++) {
00178 for (i=is; i<=ie; i++) {
00179 cc_pos(pGrid,i,j,k,&x1,&x2,&x3);
00180 pGrid->U[k][j][i].d = 1.0;
00181 pGrid->U[k][j][i].E = (1.0/Gamma - 0.1*x3)/Gamma_1;
00182 pGrid->U[k][j][i].M1 = 0.0;
00183 pGrid->U[k][j][i].M2 = 0.0;
00184 if (iprob == 1) {
00185 pGrid->U[k][j][i].M3 = amp/8.0*(1.0+cos(2.0*PI*x1/lx))*
00186 (1.0+cos(2.0*PI*x2/ly))*(1.0+cos(2.0*PI*x3/lz));
00187 }
00188 else {
00189 pGrid->U[k][j][i].M3 = amp*(ran2(&iseed) - 0.5)*
00190 (1.0+cos(2.0*PI*x3/lz));
00191 }
00192 if (x3 > 0.0) {
00193 pGrid->U[k][j][i].d = rhoh;
00194 pGrid->U[k][j][i].M3 *= rhoh;
00195 pGrid->U[k][j][i].E = (1.0/Gamma - 0.1*rhoh*x3)/Gamma_1;
00196 }
00197 pGrid->U[k][j][i].E+=0.5*SQR(pGrid->U[k][j][i].M3)/pGrid->U[k][j][i].d;
00198 #ifdef MHD
00199 switch(iprob){
00200 case 3:
00201 if (x3 <= 0.0) {
00202 pGrid->B1i[k][j][i] = b0;
00203 if (i == ie) pGrid->B1i[k][j][ie+1] = b0;
00204 pGrid->U[k][j][i].B1c = b0;
00205 }
00206 break;
00207 case 4:
00208 if (x3 <= 0.0) {
00209 pGrid->B1i[k][j][i] = b0;
00210 if (i == ie) pGrid->B1i[k][j][ie+1] = b0;
00211 pGrid->U[k][j][i].B1c = b0;
00212 pGrid->U[k][j][i].E += 0.5*b0*b0;
00213 }
00214 else {
00215 pGrid->B1i[k][j][i] = b0*cos(angle);
00216 pGrid->B2i[k][j][i] = b0*sin(angle);
00217 if (i == ie) pGrid->B1i[k][j][ie+1] = b0*cos(angle);
00218 if (j == je) pGrid->B2i[k][je+1][i] = b0*sin(angle);
00219 pGrid->U[k][j][i].B1c = b0*cos(angle);
00220 pGrid->U[k][j][i].B2c = b0*sin(angle);
00221 pGrid->U[k][j][i].E += 0.5*b0*b0;
00222 }
00223 break;
00224 case 5:
00225 if (x3 <= (-L_rot/2.0)) {
00226 pGrid->B1i[k][j][i] = b0;
00227 if (i == ie) pGrid->B1i[k][j][ie+1] = b0;
00228 pGrid->U[k][j][i].B1c = b0;
00229 pGrid->U[k][j][i].E += 0.5*b0*b0;
00230 }
00231 else if (x3 >= (L_rot/2.0)) {
00232 pGrid->B1i[k][j][i] = b0*cos(angle);
00233 pGrid->B2i[k][j][i] = b0*sin(angle);
00234 if (i == ie) pGrid->B1i[k][j][ie+1] = b0*cos(angle);
00235 if (j == je) pGrid->B2i[k][je+1][i] = b0*sin(angle);
00236 pGrid->U[k][j][i].B1c = b0*cos(angle);
00237 pGrid->U[k][j][i].B2c = b0*sin(angle);
00238 pGrid->U[k][j][i].E += 0.5*b0*b0;
00239 }
00240 else {
00241 fact = ((L_rot/2.0)+x3)/L_rot;
00242 pGrid->B1i[k][j][i] = b0*cos(fact*angle);
00243 pGrid->B2i[k][j][i] = b0*sin(fact*angle);
00244 if (i == ie) pGrid->B1i[k][j][ie+1] = b0*cos(fact*angle);
00245 if (j == je) pGrid->B2i[k][je+1][i] = b0*sin(fact*angle);
00246 pGrid->U[k][j][i].B1c = b0*cos(fact*angle);
00247 pGrid->U[k][j][i].B2c = b0*sin(fact*angle);
00248 pGrid->U[k][j][i].E += 0.5*b0*b0;
00249 }
00250
00251 break;
00252 default:
00253 pGrid->B1i[k][j][i] = b0;
00254 if (i == ie) pGrid->B1i[k][j][ie+1] = b0;
00255 pGrid->U[k][j][i].B1c = b0;
00256 pGrid->U[k][j][i].E += 0.5*b0*b0;
00257 }
00258 #endif
00259 }
00260 }
00261 }
00262
00263
00264
00265
00266
00267
00268 StaticGravPot = grav_pot3;
00269
00270 if (pDomain->Disp[2] == 0) bvals_mhd_fun(pDomain, left_x3, reflect_ix3);
00271 if (pDomain->MaxX[2] == pDomain->RootMaxX[2])
00272 bvals_mhd_fun(pDomain, right_x3, reflect_ox3);
00273
00274 }
00275
00276 return;
00277 }
00278
00279
00280
00281
00282
00283
00284
00285
00286
00287
00288
00289
00290 void problem_write_restart(MeshS *pM, FILE *fp)
00291 {
00292 return;
00293 }
00294
00295
00296
00297
00298
00299
00300 void problem_read_restart(MeshS *pM, FILE *fp)
00301 {
00302 int nl,nd;
00303
00304 if (pM->Nx[2] == 1) {
00305 StaticGravPot = grav_pot2;
00306 for (nl=0; nl<(pM->NLevels); nl++){
00307 for (nd=0; nd<(pM->DomainsPerLevel[nl]); nd++){
00308 bvals_mhd_fun(&(pM->Domain[nl][nd]), left_x2, reflect_ix2);
00309 bvals_mhd_fun(&(pM->Domain[nl][nd]), right_x2, reflect_ox2);
00310 }
00311 }
00312 }
00313
00314 if (pM->Nx[2] > 1) {
00315 StaticGravPot = grav_pot3;
00316 for (nl=0; nl<(pM->NLevels); nl++){
00317 for (nd=0; nd<(pM->DomainsPerLevel[nl]); nd++){
00318 bvals_mhd_fun(&(pM->Domain[nl][nd]), left_x3, reflect_ix3);
00319 bvals_mhd_fun(&(pM->Domain[nl][nd]), right_x3, reflect_ox3);
00320 }
00321 }
00322 }
00323
00324 return;
00325 }
00326
00327 ConsFun_t get_usr_expr(const char *expr)
00328 {
00329 return NULL;
00330 }
00331
00332 VOutFun_t get_usr_out_fun(const char *name){
00333 return NULL;
00334 }
00335
00336 void Userwork_in_loop(MeshS *pM)
00337 {
00338 }
00339
00340 void Userwork_after_loop(MeshS *pM)
00341 {
00342 }
00343
00344
00345
00346
00347
00348 #define IM1 2147483563
00349 #define IM2 2147483399
00350 #define AM (1.0/IM1)
00351 #define IMM1 (IM1-1)
00352 #define IA1 40014
00353 #define IA2 40692
00354 #define IQ1 53668
00355 #define IQ2 52774
00356 #define IR1 12211
00357 #define IR2 3791
00358 #define NTAB 32
00359 #define NDIV (1+IMM1/NTAB)
00360 #define RNMX (1.0-DBL_EPSILON)
00361
00362
00363
00364
00365
00366
00367
00368
00369
00370
00371
00372
00373
00374
00375
00376 double ran2(long int *idum)
00377 {
00378 int j;
00379 long int k;
00380 static long int idum2=123456789;
00381 static long int iy=0;
00382 static long int iv[NTAB];
00383 double temp;
00384
00385 if (*idum <= 0) {
00386 if (-(*idum) < 1) *idum=1;
00387 else *idum = -(*idum);
00388 idum2=(*idum);
00389 for (j=NTAB+7;j>=0;j--) {
00390 k=(*idum)/IQ1;
00391 *idum=IA1*(*idum-k*IQ1)-k*IR1;
00392 if (*idum < 0) *idum += IM1;
00393 if (j < NTAB) iv[j] = *idum;
00394 }
00395 iy=iv[0];
00396 }
00397 k=(*idum)/IQ1;
00398 *idum=IA1*(*idum-k*IQ1)-k*IR1;
00399 if (*idum < 0) *idum += IM1;
00400 k=idum2/IQ2;
00401 idum2=IA2*(idum2-k*IQ2)-k*IR2;
00402 if (idum2 < 0) idum2 += IM2;
00403 j=(int)(iy/NDIV);
00404 iy=iv[j]-idum2;
00405 iv[j] = *idum;
00406 if (iy < 1) iy += IMM1;
00407 if ((temp=AM*iy) > RNMX) return RNMX;
00408 else return temp;
00409 }
00410
00411 #undef IM1
00412 #undef IM2
00413 #undef AM
00414 #undef IMM1
00415 #undef IA1
00416 #undef IA2
00417 #undef IQ1
00418 #undef IQ2
00419 #undef IR1
00420 #undef IR2
00421 #undef NTAB
00422 #undef NDIV
00423 #undef RNMX
00424
00425
00426
00427
00428
00429
00430 static void reflect_ix2(GridS *pGrid)
00431 {
00432 int js = pGrid->js;
00433 int ks = pGrid->ks, ke = pGrid->ke;
00434 int i,j,k,il,iu,ku;
00435
00436 iu = pGrid->ie + nghost;
00437 il = pGrid->is - nghost;
00438
00439 for (k=ks; k<=ke; k++) {
00440 for (j=1; j<=nghost; j++) {
00441 for (i=il; i<=iu; i++) {
00442 pGrid->U[k][js-j][i] = pGrid->U[k][js+(j-1)][i];
00443 pGrid->U[k][js-j][i].M2 = -pGrid->U[k][js-j][i].M2;
00444 pGrid->U[k][js-j][i].E +=
00445 pGrid->U[k][js+(j-1)][i].d*0.1*(2*j-1)*pGrid->dx2/Gamma_1;
00446 }
00447 }
00448 }
00449
00450 #ifdef MHD
00451 for (k=ks; k<=ke; k++) {
00452 for (j=1; j<=nghost; j++) {
00453 for (i=il; i<=iu; i++) {
00454 pGrid->B1i[k][js-j][i] = pGrid->B1i[k][js+(j-1)][i];
00455 }
00456 }
00457 }
00458
00459 for (k=ks; k<=ke; k++) {
00460 for (j=1; j<=nghost; j++) {
00461 for (i=il; i<=iu; i++) {
00462 pGrid->B2i[k][js-j][i] = pGrid->B2i[k][js+(j-1)][i];
00463 }
00464 }
00465 }
00466
00467 if (pGrid->Nx[2] > 1) ku=ke+1; else ku=ke;
00468 for (k=ks; k<=ku; k++) {
00469 for (j=1; j<=nghost; j++) {
00470 for (i=il; i<=iu; i++) {
00471 pGrid->B3i[k][js-j][i] = pGrid->B3i[k][js+(j-1)][i];
00472 }
00473 }
00474 }
00475 #endif
00476
00477 return;
00478 }
00479
00480
00481
00482
00483
00484
00485 static void reflect_ox2(GridS *pGrid)
00486 {
00487 int je = pGrid->je;
00488 int ks = pGrid->ks, ke = pGrid->ke, ku;
00489 int i,j,k,il,iu,jl,ju;
00490
00491 iu = pGrid->ie + nghost;
00492 il = pGrid->is - nghost;
00493
00494 if (pGrid->Nx[1] > 1){
00495 ju = pGrid->je + nghost;
00496 jl = pGrid->js - nghost;
00497 } else {
00498 ju = pGrid->je;
00499 jl = pGrid->js;
00500 }
00501
00502 for (k=ks; k<=ke; k++) {
00503 for (j=1; j<=nghost; j++) {
00504 for (i=il; i<=iu; i++) {
00505 pGrid->U[k][je+j][i] = pGrid->U[k][je-(j-1)][i];
00506 pGrid->U[k][je+j][i].M2 = -pGrid->U[k][je+j][i].M2;
00507 pGrid->U[k][je+j][i].E -=
00508 pGrid->U[k][je-(j-1)][i].d*0.1*(2*j-1)*pGrid->dx2/Gamma_1;
00509 }
00510 }
00511 }
00512
00513 #ifdef MHD
00514 for (k=ks; k<=ke; k++) {
00515 for (j=1; j<=nghost; j++) {
00516 for (i=il; i<=iu; i++) {
00517 pGrid->B1i[k][je+j][i] = pGrid->B1i[k][je-(j-1)][i];
00518 }
00519 }
00520 }
00521
00522
00523 for (k=ks; k<=ke; k++) {
00524 for (j=2; j<=nghost; j++) {
00525 for (i=il; i<=iu; i++) {
00526 pGrid->B2i[k][je+j][i] = pGrid->B2i[k][je-(j-2)][i];
00527 }
00528 }
00529 }
00530
00531 if (pGrid->Nx[2] > 1) ku=ke+1; else ku=ke;
00532 for (k=ks; k<=ku; k++) {
00533 for (j=1; j<=nghost; j++) {
00534 for (i=il; i<=iu; i++) {
00535 pGrid->B3i[k][je+j][i] = pGrid->B3i[k][je-(j-1)][i];
00536 }
00537 }
00538 }
00539 #endif
00540
00541 return;
00542 }
00543
00544
00545
00546
00547
00548
00549 static void reflect_ix3(GridS *pGrid)
00550 {
00551 int ks = pGrid->ks;
00552 int i,j,k,il,iu,jl,ju;
00553
00554 iu = pGrid->ie + nghost;
00555 il = pGrid->is - nghost;
00556 if (pGrid->Nx[1] > 1){
00557 ju = pGrid->je + nghost;
00558 jl = pGrid->js - nghost;
00559 } else {
00560 ju = pGrid->je;
00561 jl = pGrid->js;
00562 }
00563
00564 for (k=1; k<=nghost; k++) {
00565 for (j=jl; j<=ju; j++) {
00566 for (i=il; i<=iu; i++) {
00567 pGrid->U[ks-k][j][i] = pGrid->U[ks+(k-1)][j][i];
00568 pGrid->U[ks-k][j][i].M3 = -pGrid->U[ks-k][j][i].M3;
00569 pGrid->U[ks-k][j][i].E +=
00570 pGrid->U[ks+(k-1)][j][i].d*0.1*(2*k-1)*pGrid->dx3/Gamma_1;
00571 }
00572 }
00573 }
00574
00575 #ifdef MHD
00576 for (k=1; k<=nghost; k++) {
00577 for (j=jl; j<=ju; j++) {
00578 for (i=il; i<=iu; i++) {
00579 pGrid->B1i[ks-k][j][i] = pGrid->B1i[ks+(k-1)][j][i];
00580 }
00581 }
00582 }
00583 for (k=1; k<=nghost; k++) {
00584 for (j=jl; j<=ju; j++) {
00585 for (i=il; i<=iu; i++) {
00586 pGrid->B2i[ks-k][j][i] = pGrid->B2i[ks+(k-1)][j][i];
00587 }
00588 }
00589 }
00590
00591 for (k=1; k<=nghost; k++) {
00592 for (j=jl; j<=ju; j++) {
00593 for (i=il; i<=iu; i++) {
00594 pGrid->B3i[ks-k][j][i] = pGrid->B3i[ks+(k-1)][j][i];
00595 }
00596 }
00597 }
00598 #endif
00599
00600 return;
00601 }
00602
00603
00604
00605
00606
00607
00608 static void reflect_ox3(GridS *pGrid)
00609 {
00610 int ke = pGrid->ke;
00611 int i,j,k ,il,iu,jl,ju;
00612
00613 iu = pGrid->ie + nghost;
00614 il = pGrid->is - nghost;
00615 if (pGrid->Nx[1] > 1){
00616 ju = pGrid->je + nghost;
00617 jl = pGrid->js - nghost;
00618 } else {
00619 ju = pGrid->je;
00620 jl = pGrid->js;
00621 }
00622 for (k=1; k<=nghost; k++) {
00623 for (j=jl; j<=ju; j++) {
00624 for (i=il; i<=iu; i++) {
00625 pGrid->U[ke+k][j][i] = pGrid->U[ke-(k-1)][j][i];
00626 pGrid->U[ke+k][j][i].M3 = -pGrid->U[ke+k][j][i].M3;
00627 pGrid->U[ke+k][j][i].E -=
00628 pGrid->U[ke-(k-1)][j][i].d*0.1*(2*k-1)*pGrid->dx3/Gamma_1;
00629 }
00630 }
00631 }
00632
00633 #ifdef MHD
00634 for (k=1; k<=nghost; k++) {
00635 for (j=jl; j<=ju; j++) {
00636 for (i=il; i<=iu; i++) {
00637 pGrid->B1i[ke+k][j][i] = pGrid->B1i[ke-(k-1)][j][i];
00638 }
00639 }
00640 }
00641
00642 for (k=1; k<=nghost; k++) {
00643 for (j=jl; j<=ju; j++) {
00644 for (i=il; i<=iu; i++) {
00645 pGrid->B2i[ke+k][j][i] = pGrid->B2i[ke-(k-1)][j][i];
00646 }
00647 }
00648 }
00649
00650
00651 for (k=2; k<=nghost; k++) {
00652 for (j=jl; j<=ju; j++) {
00653 for (i=il; i<=iu; i++) {
00654 pGrid->B3i[ke+k][j][i] = pGrid->B3i[ke-(k-1)][j][i];
00655 }
00656 }
00657 }
00658 #endif
00659
00660 return;
00661 }
00662
00663
00664
00665
00666
00667
00668 static Real grav_pot2(const Real x1, const Real x2, const Real x3)
00669 {
00670 return 0.1*x2;
00671 }
00672
00673
00674
00675 static Real grav_pot3(const Real x1, const Real x2, const Real x3)
00676 {
00677 return 0.1*x3;
00678 }