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 #include <math.h>
00033 #include <stdio.h>
00034 #include <stdlib.h>
00035 #include "../defs.h"
00036 #include "../athena.h"
00037 #include "../globals.h"
00038 #include "prototypes.h"
00039 #include "../prototypes.h"
00040
00041 #ifndef SPECIAL_RELATIVITY
00042 #ifdef HLLE_FLUX
00043 #define HLLE_FUNCTION fluxes
00044 #endif
00045
00046 #if defined(HLLE_FLUX) || defined (ROE_FLUX)
00047
00048
00049
00050
00051
00052
00053
00054
00055
00056
00057
00058
00059
00060
00061
00062 void HLLE_FUNCTION(const Cons1DS Ul, const Cons1DS Ur,
00063 const Prim1DS Wl, const Prim1DS Wr,
00064 const Real Bxi, Cons1DS *pFlux)
00065 {
00066 Real sqrtdl,sqrtdr,isdlpdr,droe,v1roe,v2roe,v3roe,pbl=0.0,pbr=0.0;
00067 Real asq,vaxsq=0.0,qsq,cfsq,cfl,cfr,bp,bm,ct2=0.0,tmp;
00068 #ifndef ISOTHERMAL
00069 Real hroe;
00070 #endif
00071 #ifdef MHD
00072 Real b2roe,b3roe,x,y;
00073 #endif
00074 Real ev[NWAVE],al,ar;
00075 Real *pFl, *pFr, *pF;
00076 Cons1DS Fl,Fr;
00077 int n;
00078
00079
00080
00081
00082
00083
00084
00085
00086
00087
00088
00089
00090
00091
00092
00093 sqrtdl = sqrt((double)Wl.d);
00094 sqrtdr = sqrt((double)Wr.d);
00095 isdlpdr = 1.0/(sqrtdl + sqrtdr);
00096
00097 droe = sqrtdl*sqrtdr;
00098 v1roe = (sqrtdl*Wl.Vx + sqrtdr*Wr.Vx)*isdlpdr;
00099 v2roe = (sqrtdl*Wl.Vy + sqrtdr*Wr.Vy)*isdlpdr;
00100 v3roe = (sqrtdl*Wl.Vz + sqrtdr*Wr.Vz)*isdlpdr;
00101
00102
00103
00104 #ifdef MHD
00105 b2roe = (sqrtdr*Wl.By + sqrtdl*Wr.By)*isdlpdr;
00106 b3roe = (sqrtdr*Wl.Bz + sqrtdl*Wr.Bz)*isdlpdr;
00107 x = 0.5*(SQR(Wl.By - Wr.By) + SQR(Wl.Bz - Wr.Bz))/(SQR(sqrtdl + sqrtdr));
00108 y = 0.5*(Wl.d + Wr.d)/droe;
00109 pbl = 0.5*(SQR(Bxi) + SQR(Wl.By) + SQR(Wl.Bz));
00110 pbr = 0.5*(SQR(Bxi) + SQR(Wr.By) + SQR(Wr.Bz));
00111 #endif
00112
00113
00114
00115
00116
00117
00118 #ifndef ISOTHERMAL
00119 hroe = ((Ul.E + Wl.P + pbl)/sqrtdl + (Ur.E + Wr.P + pbr)/sqrtdr)*isdlpdr;
00120 #endif
00121
00122
00123
00124
00125
00126 #ifdef HYDRO
00127 #ifdef ISOTHERMAL
00128 esys_roe_iso_hyd(v1roe, v2roe, v3roe, ev, NULL, NULL);
00129 #else
00130 esys_roe_adb_hyd(v1roe, v2roe, v3roe, hroe, ev, NULL, NULL);
00131 #endif
00132 #endif
00133
00134 #ifdef MHD
00135 #ifdef ISOTHERMAL
00136 esys_roe_iso_mhd(droe,v1roe,v2roe,v3roe, Bxi,b2roe,b3roe,x,y,ev,NULL,NULL);
00137 #else
00138 esys_roe_adb_mhd(droe,v1roe,v2roe,v3roe,hroe,Bxi,b2roe,b3roe,x,y,ev,NULL,NULL);
00139 #endif
00140 #endif
00141
00142
00143
00144
00145
00146
00147 #ifdef ISOTHERMAL
00148 asq = Iso_csound2;
00149 #else
00150 asq = Gamma*Wl.P/Wl.d;
00151 #endif
00152 #ifdef MHD
00153 vaxsq = Bxi*Bxi/Wl.d;
00154 ct2 = (Ul.By*Ul.By + Ul.Bz*Ul.Bz)/Wl.d;
00155 #endif
00156 qsq = vaxsq + ct2 + asq;
00157 tmp = vaxsq + ct2 - asq;
00158 cfsq = 0.5*(qsq + sqrt((double)(tmp*tmp + 4.0*asq*ct2)));
00159 cfl = sqrt((double)cfsq);
00160
00161
00162 #ifdef ISOTHERMAL
00163 asq = Iso_csound2;
00164 #else
00165 asq = Gamma*Wr.P/Wr.d;
00166 #endif
00167 #ifdef MHD
00168 vaxsq = Bxi*Bxi/Wr.d;
00169 ct2 = (Ur.By*Ur.By + Ur.Bz*Ur.Bz)/Wr.d;
00170 #endif
00171 qsq = vaxsq + ct2 + asq;
00172 tmp = vaxsq + ct2 - asq;
00173 cfsq = 0.5*(qsq + sqrt((double)(tmp*tmp + 4.0*asq*ct2)));
00174 cfr = sqrt((double)cfsq);
00175
00176
00177 ar = MAX(ev[NWAVE-1],(Wr.Vx + cfr));
00178 al = MIN(ev[0] ,(Wl.Vx - cfl));
00179
00180 bp = MAX(ar, 0.0);
00181 bm = MIN(al, 0.0);
00182
00183
00184
00185
00186
00187 Fl.d = Ul.Mx - bm*Ul.d;
00188 Fr.d = Ur.Mx - bp*Ur.d;
00189
00190 Fl.Mx = Ul.Mx*(Wl.Vx - bm);
00191 Fr.Mx = Ur.Mx*(Wr.Vx - bp);
00192
00193 Fl.My = Ul.My*(Wl.Vx - bm);
00194 Fr.My = Ur.My*(Wr.Vx - bp);
00195
00196 Fl.Mz = Ul.Mz*(Wl.Vx - bm);
00197 Fr.Mz = Ur.Mz*(Wr.Vx - bp);
00198
00199 #ifdef ISOTHERMAL
00200 Fl.Mx += Wl.d*Iso_csound2;
00201 Fr.Mx += Wr.d*Iso_csound2;
00202 #else
00203 Fl.Mx += Wl.P;
00204 Fr.Mx += Wr.P;
00205
00206 Fl.E = Ul.E*(Wl.Vx - bm) + Wl.P*Wl.Vx;
00207 Fr.E = Ur.E*(Wr.Vx - bp) + Wr.P*Wr.Vx;
00208 #endif
00209
00210 #ifdef MHD
00211 Fl.Mx -= 0.5*(Bxi*Bxi - SQR(Wl.By) - SQR(Wl.Bz));
00212 Fr.Mx -= 0.5*(Bxi*Bxi - SQR(Wr.By) - SQR(Wr.Bz));
00213
00214 Fl.My -= Bxi*Wl.By;
00215 Fr.My -= Bxi*Wr.By;
00216
00217 Fl.Mz -= Bxi*Wl.Bz;
00218 Fr.Mz -= Bxi*Wr.Bz;
00219
00220 #ifndef ISOTHERMAL
00221 Fl.E += (pbl*Wl.Vx - Bxi*(Bxi*Wl.Vx + Wl.By*Wl.Vy + Wl.Bz*Wl.Vz));
00222 Fr.E += (pbr*Wr.Vx - Bxi*(Bxi*Wr.Vx + Wr.By*Wr.Vy + Wr.Bz*Wr.Vz));
00223 #endif
00224
00225 Fl.By = Wl.By*(Wl.Vx - bm) - Bxi*Wl.Vy;
00226 Fr.By = Wr.By*(Wr.Vx - bp) - Bxi*Wr.Vy;
00227
00228 Fl.Bz = Wl.Bz*(Wl.Vx - bm) - Bxi*Wl.Vz;
00229 Fr.Bz = Wr.Bz*(Wr.Vx - bp) - Bxi*Wr.Vz;
00230 #endif
00231
00232 #if (NSCALARS > 0)
00233 for (n=0; n<NSCALARS; n++) {
00234 Fl.s[n] = Fl.d*Wl.r[n];
00235 Fr.s[n] = Fr.d*Wr.r[n];
00236 }
00237 #endif
00238
00239 #ifdef CYLINDRICAL
00240 #ifndef ISOTHERMAL
00241 Fl.Pflux = Wl.P;
00242 Fr.Pflux = Wr.P;
00243 #ifdef MHD
00244 Fl.Pflux += pbl;
00245 Fr.Pflux += pbr;
00246 #endif
00247 #endif
00248 #endif
00249
00250
00251
00252
00253
00254 pFl = (Real *)&(Fl);
00255 pFr = (Real *)&(Fr);
00256 pF = (Real *)pFlux;
00257 tmp = 0.5*(bp + bm)/(bp - bm);
00258 for (n=0; n<(NWAVE+NSCALARS); n++){
00259 pF[n] = 0.5*(pFl[n] + pFr[n]) + (pFl[n] - pFr[n])*tmp;
00260 }
00261
00262 #ifdef CYLINDRICAL
00263 n = NWAVE+NSCALARS;
00264 pF[n] = 0.5*(pFl[n] + pFr[n]) + (pFl[n] - pFr[n])*tmp;
00265 #endif
00266
00267 return;
00268 }
00269 #endif
00270 #endif