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 #include <math.h>
00026 #include <stdio.h>
00027 #include <stdlib.h>
00028 #include "../defs.h"
00029 #include "../athena.h"
00030 #include "../globals.h"
00031 #include "prototypes.h"
00032 #include "../prototypes.h"
00033
00034 #ifdef HLLD_FLUX
00035 #ifndef SPECIAL_RELATIVITY
00036
00037 #define SMALL_NUMBER 1e-8
00038
00039 #ifndef MHD
00040 #error : The HLLD flux only works for mhd.
00041 #endif
00042
00043 #ifndef ISOTHERMAL
00044
00045
00046
00047
00048
00049
00050
00051
00052
00053
00054
00055
00056
00057 void fluxes(const Cons1DS Ul, const Cons1DS Ur,
00058 const Prim1DS Wl, const Prim1DS Wr, const Real Bxi, Cons1DS *pFlux)
00059 {
00060 Cons1DS Ulst,Uldst,Urdst,Urst;
00061 Prim1DS Wlst,Wrst;
00062 Cons1DS Fl,Fr;
00063 Real spd[5];
00064
00065 Real sdl,sdr,sdml,sdmr;
00066 Real pbl,pbr;
00067 Real cfl,cfr,cfmax;
00068 Real gpl,gpr,gpbl,gpbr;
00069 Real sqrtdl,sqrtdr;
00070 Real invsumd;
00071 Real ptl,ptr,ptst;
00072 Real vbstl,vbstr;
00073 Real Bxsig;
00074 Real Bxsq = SQR(Bxi);
00075 Real tmp;
00076 #if (NSCALARS > 0)
00077 int n;
00078 #endif
00079
00080
00081
00082
00083
00084
00085
00086
00087
00088
00089
00090
00091
00092
00093
00094
00095 pbl = 0.5*(SQR(Bxi) + SQR(Wl.By) + SQR(Wl.Bz));
00096 pbr = 0.5*(SQR(Bxi) + SQR(Wr.By) + SQR(Wr.Bz));
00097 gpl = Gamma * Wl.P;
00098 gpr = Gamma * Wr.P;
00099 gpbl = gpl + 2.0*pbl;
00100 gpbr = gpr + 2.0*pbr;
00101
00102 cfl = sqrt((gpbl + sqrt(SQR(gpbl)-4*gpl*Bxsq))/(2.0*Wl.d));
00103 cfr = sqrt((gpbr + sqrt(SQR(gpbr)-4*gpr*Bxsq))/(2.0*Wr.d));
00104 cfmax = MAX(cfl,cfr);
00105
00106 if(Wl.Vx <= Wr.Vx) {
00107 spd[0] = Wl.Vx - cfmax;
00108 spd[4] = Wr.Vx + cfmax;
00109 }
00110 else {
00111 spd[0] = Wr.Vx - cfmax;
00112 spd[4] = Wl.Vx + cfmax;
00113 }
00114
00115
00116
00117
00118
00119
00120
00121
00122 ptl = Wl.P + pbl;
00123 ptr = Wr.P + pbr;
00124
00125 Fl.d = Ul.Mx;
00126 Fl.Mx = Ul.Mx*Wl.Vx + ptl - Bxsq;
00127 Fl.My = Ul.d*Wl.Vx*Wl.Vy - Bxi*Ul.By;
00128 Fl.Mz = Ul.d*Wl.Vx*Wl.Vz - Bxi*Ul.Bz;
00129 Fl.E = Wl.Vx*(Ul.E + ptl - Bxsq) - Bxi*(Wl.Vy*Ul.By + Wl.Vz*Ul.Bz);
00130 Fl.By = Ul.By*Wl.Vx - Bxi*Wl.Vy;
00131 Fl.Bz = Ul.Bz*Wl.Vx - Bxi*Wl.Vz;
00132
00133 Fr.d = Ur.Mx;
00134 Fr.Mx = Ur.Mx*Wr.Vx + ptr - Bxsq;
00135 Fr.My = Ur.d*Wr.Vx*Wr.Vy - Bxi*Ur.By;
00136 Fr.Mz = Ur.d*Wr.Vx*Wr.Vz - Bxi*Ur.Bz;
00137 Fr.E = Wr.Vx*(Ur.E + ptr - Bxsq) - Bxi*(Wr.Vy*Ur.By + Wr.Vz*Ur.Bz);
00138 Fr.By = Ur.By*Wr.Vx - Bxi*Wr.Vy;
00139 Fr.Bz = Ur.Bz*Wr.Vx - Bxi*Wr.Vz;
00140
00141 #if (NSCALARS > 0)
00142 for (n=0; n<NSCALARS; n++) {
00143 Fl.s[n] = Fl.d*Wl.r[n];
00144 Fr.s[n] = Fr.d*Wr.r[n];
00145 }
00146 #endif
00147
00148
00149
00150
00151
00152 if(spd[0] >= 0.0){
00153 *pFlux = Fl;
00154 #if defined(CYLINDRICAL) && !defined(BAROTROPIC)
00155 pFlux->Pflux = ptl;
00156 #endif
00157 return;
00158 }
00159
00160 if(spd[4] <= 0.0){
00161 *pFlux = Fr;
00162 #if defined(CYLINDRICAL) && !defined(BAROTROPIC)
00163 pFlux->Pflux = ptr;
00164 #endif
00165 return;
00166 }
00167
00168
00169
00170
00171
00172 sdl = spd[0] - Wl.Vx;
00173 sdr = spd[4] - Wr.Vx;
00174
00175
00176 spd[2] = (sdr*Wr.d*Wr.Vx - sdl*Wl.d*Wl.Vx - ptr + ptl) /
00177 (sdr*Wr.d-sdl*Wl.d);
00178
00179 sdml = spd[0] - spd[2];
00180 sdmr = spd[4] - spd[2];
00181
00182 Ulst.d = Ul.d * sdl/sdml;
00183 Urst.d = Ur.d * sdr/sdmr;
00184 sqrtdl = sqrt(Ulst.d);
00185 sqrtdr = sqrt(Urst.d);
00186
00187
00188 spd[1] = spd[2] - fabs(Bxi)/sqrtdl;
00189 spd[3] = spd[2] + fabs(Bxi)/sqrtdr;
00190
00191
00192
00193
00194
00195 ptst = ptl + Ul.d*sdl*(sdl-sdml);
00196
00197
00198
00199 Ulst.Mx = Ulst.d * spd[2];
00200
00201
00202
00203
00204
00205
00206
00207
00208
00209 if (fabs(Ul.d*sdl*sdml/Bxsq-1.0) < SMALL_NUMBER) {
00210
00211 Ulst.My = Ulst.d * Wl.Vy;
00212 Ulst.Mz = Ulst.d * Wl.Vz;
00213
00214 Ulst.By = Ul.By;
00215 Ulst.Bz = Ul.Bz;
00216 }
00217 else {
00218
00219 tmp = Bxi*(sdl-sdml)/(Ul.d*sdl*sdml-Bxsq);
00220 Ulst.My = Ulst.d * (Wl.Vy - Ul.By*tmp);
00221 Ulst.Mz = Ulst.d * (Wl.Vz - Ul.Bz*tmp);
00222
00223
00224
00225
00226
00227
00228
00229
00230
00231
00232
00233
00234 tmp = (Ul.d*SQR(sdl)-Bxsq)/(Ul.d*sdl*sdml - Bxsq);
00235 Ulst.By = Ul.By * tmp;
00236 Ulst.Bz = Ul.Bz * tmp;
00237 }
00238 vbstl = (Ulst.Mx*Bxi+Ulst.My*Ulst.By+Ulst.Mz*Ulst.Bz)/Ulst.d;
00239
00240 Ulst.E = (sdl*Ul.E - ptl*Wl.Vx + ptst*spd[2] +
00241 Bxi*(Wl.Vx*Bxi+Wl.Vy*Ul.By+Wl.Vz*Ul.Bz - vbstl))/sdml;
00242 Wlst = Cons1D_to_Prim1D(&Ulst,&Bxi);
00243
00244
00245
00246
00247 Urst.Mx = Urst.d * spd[2];
00248
00249
00250
00251
00252
00253
00254
00255
00256
00257 if (fabs(Ur.d*sdr*sdmr/Bxsq-1.0) < SMALL_NUMBER) {
00258
00259 Urst.My = Urst.d * Wr.Vy;
00260 Urst.Mz = Urst.d * Wr.Vz;
00261
00262 Urst.By = Ur.By;
00263 Urst.Bz = Ur.Bz;
00264 }
00265 else {
00266
00267 tmp = Bxi*(sdr-sdmr)/(Ur.d*sdr*sdmr-Bxsq);
00268 Urst.My = Urst.d * (Wr.Vy - Ur.By*tmp);
00269 Urst.Mz = Urst.d * (Wr.Vz - Ur.Bz*tmp);
00270
00271
00272
00273
00274
00275
00276
00277
00278
00279
00280
00281
00282
00283 tmp = (Ur.d*SQR(sdr)-Bxsq)/(Ur.d*sdr*sdmr - Bxsq);
00284 Urst.By = Ur.By * tmp;
00285 Urst.Bz = Ur.Bz * tmp;
00286 }
00287 vbstr = (Urst.Mx*Bxi+Urst.My*Urst.By+Urst.Mz*Urst.Bz)/Urst.d;
00288
00289 Urst.E = (sdr*Ur.E - ptr*Wr.Vx + ptst*spd[2] +
00290 Bxi*(Wr.Vx*Bxi+Wr.Vy*Ur.By+Wr.Vz*Ur.Bz - vbstr))/sdmr;
00291 Wrst = Cons1D_to_Prim1D(&Urst,&Bxi);
00292
00293
00294
00295
00296 if(0.5*Bxsq/MIN(pbl,pbr) < SQR(SMALL_NUMBER)) {
00297 Uldst = Ulst;
00298 Urdst = Urst;
00299 }
00300 else {
00301 invsumd = 1.0/(sqrtdl + sqrtdr);
00302 if(Bxi > 0) Bxsig = 1;
00303 else Bxsig = -1;
00304
00305 Uldst.d = Ulst.d;
00306 Urdst.d = Urst.d;
00307
00308 Uldst.Mx = Ulst.Mx;
00309 Urdst.Mx = Urst.Mx;
00310
00311
00312 tmp = invsumd*(sqrtdl*Wlst.Vy + sqrtdr*Wrst.Vy + Bxsig*(Urst.By-Ulst.By));
00313 Uldst.My = Uldst.d * tmp;
00314 Urdst.My = Urdst.d * tmp;
00315
00316
00317 tmp = invsumd*(sqrtdl*Wlst.Vz + sqrtdr*Wrst.Vz + Bxsig*(Urst.Bz-Ulst.Bz));
00318 Uldst.Mz = Uldst.d * tmp;
00319 Urdst.Mz = Urdst.d * tmp;
00320
00321
00322 tmp = invsumd*(sqrtdl*Urst.By + sqrtdr*Ulst.By +
00323 Bxsig*sqrtdl*sqrtdr*(Wrst.Vy-Wlst.Vy));
00324 Uldst.By = Urdst.By = tmp;
00325
00326
00327 tmp = invsumd*(sqrtdl*Urst.Bz + sqrtdr*Ulst.Bz +
00328 Bxsig*sqrtdl*sqrtdr*(Wrst.Vz-Wlst.Vz));
00329 Uldst.Bz = Urdst.Bz = tmp;
00330
00331
00332 tmp = spd[2]*Bxi + (Uldst.My*Uldst.By + Uldst.Mz*Uldst.Bz)/Uldst.d;
00333 Uldst.E = Ulst.E - sqrtdl*Bxsig*(vbstl - tmp);
00334 Urdst.E = Urst.E + sqrtdr*Bxsig*(vbstr - tmp);
00335 }
00336
00337
00338
00339
00340
00341 if(spd[1] >= 0) {
00342
00343 pFlux->d = Fl.d + spd[0]*(Ulst.d - Ul.d);
00344 pFlux->Mx = Fl.Mx + spd[0]*(Ulst.Mx - Ul.Mx);
00345 pFlux->My = Fl.My + spd[0]*(Ulst.My - Ul.My);
00346 pFlux->Mz = Fl.Mz + spd[0]*(Ulst.Mz - Ul.Mz);
00347 pFlux->E = Fl.E + spd[0]*(Ulst.E - Ul.E);
00348 pFlux->By = Fl.By + spd[0]*(Ulst.By - Ul.By);
00349 pFlux->Bz = Fl.Bz + spd[0]*(Ulst.Bz - Ul.Bz);
00350 }
00351 else if(spd[2] >= 0) {
00352
00353 tmp = spd[1] - spd[0];
00354 pFlux->d = Fl.d - spd[0]*Ul.d - tmp*Ulst.d + spd[1]*Uldst.d;
00355 pFlux->Mx = Fl.Mx - spd[0]*Ul.Mx - tmp*Ulst.Mx + spd[1]*Uldst.Mx;
00356 pFlux->My = Fl.My - spd[0]*Ul.My - tmp*Ulst.My + spd[1]*Uldst.My;
00357 pFlux->Mz = Fl.Mz - spd[0]*Ul.Mz - tmp*Ulst.Mz + spd[1]*Uldst.Mz;
00358 pFlux->E = Fl.E - spd[0]*Ul.E - tmp*Ulst.E + spd[1]*Uldst.E;
00359 pFlux->By = Fl.By - spd[0]*Ul.By - tmp*Ulst.By + spd[1]*Uldst.By;
00360 pFlux->Bz = Fl.Bz - spd[0]*Ul.Bz - tmp*Ulst.Bz + spd[1]*Uldst.Bz;
00361 }
00362 else if(spd[3] > 0) {
00363
00364 tmp = spd[3] - spd[4];
00365 pFlux->d = Fr.d - spd[4]*Ur.d - tmp*Urst.d + spd[3]*Urdst.d;
00366 pFlux->Mx = Fr.Mx - spd[4]*Ur.Mx - tmp*Urst.Mx + spd[3]*Urdst.Mx;
00367 pFlux->My = Fr.My - spd[4]*Ur.My - tmp*Urst.My + spd[3]*Urdst.My;
00368 pFlux->Mz = Fr.Mz - spd[4]*Ur.Mz - tmp*Urst.Mz + spd[3]*Urdst.Mz;
00369 pFlux->E = Fr.E - spd[4]*Ur.E - tmp*Urst.E + spd[3]*Urdst.E;
00370 pFlux->By = Fr.By - spd[4]*Ur.By - tmp*Urst.By + spd[3]*Urdst.By;
00371 pFlux->Bz = Fr.Bz - spd[4]*Ur.Bz - tmp*Urst.Bz + spd[3]*Urdst.Bz;
00372 }
00373 else {
00374
00375 pFlux->d = Fr.d + spd[4]*(Urst.d - Ur.d);
00376 pFlux->Mx = Fr.Mx + spd[4]*(Urst.Mx - Ur.Mx);
00377 pFlux->My = Fr.My + spd[4]*(Urst.My - Ur.My);
00378 pFlux->Mz = Fr.Mz + spd[4]*(Urst.Mz - Ur.Mz);
00379 pFlux->E = Fr.E + spd[4]*(Urst.E - Ur.E);
00380 pFlux->By = Fr.By + spd[4]*(Urst.By - Ur.By);
00381 pFlux->Bz = Fr.Bz + spd[4]*(Urst.Bz - Ur.Bz);
00382 }
00383
00384
00385 #if (NSCALARS > 0)
00386 if (pFlux->d >= 0.0) {
00387 for (n=0; n<NSCALARS; n++) pFlux->s[n] = pFlux->d*Wl.r[n];
00388 } else {
00389 for (n=0; n<NSCALARS; n++) pFlux->s[n] = pFlux->d*Wr.r[n];
00390 }
00391 #endif
00392
00393 #if defined(CYLINDRICAL) && !defined(BAROTROPIC)
00394 pFlux->Pflux = ptst;
00395 #endif
00396 return;
00397 }
00398
00399 #else
00400
00401
00402
00403
00404
00405
00406
00407
00408
00409
00410
00411
00412
00413 void fluxes(const Cons1DS Ul, const Cons1DS Ur,
00414 const Prim1DS Wl, const Prim1DS Wr, const Real Bxi, Cons1DS *pFlux)
00415 {
00416 Cons1DS Ulst,Ucst,Urst;
00417 Cons1DS Fl,Fr;
00418 Real spd[5],idspd;
00419 Real pbl,pbr;
00420 Real cfl,cfr;
00421 Real gpl,gpr,gpbl,gpbr;
00422 Real mxhll,ustar,dhll,sqrtdhll;
00423 Real fdhll,fmxhll;
00424 Real ptl,ptr;
00425 Real Bxsq = SQR(Bxi);
00426 Real tmp,mfact,bfact,X;
00427 int n;
00428
00429
00430
00431
00432
00433
00434
00435
00436
00437
00438
00439
00440
00441
00442
00443
00444
00445
00446
00447 pbl = 0.5*(SQR(Bxi) + SQR(Wl.By) + SQR(Wl.Bz));
00448 pbr = 0.5*(SQR(Bxi) + SQR(Wr.By) + SQR(Wr.Bz));
00449 gpl = Wl.d*Iso_csound2;
00450 gpr = Wr.d*Iso_csound2;
00451 gpbl = gpl + 2.0*pbl;
00452 gpbr = gpr + 2.0*pbr;
00453
00454 cfl = sqrt((gpbl + sqrt(SQR(gpbl)-4*gpl*Bxsq))/(2.0*Wl.d));
00455 cfr = sqrt((gpbr + sqrt(SQR(gpbr)-4*gpr*Bxsq))/(2.0*Wr.d));
00456
00457 spd[0] = MIN(Wl.Vx-cfl,Wr.Vx-cfr);
00458 spd[4] = MAX(Wl.Vx+cfl,Wr.Vx+cfr);
00459
00460
00461
00462
00463
00464
00465
00466
00467
00468
00469
00470
00471
00472
00473
00474
00475
00476
00477
00478
00479
00480
00481
00482
00483
00484
00485
00486
00487
00488 ptl = gpl + pbl;
00489 ptr = gpr + pbr;
00490
00491 Fl.d = Ul.Mx;
00492 Fl.Mx = Ul.Mx*Wl.Vx + ptl - Bxsq;
00493 Fl.My = Ul.d*Wl.Vx*Wl.Vy - Bxi*Ul.By;
00494 Fl.Mz = Ul.d*Wl.Vx*Wl.Vz - Bxi*Ul.Bz;
00495 Fl.By = Ul.By*Wl.Vx - Bxi*Wl.Vy;
00496 Fl.Bz = Ul.Bz*Wl.Vx - Bxi*Wl.Vz;
00497
00498 Fr.d = Ur.Mx;
00499 Fr.Mx = Ur.Mx*Wr.Vx + ptr - Bxsq;
00500 Fr.My = Ur.d*Wr.Vx*Wr.Vy - Bxi*Ur.By;
00501 Fr.Mz = Ur.d*Wr.Vx*Wr.Vz - Bxi*Ur.Bz;
00502 Fr.By = Ur.By*Wr.Vx - Bxi*Wr.Vy;
00503 Fr.Bz = Ur.Bz*Wr.Vx - Bxi*Wr.Vz;
00504
00505 #if (NSCALARS > 0)
00506 for (n=0; n<NSCALARS; n++) {
00507 Fl.s[n] = Fl.d*Wl.r[n];
00508 Fr.s[n] = Fr.d*Wr.r[n];
00509 }
00510 #endif
00511
00512
00513
00514
00515
00516
00517 if(spd[0] >= 0.0){
00518 *pFlux = Fl;
00519 return;
00520 }
00521
00522
00523 if(spd[4] <= 0.0){
00524 *pFlux = Fr;
00525 return;
00526 }
00527
00528
00529
00530
00531
00532
00533 idspd = 1.0/(spd[4]-spd[0]);
00534
00535
00536
00537 dhll = (spd[4]*Ur.d-spd[0]*Ul.d-Fr.d+Fl.d)*idspd;
00538 sqrtdhll = sqrt(dhll);
00539
00540
00541 fdhll = (spd[4]*Fl.d-spd[0]*Fr.d+spd[4]*spd[0]*(Ur.d-Ul.d))*idspd;
00542 fmxhll = (spd[4]*Fl.Mx-spd[0]*Fr.Mx+spd[4]*spd[0]*(Ur.Mx-Ul.Mx))*idspd;
00543
00544
00545 ustar = fdhll/dhll;
00546
00547
00548
00549 mxhll = (spd[4]*Ur.Mx-spd[0]*Ul.Mx-Fr.Mx+Fl.Mx)*idspd;
00550
00551
00552 spd[1] = ustar - fabs(Bxi)/sqrtdhll;
00553 spd[3] = ustar + fabs(Bxi)/sqrtdhll;
00554
00555
00556
00557
00558
00559
00560
00561 Ulst.d = dhll;
00562
00563 Ulst.Mx = mxhll;
00564
00565 tmp = (spd[0]-spd[1])*(spd[0]-spd[3]);
00566
00567 if ((fabs(spd[0]/spd[1]-1.0) < SMALL_NUMBER)
00568 || (fabs(spd[0]/spd[3]-1.0) < SMALL_NUMBER)) {
00569
00570 Ulst.My = Ul.My;
00571 Ulst.Mz = Ul.Mz;
00572 Ulst.By = Ul.By;
00573 Ulst.Bz = Ul.Bz;
00574 } else {
00575 mfact = Bxi*(ustar-Wl.Vx)/tmp;
00576 bfact = (Ul.d*SQR(spd[0]-Wl.Vx)-Bxsq)/(dhll*tmp);
00577
00578
00579 Ulst.My = dhll*Wl.Vy-Ul.By*mfact;
00580
00581 Ulst.Mz = dhll*Wl.Vz-Ul.Bz*mfact;
00582
00583 Ulst.By = Ul.By*bfact;
00584
00585 Ulst.Bz = Ul.Bz*bfact;
00586 }
00587
00588
00589
00590 Urst.d = dhll;
00591
00592 Urst.Mx = mxhll;
00593
00594 tmp = (spd[4]-spd[1])*(spd[4]-spd[3]);
00595
00596 if ((fabs(spd[4]/spd[1]-1.0) < SMALL_NUMBER)
00597 || (fabs(spd[4]/spd[3]-1.0) < SMALL_NUMBER)) {
00598
00599 Urst.My = Ur.My;
00600 Urst.Mz = Ur.Mz;
00601 Urst.By = Ur.By;
00602 Urst.Bz = Ur.Bz;
00603 } else {
00604 mfact = Bxi*(ustar-Wr.Vx)/tmp;
00605 bfact = (Ur.d*SQR(spd[4]-Wr.Vx)-Bxsq)/(dhll*tmp);
00606
00607
00608 Urst.My = dhll*Wr.Vy-Ur.By*mfact;
00609
00610 Urst.Mz = dhll*Wr.Vz-Ur.Bz*mfact;
00611
00612 Urst.By = Ur.By*bfact;
00613
00614 Urst.Bz = Ur.Bz*bfact;
00615 }
00616
00617
00618
00619 X = sqrtdhll*SIGN(Bxi);
00620
00621 Ucst.d = dhll;
00622
00623 Ucst.Mx = mxhll;
00624
00625 Ucst.My = 0.5*(Ulst.My+Urst.My+X*(Urst.By-Ulst.By));
00626
00627 Ucst.Mz = 0.5*(Ulst.Mz+Urst.Mz+X*(Urst.Bz-Ulst.Bz));
00628
00629 Ucst.By = 0.5*(Ulst.By+Urst.By+(Urst.My-Ulst.My)/X);
00630
00631 Ucst.Bz = 0.5*(Ulst.Bz+Urst.Bz+(Urst.Mz-Ulst.Mz)/X);
00632
00633
00634
00635
00636
00637 if(spd[1] >= 0) {
00638
00639 pFlux->d = Fl.d + spd[0]*(Ulst.d - Ul.d);
00640 pFlux->Mx = Fl.Mx + spd[0]*(Ulst.Mx - Ul.Mx);
00641 pFlux->My = Fl.My + spd[0]*(Ulst.My - Ul.My);
00642 pFlux->Mz = Fl.Mz + spd[0]*(Ulst.Mz - Ul.Mz);
00643 pFlux->By = Fl.By + spd[0]*(Ulst.By - Ul.By);
00644 pFlux->Bz = Fl.Bz + spd[0]*(Ulst.Bz - Ul.Bz);
00645 }
00646 else if (spd[3] <= 0) {
00647
00648 pFlux->d = Fr.d + spd[4]*(Urst.d - Ur.d);
00649 pFlux->Mx = Fr.Mx + spd[4]*(Urst.Mx - Ur.Mx);
00650 pFlux->My = Fr.My + spd[4]*(Urst.My - Ur.My);
00651 pFlux->Mz = Fr.Mz + spd[4]*(Urst.Mz - Ur.Mz);
00652 pFlux->By = Fr.By + spd[4]*(Urst.By - Ur.By);
00653 pFlux->Bz = Fr.Bz + spd[4]*(Urst.Bz - Ur.Bz);
00654 }
00655 else {
00656
00657 pFlux->d = dhll*ustar;
00658 pFlux->Mx = fmxhll;
00659 pFlux->My = Ucst.My*ustar - Bxi*Ucst.By;
00660 pFlux->Mz = Ucst.Mz*ustar - Bxi*Ucst.Bz;
00661 pFlux->By = Ucst.By*ustar - Bxi*Ucst.My/Ucst.d;
00662 pFlux->Bz = Ucst.Bz*ustar - Bxi*Ucst.Mz/Ucst.d;
00663 }
00664
00665
00666 #if (NSCALARS > 0)
00667 if (pFlux->d >= 0.0) {
00668 for (n=0; n<NSCALARS; n++) pFlux->s[n] = pFlux->d*Wl.r[n];
00669 } else {
00670 for (n=0; n<NSCALARS; n++) pFlux->s[n] = pFlux->d*Wr.r[n];
00671 }
00672 #endif
00673
00674 return;
00675 }
00676
00677 #endif
00678 #endif
00679 #endif