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