00001 #include "copyright.h"
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016 #include <float.h>
00017 #include <math.h>
00018 #include <stdio.h>
00019 #include <stdlib.h>
00020 #include "defs.h"
00021 #include "athena.h"
00022 #include "globals.h"
00023 #include "prototypes.h"
00024
00025
00026
00027
00028
00029
00030
00031
00032 static double ran2(long int *idum);
00033 static void pbc_ix1(Grid *pGrid);
00034 static void pbc_ox1(Grid *pGrid);
00035
00036
00037
00038
00039 void problem(Grid *pGrid, Domain *pDomain)
00040 {
00041 int i,is,ie,j,js,je,ks,nx1,nx2,iprob;
00042 long int iseed = -1;
00043 Real d0,p0,B0,v0,x1,x2,x3,x1min,x1max,Lx,k0,press,amp,kp;
00044
00045 is = pGrid->is; ie = pGrid->ie;
00046 js = pGrid->js; je = pGrid->je;
00047 ks = pGrid->ks;
00048 nx1 = (ie-is)+1;
00049 nx2 = (je-js)+1;
00050 if ((nx1 == 1) || (nx2 == 1)) {
00051 ath_error("[firehose]: This problem can only be run in 2D\n");
00052 }
00053 d0 = 1.0;
00054 p0 = 1.0;
00055 B0 = par_getd("problem","B0");
00056 v0 = par_getd("problem","v0");
00057 #ifdef BRAGINSKII
00058 nu_V = par_getd("problem","nu");
00059 #endif
00060 amp = par_getd("problem","amp");
00061 kp = par_getd("problem","kp");
00062 iprob = par_getd("problem","iprob");
00063
00064
00065 x1min = par_getd("grid","x1min");
00066 x1max = par_getd("grid","x1max");
00067 Lx = x1max - x1min;
00068 k0 = 2.0*PI/Lx;
00069
00070
00071
00072
00073 if (iprob == 1) {
00074 for (j=js; j<=je; j++) {
00075 for (i=is; i<=ie; i++) {
00076
00077 cc_pos(pGrid,i,j,ks,&x1,&x2,&x3);
00078
00079 pGrid->U[ks][j][i].d = d0;
00080 pGrid->U[ks][j][i].M1 = -d0*v0*cos(k0*x2)*sin(k0*x1);
00081 pGrid->U[ks][j][i].M2 = d0*v0*sin(k0*x2)*cos(k0*x1)+amp*cos(kp*x1);
00082 pGrid->U[ks][j][i].M3 = 0.0;
00083 #ifdef MHD
00084 pGrid->B1i[ks][j][i] = B0;
00085 pGrid->B2i[ks][j][i] = 0.0;
00086 pGrid->B3i[ks][j][i] = 0.0;
00087 pGrid->U[ks][j][i].B1c = B0;
00088 pGrid->U[ks][j][i].B2c = 0.0;
00089 pGrid->U[ks][j][i].B3c = 0.0;
00090 #endif
00091 press = p0 + 0.5*d0*v0*v0*(cos(k0*x1)*cos(k0*x1) + cos(k0*x2)*cos(k0*x2));
00092 pGrid->U[ks][j][i].E = press/Gamma_1
00093 #ifdef MHD
00094 + 0.5*(SQR(pGrid->U[ks][j][i].B1c) + SQR(pGrid->U[ks][j][i].B2c)
00095 + SQR(pGrid->U[ks][j][i].B3c))
00096 #endif
00097 + 0.5*(SQR(pGrid->U[ks][j][i].M1) + SQR(pGrid->U[ks][j][i].M2)
00098 + SQR(pGrid->U[ks][j][i].M3))/pGrid->U[ks][j][i].d;
00099 }
00100 }
00101 #ifdef MHD
00102
00103 for (j=js; j<=je; j++) {
00104 pGrid->B1i[ks][j][ie+1] = B0;
00105 }
00106 for (i=is; i<=ie; i++) {
00107 pGrid->B2i[ks][je+1][i] = 0.0;
00108 }
00109 #endif
00110 }
00111
00112
00113
00114
00115 if (iprob == 2) {
00116 for (j=js; j<=je; j++) {
00117 for (i=is; i<=ie; i++) {
00118
00119 cc_pos(pGrid,i,j,ks,&x1,&x2,&x3);
00120
00121 pGrid->U[ks][j][i].d = d0;
00122 pGrid->U[ks][j][i].M1 = d0*amp*(ran2(&iseed) - 0.5);
00123 pGrid->U[ks][j][i].M2 = d0*(amp*(ran2(&iseed) - 0.5) - 0.015*x1);
00124 pGrid->U[ks][j][i].M3 = 0.0;
00125 #ifdef MHD
00126 pGrid->B1i[ks][j][i] = B0;
00127 pGrid->B2i[ks][j][i] = B0;
00128 pGrid->B3i[ks][j][i] = 0.0;
00129 pGrid->U[ks][j][i].B1c = B0;
00130 pGrid->U[ks][j][i].B2c = B0;
00131 pGrid->U[ks][j][i].B3c = 0.0;
00132 #endif
00133 pGrid->U[ks][j][i].E = 0.1/Gamma_1
00134 #ifdef MHD
00135 + 0.5*(SQR(pGrid->U[ks][j][i].B1c) + SQR(pGrid->U[ks][j][i].B2c)
00136 + SQR(pGrid->U[ks][j][i].B3c))
00137 #endif
00138 + 0.5*(SQR(pGrid->U[ks][j][i].M1) + SQR(pGrid->U[ks][j][i].M2)
00139 + SQR(pGrid->U[ks][j][i].M3))/pGrid->U[ks][j][i].d;
00140 }
00141 }
00142 #ifdef MHD
00143
00144 for (j=js; j<=je; j++) {
00145 pGrid->B1i[ks][j][ie+1] = B0;
00146 }
00147 for (i=is; i<=ie; i++) {
00148 pGrid->B2i[ks][je+1][i] = B0;
00149 }
00150 #endif
00151
00152 set_bvals_mhd_fun(left_x1, pbc_ix1);
00153 set_bvals_mhd_fun(right_x1, pbc_ox1);
00154 }
00155
00156 return;
00157 }
00158
00159
00160
00161
00162
00163
00164
00165
00166
00167
00168
00169
00170 void problem_write_restart(Grid *pG, Domain *pD, FILE *fp)
00171 {
00172 return;
00173 }
00174
00175 void problem_read_restart(Grid *pG, Domain *pD, FILE *fp)
00176 {
00177 return;
00178 }
00179
00180 Gasfun_t get_usr_expr(const char *expr)
00181 {
00182 return NULL;
00183 }
00184
00185 VGFunout_t get_usr_out_fun(const char *name){
00186 return NULL;
00187 }
00188
00189 #ifdef PARTICLES
00190 PropFun_t get_usr_par_prop(const char *name)
00191 {
00192 return NULL;
00193 }
00194
00195 void gasvshift(const Real x1, const Real x2, const Real x3,
00196 Real *u1, Real *u2, Real *u3)
00197 {
00198 return;
00199 }
00200
00201 void Userforce_particle(Vector *ft, const Real x1, const Real x2, const Real x3,
00202 Real *w1, Real *w2, Real *w3)
00203 {
00204 return;
00205 }
00206 #endif
00207
00208 void Userwork_in_loop(Grid *pGrid, Domain *pDomain)
00209 {
00210 }
00211
00212 void Userwork_after_loop(Grid *pGrid, Domain *pDomain)
00213 {
00214 }
00215
00216
00217
00218 #define IM1 2147483563
00219 #define IM2 2147483399
00220 #define AM (1.0/IM1)
00221 #define IMM1 (IM1-1)
00222 #define IA1 40014
00223 #define IA2 40692
00224 #define IQ1 53668
00225 #define IQ2 52774
00226 #define IR1 12211
00227 #define IR2 3791
00228 #define NTAB 32
00229 #define NDIV (1+IMM1/NTAB)
00230 #define RNMX (1.0-DBL_EPSILON)
00231
00232
00233
00234
00235
00236
00237
00238
00239
00240
00241
00242
00243
00244
00245 double ran2(long int *idum)
00246 {
00247 int j;
00248 long int k;
00249 static long int idum2=123456789;
00250 static long int iy=0;
00251 static long int iv[NTAB];
00252 double temp;
00253
00254 if (*idum <= 0) {
00255 if (-(*idum) < 1) *idum=1;
00256 else *idum = -(*idum);
00257 idum2=(*idum);
00258 for (j=NTAB+7;j>=0;j--) {
00259 k=(*idum)/IQ1;
00260 *idum=IA1*(*idum-k*IQ1)-k*IR1;
00261 if (*idum < 0) *idum += IM1;
00262 if (j < NTAB) iv[j] = *idum;
00263 }
00264 iy=iv[0];
00265 }
00266 k=(*idum)/IQ1;
00267 *idum=IA1*(*idum-k*IQ1)-k*IR1;
00268 if (*idum < 0) *idum += IM1;
00269 k=idum2/IQ2;
00270 idum2=IA2*(idum2-k*IQ2)-k*IR2;
00271 if (idum2 < 0) idum2 += IM2;
00272 j=(int)(iy/NDIV);
00273 iy=iv[j]-idum2;
00274 iv[j] = *idum;
00275 if (iy < 1) iy += IMM1;
00276 if ((temp=AM*iy) > RNMX) return RNMX;
00277 else return temp;
00278 }
00279
00280 #undef IM1
00281 #undef IM2
00282 #undef AM
00283 #undef IMM1
00284 #undef IA1
00285 #undef IA2
00286 #undef IQ1
00287 #undef IQ2
00288 #undef IR1
00289 #undef IR2
00290 #undef NTAB
00291 #undef NDIV
00292 #undef RNMX
00293
00294
00295
00296
00297
00298
00299 static void pbc_ix1(Grid *pGrid)
00300 {
00301 int is = pGrid->is, ie = pGrid->ie;
00302 int js = pGrid->js, je = pGrid->je;
00303 int ks = pGrid->ks, ke = pGrid->ke;
00304 int i,j,k;
00305 #ifdef MHD
00306 int ju, ku;
00307 #endif
00308 Real xmin,xmax,Lx;
00309 xmin = par_getd("grid","x1min");
00310 xmax = par_getd("grid","x1max");
00311 Lx = xmax - xmin;
00312
00313 for (k=ks; k<=ke; k++) {
00314 for (j=js; j<=je; j++) {
00315 for (i=1; i<=nghost; i++) {
00316 pGrid->U[k][j][is-i] = pGrid->U[k][j][ie-(i-1)];
00317 pGrid->U[k][j][is-i].M2 =
00318 pGrid->U[k][j][is-i].M2 + pGrid->U[k][j][is-i].d*0.015*Lx;
00319 #ifdef SPECIAL_RELATIVITY
00320 pGrid->W[k][j][is-i] = pGrid->W[k][j][ie-(i-1)];
00321 #endif
00322 }
00323 }
00324 }
00325
00326 #ifdef MHD
00327 for (k=ks; k<=ke; k++) {
00328 for (j=js; j<=je; j++) {
00329 for (i=1; i<=nghost; i++) {
00330 pGrid->B1i[k][j][is-i] = pGrid->B1i[k][j][ie-(i-1)];
00331 }
00332 }
00333 }
00334
00335 if (pGrid->Nx2 > 1) ju=je+1; else ju=je;
00336 for (k=ks; k<=ke; k++) {
00337 for (j=js; j<=ju; j++) {
00338 for (i=1; i<=nghost; i++) {
00339 pGrid->B2i[k][j][is-i] = pGrid->B2i[k][j][ie-(i-1)];
00340 }
00341 }
00342 }
00343
00344 if (pGrid->Nx3 > 1) ku=ke+1; else ku=ke;
00345 for (k=ks; k<=ku; k++) {
00346 for (j=js; j<=je; j++) {
00347 for (i=1; i<=nghost; i++) {
00348 pGrid->B3i[k][j][is-i] = pGrid->B3i[k][j][ie-(i-1)];
00349 }
00350 }
00351 }
00352 #endif
00353
00354 return;
00355 }
00356
00357
00358
00359
00360
00361
00362 static void pbc_ox1(Grid *pGrid)
00363 {
00364 int is = pGrid->is, ie = pGrid->ie;
00365 int js = pGrid->js, je = pGrid->je;
00366 int ks = pGrid->ks, ke = pGrid->ke;
00367 int i,j,k;
00368 #ifdef MHD
00369 int ju, ku;
00370 #endif
00371 Real xmin,xmax,Lx;
00372 xmin = par_getd("grid","x1min");
00373 xmax = par_getd("grid","x1max");
00374 Lx = xmax - xmin;
00375
00376 for (k=ks; k<=ke; k++) {
00377 for (j=js; j<=je; j++) {
00378 for (i=1; i<=nghost; i++) {
00379 pGrid->U[k][j][ie+i] = pGrid->U[k][j][is+(i-1)];
00380 pGrid->U[k][j][ie+i].M2 =
00381 pGrid->U[k][j][ie+i].M2 - pGrid->U[k][j][ie+i].d*0.015*Lx;
00382 #ifdef SPECIAL_RELATIVITY
00383 pGrid->W[k][j][ie+i] = pGrid->W[k][j][is+(i-1)];
00384 #endif
00385 }
00386 }
00387 }
00388
00389 #ifdef MHD
00390
00391 for (k=ks; k<=ke; k++) {
00392 for (j=js; j<=je; j++) {
00393 for (i=2; i<=nghost; i++) {
00394 pGrid->B1i[k][j][ie+i] = pGrid->B1i[k][j][is+(i-1)];
00395 }
00396 }
00397 }
00398
00399 if (pGrid->Nx2 > 1) ju=je+1; else ju=je;
00400 for (k=ks; k<=ke; k++) {
00401 for (j=js; j<=ju; j++) {
00402 for (i=1; i<=nghost; i++) {
00403 pGrid->B2i[k][j][ie+i] = pGrid->B2i[k][j][is+(i-1)];
00404 }
00405 }
00406 }
00407
00408 if (pGrid->Nx3 > 1) ku=ke+1; else ku=ke;
00409 for (k=ks; k<=ku; k++) {
00410 for (j=js; j<=je; j++) {
00411 for (i=1; i<=nghost; i++) {
00412 pGrid->B3i[k][j][ie+i] = pGrid->B3i[k][j][is+(i-1)];
00413 }
00414 }
00415 }
00416 #endif
00417
00418 return;
00419 }