00001 #include "copyright.h"
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011 #include <float.h>
00012 #include <math.h>
00013 #include <stdio.h>
00014 #include <stdlib.h>
00015 #include "defs.h"
00016 #include "athena.h"
00017 #include "globals.h"
00018 #include "prototypes.h"
00019
00020
00021
00022
00023
00024
00025 static double ran2(long int *idum);
00026
00027
00028
00029
00030
00031 void problem(DomainS *pDomain)
00032 {
00033 GridS *pGrid = pDomain->Grid;
00034 int i=0,j=0,k=0;
00035 int is,ie,js,je,ks,ke,iprob;
00036 Real amp,drat,vflow,b0,a,sigma,x1,x2,x3;
00037 long int iseed = -1;
00038
00039 is = pGrid->is; ie = pGrid->ie;
00040 js = pGrid->js; je = pGrid->je;
00041 ks = pGrid->ks; ke = pGrid->ke;
00042
00043
00044
00045 iprob = par_geti("problem","iprob");
00046 vflow = par_getd("problem","vflow");
00047 drat = par_getd("problem","drat");
00048 amp = par_getd("problem","amp");
00049 #ifdef MHD
00050 b0 = par_getd("problem","b0");
00051 #endif
00052
00053
00054
00055 if (iprob == 1) {
00056 for (k=ks; k<=ke; k++) {
00057 for (j=js; j<=je; j++) {
00058 for (i=is; i<=ie; i++) {
00059 cc_pos(pGrid,i,j,k,&x1,&x2,&x3);
00060 pGrid->U[k][j][i].d = 1.0;
00061 pGrid->U[k][j][i].M1 = vflow + amp*(ran2(&iseed) - 0.5);
00062 pGrid->U[k][j][i].M2 = amp*(ran2(&iseed) - 0.5);
00063 pGrid->U[k][j][i].M3 = 0.0;
00064 if (fabs(x2) < 0.25) {
00065 pGrid->U[k][j][i].d = drat;
00066 pGrid->U[k][j][i].M1 = -drat*(vflow + amp*(ran2(&iseed) - 0.5));
00067 pGrid->U[k][j][i].M2 = drat*amp*(ran2(&iseed) - 0.5);
00068 }
00069
00070 #ifndef BAROTROPIC
00071 pGrid->U[k][j][i].E = 2.5/Gamma_1
00072 + 0.5*(SQR(pGrid->U[k][j][i].M1) + SQR(pGrid->U[k][j][i].M2)
00073 + SQR(pGrid->U[k][j][i].M3))/pGrid->U[k][j][i].d;
00074 #endif
00075 #ifdef MHD
00076 pGrid->B1i[k][j][i] = b0;
00077 pGrid->U[k][j][i].B1c = b0;
00078 #ifndef BAROTROPIC
00079 pGrid->U[k][j][i].E += 0.5*b0*b0;
00080 #endif
00081 #endif
00082 }
00083 #ifdef MHD
00084 pGrid->B1i[k][j][ie+1] = b0;
00085 #endif
00086 }
00087 }
00088 }
00089
00090
00091
00092
00093
00094 if (iprob == 2) {
00095 a = 0.05;
00096 sigma = 0.2;
00097 for (k=ks; k<=ke; k++) {
00098 for (j=js; j<=je; j++) {
00099 for (i=is; i<=ie; i++) {
00100 cc_pos(pGrid,i,j,k,&x1,&x2,&x3);
00101 pGrid->U[k][j][i].d = 1.0;
00102 pGrid->U[k][j][i].M1 = vflow*tanh(x2/a);
00103 pGrid->U[k][j][i].M2 = amp*sin(2.0*PI*x1)*exp(-(x2*x2)/(sigma*sigma));
00104 pGrid->U[k][j][i].M3 = 0.0;
00105 #ifndef BAROTROPIC
00106 pGrid->U[k][j][i].E = 1.0/Gamma_1
00107 + 0.5*(SQR(pGrid->U[k][j][i].M1) + SQR(pGrid->U[k][j][i].M2)
00108 + SQR(pGrid->U[k][j][i].M3))/pGrid->U[k][j][i].d;
00109 #endif
00110 #ifdef MHD
00111 pGrid->B1i[k][j][i] = b0;
00112 pGrid->U[k][j][i].B1c = b0;
00113 #ifndef BAROTROPIC
00114 pGrid->U[k][j][i].E += 0.5*b0*b0;
00115 #endif
00116 #endif
00117
00118 #if (NSCALARS > 0)
00119 pGrid->U[k][j][i].s[0] = 0.0;
00120 if (x2 > 0) pGrid->U[k][j][i].s[0] = 1.0;
00121 #endif
00122 }
00123 #ifdef MHD
00124 pGrid->B1i[k][j][ie+1] = b0;
00125 #endif
00126 }
00127 }
00128 }
00129
00130
00131
00132 #ifdef OHMIC
00133 eta_Ohm = par_getd("problem","eta");
00134 #endif
00135 #ifdef NAVIER_STOKES
00136 nu_V = par_getd("problem","nu");
00137 #endif
00138 #ifdef BRAGINSKII
00139 nu_V = par_getd("problem","nu");
00140 #endif
00141
00142 }
00143
00144
00145
00146
00147
00148
00149
00150
00151
00152
00153
00154
00155 void problem_write_restart(MeshS *pM, FILE *fp)
00156 {
00157 return;
00158 }
00159
00160 void problem_read_restart(MeshS *pM, FILE *fp)
00161 {
00162 #ifdef OHMIC
00163 eta_Ohm = par_getd("problem","eta");
00164 #endif
00165 #ifdef NAVIER_STOKES
00166 nu_V = par_getd("problem","nu");
00167 #endif
00168 #ifdef BRAGINSKII
00169 nu_V = par_getd("problem","nu");
00170 #endif
00171 return;
00172 }
00173
00174 #if (NSCALARS > 0)
00175
00176
00177 static Real color(const GridS *pG, const int i, const int j, const int k)
00178 {
00179 return pG->U[k][j][i].s[0]/pG->U[k][j][i].d;
00180 }
00181 #endif
00182
00183 ConsFun_t get_usr_expr(const char *expr)
00184 {
00185 #if (NSCALARS > 0)
00186 if(strcmp(expr,"color")==0) return color;
00187 #endif
00188 return NULL;
00189 }
00190
00191 VOutFun_t get_usr_out_fun(const char *name){
00192 return NULL;
00193 }
00194
00195 void Userwork_in_loop(MeshS *pM)
00196 {
00197 return;
00198 }
00199
00200 void Userwork_after_loop(MeshS *pM)
00201 {
00202 return;
00203 }
00204
00205
00206
00207
00208
00209 #define IM1 2147483563
00210 #define IM2 2147483399
00211 #define AM (1.0/IM1)
00212 #define IMM1 (IM1-1)
00213 #define IA1 40014
00214 #define IA2 40692
00215 #define IQ1 53668
00216 #define IQ2 52774
00217 #define IR1 12211
00218 #define IR2 3791
00219 #define NTAB 32
00220 #define NDIV (1+IMM1/NTAB)
00221 #define RNMX (1.0-DBL_EPSILON)
00222
00223
00224
00225
00226
00227
00228
00229
00230
00231
00232
00233
00234
00235
00236 double ran2(long int *idum){
00237 int j;
00238 long int k;
00239 static long int idum2=123456789;
00240 static long int iy=0;
00241 static long int iv[NTAB];
00242 double temp;
00243
00244 if (*idum <= 0) {
00245 if (-(*idum) < 1) *idum=1;
00246 else *idum = -(*idum);
00247 idum2=(*idum);
00248 for (j=NTAB+7;j>=0;j--) {
00249 k=(*idum)/IQ1;
00250 *idum=IA1*(*idum-k*IQ1)-k*IR1;
00251 if (*idum < 0) *idum += IM1;
00252 if (j < NTAB) iv[j] = *idum;
00253 }
00254 iy=iv[0];
00255 }
00256 k=(*idum)/IQ1;
00257 *idum=IA1*(*idum-k*IQ1)-k*IR1;
00258 if (*idum < 0) *idum += IM1;
00259 k=idum2/IQ2;
00260 idum2=IA2*(idum2-k*IQ2)-k*IR2;
00261 if (idum2 < 0) idum2 += IM2;
00262 j=(int)(iy/NDIV);
00263 iy=iv[j]-idum2;
00264 iv[j] = *idum;
00265 if (iy < 1) iy += IMM1;
00266 if ((temp=AM*iy) > RNMX) return RNMX;
00267 else return temp;
00268 }
00269
00270 #undef IM1
00271 #undef IM2
00272 #undef AM
00273 #undef IMM1
00274 #undef IA1
00275 #undef IA2
00276 #undef IQ1
00277 #undef IQ2
00278 #undef IR1
00279 #undef IR2
00280 #undef NTAB
00281 #undef NDIV
00282 #undef RNMX