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 #include <math.h>
00031 #include <stdio.h>
00032 #include <stdlib.h>
00033 #include <string.h>
00034 #include "defs.h"
00035 #include "athena.h"
00036 #include "globals.h"
00037 #include "prototypes.h"
00038
00039
00040
00041
00042 void problem(DomainS *pDomain)
00043 {
00044 GridS *pGrid = pDomain->Grid;
00045 int i=0,j=0,k=0;
00046 int is,ie,js,je,ks,ke,nx1,nx2,nx3,iprob;
00047 Real x1c,x2c,x3c,x1f,x2f,x3f;
00048 Real x1size,x2size,x3size,lambda=0.0,ang_2=0.0,sin_a2=0.0,cos_a2=1.0,x,y;
00049 Real rad,amp,vflow,drat,diag;
00050 Real ***az,***ay,***ax;
00051 #ifdef MHD
00052 int ku;
00053 #endif
00054 #if (NSCALARS > 0)
00055 int n;
00056 #endif
00057
00058 is = pGrid->is; ie = pGrid->ie;
00059 js = pGrid->js; je = pGrid->je;
00060 ks = pGrid->ks; ke = pGrid->ke;
00061 nx1 = (ie-is)+1 + 2*nghost;
00062 nx2 = (je-js)+1 + 2*nghost;
00063 nx3 = (ke-ks)+1 + 2*nghost;
00064
00065 if (((je-js) == 0)) {
00066 ath_error("[field_loop]: This problem can only be run in 2D or 3D\n");
00067 }
00068
00069 if ((ay = (Real***)calloc_3d_array(nx3, nx2, nx1, sizeof(Real))) == NULL) {
00070 ath_error("[field_loop]: Error allocating memory for vector pot\n");
00071 }
00072 if ((az = (Real***)calloc_3d_array(nx3, nx2, nx1, sizeof(Real))) == NULL) {
00073 ath_error("[field_loop]: Error allocating memory for vector pot\n");
00074 }
00075 if ((ax = (Real***)calloc_3d_array(nx3, nx2, nx1, sizeof(Real))) == NULL) {
00076 ath_error("[field_loop]: Error allocating memory for vector pot\n");
00077 }
00078
00079
00080
00081 rad = par_getd("problem","rad");
00082 amp = par_getd("problem","amp");
00083 vflow = par_getd("problem","vflow");
00084 drat = par_getd_def("problem","drat",1.0);
00085 iprob = par_getd("problem","iprob");
00086 #ifdef OHMIC
00087 eta_R = par_getd("problem","eta");
00088 #endif
00089 #ifdef ISOTROPIC_CONDUCTION
00090 kappa_T = par_getd("problem","kappa");
00091 #endif
00092 #ifdef ANISOTROPIC_CONDUCTION
00093 chi_C = par_getd("problem","chi");
00094 #endif
00095
00096
00097
00098
00099 if(iprob == 4){
00100 x1size = pDomain->RootMaxX[0] - pDomain->RootMinX[0];
00101 x3size = pDomain->RootMaxX[2] - pDomain->RootMinX[2];
00102
00103
00104
00105
00106
00107
00108 if(x1size == x3size){
00109 ang_2 = PI/4.0;
00110 cos_a2 = sin_a2 = sqrt(0.5);
00111 }
00112 else{
00113 ang_2 = atan(x1size/x3size);
00114 sin_a2 = sin(ang_2);
00115 cos_a2 = cos(ang_2);
00116 }
00117
00118 if (cos_a2 >= sin_a2) {
00119 lambda = x1size*cos_a2;
00120 } else {
00121 lambda = x3size*sin_a2;
00122 }
00123 }
00124
00125
00126
00127 for (k=ks; k<=ke+1; k++) {
00128 for (j=js; j<=je+1; j++) {
00129 for (i=is; i<=ie+1; i++) {
00130 cc_pos(pGrid,i,j,k,&x1c,&x2c,&x3c);
00131 x1f = x1c - 0.5*pGrid->dx1;
00132 x2f = x2c - 0.5*pGrid->dx2;
00133 x3f = x3c - 0.5*pGrid->dx3;
00134
00135
00136
00137 if(iprob==1) {
00138 ax[k][j][i] = 0.0;
00139 ay[k][j][i] = 0.0;
00140 if ((x1f*x1f + x2f*x2f) < rad*rad) {
00141 az[k][j][i] = amp*(rad - sqrt(x1f*x1f + x2f*x2f));
00142 } else {
00143 az[k][j][i] = 0.0;
00144 }
00145 }
00146
00147
00148
00149 if(iprob==2) {
00150 if ((x2f*x2f + x3f*x3f) < rad*rad) {
00151 ax[k][j][i] = amp*(rad - sqrt(x2f*x2f + x3f*x3f));
00152 } else {
00153 ax[k][j][i] = 0.0;
00154 }
00155 ay[k][j][i] = 0.0;
00156 az[k][j][i] = 0.0;
00157 }
00158
00159
00160
00161 if(iprob==3) {
00162 if ((x1f*x1f + x3f*x3f) < rad*rad) {
00163 ay[k][j][i] = amp*(rad - sqrt(x1f*x1f + x3f*x3f));
00164 } else {
00165 ay[k][j][i] = 0.0;
00166 }
00167 ax[k][j][i] = 0.0;
00168 az[k][j][i] = 0.0;
00169 }
00170
00171
00172
00173
00174
00175
00176
00177
00178
00179
00180
00181
00182
00183 if(iprob==4) {
00184 x = x1c*cos_a2 + x3f*sin_a2;
00185 y = x2f;
00186
00187 while(x > 0.5*lambda) x -= lambda;
00188 while(x < -0.5*lambda) x += lambda;
00189 if ((x*x + y*y) < rad*rad) {
00190 ax[k][j][i] = amp*(rad - sqrt(x*x + y*y))*(-sin_a2);
00191 } else {
00192 ax[k][j][i] = 0.0;
00193 }
00194
00195 ay[k][j][i] = 0.0;
00196
00197 x = x1f*cos_a2 + x3c*sin_a2;
00198 y = x2f;
00199
00200 while(x > 0.5*lambda) x -= lambda;
00201 while(x < -0.5*lambda) x += lambda;
00202 if ((x*x + y*y) < rad*rad) {
00203 az[k][j][i] = amp*(rad - sqrt(x*x + y*y))*(cos_a2);
00204 } else {
00205 az[k][j][i] = 0.0;
00206 }
00207 }
00208
00209
00210
00211 if(iprob==5) {
00212 ax[k][j][i] = 0.0;
00213 if ((x1f*x1f + x2c*x2c + x3f*x3f) < rad*rad) {
00214 ay[k][j][i] = amp*(rad - sqrt(x1f*x1f + x2c*x2c + x3f*x3f));
00215 } else {
00216 ay[k][j][i] = 0.0;
00217 }
00218 if ((x1f*x1f + x2f*x2f + x3c*x3c) < rad*rad) {
00219 az[k][j][i] = amp*(rad - sqrt(x1f*x1f + x2f*x2f + x3c*x3c));
00220 } else {
00221 az[k][j][i] = 0.0;
00222 }
00223 }
00224
00225 }}}
00226
00227
00228
00229
00230
00231 x1size = pDomain->RootMaxX[0] - pDomain->RootMinX[0];
00232 x2size = pDomain->RootMaxX[1] - pDomain->RootMinX[1];
00233 x3size = pDomain->RootMaxX[2] - pDomain->RootMinX[2];
00234 diag = sqrt(x1size*x1size + x2size*x2size + x3size*x3size);
00235 for (k=ks; k<=ke; k++) {
00236 for (j=js; j<=je; j++) {
00237 for (i=is; i<=ie; i++) {
00238 pGrid->U[k][j][i].d = 1.0;
00239 pGrid->U[k][j][i].M1 = pGrid->U[k][j][i].d*vflow*x1size/diag;
00240 pGrid->U[k][j][i].M2 = pGrid->U[k][j][i].d*vflow*x2size/diag;
00241 pGrid->U[k][j][i].M3 = pGrid->U[k][j][i].d*vflow*x3size/diag;
00242
00243 cc_pos(pGrid,i,j,k,&x1c,&x2c,&x3c);
00244 if ((x1c*x1c + x2c*x2c + x3c*x3c) < rad*rad) {
00245 pGrid->U[k][j][i].d = drat;
00246 pGrid->U[k][j][i].M1 = pGrid->U[k][j][i].d*vflow*x1size/diag;
00247 pGrid->U[k][j][i].M2 = pGrid->U[k][j][i].d*vflow*x2size/diag;
00248 pGrid->U[k][j][i].M3 = pGrid->U[k][j][i].d*vflow*x3size/diag;
00249 }
00250 #if (NSCALARS > 0)
00251 for (n=0; n<NSCALARS; n++) pGrid->U[k][j][i].s[n] = 0.0;
00252 if ((x1c*x1c + x2c*x2c + x3c*x3c) < rad*rad) {
00253 for (n=0; n<NSCALARS; n++) pGrid->U[k][j][i].s[n] = 1.0;
00254 }
00255 #endif
00256 }}}
00257
00258
00259
00260 #ifdef MHD
00261 for (k=ks; k<=ke; k++) {
00262 for (j=js; j<=je; j++) {
00263 for (i=is; i<=ie+1; i++) {
00264 pGrid->B1i[k][j][i] = (az[k][j+1][i] - az[k][j][i])/pGrid->dx2 -
00265 (ay[k+1][j][i] - ay[k][j][i])/pGrid->dx3;
00266 }}}
00267 for (k=ks; k<=ke; k++) {
00268 for (j=js; j<=je+1; j++) {
00269 for (i=is; i<=ie; i++) {
00270 pGrid->B2i[k][j][i] = (ax[k+1][j][i] - ax[k][j][i])/pGrid->dx3 -
00271 (az[k][j][i+1] - az[k][j][i])/pGrid->dx1;
00272 }}}
00273 if (ke > ks) {
00274 ku = ke+1;
00275 } else {
00276 ku = ke;
00277 }
00278 for (k=ks; k<=ku; k++) {
00279 for (j=js; j<=je; j++) {
00280 for (i=is; i<=ie; i++) {
00281 pGrid->B3i[k][j][i] = (ay[k][j][i+1] - ay[k][j][i])/pGrid->dx1 -
00282 (ax[k][j+1][i] - ax[k][j][i])/pGrid->dx2;
00283 }}}
00284 #endif
00285
00286
00287
00288 #if defined MHD || !defined ISOTHERMAL
00289 for (k=ks; k<=ke; k++) {
00290 for (j=js; j<=je; j++) {
00291 for (i=is; i<=ie; i++) {
00292 #ifdef MHD
00293 pGrid->U[k][j][i].B1c = 0.5*(pGrid->B1i[k][j][i ] +
00294 pGrid->B1i[k][j][i+1]);
00295 pGrid->U[k][j][i].B2c = 0.5*(pGrid->B2i[k][j ][i] +
00296 pGrid->B2i[k][j+1][i]);
00297 if (ke > ks)
00298 pGrid->U[k][j][i].B3c = 0.5*(pGrid->B3i[k ][j][i] +
00299 pGrid->B3i[k+1][j][i]);
00300 else
00301 pGrid->U[k][j][i].B3c = pGrid->B3i[k][j][i];
00302 #endif
00303
00304 #ifndef ISOTHERMAL
00305 pGrid->U[k][j][i].E = 1.0/Gamma_1
00306 #ifdef MHD
00307 + 0.5*(SQR(pGrid->U[k][j][i].B1c) + SQR(pGrid->U[k][j][i].B2c)
00308 + SQR(pGrid->U[k][j][i].B3c))
00309 #endif
00310 + 0.5*(SQR(pGrid->U[k][j][i].M1) + SQR(pGrid->U[k][j][i].M2)
00311 + SQR(pGrid->U[k][j][i].M3))/pGrid->U[k][j][i].d;
00312 #endif
00313 }
00314 }
00315 }
00316 #endif
00317
00318 free_3d_array((void***)az);
00319 free_3d_array((void***)ay);
00320 free_3d_array((void***)ax);
00321 }
00322
00323
00324
00325
00326
00327
00328
00329
00330
00331
00332
00333
00334
00335
00336
00337
00338 void problem_write_restart(MeshS *pM, FILE *fp)
00339 {
00340 return;
00341 }
00342
00343 void problem_read_restart(MeshS *pM, FILE *fp)
00344 {
00345 return;
00346 }
00347
00348 #ifdef MHD
00349
00350
00351
00352
00353 static Real current(const GridS *pG, const int i, const int j, const int k)
00354 {
00355 return ((pG->B2i[k][j][i]-pG->B2i[k][j][i-1])/pG->dx1 -
00356 (pG->B1i[k][j][i]-pG->B1i[k][j-1][i])/pG->dx2);
00357 }
00358
00359
00360
00361 static Real Bp2(const GridS *pG, const int i, const int j, const int k)
00362 {
00363 return (pG->U[k][j][i].B1c*pG->U[k][j][i].B1c +
00364 pG->U[k][j][i].B2c*pG->U[k][j][i].B2c);
00365 }
00366
00367
00368
00369 static Real divB(const GridS *pG, const int i, const int j, const int k)
00370 {
00371 Real qa;
00372 if (pG->Nx[2] > 1) {
00373 qa = (pG->B1i[k][j][i+1]-pG->B1i[k][j][i])/pG->dx1 +
00374 (pG->B2i[k][j+1][i]-pG->B2i[k][j][i])/pG->dx2 +
00375 (pG->B3i[k+1][j][i]-pG->B3i[k][j][i])/pG->dx3;
00376 } else {
00377 qa = (pG->B1i[k][j][i+1]-pG->B1i[k][j][i])/pG->dx1 +
00378 (pG->B2i[k][j+1][i]-pG->B2i[k][j][i])/pG->dx2;
00379 }
00380 return qa;
00381 }
00382 #endif
00383
00384 #if (NSCALARS > 0)
00385
00386
00387
00388
00389
00390
00391
00392
00393
00394
00395
00396
00397 static Real Temperature(const GridS *pG, const int i, const int j, const int k)
00398 {
00399 Real Temp;
00400 Temp = pG->U[k][j][i].E - (0.5/pG->U[k][j][i].d)*(
00401 SQR(pG->U[k][j][i].M1) + SQR(pG->U[k][j][i].M2) + SQR(pG->U[k][j][i].M3));
00402 #ifdef MHD
00403 Temp -= 0.5*(SQR(pG->U[k][j][i].B1c) + SQR(pG->U[k][j][i].B2c)
00404 + SQR(pG->U[k][j][i].B3c));
00405 #endif
00406 Temp *= Gamma_1/pG->U[k][j][i].d;
00407 return Temp;
00408 }
00409 #endif
00410
00411 ConsFun_t get_usr_expr(const char *expr)
00412 {
00413 #ifdef MHD
00414 if(strcmp(expr,"J3")==0) return current;
00415 else if(strcmp(expr,"Bp2")==0) return Bp2;
00416 else if(strcmp(expr,"divB")==0) return divB;
00417 #endif
00418 #if (NSCALARS > 0)
00419 if(strcmp(expr,"color")==0) return color;
00420 #endif
00421 if(strcmp(expr,"Temperature")==0) return Temperature;
00422 return NULL;
00423 }
00424
00425 VOutFun_t get_usr_out_fun(const char *name){
00426 return NULL;
00427 }
00428
00429 void Userwork_in_loop(MeshS *pM)
00430 {
00431 }
00432
00433 void Userwork_after_loop(MeshS *pM)
00434 {
00435 }