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