HARM
harm and utilities
 All Data Structures Files Functions Variables Typedefs Macros Pages
fluxcompute.c
Go to the documentation of this file.
1 
12 #include "decs.h"
13 
14 
15 
16 static void diag_fluxdump_1(int dir, int i, int j, int k, FTYPE *p_l, FTYPE *p_r, FTYPE *F_l, FTYPE *F_r);
17 
18 
25 int flux_compute_general(int i, int j, int k, int dir, struct of_geom *ptrgeom, FTYPE CUf, FTYPE *p_c, FTYPE *p_l, FTYPE *p_r, FTYPE *F, FTYPE *ctopallptr)
26 {
27  int flux_compute(int i, int j, int k, int dir, struct of_geom *geom, FTYPE *cminmax_l, FTYPE *cminmax_r, FTYPE *cminmax, FTYPE ctopmhd, FTYPE *cminmaxrad_l, FTYPE *cminmaxrad_r, FTYPE *cminmaxrad, FTYPE ctoprad, FTYPE CUf, FTYPE *p_l, FTYPE *p_r, FTYPE *U_l, FTYPE *U_r, FTYPE *F_l, FTYPE *F_r, FTYPE *F);
28  int p2SFUevolve(int dir, int isleftright, FTYPE *p, struct of_geom *geom, struct of_state **ptrstate, FTYPE *F, FTYPE *U);
29  FTYPE cminmax_l[NUMCS], cminmax_r[NUMCS], cminmax[NUMCS];
30  FTYPE cminmaxrad_l[NUMCS], cminmaxrad_r[NUMCS], cminmaxrad[NUMCS];
31  FTYPE cminmaxrad2_l[NUMCS], cminmaxrad2_r[NUMCS], cminmaxrad2[NUMCS];
32  FTYPE ctopmhd,ctoprad,ctoprad2;
33  struct of_state state_c, state_l, state_r;
34  struct of_state *ptrstate_c, *ptrstate_l, *ptrstate_r;
35  FTYPE F_c[NPR], F_l[NPR], F_r[NPR];
36  FTYPE U_c[NPR], U_l[NPR], U_r[NPR];
37  int pl,pliter;
38  int jj,kk;
39 
40 
41  // default
42  ptrstate_c = &state_c;
43  ptrstate_l = &state_l;
44  ptrstate_r = &state_r;
45 
46 
47 
49  //
50  // setup flux calculation based upon interpolated primitives
51  //
53  MYFUN(p2SFUevolve(dir, ISLEFT, p_l, ptrgeom, &ptrstate_l, F_l, U_l),"step_ch.c:fluxcalc()", "p2SFUevolve()", 1);
54  MYFUN(p2SFUevolve(dir, ISRIGHT, p_r, ptrgeom, &ptrstate_r, F_r, U_r),"step_ch.c:fluxcalc()", "p2SFUevolve()", 2);
55 
56 
57 #if(FLUXDUMP==2)
58  FTYPE FPAKE[NPR], FPAKE_l[NPR], FPAKE_r[NPR], UPAKE_l[NPR], UPAKE_r[NPR];
59  FTYPE FEN[NPR], FEN_l[NPR], FEN_r[NPR], UEN_l[NPR], UEN_r[NPR];
60  FTYPE FEM[NPR], FEM_l[NPR], FEM_r[NPR], UEM_l[NPR], UEM_r[NPR];
61  FTYPE FRAD[NPR], FRAD_l[NPR], FRAD_r[NPR], URAD_l[NPR], URAD_r[NPR];
62  int p2SFUevolveall(int dir, int isleftright, FTYPE *p, struct of_geom *ptrgeom, struct of_state **ptrstate, FTYPE *FPAKE, FTYPE *FEN, FTYPE *FEM, FTYPE *FRAD, FTYPE *UPAKE, FTYPE *UEN, FTYPE *UEM, FTYPE *URAD);
63  MYFUN(p2SFUevolveall(dir, ISLEFT, p_l, ptrgeom, &ptrstate_l, FPAKE_l, FEN_l, FEM_l, FRAD_l, UPAKE_l, UEN_l, UEM_l, URAD_l),"step_ch.c:fluxcalc()", "p2SFUevolve()", 1);
64  MYFUN(p2SFUevolveall(dir, ISRIGHT, p_r, ptrgeom, &ptrstate_r, FPAKE_r, FEN_r, FEM_r, FRAD_r, UPAKE_r, UEN_r, UEM_r, URAD_r),"step_ch.c:fluxcalc()", "p2SFUevolve()", 2);
65 #endif
66 
67 
68 
69  // usually "always" need cminmax_l cminmax_r and always need ctop
70  get_wavespeeds(dir, ptrgeom, p_l, p_r, U_l, U_r, F_l, F_r, ptrstate_l, ptrstate_r, cminmax_l, cminmax_r, cminmax, &ctopmhd, cminmaxrad_l, cminmaxrad_r, cminmaxrad, &ctoprad, cminmaxrad2_l, cminmaxrad2_r, cminmaxrad2, &ctoprad2);
71 
73  // choose ctopall that will go into choosing timestep
74  if(EOMRADTYPE!=EOMRADNONE){
75  // get maximum over both mhd and radiation
76  *ctopallptr=MAX(ctopmhd,ctoprad2); // choose ctoprad2 for setting timestep
77 
78  if(FORCESOLVEL){
79  *ctopallptr=1.0/sqrt(ptrgeom->gcov[GIND(dir,dir)]);
80  }
81 
82  }
83  else *ctopallptr=ctopmhd;
84 
85 
86 
87  // dualfprintf(fail_file,"ctopmhdprim=%g ctopradprim=%g\n",ctopmhd*sqrt(ptrgeom->gcov[GIND(dir,dir)]),ctoprad*sqrt(ptrgeom->gcov[GIND(dir,dir)]));
88 
89 
90 
91  // cminmaxrad and ctoprad are used for fluxes (not "2" version that's only for setting timestep)
92  MYFUN(flux_compute(i, j, k, dir, ptrgeom, cminmax_l,cminmax_r, cminmax, ctopmhd, cminmaxrad_l,cminmaxrad_r, cminmaxrad, ctoprad, CUf, p_l, p_r, U_l, U_r, F_l, F_r, F),"step_ch.c:fluxcalc()", "flux_compute", 1);
93 
94 
95 #if(FLUXDUMP==2)
96  if(dir==1){ // only for radial fluxes
97  // no accounting for RK, so only correct for RK2 midpoint method right now
98  PLOOP(pliter,pl) GLOBALMACP0A1(fluxdump,i,j,k,0*NPR + pl)=F[pl];
99 
100 
101  MYFUN(flux_compute(i, j, k, dir, ptrgeom, cminmax_l,cminmax_r, cminmax, ctopmhd, cminmaxrad_l,cminmaxrad_r, cminmaxrad, ctoprad, CUf, p_l, p_r, UPAKE_l, UPAKE_r, FPAKE_l, FPAKE_r, FPAKE),"step_ch.c:fluxcalc()", "flux_compute", 1);
102  PLOOP(pliter,pl) GLOBALMACP0A1(fluxdump,i,j,k,1*NPR + pl)=FPAKE[pl];
103 
104  MYFUN(flux_compute(i, j, k, dir, ptrgeom, cminmax_l,cminmax_r, cminmax, ctopmhd, cminmaxrad_l,cminmaxrad_r, cminmaxrad, ctoprad, CUf, p_l, p_r, UEN_l, UEN_r, FEN_l, FEN_r, FEN),"step_ch.c:fluxcalc()", "flux_compute", 1);
105  PLOOP(pliter,pl) GLOBALMACP0A1(fluxdump,i,j,k,2*NPR + pl)=FEN[pl];
106 
107  MYFUN(flux_compute(i, j, k, dir, ptrgeom, cminmax_l,cminmax_r, cminmax, ctopmhd, cminmaxrad_l,cminmaxrad_r, cminmaxrad, ctoprad, CUf, p_l, p_r, UEM_l, UEM_r, FEM_l, FEM_r, FEM),"step_ch.c:fluxcalc()", "flux_compute", 1);
108  PLOOP(pliter,pl) GLOBALMACP0A1(fluxdump,i,j,k,3*NPR + pl)=FEM[pl];
109 
110  if(EOMRADTYPE!=EOMRADNONE){
111  MYFUN(flux_compute(i, j, k, dir, ptrgeom, cminmax_l,cminmax_r, cminmax, ctopmhd, cminmaxrad_l,cminmaxrad_r, cminmaxrad, ctoprad, CUf, p_l, p_r, URAD_l, URAD_r, FRAD_l, FRAD_r, FRAD),"step_ch.c:fluxcalc()", "flux_compute", 1);
112  PLOOP(pliter,pl) GLOBALMACP0A1(fluxdump,i,j,k,4*NPR + pl)=FRAD[pl];
113  }
114  }
115 #endif
116 
117 #if(FLUXDUMP==1)
118  diag_fluxdump_1(dir,i,j,k,p_l, p_r, F_l, F_r);
119 #endif
120 
121  return(0);
122 }
123 
124 
125 
129 int flux_compute_splitmaem(int i, int j, int k, int dir, struct of_geom *ptrgeom, FTYPE CUf, FTYPE *p_c, FTYPE *p_l, FTYPE *p_r, FTYPE *F, FTYPE *FEM, FTYPE *ctopallptr)
130 {
131  int flux_compute(int i, int j, int k, int dir, struct of_geom *geom, FTYPE *cminmax_l, FTYPE *cminmax_r, FTYPE *cminmax, FTYPE ctopmhd, FTYPE *cminmaxrad_l, FTYPE *cminmaxrad_r, FTYPE *cminmaxrad, FTYPE ctoprad, FTYPE CUf, FTYPE *p_l, FTYPE *p_r, FTYPE *U_l, FTYPE *U_r, FTYPE *F_l, FTYPE *F_r, FTYPE *F);
132  int p2SFUevolve_splitmaem(int dir, int isleftright, FTYPE *p, struct of_geom *geom, struct of_state **ptrstate, FTYPE *F, FTYPE *FEM, FTYPE *U, FTYPE *UEM);
133  FTYPE cminmax_l[NUMCS], cminmax_r[NUMCS], cminmax[NUMCS];
134  FTYPE cminmaxrad_l[NUMCS], cminmaxrad_r[NUMCS], cminmaxrad[NUMCS];
135  FTYPE cminmaxrad2_l[NUMCS], cminmaxrad2_r[NUMCS], cminmaxrad2[NUMCS];
136  FTYPE ctopmhd,ctoprad,ctoprad2;
137  struct of_state state_c, state_l, state_r;
138  struct of_state *ptrstate_c, *ptrstate_l, *ptrstate_r;
139  FTYPE F_c[NPR], F_l[NPR], F_r[NPR];
140  FTYPE U_c[NPR], U_l[NPR], U_r[NPR];
141  FTYPE FEM_c[NPR], FEM_l[NPR], FEM_r[NPR];
142  FTYPE UEM_c[NPR], UEM_l[NPR], UEM_r[NPR];
143  int pl,pliter;
144  int jj,kk;
145 
146 
147  // default
148  ptrstate_c = &state_c;
149  ptrstate_l = &state_l;
150  ptrstate_r = &state_r;
151 
152 
154  //
155  // setup flux calculation based upon interpolated primitives
156  //
158  MYFUN(p2SFUevolve_splitmaem(dir, ISLEFT, p_l, ptrgeom, &ptrstate_l, F_l, FEM_l, U_l, UEM_l),"step_ch.c:fluxcalc()", "p2SFUevolve()", 1);
159  MYFUN(p2SFUevolve_splitmaem(dir, ISRIGHT, p_r, ptrgeom, &ptrstate_r, F_r, FEM_r, U_r, UEM_r),"step_ch.c:fluxcalc()", "p2SFUevolve()", 2);
160 
161 
162  // usually "always" need cminmax_l cminmax_r and always need ctop
163  get_wavespeeds(dir, ptrgeom, p_l, p_r, U_l, U_r, F_l, F_r, ptrstate_l, ptrstate_r, cminmax_l, cminmax_r, cminmax, &ctopmhd, cminmaxrad_l, cminmaxrad_r, cminmaxrad, &ctoprad, cminmaxrad2_l, cminmaxrad2_r, cminmaxrad2, &ctoprad2);
164 
165 
167  // choose ctopall that will go into choosing timestep
168  if(EOMRADTYPE!=EOMRADNONE){
169  // get maximum over both mhd and radiation for setting dt later
170  *ctopallptr=MAX(ctopmhd,ctoprad2); // choose ctoprad2 for setting timestep
171  }
172  else *ctopallptr=ctopmhd;
173 
174 
175  // cminmaxrad and ctoprad are used for fluxes (not "2" version that's only for setting timestep)
176 
177  // GODMARK:
178  // below assumes flux_compute is linear in FMA and FEM, which it is right now
179  // if used Riemann solver, would need to have it output FMA and FEM parts
180  MYFUN(flux_compute(i, j, k, dir, ptrgeom, cminmax_l,cminmax_r, cminmax, ctopmhd, cminmaxrad_l,cminmaxrad_r, cminmaxrad, ctoprad, CUf, p_l, p_r, U_l, U_r, F_l, F_r, F),"step_ch.c:fluxcalc()", "flux_compute", 1);
181 
182  MYFUN(flux_compute(i, j, k, dir, ptrgeom, cminmax_l,cminmax_r, cminmax, ctopmhd, cminmaxrad_l,cminmaxrad_r, cminmaxrad, ctoprad, CUf, p_l, p_r, UEM_l, UEM_r, FEM_l, FEM_r, FEM),"step_ch.c:fluxcalc()", "flux_compute", 1);
183 
184 
185  // Note that if splitmaem==1 then pressure term was moved as separate quasi-flux/conserved term. flux remains, and conserved quantity not used
186 
187 
188 #if(FLUXDUMP==1)
189  diag_fluxdump_1(dir,i,j,k,p_l, p_r, F_l, F_r);
190 #endif
191 
192 
193  return(0);
194 }
195 
196 
197 
199 int p2SFUevolve(int dir, int isleftright, FTYPE *p, struct of_geom *ptrgeom, struct of_state **ptrstate, FTYPE *F, FTYPE *U)
200 {
201  MYFUN(get_stateforfluxcalc(dir,isleftright, p, ptrgeom, ptrstate),"flux.c:p2SFUevolve", "get_state()", 1);
202 
203 
204  // DEBUG:
205  // if(ptrgeom->i==26 && ptrgeom->j==40 && dir==1){
206  // dualfprintf(fail_file,"NORMAL: INp2SFUevolve: gdet=%21.15g : B1=%21.15g B2=%21.15g B3=%21.15g uu0=%21.15g uu1=%21.15g uu2=%21.15g uu3=%21.15g\n",ptrgeom->gdet,p[B1],p[B2],p[B3],(*ptrstate)->ucon[TT],(*ptrstate)->ucon[1],(*ptrstate)->ucon[2],(*ptrstate)->ucon[3]);
207  // }
208 
209 
210  MYFUN(primtoflux(UEVOLVE,p, *ptrstate, dir, ptrgeom, F, NULL),"flux.c:p2SFUevolve()","primtoflux_calc() dir=1/2 l", 1);
211  MYFUN(primtoflux(UEVOLVE,p, *ptrstate, TT, ptrgeom, U, NULL),"flux.c:p2SFUevolve()", "primtoflux_calc() dir=l0", 1);
212 
213  return(0);
214 }
215 
216 
217 // p2SFUevolve() for each physical term
218 int p2SFUevolveall(int dir, int isleftright, FTYPE *p, struct of_geom *ptrgeom, struct of_state **ptrstate, FTYPE *FPAKE, FTYPE *FEN, FTYPE *FEM, FTYPE *FRAD, FTYPE *UPAKE, FTYPE *UEN, FTYPE *UEM, FTYPE *URAD)
219 {
220 
221 
222  // DEBUG:
223  // if(ptrgeom->i==26 && ptrgeom->j==40 && dir==1){
224  // dualfprintf(fail_file,"NORMAL: INp2SFUevolve: gdet=%21.15g : B1=%21.15g B2=%21.15g B3=%21.15g uu0=%21.15g uu1=%21.15g uu2=%21.15g uu3=%21.15g\n",ptrgeom->gdet,p[B1],p[B2],p[B3],(*ptrstate)->ucon[TT],(*ptrstate)->ucon[1],(*ptrstate)->ucon[2],(*ptrstate)->ucon[3]);
225  // }
226 
227  MYFUN(get_stateforfluxcalc(dir,isleftright, p, ptrgeom, ptrstate),"flux.c:p2SFUevolve", "get_state()", 1);
228 
229  int pliter,pl,jj;
230  struct of_state qtemp;
231  FTYPE ptemp[NPR];
232 
233  PLOOP(pliter,pl) ptemp[pl] = p[pl]; qtemp=**ptrstate;
234  ptemp[UU]=qtemp.pressure=qtemp.entropy=SMALL; //ug
235  ptemp[B1]=ptemp[B2]=ptemp[B3]=0.0; DLOOPA(jj) qtemp.bcon[jj]=qtemp.bcov[jj]=0.0; qtemp.bsq=0.0; // b
236  if(URAD0>=0) ptemp[URAD0]=SMALL; // urad
237  MYFUN(primtoflux(UEVOLVE,ptemp, &qtemp, dir, ptrgeom, FPAKE,NULL),"flux.c:p2SFUevolve()","primtoflux_calc() dir=1/2 l", 1);
238  MYFUN(primtoflux(UEVOLVE,ptemp, &qtemp, TT, ptrgeom, UPAKE,NULL),"flux.c:p2SFUevolve()", "primtoflux_calc() dir=l0", 1);
239 
240  ptemp[RHO]=SMALL; // rho
241  ptemp[B1]=ptemp[B2]=ptemp[B3]=0.0; DLOOPA(jj) qtemp.bcon[jj]=qtemp.bcov[jj]=0.0; qtemp.bsq=0.0; // b
242  if(URAD0>=0) ptemp[URAD0]=SMALL; // urad
243  MYFUN(primtoflux(UEVOLVE,ptemp, &qtemp, dir, ptrgeom, FEN,NULL),"flux.c:p2SFUevolve()","primtoflux_calc() dir=1/2 l", 1);
244  MYFUN(primtoflux(UEVOLVE,ptemp, &qtemp, TT, ptrgeom, UEN,NULL),"flux.c:p2SFUevolve()", "primtoflux_calc() dir=l0", 1);
245 
246  ptemp[RHO]=SMALL; // rho
247  ptemp[UU]=qtemp.pressure=qtemp.entropy=SMALL; //ug
248  if(URAD0>=0) ptemp[URAD0]=SMALL; // urad
249  MYFUN(primtoflux(UEVOLVE,ptemp, &qtemp, dir, ptrgeom, FEM,NULL),"flux.c:p2SFUevolve()","primtoflux_calc() dir=1/2 l", 1);
250  MYFUN(primtoflux(UEVOLVE,ptemp, &qtemp, TT, ptrgeom, UEM,NULL),"flux.c:p2SFUevolve()", "primtoflux_calc() dir=l0", 1);
251 
252  if(EOMRADTYPE!=EOMRADNONE){
253  ptemp[RHO]=SMALL; // rho
254  ptemp[UU]=qtemp.pressure=qtemp.entropy=SMALL; //ug
255  ptemp[B1]=ptemp[B2]=ptemp[B3]=0.0; DLOOPA(jj) qtemp.bcon[jj]=qtemp.bcov[jj]=0.0; qtemp.bsq=0.0; // b
256  MYFUN(primtoflux(UEVOLVE,ptemp, &qtemp, dir, ptrgeom, FRAD,NULL),"flux.c:p2SFUevolve()","primtoflux_calc() dir=1/2 l", 1);
257  MYFUN(primtoflux(UEVOLVE,ptemp, &qtemp, TT, ptrgeom, URAD,NULL),"flux.c:p2SFUevolve()", "primtoflux_calc() dir=l0", 1);
258  }
259 
260 
261  return(0);
262 }
263 
264 
266 int p2SFUevolve_splitmaem(int dir, int isleftright, FTYPE *p, struct of_geom *ptrgeom, struct of_state **ptrstate, FTYPE *F, FTYPE *FEM, FTYPE *U, FTYPE *UEM)
267 {
268 
269 #if(ROEAVERAGEDWAVESPEED && ((USESTOREDSPEEDSFORFLUX==0)||(STOREWAVESPEEDS==0)) )
270  dualfprintf(fail_file,"Only supposed to use this splitting of F and U pulling out particular linear terms if not assuming U is not true conserved quantity\n");
271  myexit(72679262);
272 #endif
273 
274  MYFUN(get_stateforfluxcalc(dir,isleftright, p, ptrgeom, ptrstate),"flux.c:p2SFUevolve", "get_state()", 1);
275 
276  MYFUN(primtoflux_splitmaem(UEVOLVE, p, *ptrstate, dir, dir, ptrgeom, F, FEM),"flux.c:p2SFUevolve_splitmaem()","primtoflux_ma() dir=1/2 l", 1);
277  MYFUN(primtoflux_splitmaem(UEVOLVE, p, *ptrstate, dir, TT, ptrgeom, U, UEM),"flux.c:p2SFUevolve_splitmaem()", "primtoflux_ma() dir=l0", 1);
278 
279 
280 
281 
282  return(0);
283 }
284 
285 
286 
287 
289 void diag_fluxdump_1(int dir, int i, int j, int k, FTYPE *p_l, FTYPE *p_r, FTYPE *F_l, FTYPE *F_r)
290 {
291  int pl,pliter;
292 
293 #if(FLUXDUMP==1)
294  PLOOP(pliter,pl) GLOBALMACP0A1(fluxdump,i,j,k,4*NPR + (dir-1)*NPR*5 + NPR*1 + pl)=F_l[pl];
295  PLOOP(pliter,pl) GLOBALMACP0A1(fluxdump,i,j,k,4*NPR + (dir-1)*NPR*5 + NPR*2 + pl)=F_r[pl];
296  PLOOP(pliter,pl) GLOBALMACP0A1(fluxdump,i,j,k,4*NPR + (dir-1)*NPR*5 + NPR*3 + pl)=p_l[pl];
297  PLOOP(pliter,pl) GLOBALMACP0A1(fluxdump,i,j,k,4*NPR + (dir-1)*NPR*5 + NPR*4 + pl)=p_r[pl];
298 
299  if(0 && dir==2 && j==0 && i==N1/2){ // just pick one point on polar axis
300  PLOOP(pliter,pl){
301  // dualfprintf(fail_file,"n=%ld sp=%d pl=%d :: F_l=%21.15g F_r=%21.15g U_l=%21.15g U_r=%21.15g p_l=%21.15g p_r=%21.15g\n",nstep,steppart,pl,F_l[pl],F_r[pl],U_l[pl],U_r[pl],p_l[pl],p_r[pl]);
302  }
303  }
304 #endif
305 }
306 
307 
308 
309 
310 
311 
313 #define HLLCOMPUTE(cmin,cmax,U_l,U_r,F_l,F_r) (( (cmax * F_l + cmin * F_r) - (FLUXDISSIPATION)*(cmax * cmin) * (U_r - U_l) ) / (cmax + cmin) )
314 
316 #define LAXFCOMPUTE(ctop,U_l,U_r,F_l,F_r) (0.5 * ( (F_l + F_r) - (FLUXDISSIPATION)*ctop * (U_r - U_l) ) )
317 //#define LAXFCOMPUTE(ctop,U_l,U_r,F_l,F_r) (0.5 * ( (F_l + F_r) ) )
318 //#define LAXFCOMPUTE(ctop,U_l,U_r,F_l,F_r) (0.5 * ( (F_l + F_r) - (FLUXDISSIPATION)*fabs(ctop) * fabs(U_r - U_l) ) )
319 //#define LAXFCOMPUTE(ctop,U_l,U_r,F_l,F_r) (0.5 * ( (F_l + F_r) - (FLUXDISSIPATION)*1E-5 * (U_r - U_l) ) )
320 
322 #define LFCOMPUTE(crus,U_l,U_r,F_l,F_r) (0.5*( (F_l+F_r) - (FLUXDISSIPATION)*crus*(U_r - U_l)))
323 
325 #define UHALF(crus,U_l, U_r, F_l, F_r) (0.5*((U_l+U_r)-(FLUXDISSIPATION)*(F_r-F_l)/crus))
326 
328 #define LWCOMPUTE(crus,U_l,U_r,F_l,Fhalf,F_r) (Fhalf)
329 
331 #define FORCECOMPUTE(crus,U_l,U_r,F_l,Fhalf,F_r) (0.25*( (F_l+F_r) + 2.0*Fhalf - (FLUXDISSIPATION)*crus*(U_r - U_l)))
332 
334 #define GFORCECOMPUTE(cg,crus,U_l,U_r,F_l,Fhalf,F_r) ( (LWCOMPUTE(crus,U_l,U_r,F_l,Fhalf,F_r) + cg*LFCOMPUTE(crus,U_l,U_r,F_l,F_r))/(1.0+cg) )
335 
336 
338 #define HLLLAXF1(cmin,cmax,ctop,f_s, U_l, U_r, F_l, F_r) ( HLLCOMPUTE(cmin,cmax,U_l,U_r,F_l,F_r)*f_s + LAXFCOMPUTE(ctop,U_l,U_r,F_l,F_r)*(1.0-f_s) )
339 
340 
341 
342 
343 
346 int flux_compute(int i, int j, int k, int dir, struct of_geom *geom, FTYPE *cminmax_l, FTYPE *cminmax_r, FTYPE *cminmax, FTYPE ctopmhd, FTYPE *cminmaxrad_l, FTYPE *cminmaxrad_r, FTYPE *cminmaxrad, FTYPE ctoprad, FTYPE CUf, FTYPE *p_l, FTYPE *p_r, FTYPE *U_l, FTYPE *U_r, FTYPE *F_l, FTYPE *F_r, FTYPE *F)
347 {
348  int pl,pliter;
349  void choose_flux(int fluxmethodlocal, int i, int j, int k, int pl, FTYPE *laxffrac,FTYPE *hllfrac);
350  FTYPE laxffrac[NPR],hllfrac[NPR];
351  int cminmax_calc(FTYPE cmin_l,FTYPE cmin_r,FTYPE cmax_l,FTYPE cmax_r,FTYPE *cmin,FTYPE *cmax,FTYPE *ctop);
352  int forceflux_compute(int dir,struct of_geom *geom, FTYPE *cmin, FTYPE *cmax, FTYPE *ctop, FTYPE *cforce, FTYPE *p_l, FTYPE *p_r, FTYPE *U_l, FTYPE *U_r,FTYPE *F_l,FTYPE *F_r, FTYPE *F);
353  int mustaflux_compute(int dir,struct of_geom *geom, FTYPE *cmin_l, FTYPE *cmin_r, FTYPE *cmax_l, FTYPE *cmax_r, FTYPE *cmin, FTYPE *cmax, FTYPE *ctop, FTYPE *cforce, FTYPE *p_l, FTYPE *p_r, FTYPE *U_l, FTYPE *U_r,FTYPE *F_l,FTYPE *F_r, FTYPE *F);
354  int hllflux_compute(int dir,struct of_geom *geom, FTYPE *cmin, FTYPE *cmax, FTYPE *ctop, FTYPE *p_l, FTYPE *p_r, FTYPE *U_l, FTYPE *U_r,FTYPE *F_l,FTYPE *F_r, FTYPE *F);
355  FTYPE crus[NPR];
356  FTYPE cforce[NPR];
357  FTYPE cmin_l[NPR], cmax_l[NPR], cmin_r[NPR], cmax_r[NPR], cmax[NPR], cmin[NPR];
358  FTYPE ctop[NPR];
359  FTYPE dPoP,f_s;
360 
361  // defaults
362  int fluxmethodlocal[NPR];
363  PALLLOOP(pl) fluxmethodlocal[pl]=fluxmethod[pl];
364  FTYPE cmaxfactor=1.0;
365  FTYPE cminfactor=1.0;
366 #if(OUTERRADIALSUPERFAST)
367  // force effective superfast condition on outer radial boundaries
368  if(ISBLACKHOLEMCOORD(MCOORD)){
369  if(dir==1 && startpos[1]+i==totalsize[1]){ PALLLOOP(pl) fluxmethodlocal[pl]=HLLFLUX; cminfactor=0.0; }
370  if(dir==1 && startpos[1]+i==0){ PALLLOOP(pl) fluxmethodlocal[pl]=HLLFLUX; cmaxfactor=0.0; }
371  }
372 #endif
373 
374  // assign cmin/cmax/ctop for each conserved quantity
375  PLOOP(pliter,pl){
376  if(RADFULLPL(pl)==0){
377  // MHD terms
378  cmin_l[pl] = cminmax_l[CMIN];
379  cmax_l[pl] = cminmax_l[CMAX];
380  cmin_r[pl] = cminmax_r[CMIN];
381  cmax_r[pl] = cminmax_r[CMAX];
382  cmin[pl] = cminmax[CMIN]*cminfactor;
383  cmax[pl] = cminmax[CMAX]*cmaxfactor;
384  ctop[pl] = ctopmhd;
385 
386  }
387  else{
388  // radiation terms
389  cmin_l[pl] = cminmaxrad_l[CMIN];
390  cmax_l[pl] = cminmaxrad_l[CMAX];
391  cmin_r[pl] = cminmaxrad_r[CMIN];
392  cmax_r[pl] = cminmaxrad_r[CMAX];
393  cmin[pl] = cminmaxrad[CMIN]*cminfactor;
394  cmax[pl] = cminmaxrad[CMAX]*cmaxfactor;
395  ctop[pl] = ctoprad;
396  }
397  // if(steppart==0) crus[pl]=dx[dir]/(dt*0.5);
398  // else crus[pl]=dx[dir]/(dt);
399  crus[pl]=dx[dir]/(dt*CUf);
400  // crus[pl]=dx[dir]/(dt);
401  // dualfprintf(fail_file,"CUf=%21.15g dt=%21.15g dx[%d]=%21.15g\n",CUf,dt,dir,dx[dir]);
402  // crus[pl]=dx[dir]/(dt);
403  }
404 
405 
407  //
408  // store non-dissipative contribution
409  //
411  int plsp;
412  if(NSPECIAL>=1){
413  plsp=NPR;
414  F[plsp] = 0.5 * (F_l[SPECIALPL1] + F_r[SPECIALPL1]);
415  }
416  if(NSPECIAL>=2){
417  plsp=NPR+1;
418  F[plsp] = 0.5 * (F_l[SPECIALPL2] + F_r[SPECIALPL2]);
419  }
420  if(NSPECIAL>=3){
421  plsp=NPR+2;
422  F[plsp] = 0.5 * (F_l[SPECIALPL3] + F_r[SPECIALPL3]);
423  }
424  if(NSPECIAL>=4){
425  plsp=NPR+3;
426  F[plsp] = 0.5 * (F_l[SPECIALPL4] + F_r[SPECIALPL4]);
427  }
428  if(NSPECIAL>=5){
429  plsp=NPR+4;
430  F[plsp] = 0.5 * (F_l[SPECIALPL5] + F_r[SPECIALPL5]);
431  }
432  if(NSPECIAL>=6){
433  plsp=NPR+5;
434  F[plsp] = 0.5 * (F_l[SPECIALPL6] + F_r[SPECIALPL6]);
435  }
436 
437 
439  //
440  // store full non-dissipative + dissipative contributions
441  //
443  if(fluxmethodlocal[RHO]==LAXFFLUX || fluxmethodlocal[RHO]==HLLFLUX){
444  PLOOP(pliter,pl){
445  if(fluxmethodlocal[pl]==LAXFFLUX) F[pl] = LAXFCOMPUTE(ctop[pl],U_l[pl],U_r[pl],F_l[pl],F_r[pl]);
446  else hllflux_compute(dir,geom,cmin,cmax,ctop,p_l,p_r,U_l,U_r,F_l,F_r,F);
447  }
448  }
449  else if(fluxmethodlocal[RHO]==HLLFLUX){
451  //
452  // decide which flux formula to use
453  //
455  // PLOOP(pliter,pl) choose_flux(fluxmethodlocal, i,j,k,pl,laxffrac,hllfrac);
456 
457 
459  //
460  // Decide if using different flux calculation on boundary zones
461  //
463  //#if(HLLBOUNDARY)
464  // if(
465  // ((dir == 3) && ( ((startpos[3]+k == 0)&&(BCtype[X3DN]==OUTFLOW)) || ((startpos[3]+k == totalsize[3])&&(BCtype[X3UP]==OUTFLOW)) ))||
466  // ((dir == 2) && ( ((startpos[2]+j == 0)&&(BCtype[X2DN]==POLARAXIS)) || ((startpos[2]+j == totalsize[2])&&(BCtype[X2UP]==POLARAXIS)) ))||
467  // ((dir == 1) && ( ((startpos[1]+i == 0)&&(BCtype[X1DN]==OUTFLOW)) || ((startpos[1]+i == totalsize[1])&&(BCtype[X1UP]==OUTFLOW)) ))
468  // )
469  // {
470  // PLOOP(pliter,pl){
471  // hllfrac[pl]=1.0;
472  // laxffrac[pl]=0.0;
473  // }
474  // }
475  //#endif
476 
477  hllflux_compute(dir,geom,cmin,cmax,ctop,p_l,p_r,U_l,U_r,F_l,F_r,F);
478 
479  }
480  else if(fluxmethodlocal[RHO]==FORCEFLUX){
481 
482  // normal Rusanov speed
483  //cforce=crus;
484  // LAXF speed
485  PLOOP(pliter,pl) cforce[pl]=ctop[pl];
486 
487 
488  forceflux_compute(dir,geom,cmin,cmax,ctop,cforce,p_l,p_r,U_l,U_r,F_l,F_r,F);
489 
490  }
491  else if(fluxmethodlocal[RHO]==MUSTAFLUX){
492 
493  PLOOP(pliter,pl) {
494  //cforce[pl]=crus[pl];
495  cforce[pl]=ctop[pl];
496  }
497  mustaflux_compute(dir,geom,cmin_l,cmin_r,cmax_l,cmax_r,cmin,cmax,ctop,cforce,p_l,p_r,U_l,U_r,F_l,F_r,F);
498 
499  }
500  else if(fluxmethodlocal[RHO]==HLLLAXF1FLUX){
501 
502 #define MINDPOP (0.0)
503 #define MAXDPOP (0.3)
504 
505  dPoP = fabs( (p_r[UU]-p_l[UU])/(0.5*(p_r[UU]+p_l[UU])) );
506  f_s = (1.0-0.0)/(MAXDPOP-MINDPOP) * (dPoP - MINDPOP) + 0.0;
507  if(f_s>1.0) f_s=1.0;
508  if(f_s<0.0) f_s=0.0;
509 
510  PLOOP(pliter,pl){
511  if(cmax[pl]+cmin[pl]!=0.0){
512  F[pl] = HLLLAXF1(cmin[pl],cmax[pl],ctop[pl],f_s, U_l[pl], U_r[pl], F_l[pl], F_r[pl]);
513  }
514  else{
515  F[pl] = LAXFCOMPUTE(ctop[pl],U_l[pl],U_r[pl],F_l[pl],F_r[pl]);
516  }
517  }
518 
519 
520 
521  //PLOOP(pliter,pl) dualfprintf(fail_file,"p_l[%d]=%g p_r[%d]=%g U_l[%d]=%g U_r[%d]=%g F_l[%d]=%g F_r[%d]=%g :: F[%d]=%g\n",pl,p_l[pl],pl,p_r[pl],pl,U_l[pl],pl,U_r[pl],pl,F_l[pl],pl,F_r[pl],pl,F[pl]);
522 
523  // dualfprintf(fail_file,"cmin=%g cmax=%g ctop=%g\n",cmin,cmax,ctop);
524 
525  // dualfprintf(fail_file,"i=%d dPoP=%g f_s=%g :: %g %g\n",geom->i,dPoP,f_s,p_l[UU],p_r[UU]);
526 
527 
528  }
529 
530 
531 
532 
533 
534 
535 
536 
537 
538  return(0);
539 }
540 
541 
542 
544 #define USE_CORRECTED_STATES 0
545 
547 int hllflux_compute(int dir,struct of_geom *geom, FTYPE *cmin, FTYPE *cmax, FTYPE *ctop, FTYPE *p_l, FTYPE *p_r, FTYPE *U_l, FTYPE *U_r,FTYPE *F_l,FTYPE *F_r, FTYPE *F)
548 {
549  int pl,pliter;
550  FTYPE vmin[NPR],vmax[NPR];
551  FTYPE cminreal[NPR],cmaxreal[NPR];
552 
553  PLOOP(pliter,pl) {
554 
555 
556 #if(USE_CORRECTED_STATES)
557  // get vmin/vmax
558  vmin[pl]=min(p_l[UU+dir],p_r[UU+dir]);
559  vmax[pl]=max(p_l[UU+dir],p_r[UU+dir]);
560  vmin[pl] = fabs(max(0., -vmin[pl]));
561  vmax[pl] = fabs(max(0., vmax[pl]));
562 
563  // get HLL+correction
564  if(cmax[pl]+cmin[pl]!=0.0){
565  if( pl == UU + dir && (vmin[pl]+vmax[pl]!=0.0) ) {
566  // if(vmin[pl]+vmax[pl]!=0.0){
567  cminreal[pl] = vmin[pl];
568  cmaxreal[pl] = vmax[pl];
569  }
570  else {
571  cminreal[pl] = cmin[pl];
572  cmaxreal[pl] = cmax[pl];
573  }
574  F[pl] = HLLCOMPUTE(cminreal[pl],cmaxreal[pl],U_l[pl],U_r[pl],F_l[pl],F_r[pl]);
575  }
576  else F[pl] = LAXFCOMPUTE(ctop[pl],U_l[pl],U_r[pl],F_l[pl],F_r[pl]);
577 
578 #else
579 
580  if(cmax[pl]+cmin[pl]!=0.0){
581  F[pl] = HLLCOMPUTE(cmin[pl],cmax[pl],U_l[pl],U_r[pl],F_l[pl],F_r[pl]);
582  }
583  else F[pl] = LAXFCOMPUTE(ctop[pl],U_l[pl],U_r[pl],F_l[pl],F_r[pl]);
584 
585 #endif
586 
587  } // over pl's
588 
589  return(0);
590 
591 }
592 
593 
594 #define MUSTACOEF (1.0)
595 
596 
597 
598 //#define MUSTACOEF (1.0)
599 // HLL+0.9 gives nice and smooth Noh, but blurs stationary contact
600 
601 //#define MUSTACOEF (0.49)
602 // HLL+0.49 seems no better than without MUSTA
603 
604 //#define MUSTACOEF (0.4) // ok test2 velocity with MUSTAHLL
605 
606 //#define MUSTACOEF (0.9) // ok test2 velocity with MUSTAHLL
607 
608 //#define MUSTACOEF (0.96) // ok test2 velocity with MUSTAHLL
609 
610 //#define MUSTACOEF (0.99) // death with MUSTAHLL
611 
612 //#define MUSTACOEF (0.999) // death with MUSTAHLL
613 
614 //#define MUSTACOEF (0.9990)
615 // 0.9 gives nice and smooth Noh, but blurs stationary contact
616 // 0.99 gives oscillatory Noh and blurs stationary a bit
617 // 0.9990 gives good double ratefaction as long as a2c is turned off, and for a2c on or off a nearly stationary contact and decent Noh
618 
619 
620 //#define MUSTACOEF (2.0)
621 // MUSTAFLUX==MUSTAFORCE w/ MUSTACOEF==2.0 has stationary contact
622 // does ok with rarefaction, but fails a bit like MUSTAHLL w/ MUSTACOEF=1.0
623 // bit more oscillatory in flat region for Noh
624 
626 int forceflux_compute(int dir,struct of_geom *geom, FTYPE *cmin, FTYPE *cmax, FTYPE *ctop, FTYPE *cforce, FTYPE *p_l, FTYPE *p_r, FTYPE *U_l, FTYPE *U_r,FTYPE *F_l,FTYPE *F_r, FTYPE *F)
627 {
628  int pl,pliter;
629  FTYPE umid[NPR],pmid[NPR],fmid[NPR];
630  struct of_state state;
631  struct of_state *ptrstate;
632  int doforceflux;
633  FTYPE vmin[NPR],vmax[NPR];
634  FTYPE vminorig[NPR], vmaxorig[NPR];
635  FTYPE cminreal[NPR], cmaxreal[NPR];
636  int hllflux_compute(int dir,struct of_geom *geom, FTYPE *cmin, FTYPE *cmax, FTYPE *ctop, FTYPE *p_l, FTYPE *p_r, FTYPE *U_l, FTYPE *U_r,FTYPE *F_l,FTYPE *F_r, FTYPE *F);
637  struct of_newtonstats newtonstats; setnewtonstatsdefault(&newtonstats);
638  int showmessages=1;
639  int allowlocalfailurefixandnoreport=1;
640 
642  //
643  // get middle umid
644  //
646  PLOOP(pliter,pl) umid[pl] = UHALF(cforce[pl],U_l[pl],U_r[pl],F_l[pl],F_r[pl]);
647 
648  // PLOOP(pliter,pl) dualfprintf(fail_file,"%d : cforce=%21.15g : U_l[%d]=%21.15g U_r[%d]=%21.15g umid[%d]=%21.15g\n",geom->i,cforce,pl,U_l[pl],pl,U_r[pl],pl,umid[pl]);
649 
650  // set guess for inversion
651  PLOOP(pliter,pl) pmid[pl]=0.5*(p_l[pl]+p_r[pl]);
652  // get primitive pmid(umid)
653  int eomtype=EOMDEFAULT;
654  FTYPE dissmeasure=-1.0; // assume energy try ok
655  int whichcap=CAPTYPEBASIC;
656  int whichmethod=MODEDEFAULT;
657  int modprim=0;
658  int checkoninversiongas=CHECKONINVERSION;
659  int checkoninversionrad=CHECKONINVERSIONRAD;
660  MYFUN(Utoprimgen(showmessages,checkoninversiongas,checkoninversionrad,allowlocalfailurefixandnoreport, 0, &eomtype,whichcap,whichmethod,modprim,EVOLVEUTOPRIM,UEVOLVE, umid, NULL, geom, dissmeasure, pmid, pmid,&newtonstats),"flux.c:flux_compute()", "Utoprimgen", 1);
661  doforceflux=1;
662  if(GLOBALMACP0A1(pflag,geom->i,geom->j,geom->k,FLAGUTOPRIMFAIL)){
663  if(debugfail>=1) dualfprintf(fail_file,"Failed to find inversion for FORCEFLUX, trying p_l : nstep=%ld t=%21.15g i=%d j=%d k=%d\n",nstep,t,geom->i,geom->j,geom->k);
664  PLOOP(pliter,pl) pmid[pl]=p_l[pl];
665  // get primitive pmid(umid)
666  MYFUN(Utoprimgen(showmessages,checkoninversiongas,checkoninversionrad,allowlocalfailurefixandnoreport, 0, &eomtype,whichcap,whichmethod,modprim,EVOLVEUTOPRIM,UEVOLVE, umid, NULL, geom, dissmeasure, pmid, pmid,&newtonstats),"flux.c:flux_compute()", "Utoprimgen", 1);
667  if(GLOBALMACP0A1(pflag,geom->i,geom->j,geom->k,FLAGUTOPRIMFAIL)){
668  if(debugfail>=1) dualfprintf(fail_file,"Failed to find inversion for FORCEFLUX, trying p_r : nstep=%ld t=%21.15g i=%d j=%d k=%d\n",nstep,t,geom->i,geom->j,geom->k);
669  PLOOP(pliter,pl) pmid[pl]=p_r[pl];
670  // get primitive pmid(umid)
671  MYFUN(Utoprimgen(showmessages,checkoninversiongas,checkoninversionrad,allowlocalfailurefixandnoreport, 0, &eomtype,whichcap,whichmethod,modprim,EVOLVEUTOPRIM,UEVOLVE, umid, NULL, geom, dissmeasure, pmid, pmid,&newtonstats),"flux.c:flux_compute()", "Utoprimgen", 1);
672  if(GLOBALMACP0A1(pflag,geom->i,geom->j,geom->k,FLAGUTOPRIMFAIL)){
673  if(debugfail>=1) dualfprintf(fail_file,"No initial guess worked, rejecting FORCEFLUX method : nstep=%ld t=%21.15g i=%d j=%d k=%d\n",nstep,t,geom->i,geom->j,geom->k);
674  doforceflux=0;
675  }
676  }
677  }
678 
679  // PLOOP(pliter,pl) dualfprintf(fail_file,"%d : pmid[%d]=%21.15g p_l[%d]=%21.15g p_r[%d]=%21.15g\n",geom->i,pl,pmid[pl],pl,p_l[pl],pl,p_r[pl]);
680 
682  //
683  // get flux for pmid Fmid(pmid)
684  //
686  if(doforceflux){
687  ptrstate=&state; // default
688  MYFUN(get_stateforfluxcalc(dir,ISMIDDLE, pmid, geom, &ptrstate),"flux.c:flux_compute()", "get_state()", 1);
689  MYFUN(primtoflux(UEVOLVE, pmid, ptrstate, dir, geom, fmid, NULL),"flux.c:flux_compute()","primtoflux_calc() dir=1/2 l", 1);
690 
691  // compute FORCE flux
692  // PLOOP(pliter,pl) F[pl] = FORCECOMPUTE(cforce,U_l[pl],U_r[pl],F_l[pl],fmid[pl],F_r[pl]);
693  PLOOP(pliter,pl) F[pl] = GFORCECOMPUTE(MUSTACOEF,cforce[pl],U_l[pl],U_r[pl],F_l[pl],fmid[pl],F_r[pl]);
694 
695  // PLOOP(pliter,pl) dualfprintf(fail_file,"%d : F[%d]=%21.15g Fhll[%d]=%21.15g\n",geom->i,pl,F[pl],pl,HLLCOMPUTE(cmin,cmax,U_l[pl],U_r[pl],F_l[pl],F_r[pl]));
696  }
697  else{
698  // reduce to HLL
699  hllflux_compute(dir,geom,cmin,cmax,ctop,p_l,p_r,U_l,U_r,F_l,F_r,F);
700  }
701 
702  return(0);
703 }
704 
705 
706 
708 #define NUMMUSTAITERS 1
709 
711 #define DOMULTICELL 0
712 
713 #define NUMMULTIMUSTAITERS 2 // 1=expensive way to do no musta, 2 minimum for interesting calculation
714 
715 
716 
717 #define NUMLOCALCELLS 2
718 #define MULTIMUSTACOEF (0.9)
719 
720 
721 #define MUSTAFORCE 0
722 #define MUSTALAXF 1
723 #define MUSTAHLL 2
724 
725 #define WHICHFLUX MUSTAHLL
726 //#define WHICHFLUX MUSTALAXF
727 //#define WHICHFLUX MUSTAFORCE
728 // MUSTAFORCE worse than MUSTAHLL. Contacts are blurred (at one point this wasn't true!)
729 // MUSTAFORCE doesn't fail on rarefaction however, while MUSTAHLL does a bit
730 // MUSTAFORCE behaves better with a2c in FV method
731 
732 
733 // whether to limit MUSTA by some fraction of HLL
734 #define HLLBOUNDMUSTA 0
735 
737 int mustaflux_compute(int dir,struct of_geom *geom, FTYPE *cmin_l, FTYPE *cmin_r, FTYPE *cmax_l, FTYPE *cmax_r, FTYPE *cmin, FTYPE *cmax, FTYPE *ctop, FTYPE *cforce, FTYPE *p_l, FTYPE *p_r, FTYPE *U_l, FTYPE *U_r,FTYPE *F_l,FTYPE *F_r, FTYPE *F)
738 {
739  int musta1flux_compute(int dir,struct of_geom *geom, FTYPE *cmin_l, FTYPE *cmin_r, FTYPE *cmax_l, FTYPE *cmax_r, FTYPE *cmin, FTYPE *cmax, FTYPE *ctop, FTYPE *cforce, FTYPE *p_l, FTYPE *p_r, FTYPE *U_l, FTYPE *U_r,FTYPE *F_l,FTYPE *F_r, FTYPE *F);
740  int musta2flux_compute(int dir,struct of_geom *geom, FTYPE *cmin_l, FTYPE *cmin_r, FTYPE *cmax_l, FTYPE *cmax_r, FTYPE *cmin, FTYPE *cmax, FTYPE *ctop, FTYPE *cforce, FTYPE *p_l, FTYPE *p_r, FTYPE *U_l, FTYPE *U_r,FTYPE *F_l,FTYPE *F_r, FTYPE *F);
741  int hllflux_compute(int dir,struct of_geom *geom, FTYPE *cmin, FTYPE *cmax, FTYPE *ctop, FTYPE *p_l, FTYPE *p_r, FTYPE *U_l, FTYPE *U_r,FTYPE *F_l,FTYPE *F_r, FTYPE *F);
742  int domustaflux;
743  FTYPE cmaxorig[NPR],cminorig[NPR],ctoporig[NPR];
744  FTYPE Uorig_l[NPR],Uorig_r[NPR],Forig_l[NPR],Forig_r[NPR],porig_l[NPR],porig_r[NPR];
745  FTYPE Forig[NPR];
746  FTYPE Fother[NPR];
747  int pl,pliter;
748  FTYPE shockstrength;
749  FTYPE rarestrength;
750  FTYPE contactstrength;
751  struct of_newtonstats newtonstats; setnewtonstatsdefault(&newtonstats);
752  int showmessages=1;
753  int allowlocalfailurefixandnoreport=1;
754 
755 
756 
757  if(DOEVOLVEUU){
758  // check for strong shock that MUSTA can't handle
759  shockstrength=fabs((p_r[UU]-p_l[UU])/max(min(fabs(p_l[UU]),fabs(p_r[UU])),SMALL));
760 
761  // contact strength
762  contactstrength=fabs((p_r[RHO]-p_l[RHO])/max(min(fabs(p_l[RHO]),fabs(p_r[RHO])),SMALL));
763  }
764  else{
765  shockstrength=contactstrength=0;
766  }
767 
768  // rarefaction strength
769  rarestrength=fabs((p_r[UU+dir]-p_l[UU+dir])/max(min(fabs(p_l[UU+dir]),fabs(p_r[UU+dir])),SMALL));
770 
771 
772  // keep original left/right states in case of musta failure
773  PLOOP(pliter,pl){
774 
775  cmaxorig[pl]=cmax[pl];
776  cminorig[pl]=cmin[pl];
777  ctoporig[pl]=ctop[pl];
778 
779  Uorig_l[pl]=U_l[pl];
780  Uorig_r[pl]=U_r[pl];
781  Forig_l[pl]=F_l[pl];
782  Forig_r[pl]=F_r[pl];
783  porig_l[pl]=p_l[pl];
784  porig_r[pl]=p_r[pl];
785  }
786 
787 
788 
789  if(DOMULTICELL==0){
790  domustaflux=musta1flux_compute(dir,geom,cmin_l,cmin_r,cmax_l,cmax_r,cmin,cmax,ctop,cforce,p_l,p_r,U_l,U_r,F_l,F_r,F);
791  }
792  else if(DOMULTICELL==1){
793  domustaflux=musta2flux_compute(dir,geom,cmin_l,cmin_r,cmax_l,cmax_r,cmin,cmax,ctop,cforce,p_l,p_r,U_l,U_r,F_l,F_r,F);
794  }
795 
797  //
798  // DONE WITH MUSTA. F is final musta solution or musta failed
799  //
801 
802  // check if musta failed during inversion
803  if(!domustaflux){
804  // reduce to HLL
805  hllflux_compute(dir,geom,cminorig,cmaxorig,ctoporig,porig_l,porig_r,Uorig_l,Uorig_r,Forig_l,Forig_r,F);
806 
807  // if(debugfail>=1) dualfprintf(fail_file,"DIDNOTDO_MUSTAFLUX: nstep=%ld t=%21.15g i=%d j=%d k=%d\n",nstep,t,geom->i,geom->j,geom->k);
808  }
809  else{
810  // if(debugfail>=1) dualfprintf(fail_file,"DIDDO_MUSTAFLUX: nstep=%ld t=%21.15g i=%d j=%d k=%d\n",nstep,t,geom->i,geom->j,geom->k);
811  }
812 
813 #if(0)
814  // PLOOP(pliter,pl) F[pl] = HLLCOMPUTE(cmin,cmax,U_l[pl],U_r[pl],F_l[pl],F_r[pl]);
815  hllflux_compute(dir,geom,cminorig,cmaxorig,ctoporig,porig_l,porig_r,Uorig_l,Uorig_r,Forig_l,Forig_r,Fother);
816  // average HLL and MUSTAHLL
817  PLOOP(pliter,pl) F[pl] = 0.5*(Fother[pl]+F[pl]);
818 #endif
819 
820 
821 
822 #if(HLLBOUNDMUSTA)
823  // don't allow for a change in the flux by more than 200%
824  hllflux_compute(dir,geom,cminorig,cmaxorig,ctoporig,porig_l,porig_r,Uorig_l,Uorig_r,Forig_l,Forig_r,Fother);
825  // PLOOP(pliter,pl) if( (Fother[pl]-F[pl])/(fabs(Fother[pl])+SMALL)>0.9) F[pl]=Fother[pl];
826  PLOOP(pliter,pl) if( fabs(Fother[pl]-F[pl])/(fabs(Fother[pl])+SMALL)>2.0) F[pl]=Fother[pl];
827 
828 #endif
829 
830 
831 
832  // check sign of flux and don't allow to switch
833 #if(0)
834  hllflux_compute(dir,geom,cminorig,cmaxorig,ctoporig,porig_l,porig_r,Uorig_l,Uorig_r,Forig_l,Forig_r,Fother);
835  PLOOP(pliter,pl) if( ((Fother[pl]<0.0)&&(F[pl]>0.0))||((Fother[pl]>0.0)&&(F[pl]<0.0)) ) F[pl]=Fother[pl];
836 #endif
837 
838 
839 
840 #if(0)
841  hllflux_compute(dir,geom,cminorig,cmaxorig,ctoporig,porig_l,porig_r,Uorig_l,Uorig_r,Forig_l,Forig_r,Fother);
842  //F[RHO]=Fother[RHO];
843  if( (shockstrength>0.5)||(rarestrength>0.5)){
844  F[UU]=Fother[UU];
845  F[U1]=Fother[U1];// leads to inconsistent rho*v with F[RHO]
846  F[U2]=Fother[U2];
847  F[U3]=Fother[U3];
848  }
849  // }
850  // F[B1]=Fother[B1];
851  // F[B2]=Fother[B2];
852  // F[B3]=Fother[B3];
853 
854 
855 #endif
856 
857 #if(1) // not too bad, but generates problem at moving contact.
858  hllflux_compute(dir,geom,cminorig,cmaxorig,ctoporig,porig_l,porig_r,Uorig_l,Uorig_r,Forig_l,Forig_r,Fother);
859 
860  F[UU]=Fother[UU];
861  F[U1]=Fother[U1];// leads to inconsistent rho*v with F[RHO]
862  F[U2]=Fother[U2];
863  F[U3]=Fother[U3];
864 #endif
865 
866 
867  return(0);
868 }
869 
870 
871 
873 int musta1flux_compute(int dir,struct of_geom *geom, FTYPE *cmin_l, FTYPE *cmin_r, FTYPE *cmax_l, FTYPE *cmax_r, FTYPE *cmin, FTYPE *cmax, FTYPE *ctop, FTYPE *cforce, FTYPE *p_l, FTYPE *p_r, FTYPE *U_l, FTYPE *U_r,FTYPE *F_l,FTYPE *F_r, FTYPE *F)
874 {
875  int forceflux_compute(int dir,struct of_geom *geom, FTYPE *cmin, FTYPE *cmax, FTYPE *ctop, FTYPE *cforce, FTYPE *p_l, FTYPE *p_r, FTYPE *U_l, FTYPE *U_r,FTYPE *F_l,FTYPE *F_r, FTYPE *F);
876  int hllflux_compute(int dir,struct of_geom *geom, FTYPE *cmin, FTYPE *cmax, FTYPE *ctop, FTYPE *p_l, FTYPE *p_r, FTYPE *U_l, FTYPE *U_r,FTYPE *F_l,FTYPE *F_r, FTYPE *F);
877  int pl,pliter;
878  FTYPE umid[NPR],plnew[NPR],prnew[NPR],fmid[NPR];
879  struct of_state state, state_l, state_r;
880  struct of_state *ptrstate, *ptrstate_l, *ptrstate_r;
881  int domustaflux;
882  int mustaloop;
883  FTYPE cmusta[NPR];
884  FTYPE cmineach_l[NUMEOMSETS],cmaxeach_l[NUMEOMSETS];
885  FTYPE cmineach_r[NUMEOMSETS],cmaxeach_r[NUMEOMSETS];
886  int ignorecourant;
887  int cminmax_calc(FTYPE cmin_l,FTYPE cmin_r,FTYPE cmax_l,FTYPE cmax_r,FTYPE *cmin,FTYPE *cmax,FTYPE *ctop);
888  FTYPE Fother[NPR];
889  FTYPE correctionl[NPR],correctionr[NPR],mymustacoef[NPR],mymustacoeffinal;
890  // FTYPE correctionl,correctionr,mymustacoef;
891  FTYPE Ucl,Ucr;
892  int pllargest;
893  FTYPE fracl[NPR],fracr[NPR];
894  FTYPE mustaf;
895  struct of_newtonstats newtonstats; setnewtonstatsdefault(&newtonstats);
896  FTYPE cminmhd,cmaxmhd,ctopmhd;
897  FTYPE cminrad,cmaxrad,ctoprad;
898  FTYPE cminrad2,cmaxrad2,ctoprad2;
899  int showmessages=1;
900  int allowlocalfailurefixandnoreport=1;
901 
902 
903 #if(EOMRADTYPE!=EOMRADNONE)
904  dualfprintf(fail_file,"musta1flux_compute not setup for radiation.");
905  myexit(1);
906 #endif
907 
908  // default
909  ptrstate = &state;
910  ptrstate_l = &state_l;
911  ptrstate_r = &state_r;
912 
913 
914  // default is to do musta unless musta fails
915  domustaflux=1;
916 
917 
918  // get predictor flux
919 #if(WHICHFLUX==MUSTAFORCE)
920  PLOOP(pliter,pl) cmusta[pl]=ctop[pl];
921  forceflux_compute(dir,geom,cmin,cmax,ctop,cforce,p_l,p_r,U_l,U_r,F_l,F_r,F);
922 #elif(WHICHFLUX==MUSTAHLL)
923  PLOOP(pliter,pl) cmusta[pl]=ctop[pl];
924  hllflux_compute(dir,geom,cmin,cmax,ctop,p_l,p_r,U_l,U_r,F_l,F_r,F);
925 #elif(WHICHFLUX==MUSTALAXF)
926  PLOOP(pliter,pl) cmusta[pl]=ctop[pl];
927  PLOOP(pliter,pl) F[pl] = LAXFCOMPUTE(ctop[pl],U_l[pl],U_r[pl],F_l[pl],F_r[pl]);
928 #endif
929 
930 
931 
932  // do musta stage loop
933  for(mustaloop=0;mustaloop<NUMMUSTAITERS;mustaloop++){
934 
935  PLOOP(pliter,pl){
936  // get wave speeds for opening Riemann fan
937 #if(WHICHFLUX==MUSTAFORCE)
938  // cmusta[pl]=cforce[pl];
939  cmusta[pl]=ctop[pl];
940 #elif(WHICHFLUX==MUSTAHLL)
941  cmusta[pl]=ctop[pl];
942 #elif(WHICHFLUX==MUSTALAXF)
943  cmusta[pl]=ctop[pl];
944 #endif
945  }
946 
947 
949  //
950  //
951 #define MUSTAVERSION 0
952 
953  // update U_l U_r (Open Riemann fan)
954 #if(MUSTAVERSION==0)
955  // does bad with Noh problem
956  PLOOP(pliter,pl){
957  U_l[pl] -= MUSTACOEF*(F[pl]-F_l[pl])/cmusta[pl];
958  U_r[pl] -= MUSTACOEF*(F_r[pl]-F[pl])/cmusta[pl];
959  }
960 #elif(MUSTAVERSION==1)
961  // does good with Noh problem
962 
963  PLOOP(pliter,pl){
964  correctionl[pl]=MUSTACOEF*(F[pl]-F_l[pl])/cmusta[pl];
965  correctionr[pl]=MUSTACOEF*(F_r[pl]-F[pl])/cmusta[pl];
966  // check for monotonicity
967  // if(fabs(U_l[pl]-U_r[pl])<fabs((U_l[pl]-correctionl[pl])-(U_r[pl]-correctionr[pl]))) return(0);
968  fracl[pl]=fabs(U_l[pl]-correctionl[pl])/(fabs(U_l[pl])+SMALL);
969  fracr[pl]=fabs(U_r[pl]-correctionr[pl])/(fabs(U_r[pl])+SMALL);
970 
971  mymustacoef[pl]=(U_l[pl]-U_r[pl])/(correctionl[pl]-correctionr[pl]);
972  // check that corrections are possibly tunable to be monotonic and so are at least in the right direction
973  }
974 
975  // find maximum of coefficients
976  mymustacoeffinal=0.0;
977  PLOOP(pliter,pl){
978  if(pl==RHO){
979  if(mymustacoeffinal<=mymustacoef[pl]){
980  mymustacoeffinal=mymustacoef[pl];
981  }
982  }
983  }
984  if(mymustacoeffinal>1.0) mymustacoeffinal=1.0;
985  if(mymustacoeffinal<0.9) mymustacoeffinal=0.9;
986 
987  // find maximum fractional change
988  // mustaf=0.0;
989  // pllargest=-1;
990  // PLOOP(pliter,pl){
991  // if(pl<=RHO){
992  // if(mustaf<=max(fracl[pl],fracr[pl])){
993  // mustaf=max(fracl[pl],fracr[pl]);
994  // pllargest=pl;
995  // }
996  // }
997  // }
998  // if((mustaf>1.0)||(mustaf<0.0)) return(0);
999  // MUSTA is now correctable
1000 
1001 
1002  // now use single musta coefficient
1003  PLOOP(pliter,pl){
1004  // mymustacoef=(mymustacoef>1.0) ? 1.0 : ((mymustacoef<0.0) ? 0.0 : mymustacoef);
1005  // mymustacoef=MUSTACOEF;
1006  U_l[pl] -= mymustacoeffinal*correctionl[pl];
1007  U_r[pl] -= mymustacoeffinal*correctionr[pl];
1008  }
1009 #elif(MUSTAVERSION==2)
1010  PLOOP(pliter,pl){
1011  correctionl[pl]=MUSTACOEF*(F[pl]-F_l[pl])/cmusta[pl];
1012  correctionr[pl]=MUSTACOEF*(F_r[pl]-F[pl])/cmusta[pl];
1013  // check for monotonicity
1014  // if(fabs(U_l[pl]-U_r[pl])<fabs((U_l[pl]-correctionl[pl])-(U_r[pl]-correctionr[pl]))) return(0);
1015 
1016  if(fabs(correctionl[pl]-correctionr[pl])!=0.0){
1017  mymustacoef[pl]=(U_l[pl]-U_r[pl])/(correctionl[pl]-correctionr[pl]);
1018  }
1019  else mymustacoef[pl]=BIG;
1020  // check that corrections are possibly tunable to be monotonic and so are at least in the right direction
1021  // if((mymustacoef[pl]>1.0)||(mymustacoef[pl]<0.0)) return(0);
1022 
1023  // if coef<0, then when using coef>0, this one will not be monotone
1024  // if(mymustacoef[pl]<0.0) return(0);
1025  // if(mymustacoef[pl]<-1.0) return(0);
1026  if(mymustacoef[pl]<0.0) return(0);
1027  }
1028  // MUSTA is now correctable
1029 
1030 #define MUSTACOEFTYPE 1
1031 
1032 #if(MUSTACOEFTYPE==0)
1033  if(1){
1034  PLOOP(pliter,pl){
1035  dualfprintf(fail_file,"nstep=%ld steppart=%d i=%d coef[%d]=%21.15g\n",nstep,steppart,geom->i,pl,mymustacoef[pl]);
1036  }
1037  }
1038  if(1){
1039  // find minimum of coefficients
1040  mymustacoeffinal=BIG;
1041  PLOOP(pliter,pl){
1042  if(mymustacoeffinal>mymustacoef[pl]) mymustacoeffinal=mymustacoef[pl];
1043  }
1044  if(mymustacoeffinal>1.0) return(0);
1045  // if(fabs(mymustacoef[RHO]-1.0)>2.0) return(0);
1046  // if(mymustacoeffinal>1.0) mymustacoeffinal=1.0;
1047  PLOOP(pliter,pl){
1048  mymustacoef[pl]=mymustacoeffinal;
1049  }
1050  }
1051 #elif(MUSTACOEFTYPE==1)
1052  // find maximum below 1
1053  mymustacoeffinal=0;
1054  PLOOP(pliter,pl){
1055  if( (mymustacoeffinal<mymustacoef[pl])&&(mymustacoef[pl]<=1.0)&&(mymustacoef[pl]>0.9) ) mymustacoeffinal=mymustacoef[pl];
1056  }
1057  // dualfprintf(fail_file,"nstep=%ld steppart=%d i=%d coef=%21.15g\n",nstep,steppart,geom->i,mymustacoeffinal);
1058  // dualfprintf(fail_file,"p_l[rho]=%21.15g p_r[rho]=%21.15g\n",p_l[RHO],p_r[RHO]);
1059  // if(mymustacoeffinal>1.0) return(0);
1060  // if(mymustacoeffinal>1.0) mymustacoeffinal=1.0;
1061  PLOOP(pliter,pl){
1062  mymustacoef[pl]=mymustacoeffinal;
1063  }
1064 #elif(MUSTACOEFTYPE==2)
1065  // find good coefficients (this doesn't make much sense)
1066  PLOOP(pliter,pl){
1067  if(mymustacoef[pl]<0.0) mymustacoef[pl]=0.0;
1068  if(mymustacoef[pl]>1.0) mymustacoef[pl]=1.0;
1069  }
1070 #elif(MUSTACOEFTYPE==3)
1071  // find good coefficients
1072  PLOOP(pliter,pl){
1073  mymustacoef[pl]=MUSTACOEF;
1074  }
1075 #endif
1076 
1077 
1078  // if(mymustacoeffinal<0.9) mymustacoeffinal=0.9;
1079  // mymustacoeffinal=1.0;
1080 
1081  // now use single musta coefficient
1082  PLOOP(pliter,pl){
1083  // mymustacoef=(mymustacoef>1.0) ? 1.0 : ((mymustacoef<0.0) ? 0.0 : mymustacoef);
1084  // mymustacoef=MUSTACOEF;
1085 
1086  // get corrected value
1087  Ucl=U_l[pl]-mymustacoef[pl]*correctionl[pl];
1088  Ucr=U_r[pl]-mymustacoef[pl]*correctionr[pl];
1089 
1090  // check for monotone result, can use final or pl coefficient
1091  // checking if final point is between original values
1092  if((U_r[pl] - U_l[pl])*(Ucl - U_l[pl])<0.0) return(0);
1093  if((U_r[pl] - U_l[pl])*(Ucr - U_r[pl])>0.0) return(0);
1094 
1095  // then good
1096  U_l[pl] = Ucl;
1097  U_r[pl] = Ucr;
1098  }
1099 #elif(MUSTAVERSION==3)
1100  PLOOP(pliter,pl){
1101  correctionl=(F[pl]-F_l[pl])/cmusta[pl];
1102  correctionr=(F_r[pl]-F[pl])/cmusta[pl];
1103  mymustacoef=MUSTACOEF;
1104 
1105  // fail
1106  U_l[pl] -= mymustacoef*correctionl;
1107  U_r[pl] -= mymustacoef*correctionr;
1108 
1109  // U_l[pl] -= MUSTACOEF*correctionl;
1110  // U_r[pl] -= MUSTACOEF*correctionr;
1111 
1112  // succeed
1113  //U_l[pl] -= MUSTACOEF*(F[pl]-F_l[pl])/cmusta;
1114  //U_r[pl] -= MUSTACOEF*(F_r[pl]-F[pl])/cmusta;
1115 
1116  // U_l[pl] -= mymustacoef*(F[pl]-F_l[pl])/cmusta[pl];
1117  // U_r[pl] -= mymustacoef*(F_r[pl]-F[pl])/cmusta[pl];
1118 
1119  }
1120 #endif
1121  //
1122  //
1124 
1125 
1126 
1127 
1128  // invert to get p_l p_r so can get F_l F_r and U_l U_r
1129  // get new primitive p_l
1130  int eomtype=EOMDEFAULT;
1131  FTYPE dissmeasure=-1.0; // assume energy try ok
1132  int whichcap=CAPTYPEBASIC;
1133  int whichmethod=MODEDEFAULT;
1134  int modprim=0;
1135  int checkoninversiongas=CHECKONINVERSION;
1136  int checkoninversionrad=CHECKONINVERSIONRAD;
1137  MYFUN(Utoprimgen(showmessages,checkoninversiongas,checkoninversionrad,allowlocalfailurefixandnoreport, 0, &eomtype,whichcap,whichmethod,modprim,EVOLVEUTOPRIM,UEVOLVE, U_l, ptrstate_l, geom, dissmeasure, p_l, p_l,&newtonstats),"flux.c:mustaflux_compute()", "Utoprimgen", 1);
1138  if(GLOBALMACP0A1(pflag,geom->i,geom->j,geom->k,FLAGUTOPRIMFAIL)){
1139  if(debugfail>=1) dualfprintf(fail_file,"Failed to find inversion for MUSTAFORCEFLUX(left): nstep=%ld t=%21.15g i=%d j=%d k=%d\n",nstep,t,geom->i,geom->j,geom->k);
1140  domustaflux=0;
1141  break;
1142  }
1143  // get new primitive p_r
1144  MYFUN(Utoprimgen(showmessages,checkoninversiongas,checkoninversionrad,allowlocalfailurefixandnoreport, 0, &eomtype,whichcap,whichmethod,modprim,EVOLVEUTOPRIM,UEVOLVE, U_r, ptrstate_r, geom, dissmeasure, p_r, p_r,&newtonstats),"flux.c:mustaflux_compute()", "Utoprimgen", 1);
1145  if(GLOBALMACP0A1(pflag,geom->i,geom->j,geom->k,FLAGUTOPRIMFAIL)){
1146  if(debugfail>=1) dualfprintf(fail_file,"Failed to find inversion for MUSTAFORCEFLUX(right): nstep=%ld t=%21.15g i=%d j=%d k=%d\n",nstep,t,geom->i,geom->j,geom->k);
1147  domustaflux=0;
1148  break;
1149  }
1150 
1151  // get fluxes (state used for vchar() below so need full state)
1152  MYFUN(get_state(p_l, geom, ptrstate_l),"flux.c:flux_compute()", "get_state()", 1);
1153  MYFUN(primtoflux(UEVOLVE, p_l, ptrstate_l, dir, geom, F_l, NULL),"flux.c:flux_compute()","primtoflux_calc() dir=1/2 l", 1);
1154 
1155  MYFUN(get_state(p_r, geom, ptrstate_r),"flux.c:flux_compute()", "get_state()", 1);
1156  MYFUN(primtoflux(UEVOLVE, p_r, ptrstate_r, dir, geom, F_r, NULL),"flux.c:flux_compute()","primtoflux_calc() dir=1/2 l", 1);
1157 
1158  // now have p_l, p_r, U_l, U_r, F_l, F_r
1159 
1160 
1161 #if(1||((WHICHFLUX==MUSTAHLL)||(WHICHFLUX==MUSTALAXF)))
1162  // wave speeds for new left-right states
1163  MYFUN(vchar_each(p_l, ptrstate_l, dir, geom, &cmaxeach_l[EOMSETMHD], &cmineach_l[EOMSETMHD], &cmaxeach_l[EOMSETRAD], &cmineach_l[EOMSETRAD], &cmaxeach_l[EOMSETRADFORDT], &cmineach_l[EOMSETRADFORDT],&ignorecourant),"step_ch.c:fluxcalc()", "vchar() dir=1or2", 1);
1164  MYFUN(vchar_each(p_r, ptrstate_r, dir, geom, &cmaxeach_r[EOMSETMHD], &cmineach_r[EOMSETMHD], &cmaxeach_r[EOMSETRAD], &cmineach_r[EOMSETRAD], &cmaxeach_r[EOMSETRADFORDT], &cmineach_r[EOMSETRADFORDT],&ignorecourant),"step_ch.c:fluxcalc()", "vchar() dir=1or2", 2);
1165  cminmax_calc(cmineach_l[EOMSETMHD],cmineach_r[EOMSETMHD],cmaxeach_l[EOMSETMHD],cmaxeach_r[EOMSETMHD],&cminmhd,&cmaxmhd,&ctopmhd);
1166  cminmax_calc(cmineach_l[EOMSETRAD],cmineach_r[EOMSETRAD],cmaxeach_l[EOMSETRAD],cmaxeach_r[EOMSETRAD],&cminrad,&cmaxrad,&ctoprad);
1167  cminmax_calc(cmineach_l[EOMSETRADFORDT],cmineach_r[EOMSETRADFORDT],cmaxeach_l[EOMSETRADFORDT],cmaxeach_r[EOMSETRADFORDT],&cminrad2,&cmaxrad2,&ctoprad2);
1168  // cmin, cmax, ctop for fluxes
1169  PLOOP(pliter,pl){
1170  if(pl<URAD0 || pl>URAD3){
1171  cmin[pl]=cminmhd;
1172  cmax[pl]=cmaxmhd;
1173  ctop[pl]=ctopmhd;
1174  }
1175  else{
1176  cmin[pl]=cminrad;
1177  cmax[pl]=cmaxrad;
1178  ctop[pl]=ctoprad;
1179  }
1180  }
1181 #endif
1182 
1183 
1184  // get F
1185  // get flux using corrected left/right states
1186 #if(WHICHFLUX==MUSTAFORCE)
1187  PLOOP(pliter,pl) cforce[pl]=ctop[pl];
1188  forceflux_compute(dir,geom,cmin,cmax,ctop,cforce,p_l,p_r,U_l,U_r,F_l,F_r,F);
1189 #elif(WHICHFLUX==MUSTAHLL)
1190  hllflux_compute(dir,geom,cmin,cmax,ctop,p_l,p_r,U_l,U_r,F_l,F_r,F);
1191 #elif(WHICHFLUX==MUSTALAXF)
1192  PLOOP(pliter,pl) F[pl] = LAXFCOMPUTE(ctop[pl],U_l[pl],U_r[pl],F_l[pl],F_r[pl]);
1193 #endif
1194 
1195  }// end musta loop
1196 
1197 
1198  return(domustaflux);
1199 }
1200 
1201 
1202 
1203 
1204 
1206 int musta2flux_compute(int dir,struct of_geom *geom, FTYPE *cmin_l, FTYPE *cmin_r, FTYPE *cmax_l, FTYPE *cmax_r, FTYPE *cmin, FTYPE *cmax, FTYPE *ctop, FTYPE *cforce, FTYPE *p_l, FTYPE *p_r, FTYPE *U_l, FTYPE *U_r,FTYPE *F_l,FTYPE *F_r, FTYPE *F)
1207 {
1208  int forceflux_compute(int dir,struct of_geom *geom, FTYPE *cmin, FTYPE *cmax, FTYPE *ctop, FTYPE *cforce, FTYPE *p_l, FTYPE *p_r, FTYPE *U_l, FTYPE *U_r,FTYPE *F_l,FTYPE *F_r, FTYPE *F);
1209  int hllflux_compute(int dir,struct of_geom *geom, FTYPE *cmin, FTYPE *cmax, FTYPE *ctop, FTYPE *p_l, FTYPE *p_r, FTYPE *U_l, FTYPE *U_r,FTYPE *F_l,FTYPE *F_r, FTYPE *F);
1210  int pl,pliter;
1211  FTYPE umid[NPR],plnew[NPR],prnew[NPR],fmid[NPR];
1212  struct of_state state;
1213  struct of_state *ptrstate;
1214  int domustaflux;
1215  int mustaloop;
1216  FTYPE cmusta[NPR];
1217  FTYPE *Unow,*Fnow,*pnow;
1218  FTYPE *Unow_l, *Unow_r, *Fnow_l, *Fnow_r, *pnow_l, *pnow_r;
1219  FTYPE *cmaxnow,*cminnow,*ctopnow;
1220  FTYPE Umusta[2*(NUMLOCALCELLS+1)][NPR],Fmusta[2*(NUMLOCALCELLS+1)][NPR],Fmusta_edge[2*(NUMLOCALCELLS+1)][NPR],pmusta[2*(NUMLOCALCELLS+1)][NPR];
1221  FTYPE cmaxmusta[2*(NUMLOCALCELLS+1)][NPR],cminmusta[2*(NUMLOCALCELLS+1)][NPR];
1222  FTYPE cmaxnow_l[NPR], cmaxnow_r[NPR], cminnow_l[NPR], cminnow_r[NPR];
1223  int ignorecourant;
1224  int cminmax_calc(FTYPE cmin_l,FTYPE cmin_r,FTYPE cmax_l,FTYPE cmax_r,FTYPE *cmin,FTYPE *cmax,FTYPE *ctop);
1225  FTYPE otherF[NPR];
1226  int mustacellloop;
1227  int ms,me;
1228  struct of_newtonstats newtonstats; setnewtonstatsdefault(&newtonstats);
1229  int showmessages=1;
1230  int allowlocalfailurefixandnoreport=1;
1231 
1232 
1233  // setup multi-cell MUSTA
1234  // left state
1235  for(mustacellloop=0;mustacellloop<=NUMLOCALCELLS;mustacellloop++){
1236  PLOOP(pliter,pl){
1237  Umusta[mustacellloop][pl]=U_l[pl];
1238  Fmusta[mustacellloop][pl]=F_l[pl];
1239  pmusta[mustacellloop][pl]=p_l[pl];
1240  cmaxmusta[mustacellloop][pl]=cmax_l[pl];
1241  cminmusta[mustacellloop][pl]=cmin_l[pl];
1242  }
1243  }
1244  // right state
1245  for(mustacellloop=NUMLOCALCELLS+1;mustacellloop<=2*NUMLOCALCELLS+1;mustacellloop++){
1246  PLOOP(pliter,pl){
1247  Umusta[mustacellloop][pl]=U_r[pl];
1248  Fmusta[mustacellloop][pl]=F_r[pl];
1249  pmusta[mustacellloop][pl]=p_r[pl];
1250  cmaxmusta[mustacellloop][pl]=cmax_r[pl];
1251  cminmusta[mustacellloop][pl]=cmin_r[pl];
1252  }
1253  }
1254 
1255 
1256  // default is to do musta unless musta fails
1257  domustaflux=1;
1258 
1259 
1260 
1262  //
1263  // loop over musta stages
1264  //
1265  // final flux is located on edge, while left/right states are "cell centered"
1266  //
1268  for(mustaloop=0;mustaloop<NUMMULTIMUSTAITERS;mustaloop++){
1269 
1270  /*
1271  // range of fluxes required
1272  if(mustaloop<(NUMMULTIMUSTAITERS-1)/2){
1273  ms=NUMLOCALCELLS+1-mustaloop;
1274  me=NUMLOCALCELLS+2+mustaloop;
1275  }
1276  else if(mustaloop>=(NUMMULTIMUSTAITERS-1)/2){
1277  ms=mustaloop+2;
1278  me=NUMLOCALCELLS+1-mustaloop;
1279  }
1280  */
1281 
1283  //
1284  // get predictor edge flux
1285  //
1287  for(mustacellloop=1;mustacellloop<=2*NUMLOCALCELLS+1;mustacellloop++){
1288  // for(mustacellloop=ms;mustacellloop<=me;mustacellloop++){
1289 
1290  // pointers to memory locations GETS RETRIEVED
1291  Unow_l=Umusta[mustacellloop-1];
1292  Unow_r=Umusta[mustacellloop];
1293  Fnow_l=Fmusta[mustacellloop-1];
1294  Fnow_r=Fmusta[mustacellloop];
1295  pnow_l=pmusta[mustacellloop-1];
1296  pnow_r=pmusta[mustacellloop];
1297  // assignments GETS RETRIEVED
1298  PLOOP(pliter,pl){
1299  cmaxnow_l[pl]=cmaxmusta[mustacellloop-1][pl];
1300  cmaxnow_r[pl]=cmaxmusta[mustacellloop][pl];
1301  cminnow_l[pl]=cminmusta[mustacellloop-1][pl];
1302  cminnow_r[pl]=cminmusta[mustacellloop][pl];
1303  }
1304  // Fnow on edge GETS SET
1305  Fnow=Fmusta_edge[mustacellloop];
1306 
1307 
1308  // get wave speeds and min/max/top versions of left/right states GETS RETRIEVED
1309  // cminmax_calc(cminnow_l,cminnow_r,cmaxnow_l,cmaxnow_r,&cmin,&cmax,&ctop);
1310  // SUPERGODMARK: Not setup for radiation
1311  if(EOMRADTYPE!=EOMRADNONE){
1312  dualfprintf(fail_file,"musta2 not setup for radiation\n");
1313  myexit(1);
1314  }
1315 #if(WHICHFLUX==MUSTAFORCE)
1316  PLOOP(pliter,pl) cmusta[pl]=cforce[pl];
1317  forceflux_compute(dir,geom,cmin,cmax,ctop,cforce,pnow_l,pnow_r,Unow_l,Unow_r,Fnow_l,Fnow_r,Fnow);
1318 #elif(WHICHFLUX==MUSTAHLL)
1319  PLOOP(pliter,pl) cmusta[pl]=ctop[pl];
1320  hllflux_compute(dir,geom,cmin,cmax,ctop,pnow_l,pnow_r,Unow_l,Unow_r,Fnow_l,Fnow_r,Fnow);
1321 #elif(WHICHFLUX==MUSTALAXF)
1322  PLOOP(pliter,pl) cmusta[pl]=ctop[pl];
1323  PLOOP(pliter,pl) Fnow[pl] = LAXFCOMPUTE(ctop[pl],Unow_l[pl],Unow_r[pl],Fnow_l[pl],Fnow_r[pl]);
1324 #endif
1325 
1326 
1327 
1328 
1329 
1330 
1331  }// end loop over musta cells
1332 
1333  // Fnow in middle is final flux (NUMMULTIMUSTAITERS==1 gives same as no MUSTA)
1334  // (NUMMULTIMUSTAITERS==2 gives MUSTA-1)
1335  if(mustaloop==NUMMULTIMUSTAITERS-1) break;
1336 
1338  //
1339  // update U_l U_r (Open Riemann fan)
1340  //
1342  for(mustacellloop=1;mustacellloop<=2*NUMLOCALCELLS;mustacellloop++){
1343 
1344 
1345  // Unow, pnow, and Fnow at center GETS SET
1346  Unow=Umusta[mustacellloop];
1347  pnow=pmusta[mustacellloop];
1348  Fnow=Fmusta[mustacellloop];
1349  PLOOP(pliter,pl){
1350  cmaxnow[pl]=cmaxmusta[mustacellloop][pl];
1351  cminnow[pl]=cminmusta[mustacellloop][pl];
1352  }
1353  // Fnow_l and Fnow_r on edge GETS RETRIEVED
1354  Fnow_l=Fmusta_edge[mustacellloop];
1355  Fnow_r=Fmusta_edge[mustacellloop+1];
1356 
1357 
1358 
1359  // set wave speed
1360  PLOOP(pliter,pl){
1361 #if(WHICHFLUX==MUSTAFORCE)
1362  cmusta[pl]=cforce[pl];
1363 #elif(WHICHFLUX==MUSTAHLL)
1364  cmusta[pl]=ctop[pl];
1365 #elif(WHICHFLUX==MUSTALAXF)
1366  cmusta[pl]=ctop[pl];
1367 #endif
1368  }
1369 
1370  // Open Riemann fan
1371  PLOOP(pliter,pl) Unow[pl] -= MULTIMUSTACOEF*(Fnow_r[pl]-Fnow_l[pl])/cmusta[pl];
1372 
1373 
1374  // get new primitive pnow from Unow
1375  int eomtype=EOMDEFAULT;
1376  FTYPE dissmeasure=-1.0; // assume energy try ok
1377  int whichcap=CAPTYPEBASIC;
1378  int whichmethod=MODEDEFAULT;
1379  int modprim=0;
1380  int checkoninversiongas=CHECKONINVERSION;
1381  int checkoninversionrad=CHECKONINVERSIONRAD;
1382  MYFUN(Utoprimgen(showmessages,checkoninversiongas,checkoninversionrad,allowlocalfailurefixandnoreport, 0, &eomtype,whichcap,whichmethod,modprim,EVOLVEUTOPRIM,UEVOLVE, Unow, NULL, geom, dissmeasure, pnow, pnow,&newtonstats),"flux.c:mustaflux_compute()", "Utoprimgen", 1);
1383  if(GLOBALMACP0A1(pflag,geom->i,geom->j,geom->k,FLAGUTOPRIMFAIL)){
1384  if(debugfail>=1) dualfprintf(fail_file,"Failed to find inversion for MUSTAFORCEFLUX(right): nstep=%ld t=%21.15g i=%d j=%d k=%d\n",nstep,t,geom->i,geom->j,geom->k);
1385  domustaflux=0;
1386  break;
1387  }
1388 
1389 
1390  // get corrected centered fluxes (state needed for vchar() so need full state)
1391  ptrstate=&state; // default
1392  MYFUN(get_state(pnow, geom, ptrstate),"flux.c:flux_compute()", "get_state()", 1);
1393  MYFUN(primtoflux(UEVOLVE, pnow, ptrstate, dir, geom, Fnow, NULL),"flux.c:flux_compute()","primtoflux_calc() dir=1/2 l", 1);
1394  // now have p, U, F at center again
1395 
1396  // get min and max wave speeds
1397  MYFUN(vchar(pnow, ptrstate, dir, geom, cmaxnow, cminnow, &ignorecourant),"step_ch.c:fluxcalc()", "vchar() dir=1or2", 2);
1398  // SUPERGODMARK: Not setup for radiation
1399 
1400  }// end loop over musta cells
1401 
1402 
1403  if(domustaflux==0) break;
1404 
1406  //
1407  // Set boundary conditions (zeroth order copy -- outflow)
1408  //
1410 
1411  PLOOP(pliter,pl){
1412  Umusta[0][pl]=Umusta[1][pl];
1413  Fmusta[0][pl]=Fmusta[1][pl];
1414  pmusta[0][pl]=pmusta[1][pl];
1415 
1416  Umusta[2*NUMLOCALCELLS+1][pl]=Umusta[2*NUMLOCALCELLS][pl];
1417  Fmusta[2*NUMLOCALCELLS+1][pl]=Fmusta[2*NUMLOCALCELLS][pl];
1418  pmusta[2*NUMLOCALCELLS+1][pl]=pmusta[2*NUMLOCALCELLS][pl];
1419  }
1420 
1421  }// end loop over musta stages
1422 
1423 
1425  //
1426  // assign final edge flux as MUSTA solution
1427  //
1429 
1430  PLOOP(pliter,pl) F[pl]=Fmusta_edge[NUMLOCALCELLS+1][pl];
1431 
1432 
1433 
1434  return(domustaflux);
1435 }
1436 
1437 
1438 
1440 void choose_flux(int fluxmethodlocal, int i, int j, int k, int pl, FTYPE *laxffrac,FTYPE *hllfrac)
1441 {
1442 #if(FLUXADJUST)
1443 
1444 #if(HYDROFLUXADJUSTONLY)
1445  if(pl<B1){
1446  if(GLOBALMACP0A1(pflag,i,j,k,FLAGREALFLUX)==HLLFLUX){ hllfrac[pl]=1.0; laxffrac[pl]=0.0; }
1447  else { hllfrac[pl]=0.0; laxffrac[pl]=1.0; }
1448  }
1449  else{
1450  if(fluxmethodlocal==HLLFLUX){
1451  hllfrac[pl]=1.0;
1452  laxffrac[pl]=0.0;
1453  }
1454  else if(fluxmethodlocal==LAXFFLUX){
1455  hllfrac[pl]=0.0;
1456  laxffrac[pl]=1.0;
1457  }
1458  }
1459 #else
1460  if(GLOBALMACP0A1(pflag,i,j,k,FLAGREALFLUX)==HLLFLUX){ hllfrac[pl]=1.0; laxffrac[pl]=0.0; }
1461  else { hllfrac[pl]=0.0; laxffrac[pl]=1.0; }
1462 #endif
1463 
1464 #else
1465  if(fluxmethodlocal==HLLFLUX){
1466  hllfrac[pl]=1.0;
1467  laxffrac[pl]=0.0;
1468  }
1469  else if(fluxmethodlocal==LAXFFLUX){
1470  hllfrac[pl]=0.0;
1471  laxffrac[pl]=1.0;
1472  }
1473 #endif
1474 
1475 }
1476 
Generated on Fri May 20 2016 15:52:33 for HARM by   doxygen 1.8.3.1