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 #include <math.h>
00026 #include <stdio.h>
00027 #include <stdlib.h>
00028 #include <float.h>
00029 #include "../defs.h"
00030 #include "../athena.h"
00031 #include "../globals.h"
00032 #include "prototypes.h"
00033 #include "../prototypes.h"
00034
00035 #ifdef EXACT_FLUX
00036 #ifndef SPECIAL_RELATIVITY
00037
00038 #ifdef MHD
00039 #error : The exact flux for MHD has not been implemented.
00040 #endif
00041
00042 #if (NSCALARS > 0)
00043 #error : Passive scalars have not been implemented in the exact flux.
00044 #endif
00045
00046 #ifdef ISOTHERMAL
00047
00048 static void srder(double dm, double vl, double vr, double dmin, double dmax,
00049 double *y, double *dydx);
00050 static double rtsafe(void (*funcd)(double, double, double, double, double,
00051 double *, double *), double x1, double x2, double xacc,
00052 double vl, double vr, double dmin, double dmax);
00053
00054
00055
00056
00057
00058
00059
00060
00061
00062
00063
00064
00065
00066
00067
00068 void fluxes(const Cons1DS Ul, const Cons1DS Ur,
00069 const Prim1DS Wl, const Prim1DS Wr,
00070 const Real Bxi, Cons1DS *pF)
00071 {
00072 Real zl, zr, zm, dm, Vxm, Mxm, tmp, dmin, dmax;
00073 Real sl, sr;
00074 Real hdl, hdr;
00075 Real tll, tlr;
00076 char soln;
00077
00078 if(!(Ul.d > 0.0)||!(Ur.d > 0.0))
00079 ath_error("[exact flux]: Non-positive densities: dl = %e dr = %e\n",
00080 Ul.d, Ur.d);
00081
00082
00083
00084
00085
00086 zl = sqrt((double)Wl.d);
00087 zr = sqrt((double)Wr.d);
00088
00089
00090 soln = 0;
00091
00092
00093
00094 tmp = zl*zr*(Wl.Vx - Wr.Vx)/(2.0*Iso_csound*(zl + zr));
00095 zm = tmp + sqrt((double)(tmp*tmp + zl*zr));
00096 dm = zm*zm;
00097
00098
00099 Vxm = Wl.Vx - Iso_csound*(dm-Wl.d)/(zm*zl);
00100
00101
00102
00103 dmin = MIN(Wl.d, Wr.d);
00104 dmax = MAX(Wl.d, Wr.d);
00105 if (dm < dmax) {
00106
00107 soln = 3;
00108
00109
00110
00111 dm = zl*zr*exp((Wl.Vx-Wr.Vx)/(2.0*Iso_csound));
00112
00113
00114 Vxm = Wl.Vx - Iso_csound*log(dm/Wl.d);
00115
00116
00117
00118 if (dm > dmin) {
00119
00120
00121
00122
00123
00124
00125 if (Wl.d > Wr.d) soln = 2; else soln = 1;
00126
00127 dm = rtsafe(&srder,dmin,dmax, 2.0*DBL_EPSILON, Wl.Vx, Wr.Vx, dmin, dmax);
00128
00129
00130 if ((dm > dmin) && (dm <= dmax)) {
00131 if (Wl.d > Wr.d) {
00132
00133 Vxm = Wl.Vx - Iso_csound*log(dm/Wl.d);
00134 } else {
00135
00136 Vxm = Wr.Vx + Iso_csound*log(dm/Wr.d);
00137 }
00138 } else {
00139
00140 soln = 3;
00141
00142
00143
00144
00145
00146
00147 dm = zl*zr*exp((Wl.Vx-Wr.Vx)/(2.0*Iso_csound));
00148
00149
00150 Vxm = Wl.Vx - Iso_csound*log(dm/Wl.d);
00151 }
00152 }
00153 }
00154
00155 if (dm < 0.0)
00156 ath_error("[exact flux]: Solver finds negative density %5.4e\n", dm);
00157
00158
00159
00160
00161
00162
00163 if (soln & 2) {
00164
00165 hdl = Wl.Vx - Iso_csound;
00166 tll = Vxm - Iso_csound;
00167
00168 if (hdl >= 0.0) {
00169
00170 pF->d = Ul.Mx;
00171 pF->Mx = Ul.Mx*(Wl.Vx) + Wl.d*Iso_csound2;
00172 pF->My = Ul.My*(Wl.Vx);
00173 pF->Mz = Ul.Mz*(Wl.Vx);
00174 return;
00175 } else if (tll >= 0.0) {
00176
00177 dm = Ul.d*exp(hdl/Iso_csound);
00178 Mxm = Ul.d*Iso_csound*exp(hdl/Iso_csound);
00179 Vxm = (dm == 0.0 ? 0.0 : Mxm / dm);
00180
00181 pF->d = Mxm;
00182 pF->Mx = Mxm*Vxm + dm*Iso_csound2;
00183 pF->My = Mxm*Wl.Vy;
00184 pF->Mz = Mxm*Wl.Vz;
00185 return;
00186 }
00187 } else {
00188
00189 sl = Wl.Vx - Iso_csound*sqrt(dm)/zl;
00190
00191 if(sl >= 0.0) {
00192
00193 pF->d = Ul.Mx;
00194 pF->Mx = Ul.Mx*(Wl.Vx) + Wl.d*Iso_csound2;
00195 pF->My = Ul.My*(Wl.Vx);
00196 pF->Mz = Ul.Mz*(Wl.Vx);
00197 return;
00198 }
00199 }
00200
00201 if (soln & 1) {
00202
00203 hdr = Wr.Vx + Iso_csound;
00204 tlr = Vxm + Iso_csound;
00205
00206 if (hdr <= 0.0) {
00207
00208 pF->d = Ur.Mx;
00209 pF->Mx = Ur.Mx*(Wr.Vx) + Wr.d*Iso_csound2;
00210 pF->My = Ur.My*(Wr.Vx);
00211 pF->Mz = Ur.Mz*(Wr.Vx);
00212 return;
00213 } else if (tlr <= 0.0) {
00214
00215 tmp = dm;
00216 dm = tmp*exp(-tlr/Iso_csound);
00217 Mxm = -tmp*Iso_csound*exp(-tlr/Iso_csound);
00218 Vxm = (dm == 0.0 ? 0.0 : Mxm / dm);
00219
00220 pF->d = Mxm;
00221 pF->Mx = Mxm*Vxm + dm*Iso_csound2;
00222 pF->My = Mxm*Wr.Vy;
00223 pF->Mz = Mxm*Wr.Vz;
00224 return;
00225 }
00226 } else {
00227
00228 sr = Wr.Vx + Iso_csound*sqrt(dm)/zr;
00229
00230 if(sr <= 0.0) {
00231
00232 pF->d = Ur.Mx;
00233 pF->Mx = Ur.Mx*(Wr.Vx) + Wr.d*Iso_csound2;
00234 pF->My = Ur.My*(Wr.Vx);
00235 pF->Mz = Ur.Mz*(Wr.Vx);
00236 return;
00237 }
00238 }
00239
00240
00241
00242
00243
00244
00245 if(Vxm >= 0.0){
00246 pF->d = dm*Vxm;
00247 pF->Mx = dm*Vxm*Vxm + dm*Iso_csound2;
00248 pF->My = dm*Vxm*Wl.Vy;
00249 pF->Mz = dm*Vxm*Wl.Vz;
00250 }
00251 else{
00252 pF->d = dm*Vxm;
00253 pF->Mx = dm*Vxm*Vxm + dm*Iso_csound2;
00254 pF->My = dm*Vxm*Wr.Vy;
00255 pF->Mz = dm*Vxm*Wr.Vz;
00256 }
00257
00258 return;
00259 }
00260
00261
00262
00263
00264
00265
00266 static void srder(double dm, double vl, double vr, double dmin, double dmax,
00267 double *y, double *dydx)
00268 {
00269 *y = (vr - vl) + Iso_csound*(log(dm/dmax) + (dm-dmin)/sqrt(dm*dmin));
00270 *dydx = Iso_csound/dm*(1.0 + 0.5*(dm+dmin)/sqrt(dm*dmin));
00271
00272 return;
00273 }
00274
00275
00276
00277
00278
00279
00280
00281
00282 #define MAXIT 100
00283
00284 static double rtsafe(void (*funcd)(double, double, double, double, double,
00285 double *, double *), double x1, double x2, double xacc,
00286 double vl, double vr, double dmin, double dmax)
00287 {
00288 int j;
00289 double df,dx,dxold,f,fh,fl;
00290 double temp,xh,xl,rts;
00291
00292 (*funcd)(x1,vl,vr,dmin,dmax,&fl,&df);
00293 (*funcd)(x2,vl,vr,dmin,dmax,&fh,&df);
00294 if ((fl > 0.0 && fh > 0.0) || (fl < 0.0 && fh < 0.0)) {
00295 return 0.0;
00296 }
00297 if (fl == 0.0) return x1;
00298 if (fh == 0.0) return x2;
00299 if (fl < 0.0) {
00300 xl=x1;
00301 xh=x2;
00302 } else {
00303 xh=x1;
00304 xl=x2;
00305 }
00306 rts=0.5*(x1+x2);
00307 dxold=fabs(x2-x1);
00308 dx=dxold;
00309 (*funcd)(rts,vl,vr,dmin,dmax,&f,&df);
00310 for (j=1;j<=MAXIT;j++) {
00311 if ((((rts-xh)*df-f)*((rts-xl)*df-f) > 0.0)
00312 || (fabs(2.0*f) > fabs(dxold*df))) {
00313 dxold=dx;
00314 dx=0.5*(xh-xl);
00315 rts=xl+dx;
00316 if (xl == rts) return rts;
00317 } else {
00318 dxold=dx;
00319 dx=f/df;
00320 temp=rts;
00321 rts -= dx;
00322 if (temp == rts) return rts;
00323 }
00324 if (fabs(dx) < xacc) return rts;
00325 (*funcd)(rts,vl,vr,dmin,dmax,&f,&df);
00326 if (f < 0.0)
00327 xl=rts;
00328 else
00329 xh=rts;
00330 }
00331
00332 return rts;
00333 }
00334
00335 #undef MAXIT
00336
00337 #elif defined ADIABATIC
00338
00339
00340
00341
00342
00343
00344
00345
00346
00347 static Real guessP(const Prim1DS Wl, const Prim1DS Wr)
00348 {
00349
00350 Real al = sqrt((double) Gamma*Wl.P/Wl.d);
00351 Real ar = sqrt((double) Gamma*Wr.P/Wr.d);
00352 Real ppv;
00353 Real gl, gr, p_0;
00354 Real TOL = 1.0e-6;
00355
00356 ppv = 0.5*(Wl.P + Wr.P)-0.125*(Wr.Vx-Wl.Vx)*(Wl.d+Wr.d)*(al+ar);
00357
00358 if (ppv < 0.0)
00359 ppv = 0.0;
00360
00361 gl = sqrt((double) (2.0/(Wl.d*(Gamma+1)))/
00362 ((Gamma-1)*Wl.P/(Gamma+1) + ppv));
00363 gr = sqrt((double) (2.0/(Wr.d*(Gamma+1)))/
00364 ((Gamma-1)*Wr.P/(Gamma+1) + ppv));
00365 p_0 = (gl*Wl.P + gr*Wr.P - (Wr.Vx-Wl.Vx))/(gr + gl);
00366
00367 if (p_0 < 0.0)
00368 p_0 = TOL;
00369
00370 return p_0;
00371 }
00372
00373
00374
00375
00376
00377 static Real PFunc(const Prim1DS W, Real POld)
00378 {
00379 Real f;
00380 Real a = sqrt((double) Gamma*W.P/W.d);
00381 Real Ak, Bk;
00382
00383
00384 if (POld <= W.P)
00385 f = 2*a/(Gamma-1) * (pow((double) POld/W.P, (double)
00386 (Gamma-1)/(2*Gamma)) - 1);
00387
00388
00389 else {
00390 Ak = 2/(W.d*(Gamma+1));
00391 Bk = W.P*(Gamma-1)/(Gamma+1);
00392
00393 f = (POld - W.P) * sqrt((double) Ak/(POld + Bk));
00394 }
00395
00396 return f;
00397 }
00398
00399
00400
00401
00402
00403
00404 static Real PFuncDeriv(const Prim1DS W, Real POld)
00405 {
00406 Real fDer;
00407 Real a = sqrt((double) Gamma*W.P/W.d);
00408 Real Ak, Bk;
00409
00410
00411 if (POld <= W.P)
00412 fDer = 1/(a * W.d) * pow((double) POld/W.P, (double) -
00413 (Gamma+1)/(2*Gamma));
00414
00415
00416 else {
00417 Ak = 2/(W.d*(Gamma+1));
00418 Bk = W.P*(Gamma-1)/(Gamma+1);
00419
00420 fDer = sqrt((double) Ak/(POld + Bk))
00421 *(1.0 - 0.5*(POld - W.P)/(Bk + POld));
00422 }
00423
00424 return fDer;
00425 }
00426
00427
00428
00429
00430
00431
00432
00433 static Real getPC(const Prim1DS Wl, const Prim1DS Wr)
00434 {
00435 Real POld;
00436 Real VxDiff = Wr.Vx - Wl.Vx;
00437 int i = 0;
00438 Real TOL = 1.0e-6;
00439 Real change, p, fr, fl, frder, flder;
00440 Real MAX_ITER = 100;
00441
00442 POld = guessP(Wl, Wr);
00443
00444 while (i < MAX_ITER) {
00445
00446 fl = PFunc(Wl, POld);
00447 fr = PFunc(Wr, POld);
00448 flder = PFuncDeriv(Wl, POld);
00449 frder = PFuncDeriv(Wr, POld);
00450
00451 p = POld - (fl + fr + VxDiff)/(flder + frder);
00452
00453 change = 2.0*fabs((p-POld)/(p+POld));
00454 if (change <= TOL)
00455 return p;
00456
00457 if (p < 0.0)
00458 p = TOL;
00459
00460 POld = p;
00461 i++;
00462 }
00463
00464
00465 ath_error("[exact_flux]: Divergence in Newton-Raphson \
00466 iteration: p = %e\n", p);
00467
00468 }
00469
00470
00471
00472
00473
00474
00475
00476
00477
00478
00479
00480
00481
00482 void fluxes(const Cons1DS Ul, const Cons1DS Ur,
00483 const Prim1DS Wl, const Prim1DS Wr,
00484 const Real Bxi, Cons1DS *pF)
00485 {
00486 Real pc, fr, fl;
00487 Real dc, dcl, dcr, Vxc;
00488 Real sl, sr;
00489 Real hdl, hdr;
00490 Real tll, tlr;
00491 Real al = sqrt((double) Gamma*Wl.P/Wl.d);
00492 Real ar = sqrt((double) Gamma*Wr.P/Wr.d);
00493 Real tmp1, tmp2;
00494 Real e, V, E;
00495
00496 if(!(Ul.d > 0.0)||!(Ur.d > 0.0))
00497 ath_error("[exact flux]: Non-positive densities: dl = %e dr = %e\n",
00498 Ul.d, Ur.d);
00499
00500 pc = getPC(Wl, Wr);
00501
00502
00503
00504 fr = PFunc(Wr, pc);
00505 fl = PFunc(Wl, pc);
00506 Vxc = 0.5*(Wl.Vx + Wr.Vx) + 0.5*(fr - fl);
00507
00508
00509
00510 if (pc > Wl.P) {
00511
00512
00513 Real tmp = (Gamma - 1.0) / (Gamma + 1.0);
00514 dcl = Wl.d*(pc/Wl.P + tmp)/(tmp*pc/Wl.P + 1);
00515 }
00516 else {
00517
00518
00519 dcl = Wl.d*pow((double) (pc/Wl.P), (double) (1/Gamma));
00520 }
00521
00522 if (dcl < 0.0)
00523 ath_error("[exact flux]: Solver finds negative density %5.4e\n", dcl);
00524
00525
00526
00527 if (pc > Wr.P) {
00528
00529
00530 Real tmp = (Gamma - 1)/(Gamma + 1);
00531 dcr = Wr.d*(pc/Wr.P + tmp)/(tmp*pc/Wr.P + 1);
00532 }
00533 else {
00534
00535
00536 dcr = Wr.d*pow((double) (pc/Wr.P), (double) (1/Gamma));
00537 }
00538
00539 if (dcr < 0.0)
00540 ath_error("[exact flux]: Solver finds negative density %5.4e\n", dcr);
00541
00542
00543
00544
00545 if (pc > Wl.P) {
00546
00547
00548
00549 sl = Wl.Vx - al*sqrt((double)(pc*(Gamma+1)/(2*Gamma*Wl.P) +
00550 (Gamma-1)/(2*Gamma)));
00551 if (sl >= 0.0) {
00552
00553
00554 e = Wl.P/(Wl.d*(Gamma-1));
00555 V = Wl.Vx*Wl.Vx + Wl.Vy*Wl.Vy + Wl.Vz*Wl.Vz;
00556 E = Wl.d*(0.5*V + e);
00557
00558 pF->E = Wl.Vx*(E + Wl.P);
00559 pF->d = Ul.Mx;
00560 pF->Mx = Ul.Mx*(Wl.Vx) + Wl.P;
00561 pF->My = Ul.My*(Wl.Vx);
00562 pF->Mz = Ul.Mz*(Wl.Vx);
00563
00564 return;
00565 }
00566 }
00567 else {
00568
00569
00570 Real alc = al*pow((double)(pc/Wl.P), (double)(Gamma-1)/(2*Gamma));
00571
00572 hdl = Wl.Vx - al;
00573 tll = Vxc - alc;
00574
00575 if (hdl >= 0.0) {
00576
00577
00578 e = Wl.P/(Wl.d*(Gamma-1));
00579 V = Wl.Vx*Wl.Vx + Wl.Vy*Wl.Vy + Wl.Vz*Wl.Vz;
00580 E = Wl.d*(0.5*V + e);
00581
00582 pF->E = Wl.Vx*(E + Wl.P);
00583 pF->d = Ul.Mx;
00584 pF->Mx = Ul.Mx*(Wl.Vx) + Wl.P;
00585 pF->My = Ul.My*(Wl.Vx);
00586 pF->Mz = Ul.Mz*(Wl.Vx);
00587 return;
00588
00589 }
00590 else if (tll >= 0.0) {
00591
00592
00593
00594 tmp1 = 2/(Gamma + 1);
00595 tmp2 = (Gamma - 1)/(al*(Gamma+1));
00596
00597 dc = Wl.d*pow((double)(tmp1 + tmp2*Wl.Vx), (double)(2/(Gamma-1)));
00598 Vxc = tmp1*(al + Wl.Vx*(Gamma-1)/2);
00599 pc = Wl.P*pow((double)(tmp1+tmp2*Wl.Vx),(double)(2*Gamma/(Gamma-1)));
00600
00601 e = pc/(dc*(Gamma-1));
00602 V = Vxc*Vxc + Wl.Vy*Wl.Vy + Wl.Vz*Wl.Vz;
00603 E = dc*(0.5*V + e);
00604
00605 pF->E = Vxc*(E + pc);
00606 pF->d = dc*Vxc;
00607 pF->Mx = dc*Vxc*Vxc + pc;
00608 pF->My = dc*Vxc*Wl.Vy;
00609 pF->Mz = dc*Vxc*Wl.Vz;
00610 return;
00611 }
00612 }
00613
00614 if (pc > Wr.P) {
00615
00616
00617
00618 sr = Wr.Vx + ar*sqrt((double)(pc*(Gamma+1)/(2*Gamma*Wr.P) +
00619 (Gamma-1)/(2*Gamma)));
00620 if (sr <= 0.0) {
00621
00622
00623 e = Wr.P/(Wr.d*(Gamma-1));
00624 V = Wr.Vx*Wr.Vx + Wr.Vy*Wr.Vy + Wr.Vz*Wr.Vz;
00625 E = Wr.d*(0.5*V + e);
00626
00627 pF->E = Wr.Vx*(E + Wr.P);
00628 pF->d = Ur.Mx;
00629 pF->Mx = Ur.Mx*(Wr.Vx) + Wr.P;
00630 pF->My = Ur.My*(Wr.Vx);
00631 pF->Mz = Ur.Mz*(Wr.Vx);
00632 return;
00633 }
00634 }
00635 else {
00636
00637
00638 Real arc = ar*pow((double)(pc/Wr.P), (double)(Gamma-1)/(2*Gamma));
00639
00640 hdr = Wr.Vx + ar;
00641 tlr = Vxc + arc;
00642
00643 if (hdr <= 0.0) {
00644
00645
00646 e = Wr.P/(Wr.d*(Gamma-1));
00647 V = Wr.Vx*Wr.Vx + Wr.Vy*Wr.Vy + Wr.Vz*Wr.Vz;
00648 E = Wr.d*(0.5*V + e);
00649
00650 pF->E = Wr.Vx*(E + Wr.P);
00651 pF->d = Ur.Mx;
00652 pF->Mx = Ur.Mx*(Wr.Vx) + Wr.P;
00653 pF->My = Ur.My*(Wr.Vx);
00654 pF->Mz = Ur.Mz*(Wr.Vx);
00655 return;
00656
00657 } else if (tlr <= 0.0) {
00658
00659
00660 tmp1 = 2/(Gamma + 1);
00661 tmp2 = (Gamma - 1)/(ar*(Gamma+1));
00662
00663 dc = Wr.d*pow((double)(tmp1 - tmp2*Wr.Vx), (double)(2/(Gamma-1)));
00664 Vxc = tmp1*(-ar + Wr.Vx*(Gamma-1)/2);
00665 pc = Wr.P*pow((double)(tmp1-tmp2*Wr.Vx), (double)(2*Gamma/(Gamma-1)));
00666
00667 e = pc/(dc*(Gamma-1));
00668 V = Vxc*Vxc + Wr.Vy*Wr.Vy + Wr.Vz*Wr.Vz;
00669 E = dc*(0.5*V + e);
00670
00671 pF->E = Vxc*(E + pc);
00672 pF->d = dc*Vxc;
00673 pF->Mx = dc*Vxc*Vxc + pc;
00674 pF->My = dc*Vxc*Wr.Vy;
00675 pF->Mz = dc*Vxc*Wr.Vz;
00676 return;
00677 }
00678 }
00679
00680
00681
00682
00683
00684 if (Vxc >= 0.0) {
00685
00686 e = pc/(dcl*(Gamma-1));
00687 V = Vxc*Vxc + Wl.Vy*Wl.Vy + Wl.Vz*Wl.Vz;
00688 E = dcl*(0.5*V + e);
00689
00690 pF->E = Vxc*(E + pc);
00691 pF->d = dcl*Vxc;
00692 pF->Mx = dcl*Vxc*Vxc + pc;
00693 pF->My = dcl*Vxc*Wl.Vy;
00694 pF->Mz = dcl*Vxc*Wl.Vz;
00695 }
00696 else {
00697
00698 e = pc/(dcr*(Gamma-1));
00699 V = Vxc*Vxc + Wr.Vy*Wr.Vy + Wr.Vz*Wr.Vz;
00700 E = dcr*(0.5*V + e);
00701
00702 pF->E = Vxc*(E + pc);
00703 pF->d = dcr*Vxc;
00704 pF->Mx = dcr*Vxc*Vxc + pc;
00705 pF->My = dcr*Vxc*Wr.Vy;
00706 pF->Mz = dcr*Vxc*Wr.Vz;
00707 }
00708
00709 return;
00710 }
00711
00712
00713 #endif
00714
00715 #endif
00716 #endif