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 #include <math.h>
00036 #include "defs.h"
00037 #include "athena.h"
00038 #include "globals.h"
00039 #include "prototypes.h"
00040
00041 #if defined(SPECIAL_RELATIVITY) && defined(MHD)
00042
00043 Real Gamma_1overGamma;
00044
00045
00046
00047
00048 Prim1DS vsq1D_fix (const Cons1DS *pU, const Real *pBx);
00049
00050
00051
00052
00053 Prim1DS entropy_fix1D (const Cons1DS *pU, const Real *pBx, const Real *ent);
00054
00055
00056
00057
00058 static Real calc_func (Real Q, Real E, Real Bsq, Real Ssq, Real Vsq, Real pgas);
00059
00060
00061
00062
00063 static Real calc_dfunc (Real Q, Real Bsq, Real Msq, Real Ssq, Real d, Real Vsq,
00064 Real Gsq, Real Chi);
00065
00066
00067
00068 static Real calc_ent_func (Real rho, Real pgas, Real d, Real ent);
00069
00070
00071
00072
00073
00074 static Real calc_ent_dfunc(Real Q, Real Bsq, Real Msq, Real Ssq,
00075 Real d, Real Vsq, Real Gsq, Real Chi,
00076 Real pgas, Real rho);
00077
00078
00079
00080 static Real calc_vsq (Real Bsq, Real Msq, Real Ssq, Real Q);
00081
00082
00083
00084 static Real calc_chi (Real d, Real Vsq, Real Gsq, Real Q);
00085
00086 void FUNV2(const Real d, const Real v2,
00087 const Real p, const Real Bsq, const Real Msq, const Real Ssq,
00088 Real *W, Real *f);
00089
00090
00091 static Real tol=1.0e-10;
00092 #endif
00093
00094
00095
00096
00097
00098
00099
00100
00101
00102
00103
00104 PrimS Cons_to_Prim(const ConsS *pCons)
00105 {
00106 Cons1DS U;
00107 Prim1DS W;
00108 PrimS Prim;
00109 #if (NSCALARS > 0)
00110 int n;
00111 #endif
00112 Real Bx=0.0;
00113
00114 U.d = pCons->d;
00115 U.Mx = pCons->M1;
00116 U.My = pCons->M2;
00117 U.Mz = pCons->M3;
00118 #ifndef ISOTHERMAL
00119 U.E = pCons->E;
00120 #endif
00121 #ifdef MHD
00122 Bx = pCons->B1c;
00123 U.By = pCons->B2c;
00124 U.Bz = pCons->B3c;
00125 #endif
00126 #if (NSCALARS > 0)
00127 for (n=0; n<NSCALARS; n++) U.s[n] = pCons->s[n];
00128 #endif
00129
00130 W = Cons1D_to_Prim1D(&U, &Bx);
00131
00132 Prim.d = W.d;
00133 Prim.V1 = W.Vx;
00134 Prim.V2 = W.Vy;
00135 Prim.V3 = W.Vz;
00136 #ifndef ISOTHERMAL
00137 Prim.P = W.P;
00138 #endif
00139 #ifdef MHD
00140 Prim.B1c = Bx;
00141 Prim.B2c = W.By;
00142 Prim.B3c = W.Bz;
00143 #endif
00144 #if (NSCALARS > 0)
00145 for (n=0; n<NSCALARS; n++) Prim.r[n] = W.r[n];
00146 #endif
00147
00148 return Prim;
00149 }
00150
00151
00152
00153
00154
00155
00156
00157
00158
00159 ConsS Prim_to_Cons (const PrimS *pW)
00160 {
00161 Cons1DS U;
00162 Prim1DS W;
00163 ConsS Cons;
00164 #if (NSCALARS > 0)
00165 int n;
00166 #endif
00167 Real Bx=0.0;
00168
00169 W.d = pW->d;
00170 W.Vx = pW->V1;
00171 W.Vy = pW->V2;
00172 W.Vz = pW->V3;
00173 #ifndef ISOTHERMAL
00174 W.P = pW->P;
00175 #endif
00176 #ifdef MHD
00177 Bx = pW->B1c;
00178 W.By = pW->B2c;
00179 W.Bz = pW->B3c;
00180 #endif
00181 #if (NSCALARS > 0)
00182 for (n=0; n<NSCALARS; n++) W.r[n] = pW->r[n];
00183 #endif
00184
00185 U = Prim1D_to_Cons1D(&W, &Bx);
00186
00187 Cons.d = U.d;
00188 Cons.M1 = U.Mx;
00189 Cons.M2 = U.My;
00190 Cons.M3 = U.Mz;
00191 #ifndef ISOTHERMAL
00192 Cons.E = U.E;
00193 #endif
00194 #ifdef MHD
00195 Cons.B1c = Bx;
00196 Cons.B2c = U.By;
00197 Cons.B3c = U.Bz;
00198 #endif
00199 #if (NSCALARS > 0)
00200 for (n=0; n<NSCALARS; n++) Cons.s[n] = U.s[n];
00201 #endif
00202
00203 return Cons;
00204 }
00205
00206 #ifdef SPECIAL_RELATIVITY
00207 #ifdef MHD
00208
00209
00210
00211
00212
00213
00214
00215
00216 PrimS fix_vsq (const ConsS *pCons)
00217 {
00218 Cons1DS U;
00219 Prim1DS W;
00220 PrimS Prim;
00221 #if (NSCALARS > 0)
00222 int n;
00223 #endif
00224 Real Bx=0.0;
00225
00226 U.d = pCons->d;
00227 U.Mx = pCons->M1;
00228 U.My = pCons->M2;
00229 U.Mz = pCons->M3;
00230 #ifndef ISOTHERMAL
00231 U.E = pCons->E;
00232 #endif
00233 #ifdef MHD
00234 Bx = pCons->B1c;
00235 U.By = pCons->B2c;
00236 U.Bz = pCons->B3c;
00237 #endif
00238 #if (NSCALARS > 0)
00239 for (n=0; n<NSCALARS; n++) U.s[n] = pCons->s[n];
00240 #endif
00241
00242 W = vsq1D_fix (&U, &Bx);
00243
00244 Prim.d = W.d;
00245 Prim.V1 = W.Vx;
00246 Prim.V2 = W.Vy;
00247 Prim.V3 = W.Vz;
00248 #ifndef ISOTHERMAL
00249 Prim.P = W.P;
00250 #endif
00251 #ifdef MHD
00252 Prim.B1c = Bx;
00253 Prim.B2c = W.By;
00254 Prim.B3c = W.Bz;
00255 #endif
00256 #if (NSCALARS > 0)
00257 for (n=0; n<NSCALARS; n++) Prim.r[n] = W.r[n];
00258 #endif
00259
00260 return Prim;
00261 }
00262
00263
00264
00265
00266
00267
00268
00269
00270
00271 PrimS entropy_fix (const ConsS *pCons, const Real *ent)
00272 {
00273 Cons1DS U;
00274 Prim1DS W;
00275 PrimS Prim;
00276 #if (NSCALARS > 0)
00277 int n;
00278 #endif
00279 Real Bx=0.0;
00280
00281 U.d = pCons->d;
00282 U.Mx = pCons->M1;
00283 U.My = pCons->M2;
00284 U.Mz = pCons->M3;
00285 #ifndef ISOTHERMAL
00286 U.E = pCons->E;
00287 #endif
00288 #ifdef MHD
00289 Bx = pCons->B1c;
00290 U.By = pCons->B2c;
00291 U.Bz = pCons->B3c;
00292 #endif
00293 #if (NSCALARS > 0)
00294 for (n=0; n<NSCALARS; n++) U.s[n] = pCons->s[n];
00295 #endif
00296
00297 W = entropy_fix1D (&U, &Bx, ent);
00298
00299 Prim.d = W.d;
00300 Prim.V1 = W.Vx;
00301 Prim.V2 = W.Vy;
00302 Prim.V3 = W.Vz;
00303 #ifndef ISOTHERMAL
00304 Prim.P = W.P;
00305 #endif
00306 #ifdef MHD
00307 Prim.B1c = Bx;
00308 Prim.B2c = W.By;
00309 Prim.B3c = W.Bz;
00310 #endif
00311 #if (NSCALARS > 0)
00312 for (n=0; n<NSCALARS; n++) Prim.r[n] = W.r[n];
00313 #endif
00314
00315 return Prim;
00316 }
00317 #endif
00318
00319
00320
00321
00322
00323
00324
00325
00326
00327 PrimS check_Prim(const ConsS *pCons)
00328 {
00329 Cons1DS U;
00330 Prim1DS W;
00331 PrimS Prim;
00332 #if (NSCALARS > 0)
00333 int n;
00334 #endif
00335 Real Bx=0.0;
00336
00337 U.d = pCons->d;
00338 U.Mx = pCons->M1;
00339 U.My = pCons->M2;
00340 U.Mz = pCons->M3;
00341 #ifndef ISOTHERMAL
00342 U.E = pCons->E;
00343 #endif
00344 #ifdef MHD
00345 Bx = pCons->B1c;
00346 U.By = pCons->B2c;
00347 U.Bz = pCons->B3c;
00348 #endif
00349 #if (NSCALARS > 0)
00350 for (n=0; n<NSCALARS; n++) U.s[n] = pCons->s[n];
00351 #endif
00352
00353 #ifdef MHD
00354 W = check_Prim1D(&U, &Bx);
00355 #else
00356 W = Cons1D_to_Prim1D(&U, &Bx);
00357 #endif
00358
00359 Prim.d = W.d;
00360 Prim.V1 = W.Vx;
00361 Prim.V2 = W.Vy;
00362 Prim.V3 = W.Vz;
00363 #ifndef ISOTHERMAL
00364 Prim.P = W.P;
00365 #endif
00366 #ifdef MHD
00367 Prim.B1c = Bx;
00368 Prim.B2c = W.By;
00369 Prim.B3c = W.Bz;
00370 #endif
00371 #if (NSCALARS > 0)
00372 for (n=0; n<NSCALARS; n++) Prim.r[n] = W.r[n];
00373 #endif
00374
00375 return Prim;
00376 }
00377 #endif
00378
00379 #ifndef SPECIAL_RELATIVITY
00380
00381
00382
00383
00384
00385
00386
00387
00388
00389 Prim1DS Cons1D_to_Prim1D(const Cons1DS *pU, const Real *pBx)
00390 {
00391 Prim1DS Prim1D;
00392 #if (NSCALARS > 0)
00393 int n;
00394 #endif
00395 Real di = 1.0/pU->d;
00396
00397 Prim1D.d = pU->d;
00398 Prim1D.Vx = pU->Mx*di;
00399 Prim1D.Vy = pU->My*di;
00400 Prim1D.Vz = pU->Mz*di;
00401
00402 #ifndef ISOTHERMAL
00403 Prim1D.P = pU->E - 0.5*(SQR(pU->Mx)+SQR(pU->My)+SQR(pU->Mz))*di;
00404 #ifdef MHD
00405 Prim1D.P -= 0.5*(SQR(*pBx) + SQR(pU->By) + SQR(pU->Bz));
00406 #endif
00407 Prim1D.P *= Gamma_1;
00408 Prim1D.P = MAX(Prim1D.P,TINY_NUMBER);
00409 #endif
00410
00411 #ifdef MHD
00412 Prim1D.By = pU->By;
00413 Prim1D.Bz = pU->Bz;
00414 #endif
00415
00416 #if (NSCALARS > 0)
00417 for (n=0; n<NSCALARS; n++) Prim1D.r[n] = pU->s[n]*di;
00418 #endif
00419
00420 return Prim1D;
00421 }
00422
00423
00424
00425
00426
00427
00428
00429
00430
00431
00432 Cons1DS Prim1D_to_Cons1D(const Prim1DS *pW, const Real *pBx)
00433 {
00434 Cons1DS Cons1D;
00435 #if (NSCALARS > 0)
00436 int n;
00437 #endif
00438
00439 Cons1D.d = pW->d;
00440 Cons1D.Mx = pW->d*pW->Vx;
00441 Cons1D.My = pW->d*pW->Vy;
00442 Cons1D.Mz = pW->d*pW->Vz;
00443
00444 #ifndef ISOTHERMAL
00445 Cons1D.E = pW->P/Gamma_1 + 0.5*pW->d*(SQR(pW->Vx) +SQR(pW->Vy) +SQR(pW->Vz));
00446 #ifdef MHD
00447 Cons1D.E += 0.5*(SQR(*pBx) + SQR(pW->By) + SQR(pW->Bz));
00448 #endif
00449 #endif
00450
00451 #ifdef MHD
00452 Cons1D.By = pW->By;
00453 Cons1D.Bz = pW->Bz;
00454 #endif
00455
00456 #if (NSCALARS > 0)
00457 for (n=0; n<NSCALARS; n++) Cons1D.s[n] = pW->r[n]*pW->d;
00458 #endif
00459
00460 return Cons1D;
00461 }
00462
00463
00464
00465
00466
00467
00468
00469
00470 Real cfast(const Cons1DS *U, const Real *Bx)
00471 {
00472 Real asq;
00473 #ifndef ISOTHERMAL
00474 Real p,pb=0.0;
00475 #endif
00476 #ifdef MHD
00477 Real ctsq,casq,tmp,cfsq;
00478 #endif
00479
00480 #ifdef ISOTHERMAL
00481 asq = Iso_csound2;
00482 #else
00483 #ifdef MHD
00484 pb = 0.5*(SQR(*Bx) + SQR(U->By) + SQR(U->Bz));
00485 #endif
00486 p = Gamma_1*(U->E - pb - 0.5*(SQR(U->Mx)+SQR(U->My)+SQR(U->Mz))/U->d);
00487 asq = Gamma*p/U->d;
00488 #endif
00489
00490 #ifndef MHD
00491 return sqrt(asq);
00492 #else
00493 ctsq = (SQR(U->By) + SQR(U->Bz))/U->d;
00494 casq = SQR(*Bx)/U->d;
00495 tmp = casq + ctsq - asq;
00496 cfsq = 0.5*((asq+ctsq+casq) + sqrt(tmp*tmp + 4.0*asq*ctsq));
00497 return sqrt(cfsq);
00498 #endif
00499 }
00500 #endif
00501
00502 #if defined(SPECIAL_RELATIVITY) && defined(HYDRO)
00503
00504
00505
00506
00507
00508
00509
00510
00511
00512 Prim1DS Cons1D_to_Prim1D(const Cons1DS *U, const Real *Bx)
00513 {
00514 Prim1DS Prim1D;
00515 Real Msq, M, ME, Dsq, Gamma_1sq, denom;
00516 Real a3, a2, a1, a0;
00517 Real i1, i2, i3;
00518 Real iR, iS, iT;
00519 Real ix1=0.0, iB, iC, v, vOverM;
00520
00521
00522
00523 Msq = SQR(U->Mx) + SQR(U->My) + SQR(U->Mz);
00524 M = sqrt(Msq);
00525
00526 if (fabs(M) < TINY_NUMBER) {
00527
00528 v = 0.0;
00529
00530 } else {
00531
00532 ME = M * U->E;
00533 Dsq = SQR(U->d);
00534
00535 Gamma_1sq = SQR(Gamma_1);
00536
00537 denom = 1.0 / (Gamma_1sq * (Msq + Dsq));
00538
00539 a3 = (-2.0 * Gamma * Gamma_1 * ME) * denom;
00540 a2 = (SQR(Gamma) * SQR(U->E) + 2.0*Gamma_1*Msq - Gamma_1sq*Dsq)*denom;
00541 a1 = (-2.0 * Gamma * ME) * denom;
00542 a0 = Msq * denom;
00543
00544
00545
00546 i1 = -a2;
00547 i2 = a3 * a1 - 4.0 * a0;
00548 i3 = 4.0 * a2 * a0 - SQR(a1) - SQR(a3) * a0;
00549
00550 iR = (9.0 * i1 * i2 - 27.0 * i3 - 2.0 * SQR(i1) * i1) / 54.0;
00551 iS = (3.0 * i2 - SQR(a2)) / 9.0;
00552 iT = SQR(iR) + SQR(iS) * iS;
00553
00554
00555
00556
00557 if (iT < 0) {
00558 ix1 = 2.0*pow(sqrt(iR*iR + iT),(ONE_3RD))*cos(atan2(sqrt(-iT),iR)/3.0)
00559 - i1/3.0;
00560 } else {
00561 ix1 = pow((iR + sqrt(iT)),(ONE_3RD)) + pow((iR - sqrt(iT)),(ONE_3RD))
00562 - i1/3.0;
00563 }
00564
00565 iB = 0.5*(a3 + sqrt(SQR(a3) - 4.0*a2 + 4.0*ix1));
00566 iC = 0.5*(ix1 - sqrt(SQR(ix1) - 4.0*a0));
00567 v = (-iB + sqrt(SQR(iB) - 4.0*iC))/2.0;
00568
00569 }
00570
00571
00572
00573
00574 v = MAX(v, 0.0);
00575 v = MIN(v, 1.0 - 1.0e-15);
00576
00577 vOverM = v / M;
00578
00579 if (fabs(M) < TINY_NUMBER) {
00580 vOverM = 0.0;
00581 }
00582
00583 Prim1D.d = sqrt(1.0 - SQR(v)) * U->d;
00584
00585 Prim1D.Vx = U->Mx * vOverM;
00586 Prim1D.Vy = U->My * vOverM;
00587 Prim1D.Vz = U->Mz * vOverM;
00588
00589 Prim1D.P = Gamma_1*
00590 ((U->E - U->Mx*Prim1D.Vx - U->My*Prim1D.Vy - U->Mz*Prim1D.Vz) - Prim1D.d);
00591
00592 return Prim1D;
00593 }
00594 #endif
00595
00596 #if defined(SPECIAL_RELATIVITY) && defined(MHD)
00597
00598
00599
00600
00601
00602
00603
00604
00605
00606
00607
00608
00609
00610
00611
00612
00613 Prim1DS Cons1D_to_Prim1D(const Cons1DS *U, const Real *Bx)
00614 {
00615 Prim1DS Prim1D;
00616 Real Bsq = 0.0, Msq = 0.0, S = 0.0, Ssq = 0.0;
00617 Real Qp = 0.0, Q = 0.0, Ep = 0.0, E = 0.0, d = 0.0;
00618 Real scrh1, scrh2, tmp1, tmp2;
00619 Real Usq, Vsq, Gsq, rho, Chi, pgas, ent;
00620 Real fQ, dfQ, dQstep;
00621
00622 int nr_success;
00623 int q_incs;
00624
00625 Gamma_1overGamma = Gamma_1/Gamma;
00626
00627
00628 Bsq = SQR((*Bx)) + SQR(U->By) + SQR(U->Bz);
00629 Msq = SQR(U->Mx) + SQR(U->My) + SQR(U->Mz);
00630 S = U->Mx * (*Bx) + U->My * U->By + U->Mz * U->Bz;
00631 Ssq = SQR(S);
00632
00633
00634 E = U->E;
00635 d = U->d;
00636
00637
00638
00639 scrh1 = -4.0*(E - Bsq);
00640 scrh2 = Msq - 2.0*E*Bsq + Bsq*Bsq;
00641 Q = ( - scrh1 + sqrt(fabs(scrh1*scrh1 - 12.0*scrh2)))/6.0;
00642
00643 nr_success = 0;
00644 if (Q < 0.0) {
00645 Q = d;
00646 } else if (Q != Q) {
00647 nr_success = -4;
00648 }
00649
00650
00651
00652 dQstep = 1.0;
00653 q_incs = 0;
00654 while (nr_success == 0 && q_incs < 100000) {
00655
00656 if (fabs(dQstep) <= tol) nr_success = 1;
00657
00658
00659 Vsq = calc_vsq (Bsq,Msq,Ssq,Q);
00660 if (Vsq != Vsq){
00661 nr_success = -3;
00662 }
00663 Gsq = 1.0/(1.0-Vsq);
00664 Chi = calc_chi (d,Vsq,Gsq,Q);
00665 rho = d / sqrt(fabs(Gsq));
00666 pgas = Gamma_1*Chi/Gamma;
00667
00668
00669 fQ = calc_func (Q, E, Bsq, Ssq, Vsq, pgas);
00670 dfQ = calc_dfunc (Q, Bsq, Msq, Ssq, d, Vsq, Gsq, Chi);
00671
00672
00673
00674 if (fQ != fQ) {
00675 nr_success = -1;
00676 }
00677 if (dfQ != dfQ) {
00678 nr_success = -2;
00679 }
00680
00681
00682
00683 if (fabs(fQ) < 0.1 && q_incs == 0) {
00684 Q *= 10;
00685 Vsq = calc_vsq (Bsq,Msq,Ssq,Q);
00686 Gsq = 1.0/(1.0-Vsq);
00687 Chi = calc_chi (d,Vsq,Gsq,Q);
00688 rho = d / sqrt(fabs(Gsq));
00689 pgas = Gamma_1*Chi/Gamma;
00690
00691 fQ = calc_func (Q, E, Bsq, Ssq, Vsq, pgas);
00692 dfQ = calc_dfunc (Q, Bsq, Msq, Ssq, d, Vsq, Gsq, Chi);
00693 }
00694
00695
00696 dQstep = fQ / dfQ;
00697 if (dQstep != dQstep){
00698 nr_success = -2;
00699 }
00700 Q -= dQstep;
00701 if (Q != Q){
00702 nr_success = -1;
00703 }
00704 q_incs ++;
00705
00706
00707 }
00708
00709
00710 if (nr_success == 1){
00711
00712
00713 Vsq = calc_vsq (Bsq,Msq,Ssq,Q);
00714 Gsq = 1.0/(1.0-Vsq);
00715 Chi = calc_chi (d,Vsq,Gsq,Q);
00716 rho = d / sqrt(fabs(Gsq));
00717 pgas = Gamma_1*Chi/Gamma;
00718
00719 if (pgas < 0.0) nr_success = 2;
00720 if (Vsq > 1.0) nr_success = 3;
00721 if (Vsq < 0.0) nr_success = 4;
00722 }
00723
00724 if (nr_success == 2) {
00725
00726
00727
00728
00729 tmp1 = 1.0 / Q;
00730 tmp2 = 1.0 / (Q + Bsq);
00731 Prim1D.d = MAX(rho,1.0e-4);
00732 Prim1D.P = MAX(pgas,1.0e-5);
00733 Prim1D.Vx = (U->Mx + S*(*Bx)*tmp1)*tmp2;
00734 Prim1D.Vy = (U->My + S*U->By*tmp1)*tmp2;
00735 Prim1D.Vz = (U->Mz + S*U->Bz*tmp1)*tmp2;
00736
00737 Prim1D.By = U->By;
00738 Prim1D.Bz = U->Bz;
00739 } else if (nr_success == 3) {
00740
00741
00742 tmp1 = 1.0 / Q;
00743 tmp2 = 1.0 / (Q + Bsq);
00744 Prim1D.Vx = (U->Mx + S*(*Bx)*tmp1)*tmp2;
00745 Prim1D.Vy = (U->My + S*U->By*tmp1)*tmp2;
00746 Prim1D.Vz = (U->Mz + S*U->Bz*tmp1)*tmp2;
00747
00748 scrh1 = SQR(Prim1D.Vx) + SQR(Prim1D.Vy) + SQR(Prim1D.Vz);
00749
00750 Prim1D.Vx *= 0.9999/scrh1;
00751 Prim1D.Vy *= 0.9999/scrh1;
00752 Prim1D.Vz *= 0.9999/scrh1;
00753 Vsq = SQR(Prim1D.Vx) + SQR(Prim1D.Vy) + SQR(Prim1D.Vz);
00754
00755 Gsq = 1.0/(1.0-Vsq);
00756 Chi = calc_chi (d,Vsq,Gsq,Q);
00757 rho = d / sqrt(fabs(Gsq));
00758 pgas = Gamma_1*Chi/Gamma;
00759
00760 Prim1D.d = MAX(rho,1.0e-4);
00761 Prim1D.P = MAX(pgas,1.0e-5);
00762 Prim1D.By = U->By;
00763 Prim1D.Bz = U->Bz;
00764
00765 } else if (nr_success == 4) {
00766
00767 Prim1D.d = -1.0;
00768 Prim1D.P = 1.0;
00769 Prim1D.Vx = 1.0;
00770 Prim1D.Vy =1.0;
00771 Prim1D.Vz = 1.0;
00772
00773 Prim1D.By = U->By;
00774 Prim1D.Bz = U->Bz;
00775 } else if (nr_success == 1) {
00776
00777 tmp1 = 1.0 / Q;
00778 tmp2 = 1.0 / (Q + Bsq);
00779 Prim1D.d = MAX(rho,1.0e-4);
00780 Prim1D.P = MAX(pgas,1.0e-5);
00781 Prim1D.Vx = (U->Mx + S*(*Bx)*tmp1)*tmp2;
00782 Prim1D.Vy = (U->My + S*U->By*tmp1)*tmp2;
00783 Prim1D.Vz = (U->Mz + S*U->Bz*tmp1)*tmp2;
00784
00785 Prim1D.By = U->By;
00786 Prim1D.Bz = U->Bz;
00787 } else {
00788
00789 Prim1D.d = -2.0;
00790 Prim1D.P = 2.0;
00791 Prim1D.Vx = 2.0;
00792 Prim1D.Vy = 2.0;
00793 Prim1D.Vz = 2.0;
00794
00795 Prim1D.By = U->By;
00796 Prim1D.Bz = U->Bz;
00797 }
00798
00799 return Prim1D;
00800 }
00801
00802
00803
00804
00805
00806
00807
00808
00809
00810
00811
00812
00813
00814
00815
00816
00817
00818 Prim1DS check_Prim1D (const Cons1DS *U, const Real *Bx)
00819 {
00820 Prim1DS Prim1D;
00821 Real Bsq = 0.0, Msq = 0.0, S = 0.0, Ssq = 0.0;
00822 Real Qp = 0.0, Q = 0.0, Ep = 0.0, E = 0.0, d = 0.0;
00823 Real scrh1, scrh2, tmp1, tmp2;
00824 Real Usq, Vsq, Gsq, rho, Chi, pgas, ent;
00825 Real fQ, dfQ, dQstep;
00826
00827 int nr_success;
00828 int q_incs;
00829
00830 Gamma_1overGamma = Gamma_1/Gamma;
00831
00832
00833 Bsq = SQR((*Bx)) + SQR(U->By) + SQR(U->Bz);
00834 Msq = SQR(U->Mx) + SQR(U->My) + SQR(U->Mz);
00835 S = U->Mx * (*Bx) + U->My * U->By + U->Mz * U->Bz;
00836 Ssq = SQR(S);
00837
00838
00839 E = U->E;
00840 d = U->d;
00841
00842
00843
00844 scrh1 = -4.0*(E - Bsq);
00845 scrh2 = Msq - 2.0*E*Bsq + Bsq*Bsq;
00846 Q = ( - scrh1 + sqrt(fabs(scrh1*scrh1 - 12.0*scrh2)))/6.0;
00847
00848 nr_success = 0;
00849 if (Q < 0.0) {
00850 Q = d;
00851 } else if (Q != Q) {
00852 nr_success = 9;
00853 }
00854
00855
00856
00857 dQstep = 1.0;
00858 q_incs = 0;
00859 while (nr_success == 0 && q_incs < 1000) {
00860
00861 if (fabs(dQstep) <= tol) nr_success = 1;
00862
00863
00864 Vsq = calc_vsq (Bsq,Msq,Ssq,Q);
00865 Gsq = 1.0/(1.0-Vsq);
00866 Chi = calc_chi (d,Vsq,Gsq,Q);
00867 rho = d / sqrt(fabs(Gsq));
00868 pgas = Gamma_1*Chi/Gamma;
00869
00870
00871 fQ = calc_func (Q, E, Bsq, Ssq, Vsq, pgas);
00872 dfQ = calc_dfunc (Q, Bsq, Msq, Ssq, d, Vsq, Gsq, Chi);
00873
00874
00875 if (fQ != fQ) {
00876 nr_success = 8;
00877 }
00878 if (dfQ != dfQ) {
00879 nr_success = 7;
00880 }
00881
00882
00883
00884
00885 if (fabs(fQ) < 0.1 && q_incs == 0) {
00886 Q *= 10;
00887 Vsq = calc_vsq (Bsq,Msq,Ssq,Q);
00888 Gsq = 1.0/(1.0-Vsq);
00889 Chi = calc_chi (d,Vsq,Gsq,Q);
00890 rho = d / sqrt(fabs(Gsq));
00891 pgas = Gamma_1*Chi/Gamma;
00892
00893 fQ = calc_func (Q, E, Bsq, Ssq, Vsq, pgas);
00894 dfQ = calc_dfunc (Q, Bsq, Msq, Ssq, d, Vsq, Gsq, Chi);
00895 }
00896
00897
00898 dQstep = fQ / dfQ;
00899 if (dQstep != dQstep){
00900 nr_success = 6;
00901 }
00902 Q -= dQstep;
00903 if (Q != Q){
00904 nr_success = 5;
00905 }
00906 q_incs ++;
00907
00908
00909
00910 }
00911
00912
00913 if (nr_success == 1){
00914
00915 Vsq = calc_vsq (Bsq,Msq,Ssq,Q);
00916 Gsq = 1.0/(1.0-Vsq);
00917 Chi = calc_chi (d,Vsq,Gsq,Q);
00918 rho = d / sqrt(fabs(Gsq));
00919 pgas = Gamma_1*Chi/Gamma;
00920
00921 tmp1 = 1.0 / Q;
00922 tmp2 = 1.0 / (Q + Bsq);
00923 Prim1D.d = rho;
00924 Prim1D.P = pgas;
00925 Prim1D.Vx = (U->Mx + S*(*Bx)*tmp1)*tmp2;
00926 Prim1D.Vy = (U->My + S*U->By*tmp1)*tmp2;
00927 Prim1D.Vz = (U->Mz + S*U->Bz*tmp1)*tmp2;
00928
00929 Prim1D.By = U->By;
00930 Prim1D.Bz = U->Bz;
00931
00932 } else {
00933 Prim1D.d = -1.0;
00934 Prim1D.P = -1.0;
00935 Prim1D.Vx = 1.0;
00936 Prim1D.Vy = 1.0;
00937 Prim1D.Vz = 1.0;
00938
00939 Prim1D.By = U->By;
00940 Prim1D.Bz = U->Bz;
00941 }
00942
00943 return Prim1D;
00944 }
00945 #endif
00946
00947 #ifdef SPECIAL_RELATIVITY
00948
00949
00950
00951
00952
00953
00954
00955
00956
00957 Cons1DS Prim1D_to_Cons1D(const Prim1DS *W, const Real *Bx)
00958 {
00959 Cons1DS Cons1D;
00960 Real vsq, U0, Bsq = 0.0, vDotB = 0.0, wU0sq;
00961
00962 vsq = SQR(W->Vx) + SQR(W->Vy) + SQR(W->Vz);
00963 U0 = 1.0 / (1.0 - vsq);
00964
00965 #ifdef MHD
00966 Bsq = (*Bx)*(*Bx) + SQR(W->By) + SQR(W->Bz);
00967 vDotB = (*Bx)*W->Vx + W->By*W->Vy + W->Bz*W->Vz;
00968 #endif
00969
00970 wU0sq = (W->d + Gamma/Gamma_1 * W->P)*U0;
00971
00972 Cons1D.d = sqrt(U0) * W->d;
00973 Cons1D.Mx = wU0sq * W->Vx;
00974 Cons1D.My = wU0sq * W->Vy;
00975 Cons1D.Mz = wU0sq * W->Vz;
00976
00977 Cons1D.E = wU0sq - W->P;
00978
00979 #ifdef MHD
00980 Cons1D.Mx += Bsq*W->Vx - vDotB*(*Bx);
00981 Cons1D.My += Bsq*W->Vy - vDotB*W->By;
00982 Cons1D.Mz += Bsq*W->Vz - vDotB*W->Bz;
00983
00984 Cons1D.E += (1.0 + vsq) * Bsq/2.0 - SQR(vDotB) / 2.0;
00985
00986 Cons1D.By = W->By;
00987 Cons1D.Bz = W->Bz;
00988 #endif
00989
00990 return Cons1D;
00991 }
00992 #endif
00993
00994 #if defined(SPECIAL_RELATIVITY) && defined(MHD)
00995
00996
00997
00998
00999
01000
01001
01002
01003
01004
01005
01006
01007 Prim1DS entropy_fix1D (const Cons1DS *U, const Real *Bx, const Real *ent)
01008 {
01009 Prim1DS Prim1D;
01010 Real Bsq = 0.0, Msq = 0.0, S = 0.0, Ssq = 0.0, Ent = 0.0;
01011 Real Qp = 0.0, Q = 0.0, Ep = 0.0, E = 0.0, d = 0.0;
01012 Real scrh1, scrh2, tmp1, tmp2;
01013 Real Usq, Vsq, Gsq, rho, Chi, pgas;
01014 Real fQ, dfQ, dQstep;
01015
01016 int nr_success;
01017 int q_incs;
01018
01019 Gamma_1overGamma = Gamma_1/Gamma;
01020
01021
01022 Bsq = SQR((*Bx)) + SQR(U->By) + SQR(U->Bz);
01023 Msq = SQR(U->Mx) + SQR(U->My) + SQR(U->Mz);
01024 S = U->Mx * (*Bx) + U->My * U->By + U->Mz * U->Bz;
01025 Ssq = SQR(S);
01026
01027
01028 E = U->E;
01029 d = U->d;
01030 Ent = *ent;
01031
01032
01033
01034 scrh1 = -4.0*(E - Bsq);
01035 scrh2 = Msq - 2.0*E*Bsq + Bsq*Bsq;
01036 Q = ( - scrh1 + sqrt(fabs(scrh1*scrh1 - 12.0*scrh2)))/6.0;
01037
01038 nr_success = 0;
01039 if (Q < 0.0) {
01040 Q = d;
01041 } else if (Q != Q) {
01042 nr_success = 9;
01043 }
01044
01045
01046
01047 dQstep = 1.0;
01048 q_incs = 0;
01049 while (nr_success == 0 && q_incs < 1000) {
01050
01051 if (fabs(dQstep) <= tol) nr_success = 1;
01052
01053
01054 Vsq = calc_vsq (Bsq,Msq,Ssq,Q);
01055 Gsq = 1.0/(1.0-Vsq);
01056 Chi = calc_chi (d,Vsq,Gsq,Q);
01057 rho = d / sqrt(fabs(Gsq));
01058 pgas = Gamma_1*Chi/Gamma;
01059
01060
01061 fQ = calc_ent_func (rho, pgas, d, Ent);
01062 dfQ = calc_ent_dfunc (Q, Bsq, Msq, Ssq, d, Vsq, Gsq, Chi, pgas, rho);
01063
01064
01065
01066 if (fQ != fQ) {
01067 nr_success = 8;
01068 }
01069 if (dfQ != dfQ) {
01070 nr_success = 7;
01071 }
01072
01073
01074 dQstep = fQ / dfQ;
01075 if (dQstep != dQstep){
01076 nr_success = 6; }
01077 Q -= dQstep;
01078 if (Q != Q){
01079 nr_success = 5;
01080 }
01081 q_incs ++;
01082
01083
01084
01085 }
01086
01087
01088 if (nr_success == 1){
01089
01090 Vsq = calc_vsq (Bsq,Msq,Ssq,Q);
01091 Gsq = 1.0/(1.0-Vsq);
01092 Chi = calc_chi (d,Vsq,Gsq,Q);
01093 rho = d / sqrt(fabs(Gsq));
01094 pgas = Gamma_1*Chi/Gamma;
01095
01096 tmp1 = 1.0 / Q;
01097 tmp2 = 1.0 / (Q + Bsq);
01098 Prim1D.d = rho;
01099 Prim1D.P = pgas;
01100 Prim1D.Vx = (U->Mx + S*(*Bx)*tmp1)*tmp2;
01101 Prim1D.Vy = (U->My + S*U->By*tmp1)*tmp2;
01102 Prim1D.Vz = (U->Mz + S*U->Bz*tmp1)*tmp2;
01103
01104 Prim1D.By = U->By;
01105 Prim1D.Bz = U->Bz;
01106
01107 } else {
01108 Prim1D.d = -1.0;
01109 Prim1D.P = -1.0;
01110 Prim1D.Vx = 1.0;
01111 Prim1D.Vy = 1.0;
01112 Prim1D.Vz = 1.0;
01113
01114 Prim1D.By = U->By;
01115 Prim1D.Bz = U->Bz;
01116 }
01117
01118 return Prim1D;
01119 }
01120
01121
01122
01123
01124
01125
01126
01127
01128
01129
01130 Prim1DS vsq1D_fix (const Cons1DS *U, const Real *Bx)
01131 {
01132 Prim1DS Prim1D;
01133 Cons1DS Con1D;
01134 int k, done=0;
01135 Real Bsq = 0.0, Msq = 0.0, S = 0.0, Ssq = 0.0;
01136 Real v2, v2c, fc, f, W, dW, lor2, lor;
01137 Real fmin, fmax, v2min, v2max, p, d;
01138
01139 d = 1.0e0;
01140 p = 1.0e-1;
01141
01142
01143 Bsq = SQR((*Bx)) + SQR(U->By) + SQR(U->Bz);
01144 Msq = SQR(U->Mx) + SQR(U->My) + SQR(U->Mz);
01145 S = U->Mx * (*Bx) + U->My * U->By + U->Mz * U->Bz;
01146 Ssq = SQR(S);
01147
01148 v2max = 1.0-1.e-8;
01149 v2c = 0.95;
01150 FUNV2(d, v2c, p, Bsq, Msq, Ssq, &W, &fc);
01151 v2 = 0.96;
01152 for (k = 1; k < 100; k++){
01153 FUNV2(d, v2, p, Bsq, Msq, Ssq, &W, &f);
01154 if (done == 1) break;
01155 dW = (v2 - v2c)/(f - fc)*f;
01156 v2c = v2; fc = f;
01157 v2 -= dW;
01158 v2 = MIN(v2max,v2);
01159 v2 = MAX(v2, 0.0);
01160 if (fabs(v2) < 1.e-9) done = 1;
01161 if (fabs(f) < 1.e-9) done = 1;
01162 }
01163 if (v2 > 1.0 || k >= 100) {
01164 ath_error("[fix_vsq1D]: too many iter while fixing p , v^2 = %f\n", v2);
01165 }
01166
01167 lor2 = 1.0/(1.0 - v2);
01168 lor = sqrt(lor2);
01169
01170 Con1D = *U;
01171 Con1D.d = d;
01172 Con1D.E = W - p + 0.5*(1.0 + v2)*Bsq - 0.5*Ssq/(W*W);
01173 Prim1D = Cons1D_to_Prim1D(&Con1D, Bx);
01174 v2 = SQR(Prim1D.Vx) + SQR(Prim1D.Vy) + SQR(Prim1D.Vz);
01175
01176 return Prim1D;
01177 }
01178
01179
01180
01181
01182
01183 static Real calc_func (Real Q, Real E, Real Bsq, Real Ssq, Real Vsq, Real pgas)
01184 {
01185 Real tmp1;
01186
01187 return Q - pgas + 0.5*(1.0+Vsq)*Bsq - (0.5*Ssq/Q/Q) - E;
01188 }
01189
01190
01191
01192 static Real calc_ent_func (Real rho, Real pgas, Real d, Real ent)
01193 {
01194 Real tmp1;
01195
01196 tmp1 = pgas * pow(rho,-Gamma);
01197 return d*tmp1 - ent;
01198 }
01199
01200
01201
01202
01203 static Real calc_dfunc(Real Q, Real Bsq, Real Msq, Real Ssq, Real d, Real Vsq,
01204 Real Gsq, Real Chi)
01205 {
01206
01207 Real dp_dQ, dp_dchi, dchi_dQ, dp_drho, drho_dQ, dVsq_dQ;
01208 Real G, Qsq, Qth, scrh1;
01209
01210 Qsq = Q*Q;
01211 Qth = Qsq*Q;
01212 G = sqrt(fabs(Gsq));
01213
01214 scrh1 = Q + Bsq;
01215 dVsq_dQ = Ssq*(3.0*Q*scrh1 + Bsq*Bsq) + Msq*Qth;
01216 dVsq_dQ *= -2.0/Qth/(scrh1*scrh1*scrh1);
01217
01218
01219
01220 dchi_dQ = 1.0 - Vsq - 0.5*G*(d + 2.0*Chi*G)*dVsq_dQ;
01221 drho_dQ = -0.5*d*G*dVsq_dQ;
01222
01223
01224
01225 dp_dchi = (Gamma - 1.0)/Gamma;
01226 dp_drho = 0.0;
01227
01228 dp_dQ = dp_dchi*dchi_dQ + dp_drho*drho_dQ;
01229
01230 return 1.0 - dp_dQ +0.5*Bsq*dVsq_dQ + Ssq / Qth;
01231 }
01232
01233
01234
01235
01236
01237 static Real calc_ent_dfunc(Real Q, Real Bsq, Real Msq, Real Ssq,
01238 Real d, Real Vsq, Real Gsq, Real Chi,
01239 Real pgas, Real rho)
01240 {
01241
01242 Real dp_dQ, dp_dchi, dchi_dQ, dp_drho, drho_dQ, dVsq_dQ;
01243 Real G, Qsq, Qth, scrh1;
01244
01245 Qsq = Q*Q;
01246 Qth = Qsq*Q;
01247 G = sqrt(fabs(Gsq));
01248
01249 scrh1 = Q + Bsq;
01250 dVsq_dQ = Ssq*(3.0*Q*scrh1 + Bsq*Bsq) + Msq*Qth;
01251 dVsq_dQ *= -2.0/Qth/(scrh1*scrh1*scrh1);
01252
01253
01254
01255 dchi_dQ = 1.0 - Vsq - 0.5*G*(d + 2.0*Chi*G)*dVsq_dQ;
01256 drho_dQ = -0.5*d*G*dVsq_dQ;
01257
01258
01259
01260 dp_dchi = (Gamma - 1.0)/Gamma;
01261 dp_drho = 0.0;
01262
01263 dp_dQ = dp_dchi*dchi_dQ + dp_drho*drho_dQ;
01264
01265 return d*pow(rho,-Gamma)*dp_dQ - Gamma*pgas*pow(rho,Gamma+1.0)*drho_dQ;
01266 }
01267
01268
01269
01270 static Real calc_vsq (Real Bsq, Real Msq, Real Ssq, Real Q)
01271 {
01272
01273 Real Qsq, Ssq_Qsq, scrh1;
01274 Qsq = Q*Q;
01275 Ssq_Qsq = Ssq / Qsq;
01276 scrh1 = Q + Bsq;
01277 return (Msq + Ssq_Qsq*(scrh1 + Q))/(scrh1*scrh1);
01278 }
01279
01280
01281
01282 static Real calc_chi (Real d, Real Vsq, Real Gsq, Real Q)
01283 {
01284
01285 Real G, tmp1, tmp2;
01286 G = sqrt(fabs(Gsq));
01287 tmp1 = 1.0 - Vsq;
01288 tmp2 = Q - d*G;
01289 return tmp2*tmp1;
01290 }
01291
01292 void FUNV2(const Real d, const Real v2,
01293 const Real p, const Real Bsq, const Real Msq, const Real Ssq,
01294 Real *W, Real *f)
01295 {
01296 Real lor2, lor, W2, pg;
01297
01298 lor2 = 1.0/(1.0 - v2);
01299 lor = sqrt(lor2);
01300 pg = p*lor;
01301 *W = (d + pg*Gamma/(Gamma - 1.0))*lor;
01302 W2 = SQR(*W);
01303
01304 *f = Ssq*(2.0*(*W) + Bsq) + Msq*W2;
01305 *f /= ((*W) + Bsq)*((*W) + Bsq)*W2;
01306 *f -= v2;
01307 }
01308 #endif