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