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 #include <math.h>
00024 #include <stdio.h>
00025 #include <stdlib.h>
00026 #include "../defs.h"
00027 #include "../athena.h"
00028 #include "../globals.h"
00029 #include "prototypes.h"
00030 #include "../prototypes.h"
00031
00032
00033 Real etah=0.0;
00034
00035 #ifdef ROE_FLUX
00036
00037
00038
00039
00040 void flux_hlle(const Cons1DS Ul, const Cons1DS Ur,
00041 const Prim1DS Wl, const Prim1DS Wr,
00042 const Real Bxi, Cons1DS *pFlux);
00043
00044
00045 #define TEST_INTERMEDIATE_STATES
00046
00047
00048
00049
00050
00051
00052
00053
00054
00055
00056
00057
00058
00059 void fluxes(const Cons1DS Ul, const Cons1DS Ur,
00060 const Prim1DS Wl, const Prim1DS Wr,
00061 const Real Bxi, Cons1DS *pFlux)
00062 {
00063 Real sqrtdl,sqrtdr,isdlpdr,droe,v1roe,v2roe,v3roe,pbl=0.0,pbr=0.0;
00064 #ifndef ISOTHERMAL
00065 Real hroe;
00066 #endif
00067 #ifdef MHD
00068 Real b2roe,b3roe,x,y;
00069 #endif
00070 Real coeff[NWAVE];
00071 Real ev[NWAVE],rem[NWAVE][NWAVE],lem[NWAVE][NWAVE];
00072 Real dU[NWAVE],a[NWAVE];
00073 #ifdef TEST_INTERMEDIATE_STATES
00074 Real u_inter[NWAVE];
00075 #ifdef ADIABATIC
00076 Real p_inter=0.0;
00077 #endif
00078 #endif
00079 Real *pUl, *pUr, *pFl, *pFr, *pF;
00080 Cons1DS Fl,Fr;
00081 int n,m,hlle_flag;
00082 #ifdef CYLINDRICAL
00083 Real Eint,Emag,Ekin,coeff2[NWAVE];
00084 #endif
00085
00086 for (n=0; n<NWAVE; n++) {
00087 for (m=0; m<NWAVE; m++) {
00088 rem[n][m] = 0.0;
00089 lem[n][m] = 0.0;
00090 }
00091 }
00092
00093
00094
00095
00096
00097
00098
00099
00100
00101
00102
00103
00104
00105
00106 sqrtdl = sqrt((double)Wl.d);
00107 sqrtdr = sqrt((double)Wr.d);
00108 isdlpdr = 1.0/(sqrtdl + sqrtdr);
00109
00110 droe = sqrtdl*sqrtdr;
00111 v1roe = (sqrtdl*Wl.Vx + sqrtdr*Wr.Vx)*isdlpdr;
00112 v2roe = (sqrtdl*Wl.Vy + sqrtdr*Wr.Vy)*isdlpdr;
00113 v3roe = (sqrtdl*Wl.Vz + sqrtdr*Wr.Vz)*isdlpdr;
00114
00115
00116
00117 #ifdef MHD
00118 b2roe = (sqrtdr*Wl.By + sqrtdl*Wr.By)*isdlpdr;
00119 b3roe = (sqrtdr*Wl.Bz + sqrtdl*Wr.Bz)*isdlpdr;
00120 x = 0.5*(SQR(Wl.By - Wr.By) + SQR(Wl.Bz - Wr.Bz))/(SQR(sqrtdl + sqrtdr));
00121 y = 0.5*(Wl.d + Wr.d)/droe;
00122 pbl = 0.5*(SQR(Bxi) + SQR(Wl.By) + SQR(Wl.Bz));
00123 pbr = 0.5*(SQR(Bxi) + SQR(Wr.By) + SQR(Wr.Bz));
00124 #endif
00125
00126
00127
00128
00129
00130
00131 #ifndef ISOTHERMAL
00132 hroe = ((Ul.E + Wl.P + pbl)/sqrtdl + (Ur.E + Wr.P + pbr)/sqrtdr)*isdlpdr;
00133 #endif
00134
00135
00136
00137
00138
00139 #ifdef HYDRO
00140 #ifdef ISOTHERMAL
00141 esys_roe_iso_hyd(v1roe,v2roe,v3roe, ev,rem,lem);
00142 #else
00143 esys_roe_adb_hyd(v1roe,v2roe,v3roe,hroe,ev,rem,lem);
00144 #endif
00145 #endif
00146
00147 #ifdef MHD
00148 #ifdef ISOTHERMAL
00149 esys_roe_iso_mhd(droe,v1roe,v2roe,v3roe, Bxi,b2roe,b3roe,x,y,ev,rem,lem);
00150 #else
00151 esys_roe_adb_mhd(droe,v1roe,v2roe,v3roe,hroe,Bxi,b2roe,b3roe,x,y,ev,rem,lem);
00152 #endif
00153 #endif
00154
00155
00156
00157
00158
00159 Fl.d = Ul.Mx;
00160 Fr.d = Ur.Mx;
00161
00162 Fl.Mx = Ul.Mx*Wl.Vx;
00163 Fr.Mx = Ur.Mx*Wr.Vx;
00164
00165 Fl.My = Ul.Mx*Wl.Vy;
00166 Fr.My = Ur.Mx*Wr.Vy;
00167
00168 Fl.Mz = Ul.Mx*Wl.Vz;
00169 Fr.Mz = Ur.Mx*Wr.Vz;
00170
00171 #ifdef ISOTHERMAL
00172 Fl.Mx += Wl.d*Iso_csound2;
00173 Fr.Mx += Wr.d*Iso_csound2;
00174 #else
00175 Fl.Mx += Wl.P;
00176 Fr.Mx += Wr.P;
00177
00178 Fl.E = (Ul.E + Wl.P)*Wl.Vx;
00179 Fr.E = (Ur.E + Wr.P)*Wr.Vx;
00180 #endif
00181
00182 #ifdef MHD
00183 Fl.Mx -= 0.5*(Bxi*Bxi - SQR(Wl.By) - SQR(Wl.Bz));
00184 Fr.Mx -= 0.5*(Bxi*Bxi - SQR(Wr.By) - SQR(Wr.Bz));
00185
00186 Fl.My -= Bxi*Wl.By;
00187 Fr.My -= Bxi*Wr.By;
00188
00189 Fl.Mz -= Bxi*Wl.Bz;
00190 Fr.Mz -= Bxi*Wr.Bz;
00191
00192 #ifndef ISOTHERMAL
00193 Fl.E += (pbl*Wl.Vx - Bxi*(Bxi*Wl.Vx + Wl.By*Wl.Vy + Wl.Bz*Wl.Vz));
00194 Fr.E += (pbr*Wr.Vx - Bxi*(Bxi*Wr.Vx + Wr.By*Wr.Vy + Wr.Bz*Wr.Vz));
00195 #endif
00196
00197 Fl.By = Wl.By*Wl.Vx - Bxi*Wl.Vy;
00198 Fr.By = Wr.By*Wr.Vx - Bxi*Wr.Vy;
00199
00200 Fl.Bz = Wl.Bz*Wl.Vx - Bxi*Wl.Vz;
00201 Fr.Bz = Wr.Bz*Wr.Vx - Bxi*Wr.Vz;
00202 #endif
00203
00204 #if (NSCALARS > 0)
00205 for (n=0; n<NSCALARS; n++) {
00206 Fl.s[n] = Fl.d*Wl.r[n];
00207 Fr.s[n] = Fr.d*Wr.r[n];
00208 }
00209 #endif
00210
00211
00212
00213
00214
00215 if(ev[0] >= 0.0){
00216 *pFlux = Fl;
00217 #if defined(CYLINDRICAL) && !defined(BAROTROPIC)
00218 pFlux->Pflux = Wl.P;
00219 #ifdef MHD
00220 pFlux->Pflux += pbl;
00221 #endif
00222 #endif
00223 return;
00224 }
00225
00226 if(ev[NWAVE-1] <= 0.0){
00227 *pFlux = Fr;
00228 #if defined(CYLINDRICAL) && !defined(BAROTROPIC)
00229 pFlux->Pflux = Wl.P;
00230 #ifdef MHD
00231 pFlux->Pflux += pbr;
00232 #endif
00233 #endif
00234 return;
00235 }
00236
00237
00238
00239
00240
00241 pUr = (Real *)&(Ur);
00242 pUl = (Real *)&(Ul);
00243
00244 for (n=0; n<NWAVE; n++) dU[n] = pUr[n] - pUl[n];
00245 for (n=0; n<NWAVE; n++) {
00246 a[n] = 0.0;
00247 for (m=0; m<NWAVE; m++) a[n] += lem[n][m]*dU[m];
00248 }
00249
00250
00251
00252
00253
00254
00255
00256 hlle_flag = 0;
00257 #ifdef TEST_INTERMEDIATE_STATES
00258
00259 for (n=0; n<NWAVE; n++) u_inter[n] = pUl[n];
00260 for (n=0; n<NWAVE-1; n++) {
00261 for (m=0; m<NWAVE; m++) u_inter[m] += a[n]*rem[m][n];
00262 if(ev[n+1] > ev[n]) {
00263 if (u_inter[0] <= 0.0) {
00264 hlle_flag=1;
00265 break;
00266 }
00267 #ifdef ADIABATIC
00268 p_inter = u_inter[4] - 0.5*
00269 (SQR(u_inter[1])+SQR(u_inter[2])+SQR(u_inter[3]))/u_inter[0];
00270 #ifdef MHD
00271 p_inter -= 0.5*(SQR(u_inter[NWAVE-2])+SQR(u_inter[NWAVE-1])+SQR(Bxi));
00272 #endif
00273 if (p_inter < 0.0) {
00274 hlle_flag=2;
00275 break;
00276 }
00277 #endif
00278 }
00279 }
00280
00281 if (hlle_flag != 0) {
00282 flux_hlle(Ul,Ur,Wl,Wr,Bxi,pFlux);
00283 return;
00284 }
00285
00286 #endif
00287
00288
00289
00290
00291 pFl = (Real *)&(Fl);
00292 pFr = (Real *)&(Fr);
00293 pF = (Real *)(pFlux);
00294
00295 for (m=0; m<NWAVE; m++) {
00296 coeff[m] = 0.5*MAX(fabs(ev[m]),etah)*a[m];
00297 #ifdef CYLINDRICAL
00298 coeff2[m] = 0.5*SIGN(ev[m])*a[m];
00299 #endif
00300 }
00301 for (n=0; n<NWAVE; n++) {
00302 pF[n] = 0.5*(pFl[n] + pFr[n]);
00303 #ifdef CYLINDRICAL
00304 u_inter[n] = 0.5*(pUl[n] + pUr[n]);
00305 #endif
00306 for (m=0; m<NWAVE; m++) {
00307 pF[n] -= coeff[m]*rem[n][m];
00308 #ifdef CYLINDRICAL
00309 u_inter[n] -= coeff2[m]*rem[n][m];
00310 #endif
00311 }
00312 }
00313
00314
00315 #if (NSCALARS > 0)
00316 if (pFlux->d >= 0.0) {
00317 for (n=0; n<NSCALARS; n++) pFlux->s[n] = pFlux->d*Wl.r[n];
00318 } else {
00319 for (n=0; n<NSCALARS; n++) pFlux->s[n] = pFlux->d*Wr.r[n];
00320 }
00321 #endif
00322
00323 #if defined(CYLINDRICAL) && !defined(BAROTROPIC)
00324
00325 Emag = 0.0;
00326 #ifdef MHD
00327 Emag = 0.5*(SQR(u_inter[NWAVE-2])+SQR(u_inter[NWAVE-1])+SQR(Bxi));
00328 #endif
00329 Ekin = 0.5*(SQR(u_inter[1])+SQR(u_inter[2])+SQR(u_inter[3]))/u_inter[0];
00330 Eint = u_inter[4] - Emag - Ekin;
00331
00332 pFlux->Pflux = Eint*Gamma_1 + Emag;
00333 #endif
00334
00335 return;
00336 }
00337
00338
00339 #define HLLE_FUNCTION flux_hlle
00340 #include "hlle.c"
00341 #undef HLLE_FUNCTION
00342
00343 #endif