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 #include <math.h>
00027 #include <stdio.h>
00028 #include <stdlib.h>
00029 #include "defs.h"
00030 #include "athena.h"
00031 #include "globals.h"
00032 #include "prototypes.h"
00033
00034
00035 static ConsS ***RootSoln=NULL;
00036 static int wave_flag;
00037
00038
00039
00040
00041
00042 void problem(DomainS *pDomain)
00043 {
00044 GridS *pGrid=(pDomain->Grid);
00045 ConsS ***Soln;
00046 int i=0,j=0,k=0;
00047 int is,ie,js,je,ks,ke,n,m,nx1,nx2,nx3,wave_dir;
00048 Real amp,vflow;
00049 Real d0,p0,u0,v0,w0,h0;
00050 Real x1,x2,x3,r,ev[NWAVE],rem[NWAVE][NWAVE],lem[NWAVE][NWAVE];
00051 #ifdef MHD
00052 Real bx0,by0,bz0,xfact,yfact;
00053 #endif
00054 is = pGrid->is; ie = pGrid->ie;
00055 js = pGrid->js; je = pGrid->je;
00056 ks = pGrid->ks; ke = pGrid->ke;
00057 nx1 = (ie-is)+1 + 2*nghost;
00058 nx2 = (je-js)+1 + 2*nghost;
00059 nx3 = (ke-ks)+1 + 2*nghost;
00060
00061
00062
00063 wave_flag = par_geti("problem","wave_flag");
00064 amp = par_getd("problem","amp");
00065 vflow = par_getd("problem","vflow");
00066 wave_dir = par_getd("problem","wave_dir");
00067
00068
00069
00070
00071 if ((Soln = (ConsS***)calloc_3d_array(nx3,nx2,nx1,sizeof(ConsS)))==NULL)
00072 ath_error("[problem]: Error allocating memory for Soln\n");
00073 if (pDomain->Level == 0){
00074 if ((RootSoln = (ConsS***)calloc_3d_array(nx3,nx2,nx1,sizeof(ConsS)))
00075 == NULL) ath_error("[problem]: Error alloc memory for RootSoln\n");
00076 }
00077
00078
00079
00080
00081
00082 d0 = 1.0;
00083 #ifndef ISOTHERMAL
00084 p0 = 1.0/Gamma;
00085 u0 = vflow*sqrt(Gamma*p0/d0);
00086 #else
00087 u0 = vflow*Iso_csound;
00088 #endif
00089 v0 = 0.0;
00090 w0 = 0.0;
00091 #ifdef MHD
00092 bx0 = 1.0;
00093 by0 = sqrt(2.0);
00094 bz0 = 0.5;
00095 xfact = 0.0;
00096 yfact = 1.0;
00097 #endif
00098
00099 for (n=0; n<NWAVE; n++) {
00100 for (m=0; m<NWAVE; m++) {
00101 rem[n][m] = 0.0;
00102 lem[n][m] = 0.0;
00103 }
00104 }
00105
00106 #ifdef HYDRO
00107 #ifdef ISOTHERMAL
00108 esys_roe_iso_hyd(u0,v0,w0, ev,rem,lem);
00109 #else
00110 h0 = ((p0/Gamma_1 + 0.5*d0*(u0*u0+v0*v0+w0*w0)) + p0)/d0;
00111 esys_roe_adb_hyd(u0,v0,w0,h0,ev,rem,lem);
00112 ath_pout(0,"Ux - Cs = %e, %e\n",ev[0],rem[0][wave_flag]);
00113 ath_pout(0,"Ux = %e, %e\n",ev[1],rem[1][wave_flag]);
00114 ath_pout(0,"Ux + Cs = %e, %e\n",ev[4],rem[4][wave_flag]);
00115 #endif
00116 #endif
00117
00118 #ifdef MHD
00119 #if defined(ISOTHERMAL)
00120 esys_roe_iso_mhd(d0,u0,v0,w0,bx0,by0,bz0,xfact,yfact,ev,rem,lem);
00121 ath_pout(0,"Ux - Cf = %e, %e\n",ev[0],rem[0][wave_flag]);
00122 ath_pout(0,"Ux - Ca = %e, %e\n",ev[1],rem[1][wave_flag]);
00123 ath_pout(0,"Ux - Cs = %e, %e\n",ev[2],rem[2][wave_flag]);
00124 ath_pout(0,"Ux + Cs = %e, %e\n",ev[3],rem[3][wave_flag]);
00125 ath_pout(0,"Ux + Ca = %e, %e\n",ev[4],rem[4][wave_flag]);
00126 ath_pout(0,"Ux + Cf = %e, %e\n",ev[5],rem[5][wave_flag]);
00127 #else
00128 h0 = ((p0/Gamma_1+0.5*(bx0*bx0+by0*by0+bz0*bz0)+0.5*d0*(u0*u0+v0*v0+w0*w0))
00129 + (p0+0.5*(bx0*bx0+by0*by0+bz0*bz0)))/d0;
00130 esys_roe_adb_mhd(d0,u0,v0,w0,h0,bx0,by0,bz0,xfact,yfact,ev,rem,lem);
00131 ath_pout(0,"Ux - Cf = %e, %e\n",ev[0],rem[0][wave_flag]);
00132 ath_pout(0,"Ux - Ca = %e, %e\n",ev[1],rem[1][wave_flag]);
00133 ath_pout(0,"Ux - Cs = %e, %e\n",ev[2],rem[2][wave_flag]);
00134 ath_pout(0,"Ux = %e, %e\n",ev[3],rem[3][wave_flag]);
00135 ath_pout(0,"Ux + Cs = %e, %e\n",ev[4],rem[4][wave_flag]);
00136 ath_pout(0,"Ux + Ca = %e, %e\n",ev[5],rem[5][wave_flag]);
00137 ath_pout(0,"Ux + Cf = %e, %e\n",ev[6],rem[6][wave_flag]);
00138 #endif
00139 #endif
00140
00141
00142
00143 for (k=ks; k<=ke; k++) {
00144 for (j=js; j<=je; j++) {
00145 for (i=is; i<=ie; i++) {
00146
00147
00148
00149 cc_pos(pGrid,i,j,k,&x1,&x2,&x3);
00150 Soln[k][j][i].d = d0;
00151 #ifndef ISOTHERMAL
00152 Soln[k][j][i].E = p0/Gamma_1 + 0.5*d0*u0*u0;
00153 #ifdef MHD
00154 Soln[k][j][i].E += 0.5*(bx0*bx0 + by0*by0 + bz0*bz0);
00155 #endif
00156 #endif
00157
00158
00159
00160 switch(wave_dir){
00161
00162
00163
00164 case 1:
00165 Soln[k][j][i].d += amp*sin(2.0*PI*x1)*rem[0][wave_flag];
00166 #ifndef ISOTHERMAL
00167 Soln[k][j][i].E += amp*sin(2.0*PI*x1)*rem[4][wave_flag];
00168 #endif
00169 Soln[k][j][i].M1 = d0*vflow;
00170 Soln[k][j][i].M2 = 0.0;
00171 Soln[k][j][i].M3 = 0.0;
00172 Soln[k][j][i].M1 += amp*sin(2.0*PI*x1)*rem[1][wave_flag];
00173 Soln[k][j][i].M2 += amp*sin(2.0*PI*x1)*rem[2][wave_flag];
00174 Soln[k][j][i].M3 += amp*sin(2.0*PI*x1)*rem[3][wave_flag];
00175 #ifdef MHD
00176 Soln[k][j][i].B1c = bx0;
00177 Soln[k][j][i].B2c = by0 + amp*sin(2.0*PI*x1)*rem[NWAVE-2][wave_flag];
00178 Soln[k][j][i].B3c = bz0 + amp*sin(2.0*PI*x1)*rem[NWAVE-1][wave_flag];
00179 #endif
00180 #if (NSCALARS > 0)
00181 for (n=0; n<NSCALARS; n++)
00182 Soln[k][j][i].s[n] = amp*(1.0 + sin(2.0*PI*x1));
00183 #endif
00184 break;
00185
00186
00187
00188 case 2:
00189 Soln[k][j][i].d += amp*sin(2.0*PI*x2)*rem[0][wave_flag];
00190 #ifndef ISOTHERMAL
00191 Soln[k][j][i].E += amp*sin(2.0*PI*x2)*rem[4][wave_flag];
00192 #endif
00193 Soln[k][j][i].M1 = 0.0;
00194 Soln[k][j][i].M2 = d0*vflow;
00195 Soln[k][j][i].M3 = 0.0;
00196 Soln[k][j][i].M1 += amp*sin(2.0*PI*x2)*rem[3][wave_flag];
00197 Soln[k][j][i].M2 += amp*sin(2.0*PI*x2)*rem[1][wave_flag];
00198 Soln[k][j][i].M3 += amp*sin(2.0*PI*x2)*rem[2][wave_flag];
00199 #ifdef MHD
00200 Soln[k][j][i].B1c = bz0 + amp*sin(2.0*PI*x2)*rem[NWAVE-1][wave_flag];
00201 Soln[k][j][i].B2c = bx0;
00202 Soln[k][j][i].B3c = by0 + amp*sin(2.0*PI*x2)*rem[NWAVE-2][wave_flag];
00203 #endif
00204 #if (NSCALARS > 0)
00205 for (n=0; n<NSCALARS; n++)
00206 Soln[k][j][i].s[n] = amp*(1.0 + sin(2.0*PI*x2));
00207 #endif
00208 break;
00209
00210
00211
00212 case 3:
00213 Soln[k][j][i].d += amp*sin(2.0*PI*x3)*rem[0][wave_flag];
00214 #ifndef ISOTHERMAL
00215 Soln[k][j][i].E += amp*sin(2.0*PI*x3)*rem[4][wave_flag];
00216 #endif
00217 Soln[k][j][i].M1 = 0.0;
00218 Soln[k][j][i].M2 = 0.0;
00219 Soln[k][j][i].M3 = d0*vflow;
00220 Soln[k][j][i].M1 += amp*sin(2.0*PI*x3)*rem[2][wave_flag];
00221 Soln[k][j][i].M2 += amp*sin(2.0*PI*x3)*rem[3][wave_flag];
00222 Soln[k][j][i].M3 += amp*sin(2.0*PI*x3)*rem[1][wave_flag];
00223 #ifdef MHD
00224 Soln[k][j][i].B1c = by0 + amp*sin(2.0*PI*x3)*rem[NWAVE-2][wave_flag];
00225 Soln[k][j][i].B2c = bz0 + amp*sin(2.0*PI*x3)*rem[NWAVE-1][wave_flag];
00226 Soln[k][j][i].B3c = bx0;
00227 #endif
00228 #if (NSCALARS > 0)
00229 for (n=0; n<NSCALARS; n++)
00230 Soln[k][j][i].s[n] = amp*(1.0 + sin(2.0*PI*x3));
00231 #endif
00232 break;
00233 default:
00234 ath_error("[linear_wave1d]: wave direction %d not allowed\n",wave_dir);
00235 }
00236 }}}
00237
00238
00239
00240 for (k=ks; k<=ke; k++) {
00241 for (j=js; j<=je; j++) {
00242 for (i=is; i<=ie; i++) {
00243 pGrid->U[k][j][i].d = Soln[k][j][i].d;
00244 #ifndef ISOTHERMAL
00245 pGrid->U[k][j][i].E = Soln[k][j][i].E;
00246 #endif
00247 pGrid->U[k][j][i].M1 = Soln[k][j][i].M1;
00248 pGrid->U[k][j][i].M2 = Soln[k][j][i].M2;
00249 pGrid->U[k][j][i].M3 = Soln[k][j][i].M3;
00250 #ifdef MHD
00251 pGrid->U[k][j][i].B1c = Soln[k][j][i].B1c;
00252 pGrid->U[k][j][i].B2c = Soln[k][j][i].B2c;
00253 pGrid->U[k][j][i].B3c = Soln[k][j][i].B3c;
00254 pGrid->B1i[k][j][i] = Soln[k][j][i].B1c;
00255 pGrid->B2i[k][j][i] = Soln[k][j][i].B2c;
00256 pGrid->B3i[k][j][i] = Soln[k][j][i].B3c;
00257 #endif
00258 #if (NSCALARS > 0)
00259 for (n=0; n<NSCALARS; n++)
00260 pGrid->U[k][j][i].s[n] = Soln[k][j][i].s[n];
00261 #endif
00262 }}}
00263 #ifdef MHD
00264 for (k=ks; k<=ke; k++) {
00265 for (j=js; j<=je; j++) {
00266 pGrid->B1i[k][j][ie+1] = pGrid->B1i[k][j][ie];
00267 }
00268 }
00269 if (pGrid->Nx[1] > 1) {
00270 for (k=ks; k<=ke; k++) {
00271 for (i=is; i<=ie; i++) {
00272 pGrid->B2i[k][je+1][i] = pGrid->B2i[k][je][i];
00273 }
00274 }
00275 }
00276 if (pGrid->Nx[2] > 1) {
00277 for (j=js; j<=je; j++) {
00278 for (i=is; i<=ie; i++) {
00279 pGrid->B3i[ke+1][j][i] = pGrid->B3i[ke][j][i];
00280 }
00281 }
00282 }
00283 #endif
00284
00285
00286
00287 #ifdef SELF_GRAVITY
00288 four_pi_G = par_getd("problem","four_pi_G");
00289 grav_mean_rho = d0;
00290 #endif
00291
00292
00293
00294 #ifdef RESISTIVITY
00295 eta_Ohm = par_getd("problem","eta_O");
00296 Q_AD = par_getd_def("problem","Q_AD",0.0);
00297 d_ind = par_getd_def("problem","d_ind",0.0);
00298 #endif
00299 #ifdef NAVIER_STOKES
00300 nu_V = par_getd("problem","nu");
00301 #endif
00302
00303
00304
00305 if (pDomain->Level == 0) {
00306 for (k=ks; k<=ke; k++) {
00307 for (j=js; j<=je; j++) {
00308 for (i=is; i<=ie; i++) {
00309 RootSoln[k][j][i].d = Soln[k][j][i].d ;
00310 RootSoln[k][j][i].M1 = Soln[k][j][i].M1;
00311 RootSoln[k][j][i].M2 = Soln[k][j][i].M2;
00312 RootSoln[k][j][i].M3 = Soln[k][j][i].M3;
00313 #ifndef ISOTHERMAL
00314 RootSoln[k][j][i].E = Soln[k][j][i].E ;
00315 #endif
00316 #ifdef MHD
00317 RootSoln[k][j][i].B1c = Soln[k][j][i].B1c;
00318 RootSoln[k][j][i].B2c = Soln[k][j][i].B2c;
00319 RootSoln[k][j][i].B3c = Soln[k][j][i].B3c;
00320 #endif
00321 #if (NSCALARS > 0)
00322 for (n=0; n<NSCALARS; n++)
00323 RootSoln[k][j][i].s[n] = Soln[k][j][i].s[n];
00324 #endif
00325 }}}
00326 }
00327
00328 free_3d_array(Soln);
00329
00330 return;
00331 }
00332
00333
00334
00335
00336
00337
00338
00339
00340
00341
00342
00343
00344 void problem_write_restart(MeshS *pM, FILE *fp)
00345 {
00346 return;
00347 }
00348
00349 void problem_read_restart(MeshS *pM, FILE *fp)
00350 {
00351 #ifdef SELF_GRAVITY
00352 Real d0 = 1.0;
00353 four_pi_G = par_getd("problem","four_pi_G");
00354 grav_mean_rho = d0;
00355 #endif
00356 #ifdef RESISTIVITY
00357 eta_Ohm = par_getd_def("problem","eta_O",0.0);
00358 Q_Hall = par_getd_def("problem","Q_H",0.0);
00359 Q_AD = par_getd_def("problem","Q_AD",0.0);
00360 #endif
00361 #ifdef NAVIER_STOKES
00362 nu_V = par_getd("problem","nu");
00363 #endif
00364 return;
00365 }
00366
00367 ConsFun_t get_usr_expr(const char *expr)
00368 {
00369 return NULL;
00370 }
00371
00372 VOutFun_t get_usr_out_fun(const char *name){
00373 return NULL;
00374 }
00375
00376 #ifdef RESISTIVITY
00377
00378
00379
00380 void get_eta_user(GridS *pG, int i, int j, int k,
00381 Real *eta_O, Real *eta_H, Real *eta_A)
00382 {
00383 *eta_O = 0.0;
00384 *eta_H = 0.0;
00385 *eta_A = 0.0;
00386
00387 return;
00388 }
00389 #endif
00390
00391 void Userwork_in_loop(MeshS *pM)
00392 {
00393 }
00394
00395
00396
00397
00398
00399
00400
00401
00402
00403 void Userwork_after_loop(MeshS *pM)
00404 {
00405 GridS *pGrid;
00406 int i=0,j=0,k=0;
00407 int is,ie,js,je,ks,ke;
00408 #if (NSCALARS > 0)
00409 int n;
00410 #endif
00411
00412 Real rms_error=0.0;
00413 ConsS error,total_error;
00414 FILE *fp;
00415 char *fname;
00416 int Nx1, Nx2, Nx3, count;
00417
00418 total_error.d = 0.0;
00419 total_error.M1 = 0.0;
00420 total_error.M2 = 0.0;
00421 total_error.M3 = 0.0;
00422 #ifdef MHD
00423 total_error.B1c = 0.0;
00424 total_error.B2c = 0.0;
00425 total_error.B3c = 0.0;
00426 #endif
00427 #ifndef ISOTHERMAL
00428 total_error.E = 0.0;
00429 #endif
00430 #if (NSCALARS > 0)
00431 for (n=0; n<NSCALARS; n++) total_error.s[n] = 0.0;
00432 #endif
00433
00434
00435
00436 pGrid=pM->Domain[0][0].Grid;
00437 if (pGrid == NULL) return;
00438
00439
00440
00441 is = pGrid->is; ie = pGrid->ie;
00442 js = pGrid->js; je = pGrid->je;
00443 ks = pGrid->ks; ke = pGrid->ke;
00444 for (k=ks; k<=ke; k++) {
00445 for (j=js; j<=je; j++) {
00446 error.d = 0.0;
00447 error.M1 = 0.0;
00448 error.M2 = 0.0;
00449 error.M3 = 0.0;
00450 #ifdef MHD
00451 error.B1c = 0.0;
00452 error.B2c = 0.0;
00453 error.B3c = 0.0;
00454 #endif
00455 #ifndef ISOTHERMAL
00456 error.E = 0.0;
00457 #endif
00458 #if (NSCALARS > 0)
00459 for (n=0; n<NSCALARS; n++) error.s[n] = 0.0;
00460 #endif
00461
00462 for (i=is; i<=ie; i++) {
00463 error.d += fabs(pGrid->U[k][j][i].d - RootSoln[k][j][i].d );
00464 error.M1 += fabs(pGrid->U[k][j][i].M1 - RootSoln[k][j][i].M1);
00465 error.M2 += fabs(pGrid->U[k][j][i].M2 - RootSoln[k][j][i].M2);
00466 error.M3 += fabs(pGrid->U[k][j][i].M3 - RootSoln[k][j][i].M3);
00467 #ifdef MHD
00468 error.B1c += fabs(pGrid->U[k][j][i].B1c - RootSoln[k][j][i].B1c);
00469 error.B2c += fabs(pGrid->U[k][j][i].B2c - RootSoln[k][j][i].B2c);
00470 error.B3c += fabs(pGrid->U[k][j][i].B3c - RootSoln[k][j][i].B3c);
00471 #endif
00472 #ifndef ISOTHERMAL
00473 error.E += fabs(pGrid->U[k][j][i].E - RootSoln[k][j][i].E );
00474 #endif
00475 #if (NSCALARS > 0)
00476 for (n=0; n<NSCALARS; n++)
00477 error.s[n] += fabs(pGrid->U[k][j][i].s[n] - RootSoln[k][j][i].s[n]);;
00478 #endif
00479 }
00480
00481 total_error.d += error.d;
00482 total_error.M1 += error.M1;
00483 total_error.M2 += error.M2;
00484 total_error.M3 += error.M3;
00485 #ifdef MHD
00486 total_error.B1c += error.B1c;
00487 total_error.B2c += error.B2c;
00488 total_error.B3c += error.B3c;
00489 #endif
00490 #ifndef ISOTHERMAL
00491 total_error.E += error.E;
00492 #endif
00493 #if (NSCALARS > 0)
00494 for (n=0; n<NSCALARS; n++) total_error.s[n] += error.s[n];
00495 #endif
00496 }}
00497
00498 Nx1 = ie - is + 1;
00499 Nx2 = je - js + 1;
00500 Nx3 = ke - ks + 1;
00501
00502 count = Nx1*Nx2*Nx3;
00503
00504
00505
00506 rms_error = SQR(total_error.d) + SQR(total_error.M1) + SQR(total_error.M2)
00507 + SQR(total_error.M3);
00508 #ifdef MHD
00509 rms_error += SQR(total_error.B1c) + SQR(total_error.B2c)
00510 + SQR(total_error.B3c);
00511 #endif
00512 #ifndef ISOTHERMAL
00513 rms_error += SQR(total_error.E);
00514 #endif
00515 rms_error = sqrt(rms_error)/(double)count;
00516
00517
00518
00519 fname = ath_fname(NULL,"LinWave-errors",NULL,NULL,1,wave_flag,NULL,"dat");
00520
00521 if((fp=fopen(fname,"r")) != NULL){
00522 if((fp = freopen(fname,"a",fp)) == NULL){
00523 ath_error("[Userwork_after_loop]: Unable to reopen file.\n");
00524 return;
00525 }
00526 }
00527
00528 else{
00529 if((fp = fopen(fname,"w")) == NULL){
00530 ath_error("[Userwork_after_loop]: Unable to open file.\n");
00531 return;
00532 }
00533
00534 fprintf(fp,"# Nx1 Nx2 Nx3 RMS-Error d M1 M2 M3");
00535 #ifndef ISOTHERMAL
00536 fprintf(fp," E");
00537 #endif
00538 #ifdef MHD
00539 fprintf(fp," B1c B2c B3c");
00540 #endif
00541 #if (NSCALARS > 0)
00542 for (n=0; n<NSCALARS; n++) {
00543 fprintf(fp," S[ %d ]",n);
00544 }
00545 #endif
00546
00547 fprintf(fp,"\n#\n");
00548 }
00549
00550 fprintf(fp,"%d %d %d %e",Nx1,Nx2,Nx3,rms_error);
00551
00552 fprintf(fp," %e %e %e %e",
00553 (total_error.d/(double) count),
00554 (total_error.M1/(double)count),
00555 (total_error.M2/(double)count),
00556 (total_error.M3/(double)count));
00557 #ifndef ISOTHERMAL
00558 fprintf(fp," %e",(total_error.E/(double)count));
00559 #endif
00560 #ifdef MHD
00561 fprintf(fp," %e %e %e",
00562 (total_error.B1c/(double)count),
00563 (total_error.B2c/(double)count),
00564 (total_error.B3c/(double)count));
00565 #endif
00566 #if (NSCALARS > 0)
00567 for (n=0; n<NSCALARS; n++) {
00568 fprintf(fp," %e",total_error.s[n]/(double)count);
00569 }
00570 #endif
00571 fprintf(fp,"\n");
00572
00573 fclose(fp);
00574
00575 return;
00576 }