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 #include <float.h>
00028 #include <math.h>
00029 #include <stdio.h>
00030 #include <stdlib.h>
00031 #include <string.h>
00032 #include "defs.h"
00033 #include "athena.h"
00034 #include "prototypes.h"
00035 #include "globals.h"
00036
00037 #ifndef ISOTHERMAL
00038 #error Problem generator only works for isothermal turbulence
00039 #endif
00040
00041 #ifdef MPI_PARALLEL
00042 #include "mpi.h"
00043 #ifdef DOUBLE_PREC
00044 #define MPI_RL MPI_DOUBLE
00045 #else
00046 #define MPI_RL MPI_FLOAT
00047 #endif
00048 #endif
00049
00050
00051
00052
00053
00054
00055
00056
00057 #define OFST(i, j, k) ((k) + nx3*((j) + nx2*(i)))
00058
00059 #define KWVM(i, j, k) (sqrt(SQR(KCOMP(i,gis,gnx1))+ \
00060 SQR(KCOMP(j,gjs,gnx2))+SQR(KCOMP(k,gks,gnx3))))
00061
00062
00063
00064
00065 static struct ath_3d_fft_plan *plan;
00066
00067
00068 static ath_fft_data *fv1=NULL, *fv2=NULL, *fv3=NULL;
00069
00070
00071 static Real ***dv1=NULL, ***dv2=NULL, ***dv3=NULL;
00072
00073 static Real klow,khigh,kpeak,expo,dkx;
00074
00075
00076
00077 static Real dedt,tdrive,dtdrive;
00078
00079 static int ispect,idrive;
00080
00081 static int nx1,nx2,nx3,gnx1,gnx2,gnx3;
00082
00083 static int gis,gie,gjs,gje,gks,gke;
00084
00085 long int rseed;
00086 #ifdef MHD
00087
00088
00089 static Real beta,B0;
00090 #endif
00091
00092 static const Real rhobar = 1.0;
00093
00094
00095
00096
00097
00098 static void pspect(ath_fft_data *ampl);
00099 static void project();
00100 static inline void transform();
00101 static inline void generate();
00102 static void perturb(Grid *pGrid, Real dt);
00103
00104
00105 static void initialize(Grid *pGrid, Domain *pD);
00106
00107
00108
00109
00110
00111
00112
00113
00114 static Real hst_dEk(const Grid *pG, const int i, const int j, const int k);
00115 static Real hst_dEb(const Grid *pG, const int i, const int j, const int k);
00116
00117
00118 static double ran2(long int *idum);
00119
00120
00121
00122
00123
00124
00125
00126
00127
00128
00129
00130
00131
00132
00133
00134
00135 static void pspect(ath_fft_data *ampl)
00136 {
00137 int i,j,k;
00138 double q1,q2,q3;
00139
00140
00141 for (k=0; k<nx3; k++) {
00142 for (j=0; j<nx2; j++) {
00143 for (i=0; i<nx1; i++) {
00144 q1 = ran2(&rseed);
00145 q2 = ran2(&rseed);
00146 q3 = sqrt(-2.0*log(q1+1.0e-20))*cos(2.0*PI*q2);
00147 q1 = ran2(&rseed);
00148 ampl[OFST(i,j,k)][0] = q3*cos(2.0*PI*q1);
00149 ampl[OFST(i,j,k)][1] = q3*sin(2.0*PI*q1);
00150 }
00151 }
00152 }
00153
00154
00155
00156
00157
00158 for (k=0; k<nx3; k++) {
00159 for (j=0; j<nx2; j++) {
00160 for (i=0; i<nx1; i++) {
00161
00162 q3 = KWVM(i,j,k);
00163 if ((q3 > klow) && (q3 < khigh)) {
00164 q3 *= dkx;
00165 if (ispect == 1) {
00166
00167 ampl[OFST(i,j,k)][0] /= pow(q3,(expo+2.0)/2.0);
00168 ampl[OFST(i,j,k)][1] /= pow(q3,(expo+2.0)/2.0);
00169 } else if (ispect == 2) {
00170
00171 ampl[OFST(i,j,k)][0] *= pow(q3,3.0)*exp(-4.0*q3/kpeak);
00172 ampl[OFST(i,j,k)][1] *= pow(q3,3.0)*exp(-4.0*q3/kpeak);
00173 }
00174 } else {
00175
00176 ampl[OFST(i,j,k)][0] = 0.0;
00177 ampl[OFST(i,j,k)][1] = 0.0;
00178 }
00179 }
00180 }
00181 }
00182 ampl[0][0] = 0.0;
00183 ampl[0][1] = 0.0;
00184
00185 return;
00186 }
00187
00188
00189
00190
00191
00192
00193 static void project()
00194 {
00195 int i,j,k,m,ind;
00196 double kap[3], kapn[3], mag;
00197 ath_fft_data dot;
00198
00199
00200 for (k=0; k<nx3; k++) {
00201 kap[2] = sin(2.0*PI*(gks+k)/gnx3);
00202 for (j=0; j<nx2; j++) {
00203 kap[1] = sin(2.0*PI*(gjs+j)/gnx2);
00204 for (i=0; i<nx1; i++) {
00205 if (((gis+i)+(gjs+j)+(gks+k)) != 0) {
00206 kap[0] = sin(2.0*PI*(gis+i)/gnx1);
00207 ind = OFST(i,j,k);
00208
00209
00210 mag = sqrt(SQR(kap[0]) + SQR(kap[1]) + SQR(kap[2]));
00211 for (m=0; m<3; m++) kapn[m] = kap[m] / mag;
00212
00213
00214 dot[0] = fv1[ind][0]*kapn[0]+fv2[ind][0]*kapn[1]+fv3[ind][0]*kapn[2];
00215 dot[1] = fv1[ind][1]*kapn[0]+fv2[ind][1]*kapn[1]+fv3[ind][1]*kapn[2];
00216
00217
00218 fv1[ind][0] -= dot[0]*kapn[0];
00219 fv2[ind][0] -= dot[0]*kapn[1];
00220 fv3[ind][0] -= dot[0]*kapn[2];
00221
00222 fv1[ind][1] -= dot[1]*kapn[0];
00223 fv2[ind][1] -= dot[1]*kapn[1];
00224 fv3[ind][1] -= dot[1]*kapn[2];
00225 }
00226 }
00227 }
00228 }
00229
00230 return;
00231 }
00232
00233
00234
00235
00236
00237
00238 static inline void transform()
00239 {
00240
00241 ath_3d_fft(plan, fv1);
00242 ath_3d_fft(plan, fv2);
00243 ath_3d_fft(plan, fv3);
00244
00245
00246
00247
00248
00249 return;
00250 }
00251
00252
00253
00254
00255
00256
00257 static inline void generate()
00258 {
00259
00260 pspect(fv1);
00261 pspect(fv2);
00262 pspect(fv3);
00263
00264
00265 project();
00266
00267
00268
00269 transform();
00270
00271 return;
00272 }
00273
00274
00275
00276
00277
00278
00279
00280 static void perturb(Grid *pGrid, Real dt)
00281 {
00282 int i, is=pGrid->is, ie = pGrid->ie;
00283 int j, js=pGrid->js, je = pGrid->je;
00284 int k, ks=pGrid->ks, ke = pGrid->ke;
00285 int ind, mpierr;
00286 Real dvol, aa, b, c, s, de, qa, v1, v2, v3;
00287 Real t0, t0ij, t0i, t1, t1ij, t1i;
00288 Real t2, t2ij, t2i, t3, t3ij, t3i;
00289 Real m[4], gm[4];
00290
00291
00292 dvol = 1.0/((Real)(gnx1*gnx2*gnx3));
00293 for (k=ks; k<=ke; k++) {
00294 for (j=js; j<=je; j++) {
00295 for (i=is; i<=ie; i++) {
00296 ind = OFST(i-is,j-js,k-ks);
00297 dv1[k][j][i] = fv1[ind][0]*dvol;
00298 dv2[k][j][i] = fv2[ind][0]*dvol;
00299 dv3[k][j][i] = fv3[ind][0]*dvol;
00300 }
00301 }
00302 }
00303
00304
00305 t0 = 0.0; t1 = 0.0; t2 = 0.0; t3 = 0.0;
00306 for (k=ks; k<=ke; k++) {
00307 t0ij = 0.0; t1ij = 0.0; t2ij = 0.0; t3ij = 0.0;
00308 for (j=js; j<=je; j++) {
00309 t0i = 0.0; t1i = 0.0; t2i = 0.0; t3i = 0.0;
00310 for (i=is; i<=ie; i++) {
00311 t0i += pGrid->U[k][j][i].d;
00312
00313
00314 t1i += pGrid->U[k][j][i].d * dv1[k][j][i];
00315 t2i += pGrid->U[k][j][i].d * dv2[k][j][i];
00316 t3i += pGrid->U[k][j][i].d * dv3[k][j][i];
00317 }
00318 t0ij += t0i; t1ij += t1i; t2ij += t2i; t3ij += t3i;
00319 }
00320 t0 += t0ij; t1 += t1ij; t2 += t2ij; t3 += t3ij;
00321 }
00322
00323 #ifdef MPI_PARALLEL
00324
00325 m[0] = t0; m[1] = t1; m[2] = t2; m[3] = t3;
00326 mpierr = MPI_Allreduce(m, gm, 4, MPI_RL, MPI_SUM, MPI_COMM_WORLD);
00327 if (mpierr) ath_error("[normalize]: MPI_Allreduce error = %d\n", mpierr);
00328 t0 = gm[0]; t1 = gm[1]; t2 = gm[2]; t3 = gm[3];
00329 #endif
00330
00331
00332
00333 for (k=ks; k<=ke; k++) {
00334 for (j=js; j<=je; j++) {
00335 for (i=is; i<=ie; i++) {
00336 dv1[k][j][i] -= t1/t0;
00337 dv2[k][j][i] -= t2/t0;
00338 dv3[k][j][i] -= t3/t0;
00339 }
00340 }
00341 }
00342
00343
00344 t1 = 0.0; t2 = 0.0;
00345 for (k=ks; k<=ke; k++) {
00346 t1ij = 0.0; t2ij = 0.0;
00347 for (j=js; j<=je; j++) {
00348 t1i = 0.0; t2i = 0.0;
00349 for (i=is; i<=ie; i++) {
00350
00351
00352 v1 = dv1[k][j][i];
00353 v2 = dv2[k][j][i];
00354 v3 = dv3[k][j][i];
00355
00356 t1i += (pGrid->U[k][j][i].d)*(SQR(v1) + SQR(v2) + SQR(v3));
00357 t2i += (pGrid->U[k][j][i].M1)*v1 + (pGrid->U[k][j][i].M2)*v2 +
00358 (pGrid->U[k][j][i].M3)*v3;
00359 }
00360 t1ij += t1i; t2ij += t2i;
00361 }
00362 t1 += t1ij; t2 += t2ij;
00363 }
00364
00365 #ifdef MPI_PARALLEL
00366
00367 m[0] = t1; m[1] = t2;
00368 mpierr = MPI_Allreduce(m, gm, 2, MPI_RL, MPI_SUM, MPI_COMM_WORLD);
00369 if (mpierr) ath_error("[normalize]: MPI_Allreduce error = %d\n", mpierr);
00370 t1 = gm[0]; t2 = gm[1];
00371 #endif
00372
00373
00374 dvol = pGrid->dx1*pGrid->dx2*pGrid->dx3;
00375 if (idrive == 0) {
00376
00377 de = dedt*dt;
00378 } else {
00379
00380 de = dedt;
00381 }
00382 aa = 0.5*t1;
00383 aa = MAX(aa,1.0e-20);
00384 b = t2;
00385 c = -de/dvol;
00386 if(b >= 0.0)
00387 s = (-2.0*c)/(b + sqrt(b*b - 4.0*aa*c));
00388 else
00389 s = (-b + sqrt(b*b - 4.0*aa*c))/(2.0*aa);
00390
00391 if (isnan(s)) ath_error("[perturb]: s is NaN!\n");
00392
00393
00394 for (k=ks; k<=ke; k++) {
00395 for (j=js; j<=je; j++) {
00396 for (i=is; i<=ie; i++) {
00397 qa = s*pGrid->U[k][j][i].d;
00398 pGrid->U[k][j][i].M1 += qa*dv1[k][j][i];
00399 pGrid->U[k][j][i].M2 += qa*dv2[k][j][i];
00400 pGrid->U[k][j][i].M3 += qa*dv3[k][j][i];
00401 }
00402 }
00403 }
00404
00405 return;
00406 }
00407
00408
00409
00410
00411 static void initialize(Grid *pGrid, Domain *pD)
00412 {
00413 int i, is=pGrid->is, ie = pGrid->ie;
00414 int j, js=pGrid->js, je = pGrid->je;
00415 int k, ks=pGrid->ks, ke = pGrid->ke;
00416 int nbuf, mpierr, nx1gh, nx2gh, nx3gh;
00417 float kwv, kpara, kperp;
00418 char donedrive = 0;
00419
00420
00421
00422
00423
00424
00425
00426
00427 nx1 = (ie-is+1);
00428 nx2 = (je-js+1);
00429 nx3 = (ke-ks+1);
00430
00431
00432 gnx1 = pD->ide - pD->ids + 1;
00433 gnx2 = pD->jde - pD->jds + 1;
00434 gnx3 = pD->kde - pD->kds + 1;
00435
00436
00437 gis=is+pGrid->idisp; gie=ie+pGrid->idisp;
00438 gjs=js+pGrid->jdisp; gje=je+pGrid->jdisp;
00439 gks=ks+pGrid->kdisp; gke=ke+pGrid->kdisp;
00440
00441
00442
00443 nx1gh = nx1 + 2*nghost;
00444 nx2gh = nx2 + 2*nghost;
00445 nx3gh = nx3 + 2*nghost;
00446
00447
00448
00449
00450
00451 dtdrive = par_getd("problem","dtdrive");
00452 #ifdef MHD
00453
00454 beta = par_getd("problem","beta");
00455
00456 B0 = sqrt(2.0*Iso_csound2*rhobar/beta);
00457 #endif
00458
00459 dedt = par_getd("problem","dedt");
00460
00461
00462 ispect = par_geti("problem","ispect");
00463 if (ispect == 1) {
00464 expo = par_getd("problem","expo");
00465 } else if (ispect == 2) {
00466 kpeak = par_getd("problem","kpeak")*2.0*PI;
00467 } else {
00468 ath_error("Invalid value for ispect\n");
00469 }
00470
00471 klow = par_getd("problem","klow");
00472 khigh = par_getd("problem","khigh");
00473 dkx = 2.0*PI/(pGrid->dx1*gnx1);
00474
00475
00476 idrive = par_geti("problem","idrive");
00477 if ((idrive < 0) || (idrive > 1)) ath_error("Invalid value for idrive\n");
00478
00479 if ((idrive == 1) && (pGrid->nstep > 0)) {
00480 donedrive = 1;
00481 }
00482
00483 if (donedrive == 0) {
00484
00485 if ((dv1=(Real***)calloc_3d_array(nx3gh,nx2gh,nx1gh,sizeof(Real)))==NULL) {
00486 ath_error("[problem]: Error allocating memory for vel pert\n");
00487 }
00488 if ((dv2=(Real***)calloc_3d_array(nx3gh,nx2gh,nx1gh,sizeof(Real)))==NULL) {
00489 ath_error("[problem]: Error allocating memory for vel pert\n");
00490 }
00491 if ((dv3=(Real***)calloc_3d_array(nx3gh,nx2gh,nx1gh,sizeof(Real)))==NULL) {
00492 ath_error("[problem]: Error allocating memory for vel pert\n");
00493 }
00494 }
00495
00496
00497 plan = ath_3d_fft_quick_plan(pGrid, pD, NULL, ATH_FFT_BACKWARD);
00498
00499
00500 if (donedrive == 0) {
00501 fv1 = ath_3d_fft_malloc(plan);
00502 fv2 = ath_3d_fft_malloc(plan);
00503 fv3 = ath_3d_fft_malloc(plan);
00504 }
00505
00506
00507 dump_history_enroll(hst_dEk,"<dE_K>");
00508 dump_history_enroll(hst_dEb,"<dE_B>");
00509
00510 return;
00511 }
00512
00513
00514
00515
00516
00517
00518
00519
00520
00521 void problem(Grid *pGrid, Domain *pD)
00522 {
00523 int i, is=pGrid->is, ie = pGrid->ie;
00524 int j, js=pGrid->js, je = pGrid->je;
00525 int k, ks=pGrid->ks, ke = pGrid->ke;
00526
00527 rseed = (pGrid->my_id+1);
00528 initialize(pGrid, pD);
00529 tdrive = 0.0;
00530
00531
00532 for (k=ks-nghost; k<=ke+nghost; k++) {
00533 for (j=js-nghost; j<=je+nghost; j++) {
00534 for (i=is-nghost; i<=ie+nghost; i++) {
00535 pGrid->U[k][j][i].d = rhobar;
00536 pGrid->U[k][j][i].M1 = 0.0;
00537 pGrid->U[k][j][i].M2 = 0.0;
00538 pGrid->U[k][j][i].M3 = 0.0;
00539 }
00540 }
00541 }
00542
00543 #ifdef MHD
00544
00545 for (k=ks-nghost; k<=ke+nghost; k++) {
00546 for (j=js-nghost; j<=je+nghost; j++) {
00547 for (i=is-nghost; i<=ie+nghost; i++) {
00548 pGrid->U[k][j][i].B1c = B0;
00549 pGrid->U[k][j][i].B2c = 0.0;
00550 pGrid->U[k][j][i].B3c = 0.0;
00551 pGrid->B1i[k][j][i] = B0;
00552 pGrid->B2i[k][j][i] = 0.0;
00553 pGrid->B3i[k][j][i] = 0.0;
00554 }
00555 }
00556 }
00557 #endif
00558
00559
00560
00561
00562 generate();
00563 perturb(pGrid, dtdrive);
00564
00565
00566 if (idrive == 1) {
00567 ath_pout(0,"De-allocating driving memory.\n");
00568
00569
00570 free_3d_array(dv1);
00571 free_3d_array(dv2);
00572 free_3d_array(dv3);
00573
00574
00575 ath_3d_fft_free(fv1);
00576 ath_3d_fft_free(fv2);
00577 ath_3d_fft_free(fv3);
00578 }
00579
00580 return;
00581 }
00582
00583
00584
00585
00586
00587
00588
00589
00590
00591 void Userwork_in_loop(Grid *pGrid, Domain *pD)
00592 {
00593 int i, is=pGrid->is, ie = pGrid->ie;
00594 int j, js=pGrid->js, je = pGrid->je;
00595 int k, ks=pGrid->ks, ke = pGrid->ke;
00596 Real newtime;
00597
00598 if (isnan(pGrid->dt)) ath_error("Time step is NaN!");
00599
00600 if (idrive == 0) {
00601
00602 newtime = pGrid->time + pGrid->dt;
00603
00604 #ifndef IMPULSIVE_DRIVING
00605
00606 perturb(pGrid, pGrid->dt);
00607 #endif
00608
00609 if (newtime >= (tdrive+dtdrive)) {
00610
00611
00612
00613 while ((tdrive+dtdrive) <= newtime) tdrive += dtdrive;
00614
00615 #ifdef IMPULSIVE_DRIVING
00616
00617 perturb(pGrid, dtdrive);
00618 #endif
00619
00620
00621
00622
00623
00624
00625
00626 generate();
00627 }
00628 }
00629
00630 return;
00631 }
00632
00633
00634
00635 void Userwork_after_loop(Grid *pGrid, Domain *pD)
00636 {
00637
00638
00639 return;
00640 }
00641
00642 void problem_write_restart(Grid *pG, Domain *pD, FILE *fp)
00643 { return; }
00644
00645 void problem_read_restart(Grid *pG, Domain *pD, FILE *fp)
00646 {
00647
00648 rseed = (pG->my_id+1);
00649 initialize(pG, pD);
00650 tdrive = pG->time;
00651
00652
00653 if (idrive == 0) generate();
00654
00655 return;
00656 }
00657
00658 Gasfun_t get_usr_expr(const char *expr)
00659 { return NULL; }
00660
00661 VGFunout_t get_usr_out_fun(const char *name){
00662 return NULL;
00663 }
00664
00665
00666
00667
00668
00669
00670
00671
00672
00673
00674
00675 static Real hst_dEk(const Grid *pG, const int i, const int j, const int k)
00676 {
00677 return 0.5*(pG->U[k][j][i].M1*pG->U[k][j][i].M1 +
00678 pG->U[k][j][i].M2*pG->U[k][j][i].M2 +
00679 pG->U[k][j][i].M3*pG->U[k][j][i].M3)/pG->U[k][j][i].d;
00680 }
00681
00682
00683
00684 static Real hst_dEb(const Grid *pG, const int i, const int j, const int k)
00685 {
00686 #ifdef MHD
00687 return 0.5*((pG->U[k][j][i].B1c*pG->U[k][j][i].B1c +
00688 pG->U[k][j][i].B2c*pG->U[k][j][i].B2c +
00689 pG->U[k][j][i].B3c*pG->U[k][j][i].B3c)-B0*B0);
00690 #else
00691 return 0.0;
00692 #endif
00693 }
00694
00695
00696
00697 #define IM1 2147483563
00698 #define IM2 2147483399
00699 #define AM (1.0/IM1)
00700 #define IMM1 (IM1-1)
00701 #define IA1 40014
00702 #define IA2 40692
00703 #define IQ1 53668
00704 #define IQ2 52774
00705 #define IR1 12211
00706 #define IR2 3791
00707 #define NDIV (1+IMM1/NTAB)
00708 #define RNMX (1.0-DBL_EPSILON)
00709 #define NTAB 32
00710
00711
00712
00713
00714
00715
00716
00717
00718
00719
00720
00721
00722
00723
00724
00725
00726 double ran2(long int *idum){
00727 int j;
00728 long int k;
00729 static long int idum2=123456789;
00730 static long int iy=0;
00731 static long int iv[NTAB];
00732 double temp;
00733
00734 if (*idum <= 0) {
00735 if (-(*idum) < 1) *idum=1;
00736 else *idum = -(*idum);
00737 idum2=(*idum);
00738 for (j=NTAB+7;j>=0;j--) {
00739 k=(*idum)/IQ1;
00740 *idum=IA1*(*idum-k*IQ1)-k*IR1;
00741 if (*idum < 0) *idum += IM1;
00742 if (j < NTAB) iv[j] = *idum;
00743 }
00744 iy=iv[0];
00745 }
00746 k=(*idum)/IQ1;
00747 *idum=IA1*(*idum-k*IQ1)-k*IR1;
00748 if (*idum < 0) *idum += IM1;
00749 k=idum2/IQ2;
00750 idum2=IA2*(idum2-k*IQ2)-k*IR2;
00751 if (idum2 < 0) idum2 += IM2;
00752 j=(int)(iy/NDIV);
00753 iy=iv[j]-idum2;
00754 iv[j] = *idum;
00755 if (iy < 1) iy += IMM1;
00756 if ((temp=AM*iy) > RNMX) return RNMX;
00757 else return temp;
00758 }
00759
00760 #undef IM1
00761 #undef IM2
00762 #undef AM
00763 #undef IMM1
00764 #undef IA1
00765 #undef IA2
00766 #undef IQ1
00767 #undef IQ2
00768 #undef IR1
00769 #undef IR2
00770 #undef NTAB
00771 #undef NDIV
00772 #undef RNMX
00773
00774 #undef OFST
00775 #undef KCOMP
00776 #undef KWVM