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 #include <math.h>
00031 #include <stdio.h>
00032 #include <stdlib.h>
00033 #include "../defs.h"
00034 #include "../athena.h"
00035 #include "../globals.h"
00036 #include "prototypes.h"
00037 #include "../prototypes.h"
00038
00039 #ifdef SECOND_ORDER_CHAR
00040 #ifdef SPECIAL_RELATIVITY
00041 #error : PLM reconstruction (order=2) cannot be used for special relativity.
00042 #endif
00043
00044 static Real **pW=NULL;
00045
00046
00047
00048
00049
00050
00051
00052
00053
00054
00055
00056
00057
00058
00059
00060
00061
00062 void lr_states(const GridS *pG, const Prim1DS W[], const Real Bxc[],
00063 const Real dt, const Real dx, const int il, const int iu,
00064 Prim1DS Wl[], Prim1DS Wr[], const int dir)
00065 {
00066 int i,n,m;
00067 Real lim_slope1,lim_slope2,qa,qx;
00068 Real ev[NWAVE],rem[NWAVE][NWAVE],lem[NWAVE][NWAVE];
00069 Real dWc[NWAVE+NSCALARS],dWl[NWAVE+NSCALARS];
00070 Real dWr[NWAVE+NSCALARS],dWg[NWAVE+NSCALARS];
00071 Real dac[NWAVE+NSCALARS],dal[NWAVE+NSCALARS];
00072 Real dar[NWAVE+NSCALARS],dag[NWAVE+NSCALARS],da[NWAVE+NSCALARS];
00073 Real Wlv[NWAVE+NSCALARS],Wrv[NWAVE+NSCALARS];
00074 Real dW[NWAVE+NSCALARS],dWm[NWAVE+NSCALARS];
00075 Real *pWl, *pWr;
00076 Real qx1,qx2,C;
00077
00078
00079 Real zl,zr,zc,gamma_curv,opg,omg,beta,betai;
00080 int hllallwave_flag = 0;
00081 const Real dtodx = dt/dx;
00082 #ifdef CYLINDRICAL
00083 const Real *r=pG->r, *ri=pG->ri;
00084 #endif
00085
00086
00087 for (n=0; n<NWAVE; n++) {
00088 for (m=0; m<NWAVE; m++) {
00089 rem[n][m] = 0.0;
00090 lem[n][m] = 0.0;
00091 }
00092 }
00093 for (i=il-2; i<=iu+2; i++) pW[i] = (Real*)&(W[i]);
00094
00095
00096 for (i=il-1; i<=iu+1; i++) {
00097
00098
00099
00100
00101 #ifdef HYDRO
00102 #ifdef ISOTHERMAL
00103 esys_prim_iso_hyd(W[i].d,W[i].Vx, ev,rem,lem);
00104 #else
00105 esys_prim_adb_hyd(W[i].d,W[i].Vx,Gamma*W[i].P,ev,rem,lem);
00106 #endif
00107 #endif
00108
00109 #ifdef MHD
00110 #ifdef ISOTHERMAL
00111 esys_prim_iso_mhd(
00112 W[i].d,W[i].Vx, Bxc[i],W[i].By,W[i].Bz,ev,rem,lem);
00113 #else
00114 esys_prim_adb_mhd(
00115 W[i].d,W[i].Vx,Gamma*W[i].P,Bxc[i],W[i].By,W[i].Bz,ev,rem,lem);
00116 #endif
00117 #endif
00118
00119
00120
00121
00122
00123 #ifdef CYLINDRICAL
00124 if (dir==1) {
00125
00126 zc = 1.0/(1.0 - SQR(dx)/(12.0*r[i+1]*r[i-1]));
00127 zl = 1.0/(1.0 - SQR(dx)/(12.0*r[i ]*r[i-1]));
00128 zr = 1.0/(1.0 - SQR(dx)/(12.0*r[i+1]*r[i ]));
00129 }
00130 #endif
00131 for (n=0; n<(NWAVE+NSCALARS); n++) {
00132 dWc[n] = pW[i+1][n] - pW[i-1][n];
00133 dWl[n] = pW[i][n] - pW[i-1][n];
00134 dWr[n] = pW[i+1][n] - pW[i][n];
00135 #ifdef CYLINDRICAL
00136 if (dir==1) {
00137 dWc[n] *= zc;
00138 dWl[n] *= zl;
00139 dWr[n] *= zr;
00140 }
00141 #endif
00142 if (dWl[n]*dWr[n] > 0.0) {
00143 dWg[n] = 2.0*dWl[n]*dWr[n]/(dWl[n]+dWr[n]);
00144 } else {
00145 dWg[n] = 0.0;
00146 }
00147 }
00148
00149
00150
00151
00152 for (n=0; n<NWAVE; n++) {
00153 dac[n] = lem[n][0]*dWc[0];
00154 dal[n] = lem[n][0]*dWl[0];
00155 dar[n] = lem[n][0]*dWr[0];
00156 dag[n] = lem[n][0]*dWg[0];
00157 for (m=1; m<NWAVE; m++) {
00158 dac[n] += lem[n][m]*dWc[m];
00159 dal[n] += lem[n][m]*dWl[m];
00160 dar[n] += lem[n][m]*dWr[m];
00161 dag[n] += lem[n][m]*dWg[m];
00162 }
00163 }
00164
00165
00166
00167
00168 #if (NSCALARS > 0)
00169 for (n=NWAVE; n<(NWAVE+NSCALARS); n++) {
00170 dac[n] = dWc[n];
00171 dal[n] = dWl[n];
00172 dar[n] = dWr[n];
00173 dag[n] = dWg[n];
00174 }
00175 #endif
00176
00177
00178
00179
00180 for (n=0; n<(NWAVE+NSCALARS); n++) {
00181 da[n] = 0.0;
00182 if (dal[n]*dar[n] > 0.0) {
00183 lim_slope1 = MIN( fabs(dal[n]),fabs(dar[n]));
00184 lim_slope2 = MIN(0.5*fabs(dac[n]),fabs(dag[n]));
00185 da[n] = SIGN(dac[n])*MIN(2.0*lim_slope1,lim_slope2);
00186 }
00187 }
00188
00189
00190
00191
00192 for (n=0; n<NWAVE; n++) {
00193 dWm[n] = da[0]*rem[n][0];
00194 for (m=1; m<NWAVE; m++) {
00195 dWm[n] += da[m]*rem[n][m];
00196 }
00197 }
00198
00199 #if (NSCALARS > 0)
00200 for (n=NWAVE; n<(NWAVE+NSCALARS); n++) {
00201 dWm[n] = da[n];
00202 }
00203 #endif
00204
00205
00206
00207
00208
00209
00210
00211
00212
00213 gamma_curv = 0.0;
00214 #ifdef CYLINDRICAL
00215 if (dir==1) gamma_curv = dx/(6.0*r[i]);
00216 #endif
00217 opg = 1.0 + gamma_curv;
00218 omg = 1.0 - gamma_curv;
00219 beta = omg/opg;
00220 betai = opg/omg;
00221
00222 for (n=0; n<(NWAVE+NSCALARS); n++) {
00223 Wlv[n] = pW[i][n] - 0.5*dWm[n]*opg;
00224 Wrv[n] = pW[i][n] + 0.5*dWm[n]*omg;
00225 }
00226
00227 for (n=0; n<(NWAVE+NSCALARS); n++) {
00228 C = Wrv[n] + beta*Wlv[n];
00229 Wlv[n] = MAX(MIN(pW[i][n],pW[i-1][n]),Wlv[n]);
00230 Wlv[n] = MIN(MAX(pW[i][n],pW[i-1][n]),Wlv[n]);
00231 Wrv[n] = C - beta*Wlv[n];
00232
00233 Wrv[n] = MAX(MIN(pW[i][n],pW[i+1][n]),Wrv[n]);
00234 Wrv[n] = MIN(MAX(pW[i][n],pW[i+1][n]),Wrv[n]);
00235 Wlv[n] = (C - Wrv[n])*betai;
00236 }
00237
00238 for (n=0; n<(NWAVE+NSCALARS); n++) {
00239 dW[n] = Wrv[n] - Wlv[n];
00240 }
00241
00242
00243
00244
00245
00246
00247 pWl = (Real *) &(Wl[i+1]);
00248 pWr = (Real *) &(Wr[i]);
00249
00250 #ifndef CTU_INTEGRATOR
00251
00252 for (n=0; n<(NWAVE+NSCALARS); n++) {
00253 pWl[n] = Wrv[n];
00254 pWr[n] = Wlv[n];
00255 }
00256
00257 #elif defined(HLLE_FLUX) || defined(HLLC_FLUX) || defined(HLLD_FLUX)
00258
00259 for (n=0; n<NWAVE; n++) {
00260 pWl[n] = Wrv[n];
00261 pWr[n] = Wlv[n];
00262 }
00263
00264 #ifdef HLL_ALL_WAVE
00265 hllallwave_flag = 1;
00266 #endif
00267
00268 for (n=0; n<NWAVE; n++) {
00269 qa = 0.0;
00270 if (hllallwave_flag || ev[n] > 0.) {
00271 qx = 0.5*dtodx*ev[n];
00272 #ifdef CYLINDRICAL
00273 if (dir==1)
00274 qx *= 1.0 - dx*qx/(3.0*(ri[i+1]-dx*qx));
00275 #endif
00276 for (m=0; m<NWAVE; m++) {
00277 qa += lem[n][m]*qx*dW[m];
00278 }
00279 for (m=0; m<NWAVE; m++) pWl[m] -= qa*rem[m][n];
00280 }
00281
00282 qa = 0.0;
00283 if (hllallwave_flag || ev[n] < 0.) {
00284 qx = 0.5*dtodx*ev[n];
00285 #ifdef CYLINDRICAL
00286 if (dir==1)
00287 qx *= 1.0 - dx*qx/(3.0*(ri[i]-dx*qx));
00288 #endif
00289 for (m=0; m<NWAVE; m++) {
00290 qa += lem[n][m]*qx*dW[m];
00291 }
00292 for (m=0; m<NWAVE; m++) pWr[m] -= qa*rem[m][n];
00293 }
00294 }
00295
00296 #else
00297 qx = 0.5*MAX(ev[NWAVE-1],0.0)*dtodx;
00298 #ifdef CYLINDRICAL
00299 if (dir==1)
00300 qx *= 1.0 - dx*qx/(3.0*(ri[i+1]-dx*qx));
00301 #endif
00302 for (n=0; n<(NWAVE+NSCALARS); n++) {
00303 pWl[n] = Wrv[n] - qx*dW[n];
00304 }
00305
00306 qx = -0.5*MIN(ev[0],0.0)*dtodx;
00307 #ifdef CYLINDRICAL
00308 if (dir==1)
00309 qx *= 1.0 + dx*qx/(3.0*(ri[i]+dx*qx));
00310 #endif
00311 for (n=0; n<(NWAVE+NSCALARS); n++) {
00312 pWr[n] = Wlv[n] + qx*dW[n];
00313 }
00314
00315
00316
00317
00318
00319
00320
00321
00322 for (n=0; n<NWAVE; n++) {
00323 if (ev[n] >= 0.0) {
00324 qa = 0.0;
00325 qx1 = 0.5*dtodx*ev[NWAVE-1];
00326 qx2 = 0.5*dtodx*ev[n];
00327 #ifdef CYLINDRICAL
00328 if (dir==1) {
00329 qx1 *= 1.0 - dx*qx1/(3.0*(ri[i+1]-dx*qx1));
00330 qx2 *= 1.0 - dx*qx2/(3.0*(ri[i+1]-dx*qx2));
00331 }
00332 #endif
00333 qx = qx1 - qx2;
00334 for (m=0; m<NWAVE; m++) {
00335 qa += lem[n][m]*qx*dW[m];
00336 }
00337 for (m=0; m<NWAVE; m++) pWl[m] += qa*rem[m][n];
00338 }
00339 }
00340
00341 for (n=0; n<NWAVE; n++) {
00342 if (ev[n] <= 0.0) {
00343 qa = 0.0;
00344 qx1 = -0.5*dtodx*ev[0];
00345 qx2 = -0.5*dtodx*ev[n];
00346 #ifdef CYLINDRICAL
00347 if (dir==1) {
00348 qx1 *= 1.0 + dx*qx1/(3.0*(ri[i]+dx*qx1));
00349 qx2 *= 1.0 + dx*qx2/(3.0*(ri[i]+dx*qx2));
00350 }
00351 #endif
00352 qx = -qx1 + qx2;
00353 for (m=0; m<NWAVE; m++) {
00354 qa += lem[n][m]*qx*dW[m];
00355 }
00356 for (m=0; m<NWAVE; m++) pWr[m] += qa*rem[m][n];
00357 }
00358 }
00359
00360
00361 for (n=NWAVE; n<(NWAVE+NSCALARS); n++) {
00362 if (W[i].Vx > 0.) {
00363 pWl[n] += 0.5*dtodx*(ev[NWAVE-1]-W[i].Vx)*dW[n];
00364 } else if (W[i].Vx < 0.) {
00365 pWr[n] += 0.5*dtodx*(ev[0]-W[i].Vx)*dW[n];
00366 }
00367 }
00368
00369 #endif
00370
00371 }
00372
00373 return;
00374 }
00375
00376
00377
00378
00379
00380
00381 void lr_states_init(MeshS *pM)
00382 {
00383 int nmax,size1=0,size2=0,size3=0,nl,nd;
00384
00385
00386 for (nl=0; nl<(pM->NLevels); nl++){
00387 for (nd=0; nd<(pM->DomainsPerLevel[nl]); nd++){
00388 if (pM->Domain[nl][nd].Grid != NULL) {
00389 if (pM->Domain[nl][nd].Grid->Nx[0] > size1){
00390 size1 = pM->Domain[nl][nd].Grid->Nx[0];
00391 }
00392 if (pM->Domain[nl][nd].Grid->Nx[1] > size2){
00393 size2 = pM->Domain[nl][nd].Grid->Nx[1];
00394 }
00395 if (pM->Domain[nl][nd].Grid->Nx[2] > size3){
00396 size3 = pM->Domain[nl][nd].Grid->Nx[2];
00397 }
00398 }
00399 }
00400 }
00401
00402 size1 = size1 + 2*nghost;
00403 size2 = size2 + 2*nghost;
00404 size3 = size3 + 2*nghost;
00405 nmax = MAX((MAX(size1,size2)),size3);
00406
00407 if ((pW = (Real**)malloc(nmax*sizeof(Real*))) == NULL) goto on_error;
00408
00409 return;
00410 on_error:
00411 lr_states_destruct();
00412 ath_error("[lr_states_init]: malloc returned a NULL pointer\n");
00413 }
00414
00415
00416
00417
00418
00419 void lr_states_destruct(void)
00420 {
00421 if (pW != NULL) free(pW);
00422 return;
00423 }
00424
00425 #endif