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 #include <math.h>
00034 #include <float.h>
00035 #include "../defs.h"
00036 #include "../athena.h"
00037 #include "../globals.h"
00038 #include "prototypes.h"
00039 #include "../prototypes.h"
00040
00041
00042
00043
00044
00045
00046
00047
00048
00049
00050 #if defined(BAROTROPIC) && defined(HYDRO)
00051 void esys_prim_iso_hyd(const Real d, const Real v1,
00052 Real eigenvalues[],
00053 Real right_eigenmatrix[][4], Real left_eigenmatrix[][4])
00054 {
00055
00056
00057
00058 eigenvalues[0] = v1 - Iso_csound;
00059 eigenvalues[1] = v1;
00060 eigenvalues[2] = v1;
00061 eigenvalues[3] = v1 + Iso_csound;
00062
00063
00064
00065 right_eigenmatrix[0][0] = 1.0;
00066 right_eigenmatrix[1][0] = -Iso_csound/d;
00067
00068
00069
00070
00071
00072 right_eigenmatrix[2][1] = 1.0;
00073
00074
00075
00076
00077
00078 right_eigenmatrix[3][2] = 1.0;
00079
00080 right_eigenmatrix[0][3] = 1.0;
00081 right_eigenmatrix[1][3] = Iso_csound/d;
00082
00083
00084
00085
00086
00087 left_eigenmatrix[0][0] = 0.5;
00088 left_eigenmatrix[0][1] = -0.5*d/Iso_csound;
00089
00090
00091
00092
00093
00094 left_eigenmatrix[1][2] = 1.0;
00095
00096
00097
00098
00099
00100 left_eigenmatrix[2][3] = 1.0;
00101
00102 left_eigenmatrix[3][0] = 0.5;
00103 left_eigenmatrix[3][1] = 0.5*d/Iso_csound;
00104
00105
00106
00107 }
00108 #endif
00109
00110
00111
00112
00113
00114
00115
00116
00117
00118
00119 #if !defined(BAROTROPIC) && defined(HYDRO)
00120 void esys_prim_adb_hyd(const Real d, const Real v1, const Real rho_a2,
00121 Real eigenvalues[],
00122 Real right_eigenmatrix[][5], Real left_eigenmatrix[][5])
00123 {
00124 Real asq,a;
00125 asq = rho_a2/d;
00126 a = sqrt(asq);
00127
00128
00129
00130 eigenvalues[0] = v1 - a;
00131 eigenvalues[1] = v1;
00132 eigenvalues[2] = v1;
00133 eigenvalues[3] = v1;
00134 eigenvalues[4] = v1 + a;
00135
00136
00137
00138 right_eigenmatrix[0][0] = 1.0;
00139 right_eigenmatrix[1][0] = -a/d;
00140
00141
00142 right_eigenmatrix[4][0] = asq;
00143
00144 right_eigenmatrix[0][1] = 1.0;
00145
00146
00147
00148
00149
00150
00151
00152 right_eigenmatrix[2][2] = 1.0;
00153
00154
00155
00156
00157
00158
00159 right_eigenmatrix[3][3] = 1.0;
00160
00161
00162 right_eigenmatrix[0][4] = 1.0;
00163 right_eigenmatrix[1][4] = -right_eigenmatrix[1][0];
00164
00165
00166 right_eigenmatrix[4][4] = asq;
00167
00168
00169
00170
00171 left_eigenmatrix[0][1] = -0.5*d/a;
00172
00173
00174 left_eigenmatrix[0][4] = 0.5/asq;
00175
00176 left_eigenmatrix[1][0] = 1.0;
00177
00178
00179
00180 left_eigenmatrix[1][4] = -1.0/asq;
00181
00182
00183
00184 left_eigenmatrix[2][2] = 1.0;
00185
00186
00187
00188
00189
00190
00191 left_eigenmatrix[3][3] = 1.0;
00192
00193
00194
00195 left_eigenmatrix[4][1] = -left_eigenmatrix[0][1];
00196
00197
00198 left_eigenmatrix[4][4] = left_eigenmatrix[0][4];
00199 }
00200 #endif
00201
00202
00203
00204
00205
00206
00207
00208
00209
00210
00211 #if defined(BAROTROPIC) && defined(MHD)
00212 void esys_prim_iso_mhd(const Real d, const Real v1, const Real b1,
00213 const Real b2, const Real b3, Real eigenvalues[],
00214 Real right_eigenmatrix[][6], Real left_eigenmatrix[][6])
00215 {
00216 Real btsq,vaxsq,cfsq,cf,cssq,cs,bt,bet2,bet3,alpha_f,alpha_s;
00217 Real sqrtd,s,qf,qs,af,as,vax,norm,af_prime,as_prime;
00218 Real ct2,tsum,tdif,cf2_cs2;
00219 Real di = 1.0/d;
00220 btsq = b2*b2 + b3*b3;
00221 vaxsq = b1*b1*di;
00222
00223
00224
00225 ct2 = btsq*di;
00226 tsum = vaxsq + ct2 + (Iso_csound2);
00227 tdif = vaxsq + ct2 - (Iso_csound2);
00228 cf2_cs2 = sqrt((double)(tdif*tdif + 4.0*(Iso_csound2)*ct2));
00229
00230 cfsq = 0.5*(tsum + cf2_cs2);
00231 cf = sqrt((double)cfsq);
00232
00233 cssq = (Iso_csound2)*vaxsq/cfsq;
00234 cs = sqrt((double)cssq);
00235
00236
00237
00238 bt = sqrt(btsq);
00239 if (bt == 0.0) {
00240 bet2 = 1.0;
00241 bet3 = 0.0;
00242 } else {
00243 bet2 = b2/bt;
00244 bet3 = b3/bt;
00245 }
00246
00247
00248
00249 if ((cfsq-cssq) == 0.0) {
00250 alpha_f = 1.0;
00251 alpha_s = 0.0;
00252 } else if ( (Iso_csound2 - cssq) <= 0.0) {
00253 alpha_f = 0.0;
00254 alpha_s = 1.0;
00255 } else if ( (cfsq - Iso_csound2) <= 0.0) {
00256 alpha_f = 1.0;
00257 alpha_s = 0.0;
00258 } else {
00259 alpha_f = sqrt((Iso_csound2 - cssq)/(cfsq - cssq));
00260 alpha_s = sqrt((cfsq - Iso_csound2)/(cfsq - cssq));
00261 }
00262
00263
00264
00265 sqrtd = sqrt(d);
00266 s = SIGN(b1);
00267 qf = cf*alpha_f*s;
00268 qs = cs*alpha_s*s;
00269 af = Iso_csound*alpha_f*sqrtd;
00270 as = Iso_csound*alpha_s*sqrtd;
00271
00272
00273
00274 vax = sqrt(vaxsq);
00275 eigenvalues[0] = v1 - cf;
00276 eigenvalues[1] = v1 - vax;
00277 eigenvalues[2] = v1 - cs;
00278 eigenvalues[3] = v1 + cs;
00279 eigenvalues[4] = v1 + vax;
00280 eigenvalues[5] = v1 + cf;
00281
00282
00283
00284 right_eigenmatrix[0][0] = d*alpha_f;
00285 right_eigenmatrix[1][0] = -cf*alpha_f;
00286 right_eigenmatrix[2][0] = qs*bet2;
00287 right_eigenmatrix[3][0] = qs*bet3;
00288 right_eigenmatrix[4][0] = as*bet2;
00289 right_eigenmatrix[5][0] = as*bet3;
00290
00291
00292
00293 right_eigenmatrix[2][1] = -bet3;
00294 right_eigenmatrix[3][1] = bet2;
00295 right_eigenmatrix[4][1] = -bet3*s*sqrtd;
00296 right_eigenmatrix[5][1] = bet2*s*sqrtd;
00297
00298 right_eigenmatrix[0][2] = d*alpha_s;
00299 right_eigenmatrix[1][2] = -cs*alpha_s;
00300 right_eigenmatrix[2][2] = -qf*bet2;
00301 right_eigenmatrix[3][2] = -qf*bet3;
00302 right_eigenmatrix[4][2] = -af*bet2;
00303 right_eigenmatrix[5][2] = -af*bet3;
00304
00305 right_eigenmatrix[0][3] = right_eigenmatrix[0][2];
00306 right_eigenmatrix[1][3] = -right_eigenmatrix[1][2];
00307 right_eigenmatrix[2][3] = -right_eigenmatrix[2][2];
00308 right_eigenmatrix[3][3] = -right_eigenmatrix[3][2];
00309 right_eigenmatrix[4][3] = right_eigenmatrix[4][2];
00310 right_eigenmatrix[5][3] = right_eigenmatrix[5][2];
00311
00312
00313
00314 right_eigenmatrix[2][4] = bet3;
00315 right_eigenmatrix[3][4] = -bet2;
00316 right_eigenmatrix[4][4] = right_eigenmatrix[4][1];
00317 right_eigenmatrix[5][4] = right_eigenmatrix[5][1];
00318
00319 right_eigenmatrix[0][5] = right_eigenmatrix[0][0];
00320 right_eigenmatrix[1][5] = -right_eigenmatrix[1][0];
00321 right_eigenmatrix[2][5] = -right_eigenmatrix[2][0];
00322 right_eigenmatrix[3][5] = -right_eigenmatrix[3][0];
00323 right_eigenmatrix[4][5] = right_eigenmatrix[4][0];
00324 right_eigenmatrix[5][5] = right_eigenmatrix[5][0];
00325
00326
00327
00328 norm = 0.5/Iso_csound2;
00329 qf = norm*qf;
00330 qs = norm*qs;
00331 af_prime = norm*af*di;
00332 as_prime = norm*as*di;
00333
00334 left_eigenmatrix[0][0] = norm*alpha_f*Iso_csound2*di;
00335 left_eigenmatrix[0][1] = -norm*cf*alpha_f;
00336 left_eigenmatrix[0][2] = qs*bet2;
00337 left_eigenmatrix[0][3] = qs*bet3;
00338 left_eigenmatrix[0][4] = as_prime*bet2;
00339 left_eigenmatrix[0][5] = as_prime*bet3;
00340
00341
00342
00343 left_eigenmatrix[1][2] = -0.5*bet3;
00344 left_eigenmatrix[1][3] = 0.5*bet2;
00345 left_eigenmatrix[1][4] = -0.5*bet3*s/sqrtd;
00346 left_eigenmatrix[1][5] = 0.5*bet2*s/sqrtd;
00347
00348 left_eigenmatrix[2][0] = norm*alpha_s*Iso_csound2*di;
00349 left_eigenmatrix[2][1] = -norm*cs*alpha_s;
00350 left_eigenmatrix[2][2] = -qf*bet2;
00351 left_eigenmatrix[2][3] = -qf*bet3;
00352 left_eigenmatrix[2][4] = -af_prime*bet2;
00353 left_eigenmatrix[2][5] = -af_prime*bet3;
00354
00355 left_eigenmatrix[3][0] = left_eigenmatrix[2][0];
00356 left_eigenmatrix[3][1] = -left_eigenmatrix[2][1];
00357 left_eigenmatrix[3][2] = -left_eigenmatrix[2][2];
00358 left_eigenmatrix[3][3] = -left_eigenmatrix[2][3];
00359 left_eigenmatrix[3][4] = left_eigenmatrix[2][4];
00360 left_eigenmatrix[3][5] = left_eigenmatrix[2][5];
00361
00362
00363
00364 left_eigenmatrix[4][2] = -left_eigenmatrix[1][2];
00365 left_eigenmatrix[4][3] = -left_eigenmatrix[1][3];
00366 left_eigenmatrix[4][4] = left_eigenmatrix[1][4];
00367 left_eigenmatrix[4][5] = left_eigenmatrix[1][5];
00368
00369 left_eigenmatrix[5][0] = left_eigenmatrix[0][0];
00370 left_eigenmatrix[5][1] = -left_eigenmatrix[0][1];
00371 left_eigenmatrix[5][2] = -left_eigenmatrix[0][2];
00372 left_eigenmatrix[5][3] = -left_eigenmatrix[0][3];
00373 left_eigenmatrix[5][4] = left_eigenmatrix[0][4];
00374 left_eigenmatrix[5][5] = left_eigenmatrix[0][5];
00375 }
00376 #endif
00377
00378
00379
00380
00381
00382
00383
00384
00385
00386
00387 #if !defined(BAROTROPIC) && defined(MHD)
00388 void esys_prim_adb_mhd(const Real d, const Real v1, const Real rho_a2,
00389 const Real b1, const Real b2, const Real b3, Real eigenvalues[],
00390 Real right_eigenmatrix[][7], Real left_eigenmatrix[][7])
00391 {
00392 Real di,cfsq,cf,cssq,cs,bt,bet2,bet3,alpha_f,alpha_s;
00393 Real sqrtd,s,a,qf,qs,af,as,vax,na,af_prime,as_prime;
00394 Real tsum,tdif,cf2_cs2,ct2;
00395 Real btsq,vaxsq,asq;
00396 di = 1.0/d;
00397 btsq = b2*b2 + b3*b3;
00398 vaxsq = b1*b1*di;
00399 asq = rho_a2*di;
00400
00401
00402
00403 ct2 = btsq*di;
00404 tsum = vaxsq + ct2 + asq;
00405 tdif = vaxsq + ct2 - asq;
00406 cf2_cs2 = sqrt((double)(tdif*tdif + 4.0*asq*ct2));
00407
00408 cfsq = 0.5*(tsum + cf2_cs2);
00409 cf = sqrt((double)cfsq);
00410
00411 cssq = asq*vaxsq/cfsq;
00412 cs = sqrt((double)cssq);
00413
00414
00415
00416 bt = sqrt(btsq);
00417 if (bt == 0.0) {
00418 bet2 = 1.0;
00419 bet3 = 0.0;
00420 } else {
00421 bet2 = b2/bt;
00422 bet3 = b3/bt;
00423 }
00424
00425
00426
00427 if (cf2_cs2 == 0.0) {
00428 alpha_f = 1.0;
00429 alpha_s = 0.0;
00430 } else if ( (asq - cssq) <= 0.0) {
00431 alpha_f = 0.0;
00432 alpha_s = 1.0;
00433 } else if ( (cfsq - asq) <= 0.0) {
00434 alpha_f = 1.0;
00435 alpha_s = 0.0;
00436 } else {
00437 alpha_f = sqrt((asq - cssq)/cf2_cs2);
00438 alpha_s = sqrt((cfsq - asq)/cf2_cs2);
00439 }
00440
00441
00442
00443 sqrtd = sqrt(d);
00444 s = SIGN(b1);
00445 a = sqrt(asq);
00446 qf = cf*alpha_f*s;
00447 qs = cs*alpha_s*s;
00448 af = a*alpha_f*sqrtd;
00449 as = a*alpha_s*sqrtd;
00450
00451
00452
00453 vax = sqrt(vaxsq);
00454 eigenvalues[0] = v1 - cf;
00455 eigenvalues[1] = v1 - vax;
00456 eigenvalues[2] = v1 - cs;
00457 eigenvalues[3] = v1;
00458 eigenvalues[4] = v1 + cs;
00459 eigenvalues[5] = v1 + vax;
00460 eigenvalues[6] = v1 + cf;
00461
00462
00463
00464
00465
00466
00467 right_eigenmatrix[0][0] = d*alpha_f;
00468
00469 right_eigenmatrix[0][2] = d*alpha_s;
00470 right_eigenmatrix[0][3] = 1.0;
00471 right_eigenmatrix[0][4] = right_eigenmatrix[0][2];
00472
00473 right_eigenmatrix[0][6] = right_eigenmatrix[0][0];
00474
00475 right_eigenmatrix[1][0] = -cf*alpha_f;
00476
00477 right_eigenmatrix[1][2] = -cs*alpha_s;
00478
00479 right_eigenmatrix[1][4] = -right_eigenmatrix[1][2];
00480
00481 right_eigenmatrix[1][6] = -right_eigenmatrix[1][0];
00482
00483 right_eigenmatrix[2][0] = qs*bet2;
00484 right_eigenmatrix[2][1] = -bet3;
00485 right_eigenmatrix[2][2] = -qf*bet2;
00486
00487 right_eigenmatrix[2][4] = -right_eigenmatrix[2][2];
00488 right_eigenmatrix[2][5] = bet3;
00489 right_eigenmatrix[2][6] = -right_eigenmatrix[2][0];
00490
00491 right_eigenmatrix[3][0] = qs*bet3;
00492 right_eigenmatrix[3][1] = bet2;
00493 right_eigenmatrix[3][2] = -qf*bet3;
00494
00495 right_eigenmatrix[3][4] = -right_eigenmatrix[3][2];
00496 right_eigenmatrix[3][5] = -bet2;
00497 right_eigenmatrix[3][6] = -right_eigenmatrix[3][0];
00498
00499 right_eigenmatrix[4][0] = d*asq*alpha_f;
00500
00501 right_eigenmatrix[4][2] = d*asq*alpha_s;
00502
00503 right_eigenmatrix[4][4] = right_eigenmatrix[4][2];
00504
00505 right_eigenmatrix[4][6] = right_eigenmatrix[4][0];
00506
00507 right_eigenmatrix[5][0] = as*bet2;
00508 right_eigenmatrix[5][1] = -bet3*s*sqrtd;
00509 right_eigenmatrix[5][2] = -af*bet2;
00510
00511 right_eigenmatrix[5][4] = right_eigenmatrix[5][2];
00512 right_eigenmatrix[5][5] = right_eigenmatrix[5][1];
00513 right_eigenmatrix[5][6] = right_eigenmatrix[5][0];
00514
00515 right_eigenmatrix[6][0] = as*bet3;
00516 right_eigenmatrix[6][1] = bet2*s*sqrtd;
00517 right_eigenmatrix[6][2] = -af*bet3;
00518
00519 right_eigenmatrix[6][4] = right_eigenmatrix[6][2];
00520 right_eigenmatrix[6][5] = right_eigenmatrix[6][1];
00521 right_eigenmatrix[6][6] = right_eigenmatrix[6][0];
00522
00523
00524
00525 na = 0.5/asq;
00526 qf = na*qf;
00527 qs = na*qs;
00528 af_prime = na*af*di;
00529 as_prime = na*as*di;
00530
00531
00532 left_eigenmatrix[0][1] = -na*cf*alpha_f;
00533 left_eigenmatrix[0][2] = qs*bet2;
00534 left_eigenmatrix[0][3] = qs*bet3;
00535 left_eigenmatrix[0][4] = na*alpha_f*di;
00536 left_eigenmatrix[0][5] = as_prime*bet2;
00537 left_eigenmatrix[0][6] = as_prime*bet3;
00538
00539
00540
00541 left_eigenmatrix[1][2] = -0.5*bet3;
00542 left_eigenmatrix[1][3] = 0.5*bet2;
00543
00544 left_eigenmatrix[1][5] = -0.5*bet3*s/sqrtd;
00545 left_eigenmatrix[1][6] = 0.5*bet2*s/sqrtd;
00546
00547
00548 left_eigenmatrix[2][1] = -na*cs*alpha_s;
00549 left_eigenmatrix[2][2] = -qf*bet2;
00550 left_eigenmatrix[2][3] = -qf*bet3;
00551 left_eigenmatrix[2][4] = na*alpha_s*di;
00552 left_eigenmatrix[2][5] = -af_prime*bet2;
00553 left_eigenmatrix[2][6] = -af_prime*bet3;
00554
00555 left_eigenmatrix[3][0] = 1.0;
00556
00557
00558
00559 left_eigenmatrix[3][4] = -1.0/asq;
00560
00561
00562
00563
00564 left_eigenmatrix[4][1] = -left_eigenmatrix[2][1];
00565 left_eigenmatrix[4][2] = -left_eigenmatrix[2][2];
00566 left_eigenmatrix[4][3] = -left_eigenmatrix[2][3];
00567 left_eigenmatrix[4][4] = left_eigenmatrix[2][4];
00568 left_eigenmatrix[4][5] = left_eigenmatrix[2][5];
00569 left_eigenmatrix[4][6] = left_eigenmatrix[2][6];
00570
00571
00572
00573 left_eigenmatrix[5][2] = -left_eigenmatrix[1][2];
00574 left_eigenmatrix[5][3] = -left_eigenmatrix[1][3];
00575
00576 left_eigenmatrix[5][5] = left_eigenmatrix[1][5];
00577 left_eigenmatrix[5][6] = left_eigenmatrix[1][6];
00578
00579
00580 left_eigenmatrix[6][1] = -left_eigenmatrix[0][1];
00581 left_eigenmatrix[6][2] = -left_eigenmatrix[0][2];
00582 left_eigenmatrix[6][3] = -left_eigenmatrix[0][3];
00583 left_eigenmatrix[6][4] = left_eigenmatrix[0][4];
00584 left_eigenmatrix[6][5] = left_eigenmatrix[0][5];
00585 left_eigenmatrix[6][6] = left_eigenmatrix[0][6];
00586 }
00587 #endif