Go to the documentation of this file.00001 #include "../copyright.h"
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018 #include <math.h>
00019 #include <stdio.h>
00020 #include <stdlib.h>
00021 #include "../defs.h"
00022 #include "../athena.h"
00023 #include "../globals.h"
00024 #include "prototypes.h"
00025 #include "../prototypes.h"
00026
00027 #ifdef TWO_SHOCK_FLUX
00028
00029 #ifdef MHD
00030 #error : The 2-shock flux for MHD has not been implemented.
00031 #endif
00032
00033 #ifndef ISOTHERMAL
00034 #error : The 2-shock flux for adiabatic EOS has not been implemented.
00035 #endif
00036
00037 #if (NSCALARS > 0)
00038 #error : Passive scalars have not been implemented in the 2-shock flux.
00039 #endif
00040
00041
00042
00043
00044
00045
00046
00047
00048
00049
00050
00051
00052
00053 void fluxes(const Cons1DS Ul, const Cons1DS Ur,
00054 const Prim1DS Wl, const Prim1DS Wr,
00055 const Real Bxi, Cons1DS *pFlux)
00056 {
00057 Real zl, zc, zr, dc, Vxc, tmp;
00058 Real sl, sr;
00059 Real al, ar;
00060
00061 if(!(Ul.d > 0.0)||!(Ur.d > 0.0))
00062 ath_error("[two-shock flux]: Non-positive densities: dl = %e dr = %e\n",
00063 Ul.d, Ur.d);
00064
00065
00066
00067
00068
00069
00070
00071
00072
00073
00074
00075
00076
00077
00078 zl = sqrt((double)Wl.d);
00079 zr = sqrt((double)Wr.d);
00080
00081 tmp = zl*zr*(Wl.Vx - Wr.Vx)/(2.0*Iso_csound*(zl + zr));
00082 zc = tmp + sqrt((double)(tmp*tmp + zl*zr));
00083
00084 Vxc = (Wl.Vx*zl + Wr.Vx*zr)/(zl + zr) + Iso_csound*(zl - zr)/zc;
00085
00086
00087 sl = Wl.Vx - Iso_csound*zc/zl;
00088 sr = Wr.Vx + Iso_csound*zc/zr;
00089
00090
00091
00092
00093
00094 if(sr <= 0.0){
00095 pF->d = Ur.Mx;
00096 pF->Mx = Ur.Mx*(Wr.Vx) + Wr.d*Iso_csound2;
00097 pF->My = Ur.My*(Wr.Vx);
00098 pF->Mz = Ur.Mz*(Wr.Vx);
00099 return;
00100 }
00101
00102 if(sl >= 0.0){
00103 pF->d = Ul.Mx;
00104 pF->Mx = Ul.Mx*(Wl.Vx) + Wl.d*Iso_csound2;
00105 pF->My = Ul.My*(Wl.Vx);
00106 pF->Mz = Ul.Mz*(Wl.Vx);
00107 return;
00108 }
00109
00110
00111
00112
00113 dc = zc*zc;
00114 if(Vxc >= 0.0){
00115 pF->d = dc*Vxc;
00116 pF->Mx = dc*Vxc*Vxc + dc*Iso_csound2;
00117 pF->My = dc*Vxc*Wl.Vy;
00118 pF->Mz = dc*Vxc*Wl.Vz;
00119 }
00120 else{
00121 pF->d = dc*Vxc;
00122 pF->Mx = dc*Vxc*Vxc + dc*Iso_csound2;
00123 pF->My = dc*Vxc*Wr.Vy;
00124 pF->Mz = dc*Vxc*Wr.Vz;
00125 }
00126
00127 return;
00128 }
00129 #endif