00001 #include "copyright.h"
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019 #include <stdio.h>
00020 #include <stdlib.h>
00021 #include <string.h>
00022 #include <stdarg.h>
00023 #include <math.h>
00024 #include "defs.h"
00025 #include "athena.h"
00026 #include "prototypes.h"
00027
00028
00029
00030
00031
00032
00033
00034
00035 char *ath_strdup(const char *in)
00036 {
00037 char *out = (char *)malloc((1+strlen(in))*sizeof(char));
00038 if(out == NULL) {
00039 ath_perr(-1,"ath_strdup: failed to alloc %d\n",(int)(1+strlen(in)));
00040 return NULL;
00041 }
00042 return strcpy(out,in);
00043 }
00044
00045
00046
00047
00048
00049
00050 int ath_gcd(int a, int b)
00051 {
00052 int c;
00053 if(b>a) {c=a; a=b; b=c;}
00054 while((c=a%b)) {a=b; b=c;}
00055 return b;
00056 }
00057
00058
00059
00060
00061
00062
00063
00064 int ath_big_endian(void)
00065 {
00066 short int n = 1;
00067 char *ep = (char *)&n;
00068
00069 return (*ep == 0);
00070 }
00071
00072
00073
00074
00075
00076
00077 void ath_bswap(void *vdat, int len, int cnt)
00078 {
00079 char tmp, *dat = (char *) vdat;
00080 int k;
00081
00082 if (len==1)
00083 return;
00084 else if (len==2)
00085 while (cnt--) {
00086 tmp = dat[0]; dat[0] = dat[1]; dat[1] = tmp;
00087 dat += 2;
00088 }
00089 else if (len==4)
00090 while (cnt--) {
00091 tmp = dat[0]; dat[0] = dat[3]; dat[3] = tmp;
00092 tmp = dat[1]; dat[1] = dat[2]; dat[2] = tmp;
00093 dat += 4;
00094 }
00095 else if (len==8)
00096 while (cnt--) {
00097 tmp = dat[0]; dat[0] = dat[7]; dat[7] = tmp;
00098 tmp = dat[1]; dat[1] = dat[6]; dat[6] = tmp;
00099 tmp = dat[2]; dat[2] = dat[5]; dat[5] = tmp;
00100 tmp = dat[3]; dat[3] = dat[4]; dat[4] = tmp;
00101 dat += 8;
00102 }
00103 else {
00104 for(k=0; k<len/2; k++) {
00105 tmp = dat[k];
00106 dat[k] = dat[len-1-k];
00107 dat[len-1-k] = tmp;
00108 }
00109 }
00110 }
00111
00112
00113
00114
00115
00116
00117
00118 void ath_error(char *fmt, ...)
00119 {
00120 va_list ap;
00121 FILE *atherr = atherr_fp();
00122
00123 fprintf(atherr,"### Fatal error: ");
00124 va_start(ap, fmt);
00125 vfprintf(atherr, fmt, ap);
00126 fflush(atherr);
00127 va_end(ap);
00128
00129 #ifdef MPI_PARALLEL
00130 MPI_Abort(MPI_COMM_WORLD, 1);
00131 #endif
00132
00133 exit(EXIT_FAILURE);
00134 }
00135
00136
00137
00138
00139
00140
00141
00142 void minmax1(Real *data, int nx1, Real *dmino, Real *dmaxo)
00143 {
00144 int i;
00145 register Real dmin, dmax;
00146
00147 dmin = dmax = data[0];
00148 for (i=0; i<nx1; i++) {
00149 dmin = MIN(dmin,data[i]);
00150 dmax = MAX(dmax,data[i]);
00151 }
00152 *dmino = dmin;
00153 *dmaxo = dmax;
00154 }
00155
00156
00157
00158
00159
00160 void minmax2(Real **data, int nx2, int nx1, Real *dmino, Real *dmaxo)
00161 {
00162 int i,j;
00163 register Real dmin, dmax;
00164
00165 dmin = dmax = data[0][0];
00166 for (j=0; j<nx2; j++) {
00167 for (i=0; i<nx1; i++) {
00168 dmin = MIN(dmin,data[j][i]);
00169 dmax = MAX(dmax,data[j][i]);
00170 }
00171 }
00172 *dmino = dmin;
00173 *dmaxo = dmax;
00174 }
00175
00176
00177
00178
00179
00180
00181 void minmax3(Real ***data, int nx3, int nx2, int nx1, Real *dmino, Real *dmaxo)
00182 {
00183 int i,j,k;
00184 register Real dmin, dmax;
00185
00186 dmin = dmax = data[0][0][0];
00187 for (k=0; k<nx3; k++) {
00188 for (j=0; j<nx2; j++) {
00189 for (i=0; i<nx1; i++) {
00190 dmin = MIN(dmin,data[k][j][i]);
00191 dmax = MAX(dmax,data[k][j][i]);
00192 }
00193 }
00194 }
00195 *dmino = dmin;
00196 *dmaxo = dmax;
00197 }
00198
00199
00200
00201
00202
00203
00204
00205 void do_nothing_bc(GridS *pG)
00206 {
00207 }
00208
00209
00210
00211
00212
00213
00214
00215
00216
00217
00218 Real compute_div_b(GridS *pG)
00219 {
00220 #ifdef MHD
00221 int i,j,k,is,ie,js,je,ks,ke;
00222 Real x1,x2,x3,divB,maxdivB=0.0;
00223 Real lsf=1.0,rsf=1.0,dx2=pG->dx2;
00224
00225 is = pG->is; ie = pG->ie;
00226 js = pG->js; je = pG->je;
00227 ks = pG->ks; ke = pG->ke;
00228
00229 for (k=ks; k<=ke; k++) {
00230 for (j=js; j<=je; j++) {
00231 for (i=is; i<=ie; i++) {
00232 cc_pos(pG,i,j,k,&x1,&x2,&x3);
00233 #ifdef CYLINDRICAL
00234 rsf = (x1+0.5*pG->dx1)/x1; lsf = (x1-0.5*pG->dx1)/x1;
00235 dx2 = x1*pG->dx2;
00236 #endif
00237 divB = (rsf*pG->B1i[k][j][i+1] - lsf*pG->B1i[k][j][i])/pG->dx1;
00238 if (je > js)
00239 divB += (pG->B2i[k][j+1][i] - pG->B2i[k][j][i])/dx2;
00240 if (ke > ks)
00241 divB += (pG->B3i[k+1][j][i] - pG->B3i[k][j][i])/pG->dx3;
00242
00243 maxdivB = MAX(maxdivB,fabs(divB));
00244 }
00245 }
00246 }
00247
00248 return maxdivB;
00249
00250 #else
00251 fprintf(stderr,"[compute_div_b]: This only works for MHD!\n");
00252 exit(EXIT_FAILURE);
00253 return 0.0;
00254 #endif
00255 }
00256
00257
00258
00259
00260
00261
00262
00263
00264
00265
00266
00267
00268
00269 void compute_l1_error(const char *problem, const MeshS *pM,
00270 const ConsS ***RootSoln, const int errortype)
00271 {
00272 DomainS *pD=&(pM->Domain[0][0]);
00273 GridS *pG=pM->Domain[0][0].Grid;
00274 int i=0,j=0,k=0;
00275 #if (NSCALARS > 0)
00276 int n;
00277 #endif
00278 int is,ie,js,je,ks,ke;
00279 Real rms_error=0.0;
00280 Real x1,x2,x3,dVol,totVol;
00281 ConsS error,total_error;
00282 FILE *fp;
00283 char *fname, fnamestr[256];
00284 int Nx1,Nx2,Nx3;
00285 #if defined MPI_PARALLEL
00286 double err[8+NSCALARS], tot_err[8+NSCALARS];
00287 int mpi_err;
00288 #endif
00289
00290
00291 memset(&total_error,0.0,sizeof(ConsS));
00292 if (pG == NULL) return;
00293
00294
00295
00296 is = pG->is; ie = pG->ie;
00297 js = pG->js; je = pG->je;
00298 ks = pG->ks; ke = pG->ke;
00299
00300 for (k=ks; k<=ke; k++) {
00301 for (j=js; j<=je; j++) {
00302 memset(&error,0.0,sizeof(ConsS));
00303 for (i=is; i<=ie; i++) {
00304 dVol = 1.0;
00305 if (pG->dx1 > 0.0) dVol *= pG->dx1;
00306 if (pG->dx2 > 0.0) dVol *= pG->dx2;
00307 if (pG->dx3 > 0.0) dVol *= pG->dx3;
00308 #ifdef CYLINDRICAL
00309 cc_pos(pG,i,j,k,&x1,&x2,&x3);
00310 dVol *= x1;
00311 #endif
00312
00313
00314 error.d += dVol*fabs(pG->U[k][j][i].d - RootSoln[k][j][i].d );
00315 error.M1 += dVol*fabs(pG->U[k][j][i].M1 - RootSoln[k][j][i].M1);
00316 error.M2 += dVol*fabs(pG->U[k][j][i].M2 - RootSoln[k][j][i].M2);
00317 error.M3 += dVol*fabs(pG->U[k][j][i].M3 - RootSoln[k][j][i].M3);
00318 #ifdef MHD
00319 error.B1c += dVol*fabs(pG->U[k][j][i].B1c - RootSoln[k][j][i].B1c);
00320 error.B2c += dVol*fabs(pG->U[k][j][i].B2c - RootSoln[k][j][i].B2c);
00321 error.B3c += dVol*fabs(pG->U[k][j][i].B3c - RootSoln[k][j][i].B3c);
00322 #endif
00323 #ifndef ISOTHERMAL
00324 error.E += dVol*fabs(pG->U[k][j][i].E - RootSoln[k][j][i].E );
00325 #endif
00326 #if (NSCALARS > 0)
00327 for (n=0; n<NSCALARS; n++)
00328 error.s[n] += dVol*fabs(pG->U[k][j][i].s[n] - RootSoln[k][j][i].s[n]);;
00329 #endif
00330 }
00331
00332
00333 total_error.d += error.d;
00334 total_error.M1 += error.M1;
00335 total_error.M2 += error.M2;
00336 total_error.M3 += error.M3;
00337 #ifdef MHD
00338 total_error.B1c += error.B1c;
00339 total_error.B2c += error.B2c;
00340 total_error.B3c += error.B3c;
00341 #endif
00342 #ifndef ISOTHERMAL
00343 total_error.E += error.E;
00344 #endif
00345 #if (NSCALARS > 0)
00346 for (n=0; n<NSCALARS; n++) total_error.s[n] += error.s[n];
00347 #endif
00348 }
00349 }
00350
00351 #ifdef MPI_PARALLEL
00352
00353
00354
00355 err[0] = total_error.d;
00356 err[1] = total_error.M1;
00357 err[2] = total_error.M2;
00358 err[3] = total_error.M3;
00359 #ifdef MHD
00360 err[4] = total_error.B1c;
00361 err[5] = total_error.B2c;
00362 err[6] = total_error.B3c;
00363 #endif
00364 #ifndef ISOTHERMAL
00365 err[7] = total_error.E;
00366 #endif
00367 #if (NSCALARS > 0)
00368 for (n=0; n<NSCALARS; n++) err[8+n] = total_error.s[n];
00369 #endif
00370
00371
00372 mpi_err = MPI_Reduce(err, tot_err, (8+NSCALARS),
00373 MPI_DOUBLE, MPI_SUM, 0, pD->Comm_Domain);
00374 if(mpi_err)
00375 ath_error("[compute_l1_error]: MPI_Reduce call returned error = %d\n",
00376 mpi_err);
00377
00378
00379 if(pD->DomNumber == 0){
00380 total_error.d = tot_err[0];
00381 total_error.M1 = tot_err[1];
00382 total_error.M2 = tot_err[2];
00383 total_error.M3 = tot_err[3];
00384 #ifdef MHD
00385 total_error.B1c = tot_err[4];
00386 total_error.B2c = tot_err[5];
00387 total_error.B3c = tot_err[6];
00388 #endif
00389 #ifndef ISOTHERMAL
00390 total_error.E = tot_err[7];
00391 #endif
00392 #if (NSCALARS > 0)
00393 for (n=0; n<NSCALARS; n++) total_error.s[n] = err[8+n];
00394 #endif
00395
00396 }
00397 else return;
00398 #endif
00399
00400
00401 Nx1 = pD->Nx[0];
00402 Nx2 = pD->Nx[1];
00403 Nx3 = pD->Nx[2];
00404
00405 totVol = 1.0;
00406 if (errortype == 1) {
00407 if (pD->MaxX[0] > pD->MinX[0]) totVol *= pD->MaxX[0] - pD->MinX[0];
00408 if (pD->MaxX[1] > pD->MinX[1]) totVol *= pD->MaxX[1] - pD->MinX[1];
00409 if (pD->MaxX[2] > pD->MinX[2]) totVol *= pD->MaxX[2] - pD->MinX[2];
00410 #ifdef CYLINDRICAL
00411 totVol *= 0.5*(pD->MinX[0] + pD->MaxX[0]);
00412 #endif
00413 }
00414
00415
00416
00417
00418 rms_error = SQR(total_error.d) + SQR(total_error.M1) + SQR(total_error.M2)
00419 + SQR(total_error.M3);
00420 #ifdef MHD
00421 rms_error += SQR(total_error.B1c) + SQR(total_error.B2c)
00422 + SQR(total_error.B3c);
00423 #endif
00424 #ifndef ISOTHERMAL
00425 rms_error += SQR(total_error.E);
00426 #endif
00427
00428 rms_error = sqrt(rms_error)/totVol;
00429
00430
00431 sprintf(fnamestr,"%s-errors",problem);
00432 fname = ath_fname(NULL,fnamestr,NULL,NULL,1,0,NULL,"dat");
00433
00434
00435 if((fp=fopen(fname,"r")) != NULL){
00436 if((fp = freopen(fname,"a",fp)) == NULL){
00437 ath_error("[compute_l1_error]: Unable to reopen file.\n");
00438 free(fname);
00439 return;
00440 }
00441 }
00442
00443 else{
00444 if((fp = fopen(fname,"w")) == NULL){
00445 ath_error("[compute_l1_error]: Unable to open file.\n");
00446 free(fname);
00447 return;
00448 }
00449
00450 fprintf(fp,"# Nx1 Nx2 Nx3 RMS-Error d M1 M2 M3");
00451 #ifndef ISOTHERMAL
00452 fprintf(fp," E");
00453 #endif
00454 #ifdef MHD
00455 fprintf(fp," B1c B2c B3c");
00456 #endif
00457 #if (NSCALARS > 0)
00458 for (n=0; n<NSCALARS; n++) {
00459 fprintf(fp," S[ %d ]",n);
00460 }
00461 #endif
00462 fprintf(fp,"\n#\n");
00463 }
00464
00465 fprintf(fp,"%d %d %d %e",Nx1,Nx2,Nx3,rms_error);
00466
00467 fprintf(fp," %e %e %e %e",
00468 (total_error.d /totVol),
00469 (total_error.M1/totVol),
00470 (total_error.M2/totVol),
00471 (total_error.M3/totVol));
00472
00473 #ifndef ISOTHERMAL
00474 fprintf(fp," %e",total_error.E/totVol);
00475 #endif
00476
00477 #ifdef MHD
00478 fprintf(fp," %e %e %e",
00479 (total_error.B1c/totVol),
00480 (total_error.B2c/totVol),
00481 (total_error.B3c/totVol));
00482 #endif
00483 #if (NSCALARS > 0)
00484 for (n=0; n<NSCALARS; n++) {
00485 fprintf(fp," %e",total_error.s[n]/totVol);
00486 }
00487 #endif
00488
00489 fprintf(fp,"\n");
00490
00491 fclose(fp);
00492 free(fname);
00493
00494 return;
00495 }
00496
00497
00498
00499
00500
00501
00502
00503
00504
00505
00506
00507
00508
00509
00510
00511 int sign_change(Real (*func)(const Real,const Real), const Real a0,
00512 const Real b0, const Real x, Real *a, Real *b) {
00513 const int kmax=20;
00514 int k, n, i;
00515 Real delta, fk, fkp1;
00516
00517 for (k=1; k<=kmax; k++) {
00518 n = pow(2,k);
00519 delta = (b0-a0)/(n-1);
00520 *a = a0;
00521 fk = func(x,*a);
00522 for (i=1; i<n; i++) {
00523 *b = *a + delta;
00524 fkp1 = func(x,*b);
00525 if (fkp1*fk < 0)
00526 return 1;
00527 *a = *b;
00528 fk = fkp1;
00529 }
00530 }
00531
00532 return 0;
00533 }
00534
00535
00536
00537
00538
00539
00540
00541 int bisection(Real (*func)(const Real,const Real), const Real a0, const Real b0,
00542 const Real x, Real *root)
00543 {
00544 const Real tol = 1.0E-10;
00545 const int maxiter = 400;
00546 Real a=a0, b=b0, c, fa, fb, fc;
00547 int i;
00548
00549 fa = func(x,a);
00550 fb = func(x,b);
00551 if (fabs(fa) < tol) {
00552 *root = a;
00553 return 1;
00554 }
00555 if (fabs(fb) < tol) {
00556 *root = b;
00557 return 1;
00558 }
00559
00560
00561 for (i = 0; i < maxiter; i++) {
00562 c = 0.5*(a+b);
00563
00564 #ifdef MYDEBUG
00565 printf("c = %f\n", c);
00566 #endif
00567 if (fabs((b-a)/c) < tol) {
00568 #ifdef MYDEBUG
00569 printf("Bisection converged within tolerance of %f!\n", eps);
00570 #endif
00571 *root = c;
00572 return 1;
00573 }
00574 fc = func(x,c);
00575 if (fa*fc < 0) {
00576 b = c;
00577 fb = fc;
00578 }
00579 else if (fc*fb < 0) {
00580 a = c;
00581 fa = fc;
00582 }
00583 else if (fc == 0) {
00584 *root = c;
00585 return 1;
00586 }
00587 else {
00588 ath_error("[bisection]: There is no single root in (%f,%f) for x = %13.10f!!\n", a, b,x);
00589 *root = c;
00590 return 0;
00591 }
00592 }
00593
00594 ath_error("[bisection]: Bisection did not converge in %d iterations for x = %13.10f!!\n", maxiter,x);
00595 *root = c;
00596 return 0;
00597 }
00598
00599
00600
00601
00602
00603
00604
00605
00606
00607
00608
00609
00610
00611
00612
00613
00614
00615
00616
00617 Real trapzd(Real (*func)(Real), const Real a, const Real b, const int n,
00618 const Real s)
00619 {
00620 Real x,tnm,sum,dx;
00621 int it,j;
00622
00623 if (n == 1) {
00624 return 0.5*(b-a)*(func(a)+func(b));
00625 }
00626 else {
00627 for (it=1,j=1; j<n-1; j++) it <<= 1;
00628 tnm = it;
00629 dx = (b-a)/tnm;
00630 x = a + 0.5*dx;
00631 for (sum=0.0,j=1; j<=it; j++,x+=dx) sum += func(x);
00632 return 0.5*(s+(b-a)*sum/tnm);
00633 }
00634 }
00635
00636
00637
00638
00639
00640
00641
00642
00643
00644
00645
00646
00647 #define EPS 1.0e-8
00648 #define JMAX 20
00649
00650 Real qsimp(Real (*func)(Real), const Real a, const Real b)
00651 {
00652 int j;
00653 Real s,st,ost,os;
00654
00655 ost = os = -1.0e30;
00656 for (j=1; j<=JMAX; j++) {
00657 st = trapzd(func,a,b,j,ost);
00658 s = (4.0*st-ost)/3.0;
00659 if (j > 5)
00660 if (fabs(s-os) < EPS*fabs(os) || (s == 0.0 && os == 0.0)) return s;
00661 os=s;
00662 ost=st;
00663 }
00664
00665 ath_error("[qsimp]: Too many steps!\n");
00666 return 0.0;
00667 }
00668
00669
00670
00671
00672
00673
00674
00675
00676
00677 static Real xsav,ysav,zsav,xmin,xmax,ymin,ymax,zmin,zmax;
00678 static Real (*nrfunc)(Real,Real,Real);
00679
00680
00681
00682
00683
00684
00685
00686
00687
00688 Real avg1d(Real (*func)(Real, Real, Real), const GridS *pG,
00689 const int i, const int j, const int k)
00690 {
00691 Real x1,x2,x3,dvol=pG->dx1;
00692 Real fx(Real x);
00693
00694 nrfunc=func;
00695 cc_pos(pG,i,j,k,&x1,&x2,&x3);
00696 xmin = x1 - 0.5*pG->dx1; xmax = x1 + 0.5*pG->dx1;
00697
00698 ysav = x2;
00699 zsav = x3;
00700 #ifdef CYLINDRICAL
00701 dvol *= x1;
00702 #endif
00703
00704 return qsimp(fx,xmin,xmax)/dvol;
00705 }
00706
00707
00708
00709
00710
00711
00712
00713
00714
00715 Real avg2d(Real (*func)(Real, Real, Real), const GridS *pG,
00716 const int i, const int j, const int k)
00717 {
00718 Real x1,x2,x3,dvol=pG->dx1*pG->dx2;
00719 Real fy(Real y);
00720
00721 nrfunc=func;
00722 cc_pos(pG,i,j,k,&x1,&x2,&x3);
00723 xmin = x1 - 0.5*pG->dx1; xmax = x1 + 0.5*pG->dx1;
00724 ymin = x2 - 0.5*pG->dx2; ymax = x2 + 0.5*pG->dx2;
00725
00726 zsav = x3;
00727 #ifdef CYLINDRICAL
00728 dvol *= x1;
00729 #endif
00730
00731 return qsimp(fy,ymin,ymax)/dvol;
00732 }
00733
00734
00735
00736
00737
00738
00739
00740
00741
00742 Real avg3d(Real (*func)(Real, Real, Real), const GridS *pG,
00743 const int i, const int j, const int k)
00744 {
00745 Real x1,x2,x3,dvol=pG->dx1*pG->dx2*pG->dx3;
00746 Real fz(Real z);
00747
00748 nrfunc=func;
00749 cc_pos(pG,i,j,k,&x1,&x2,&x3);
00750 xmin = x1 - 0.5*pG->dx1; xmax = x1 + 0.5*pG->dx1;
00751 ymin = x2 - 0.5*pG->dx2; ymax = x2 + 0.5*pG->dx2;
00752 zmin = x3 - 0.5*pG->dx3; zmax = x3 + 0.5*pG->dx3;
00753
00754 #ifdef CYLINDRICAL
00755 dvol *= x1;
00756 #endif
00757
00758 return qsimp(fz,zmin,zmax)/dvol;
00759 }
00760
00761 Real avgXZ(Real (*func)(Real, Real, Real), const GridS *pG, const int i, const int j, const int k) {
00762 Real x1,x2,x3;
00763
00764 Real fXZ(Real z);
00765
00766 nrfunc=func;
00767 cc_pos(pG,i,j,k,&x1,&x2,&x3);
00768 xmin = x1 - 0.5*pG->dx1; xmax = x1 + 0.5*pG->dx1;
00769 zmin = x3 - 0.5*pG->dx3; zmax = x3 + 0.5*pG->dx3;
00770
00771 ysav = x2;
00772 return qsimp(fXZ,zmin,zmax)/(x1*pG->dx1*pG->dx3);
00773
00774 }
00775
00776 Real fz(Real z)
00777 {
00778 Real fy(Real y);
00779
00780 zsav = z;
00781 return qsimp(fy,ymin,ymax);
00782 }
00783
00784 Real fy(Real y)
00785 {
00786 Real fx(Real x);
00787
00788 ysav = y;
00789 return qsimp(fx,xmin,xmax);
00790 }
00791
00792 Real fx(Real x)
00793 {
00794 #ifdef CYLINDRICAL
00795 return x*nrfunc(x,ysav,zsav);
00796 #else
00797 return nrfunc(x,ysav,zsav);
00798 #endif
00799 }
00800
00801 Real fXZ(Real z) {
00802 Real fx(Real x);
00803
00804 zsav = z;
00805 return qsimp(fx,xmin,xmax);
00806
00807 }
00808
00809
00810
00811
00812
00813
00814
00815
00816
00817
00818 static Real (*a1func)(Real,Real,Real);
00819 static Real (*a2func)(Real,Real,Real);
00820 static Real (*a3func)(Real,Real,Real);
00821
00822
00823
00824
00825
00826
00827
00828
00829
00830
00831 Real vecpot2b1i(Real (*A2)(Real,Real,Real), Real (*A3)(Real,Real,Real),
00832 const GridS *pG, const int i, const int j, const int k)
00833 {
00834 Real x1,x2,x3,b1i=0.0,lsf=1.0,rsf=1.0,dx2=pG->dx2;
00835 Real f2(Real y);
00836 Real f3(Real z);
00837
00838 a2func = A2;
00839 a3func = A3;
00840 cc_pos(pG,i,j,k,&x1,&x2,&x3);
00841 xmin = x1 - 0.5*pG->dx1; xmax = x1 + 0.5*pG->dx1;
00842 ymin = x2 - 0.5*pG->dx2; ymax = x2 + 0.5*pG->dx2;
00843 zmin = x3 - 0.5*pG->dx3; zmax = x3 + 0.5*pG->dx3;
00844
00845 xsav = xmin;
00846 #ifdef CYLINDRICAL
00847 lsf = xmin; rsf = xmin;
00848 dx2 = xmin*pG->dx2;
00849 #endif
00850
00851 if (A2 != NULL) {
00852 if (ymin == ymax)
00853 b1i += rsf*A2(xmin,ymin,zmin) - lsf*A2(xmin,ymin,zmax);
00854 else {
00855 zsav = zmin;
00856 b1i += rsf*qsimp(f2,ymin,ymax);
00857 zsav = zmax;
00858 b1i -= lsf*qsimp(f2,ymin,ymax);
00859 }
00860 }
00861 if (A3 != NULL) {
00862 if (zmin == zmax)
00863 b1i += A3(xmin,ymax,zmin) - A3(xmin,ymin,zmin);
00864 else {
00865 ysav = ymax;
00866 b1i += qsimp(f3,zmin,zmax);
00867 ysav = ymin;
00868 b1i -= qsimp(f3,zmin,zmax);
00869 }
00870 }
00871
00872 if (pG->dx2 > 0.0) b1i /= dx2;
00873 if (pG->dx3 > 0.0) b1i /= pG->dx3;
00874
00875 return b1i;
00876 }
00877
00878
00879
00880
00881
00882
00883
00884
00885
00886
00887 Real vecpot2b2i(Real (*A1)(Real,Real,Real), Real (*A3)(Real,Real,Real),
00888 const GridS *pG, const int i, const int j, const int k)
00889 {
00890 Real x1,x2,x3,b2i=0.0;
00891 Real f1(Real x);
00892 Real f3(Real z);
00893
00894 a1func = A1;
00895 a3func = A3;
00896 cc_pos(pG,i,j,k,&x1,&x2,&x3);
00897 xmin = x1 - 0.5*pG->dx1; xmax = x1 + 0.5*pG->dx1;
00898 ymin = x2 - 0.5*pG->dx2; ymax = x2 + 0.5*pG->dx2;
00899 zmin = x3 - 0.5*pG->dx3; zmax = x3 + 0.5*pG->dx3;
00900
00901 ysav = ymin;
00902
00903 if (A1 != NULL) {
00904 if (xmin == xmax)
00905 b2i += A1(xmin,ymin,zmax) - A1(xmin,ymin,zmin);
00906 else {
00907 zsav = zmax;
00908 b2i += qsimp(f1,xmin,xmax);
00909 zsav = zmin;
00910 b2i -= qsimp(f1,xmin,xmax);
00911 }
00912 }
00913 if (A3 != NULL) {
00914 if (zmin == zmax)
00915 b2i += A3(xmin,ymin,zmin) - A3(xmax,ymin,zmin);
00916 else {
00917 xsav = xmin;
00918 b2i += qsimp(f3,zmin,zmax);
00919 xsav = xmax;
00920 b2i -= qsimp(f3,zmin,zmax);
00921 }
00922 }
00923
00924 if (pG->dx1 > 0.0) b2i /= pG->dx1;
00925 if (pG->dx3 > 0.0) b2i /= pG->dx3;
00926
00927 return b2i;
00928 }
00929
00930
00931
00932
00933
00934
00935
00936
00937
00938
00939 Real vecpot2b3i(Real (*A1)(Real,Real,Real), Real (*A2)(Real,Real,Real),
00940 const GridS *pG, const int i, const int j, const int k)
00941 {
00942 Real x1,x2,x3,b3i=0.0,lsf=1.0,rsf=1.0,dx2=pG->dx2;
00943 Real f1(Real x);
00944 Real f2(Real y);
00945
00946 a1func = A1;
00947 a2func = A2;
00948 cc_pos(pG,i,j,k,&x1,&x2,&x3);
00949 xmin = x1 - 0.5*pG->dx1; xmax = x1 + 0.5*pG->dx1;
00950 ymin = x2 - 0.5*pG->dx2; ymax = x2 + 0.5*pG->dx2;
00951 zmin = x3 - 0.5*pG->dx3; zmax = x3 + 0.5*pG->dx3;
00952
00953 zsav = zmin;
00954 #ifdef CYLINDRICAL
00955 rsf = xmax; lsf = xmin;
00956 dx2 = x1*pG->dx2;
00957 #endif
00958
00959 if (A1 != NULL) {
00960 if (xmin == xmax)
00961 b3i += A1(xmin,ymin,zmin) - A1(xmin,ymax,zmin);
00962 else {
00963 ysav = ymin;
00964 b3i += qsimp(f1,xmin,xmax);
00965 ysav = ymax;
00966 b3i -= qsimp(f1,xmin,xmax);
00967 }
00968 }
00969 if (A2 != NULL) {
00970 if (ymin == ymax)
00971 b3i += rsf*A2(xmax,ymin,zmin) - lsf*A2(xmin,ymin,zmin);
00972 else {
00973 xsav = xmax;
00974 b3i += rsf*qsimp(f2,ymin,ymax);
00975 xsav = xmin;
00976 b3i -= lsf*qsimp(f2,ymin,ymax);
00977 }
00978 }
00979
00980 if (pG->dx1 > 0.0) b3i /= pG->dx1;
00981 if (pG->dx2 > 0.0) b3i /= dx2;
00982
00983 return b3i;
00984 }
00985
00986 Real f1(Real x)
00987 {
00988 return a1func(x,ysav,zsav);
00989 }
00990
00991 Real f2(Real y)
00992 {
00993 return a2func(xsav,y,zsav);
00994 }
00995
00996 Real f3(Real z)
00997 {
00998 return a3func(xsav,ysav,z);
00999 }
01000
01001
01002 #if defined(PARTICLES) || defined(CHEMISTRY)
01003
01004
01005
01006
01007
01008
01009
01010
01011
01012 void ludcmp(Real **a, int n, int *indx, Real *d)
01013 {
01014 int i,imax,j,k;
01015 Real big,dum,sum,temp;
01016 Real *rowscale;
01017
01018 rowscale = (Real*)calloc_1d_array(n, sizeof(Real));
01019 *d=1.0;
01020
01021 for (i=0;i<n;i++)
01022 {
01023 big=0.0;
01024 for (j=0;j<n;j++)
01025 if ((temp=fabs(a[i][j])) > big) big=temp;
01026 if (big == 0.0) ath_error("[LUdecomp]:Input matrix is singular!");
01027 rowscale[i]=1.0/big;
01028 }
01029
01030 for (j=0;j<n;j++) {
01031
01032 for (i=0;i<j;i++) {
01033 sum=a[i][j];
01034 for (k=0;k<i;k++) sum -= a[i][k]*a[k][j];
01035 a[i][j]=sum;
01036 }
01037
01038 big=0.0;
01039 for (i=j;i<n;i++) {
01040 sum=a[i][j];
01041 for (k=0;k<j;k++)
01042 sum -= a[i][k]*a[k][j];
01043 a[i][j]=sum;
01044
01045 if ( (dum=rowscale[i]*fabs(sum)) >= big) {
01046 big=dum;
01047 imax=i;
01048 }
01049 }
01050
01051 if (j != imax) {
01052 for (k=0;k<n;k++) {
01053 dum=a[imax][k];
01054 a[imax][k]=a[j][k];
01055 a[j][k]=dum;
01056 }
01057 *d = -(*d);
01058 rowscale[imax]=rowscale[j];
01059 }
01060 indx[j]=imax;
01061
01062 if (a[j][j] == 0.0) a[j][j]=TINY_NUMBER;
01063 dum=1.0/(a[j][j]);
01064 for (i=j+1;i<n;i++) a[i][j] *= dum;
01065 }
01066 free(rowscale);
01067 }
01068
01069
01070
01071
01072
01073
01074
01075
01076
01077 void lubksb(Real **a, int n, int *indx, Real b[])
01078 {
01079 int i,ii=-1,ip,j;
01080 Real sum;
01081
01082 for (i=0;i<n;i++) {
01083 ip=indx[i];
01084 sum=b[ip];
01085 b[ip]=b[i];
01086 if (ii>=0)
01087 for (j=ii;j<=i-1;j++) sum -= a[i][j]*b[j];
01088 else if (sum) ii=i;
01089 b[i]=sum;
01090 }
01091
01092 for (i=n-1;i>=0;i--) {
01093 sum=b[i];
01094 for (j=i+1;j<n;j++) sum -= a[i][j]*b[j];
01095 b[i]=sum/a[i][i];
01096 }
01097 }
01098
01099
01100
01101
01102
01103
01104
01105
01106 void InverseMatrix(Real **a, int n, Real **b)
01107 {
01108 int i,j,*indx;
01109 Real *col,d;
01110
01111 indx = (int*)calloc_1d_array(n, sizeof(int));
01112 col = (Real*)calloc_1d_array(n, sizeof(Real));
01113
01114 ludcmp(a,n,indx,&d);
01115
01116 for (j=0; j<n; j++) {
01117 for (i=0; i<n; i++) col[i]=0.0;
01118 col[j]=1.0;
01119 lubksb(a, n, indx, col);
01120 for (i=0; i<n; i++) b[i][j] = col[i];
01121 }
01122
01123 return;
01124 }
01125
01126
01127
01128
01129 void MatrixMult(Real **a, Real **b, int m, int n, int l, Real **c)
01130 {
01131 int i, j, k;
01132 for (i=0; i<m; i++)
01133 for (j=0; j<l; j++)
01134 {
01135 c[i][j] = 0.0;
01136 for (k=0; k<n; k++) c[i][j] += a[i][k] * b[k][j];
01137 }
01138 }
01139 #endif