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 "../defs.h"
00029 #include "../athena.h"
00030 #include "../globals.h"
00031 #include "prototypes.h"
00032 #include "../prototypes.h"
00033
00034 #ifdef HLLC_FLUX
00035 #ifdef SPECIAL_RELATIVITY
00036
00037 #if (NSCALARS > 0)
00038 #error : The SR HLLC flux does not work with passive scalars.
00039 #endif
00040
00041 #ifdef MHD
00042 void flux_LR(Cons1DS U, Prim1DS W, Cons1DS *flux, Real Bx, Real* p);
00043 void getMaxSignalSpeeds_pluto(const Prim1DS Wl, const Prim1DS Wr,
00044 const Real Bx, Real* low, Real* high);
00045 void getMaxSignalSpeeds_echo(const Prim1DS Wl, const Prim1DS Wr,
00046 const Real Bx, Real* low, Real* high);
00047 void getVChar_pluto(const Prim1DS W, const Real Bx, Real* lml, Real* lmr);
00048 void getVChar_echo (const Prim1DS W, const Real Bx, Real* lml, Real* lmr);
00049
00050
00051
00052
00053 int QUARTIC (Real b, Real c, Real d, Real e, Real z[]);
00054
00055
00056
00057 int CUBIC(Real b, Real c, Real d, Real z[]);
00058 #endif
00059
00060 #ifdef HYDRO
00061
00062
00063
00064
00065
00066
00067
00068
00069
00070
00071
00072 void fluxes(const Cons1DS Ul, const Cons1DS Ur,
00073 const Prim1DS Wl, const Prim1DS Wr, const Real Bxi, Cons1DS *pFlux)
00074 {
00075 Cons1DS Fl,Fr,Fhll,Uhll,Usl,Usr;
00076 Real rhl, rhr, csl, csr, cslsq, csrsq, vsql, vsqr, gammasql, gammasqr;
00077 Real ssl, ssr, radl, radr, lmdapl, lmdapr, lmdaml, lmdamr, lmdatlmda;
00078 Real lmdal,lmdar;
00079 Real lmdas;
00080 Real ovlrmll;
00081 Real a,b,c,quad,rad;
00082 Real den,ps;
00083
00084
00085
00086
00087
00088 rhl = Wl.d + Wl.P * Gamma / Gamma_1;
00089 rhr = Wr.d + Wr.P * Gamma / Gamma_1;
00090
00091 csl = sqrt(Gamma * Wl.P / rhl);
00092 csr = sqrt(Gamma * Wr.P / rhr);
00093
00094 cslsq = SQR(csl);
00095 csrsq = SQR(csr);
00096
00097 vsql = SQR(Wl.Vx) + SQR(Wl.Vy) + SQR(Wl.Vz);
00098 vsqr = SQR(Wr.Vx) + SQR(Wr.Vy) + SQR(Wr.Vz);
00099
00100 gammasql = 1.0 / (1.0 - vsql);
00101 gammasqr = 1.0 / (1.0 - vsqr);
00102
00103 ssl = cslsq / ( gammasql * (1.0 - cslsq) );
00104 ssr = csrsq / ( gammasqr * (1.0 - csrsq) );
00105
00106 radl = sqrt( ssl*(1.0-SQR(Wl.Vx)+ssl) );
00107 radr = sqrt( ssr*(1.0-SQR(Wr.Vx)+ssr) );
00108
00109 lmdapl = (Wl.Vx + radl) / (1.0 + ssl);
00110 lmdapr = (Wr.Vx + radr) / (1.0 + ssr);
00111 lmdaml = (Wl.Vx - radl) / (1.0 + ssl);
00112 lmdamr = (Wr.Vx - radr) / (1.0 + ssr);
00113
00114 lmdal = MIN(lmdaml, lmdamr);
00115 lmdar = MAX(lmdapl, lmdapr);
00116
00117
00118
00119
00120
00121
00122 Fl.d = Ul.d * Wl.Vx;
00123 Fl.Mx = Ul.Mx * Wl.Vx + Wl.P;
00124 Fl.My = Ul.My * Wl.Vx;
00125 Fl.Mz = Ul.Mz * Wl.Vx;
00126 Fl.E = Ul.Mx;
00127
00128 Fr.d = Ur.d * Wr.Vx;
00129 Fr.Mx = Ur.Mx * Wr.Vx + Wr.P;
00130 Fr.My = Ur.My * Wr.Vx;
00131 Fr.Mz = Ur.Mz * Wr.Vx;
00132 Fr.E = Ur.Mx;
00133
00134
00135
00136
00137
00138
00139 ovlrmll = 1.0 / ( lmdar - lmdal );
00140 lmdatlmda = lmdal*lmdar;
00141
00142 Fhll.d = (lmdar*Fl.d - lmdal*Fr.d + lmdatlmda * (Ur.d - Ul.d) ) * ovlrmll;
00143 Fhll.Mx = (lmdar*Fl.Mx - lmdal*Fr.Mx + lmdatlmda * (Ur.Mx - Ul.Mx)) * ovlrmll;
00144 Fhll.My = (lmdar*Fl.My - lmdal*Fr.My + lmdatlmda * (Ur.My - Ul.My)) * ovlrmll;
00145 Fhll.Mz = (lmdar*Fl.Mz - lmdal*Fr.Mz + lmdatlmda * (Ur.Mz - Ul.Mz)) * ovlrmll;
00146 Fhll.E = (lmdar*Fl.E - lmdal*Fr.E + lmdatlmda * (Ur.E - Ul.E )) * ovlrmll;
00147
00148 Uhll.d = (lmdar * Ur.d - lmdal * Ul.d + Fl.d - Fr.d ) * ovlrmll;
00149 Uhll.Mx = (lmdar * Ur.Mx - lmdal * Ul.Mx + Fl.Mx - Fr.Mx) * ovlrmll;
00150 Uhll.My = (lmdar * Ur.My - lmdal * Ul.My + Fl.My - Fr.My) * ovlrmll;
00151 Uhll.Mz = (lmdar * Ur.Mz - lmdal * Ul.Mz + Fl.Mz - Fr.Mz) * ovlrmll;
00152 Uhll.E = (lmdar * Ur.E - lmdal * Ul.E + Fl.E - Fr.E ) * ovlrmll;
00153
00154
00155
00156
00157
00158
00159
00160
00161 a = Fhll.E;
00162 b = -(Uhll.E + Fhll.Mx);
00163 c = Uhll.Mx;
00164
00165
00166 quad = -0.5*(b + SIGN(b)*sqrt(b*b - 4.0*a*c));
00167 lmdas = c/quad;
00168
00169
00170
00171
00172
00173 if( lmdal >= 0.0){
00174
00175 pFlux->d = Fl.d;
00176 pFlux->Mx = Fl.Mx;
00177 pFlux->My = Fl.My;
00178 pFlux->Mz = Fl.Mz;
00179 pFlux->E = Fl.E;
00180
00181 return;
00182 }
00183 else if( lmdas >= 0.0){
00184
00185
00186 ps = -Fhll.E*lmdas + Fhll.Mx;
00187
00188
00189 den = 1.0 / (lmdal - lmdas);
00190
00191 Usl.d = Ul.d * (lmdal - Wl.Vx) * den;
00192 Usl.Mx = (Ul.Mx * (lmdal - Wl.Vx) + ps - Wl.P) * den;
00193 Usl.My = Ul.My * (lmdal - Wl.Vx) * den;
00194 Usl.Mz = Ul.Mz * (lmdal - Wl.Vx) * den;
00195 Usl.E = (Ul.E * (lmdal - Wl.Vx) + ps * lmdas - Wl.P * Wl.Vx) * den;
00196
00197
00198
00199 pFlux->d = lmdal*(Usl.d - Ul.d ) + Fl.d;
00200 pFlux->Mx = lmdal*(Usl.Mx - Ul.Mx) + Fl.Mx;
00201 pFlux->My = lmdal*(Usl.My - Ul.My) + Fl.My;
00202 pFlux->Mz = lmdal*(Usl.Mz - Ul.Mz) + Fl.Mz;
00203 pFlux->E = lmdal*(Usl.E - Ul.E ) + Fl.E;
00204
00205 return;
00206 }
00207 else if( lmdar >= 0.0){
00208
00209
00210 ps = -Fhll.E*lmdas + Fhll.Mx;
00211
00212
00213 den = 1.0 / (lmdar - lmdas);
00214
00215 Usr.d = Ur.d * (lmdar - Wr.Vx) * den;
00216 Usr.Mx = (Ur.Mx * (lmdar - Wr.Vx) + ps - Wr.P) * den;
00217 Usr.My = Ur.My * (lmdar - Wr.Vx) * den;
00218 Usr.Mz = Ur.Mz * (lmdar - Wr.Vx) * den;
00219 Usr.E = (Ur.E * (lmdar - Wr.Vx) + ps * lmdas - Wr.P * Wr.Vx) * den;
00220
00221
00222
00223 pFlux->d = lmdar*(Usr.d - Ur.d ) + Fr.d;
00224 pFlux->Mx = lmdar*(Usr.Mx - Ur.Mx) + Fr.Mx;
00225 pFlux->My = lmdar*(Usr.My - Ur.My) + Fr.My;
00226 pFlux->Mz = lmdar*(Usr.Mz - Ur.Mz) + Fr.Mz;
00227 pFlux->E = lmdar*(Usr.E - Ur.E ) + Fr.E;
00228
00229 return;
00230 }
00231 else{
00232
00233 pFlux->d = Fr.d;
00234 pFlux->Mx = Fr.Mx;
00235 pFlux->My = Fr.My;
00236 pFlux->Mz = Fr.Mz;
00237 pFlux->E = Fr.E;
00238
00239 return;
00240 }
00241
00242
00243 }
00244
00245 #endif
00246
00247 #ifdef MHD
00248
00249
00250
00251
00252
00253
00254
00255
00256
00257
00258
00259 void fluxes(const Cons1DS Ul, const Cons1DS Ur,
00260 const Prim1DS Wl, const Prim1DS Wr, const Real Bxi, Cons1DS *pFlux)
00261 {
00262 Cons1DS Fl,Fr,Fhll,Uhll,Usl,Usr;
00263 Prim1DS Wsl,Wsr,Whll;
00264 Real Sl, Sr, Pl, Pr;
00265 Real Sla, Sra;
00266 Real dS_1, scrh;
00267 Real Bx, Bys, Bzs;
00268 Real BtFBt, Bt2, FBt2;
00269 Real a, b, c;
00270 Real ps, vxs, vys, vzs, gammas_2, vBs, V2l, V2r;
00271 Real vxl, vxr, alpha_l, alpha_r;
00272 int switch_to_hll,wave_speed_fail;
00273
00274 wave_speed_fail = 0;
00275 switch_to_hll = 0;
00276
00277
00278
00279
00280 getMaxSignalSpeeds_pluto(Wl,Wr,Bxi,&Sl,&Sr);
00281
00282 if (Sl != Sl) {
00283 wave_speed_fail = 1;
00284 printf("[hllc_sr_mhd]: NaN in Sl %10.4e %10.4e\n",Sl,Sr);
00285 Sl = -1.0;
00286 Sr = 1.0;
00287 }
00288
00289 if (Sr != Sr) {
00290 wave_speed_fail = 1;
00291 printf("[hllc_sr_mhd]: NaN in Sr %10.4e %10.4e\n",Sl,Sr);
00292 Sl = -1.0;
00293 Sr = 1.0;
00294 }
00295
00296 if (Sl < -1.0) {
00297 wave_speed_fail = 1;
00298 printf("[hllc_sr_mhd]: Superluminal Sl %10.4e %10.4e\n",Sl,Sr);
00299 Sl = -1.0;
00300 Sr = 1.0;
00301 }
00302 if (Sr > 1.0) {
00303 wave_speed_fail = 1;
00304 printf("[hllc_sr_mhd]: Superluminal Sr %10.4e %10.4e\n",Sl,Sr);
00305 Sl = -1.0;
00306 Sr = 1.0;
00307 }
00308
00309
00310
00311
00312 if (wave_speed_fail){
00313 getMaxSignalSpeeds_echo (Wl,Wr,Bxi,&Sla,&Sra);
00314
00315 if (Sla != Sla) {
00316 switch_to_hll = 1;
00317 printf("[hllc_sr_mhd]: NaN in Sl %10.4e %10.4e\n",Sl,Sr);
00318 Sla = -1.0;
00319 Sra = 1.0;
00320 }
00321
00322 if (Sra != Sra) {
00323 switch_to_hll = 1;
00324 printf("[hllc_sr_mhd]: NaN in Sr %10.4e %10.4e\n",Sl,Sr);
00325 Sla = -1.0;
00326 Sra = 1.0;
00327 }
00328
00329 if (Sla < -1.0) {
00330 switch_to_hll = 1;
00331 printf("[hllc_sr_mhd]: Superluminal Sl %10.4e %10.4e\n",Sl,Sr);
00332 Sla = -1.0;
00333 Sra = 1.0;
00334 }
00335 if (Sra > 1.0) {
00336 switch_to_hll = 1;
00337 printf("[hllc_sr_mhd]: Superluminal Sr %10.4e %10.4e\n",Sl,Sr);
00338 Sla = -1.0;
00339 Sra = 1.0;
00340 }
00341
00342 Sl = Sla;
00343 Sr = Sra;
00344
00345 }
00346
00347
00348 flux_LR(Ul,Wl,&Fl,Bxi,&Pl);
00349 flux_LR(Ur,Wr,&Fr,Bxi,&Pr);
00350
00351
00352
00353
00354 dS_1 = 1.0/(Sr - Sl);
00355
00356 Uhll.d = (Sr*Ur.d - Sl*Ul.d + Fl.d - Fr.d ) * dS_1;
00357 Uhll.Mx = (Sr*Ur.Mx - Sl*Ul.Mx + Fl.Mx - Fr.Mx) * dS_1;
00358 Uhll.My = (Sr*Ur.My - Sl*Ul.My + Fl.My - Fr.My) * dS_1;
00359 Uhll.Mz = (Sr*Ur.Mz - Sl*Ul.Mz + Fl.Mz - Fr.Mz) * dS_1;
00360 Uhll.E = (Sr*Ur.E - Sl*Ul.E + Fl.E - Fr.E ) * dS_1;
00361 Uhll.By = (Sr*Ur.By - Sl*Ul.By + Fl.By - Fr.By) * dS_1;
00362 Uhll.Bz = (Sr*Ur.Bz - Sl*Ul.Bz + Fl.Bz - Fr.Bz) * dS_1;
00363
00364 Fhll.d = (Sr*Fl.d - Sl*Fr.d + Sl*Sr*(Ur.d - Ul.d )) * dS_1;
00365 Fhll.Mx = (Sr*Fl.Mx - Sl*Fr.Mx + Sl*Sr*(Ur.Mx - Ul.Mx)) * dS_1;
00366 Fhll.My = (Sr*Fl.My - Sl*Fr.My + Sl*Sr*(Ur.My - Ul.My)) * dS_1;
00367 Fhll.Mz = (Sr*Fl.Mz - Sl*Fr.Mz + Sl*Sr*(Ur.Mz - Ul.Mz)) * dS_1;
00368 Fhll.E = (Sr*Fl.E - Sl*Fr.E + Sl*Sr*(Ur.E - Ul.E )) * dS_1;
00369 Fhll.By = (Sr*Fl.By - Sl*Fr.By + Sl*Sr*(Ur.By - Ul.By)) * dS_1;
00370 Fhll.Bz = (Sr*Fl.Bz - Sl*Fr.Bz + Sl*Sr*(Ur.Bz - Ul.Bz)) * dS_1;
00371
00372
00373 if (switch_to_hll) {
00374 pFlux->d = Fhll.d;
00375 pFlux->Mx = Fhll.Mx;
00376 pFlux->My = Fhll.My;
00377 pFlux->Mz = Fhll.Mz;
00378 pFlux->E = Fhll.E;
00379 pFlux->By = Fhll.By;
00380 pFlux->Bz = Fhll.Bz;
00381
00382 return;
00383 }
00384
00385
00386
00387
00388
00389 if(Sl >= 0.0){
00390 pFlux->d = Fl.d;
00391 pFlux->Mx = Fl.Mx;
00392 pFlux->My = Fl.My;
00393 pFlux->Mz = Fl.Mz;
00394 pFlux->E = Fl.E;
00395 pFlux->By = Fl.By;
00396 pFlux->Bz = Fl.Bz;
00397
00398 return;
00399
00400 }
00401 else if(Sr <= 0.0){
00402 pFlux->d = Fr.d;
00403 pFlux->Mx = Fr.Mx;
00404 pFlux->My = Fr.My;
00405 pFlux->Mz = Fr.Mz;
00406 pFlux->E = Fr.E;
00407 pFlux->By = Fr.By;
00408 pFlux->Bz = Fr.Bz;
00409
00410 return;
00411 }
00412 else {
00413
00414 switch_to_hll = 0;
00415
00416
00417 vxl = Wl.Vx;
00418 vxr = Wr.Vx;
00419
00420 Bx = Bxi;
00421 Bys = Uhll.By;
00422 Bzs = Uhll.Bz;
00423
00424 if (fabs(Bx) < 1.0e-12) {
00425 a = Fhll.E;
00426 b = - (Fhll.Mx + Uhll.E);
00427 c = Uhll.Mx;
00428 } else {
00429 BtFBt = Uhll.By*Fhll.By + Uhll.Bz*Fhll.Bz;
00430 Bt2 = Uhll.By*Uhll.By + Uhll.Bz*Uhll.Bz;
00431 FBt2 = Fhll.By*Fhll.By + Fhll.Bz*Fhll.Bz;
00432
00433 a = Fhll.E - BtFBt;
00434 b = Bt2 + FBt2 - (Fhll.Mx + Uhll.E);
00435 c = Uhll.Mx - BtFBt;
00436 }
00437
00438 if (fabs(a) > 1.e-12){
00439 scrh = 1.0 + sqrt(1.0 - 4.0*a*c/(b*b));
00440 if (scrh != scrh) {
00441 switch_to_hll = 1;
00442 }
00443 vxs = - 2.0*c/(b*scrh);
00444 } else {
00445 vxs = -c/b;
00446 }
00447 if ((vxs != vxs || vxs > 1.0) && (switch_to_hll == 0)) {
00448 switch_to_hll = 1;
00449 }
00450
00451 if (switch_to_hll) {
00452 pFlux->d = Fhll.d;
00453 pFlux->Mx = Fhll.Mx;
00454 pFlux->My = Fhll.My;
00455 pFlux->Mz = Fhll.Mz;
00456 pFlux->E = Fhll.E;
00457 pFlux->By = Fhll.By;
00458 pFlux->Bz = Fhll.Bz;
00459
00460 if (pFlux->d != pFlux->d) {
00461 }
00462
00463 return;
00464 }
00465
00466 if (fabs(Bx) < 1.0e-12) {
00467
00468
00469
00470
00471
00472
00473 ps = Fhll.Mx - Fhll.E*vxs;
00474
00475 if (ps < 0) {
00476 switch_to_hll = 1;
00477 } else {
00478 alpha_l = (Sl - vxl)/(Sl - vxs);
00479 alpha_r = (Sr - vxr)/(Sr - vxs);
00480
00481 Usl.d = Ul.d*alpha_l;
00482 Usr.d = Ur.d*alpha_r;
00483
00484 Usl.E = (Sl*Ul.E - Fl.E + ps*vxs)/(Sl - vxs);
00485 Usr.E = (Sr*Ur.E - Fr.E + ps*vxs)/(Sr - vxs);
00486
00487 Usl.Mx = (Usl.E + ps)*vxs;
00488 Usr.Mx = (Usr.E + ps)*vxs;
00489 Usl.My = Ul.My*alpha_l;
00490 Usr.My = Ur.My*alpha_r;
00491 Usl.Mz = Ul.Mz*alpha_l;
00492 Usr.Mz = Ur.Mz*alpha_r;
00493
00494 Usl.By = Ul.By*alpha_l;
00495 Usr.By = Ur.By*alpha_r;
00496 Usl.Bz = Ul.Bz*alpha_l;
00497 Usr.Bz = Ur.Bz*alpha_r;
00498 }
00499
00500 } else {
00501
00502 vys = (Bys*vxs - Fhll.By)/Bx;
00503 vzs = (Bzs*vxs - Fhll.Bz)/Bx;
00504
00505 gammas_2 = vxs*vxs + vys*vys + vzs*vzs;
00506 gammas_2 = 1.0 - gammas_2;
00507 vBs = vxs*Bx + vys*Bys + vzs*Bzs;
00508
00509 ps = (Bx*vBs - Fhll.E)*vxs + (Bx*Bx*gammas_2) + Fhll.Mx;
00510 if (ps < 0) {
00511 switch_to_hll = 1;
00512 } else {
00513
00514 alpha_l = (Sl - vxl)/(Sl - vxs);
00515 alpha_r = (Sr - vxr)/(Sr - vxs);
00516
00517 if (alpha_l != alpha_l) {
00518 switch_to_hll = 1;
00519 }
00520
00521 if (alpha_r != alpha_r) {
00522 switch_to_hll = 1;
00523
00524 }
00525
00526 if (switch_to_hll == 0){
00527
00528 Usl.d = Ul.d*alpha_l;
00529 Usr.d = Ur.d*alpha_r;
00530
00531 Usl.E = (Sl*Ul.E - Fl.E + ps*vxs - vBs*Bx)/(Sl - vxs);
00532 Usr.E = (Sr*Ur.E - Fr.E + ps*vxs - vBs*Bx)/(Sr - vxs);
00533
00534 Usl.Mx = (Usl.E + ps)*vxs - vBs*Bx;
00535 Usr.Mx = (Usr.E + ps)*vxs - vBs*Bx;
00536 Usl.My = (Sl*Ul.My - Fl.My - Bx*(Bys*gammas_2 + vBs*vys))/(Sl - vxs);
00537 Usr.My = (Sr*Ur.My - Fr.My - Bx*(Bys*gammas_2 + vBs*vys))/(Sr - vxs);
00538 Usl.Mz = (Sl*Ul.Mz - Fl.Mz - Bx*(Bzs*gammas_2 + vBs*vzs))/(Sl - vxs);
00539 Usr.Mz = (Sr*Ur.Mz - Fr.Mz - Bx*(Bzs*gammas_2 + vBs*vzs))/(Sr - vxs);
00540
00541 Usl.By = Usr.By = Bys;
00542 Usl.Bz = Usr.Bz = Bzs;
00543 }
00544 }
00545 }
00546
00547 if (switch_to_hll) {
00548 pFlux->d = Fhll.d;
00549 pFlux->Mx = Fhll.Mx;
00550 pFlux->My = Fhll.My;
00551 pFlux->Mz = Fhll.Mz;
00552 pFlux->E = Fhll.E;
00553 pFlux->By = Fhll.By;
00554 pFlux->Bz = Fhll.Bz;
00555
00556 return;
00557 }
00558
00559
00560 if (vxs > 0.0) {
00561 pFlux->d = Fl.d + Sl*(Usl.d - Ul.d );
00562 pFlux->Mx = Fl.Mx + Sl*(Usl.Mx - Ul.Mx);
00563 pFlux->My = Fl.My + Sl*(Usl.My - Ul.My);
00564 pFlux->Mz = Fl.Mz + Sl*(Usl.Mz - Ul.Mz);
00565 pFlux->E = Fl.E + Sl*(Usl.E - Ul.E );
00566 pFlux->By = Fl.By + Sl*(Usl.By - Ul.By);
00567 pFlux->Bz = Fl.Bz + Sl*(Usl.Bz - Ul.Bz);
00568
00569 if (pFlux->d != pFlux->d) {
00570
00571 pFlux->d = Fhll.d;
00572 pFlux->Mx = Fhll.Mx;
00573 pFlux->My = Fhll.My;
00574 pFlux->Mz = Fhll.Mz;
00575 pFlux->E = Fhll.E;
00576 pFlux->By = Fhll.By;
00577 pFlux->Bz = Fhll.Bz;
00578
00579 }
00580
00581 return;
00582
00583 } else {
00584 pFlux->d = Fr.d + Sr*(Usr.d - Ur.d );
00585 pFlux->Mx = Fr.Mx + Sr*(Usr.Mx - Ur.Mx);
00586 pFlux->My = Fr.My + Sr*(Usr.My - Ur.My);
00587 pFlux->Mz = Fr.Mz + Sr*(Usr.Mz - Ur.Mz);
00588 pFlux->E = Fr.E + Sr*(Usr.E - Ur.E );
00589 pFlux->By = Fr.By + Sr*(Usr.By - Ur.By);
00590 pFlux->Bz = Fr.Bz + Sr*(Usr.Bz - Ur.Bz);
00591
00592 if (pFlux->d != pFlux->d) {
00593
00594 pFlux->d = Fhll.d;
00595 pFlux->Mx = Fhll.Mx;
00596 pFlux->My = Fhll.My;
00597 pFlux->Mz = Fhll.Mz;
00598 pFlux->E = Fhll.E;
00599 pFlux->By = Fhll.By;
00600 pFlux->Bz = Fhll.Bz;
00601
00602 }
00603
00604 return;
00605 }
00606
00607 }
00608
00609 }
00610
00611
00612
00613
00614 void entropy_flux (const Cons1DS Ul, const Cons1DS Ur,
00615 const Prim1DS Wl, const Prim1DS Wr, const Real Bx, Real *pFlux)
00616 {
00617 Real Fl, Fr;
00618 Real USl, USr;
00619 Real WSl, WSr;
00620 Real Uhll, Fhll;
00621 Real Pl, Pr;
00622 Real Sl, Sr;
00623 Real Sla, Sra;
00624 Real dS_1;
00625 int wave_speed_fail;
00626
00627 wave_speed_fail = 0;
00628
00629
00630
00631
00632 getMaxSignalSpeeds_pluto(Wl,Wr,Bx,&Sl,&Sr);
00633
00634 if (Sl != Sl) {
00635 wave_speed_fail = 1;
00636 printf("[hlle_sr_mhd]: NaN in Sl %10.4e %10.4e\n",Sl,Sr);
00637 Sl = -1.0;
00638 Sr = 1.0;
00639 }
00640
00641 if (Sr != Sr) {
00642 wave_speed_fail = 1;
00643 printf("[hlle_sr_mhd]: NaN in Sr %10.4e %10.4e\n",Sl,Sr);
00644 Sl = -1.0;
00645 Sr = 1.0;
00646 }
00647
00648 if (Sl < -1.0) {
00649 wave_speed_fail = 1;
00650 printf("[hlle_sr_mhd]: Superluminal Sl %10.4e %10.4e\n",Sl,Sr);
00651 Sl = -1.0;
00652 Sr = 1.0;
00653 }
00654 if (Sr > 1.0) {
00655 wave_speed_fail = 1;
00656 printf("[hlle_sr_mhd]: Superluminal Sr %10.4e %10.4e\n",Sl,Sr);
00657 Sl = -1.0;
00658 Sr = 1.0;
00659 }
00660
00661
00662
00663
00664 if (wave_speed_fail){
00665 getMaxSignalSpeeds_echo (Wl,Wr,Bx,&Sla,&Sra);
00666
00667 if (Sla != Sla) {
00668 printf("[hlle_sr_mhd]: NaN in Sl %10.4e %10.4e\n",Sl,Sr);
00669 Sla = -1.0;
00670 Sra = 1.0;
00671 }
00672
00673 if (Sra != Sra) {
00674 printf("[hlle_sr_mhd]: NaN in Sr %10.4e %10.4e\n",Sl,Sr);
00675 Sla = -1.0;
00676 Sra = 1.0;
00677 }
00678
00679 if (Sla < -1.0) {
00680 printf("[hlle_sr_mhd]: Superluminal Sl %10.4e %10.4e\n",Sl,Sr);
00681 Sla = -1.0;
00682 Sra = 1.0;
00683 }
00684 if (Sra > 1.0) {
00685 printf("[hlle_sr_mhd]: Superluminal Sr %10.4e %10.4e\n",Sl,Sr);
00686 Sla = -1.0;
00687 Sra = 1.0;
00688 }
00689
00690 Sl = Sla;
00691 Sr = Sra;
00692
00693 }
00694
00695
00696 WSl = Wl.P*pow(Wl.d,1.0-Gamma);
00697 WSr = Wr.P*pow(Wr.d,1.0-Gamma);
00698 USl = WSl * Ul.d/Wl.d;
00699 USr = WSr * Ur.d/Wr.d;
00700 Fl = USl * Wl.Vx;
00701 Fr = USr * Wr.Vx;
00702
00703 if(Sl >= 0.0){
00704 *pFlux = Fl;
00705 return;
00706 }
00707 else if(Sr <= 0.0){
00708 *pFlux = Fr;
00709 return;
00710 }
00711 else{
00712
00713 dS_1 = 1.0/(Sr - Sl);
00714 *pFlux = (Sr*Fl - Sl*Fr + Sl*Sr*(USr - USl)) * dS_1;
00715 return;
00716 }
00717 }
00718
00719
00720
00721 void flux_LR(Cons1DS U, Prim1DS W, Cons1DS *flux, Real Bx, Real* p)
00722 {
00723 Real wtg2, pt, g, g2, g_1,g_2, h, gmmr, theta;
00724 Real bx, by, bz, vB, b2, Bmag2;
00725
00726
00727
00728 theta = W.P/W.d;
00729 gmmr = Gamma / Gamma_1;
00730
00731 h = 1.0 + gmmr*theta;
00732
00733
00734 g = U.d/W.d;
00735 g2 = SQR(g);
00736 g_1 = 1.0/g;
00737 g_2 = 1.0/g2;
00738
00739 pt = W.P;
00740 wtg2 = W.d*h*g2;
00741
00742 vB = W.Vx*Bx + W.Vy*W.By + W.Vz*W.Bz;
00743 Bmag2 = SQR(Bx) + SQR(W.By) + SQR(W.Bz);
00744
00745 bx = g*( Bx*g_2 + vB*W.Vx);
00746 by = g*(W.By*g_2 + vB*W.Vy);
00747 bz = g*(W.Bz*g_2 + vB*W.Vz);
00748
00749 b2 = Bmag2*g_2 + vB*vB;
00750
00751 pt += 0.5*b2;
00752 wtg2 += b2*g2;
00753
00754 flux->d = U.d*W.Vx;
00755 flux->Mx = wtg2*W.Vx*W.Vx + pt;
00756 flux->My = wtg2*W.Vy*W.Vx;
00757 flux->Mz = wtg2*W.Vz*W.Vx;
00758 flux->E = U.Mx;
00759
00760 flux->Mx -= bx*bx;
00761 flux->My -= by*bx;
00762 flux->Mz -= bz*bx;
00763 flux->By = W.Vx*W.By - Bx*W.Vy;
00764 flux->Bz = W.Vx*W.Bz - Bx*W.Vz;
00765
00766 *p = pt;
00767 }
00768
00769
00770
00771
00772
00773 void getMaxSignalSpeeds_pluto(const Prim1DS Wl, const Prim1DS Wr,
00774 const Real Bx, Real* low, Real* high)
00775 {
00776
00777 Real lml,lmr;
00778 Real lpl,lpr;
00779 Real al,ar;
00780
00781 getVChar_pluto(Wl,Bx,&lml,&lpl);
00782 getVChar_pluto(Wr,Bx,&lmr,&lpr);
00783
00784 *low = MIN(lml, lmr);
00785 *high = MAX(lpl, lpr);
00786 }
00787
00788
00789
00790
00791 void getVChar_pluto(const Prim1DS W, const Real Bx, Real* lm, Real* lp)
00792 {
00793 Real rhoh,vsq,bsq;
00794 Real cssq,vasq,asq;
00795 Real Vx2,gamma,gamma2;
00796 Real Bx2,Bsq,vDotB,vDotBsq,b0,bx;
00797 Real w_1,a0,a1,a2,a3,a4,Q;
00798 Real scrh,scrh1,scrh2,eps2;
00799 Real a2_w,one_m_eps2,lambda[5];
00800 int iflag;
00801
00802 rhoh = W.d + (Gamma/Gamma_1) * (W.P);
00803
00804 Vx2 = SQR(W.Vx);
00805 vsq = Vx2 + SQR(W.Vy) + SQR(W.Vz);
00806 if (vsq > 1.0){
00807
00808 *lm = -1.0;
00809 *lp = 1.0;
00810 return;
00811 }
00812 gamma = 1.0 / sqrt(1 - vsq);
00813 gamma2 = SQR(gamma);
00814
00815 Bx2 = SQR(Bx);
00816 Bsq = Bx2 + SQR(W.By) + SQR(W.Bz);
00817 vDotB = W.Vx*Bx + W.Vy*W.By + W.Vz*W.Bz;
00818 vDotBsq = SQR(vDotB);
00819 b0 = gamma * vDotB;
00820 bx = Bx/gamma2 + W.Vx*vDotB;
00821 bsq = Bsq / gamma2 + SQR(vDotB);
00822
00823 cssq = (Gamma * W.P) / (rhoh);
00824 vasq = bsq / (rhoh + bsq);
00825
00826 if (bsq < 0.0) bsq = 0.0;
00827 if (cssq < 0.0) cssq = 0.0;
00828 if (cssq > 1.0) cssq = 1.0;
00829 if (vasq > 1.0) bsq = rhoh + bsq;
00830
00831 if (vsq < 1.0e-12) {
00832 w_1 = 1.0/(rhoh + bsq);
00833 eps2 = cssq + bsq*w_1*(1.0 - cssq);
00834 a0 = cssq*Bx*Bx*w_1;
00835 a1 = - a0 - eps2;
00836 scrh = a1*a1 - 4.0*a0;
00837 if (scrh < 0.0) scrh = 0.0;
00838
00839 scrh = sqrt(0.5*(-a1 + sqrt(scrh)));
00840 *lp = scrh;
00841 *lm = -scrh;
00842 return;
00843 }
00844
00845 w_1 = 1.0/(rhoh + bsq);
00846
00847 if (Bx < 1.0e-14) {
00848
00849 eps2 = cssq + bsq*w_1*(1.0 - cssq);
00850
00851 scrh1 = (1.0 - eps2)*gamma2;
00852 scrh2 = cssq*vDotBsq*w_1 - eps2;
00853
00854 a2 = scrh1 - scrh2;
00855 a1 = -2.0*W.Vx*scrh1;
00856 a0 = Vx2*scrh1 + scrh2;
00857
00858 *lp = 0.5*(-a1 + sqrt(a1*a1 - 4.0*a2*a0))/a2;
00859 *lm = 0.5*(-a1 - sqrt(a1*a1 - 4.0*a2*a0))/a2;
00860
00861 return;
00862 }
00863
00864 scrh1 = bx;
00865 scrh2 = scrh1*scrh1;
00866
00867 a2_w = cssq*w_1;
00868 eps2 = (cssq*rhoh + bsq)*w_1;
00869 one_m_eps2 = gamma2*rhoh*(1.0 - cssq)*w_1;
00870
00871
00872
00873
00874
00875 scrh = 2.0*(a2_w*vDotB*scrh1 - eps2*W.Vx);
00876 a4 = one_m_eps2 - a2_w*vDotBsq + eps2;
00877 a3 = - 4.0*W.Vx*one_m_eps2 + scrh;
00878 a2 = 6.0*Vx2*one_m_eps2 + a2_w*(vDotBsq - scrh2) + eps2*(Vx2 - 1.0);
00879 a1 = - 4.0*W.Vx*Vx2*one_m_eps2 - scrh;
00880 a0 = Vx2*Vx2*one_m_eps2 + a2_w*scrh2 - eps2*Vx2;
00881
00882 if (a4 < 1.e-12){
00883
00884 printf("[MAX_CH_SPEED]: Can not divide by a4 in MAX_CH_SPEED\n");
00885
00886 *lm = -1.0;
00887 *lp = 1.0;
00888
00889 return;
00890 }
00891
00892 scrh = 1.0/a4;
00893
00894 a3 *= scrh;
00895 a2 *= scrh;
00896 a1 *= scrh;
00897 a0 *= scrh;
00898 iflag = QUARTIC(a3, a2, a1, a0, lambda);
00899
00900 if (iflag){
00901 printf ("Can not find max speed:\n");
00902
00903 printf("QUARTIC: f(x) = %12.6e + x*(%12.6e + x*(%12.6e ",
00904 a0*a4, a1*a4, a2*a4);
00905 printf("+ x*(%12.6e + x*%12.6e)))\n", a3*a4, a4);
00906 printf("[MAX_CH_SPEED]: Failed to find wave speeds");
00907
00908 *lm = -1.0;
00909 *lp = 1.0;
00910
00911 return;
00912 }
00913
00914 *lp = MIN(1.0,MAX(lambda[3], lambda[2]));
00915 *lp = MIN(1.0,MAX(*lp, lambda[1]));
00916 *lp = MIN(1.0,MAX(*lp, lambda[0]));
00917
00918 *lm = MAX(-1.0,MIN(lambda[3], lambda[2]));
00919 *lm = MAX(-1.0,MIN(*lm, lambda[1]));
00920 *lm = MAX(-1.0,MIN(*lm, lambda[0]));
00921
00922 return;
00923
00924 }
00925
00926
00927
00928
00929
00930 void getMaxSignalSpeeds_echo (const Prim1DS Wl, const Prim1DS Wr,
00931 const Real Bx, Real* low, Real* high)
00932 {
00933
00934 Real lml,lmr;
00935 Real lpl,lpr;
00936 Real al,ar;
00937
00938 getVChar_echo(Wl,Bx,&lml,&lpl);
00939 getVChar_echo(Wr,Bx,&lmr,&lpr);
00940
00941 *low = MIN(lml, lmr);
00942 *high = MAX(lpl, lpr);
00943 }
00944
00945
00946
00947
00948 void getVChar_echo(const Prim1DS W, const Real Bx, Real* lm, Real* lp)
00949 {
00950 Real rhoh,vsq,bsq;
00951 Real cssq,vasq,asq;
00952 Real gamma,gamma2;
00953 Real Bsq,vDotB,b0,bx;
00954 Real tmp1,tmp2,tmp3,tmp4,tmp5;
00955 Real vm,vp;
00956
00957 rhoh = W.d + (Gamma/Gamma_1) * (W.P);
00958
00959 vsq = SQR(W.Vx) + SQR(W.Vy) + SQR(W.Vz);
00960 gamma = 1.0 / sqrt(1 - vsq);
00961 gamma2 = SQR(gamma);
00962
00963 Bsq = SQR(Bx) + SQR(W.By) + SQR(W.Bz);
00964 vDotB = W.Vx*Bx + W.Vy*W.By + W.Vz*W.Bz;
00965 b0 = gamma * vDotB;
00966 bx = Bx/gamma2 + W.Vx*vDotB;
00967 bsq = Bsq / gamma2 + SQR(vDotB);
00968
00969 cssq = (Gamma * W.P) / (rhoh);
00970 vasq = bsq / (rhoh + bsq);
00971 asq = cssq + vasq - (cssq*vasq);
00972
00973 if (cssq < 0.0) cssq = 0.0;
00974 if (vasq > 0.0) vasq = 0.0;
00975 if (asq < 0.0) asq = 0.0;
00976 if (cssq > 1.0) cssq = 1.0;
00977 if (vasq > 1.0) vasq = 1.0;
00978 if (asq > 1.0) asq = 1.0;
00979
00980 tmp1 = (1.0 - asq);
00981 tmp2 = (1.0 - vsq);
00982 tmp3 = (1.0 - vsq*asq);
00983 tmp4 = SQR(W.Vx);
00984 tmp5 = 1.0 / tmp3;
00985
00986 vm = tmp1*W.Vx - sqrt(asq*tmp2*(tmp3 - tmp1*tmp4));
00987 vp = tmp1*W.Vx + sqrt(asq*tmp2*(tmp3 - tmp1*tmp4));
00988 vm *=tmp5;
00989 vp *=tmp5;
00990
00991 if (vp > vm) {
00992 *lm = vm;
00993 *lp = vp;
00994 } else {
00995 *lm = vp;
00996 *lp = vm;
00997 }
00998 }
00999
01000
01001
01002
01003
01004
01005
01006
01007
01008
01009
01010
01011
01012
01013
01014
01015
01016
01017
01018
01019
01020
01021
01022
01023
01024
01025
01026
01027
01028
01029
01030
01031
01032 int QUARTIC (Real b, Real c, Real d,
01033 Real e, Real z[])
01034 {
01035 int n, ifail;
01036 Real b2, f, g, h;
01037 Real a2, a1, a0, u[4];
01038 Real p, q, r, s;
01039 static Real three_256 = 3.0/256.0;
01040 static Real one_64 = 1.0/64.0;
01041
01042 b2 = b*b;
01043
01044 f = c - b2*0.375;
01045 g = d + b2*b*0.125 - b*c*0.5;
01046 h = e - b2*b2*three_256 + 0.0625*b2*c - 0.25*b*d;
01047
01048 a2 = 0.5*f;
01049 a1 = (f*f - 4.0*h)*0.0625;
01050 a0 = -g*g*one_64;
01051
01052 ifail = CUBIC(a2, a1, a0, u);
01053
01054 if (ifail)return(1);
01055
01056 if (u[1] < 1.e-14){
01057
01058 p = sqrt(u[2]);
01059 s = 0.25*b;
01060 z[0] = z[2] = - p - s;
01061 z[1] = z[3] = + p - s;
01062
01063 }else{
01064
01065 p = sqrt(u[1]);
01066 q = sqrt(u[2]);
01067
01068 r = -0.125*g/(p*q);
01069 s = 0.25*b;
01070
01071 z[0] = - p - q + r - s;
01072 z[1] = p - q - r - s;
01073 z[2] = - p + q - r - s;
01074 z[3] = p + q + r - s;
01075
01076 }
01077
01078
01079
01080
01081
01082
01083 for (n = 0; n < 4; n++){
01084 s = e + z[n]*(d + z[n]*(c + z[n]*(b + z[n])));
01085 if (s != s) {
01086 printf ("Nan found in QUARTIC \n");
01087 return(1);
01088 }
01089 if (fabs(s) > 1.e-6) {
01090 printf ("Solution does not satisfy f(z) = 0; f(z) = %12.6e\n",s);
01091 return(1);
01092 }
01093 }
01094
01095 return(0);
01096
01097
01098
01099 }
01100
01101
01102
01103
01104
01105
01106
01107
01108
01109
01110
01111
01112
01113
01114
01115
01116
01117
01118
01119
01120
01121
01122
01123
01124
01125
01126
01127
01128
01129
01130
01131
01132 int CUBIC(Real b, Real c, Real d, Real z[])
01133 {
01134 Real b2, g2;
01135 Real f, g, h;
01136 Real i, i2, j, k, m, n, p;
01137 static Real one_3 = 1.0/3.0, one_27=1.0/27.0;
01138
01139 b2 = b*b;
01140
01141
01142
01143
01144
01145
01146
01147
01148 f = c*(1.0 - 1.e-16) - b2*one_3;
01149 g = b*(2.0*b2 - 9.0*c)*one_27 + d;
01150 g2 = g*g;
01151 i2 = -f*f*f*one_27;
01152 h = g2*0.25 - i2;
01153
01154
01155
01156
01157
01158
01159
01160 if (h > 1.e-12){
01161 printf ("Only one real root (%12.6e)!\n", h);
01162 }
01163 if (i2 < 0.0){
01164
01165
01166
01167
01168 i2 = 0.0;
01169 }
01170
01171
01172
01173
01174
01175 i = sqrt(i2);
01176 j = pow(i, one_3);
01177 k = -0.5*g/i;
01178
01179
01180
01181
01182 k = (k < -1.0 ? -1.0:k);
01183 k = (k > 1.0 ? 1.0:k);
01184
01185 k = acos(k)*one_3;
01186
01187 m = cos(k);
01188 n = sqrt(3.0)*sin(k);
01189 p = -b*one_3;
01190
01191 z[0] = -j*(m + n) + p;
01192 z[1] = -j*(m - n) + p;
01193 z[2] = 2.0*j*m + p;
01194
01195
01196
01197
01198
01199
01200
01201
01202
01203
01204
01205
01206 return(0);
01207 }
01208
01209 #endif
01210
01211 #endif
01212 #endif