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