HARM
harm and utilities
 All Data Structures Files Functions Variables Typedefs Macros Pages
interpline.c
Go to the documentation of this file.
1 
2 
3 
10 #include "decs.h"
11 
12 
13 
14 
15 
16 
17 
18 #if(MERGEDC2EA2CMETHOD)
19 #define LINETYPEDEFINESOTHER FTYPE a_youtpolycoef[NPR2INTERP][MAXSPACEORDER][NBIGM];
20 #error "SUPERGODMARK: Should only need this if doing more than 3 point stencil -- fix"
21 #else
22 #define LINETYPEDEFINESOTHER FTYPE a_youtpolycoef[1][1][1]; //VS complained about defining zero-size object
23 #endif
24 
25 
27 #define LINETYPEDEFINES1 \
28  int bs,be,ps,pe; \
29  int di,dj,dk; \
30  int is,ie,js,je,ks,ke; \
31  int preforder, whichreduce; \
32  int pl,pliter; \
33  int recontype; \
34  extern void pass_1d_line_weno(int whichquantity, int dir, int do_weight_or_recon, int recontype, int whichreduce, int preforder, int pl, int bs, int ps, int pe, int be, int *minorder, int *maxorder, int *shift, FTYPE (*shockindicator)[NBIGM], FTYPE *stiffindicator, FTYPE (*Vline)[NBIGM], FTYPE (*Pline)[NBIGM], FTYPE (*df)[NBIGM], FTYPE (*dP)[NBIGM], FTYPE (*etai)[NBIGM], FTYPE (*monoindicator)[NBIGM], FTYPE (*yprim)[2][NBIGM], FTYPE (*ystencilvar)[NBIGM], FTYPE (*yin)[NBIGM], FTYPE (*yout)[NBIGM], FTYPE (*youtpolycoef)[NBIGM], struct of_trueijkp *trueijkp); \
35  extern void pass_1d_line_paraline(int whichquantity, int dir, int do_weight_or_recon, int recontype, int whichreduce, int preforder, int pl, int bs, int ps, int pe, int be, int *minorder, int *maxorder, int *shift, FTYPE (*shockindicator)[NBIGM], FTYPE *stiffindicator, FTYPE (*Vline)[NBIGM], FTYPE (*Pline)[NBIGM], FTYPE (*df)[NBIGM], FTYPE (*dP)[NBIGM], FTYPE (*etai)[NBIGM], FTYPE (*monoindicator)[NBIGM], FTYPE (*yprim)[2][NBIGM], FTYPE (*ystencilvar)[NBIGM], FTYPE (*yin)[NBIGM], FTYPE (*yout)[NBIGM], FTYPE (*youtpolycoef)[NBIGM], struct of_trueijkp *trueijkp); \
36  extern void pass_1d_line_multipl_weno(int MULTIPLTYPE, int whichquantity, int dir, int do_weight_or_recon, int recontype, int whichreduce, int preforder, int bs, int ps, int pe, int be, int *minorder, int *maxorder, int *shift, FTYPE (*shockindicator)[NBIGM], FTYPE *stiffindicator, FTYPE (*Vline)[NBIGM], FTYPE (*Pline)[NBIGM], FTYPE (*df)[NUMDFS][NBIGM], FTYPE (*dP)[NBIGM], FTYPE (*etai)[NUMTRUEEOMSETS][NBIGM], FTYPE (*monoindicator)[NUMMONOINDICATORS][NBIGM], FTYPE (*yprim)[2][NBIGM], FTYPE (*ystecilvar)[2][NBIGM], FTYPE (*yin)[2][NBIGM], FTYPE (*yout)[2][NBIGM], FTYPE (*youtpolycoef)[MAXSPACEORDER][NBIGM], struct of_trueijkp *trueijkp); \
37  extern void pass_1d_line_multipl_paraline(int MULTIPLTYPE, int whichquantity, int dir, int do_weight_or_recon, int recontype, int whichreduce, int preforder, int bs, int ps, int pe, int be, int *minorder, int *maxorder, int *shift, FTYPE (*shockindicator)[NBIGM], FTYPE *stiffindicator, FTYPE (*Vline)[NBIGM], FTYPE (*Pline)[NBIGM], FTYPE (*df)[NUMDFS][NBIGM], FTYPE (*dP)[NBIGM], FTYPE (*etai)[NUMTRUEEOMSETS][NBIGM], FTYPE (*monoindicator)[NUMMONOINDICATORS][NBIGM], FTYPE (*yprim)[2][NBIGM], FTYPE (*ystecilvar)[2][NBIGM], FTYPE (*yin)[2][NBIGM], FTYPE (*yout)[2][NBIGM], FTYPE (*youtpolycoef)[MAXSPACEORDER][NBIGM], struct of_trueijkp *trueijkp); \
38  int intdir; \
39  int withshifts; \
40  int nprlocalstart,nprlocalend; \
41  int nprlocallist[MAXNPR]; \
42  int pllocal; \
43  int numprims; \
44  int stencilvarisnull; \
45  int doingweno; // last thing has no semicolon
46 
47 
48 
49 #define GEN3MAC(numypl,pl,numywhich,which,numyi,i) (i + numyi*which + numyi*numywhich*pl)
50 
52 #define GEN2MAC(numypl,pl,numyi,i) (i + numyi*pl)
53 
54 #define GEN1MAC(numyi,i) (i)
55 
57 #define YINOUTMAC(numypl,pl,numywhich,which,numyi,i) GEN3MAC(numypl,pl,numywhich,which,numyi,i)
58 
59 
60 
64 #define LINETYPEDEFINESMEMORY \
65  FTYPE a_yin[NPR2INTERP][2][NBIGM]; \
66  FTYPE a_yout[NPR2INTERP][2][NBIGM]; \
67  \
68  FTYPE a_shockindicator[NUMTRUEEOMSETS][NBIGM]; \
69  FTYPE a_stiffindicator[NBIGM]; \
70  \
71  FTYPE a_yprim[NPR2INTERP][2][NBIGM]; \
72  FTYPE a_ystencilvar[NPR2INTERP][2][NBIGM]; \
73  FTYPE a_df[NPR2INTERP][NUMDFS][NBIGM]; \
74  FTYPE a_dfformono[NPR2INTERP][NUMDFS][NBIGM]; \
75  FTYPE a_drho[NUMDFS][NBIGM]; \
76  FTYPE a_dP[NUMDFS][NBIGM]; \
77  FTYPE a_etai[NPR2INTERP][NUMTRUEEOMSETS][NBIGM]; \
78  FTYPE a_Pline[NUMTRUEEOMSETS][NBIGM]; \
79  FTYPE a_Vline[NUMTRUEEOMSETS][NBIGM]; \
80  int a_shift[NBIGM]; \
81  FTYPE a_monoindicator[NPR2INTERP][NUMMONOINDICATORS][NBIGM]; \
82  int a_minorder[NBIGM]; \
83  int a_maxorder[NBIGM];
84 
85 
86 
89 #define LINETYPEDEFINEPOINTERS \
90  FTYPE (*yin)[2][NBIGM]; \
91  FTYPE (*yout)[2][NBIGM]; \
92  \
93  FTYPE (*shockindicator)[NBIGM]; \
94  FTYPE *stiffindicator; \
95  \
96  FTYPE (*yprim)[2][NBIGM]; \
97  FTYPE (*ystencilvar)[2][NBIGM]; \
98  FTYPE (*df)[NUMDFS][NBIGM]; \
99  FTYPE (*dfformono)[NUMDFS][NBIGM]; \
100  FTYPE (*drho)[NBIGM]; \
101  FTYPE (*dP)[NBIGM]; \
102  FTYPE (*etai)[NUMTRUEEOMSETS][NBIGM]; \
103  FTYPE (*Pline)[NBIGM]; \
104  FTYPE (*Vline)[NBIGM]; \
105  int *shift; \
106  FTYPE (*monoindicator)[NUMMONOINDICATORS][NBIGM]; \
107  int *minorder; \
108  int *maxorder; \
109  FTYPE (*youtpolycoef)[MAXSPACEORDER][NBIGM];
110 
111 
112 
113 
114 
115 #define LINETYPEDEFINES2 \
116  int reallim; \
117  int pl,pliter; \
118  int yiter; \
119  int dfiter; \
120  int counter; \
121  int pllocal; \
122  struct of_trueijkp trueijkp; \
123  extern int choose_limiter(int dir, int i, int j, int k, int pl); \
124  extern void compute_monotonicity_line(int recontype, int whichreduce, int preforder, int pl, int bs, int ps, int pe, int be, int *minorder, int *maxorder, int *shift, FTYPE (*shockindicator)[NBIGM], FTYPE (*df)[NBIGM], FTYPE (*monoindicator)[NBIGM] , FTYPE *yin, FTYPE (*yout)[NBIGM], FTYPE (*youtpolycoef)[NBIGM]); \
125  extern void compute_monotonicity_line_multipl(int stencilvarisnull, int MULTIPLTYPE, int whichquantity, int dir, int do_weight_or_recon, int recontype, int whichreduce, int preforder, int bs, int ps, int pe, int be, int *minorder, int *maxorder, int *shift, FTYPE (*shockindicator)[NBIGM], FTYPE *stiffindicator, FTYPE (*Vline)[NBIGM], FTYPE (*Pline)[NBIGM], FTYPE (*df)[NUMDFS][NBIGM], FTYPE (*dP)[NBIGM], FTYPE (*etai)[NUMTRUEEOMSETS][NBIGM], FTYPE (*monoindicator)[NUMMONOINDICATORS][NBIGM], FTYPE (*yprim)[2][NBIGM], FTYPE (*ystencilvar)[2][NBIGM], FTYPE (*yin)[2][NBIGM], FTYPE (*yout)[2][NBIGM], FTYPE (*youtpolycoef)[MAXSPACEORDER][NBIGM]); \
126  extern void compute_monotonicity_line_indicatoronly(int recontype, int whichreduce, int preforder, int pl, int bs, int ps, int pe, int be, int *minorder, int *maxorder, int *shift, FTYPE (*shockindicator)[NBIGM], FTYPE (*df)[NBIGM], FTYPE (*monoindicator)[NBIGM] , FTYPE *yin, FTYPE (*yout)[NBIGM], FTYPE (*youtpolycoef)[NBIGM]); \
127  extern void compute_monotonicity_line_valueonly(int recontype, int whichreduce, int preforder, int pl, int bs, int ps, int pe, int be, int *minorder, int *maxorder, int *shift, FTYPE (*shockindicator)[NBIGM], FTYPE (*df)[NBIGM], FTYPE (*monoindicator)[NBIGM] , FTYPE *yin, FTYPE (*yout)[NBIGM], FTYPE (*youtpolycoef)[NBIGM])
128 
129 // last thing has no semicolon
130 
131 
132 
133 
135 #define LINETYPESHIFTS \
136  {yin =(FTYPE (*)[2][NBIGM]) (&(a_yin[0][0][NBIGBND])); \
137  yout=(FTYPE (*)[2][NBIGM]) (&(a_yout[0][0][NBIGBND])); \
138  \
139  shockindicator=(FTYPE (*)[NBIGM]) (&(a_shockindicator[0][NBIGBND])); \
140  stiffindicator=(FTYPE (*)) (&(a_stiffindicator[NBIGBND])); \
141  \
142  yprim =(FTYPE (*)[2][NBIGM]) (&(a_yprim[0][0][NBIGBND])); \
143  ystencilvar =(FTYPE (*)[2][NBIGM]) (&(a_ystencilvar[0][0][NBIGBND])); \
144  df=(FTYPE (*)[NUMDFS][NBIGM]) (&(a_df[0][0][NBIGBND])); \
145  dfformono=(FTYPE (*)[NUMDFS][NBIGM]) (&(a_dfformono[0][0][NBIGBND])); \
146  drho=(FTYPE (*)[NBIGM]) (&(a_drho[0][NBIGBND])); \
147  dP=(FTYPE (*)[NBIGM]) (&(a_dP[0][NBIGBND])); \
148  etai=(FTYPE (*)[NUMTRUEEOMSETS][NBIGM]) (&(a_etai[0][0][NBIGBND])); \
149  Pline=(FTYPE (*)[NBIGM]) (&(a_Pline[0][NBIGBND])); \
150  Vline=(FTYPE (*)[NBIGM]) (&(a_Vline[0][NBIGBND])); \
151  shift=(int (*)) (&(a_shift[NBIGBND])); \
152  monoindicator=(FTYPE (*)[NUMMONOINDICATORS][NBIGM]) (&(a_monoindicator[0][0][NBIGBND])); \
153  minorder=(int (*)) (&(a_minorder[NBIGBND])); \
154  maxorder=(int (*)) (&(a_maxorder[NBIGBND])); \
155  youtpolycoef=(FTYPE (*)[MAXSPACEORDER][NBIGM]) (&(a_youtpolycoef[0][0][NBIGBND]));}
156 
157 
158 
159 #define LINETYPEPOINTERSFUNCTIONDECLARE \
160  FTYPE (*yin)[2][NBIGM], \
161  FTYPE (*yout)[2][NBIGM], \
162  \
163  FTYPE (*shockindicator)[NBIGM], \
164  FTYPE *stiffindicator, \
165  \
166  FTYPE (*yprim)[2][NBIGM], \
167  FTYPE (*ystencilvar)[2][NBIGM], \
168  FTYPE (*df)[NUMDFS][NBIGM], \
169  FTYPE (*dfformono)[NUMDFS][NBIGM], \
170  FTYPE (*drho)[NBIGM], \
171  FTYPE (*dP)[NBIGM], \
172  FTYPE (*etai)[NUMTRUEEOMSETS][NBIGM], \
173  FTYPE (*Pline)[NBIGM], \
174  FTYPE (*Vline)[NBIGM], \
175  int *shift, \
176  FTYPE (*monoindicator)[NUMMONOINDICATORS][NBIGM], \
177  int *minorder, \
178  int *maxorder, \
179  FTYPE (*youtpolycoef)[MAXSPACEORDER][NBIGM]
180 // no comma on end
181 
182 
183 
184 #define LINETYPEPOINTERSPASSFUNCTIONARG \
185  yin, \
186  yout, \
187  \
188  shockindicator, \
189  stiffindicator, \
190  \
191  yprim , \
192  ystencilvar, \
193  df, \
194  dfformono, \
195  drho, \
196  dP, \
197  etai, \
198  Pline, \
199  Vline, \
200  shift, \
201  monoindicator, \
202  minorder, \
203  maxorder, \
204  youtpolycoef
205 // no comma on end
206 
207 
208 
209 
210 
211 
212 
213 static void slope_lim_linetype_c2e_perijk_precomputeindicators(
214  LINETYPEPOINTERSFUNCTIONDECLARE,
215  int i, int j, int k,
216  int realisinterp, int whichprimtype, int interporflux,
217  int nprlocalstart, int nprlocalend, int *nprlocallist, int numprims,
218  int dir, int intdir,
219  int is, int ie, int js, int je, int ks, int ke, int di, int dj, int dk, int bs, int ps, int pe, int be,
220  int recontype, int doingweno,
221  void (*pass_1d_line)(int whichquantity, int dir, int do_weight_or_recon, int recontype, int whichreduce, int preforder, int pl, int bs, int ps, int pe, int be, int *minorder, int *maxorder, int *shift, FTYPE (*shockindicator)[NBIGM], FTYPE *stiffindicator, FTYPE (*Vline)[NBIGM], FTYPE (*Pline)[NBIGM], FTYPE (*df)[NBIGM], FTYPE (*dP)[NBIGM], FTYPE (*etai)[NBIGM], FTYPE (*monoindicator)[NBIGM], FTYPE (*yprim)[2][NBIGM], FTYPE (*ystencilvar)[NBIGM], FTYPE (*yin)[NBIGM], FTYPE (*yout)[NBIGM], FTYPE (*youtpolycoef)[NBIGM], struct of_trueijkp *trueijkp),
222  void (*pass_1d_line_multipl)(int MULTIPLTYPE, int whichquantity, int dir, int do_weight_or_recon, int recontype, int whichreduce, int preforder, int bs, int ps, int pe, int be, int *minorder, int *maxorder, int *shift, FTYPE (*shockindicator)[NBIGM], FTYPE *stiffindicator, FTYPE (*Vline)[NBIGM], FTYPE (*Pline)[NBIGM], FTYPE (*df)[NUMDFS][NBIGM], FTYPE (*dP)[NBIGM], FTYPE (*etai)[NUMTRUEEOMSETS][NBIGM], FTYPE (*monoindicator)[NUMMONOINDICATORS][NBIGM], FTYPE (*yprim)[2][NBIGM], FTYPE (*ystecilvar)[2][NBIGM], FTYPE (*yin)[2][NBIGM], FTYPE (*yout)[2][NBIGM], FTYPE (*youtpolycoef)[MAXSPACEORDER][NBIGM], struct of_trueijkp *trueijkp),
223  int stencilvarisnull, int preforder, int whichreduce,
224  int idel, int jdel, int kdel,
225  FTYPE (*primreal)[NSTORE2][NSTORE3][NPR], FTYPE (*stencilvar)[NSTORE2][NSTORE3][NPR2INTERP], FTYPE (*p2interp)[NSTORE2][NSTORE3][NPR2INTERP], FTYPE (*pleft)[NSTORE2][NSTORE3][NPR2INTERP], FTYPE (*pright)[NSTORE2][NSTORE3][NPR2INTERP]);
226 
227 
228 static void slope_lim_linetype_c2e_perijk(
229  LINETYPEPOINTERSFUNCTIONDECLARE,
230  int i, int j, int k,
231  int realisinterp, int whichprimtype, int interporflux,
232  int nprlocalstart, int nprlocalend, int *nprlocallist, int numprims,
233  int dir, int intdir,
234  int is, int ie, int js, int je, int ks, int ke, int di, int dj, int dk, int bs, int ps, int pe, int be,
235  int recontype, int doingweno,
236  void (*pass_1d_line)(int whichquantity, int dir, int do_weight_or_recon, int recontype, int whichreduce, int preforder, int pl, int bs, int ps, int pe, int be, int *minorder, int *maxorder, int *shift, FTYPE (*shockindicator)[NBIGM], FTYPE *stiffindicator, FTYPE (*Vline)[NBIGM], FTYPE (*Pline)[NBIGM], FTYPE (*df)[NBIGM], FTYPE (*dP)[NBIGM], FTYPE (*etai)[NBIGM], FTYPE (*monoindicator)[NBIGM], FTYPE (*yprim)[2][NBIGM], FTYPE (*ystencilvar)[NBIGM], FTYPE (*yin)[NBIGM], FTYPE (*yout)[NBIGM], FTYPE (*youtpolycoef)[NBIGM], struct of_trueijkp *trueijkp),
237  void (*pass_1d_line_multipl)(int MULTIPLTYPE, int whichquantity, int dir, int do_weight_or_recon, int recontype, int whichreduce, int preforder, int bs, int ps, int pe, int be, int *minorder, int *maxorder, int *shift, FTYPE (*shockindicator)[NBIGM], FTYPE *stiffindicator, FTYPE (*Vline)[NBIGM], FTYPE (*Pline)[NBIGM], FTYPE (*df)[NUMDFS][NBIGM], FTYPE (*dP)[NBIGM], FTYPE (*etai)[NUMTRUEEOMSETS][NBIGM], FTYPE (*monoindicator)[NUMMONOINDICATORS][NBIGM], FTYPE (*yprim)[2][NBIGM], FTYPE (*ystecilvar)[2][NBIGM], FTYPE (*yin)[2][NBIGM], FTYPE (*yout)[2][NBIGM], FTYPE (*youtpolycoef)[MAXSPACEORDER][NBIGM], struct of_trueijkp *trueijkp),
238  int stencilvarisnull, int preforder, int whichreduce,
239  int idel, int jdel, int kdel,
240  FTYPE (*primreal)[NSTORE2][NSTORE3][NPR], FTYPE (*stencilvar)[NSTORE2][NSTORE3][NPR2INTERP], FTYPE (*p2interp)[NSTORE2][NSTORE3][NPR2INTERP], FTYPE (*pleft)[NSTORE2][NSTORE3][NPR2INTERP], FTYPE (*pright)[NSTORE2][NSTORE3][NPR2INTERP]);
241 
242 
243 static void slope_lim_linetype_perijk(
244  LINETYPEPOINTERSFUNCTIONDECLARE,
245  int i, int j, int k,
246  int realisinterp, int whichprimtype, int interporflux,
247  int nprlocalstart, int nprlocalend, int *nprlocallist, int numprims,
248  int dir, int intdir,
249  int is, int ie, int js, int je, int ks, int ke, int di, int dj, int dk, int bs, int ps, int pe, int be,
250  int recontype, int doingweno,
251  void (*pass_1d_line)(int whichquantity, int dir, int do_weight_or_recon, int recontype, int whichreduce, int preforder, int pl, int bs, int ps, int pe, int be, int *minorder, int *maxorder, int *shift, FTYPE (*shockindicator)[NBIGM], FTYPE *stiffindicator, FTYPE (*Vline)[NBIGM], FTYPE (*Pline)[NBIGM], FTYPE (*df)[NBIGM], FTYPE (*dP)[NBIGM], FTYPE (*etai)[NBIGM], FTYPE (*monoindicator)[NBIGM], FTYPE (*yprim)[2][NBIGM], FTYPE (*ystencilvar)[NBIGM], FTYPE (*yin)[NBIGM], FTYPE (*yout)[NBIGM], FTYPE (*youtpolycoef)[NBIGM], struct of_trueijkp *trueijkp),
252  void (*pass_1d_line_multipl)(int MULTIPLTYPE, int whichquantity, int dir, int do_weight_or_recon, int recontype, int whichreduce, int preforder, int bs, int ps, int pe, int be, int *minorder, int *maxorder, int *shift, FTYPE (*shockindicator)[NBIGM], FTYPE *stiffindicator, FTYPE (*Vline)[NBIGM], FTYPE (*Pline)[NBIGM], FTYPE (*df)[NUMDFS][NBIGM], FTYPE (*dP)[NBIGM], FTYPE (*etai)[NUMTRUEEOMSETS][NBIGM], FTYPE (*monoindicator)[NUMMONOINDICATORS][NBIGM], FTYPE (*yprim)[2][NBIGM], FTYPE (*ystecilvar)[2][NBIGM], FTYPE (*yin)[2][NBIGM], FTYPE (*yout)[2][NBIGM], FTYPE (*youtpolycoef)[MAXSPACEORDER][NBIGM], struct of_trueijkp *trueijkp),
253  int stencilvarisnull, int preforder, int whichreduce,
254  int weightsplittype, int whichquantity,
255  int idel, int jdel, int kdel,
256  FTYPE (*primreal)[NSTORE2][NSTORE3][NPR], FTYPE (*stencilvar)[NSTORE2][NSTORE3][NPR2INTERP], FTYPE (*p2interpm)[NSTORE2][NSTORE3][NPR], FTYPE (*p2interpp)[NSTORE2][NSTORE3][NPR], FTYPE (*pleft)[NSTORE2][NSTORE3][NPR], FTYPE (*pright)[NSTORE2][NSTORE3][NPR]);
257 
258 
259 static void assign_eno_result(int interporflux, int pl, int bs, int ps, int pe, int be, int i, int j, int k, int idel, int jdel, int kdel, FTYPE (*yout)[NBIGM], FTYPE (*result0)[NSTORE2][NSTORE3][NPR], FTYPE (*result1)[NSTORE2][NSTORE3][NPR]);
260 
261 
263 static void get_1d_line_c2e_multipl(int whichquantity, int dir, int interporflux, int bs, int ps, int pe, int be, int i, int j, int k, FTYPE (*p2interp)[NSTORE2][NSTORE3][NPR2INTERP], FTYPE (*yin)[2][NBIGM], struct of_trueijkp *trueijkp);
264 static void get_1d_line_c2e(int dir, int interporflux, int pl, int bs, int ps, int pe, int be, int i, int j, int k, FTYPE (*p2interp)[NSTORE2][NSTORE3][NPR2INTERP],FTYPE *yin, struct of_trueijkp *trueijkp);
265 static void get_1d_line_shockarray(int dir, int interporflux, int bs, int ps, int pe, int be, int i, int j, int k, FTYPE (*shockarray)[NSTORE1][NSTORE2][NSTORE3], FTYPE (*shockindicator)[NBIGM], struct of_trueijkp *trueijkp);
266 
268 static void assign_eno_result_c2e_multipl(int whichprimtype, int interporflux, int bs, int ps, int pe, int be, int i, int j, int k, int idel, int jdel, int kdel, FTYPE (*yout)[2][NBIGM], FTYPE (*result0)[NSTORE2][NSTORE3][NPR2INTERP], FTYPE (*result1)[NSTORE2][NSTORE3][NPR2INTERP]);
269 static void assign_eno_result_c2e(int interporflux, int pl, int bs, int ps, int pe, int be, int i, int j, int k, int idel, int jdel, int kdel, FTYPE (*yout)[NBIGM], FTYPE (*result0)[NSTORE2][NSTORE3][NPR2INTERP], FTYPE (*result1)[NSTORE2][NSTORE3][NPR2INTERP]);
271 static void get_1d_line_c2e_gammaeffhydro(int dir, int interporflux, int pl, int bs, int ps, int pe, int be, int i, int j, int k, int idel, int jdel, int kdel, FTYPE (*p2interp)[NSTORE2][NSTORE3][NPR2INTERP], FTYPE *yin, struct of_trueijkp *trueijkp);
272 
274 static void get_df_line_gen_new(int realisinterp, int doingweno, int whichprimtype, int interporflux, int recontype, int dir, int whichreduce, int preforder, int bs, int ps, int pe, int be, int *minorder, int *maxorder, int *shift, FTYPE (*yprim)[2][NBIGM], FTYPE (*yin)[2][NBIGM], FTYPE (*df)[NUMDFS][NBIGM], FTYPE (**drhoptr)[NBIGM], FTYPE (**dPptr)[NBIGM], FTYPE *stiffindicator, FTYPE (*Vline)[NBIGM], FTYPE (*Pline)[NBIGM], struct of_trueijkp *trueijkp);
275 static void get_df_line_paraline(int realisinterp, int doingweno, int whichprimtype, int interporflux, int recontype, int dir, int whichreduce, int preforder, int bs, int ps, int pe, int be, int *minorder, int *maxorder, int *shift, FTYPE (*yprim)[2][NBIGM], FTYPE (*yin)[2][NBIGM], FTYPE (*df)[NUMDFS][NBIGM], FTYPE (**drhoptr)[NBIGM], FTYPE (**dPptr)[NBIGM], FTYPE *stiffindicator, FTYPE (*Vline)[NBIGM], FTYPE (*Pline)[NBIGM], struct of_trueijkp *trueijkp);
276 static int compute_df_line(int doingweno,int interporflux, int recontype, int whichreduce, int preforder, int pl, int bs, int ps, int pe, int be, int *minorder, int *maxorder, int *shift, FTYPE *yin, FTYPE (*df)[NBIGM]);
277 static int compute_df_line_paraline(int doingweno,int interporflux, int recontype, int whichreduce, int preforder, int pl, int bs, int ps, int pe, int be, int *minorder, int *maxorder, int *shift, FTYPE *yin, FTYPE (*df)[NBIGM]);
278 static int compute_df_line_new(int doingweno, int whichprimtype, int interporflux, int recontype, int dir, int whichreduce, int preforder, int bs, int ps, int pe, int be, int *minorder, int *maxorder, int *shift, FTYPE (*yprim)[2][NBIGM], FTYPE *df, FTYPE *stiffindicator, FTYPE (*Vline)[NBIGM], FTYPE (*Pline)[NBIGM], struct of_trueijkp *trueijkp);
279 static int compute_df_line_formono(int doingweno,int interporflux, int recontype, int whichreduce, int preforder, int pl, int bs, int ps, int pe, int be, int *minorder, int *maxorder, int *shift, FTYPE *yin, FTYPE (*df)[NBIGM]);
280 static void get_1d_line(int dir, int interporflux, int pl, int bs, int ps, int pe, int be, int i, int j, int k, FTYPE (*p2interpm)[NSTORE2][NSTORE3][NPR],FTYPE (*p2interpp)[NSTORE2][NSTORE3][NPR], FTYPE (*yin)[NBIGM], struct of_trueijkp *trueijkp);
281 static int get_V_and_P(int whichprimtype, int interporflux, int dir, int bs, int ps, int pe, int be, int i, int j, int k, int idel, int jdel, int kdel, FTYPE (*yin)[2][NBIGM], FTYPE (*Vline)[NBIGM], FTYPE (*Pline)[NBIGM], struct of_trueijkp *trueijkp);
282 static int get_shock_indicator(int whichprimtype, int interporflux, int dir, int bs, int ps, int pe, int be, int i, int j, int k, int idel, int jdel, int kdel, FTYPE (*yin)[2][NBIGM], FTYPE (*Vline)[NBIGM], FTYPE (*Pline)[NBIGM], FTYPE (*shockindicator)[NBIGM], struct of_trueijkp *trueijkp);
283 static int get_contact_indicator(int realisinterp, int whichprimtype, int interporflux, int dir, int bs, int ps, int pe, int be, int i, int j, int k, int idel, int jdel, int kdel, FTYPE (*yin)[2][NBIGM], FTYPE (*Vline)[NBIGM], FTYPE (*Pline)[NBIGM], FTYPE (*etai)[NUMTRUEEOMSETS][NBIGM]);
284 static void causal_shift_order(int whichquantitiy, int interporflux, int dir, int preforder, int bs, int ps, int pe, int be, int i, int j, int k, int idel, int jdel, int kdel, int *shift, int *minorder, int *maxorder);
285 
286 
287 static void set_preforder(int dir, int interporflux, int *preforder, int*whichreduce);
288 static int get_recon_type(int interporflux);
289 static void (*pass_1d_line)(int whichquantity, int dir, int do_weight_or_recon, int recontype, int whichreduce, int preforder, int pl, int bs, int ps, int pe, int be, int *minorder, int *maxorder, int *shift, FTYPE (*shockindicator)[NBIGM], FTYPE *stiffindicator, FTYPE (*Vline)[NBIGM], FTYPE (*Pline)[NBIGM], FTYPE (*df)[NBIGM], FTYPE (*dP)[NBIGM], FTYPE (*etai)[NBIGM], FTYPE (*monoindicator)[NBIGM], FTYPE (*yprim)[2][NBIGM], FTYPE (*ystencilvar)[NBIGM], FTYPE (*yin)[NBIGM], FTYPE (*yout)[NBIGM], FTYPE (*youtpolycoef)[NBIGM], struct of_trueijkp *trueijkp);
290 static void (*pass_1d_line_multipl)(int MULTIPLTYPE, int whichquantity, int dir, int do_weight_or_recon, int recontype, int whichreduce, int preforder, int bs, int ps, int pe, int be, int *minorder, int *maxorder, int *shift, FTYPE (*shockindicator)[NBIGM], FTYPE *stiffindicator, FTYPE (*Vline)[NBIGM], FTYPE (*Pline)[NBIGM], FTYPE (*df)[NUMDFS][NBIGM], FTYPE (*dP)[NBIGM], FTYPE (*etai)[NUMTRUEEOMSETS][NBIGM], FTYPE (*monoindicator)[NUMMONOINDICATORS][NBIGM], FTYPE (*yprim)[2][NBIGM], FTYPE (*ystecilvar)[2][NBIGM], FTYPE (*yin)[2][NBIGM], FTYPE (*yout)[2][NBIGM], FTYPE (*youtpolycoef)[MAXSPACEORDER][NBIGM], struct of_trueijkp *trueijkp);
291 
292 
293 static void set_interp_loop(int withshifts, int interporflux, int dir, int loc, int continuous, int *intdir, int *is, int *ie, int *js, int *je, int *ks, int *ke, int *di, int *dj, int *dk, int *bs, int *ps, int *pe, int *be);
294 static void set_interp_loop_expanded(int withshifts, int interporflux, int dir, int loc, int continuous, int *intdir, int *is, int *ie, int *js, int *je, int *ks, int *ke, int *di, int *dj, int *dk, int *bs, int *ps, int *pe, int *be);
295 
296 
297 
298 
299 
300 
302 void slope_lim_linetype_c2e(int realisinterp, int whichprimtype, int interporflux, int dir, int idel, int jdel, int kdel, FTYPE (*primreal)[NSTORE2][NSTORE3][NPR], FTYPE (*stencilvar)[NSTORE2][NSTORE3][NPR2INTERP], FTYPE (*p2interp)[NSTORE2][NSTORE3][NPR2INTERP], FTYPE (*pleft)[NSTORE2][NSTORE3][NPR2INTERP], FTYPE (*pright)[NSTORE2][NSTORE3][NPR2INTERP])
303 {
304 
305  // VS doesn't compile unless this after all function/variable definitions
306  // VS also wants no semicolons back-to-back?
307  LINETYPEDEFINES1;
308 
309 
310 
311  // RESIDUAL NOTES:
312  // DEATHMARK
313  // yin is of size yin[bf-bs+1] + yin[bs] is starting point for yin data. Size: yout[2][pf-ps+1] yout[0/1][0] is first data point
314  // shift is same size as yout. It indicates the location of the center of the stencil w.r.t. the point of interest (i). So a shift if 2 would mean center of stencil is at i+2.
315  // primitives of some form are interpolated here and their dimensions are in standard [i][j][k] form.
316 
317 
319  //
320  // Define which quantities (pl) to operate on
321  //
323 
324  setup_nprlocalist(whichprimtype,&nprlocalstart,&nprlocalend,nprlocallist,&numprims);
325 
326 
328  //
329  // Define loop over starting positions and range of loop for each starting position
330  //
332 
333 
334  withshifts=1; // need with shifts since SUPERGENLOOP below has no shifts and shouldn't have shifts since (e.g.) for dir=1 ie=is and shifts would split the loop into 3D type loop instead of over starting positions
335 
336  int loc =CENT;
337  int continuous=0;
338  set_interp_loop_gen(withshifts, interporflux, dir, loc, continuous, &intdir, &is, &ie, &js, &je, &ks, &ke, &di, &dj, &dk, &bs, &ps, &pe, &be);
339 
341  //
342  // get reconstruction type (c2e, a2c, c2a)
343  //
345 
346  recontype=get_recon_type(interporflux);
347 
348 
350  //
351  // determine which pass_1d_line and pass_1d_line_multipl to use
352  //
353  // assume recontype is CVT_C2E
354  //
356 
357  doingweno=0;
358  if(
359  WENOINTERPTYPE(lim[dir])
360  ){
361  doingweno=1;
362  pass_1d_line=&pass_1d_line_weno;
363  pass_1d_line_multipl=&pass_1d_line_multipl_weno;
364  }
365  else if(
366  lim[dir]==PARALINE
367  ){
368  pass_1d_line=&pass_1d_line_paraline;
369  pass_1d_line_multipl=&pass_1d_line_multipl_paraline;
370  }
371  else{
372  dualfprintf(fail_file,"No such defined pass_1d_line for lim[dir=%d]=%d\n",dir,lim[dir]);
373  myexit(91758726);
374  }
375 
377  //
378  // set stencil variable indicator
379  //
381 
382  if(stencilvar==NULL) stencilvarisnull=1;
383  else stencilvarisnull=0;
384 
385 
387  //
388  // set preferred order
389  //
391 
392  set_preforder(dir, interporflux, &preforder, &whichreduce);
393 
394 
396  //
397  // LOOP OVER STARTING POSITIONS: this loop is over starting positions only
398  //
400  // SUPERGENLOOP(i,j,k,is,ie,js,je,ks,ke,di,dj,dk){
401 
402 #if(STOREFLUXSTATE)
403  // don't need to get EOS if already stored since only access already-stored data (i.e. we don't recompute state)
404 #pragma omp parallel
405 #else
406 #pragma omp parallel OPENMPGLOBALPRIVATEFORSTATEANDGEOMINTERP
407 #endif
408  {
409  OPENMP3DLOOPVARSDEFINE;
410  int i,j,k;
411  LINETYPEDEFINEPOINTERS; // must be within parallel region
412  LINETYPEDEFINESOTHER; // must be within parallel region
413  LINETYPEDEFINESMEMORY; // must be within parallel region
414 
415 
417  // perform pointer shifts
419  LINETYPESHIFTS;
420 
421 
422  OPENMP3DLOOPSETUPSUPERGEN(is,ie,js,je,ks,ke,di,dj,dk);
423 #pragma omp for schedule(OPENMPSCHEDULE(),OPENMPCHUNKSIZE(blocksize))
424  OPENMP3DLOOPBLOCK{
425  OPENMP3DLOOPBLOCK2IJK(i,j,k);
426 
427 
428 #if(STORESHOCKINDICATOR)
429  // avoids expensive get_V_and_P(), Ficalc, etc. since not necessary to do this for each call to this function since just diffusive-type correction. Just ensure average (or use max or min with fabs() ) original quantity to right location when finally used.
430  slope_lim_linetype_c2e_perijk_precomputeindicators(
431  LINETYPEPOINTERSPASSFUNCTIONARG,
432  i,j,k,
433  realisinterp, whichprimtype, interporflux,
434  nprlocalstart, nprlocalend, nprlocallist, numprims,
435  dir, intdir,
436  is, ie, js, je, ks, ke, di, dj, dk, bs, ps, pe, be,
437  recontype, doingweno,
438  pass_1d_line, pass_1d_line_multipl,
439  stencilvarisnull, preforder, whichreduce,
440  idel, jdel, kdel,
441  primreal, stencilvar, p2interp, pleft, pright);
442 #else
443  slope_lim_linetype_c2e_perijk(
444  LINETYPEPOINTERSPASSFUNCTIONARG,
445  i,j,k,
446  realisinterp, whichprimtype, interporflux,
447  nprlocalstart, nprlocalend, nprlocallist, numprims,
448  dir, intdir,
449  is, ie, js, je, ks, ke, di, dj, dk, bs, ps, pe, be,
450  recontype, doingweno,
451  pass_1d_line, pass_1d_line_multipl,
452  stencilvarisnull, preforder, whichreduce,
453  idel, jdel, kdel,
454  primreal, stencilvar, p2interp, pleft, pright);
455 #endif
456 
457  } // done with main loop over starting points
458  }// end parallel region
459 
460 
461 }
462 
463 
464 
465 
467 void slope_lim_linetype_c2e_perijk_precomputeindicators(
468  LINETYPEPOINTERSFUNCTIONDECLARE,
469  int i, int j, int k,
470  int realisinterp, int whichprimtype, int interporflux,
471  int nprlocalstart, int nprlocalend, int *nprlocallist, int numprims,
472  int dir, int intdir,
473  int is, int ie, int js, int je, int ks, int ke, int di, int dj, int dk, int bs, int ps, int pe, int be,
474  int recontype, int doingweno,
475  void (*pass_1d_line)(int whichquantity, int dir, int do_weight_or_recon, int recontype, int whichreduce, int preforder, int pl, int bs, int ps, int pe, int be, int *minorder, int *maxorder, int *shift, FTYPE (*shockindicator)[NBIGM], FTYPE *stiffindicator, FTYPE (*Vline)[NBIGM], FTYPE (*Pline)[NBIGM], FTYPE (*df)[NBIGM], FTYPE (*dP)[NBIGM], FTYPE (*etai)[NBIGM], FTYPE (*monoindicator)[NBIGM], FTYPE (*yprim)[2][NBIGM], FTYPE (*ystencilvar)[NBIGM], FTYPE (*yin)[NBIGM], FTYPE (*yout)[NBIGM], FTYPE (*youtpolycoef)[NBIGM], struct of_trueijkp *trueijkp),
476  void (*pass_1d_line_multipl)(int MULTIPLTYPE, int whichquantity, int dir, int do_weight_or_recon, int recontype, int whichreduce, int preforder, int bs, int ps, int pe, int be, int *minorder, int *maxorder, int *shift, FTYPE (*shockindicator)[NBIGM], FTYPE *stiffindicator, FTYPE (*Vline)[NBIGM], FTYPE (*Pline)[NBIGM], FTYPE (*df)[NUMDFS][NBIGM], FTYPE (*dP)[NBIGM], FTYPE (*etai)[NUMTRUEEOMSETS][NBIGM], FTYPE (*monoindicator)[NUMMONOINDICATORS][NBIGM], FTYPE (*yprim)[2][NBIGM], FTYPE (*ystecilvar)[2][NBIGM], FTYPE (*yin)[2][NBIGM], FTYPE (*yout)[2][NBIGM], FTYPE (*youtpolycoef)[MAXSPACEORDER][NBIGM], struct of_trueijkp *trueijkp),
477  int stencilvarisnull, int preforder, int whichreduce,
478  int idel, int jdel, int kdel,
479  FTYPE (*primreal)[NSTORE2][NSTORE3][NPR], FTYPE (*stencilvar)[NSTORE2][NSTORE3][NPR2INTERP], FTYPE (*p2interp)[NSTORE2][NSTORE3][NPR2INTERP], FTYPE (*pleft)[NSTORE2][NSTORE3][NPR2INTERP], FTYPE (*pright)[NSTORE2][NSTORE3][NPR2INTERP])
480 {
481 
482  // VS doesn't compile unless this after all function/variable definitions
483  // VS also wants no semicolons back-to-back?
484  LINETYPEDEFINES2;
485 
486 
487 
488 
489 
490 
491  // RESIDUAL NOTES:
492  // DEATHMARK
493  // yin is of size yin[bf-bs+1] + yin[bs] is starting point for yin data. Size: yout[2][pf-ps+1] yout[0/1][0] is first data point
494  // shift is same size as yout. It indicates the location of the center of the stencil w.r.t. the point of interest (i). So a shift if 2 would mean center of stencil is at i+2.
495  // primitives of some form are interpolated here and their dimensions are in standard [i][j][k] form.
496 
497 
498 
499  if(doingweno){ // GODMARK: shift,minorder,maxorder only used by weno right now
501  //
502  // Figure out shifting of stencil and order of stencil to more closely match with causality.
503  // If no stencil shifting, then set shift and min/max order to default
504  //
506  causal_shift_order(whichprimtype, interporflux, dir, preforder, bs, ps, pe, be, i, j, k, idel, jdel ,kdel, shift, minorder, maxorder);
507  }
508 
509  // subloop loops from starting position to end of that line
510 
511 
513  //
514  // GET 1D LINES [assume shock indicator (or other things) already computed and exists as part of "primitives"]
515  //
517 
518  // for this we don't need realisinterp and should never use yprim or primreal
519  get_1d_line_c2e_multipl(ENOPRIMITIVE, dir, interporflux, bs, ps, pe, be, i, j, k, p2interp, yin, &trueijkp);
520 
521  // get shock indicator from stored array
522  if(whichreduce == WENO_REDUCE_TYPE_PPM || CONTACTINDICATOR || COMPUTEDRHODP|| SHOCKINDICATOR ){
523 #if(STORESHOCKINDICATOR==0)
524  dualfprintf(fail_file,"If using simple c2e method, need to set STORESHOCKINDICATOR=1\n");
525  myexit(394834);
526 #endif
527 
528  get_1d_line_shockarray(dir, interporflux, bs, ps, pe, be, i, j, k, GLOBALPOINT(shockindicatorarray), shockindicator,&trueijkp);
529  }
530 
531 
532 
533 #if(CONTACTINDICATOR)
534  dualfprintf(fail_file,"If using simple c2e method with contacts, need to setup flux.c's compute_and_store_fluxstatecent().\n");
535  myexit(394834);
536 #endif
537 
538 
539 
540  if(PARALINEUSESMONO||doingweno){ // let paraline use mono (which needs ystencilvar)
542  //
543  // 1D GET LINE of stencilvar
544  //
546  if(!stencilvarisnull){
547  get_1d_line_c2e_multipl(ENOPRIMITIVE, dir, interporflux, bs, ps, pe, be, i, j, k, stencilvar, ystencilvar,&trueijkp);
548  }
549  else{
550  // assign rather than pointer assign since yin//ystencilvar may be modified and don't want to have to worry if was
551  // notice only requesting yin[0], not yin[1], so not good for FLUXSPLIT-- GODMARK
552  NUMPRIMLOOP(pliter,pl) for(yiter=bs;yiter<=be;yiter++) ystencilvar[pl][0][yiter]=yin[pl][0][yiter];
553  }
554  }
555 
556 
557 
558  if(doingweno){
560  //
561  // get df's for contact indicator or in general (use this to get stiffindicator too)
562  //
563  // gets Vline and Pline now too if doingweno==1
565  get_df_line_gen_new(realisinterp, doingweno,whichprimtype,interporflux,recontype,dir,whichreduce,preforder, bs, ps, pe, be, minorder, maxorder, shift, yprim, yin, df, &drho, &dP,stiffindicator,Vline,Pline,&trueijkp);
566  }
567  else{
568  get_df_line_paraline(realisinterp, doingweno,whichprimtype,interporflux,recontype,dir,whichreduce,preforder, bs, ps, pe, be, minorder, maxorder, shift, yprim, yin, df, &drho, &dP,stiffindicator,Vline,Pline,&trueijkp);
569  }
570 
571 
572 
573  if(PARALINEUSESMONO||doingweno){ // let paraline use mono (which needs dfformono[][][])
575  //
576  // 1D GET LINE of mono df's if needed
577  //
578  // Point is that MONO (e.g. SMONO) uses df to compute if cusp exists (even if DO_SMONO_???==0).
579  // Cusp sets monoindicator that sets effective stencil, so need to use stencilvar for this
580  //
582  if(!stencilvarisnull){
583  // notice only requesting yin[0], not yin[1], so not good for FLUXSPLIT-- GODMARK
584  NUMPRIMLOOP(pliter,pl) compute_df_line_formono(doingweno,interporflux, recontype, whichreduce, preforder, pl, bs, ps, pe, be, minorder, maxorder, shift, yin[pl][0], dfformono[pl]);
585  }
586  else{
587  // assign rather than pointer assign since yin//ystencilvar may be modified and don't want to have to worry if was
588  // should be expensive
589  NUMPRIMLOOP(pliter,pl) for(yiter=bs;yiter<=be;yiter++) for(dfiter=0;dfiter<NUMDFS;dfiter++) dfformono[pl][dfiter][yiter]=df[pl][dfiter][yiter];
590  }
591 
592 
594  //
595  // Compute monotonicity indicator and if monotonic or rough then compute interface value
596  // Assume only needed by WENO routines
597  //
599 
600  if(stencilvar==NULL){
601  NUMPRIMLOOP(pliter,pl) compute_monotonicity_line(recontype, whichreduce, preforder, pl, bs, ps, pe, be, minorder, maxorder, shift, shockindicator, dfformono[pl], monoindicator[pl], yin[pl][0], yout[pl], youtpolycoef[pl]);
602  }
603  else{
604  // then split indicator from value
605  NUMPRIMLOOP(pliter,pl) compute_monotonicity_line_indicatoronly(recontype, whichreduce, preforder, pl, bs, ps, pe, be, minorder, maxorder, shift, shockindicator, dfformono[pl], monoindicator[pl], ystencilvar[pl][0], yout[pl], youtpolycoef[pl]);
606 
607  NUMPRIMLOOP(pliter,pl) compute_monotonicity_line_valueonly(recontype, whichreduce, preforder, pl, bs, ps, pe, be, minorder, maxorder, shift, shockindicator, dfformono[pl], monoindicator[pl], yin[pl][0], yout[pl], youtpolycoef[pl]);
608  }
609  }
610 
611 
612 
613 
615  //
616  // PASS 1D LINE
617  //
619 
620  // assume normally want all pl's
621  pass_1d_line_multipl( DO_SPLITC2E, ENOPRIMITIVE, dir, ALL_CALC, recontype, whichreduce, preforder, bs, ps, pe, be, minorder, maxorder, shift, shockindicator, stiffindicator, Vline, Pline, df, dP, etai, monoindicator, yprim, ystencilvar, yin, yout, youtpolycoef,&trueijkp);
622 
624  //
625  // now have 1D line of point data corresponding to left and right interpolated values
626  //
627  // Assign result back to pleft/pright (or just pleft if one result)
628  //
630  assign_eno_result_c2e_multipl(ENOPRIMITIVE, recontype, bs, ps, pe, be, i, j, k, idel, jdel, kdel, yout, pleft, pright);
631 
632 }
633 
634 
635 
636 
638 void slope_lim_linetype_c2e_perijk(
639  LINETYPEPOINTERSFUNCTIONDECLARE,
640  int i, int j, int k,
641  int realisinterp, int whichprimtype, int interporflux,
642  int nprlocalstart, int nprlocalend, int *nprlocallist, int numprims,
643  int dir, int intdir,
644  int is, int ie, int js, int je, int ks, int ke, int di, int dj, int dk, int bs, int ps, int pe, int be,
645  int recontype, int doingweno,
646  void (*pass_1d_line)(int whichquantity, int dir, int do_weight_or_recon, int recontype, int whichreduce, int preforder, int pl, int bs, int ps, int pe, int be, int *minorder, int *maxorder, int *shift, FTYPE (*shockindicator)[NBIGM], FTYPE *stiffindicator, FTYPE (*Vline)[NBIGM], FTYPE (*Pline)[NBIGM], FTYPE (*df)[NBIGM], FTYPE (*dP)[NBIGM], FTYPE (*etai)[NBIGM], FTYPE (*monoindicator)[NBIGM], FTYPE (*yprim)[2][NBIGM], FTYPE (*ystencilvar)[NBIGM], FTYPE (*yin)[NBIGM], FTYPE (*yout)[NBIGM], FTYPE (*youtpolycoef)[NBIGM], struct of_trueijkp *trueijkp),
647  void (*pass_1d_line_multipl)(int MULTIPLTYPE, int whichquantity, int dir, int do_weight_or_recon, int recontype, int whichreduce, int preforder, int bs, int ps, int pe, int be, int *minorder, int *maxorder, int *shift, FTYPE (*shockindicator)[NBIGM], FTYPE *stiffindicator, FTYPE (*Vline)[NBIGM], FTYPE (*Pline)[NBIGM], FTYPE (*df)[NUMDFS][NBIGM], FTYPE (*dP)[NBIGM], FTYPE (*etai)[NUMTRUEEOMSETS][NBIGM], FTYPE (*monoindicator)[NUMMONOINDICATORS][NBIGM], FTYPE (*yprim)[2][NBIGM], FTYPE (*ystecilvar)[2][NBIGM], FTYPE (*yin)[2][NBIGM], FTYPE (*yout)[2][NBIGM], FTYPE (*youtpolycoef)[MAXSPACEORDER][NBIGM], struct of_trueijkp *trueijkp),
648  int stencilvarisnull, int preforder, int whichreduce,
649  int idel, int jdel, int kdel,
650  FTYPE (*primreal)[NSTORE2][NSTORE3][NPR], FTYPE (*stencilvar)[NSTORE2][NSTORE3][NPR2INTERP], FTYPE (*p2interp)[NSTORE2][NSTORE3][NPR2INTERP], FTYPE (*pleft)[NSTORE2][NSTORE3][NPR2INTERP], FTYPE (*pright)[NSTORE2][NSTORE3][NPR2INTERP])
651 {
652 
653  // VS doesn't compile unless this after all function/variable definitions
654  // VS also wants no semicolons back-to-back?
655  LINETYPEDEFINES2;
656 
657 
658 
659 
660 
661 
662  // RESIDUAL NOTES:
663  // DEATHMARK
664  // yin is of size yin[bf-bs+1] + yin[bs] is starting point for yin data. Size: yout[2][pf-ps+1] yout[0/1][0] is first data point
665  // shift is same size as yout. It indicates the location of the center of the stencil w.r.t. the point of interest (i). So a shift if 2 would mean center of stencil is at i+2.
666  // primitives of some form are interpolated here and their dimensions are in standard [i][j][k] form.
667 
668 
669 
670  if(doingweno){ // GODMARK: shift,minorder,maxorder only used by weno right now
672  //
673  // Figure out shifting of stencil and order of stencil to more closely match with causality.
674  // If no stencil shifting, then set shift and min/max order to default
675  //
677  causal_shift_order(whichprimtype, interporflux, dir, preforder, bs, ps, pe, be, i, j, k, idel, jdel ,kdel, shift, minorder, maxorder);
678  }
679 
680  // subloop loops from starting position to end of that line
681 
682 
684  //
685  // GET 1D LINES
686  //
688 
689  if(whichreduce == WENO_REDUCE_TYPE_PPM || CONTACTINDICATOR || COMPUTEDRHODP|| SHOCKINDICATOR ){
690  if(realisinterp){
691  // means primreal=p2interp
692  // then get all primitives and assign to normal yin
693  NUMPRIMLOOP(pliter,pl) get_1d_line_c2e(dir, interporflux, pl, bs, ps, pe, be, i, j, k, primreal, yin[pl][0],&trueijkp);
694  // below means can't modify yin without having changed yprim (not general GODMARK)
695  yprim=yin;
696  }
697  else{
698  // means primreal!=p2interp
699  // then need to get separate p2interp line
700  get_1d_line_c2e_multipl(ENOPRIMITIVE, dir, interporflux, bs, ps, pe, be, i, j, k, p2interp, yin, &trueijkp);
701 
702  // get prim line
703  if(whichreduce == WENO_REDUCE_TYPE_PPM || SHOCKINDICATOR ){
704  // then need to get primitives separately from interpolated quantities
705  PALLREALLOOP(pl) get_1d_line_c2e(dir, interporflux, pl, bs, ps, pe, be, i, j, k, primreal, yprim[pl][0],&trueijkp);
706  }
707  else if(CONTACTINDICATOR||COMPUTEDRHODP){
708  // then only need RHO and UU
709  get_1d_line_c2e(dir, interporflux, RHO, bs, ps, pe, be, i, j, k, primreal, yprim[RHO][0],&trueijkp);
710  get_1d_line_c2e(dir, interporflux, UU, bs, ps, pe, be, i, j, k, primreal, yprim[UU][0],&trueijkp);
711  }
712 
713 #if(VSQ!=-100)
714  //put the effective hydro value of $\gamma^2$ into the VSQ element of the primitive array that will be passed to reconstruction
715  get_1d_line_c2e_gammaeffhydro(dir, interporflux, VSQ, bs, ps, pe, be, i, j, k, idel, jdel, kdel, p2interp, yprim[VSQ][0],&trueijkp);
716 #endif
717 
718  }
719  }
720  else{
721  // then no ppm reduce or contact indicator, so just get normal line
723  //
724  // 1D GET LINE (interpolated quantity is never the same as the prims4shock quantity since not doing c2e here)
725  //
727  get_1d_line_c2e_multipl(ENOPRIMITIVE, dir, interporflux, bs, ps, pe, be, i, j, k, p2interp, yin, &trueijkp);
728  }
729 
730 
731 
732  if(PARALINEUSESMONO||doingweno){ // let paraline use mono (which needs ystencilvar)
734  //
735  // 1D GET LINE of stencilvar
736  //
738  if(!stencilvarisnull){
739  get_1d_line_c2e_multipl(ENOPRIMITIVE, dir, interporflux, bs, ps, pe, be, i, j, k, stencilvar, ystencilvar,&trueijkp);
740  }
741  else{
742  // assign rather than pointer assign since yin//ystencilvar may be modified and don't want to have to worry if was
743  // notice only requesting yin[0], not yin[1], so not good for FLUXSPLIT-- GODMARK
744  NUMPRIMLOOP(pliter,pl) for(yiter=bs;yiter<=be;yiter++) ystencilvar[pl][0][yiter]=yin[pl][0][yiter];
745  }
746  }
747 
748 
749 
750  if(doingweno){
752  //
753  // get df's for contact indicator or in general (use this to get stiffindicator too)
754  //
755  // gets Vline and Pline now too if doingweno==1
757  get_df_line_gen_new(realisinterp, doingweno,whichprimtype,interporflux,recontype,dir,whichreduce,preforder, bs, ps, pe, be, minorder, maxorder, shift, yprim, yin, df, &drho, &dP,stiffindicator,Vline,Pline,&trueijkp);
758  }
759  else{
760  get_df_line_paraline(realisinterp, doingweno,whichprimtype,interporflux,recontype,dir,whichreduce,preforder, bs, ps, pe, be, minorder, maxorder, shift, yprim, yin, df, &drho, &dP,stiffindicator,Vline,Pline,&trueijkp);
761  }
762 
763 
764  if(PARALINEUSESMONO||doingweno){ // let paraline use mono (which needs dfformono[][][])
766  //
767  // 1D GET LINE of mono df's if needed
768  //
769  // Point is that MONO (e.g. SMONO) uses df to compute if cusp exists (even if DO_SMONO_???==0).
770  // Cusp sets monoindicator that sets effective stencil, so need to use stencilvar for this
771  //
773  if(!stencilvarisnull){
774  // notice only requesting yin[0], not yin[1], so not good for FLUXSPLIT-- GODMARK
775  NUMPRIMLOOP(pliter,pl) compute_df_line_formono(doingweno,interporflux, recontype, whichreduce, preforder, pl, bs, ps, pe, be, minorder, maxorder, shift, yin[pl][0], dfformono[pl]);
776  }
777  else{
778  // assign rather than pointer assign since yin//ystencilvar may be modified and don't want to have to worry if was
779  // should be expensive
780  NUMPRIMLOOP(pliter,pl) for(yiter=bs;yiter<=be;yiter++) for(dfiter=0;dfiter<NUMDFS;dfiter++) dfformono[pl][dfiter][yiter]=df[pl][dfiter][yiter];
781  }
782  }
783 
784 
785  if(!doingweno){ // need this because get_df_line_gen_new() didn't get Vline,Pline if doingweno==0
787  //
788  // 1D P and V
789  //
791  if(SHOCKINDICATOR || CONTACTINDICATOR){
792  get_V_and_P(whichprimtype, interporflux, dir, bs, ps, pe, be, i, j, k, idel, jdel, kdel, yprim, Vline, Pline,&trueijkp);
793  }
794  }
795 
796 
797 
798  if(whichreduce == WENO_REDUCE_TYPE_PPM || SHOCKINDICATOR ){
800  //
801  // 1D GET SHOCK INDICATOR
802  //
804  //use primitive values that correspond to the quantities being interpolated
805  get_shock_indicator(whichprimtype, interporflux, dir, bs, ps, pe, be, i, j, k, idel, jdel ,kdel, yprim, Vline, Pline, shockindicator,&trueijkp);
806  }
807 
808 
809  if(CONTACTINDICATOR){
811  //
812  // 1D GET CONTACT INDICATOR
813  //
815  //use primitive values that correspond to the quantities being interpolated
816  get_contact_indicator(realisinterp, whichprimtype, interporflux, dir, bs, ps, pe, be, i, j, k, idel, jdel ,kdel, yin, Vline, Pline, etai);
817  }
818 
819 
820 
821 
822  if(PARALINEUSESMONO||doingweno){ // let paraline use mono
824  //
825  // Compute monotonicity indicator and if monotonic or rough then compute interface value
826  // Assume only needed by WENO routines
827  //
829 
830  if(stencilvar==NULL){
831  NUMPRIMLOOP(pliter,pl) compute_monotonicity_line(recontype, whichreduce, preforder, pl, bs, ps, pe, be, minorder, maxorder, shift, shockindicator, dfformono[pl], monoindicator[pl], yin[pl][0], yout[pl], youtpolycoef[pl]);
832  }
833  else{
834  // then split indicator from value
835  NUMPRIMLOOP(pliter,pl) compute_monotonicity_line_indicatoronly(recontype, whichreduce, preforder, pl, bs, ps, pe, be, minorder, maxorder, shift, shockindicator, dfformono[pl], monoindicator[pl], ystencilvar[pl][0], yout[pl], youtpolycoef[pl]);
836 
837  NUMPRIMLOOP(pliter,pl) compute_monotonicity_line_valueonly(recontype, whichreduce, preforder, pl, bs, ps, pe, be, minorder, maxorder, shift, shockindicator, dfformono[pl], monoindicator[pl], yin[pl][0], yout[pl], youtpolycoef[pl]);
838  }
839  }
840 
841 
842 
843 
845  //
846  // PASS 1D LINE
847  //
849 
850  // assume normally want all pl's
851  pass_1d_line_multipl( DO_SPLITC2E, ENOPRIMITIVE, dir, ALL_CALC, recontype, whichreduce, preforder, bs, ps, pe, be, minorder, maxorder, shift, shockindicator, stiffindicator, Vline, Pline, df, dP, etai, monoindicator, yprim, ystencilvar, yin, yout, youtpolycoef,&trueijkp);
852 
854  //
855  // now have 1D line of point data corresponding to left and right interpolated values
856  //
857  // Assign result back to pleft/pright (or just pleft if one result)
858  //
860  assign_eno_result_c2e_multipl(ENOPRIMITIVE, recontype, bs, ps, pe, be, i, j, k, idel, jdel, kdel, yout, pleft, pright);
861 
862 }
863 
864 
865 
866 
867 
868 
869 
870 
871 
872 
876 void slope_lim_linetype(int whichquantity, int interporflux, int dir, int idel, int jdel, int kdel, FTYPE (*primreal)[NSTORE2][NSTORE3][NPR], FTYPE (*stencilvar)[NSTORE2][NSTORE3][NPR], FTYPE (*p2interpm)[NSTORE2][NSTORE3][NPR], FTYPE (*p2interpp)[NSTORE2][NSTORE3][NPR], FTYPE (*pleft)[NSTORE2][NSTORE3][NPR], FTYPE (*pright)[NSTORE2][NSTORE3][NPR])
877 {
878  int whichprimtype=whichquantity;
879  int realisinterp;
880  int weightsplittype;
881  // VS doesn't compile unless this after all function/variable definitions
882  // VS also wants no semicolons back-to-back?
883  LINETYPEDEFINES1;
884 
885 
886 
888  //
889  // Define which quantities (pl) to operate on
890  //
892 
893  setup_nprlocalist(whichprimtype,&nprlocalstart,&nprlocalend,nprlocallist,&numprims);
894 
895  // set realisinterp *after* setting nprlists
896  // below not generaly correct if doing c2e, but always doing c2e with slope_lim_linetype_c2e() anyways -- GODMARK
897  if(interporflux==ENOINTERPTYPE) set_normal_realisinterp(&realisinterp);
898  else realisinterp=0;
899 
900 
901 
902 
903 
904  // RESIDUAL NOTES:
905  // DEATHMARK
906  // yin is of size yin[bf-bs+1] + yin[bs] is starting point for yin data. Size: yout[2][pf-ps+1] yout[0/1][0] is first data point
907  // shift is same size as yout. It indicates the location of the center of the stencil w.r.t. the point of interest (i). So a shift if 2 would mean center of stencil is at i+2.
908  // primitives of some form are interpolated here and their dimensions are in standard [i][j][k] form.
909 
910 
912  //
913  // Define loop over starting positions and range of loop for each starting position
914  //
916 
917  withshifts=1; // need with shifts since SUPERGENLOOP below has no shifts and shouldn't have shifts since (e.g.) for dir=1 ie=is and shifts would split the loop into 3D type loop instead of over starting positions
918 
919  int loc=CENT;
920  int continuous=0;
921  set_interp_loop_gen(withshifts,interporflux, dir, continuous, loc, &intdir, &is, &ie, &js, &je, &ks, &ke, &di, &dj, &dk, &bs, &ps, &pe, &be);
922 
924  //
925  // get reconstruction type (c2e, a2c, c2a)
926  //
928 
929  recontype=get_recon_type(interporflux);
930 
932  //
933  // determine which pass_1d_line and pass_1d_line_multipl to use
934  //
936 
937  doingweno=0;
938  if(
939  recontype==CVT_C2E && WENOINTERPTYPE(lim[dir]) ||
940  recontype==CVT_C2A && WENOINTERPTYPE(avgscheme[dir]) ||
941  recontype==CVT_A2C && WENOINTERPTYPE(avgscheme[dir])
942  ){
943  doingweno=1;
944  pass_1d_line=&pass_1d_line_weno;
945  pass_1d_line_multipl=&pass_1d_line_multipl_weno;
946  }
947  else if(
948  recontype==CVT_C2E && lim[dir]==PARALINE ||
949  recontype==CVT_C2A && avgscheme[dir]==PARALINE ||
950  recontype==CVT_A2C && avgscheme[dir]==PARALINE
951  ){
952  pass_1d_line=&pass_1d_line_paraline;
953  pass_1d_line_multipl=&pass_1d_line_multipl_paraline;
954  }
955  else{
956  dualfprintf(fail_file,"No such defined pass_1d_line for lim[dir=%d]=%d avgscheme[1,2,3]=%d %d %d\n",dir,lim[dir],avgscheme[1],avgscheme[2],avgscheme[3]);
957  myexit(91758727);
958  }
959 
960 
961 
963  //
964  // set stencil variable indicator
965  //
967 
968  if(stencilvar==NULL) stencilvarisnull=1;
969  else stencilvarisnull=0;
970 
971 
973  //
974  // set stencil weight reduction type (some arbitrariness for user definition)
975  //
977 
978  higherorder_set(whichquantity, recontype, &weightsplittype);
979 
981  //
982  // set preferred order
983  //
985 
986  set_preforder(dir, interporflux, &preforder, &whichreduce);
987 
988 
990  //
991  // LOOP OVER STARTING POSITIONS: this loop is over starting positions only
992  //
994  // SUPERGENLOOP(i,j,k,is,ie,js,je,ks,ke,di,dj,dk){
995 
996 #if(STOREFLUXSTATE)
997  // don't need to get EOS if already stored since only access already-stored data (i.e. we don't recompute state)
998 #pragma omp parallel
999 #else
1000 #pragma omp parallel OPENMPGLOBALPRIVATEFORSTATEANDGEOMINTERP
1001 #endif
1002  {
1003  OPENMP3DLOOPVARSDEFINE;
1004  int i,j,k;
1005  LINETYPEDEFINEPOINTERS; // must be within parallel region
1006  LINETYPEDEFINESOTHER; // must be within parallel region
1007  LINETYPEDEFINESMEMORY; // must be within parallel region
1008 
1009 
1011  // perform pointer shifts
1013  LINETYPESHIFTS;
1014 
1015 
1016  OPENMP3DLOOPSETUPSUPERGEN(is,ie,js,je,ks,ke,di,dj,dk);
1017 #pragma omp for schedule(OPENMPSCHEDULE(),OPENMPCHUNKSIZE(blocksize))
1018  OPENMP3DLOOPBLOCK{
1019  OPENMP3DLOOPBLOCK2IJK(i,j,k);
1020 
1021  slope_lim_linetype_perijk(
1022  LINETYPEPOINTERSPASSFUNCTIONARG,
1023  i,j,k,
1024  realisinterp, whichprimtype, interporflux,
1025  nprlocalstart, nprlocalend, nprlocallist, numprims,
1026  dir, intdir,
1027  is, ie, js, je, ks, ke, di, dj, dk, bs, ps, pe, be,
1028  recontype, doingweno,
1029  pass_1d_line, pass_1d_line_multipl,
1030  stencilvarisnull, preforder, whichreduce,
1031  weightsplittype, whichquantity,
1032  idel, jdel, kdel,
1033  primreal, stencilvar, p2interpm, p2interpp, pleft, pright);
1034 
1035 
1036  } // done with main loop over starting points
1037  }// end parallel region
1038 
1039 }
1040 
1041 
1042 
1043 
1044 
1045 
1046 
1047 
1049 void slope_lim_linetype_perijk(
1050  LINETYPEPOINTERSFUNCTIONDECLARE,
1051  int i, int j, int k,
1052  int realisinterp, int whichprimtype, int interporflux,
1053  int nprlocalstart, int nprlocalend, int *nprlocallist, int numprims,
1054  int dir, int intdir,
1055  int is, int ie, int js, int je, int ks, int ke, int di, int dj, int dk, int bs, int ps, int pe, int be,
1056  int recontype, int doingweno,
1057  void (*pass_1d_line)(int whichquantity, int dir, int do_weight_or_recon, int recontype, int whichreduce, int preforder, int pl, int bs, int ps, int pe, int be, int *minorder, int *maxorder, int *shift, FTYPE (*shockindicator)[NBIGM], FTYPE *stiffindicator, FTYPE (*Vline)[NBIGM], FTYPE (*Pline)[NBIGM], FTYPE (*df)[NBIGM], FTYPE (*dP)[NBIGM], FTYPE (*etai)[NBIGM], FTYPE (*monoindicator)[NBIGM], FTYPE (*yprim)[2][NBIGM], FTYPE (*ystencilvar)[NBIGM], FTYPE (*yin)[NBIGM], FTYPE (*yout)[NBIGM], FTYPE (*youtpolycoef)[NBIGM], struct of_trueijkp *trueijkp),
1058  void (*pass_1d_line_multipl)(int MULTIPLTYPE, int whichquantity, int dir, int do_weight_or_recon, int recontype, int whichreduce, int preforder, int bs, int ps, int pe, int be, int *minorder, int *maxorder, int *shift, FTYPE (*shockindicator)[NBIGM], FTYPE *stiffindicator, FTYPE (*Vline)[NBIGM], FTYPE (*Pline)[NBIGM], FTYPE (*df)[NUMDFS][NBIGM], FTYPE (*dP)[NBIGM], FTYPE (*etai)[NUMTRUEEOMSETS][NBIGM], FTYPE (*monoindicator)[NUMMONOINDICATORS][NBIGM], FTYPE (*yprim)[2][NBIGM], FTYPE (*ystecilvar)[2][NBIGM], FTYPE (*yin)[2][NBIGM], FTYPE (*yout)[2][NBIGM], FTYPE (*youtpolycoef)[MAXSPACEORDER][NBIGM], struct of_trueijkp *trueijkp),
1059  int stencilvarisnull, int preforder, int whichreduce,
1060  int weightsplittype, int whichquantity,
1061  int idel, int jdel, int kdel,
1062  FTYPE (*primreal)[NSTORE2][NSTORE3][NPR], FTYPE (*stencilvar)[NSTORE2][NSTORE3][NPR2INTERP], FTYPE (*p2interpm)[NSTORE2][NSTORE3][NPR], FTYPE (*p2interpp)[NSTORE2][NSTORE3][NPR], FTYPE (*pleft)[NSTORE2][NSTORE3][NPR], FTYPE (*pright)[NSTORE2][NSTORE3][NPR])
1063 {
1064  extern int apply_bc_line(int nprlocalstart, int nprlocalend, int*nprlocallist, int doinverse, int iter, int recontype, int bs, int be, FTYPE (*yin)[2][NBIGM],FTYPE (*yout)[2][NBIGM], FTYPE (*youtpolycoef)[MAXSPACEORDER][NBIGM]);
1065 
1066 
1067  int monoi;
1068  // VS doesn't compile unless this after all function/variable definitions
1069  // VS also wants no semicolons back-to-back?
1070  LINETYPEDEFINES2;
1071 
1072 
1073 
1074 
1075 
1076 
1077  if(weightsplittype==CONSTANT_ALL_WEIGHTS){
1078 
1079  if(DOMONOINTERP==NOMONOINTERP){
1080  dualfprintf(fail_file,"Cannot do weightsplittype==CONSTANT_ALL_WEIGHTS if DOMONOINTERP==%d\n",DOMONOINTERP);
1081  myexit(2469238463);
1082  }
1083 
1084  // then just do mono with fixed weights
1085 
1086  // get line (don't need yprim,df,shockindicator,etc.)
1087  NUMPRIMLOOP(pliter,pl) get_1d_line(dir, interporflux, pl, bs, ps, pe, be, i, j, k, p2interpm, p2interpp, yin[pl],&trueijkp);
1088 
1089  // user-defined modification of the line -- usually adjusting line data so interpolation can be higher-order
1090  apply_bc_line(nprlocalstart,nprlocalend,nprlocallist,0,trueijkp.iter,recontype,bs,be,yin,NULL,NULL);
1091 
1092  NUMPRIMLOOP(pliter,pl){
1093  // set all weights equal
1094  for(monoi=bs;monoi<=be;monoi++){
1095  monoindicator[pl][MONOINDTYPE][monoi]=1;
1096  monoindicator[pl][MONOLEFTSET][monoi]=1; // or center
1097  monoindicator[pl][MONORIGHTSET][monoi]=1;
1098  }
1099 
1100  // compute MONO reconstructions given weights for all "pl"
1101  // note that shift, shockindicator, and df not used
1102  compute_monotonicity_line_valueonly(recontype, whichreduce, preforder, pl, bs, ps, pe, be, minorder, maxorder, shift, shockindicator, df[pl], monoindicator[pl], yin[pl][0], yout[pl], youtpolycoef[pl]);
1103  }
1104 
1105 
1106 
1107 
1108 
1109  }
1110  else{
1111 
1112 
1113 
1115  //
1116  // Figure out shifting of stencil and order of stencil to more closely match with causality.
1117  // If no stencil shifting, then set shift and min/max order to default
1118  //
1120  causal_shift_order(whichquantity, interporflux, dir, preforder, bs, ps, pe, be, i, j, k, idel, jdel ,kdel, shift, minorder, maxorder);
1121 
1122 
1123 
1125  //
1126  // GET 1D LINES
1127  //
1129  // subloop loops from starting position to end of that line
1130 
1131  if(whichreduce == WENO_REDUCE_TYPE_PPM || CONTACTINDICATOR || COMPUTEDRHODP|| SHOCKINDICATOR ){
1132  if((interporflux==ENOINTERPTYPE)&&(realisinterp)){
1133  // means primreal=p2interpm
1134  // then get all primitives and assign to normal yin
1135  NUMPRIMLOOP(pliter,pl) get_1d_line(dir, interporflux, pl, bs, ps, pe, be, i, j, k, primreal, NULL, yin[pl],&trueijkp);
1136  yprim=yin;
1137  }
1138  else{
1139  // means primreal!=p2interpm
1140  // then need to get separate p2interp line
1141  NUMPRIMLOOP(pliter,pl) get_1d_line(dir, interporflux, pl, bs, ps, pe, be, i, j, k, p2interpm, p2interpp, yin[pl],&trueijkp);
1142 
1143  // get prim line
1144  if(whichreduce == WENO_REDUCE_TYPE_PPM || SHOCKINDICATOR ){
1145  // then need to get primitives separately from interpolated quantities
1146  PALLREALLOOP(pl) get_1d_line(dir, interporflux, pl, bs, ps, pe, be, i, j, k, primreal, NULL, yprim[pl],&trueijkp);
1147  }
1148  else if(CONTACTINDICATOR||COMPUTEDRHODP){
1149  // then only need RHO and UU
1150  get_1d_line(dir, interporflux, RHO, bs, ps, pe, be, i, j, k, primreal, NULL, yprim[RHO],&trueijkp);
1151  get_1d_line(dir, interporflux, UU, bs, ps, pe, be, i, j, k, primreal, NULL, yprim[UU],&trueijkp);
1152  }
1153  }
1154  }
1155  else{
1156  // then no ppm reduce or contact indicator
1158  //
1159  // 1D GET LINE (interpolated quantity is never the same as the prims4shock quantity since not doing c2e here)
1160  //
1162  NUMPRIMLOOP(pliter,pl) get_1d_line(dir, interporflux, pl, bs, ps, pe, be, i, j, k, p2interpm, p2interpp, yin[pl],&trueijkp);
1163  }
1164 
1165 
1166 
1168  //
1169  // 1D GET LINE of stencilvar
1170  //
1172  if(!stencilvarisnull){
1173  NUMPRIMLOOP(pliter,pl) get_1d_line(dir, interporflux, pl, bs, ps, pe, be, i, j, k, stencilvar, stencilvar, ystencilvar[pl],&trueijkp);
1174  }
1175  else{
1176  // assign rather than pointer assign since yin//ystencilvar may be modified and don't want to have to worry if was
1177  NUMPRIMLOOP(pliter,pl) for(yiter=bs;yiter<=be;yiter++){
1178  ystencilvar[pl][0][yiter]=yin[pl][0][yiter];
1179  ystencilvar[pl][1][yiter]=yin[pl][1][yiter];
1180  }
1181  }
1182 
1183 
1184 
1185  // user-defined modification of the line -- usually adjusting line data so interpolation can be higher-order
1186  apply_bc_line(nprlocalstart,nprlocalend,nprlocallist,0,trueijkp.iter,recontype,bs,be,yin,NULL,NULL);
1187 
1188 
1190  //
1191  // get df's for contact indicator or in general
1192  // GODMARK: not done for flux splitting method (i.e. assume only 1 input always)
1193  // gets Vline and Pline now too if doingweno==1
1194  //
1196  get_df_line_gen_new(realisinterp, doingweno,whichprimtype,interporflux,recontype,dir,whichreduce,preforder, bs, ps, pe, be, minorder, maxorder, shift, yprim, yin, df, &drho, &dP,stiffindicator,Vline,Pline,&trueijkp);
1197 
1198 
1199 
1201  //
1202  // 1D GET LINE of mono df's if needed
1203  //
1205  if(!stencilvarisnull){
1206  // notice only requesting yin[0], not yin[1], so not good for FLUXSPLIT-- GODMARK
1207  NUMPRIMLOOP(pliter,pl) compute_df_line_formono(doingweno,interporflux, recontype, whichreduce, preforder, pl, bs, ps, pe, be, minorder, maxorder, shift, ystencilvar[pl][0], dfformono[pl]);
1208  }
1209  else{
1210  // assign rather than pointer assign since yin//ystencilvar may be modified and don't want to have to worry if was
1211  // should be expensive
1212  NUMPRIMLOOP(pliter,pl) for(yiter=bs;yiter<=be;yiter++) for(dfiter=0;dfiter<NUMDFS;dfiter++) dfformono[pl][dfiter][yiter]=df[pl][dfiter][yiter];
1213  }
1214 
1215 
1216 
1218  //
1219  // 1D P and V
1220  //
1222  if(!doingweno){ // need this because get_df_line_gen_new() didn't get Vline,Pline if doingweno==0
1223  if(SHOCKINDICATOR || CONTACTINDICATOR){
1224  get_V_and_P(whichprimtype, interporflux, dir, bs, ps, pe, be, i, j, k, idel, jdel, kdel, yprim, Vline, Pline,&trueijkp);
1225  }
1226  }
1227 
1229  //
1230  // 1D GET SHOCK INDICATOR
1231  //
1233 
1234  if(whichreduce == WENO_REDUCE_TYPE_PPM|| SHOCKINDICATOR ){
1235  //use primitive values that correspond to the quantities being interpolated
1236  get_shock_indicator(whichprimtype, interporflux, dir, bs, ps, pe, be, i, j, k, idel, jdel ,kdel, yprim, Vline, Pline, shockindicator,&trueijkp);
1237  }
1238 
1240  //
1241  // 1D GET CONTACT INDICATOR
1242  //
1244  if(CONTACTINDICATOR){
1245  //use primitive values that correspond to the quantities being interpolated
1246  get_contact_indicator(realisinterp, whichprimtype, interporflux, dir, bs, ps, pe, be, i, j, k, idel, jdel ,kdel, yin, Vline, Pline, etai);
1247  }
1248 
1249 
1250 
1251 
1253  //
1254  // Compute monotonicity indicator and if monotonic or rough then compute interface value
1255  // GODMARK: not done for flux splitting method (i.e. assume only 1 input always)
1256  //
1258 
1259  if(PARALINEUSESMONO||doingweno){ // let paraline use mono
1260  // note dfformono used, not df
1261  compute_monotonicity_line_multipl(stencilvarisnull,weightsplittype,whichquantity,dir,ALL_CALC,recontype, whichreduce, preforder, bs, ps, pe, be, minorder, maxorder, shift, shockindicator, stiffindicator, Vline, Pline, dfformono, dP, etai, monoindicator, yprim, ystencilvar, yin, yout, youtpolycoef);
1262  }
1263 
1264 
1265 
1266 
1267 #if(0)
1268 
1269  if(1||crapdebug){
1270  NUMPRIMLOOP(pliter,pl){
1271  if(recontype==CVT_C2E && (monoindicator[pl][0][0]!=monoindicator[pl][0][0+N1*idel+N2*jdel+N3*kdel] || monoindicator[pl][1][0]!=monoindicator[pl][1][0+N1*idel+N2*jdel+N3*kdel] || monoindicator[pl][2][0]!=monoindicator[pl][2][0+N1*idel+N2*jdel+N3*kdel])
1272  ||
1273  (monoindicator[pl][0][0]!=monoindicator[pl][0][0+N1*idel+N2*jdel+N3*kdel] || monoindicator[pl][1][0]!=monoindicator[pl][1][0+N1*idel+N2*jdel+N3*kdel])
1274  ){
1275  dualfprintf(fail_file,"Detected asymmetry: recontype=%d dir=%d ps=%d pe=%d\n",recontype,dir,ps,pe);
1276  dualfprintf(fail_file,"i=%d j=%d k=%d di=%d dj=%d dk=%d\n",i,j,k,di,dj,dk);
1277  dualfprintf(fail_file,"%21.15g %21.15g %21.15g %21.15g %21.15g %21.15g\n",monoindicator[pl][0][0],monoindicator[pl][0][0+N1*idel+N2*jdel+N3*kdel], monoindicator[pl][1][0],monoindicator[pl][1][0+N1*idel+N2*jdel+N3*kdel],monoindicator[pl][2][0],monoindicator[pl][2][0+N1*idel+N2*jdel+N3*kdel]);
1278  }
1279  }
1280  }
1281 #endif
1282 
1283 
1284 
1285 
1286 
1287 
1289  //
1290  // PASS 1D LINE
1291  //
1293 
1294  // assume default behavior is to send all pl's
1295  pass_1d_line_multipl(weightsplittype, whichquantity, dir, ALL_CALC, recontype, whichreduce, preforder, bs, ps, pe, be, minorder, maxorder, shift, shockindicator, stiffindicator, Vline, Pline, df, dP, etai, monoindicator, yprim, ystencilvar, yin, yout, youtpolycoef,&trueijkp);
1296 
1297 
1298  // user-defined de-modification of the line -- usually adjusting line data so interpolation can be higher-order
1299  apply_bc_line(nprlocalstart,nprlocalend,nprlocallist,1,trueijkp.iter,recontype,bs,be,yin,yout,youtpolycoef);
1300 
1301 
1302  }// end if doing normal variation of weights
1303 
1304 
1305 
1306 
1307 
1309  //
1310  // now have 1D line of point data corresponding to left and right interpolated values
1311  //
1312  // Assign result back to pleft/pright (or just pleft if one result)
1313  //
1315  NUMPRIMLOOP(pliter,pl){
1316 #if(0)
1317  if(crapdebug){
1318  dualfprintf(fail_file,"pl=%d\n",pl);
1319  }
1320  // if(crapdebug==0 && pl!=3) assign_eno_result(recontype, pl, bs, ps, pe, be, i, j, k, idel, jdel, kdel, yout[pl], pleft, pright);
1321  //assign_eno_result(recontype, pl, bs, ps, pe, be, i, j, k, idel, jdel, kdel, yout[pl], pleft, pright);
1322  if(crapdebug==0 || crapdebug==1 || pl!=3 ) assign_eno_result(recontype, pl, bs, ps, pe, be, i, j, k, idel, jdel, kdel, yout[pl], pleft, pright);
1323  else dualfprintf(fail_file,"crapdebug=%d\n",crapdebug);
1324 #endif
1325 
1326  assign_eno_result(recontype, pl, bs, ps, pe, be, i, j, k, idel, jdel, kdel, yout[pl], pleft, pright);
1327 
1328  }
1329 
1330 
1331 
1332 }
1333 
1334 
1335 
1336 
1337 
1338 
1339 
1340 
1341 
1342 
1344 static int get_recon_type(int interporflux)
1345 {
1346 
1347  // c2a stuff
1348  if(interporflux==ENOFLUXAVG1TYPE|| interporflux==ENOFLUXAVG2TYPE || interporflux==ENOFLUXAVG3TYPE || interporflux==ENOCENT2AVGTYPE || interporflux==ENOQUASIFIELDFLUXRECONTYPE) return(CVT_C2A);
1349  // c2e stuff
1350  else if(interporflux==ENOFLUXSPLITTYPE || interporflux==ENOINTERPTYPE || interporflux==ENOINTERPTYPE4EMF) return(CVT_C2E);
1351  // a2c stuff
1352  else if(interporflux==ENOFLUXRECONTYPE || interporflux==ENOFLUXRECONTYPEGHOSTACTIVE || interporflux==ENOAVG2CENTTYPE) return(CVT_A2C);
1353  else{
1354  dualfprintf(fail_file,"No such interporflux=%d\n", interporflux);
1355  myexit(177);
1356  }
1357 
1358  return(-100);
1359 
1360 }
1361 
1362 
1363 
1365 static void set_preforder(int dir, int interporflux, int *preforder, int*whichreduce)
1366 {
1367 
1368  if(interporflux==ENOINTERPTYPE || interporflux==ENOINTERPTYPE4EMF){
1369  // reallim=choose_limiter(dir,i,j,k,pl); // starting point chooses limiter type
1370  *preforder=interporder[lim[dir]]; // get order of scheme
1371 
1372  if( lim[dir] == WENO5FLAT ) { //correct the order for WENO5FLAT: in this case need more points than the order for stencil reduction
1373  *whichreduce = WENO_REDUCE_TYPE_PPM;
1374  *preforder = 5;
1375  }
1376  else if( WENOBNDPINTERPTYPE(lim[dir]) ) { //correct the order for WENO5FLAT: in this case need more points than the order for stencil reduction
1377  *whichreduce = WENO_REDUCE_TYPE_DEFAULT;
1378  *preforder = 5; //correct the order of the scheme because the number of points passed to it is larger than its order (the order is 5)
1379  }
1380  else if( lim[dir] == PARALINE ) { // true order
1381  *whichreduce = WENO_REDUCE_TYPE_DEFAULT;
1382  *preforder = 5; //correct the order of the scheme because the number of points passed to it is larger than its order (the order is 5)
1383  }
1384  else {
1385  *whichreduce = WENO_REDUCE_TYPE_DEFAULT;
1386  }
1387  }
1388  else{
1389  *preforder=interporder[avgscheme[dir]];
1390 
1391  if( avgscheme[dir] == WENO5FLAT ) { //correct the order for WENO5FLAT: in this case need more points than the order for stencil reduction
1392  *whichreduce = WENO_REDUCE_TYPE_PPM;
1393  *preforder -= 2;
1394  }
1395  else if( WENOBNDPINTERPTYPE(avgscheme[dir]) ) { //correct the order for WENO5FLAT: in this case need more points than the order for stencil reduction
1396  *whichreduce = WENO_REDUCE_TYPE_DEFAULT;
1397  *preforder = 5; //correct the order of the scheme because the number of points passed to it is larger than its order (the order is 5)
1398  }
1399  else {
1400  *whichreduce = WENO_REDUCE_TYPE_DEFAULT;
1401  }
1402  }
1403 
1404 
1405 }
1406 
1407 
1408 
1409 
1410 
1411 
1412 
1413 
1416 static void get_df_line_gen(int realisinterp, int doingweno, int whichprimtype, int interporflux, int recontype, int dir, int whichreduce, int preforder, int bs, int ps, int pe, int be, int *minorder, int *maxorder, int *shift, FTYPE (*yprim)[2][NBIGM], FTYPE (*yin)[2][NBIGM], FTYPE (*df)[NUMDFS][NBIGM], FTYPE (**drhoptr)[NBIGM], FTYPE (**dPptr)[NBIGM], FTYPE *stiffindicator)
1417 {
1418  int pl,pliter;
1419  int nprlocalstart,nprlocalend;
1420  int nprlocallist[MAXNPR];
1421  int pllocal;
1422  int numprims;
1423 
1424 
1426  //
1427  // Define which quantities (pl) to operate on
1428  //
1430 
1431  setup_nprlocalist(whichprimtype,&nprlocalstart,&nprlocalend,nprlocallist,&numprims);
1432 
1433 
1434 
1435  if(CONTACTINDICATOR||COMPUTEDRHODP){
1436  if((interporflux==ENOINTERPTYPE)&&(realisinterp)){
1437  // then don't need to get separate drho and dP
1438 
1439  // then get all df's
1440  NUMPRIMLOOP(pliter,pl) compute_df_line(doingweno,interporflux,recontype, whichreduce,preforder, pl, bs, ps, pe, be, minorder, maxorder, shift, yin[pl][0], df[pl]);
1441 
1442  // then primitives are same as interpolated quantities
1443  // assign drho and dP. This overwrites original pointer that pointed to independent memory for drho and dP
1444  *drhoptr=df[RHO];
1445  *dPptr=df[UU];
1446  }
1447  else{
1448  // then need to compute separate drho and dP
1450  // Compute differentials needed for contactindicator
1452  compute_df_line(doingweno,interporflux,recontype, whichreduce,preforder, RHO, bs, ps, pe, be, minorder, maxorder, shift, yprim[RHO][0], *drhoptr);
1453  compute_df_line(doingweno,interporflux,recontype, whichreduce,preforder, UU, bs, ps, pe, be, minorder, maxorder, shift, yprim[UU][0], *dPptr);
1454 
1455  // then get all p2interp df's
1456  NUMPRIMLOOP(pliter,pl) compute_df_line(doingweno,interporflux,recontype, whichreduce,preforder, pl, bs, ps, pe, be, minorder, maxorder, shift, yin[pl][0], df[pl]);
1457  // if(interporflux==ENOFLUXSPLITTYPE) NUMPRIMLOOP(pliter,pl) compute_df_line(doingweno,interporflux,recontype, whichreduce,preforder, pl, bs, ps, pe, be, minorder, maxorder, shift, yin[pl][1], dfp[pl]);
1458  }
1459  }
1460  else{
1461  // then don't need drho or dP at all, just get p2interp df's
1462  // then get all p2interp df's
1463  NUMPRIMLOOP(pliter,pl) compute_df_line(doingweno,interporflux,recontype, whichreduce,preforder, pl, bs, ps, pe, be, minorder, maxorder, shift, yin[pl][0], df[pl]);
1464  // if(interporflux==ENOFLUXSPLITTYPE) NUMPRIMLOOP(pliter,pl) compute_df_line(doingweno,interporflux,recontype, whichreduce,preforder, pl, bs, ps, pe, be, minorder, maxorder, shift, yin[pl][1], dfp[pl]);
1465  }
1466 }
1467 
1468 
1469 
1470 
1471 
1472 
1473 
1475 static void get_df_line_gen_new(int realisinterp, int doingweno, int whichprimtype, int interporflux, int recontype, int dir, int whichreduce, int preforder, int bs, int ps, int pe, int be, int *minorder, int *maxorder, int *shift, FTYPE (*yprim)[2][NBIGM], FTYPE (*yin)[2][NBIGM], FTYPE (*df)[NUMDFS][NBIGM], FTYPE (**drhoptr)[NBIGM], FTYPE (**dPptr)[NBIGM], FTYPE *stiffindicator, FTYPE (*Vline)[NBIGM], FTYPE (*Pline)[NBIGM], struct of_trueijkp *trueijkp)
1476 {
1477  int pl,pliter;
1478  int nprlocalstart,nprlocalend;
1479  int nprlocallist[MAXNPR];
1480  int pllocal;
1481  int numprims;
1482 
1483 
1484 
1485 
1487  //
1488  // Define which quantities (pl) to operate on
1489  //
1491 
1492  setup_nprlocalist(whichprimtype,&nprlocalstart,&nprlocalend,nprlocallist,&numprims);
1493 
1494 
1495 
1496 
1497 
1498  if(SHOCKINDICATOR || CONTACTINDICATOR||COMPUTEDRHODP){
1499  if((interporflux==ENOINTERPTYPE)&&(realisinterp)){
1500  // then don't need to get separate drho and dP
1501 
1502 
1503  // then get all df's
1504  NUMPRIMLOOP(pliter,pl) compute_df_line(doingweno,interporflux,recontype,whichreduce,preforder, pl, bs, ps, pe, be, minorder, maxorder, shift, yin[pl][0], df[pl]);
1505 
1506  // weno uses Lorentz factor-based differences (gets stiffindicator as well)
1507  if(doingweno) compute_df_line_new(doingweno, whichprimtype, interporflux, recontype, dir, whichreduce,preforder, bs, ps, pe, be, minorder, maxorder, shift, yprim, (*dPptr)[0], stiffindicator, Vline, Pline,trueijkp);
1508 
1509  // then primitives are same as interpolated quantities
1510  // assign drho and dP. This overwrites original pointer that pointed to independent memory for drho and dP
1511  *drhoptr=df[RHO];
1512  // *dPptr=df[UU];
1513  }
1514  else{
1515  // then need to compute separate drho and dP
1517  // Compute differentials needed for contactindicator
1519 
1520  // then get all p2interp df's
1521  NUMPRIMLOOP(pliter,pl) compute_df_line(doingweno,interporflux,recontype,whichreduce,preforder, pl, bs, ps, pe, be, minorder, maxorder, shift, yin[pl][0], df[pl]);
1522  // if(interporflux==ENOFLUXSPLITTYPE) NUMPRIMLOOP(pliter,pl) compute_df_line(doingweno,interporflux,recontype,whichreduce,preforder, pl, bs, ps, pe, be, minorder, maxorder, shift, yin[pl][1], dfp[pl]);
1523 
1524  compute_df_line(doingweno,interporflux,recontype,whichreduce,preforder, RHO, bs, ps, pe, be, minorder, maxorder, shift, yprim[RHO][0], *drhoptr);
1525 
1526  if(doingweno) compute_df_line_new(doingweno,whichprimtype,interporflux,recontype,dir, whichreduce,preforder, bs, ps, pe, be, minorder, maxorder, shift, yprim, (*dPptr)[0], stiffindicator, Vline, Pline,trueijkp);
1527 
1528  }
1529  }
1530  else{
1531  // then don't need drho or dP at all, just get p2interp df's
1532  // then get all p2interp df's
1533  NUMPRIMLOOP(pliter,pl) compute_df_line(doingweno,interporflux,recontype,whichreduce,preforder, pl, bs, ps, pe, be, minorder, maxorder, shift, yin[pl][0], df[pl]);
1534  // if(interporflux==ENOFLUXSPLITTYPE) NUMPRIMLOOP(pliter,pl) compute_df_line(doingweno,interporflux,recontype,whichreduce,preforder, pl, bs, ps, pe, be, minorder, maxorder, shift, yin[pl][1], dfp[pl]);
1535  }
1536 }
1537 
1538 
1539 
1541 static void get_df_line_paraline(int realisinterp, int doingweno, int whichprimtype, int interporflux, int recontype, int dir, int whichreduce, int preforder, int bs, int ps, int pe, int be, int *minorder, int *maxorder, int *shift, FTYPE (*yprim)[2][NBIGM], FTYPE (*yin)[2][NBIGM], FTYPE (*df)[NUMDFS][NBIGM], FTYPE (**drhoptr)[NBIGM], FTYPE (**dPptr)[NBIGM], FTYPE *stiffindicator, FTYPE (*Vline)[NBIGM], FTYPE (*Pline)[NBIGM], struct of_trueijkp *trueijkp)
1542 {
1543  int pl,pliter;
1544  int nprlocalstart,nprlocalend;
1545  int nprlocallist[MAXNPR];
1546  int pllocal;
1547  int numprims;
1548 
1549 
1550 
1551 
1553  //
1554  // Define which quantities (pl) to operate on
1555  //
1557 
1558  setup_nprlocalist(whichprimtype,&nprlocalstart,&nprlocalend,nprlocallist,&numprims);
1559 
1560 
1561 
1562  // then don't need drho or dP at all, just get p2interp df's
1563  // then get all p2interp df's
1564  NUMPRIMLOOP(pliter,pl) compute_df_line_paraline(doingweno,interporflux,recontype,whichreduce,preforder, pl, bs, ps, pe, be, minorder, maxorder, shift, yin[pl][0], df[pl]);
1565  // if(interporflux==ENOFLUXSPLITTYPE) NUMPRIMLOOP(pliter,pl) compute_df_line(doingweno,interporflux,recontype,whichreduce,preforder, pl, bs, ps, pe, be, minorder, maxorder, shift, yin[pl][1], dfp[pl]);
1566 
1567 
1568 }
1569 
1570 
1571 
1572 
1573 
1575 void set_interp_loop_gen(int withshifts, int interporflux, int dir, int loc, int continuous, int *intdir, int *is, int *ie, int *js, int *je, int *ks, int *ke, int *di, int *dj, int *dk, int *bs, int *ps, int *pe, int *be)
1576 {
1577 
1578  if(useghostplusactive){
1579  // now assume active + active/ghost + ghost layers so only bound primitives during calculation
1580  set_interp_loop_expanded(withshifts, interporflux, dir, loc, continuous, intdir, is, ie, js, je, ks, ke, di, dj, dk, bs, ps, pe, be);
1581  }
1582  else{
1583  // straight-forward average or de-average along dir and just one ghost layer (original method)
1584  // GODMARK: later should convert ENOFLUXRECON method to have expanded ghost layer and ghost+active layer
1585  set_interp_loop(withshifts, interporflux, dir, loc, continuous, intdir, is, ie, js, je, ks, ke, di, dj, dk, bs, ps, pe, be);
1586  }
1587 
1588 }
1589 
1590 
1591 
1596 static void set_interp_loop(int withshifts, int interporflux, int dir, int loc, int continuous, int *intdir, int *is, int *ie, int *js, int *je, int *ks, int *ke, int *di, int *dj, int *dk, int *bs, int *ps, int *pe, int *be)
1597 {
1598 
1600  //
1601  //
1602  // DIRECTION OF LINE EXTRACTION and INTERPOLATION = 1
1603  //
1604  //
1606  // determine range of outer loop and range to feed to eno scheme
1607  if(dir==1){
1608  *intdir=dir;
1609 
1610  *is=*ie=0; // anything so di=1 iterates, next subloop overwrites i actually useds
1611 
1612  if((interporflux==ENOFLUXSPLITTYPE)||(interporflux==ENOINTERPTYPE || interporflux==ENOINTERPTYPE4EMF)){
1613  // then centered quantity, and need from -1 .. N
1614  *bs=INFULL1;
1615  *be=OUTFULL1;
1616 
1617  *ps=0 -SHIFT1;
1618  *pe=N1-1 +SHIFT1;
1619  }
1620  else if(interporflux==ENOFLUXRECONTYPE || interporflux==ENOFLUXRECONTYPEGHOSTACTIVE || interporflux==ENOQUASIFIELDFLUXRECONTYPE){
1621  // then edge quantity and need from 0 .. N
1622  *bs=INFULL1;
1623  *be=OUTFULL1;
1624 
1625  *ps=0;
1626  *pe=N1-1+SHIFT1;
1627 
1628  if(interporflux==ENOFLUXRECONTYPEGHOSTACTIVE){
1629  *ps -= (interporder[avgscheme[dir]]-1)/2;
1630  *pe += (interporder[avgscheme[dir]]-1)/2;
1631  }
1632 
1633 
1634  }
1635 
1636  *js=fluxloop[dir][FJS]; // 0
1637  *je=fluxloop[dir][FJE]; // N2-1;
1638 
1639  *ks=fluxloop[dir][FKS]; //0;
1640  *ke=fluxloop[dir][FKE]; // N3-1;
1641 
1642  *di=1;
1643  *dj=1;
1644  *dk=1;
1645 
1646  if(withshifts){
1647  // shift due to grid sectionion SECTIONMARK
1648  *bs += SHIFTX1DN;
1649  *be += SHIFTX1UP;
1650  *ps += SHIFTX1DN;
1651  *pe += SHIFTX1UP;
1652  *js += SHIFTX2DN;
1653  *je += SHIFTX2UP;
1654  *ks += SHIFTX3DN;
1655  *ke += SHIFTX3UP;
1656  }
1657 
1658  }
1660  //
1661  //
1662  // DIRECTION OF LINE EXTRACTION and INTERPOLATION = 2
1663  //
1664  //
1666  else if(dir==2){
1667  *intdir=dir;
1668 
1669  *is=fluxloop[dir][FIS]; //0;
1670  *ie=fluxloop[dir][FIE]; //N1-1;
1671 
1672 
1673  *js=*je=0; // anything so dj=1 iterates, next subloop overwrites j actually useds
1674 
1675  if((interporflux==ENOFLUXSPLITTYPE)||(interporflux==ENOINTERPTYPE) || interporflux==ENOINTERPTYPE4EMF){
1676  // then centered quantity, and need from -1 .. N
1677  *bs=INFULL2;
1678  *be=OUTFULL2;
1679 
1680  *ps=0 -SHIFT2;
1681  *pe=N2-1 +SHIFT2;
1682  }
1683  else if(interporflux==ENOFLUXRECONTYPE || interporflux==ENOFLUXRECONTYPEGHOSTACTIVE || interporflux==ENOQUASIFIELDFLUXRECONTYPE){
1684  // then edge quantity and need from 0 .. N
1685  *bs=INFULL2;
1686  *be=OUTFULL2;
1687 
1688  *ps=0;
1689  *pe=N2-1+SHIFT2;
1690 
1691  if(interporflux==ENOFLUXRECONTYPEGHOSTACTIVE){
1692  *ps -= (interporder[avgscheme[dir]]-1)/2;
1693  *pe += (interporder[avgscheme[dir]]-1)/2;
1694  }
1695 
1696  }
1697 
1698  *ks=fluxloop[dir][FKS]; // 0;
1699  *ke=fluxloop[dir][FKE]; // N3-1;
1700 
1701  *di=1;
1702  *dj=1;
1703  *dk=1;
1704 
1705  if(withshifts){
1706  // shift due to grid sectionion SECTIONMARK
1707  *is += SHIFTX1DN;
1708  *ie += SHIFTX1UP;
1709  *bs += SHIFTX2DN;
1710  *be += SHIFTX2UP;
1711  *ps += SHIFTX2DN;
1712  *pe += SHIFTX2UP;
1713  *ks += SHIFTX3DN;
1714  *ke += SHIFTX3UP;
1715  }
1716 
1717  }
1719  //
1720  //
1721  // DIRECTION OF LINE EXTRACTION and INTERPOLATION = 3
1722  //
1723  //
1725  else if(dir==3){
1726  *intdir=dir;
1727 
1728  *is=fluxloop[dir][FIS]; // 0;
1729  *ie=fluxloop[dir][FIE]; // N1-1;
1730 
1731  *js=fluxloop[dir][FJS]; // 0;
1732  *je=fluxloop[dir][FJE]; // N2-1;
1733 
1734  *ks=*ke=0; // anything so dk=1 iterates, next subloop overwrites k actually used
1735 
1736 
1737  if((interporflux==ENOFLUXSPLITTYPE)||(interporflux==ENOINTERPTYPE) || interporflux==ENOINTERPTYPE4EMF){
1738  // then centered quantity, and need from -1 .. N
1739  *bs=INFULL3;
1740  *be=OUTFULL3;
1741 
1742  *ps=0 -SHIFT3;
1743  *pe=N3-1 +SHIFT3;
1744  }
1745  else if(interporflux==ENOFLUXRECONTYPE || interporflux==ENOFLUXRECONTYPEGHOSTACTIVE || interporflux==ENOQUASIFIELDFLUXRECONTYPE){
1746  // then edge quantity and need from 0 .. N
1747  *bs=INFULL3;
1748  *be=OUTFULL3;
1749 
1750  *ps=0;
1751  *pe=N3-1+SHIFT3;
1752 
1753  if(interporflux==ENOFLUXRECONTYPEGHOSTACTIVE){
1754  *ps -= (interporder[avgscheme[dir]]-1)/2;
1755  *pe += (interporder[avgscheme[dir]]-1)/2;
1756  }
1757 
1758  }
1759 
1760  *di=1;
1761  *dj=1;
1762  *dk=1;
1763 
1764  if(withshifts){
1765  // shift due to grid sectionion SECTIONMARK
1766  *is += SHIFTX1DN;
1767  *ie += SHIFTX1UP;
1768  *js += SHIFTX2DN;
1769  *je += SHIFTX2UP;
1770  *bs += SHIFTX3DN;
1771  *be += SHIFTX3UP;
1772  *ps += SHIFTX3DN;
1773  *pe += SHIFTX3UP;
1774  }
1775 
1776 
1777  }
1778 
1779 }
1780 
1781 
1782 
1783 
1784 
1789 static void set_interp_loop_expanded(int withshifts, int interporflux, int dir, int loc, int continuous, int *intdir, int *is, int *ie, int *js, int *je, int *ks, int *ke, int *di, int *dj, int *dk, int *bs, int *ps, int *pe, int *be)
1790 {
1791  int dir_exception[NDIM];
1792  int *myUconsloop;
1793 
1794 
1795  if(interporflux==ENOINTERPTYPE4EMF){
1796  myUconsloop=emfUconsloop;
1797  }
1798  else{
1799  myUconsloop=Uconsloop;
1800  }
1801 
1802  // if((interporflux==ENOFLUXAVG1TYPE)||(interporflux==ENOFLUXAVG2TYPE)||(interporflux==ENOFLUXAVG3TYPE) ){
1803 
1804 
1805  // ENOFLUXAVG?TYPE is ?=direction of integration not direction of flux
1806  if(
1807  ((interporflux==ENOFLUXAVG1TYPE)&&(dir==1))||
1808  ((interporflux==ENOFLUXAVG2TYPE)&&(dir==2))||
1809  ((interporflux==ENOFLUXAVG3TYPE)&&(dir==3))
1810  ){
1811  dualfprintf(fail_file,"No such method with interporflux=%d and dir=%d\n",interporflux,dir);
1812  myexit(1);
1813  }
1814 
1815  dir_exception[1] = (interporflux==ENOFLUXAVG2TYPE) || (interporflux==ENOFLUXAVG3TYPE);
1816  dir_exception[2] = (interporflux==ENOFLUXAVG1TYPE) || (interporflux==ENOFLUXAVG3TYPE);
1817  dir_exception[3] = (interporflux==ENOFLUXAVG1TYPE) || (interporflux==ENOFLUXAVG2TYPE);
1818 
1819  // determine range of outer loop and range to feed to eno scheme
1820  // Effectively below corresponds to 3 cases:
1821  // For example, for dir==1:
1822  // 1) not an AVG?TYPE : enter if dir==1
1823  // 2) AVG<dir>type : enter
1824  // 3) AVG<other dir>type : do not enter
1825 
1827  //
1828  //
1829  // DIRECTION OF LINE EXTRACTION and INTERPOLATION = 1
1830  //
1831  //
1833  if( ( (!dir_exception[1]) && (dir==1) ) || (interporflux==ENOFLUXAVG1TYPE) ){
1834  *intdir=1; // not generally true that intdir=dir
1835 
1836  *is=*ie=0; // anything so di=1 iterates, next subloop overwrites i actually used
1837 
1838  // input and output at different location
1839 
1840  if((interporflux==ENOFLUXSPLITTYPE)||(interporflux==ENOINTERPTYPE) || interporflux==ENOINTERPTYPE4EMF){
1841  // then centered quantity, and need from -1 .. N
1842  *bs=INFULL1;
1843  *be=OUTFULL1;
1844 
1845  *ps=myUconsloop[FIS]-SHIFT1;
1846  *pe=myUconsloop[FIE]+SHIFT1;
1847 
1848  *js=fluxloop[dir][FJS];
1849  *je=fluxloop[dir][FJE];
1850 
1851  *ks=fluxloop[dir][FKS];
1852  *ke=fluxloop[dir][FKE];
1853  }
1854  else if(interporflux==ENOFLUXRECONTYPE || interporflux==ENOFLUXRECONTYPEGHOSTACTIVE || interporflux==ENOQUASIFIELDFLUXRECONTYPE){
1855  // input and output at same location
1856 
1857  // then edge quantity and need from 0 .. N
1858  *bs=INFULL1;
1859  *be=OUTFULL1;
1860 
1861  // *ps=myUconsloop[FIS];
1862  // *pe=myUconsloop[FIE]+1;
1863  *ps=0;
1864  *pe=N1-1+SHIFT1;
1865 
1866  if(interporflux==ENOFLUXRECONTYPEGHOSTACTIVE){
1867  *ps -= (interporder[avgscheme[dir]]-1)/2;
1868  *pe += (interporder[avgscheme[dir]]-1)/2;
1869  }
1870 
1871 
1872  *js=fluxloop[dir][FJS];
1873  *je=fluxloop[dir][FJE];
1874 
1875  *ks=fluxloop[dir][FKS];
1876  *ke=fluxloop[dir][FKE];
1877 
1878  }
1879  else if( interporflux==ENOAVG2CENTTYPE ){
1880  // input and output at same location
1881 
1882  // then edge quantity and need from 0 .. N
1883  *bs=myUconsloop[FIS];
1884  *be=myUconsloop[FIE];
1885 
1886  *ps=0;
1887  *pe=N1-1;
1888 
1889  *js=fluxloop[dir][FJS];
1890  *je=fluxloop[dir][FJE];
1891 
1892  *ks=fluxloop[dir][FKS];
1893  *ke=fluxloop[dir][FKE];
1894 
1895 
1896  }
1897  else if( interporflux==ENOCENT2AVGTYPE ){
1898  *bs=INFULL1;
1899  *be=OUTFULL1;
1900 
1901  *ps=myUconsloop[FIS];
1902  *pe=myUconsloop[FIE];
1903 
1904  *js=fluxloop[dir][FJS];
1905  *je=fluxloop[dir][FJE];
1906 
1907  *ks=fluxloop[dir][FKS];
1908  *ke=fluxloop[dir][FKE];
1909 
1910  }
1911  else if(interporflux==ENOFLUXAVG1TYPE){
1912 
1913  // then edge quantity and need from 0 .. N and have ghost zones out to ijkminmax[1][0]+1
1914  *bs=INFULL1;
1915  *be=OUTFULL1;
1916 
1917  *ps=myUconsloop[FIS];
1918  *pe=myUconsloop[FIE];
1919 
1920  if(dir==2){
1921  // This method doesn't need "boundary" direction of "fluxes" since they don't actually exist
1922  *js=myUconsloop[FJS];
1923  *je=myUconsloop[FJE]+SHIFT2;
1924  }
1925  else{
1926  // GODMARK: needed?
1927  *js=myUconsloop[FJS];
1928  *je=myUconsloop[FJE];
1929  // *js=fluxloop[dir][FJS];
1930  // *je=fluxloop[dir][FJE];
1931  }
1932 
1933  if(dir==3){
1934  // This method doesn't need "boundary" direction of "fluxes" since they don't actually exist
1935  *ks=myUconsloop[FKS];
1936  *ke=myUconsloop[FKE]+SHIFT3;
1937  }
1938  else{
1939  // GODMARK: needed?
1940  *ks=myUconsloop[FKS];
1941  *ke=myUconsloop[FKE];
1942  //*ks=fluxloop[dir][FKS];
1943  // *ke=fluxloop[dir][FKE];
1944  }
1945 
1946  }
1947 
1948 
1949  *di=1;
1950  *dj=1;
1951  *dk=1;
1952 
1953  if(withshifts){
1954  // shift due to grid sectionion SECTIONMARK
1955  *bs += SHIFTX1DN;
1956  *be += SHIFTX1UP;
1957  *ps += SHIFTX1DN;
1958  *pe += SHIFTX1UP;
1959  *js += SHIFTX2DN;
1960  *je += SHIFTX2UP;
1961  *ks += SHIFTX3DN;
1962  *ke += SHIFTX3UP;
1963  }
1964 
1965  }
1967  //
1968  //
1969  // DIRECTION OF LINE EXTRACTION and INTERPOLATION = 2
1970  //
1971  //
1973  else if( ( (!dir_exception[2]) && (dir==2) ) || (interporflux==ENOFLUXAVG2TYPE) ){
1974  *intdir=2; // not generally true that intdir=dir
1975 
1976  *js=*je=0; // anything so dj=1 iterates, next subloop overwrites j actually used
1977 
1978  if((interporflux==ENOFLUXSPLITTYPE)||(interporflux==ENOINTERPTYPE) || interporflux==ENOINTERPTYPE4EMF){
1979  // then centered quantity, and need from -1 .. N
1980  *bs=INFULL2;
1981  *be=OUTFULL2;
1982 
1983  *ps=myUconsloop[FJS]-SHIFT2;
1984  *pe=myUconsloop[FJE]+SHIFT2;
1985 
1986  *is=INFULL1;
1987  *ie=OUTFULL1;
1988 
1989  *ks=fluxloop[dir][FKS];
1990  *ke=fluxloop[dir][FKE];
1991  }
1992  else if(interporflux==ENOFLUXRECONTYPE || interporflux==ENOFLUXRECONTYPEGHOSTACTIVE || interporflux==ENOQUASIFIELDFLUXRECONTYPE){
1993  // then edge quantity and need from 0 .. N
1994  *bs=INFULL2;
1995  *be=OUTFULL2;
1996 
1997  // *ps=myUconsloop[FJS];
1998  // *pe=myUconsloop[FJE]+1;
1999  *ps=0;
2000  *pe=N2-1+SHIFT2;
2001 
2002  if(interporflux==ENOFLUXRECONTYPEGHOSTACTIVE){
2003  *ps -= (interporder[avgscheme[dir]]-1)/2;
2004  *pe += (interporder[avgscheme[dir]]-1)/2;
2005  }
2006 
2007 
2008  *is=INFULL1;
2009  *ie=OUTFULL1;
2010 
2011  *ks=fluxloop[dir][FKS];
2012  *ke=fluxloop[dir][FKE];
2013  }
2014  else if( interporflux==ENOAVG2CENTTYPE ){
2015  // input and output at same location
2016 
2017  // then edge quantity and need from 0 .. N
2018  *bs=myUconsloop[FJS]; //atch correct SASMARK; ALSO need to set up the same thing for other dimensions-- not sure if this is enough
2019  *be=myUconsloop[FJE];
2020 
2021  *ps=0;
2022  *pe=N2 - 1;
2023 
2024  *is=INFULL1;
2025  *ie=OUTFULL1;
2026 
2027  *ks=fluxloop[dir][FKS];
2028  *ke=fluxloop[dir][FKE];
2029  }
2030  else if(interporflux==ENOCENT2AVGTYPE){
2031  // then edge quantity and need from 0 .. N
2032  *bs=INFULL2;
2033  *be=OUTFULL2;
2034 
2035  *ps=myUconsloop[FJS];
2036  *pe=myUconsloop[FJE];
2037 
2038  *is=INFULL1;
2039  *ie=OUTFULL1;
2040 
2041  *ks=fluxloop[dir][FKS];
2042  *ke=fluxloop[dir][FKE];
2043  }
2044  else if(interporflux==ENOFLUXAVG2TYPE){
2045  //is bs and be to be corrected for the shock indicator? SASMARK
2046  // then edge quantity and need from 0 .. N
2047  *bs=INFULL2;
2048  *be=OUTFULL2;
2049 
2050  *ps=myUconsloop[FJS];
2051  *pe=myUconsloop[FJE];
2052 
2053  if(dir==1){
2054  // This method doesn't need "boundary" direction of "fluxes" since they don't actually exist
2055  *is=myUconsloop[FIS];
2056  *ie=myUconsloop[FIE]+SHIFT1;
2057  }
2058  else{
2059  *is=myUconsloop[FIS];
2060  *ie=myUconsloop[FIE];
2061  // GODMARK: needed?
2062  // *is=INFULL1;
2063  // *ie=OUTFULL1;
2064  }
2065 
2066  if(dir==3){
2067  // This method doesn't need "boundary" direction of "fluxes" since they don't actually exist
2068  *ks=myUconsloop[FKS];
2069  *ke=myUconsloop[FKE]+SHIFT3;
2070  }
2071  else{
2072  // GODMARK: needed?
2073  *ks=myUconsloop[FKS];
2074  *ke=myUconsloop[FKE];
2075  // *ks=fluxloop[dir][FKS];
2076  // *ke=fluxloop[dir][FKE];
2077  }
2078 
2079  }
2080 
2081 
2082  *di=1;
2083  *dj=1;
2084  *dk=1;
2085 
2086  if(withshifts){
2087  // shift due to grid sectionion SECTIONMARK
2088  *is += SHIFTX1DN;
2089  *ie += SHIFTX1UP;
2090  *bs += SHIFTX2DN;
2091  *be += SHIFTX2UP;
2092  *ps += SHIFTX2DN;
2093  *pe += SHIFTX2UP;
2094  *ks += SHIFTX3DN;
2095  *ke += SHIFTX3UP;
2096  }
2097 
2098  }
2100  //
2101  //
2102  // DIRECTION OF LINE EXTRACTION and INTERPOLATION = 3
2103  //
2104  //
2106  else if( ( (!dir_exception[3]) && (dir==3) ) || (interporflux==ENOFLUXAVG3TYPE) ){
2107  *intdir=3; // not generally true that intdir=dir
2108 
2109 
2110  *ks=*ke=0; // anything so dk=1 iterates, next subloop overwrites k actually used
2111 
2112 
2113  if((interporflux==ENOFLUXSPLITTYPE)||(interporflux==ENOINTERPTYPE) || interporflux==ENOINTERPTYPE4EMF){
2114  // then centered quantity, and need from -1 .. N
2115  *bs=INFULL3;
2116  *be=OUTFULL3;
2117 
2118  *ps=myUconsloop[FKS]-SHIFT3;
2119  *pe=myUconsloop[FKE]+SHIFT3;
2120 
2121  *is=INFULL1;
2122  *ie=OUTFULL1;
2123 
2124  *js=fluxloop[dir][FJS];
2125  *je=fluxloop[dir][FJE];
2126 
2127  }
2128  else if(interporflux==ENOFLUXRECONTYPE || interporflux==ENOFLUXRECONTYPEGHOSTACTIVE || interporflux==ENOQUASIFIELDFLUXRECONTYPE){
2129  // then edge quantity and need from 0 .. N
2130  *bs=INFULL3;
2131  *be=OUTFULL3;
2132 
2133  // *ps=myUconsloop[FKS];
2134  // *pe=myUconsloop[FKE]+1;
2135  *ps=0;
2136  *pe=N3-1+SHIFT3;
2137 
2138  if(interporflux==ENOFLUXRECONTYPEGHOSTACTIVE){
2139  *ps -= (interporder[avgscheme[dir]]-1)/2;
2140  *pe += (interporder[avgscheme[dir]]-1)/2;
2141  }
2142 
2143 
2144  *is=INFULL1;
2145  *ie=OUTFULL1;
2146 
2147  *js=fluxloop[dir][FJS];
2148  *je=fluxloop[dir][FJE];
2149 
2150  }
2151  else if( interporflux==ENOAVG2CENTTYPE ){
2152  // input and output at same location
2153 
2154  // then edge quantity and need from 0 .. N
2155  *bs=myUconsloop[FKS]; //atch correct SASMARK; ALSO need to set up the same thing for other dimensions-- not sure if this is enough
2156  *be=myUconsloop[FKE];
2157 
2158  *ps=0;
2159  *pe=N3 - 1;
2160 
2161  *is=INFULL1;
2162  *ie=OUTFULL1;
2163 
2164  *js=fluxloop[dir][FJS];
2165  *je=fluxloop[dir][FJE];
2166  }
2167  else if(interporflux==ENOCENT2AVGTYPE){
2168  // then edge quantity and need from 0 .. N
2169  *bs=INFULL3;
2170  *be=OUTFULL3;
2171 
2172  *ps=myUconsloop[FKS];
2173  *pe=myUconsloop[FKE];
2174 
2175  *is=INFULL1;
2176  *ie=OUTFULL1;
2177 
2178  *js=fluxloop[dir][FJS];
2179  *je=fluxloop[dir][FJE];
2180 
2181  }
2182  else if(interporflux==ENOFLUXAVG3TYPE){
2183 
2184  // then edge quantity and need from 0 .. N
2185  *bs=INFULL3;
2186  *be=OUTFULL3;
2187 
2188  *ps=myUconsloop[FKS];
2189  *pe=myUconsloop[FKE];
2190 
2191  if(dir==1){
2192  // This method doesn't need "boundary" direction of "fluxes" since they don't actually exist
2193  *is=myUconsloop[FIS];
2194  *ie=myUconsloop[FIE]+SHIFT1;
2195  }
2196  else{
2197  // GODMARK: needed?
2198  *is=myUconsloop[FIS];
2199  *ie=myUconsloop[FIE];
2200  // *is=INFULL1;
2201  // *ie=OUTFULL1;
2202  }
2203 
2204  if(dir==2){
2205  // This method doesn't need "boundary" direction of "fluxes" since they don't actually exist
2206  *js=myUconsloop[FJS];
2207  *je=myUconsloop[FJE]+SHIFT2;
2208  }
2209  else{
2210  *js=myUconsloop[FJS];
2211  *je=myUconsloop[FJE];
2212  // GODMARK: needed?
2213  // *js=fluxloop[dir][FJS];
2214  // *je=fluxloop[dir][FJE];
2215  }
2216 
2217  }
2218 
2219  *di=1;
2220  *dj=1;
2221  *dk=1;
2222 
2223  if(withshifts){
2224  // shift due to grid sectionion SECTIONMARK
2225  *is += SHIFTX1DN;
2226  *ie += SHIFTX1UP;
2227  *js += SHIFTX2DN;
2228  *je += SHIFTX2UP;
2229  *bs += SHIFTX3DN;
2230  *be += SHIFTX3UP;
2231  *ps += SHIFTX3DN;
2232  *pe += SHIFTX3UP;
2233  }
2234 
2235  }
2236 
2237 
2238  // dualfprintf(fail_file,"pspe: %d %d %d %d %d %d %d %d %d %d %d %d %d %d %d\n",dir, interporflux, *is, *ie, *js, *je, *ks, *ke, *di, *dj, *dk, *bs, *ps, *pe, *be);
2239 
2240 
2241 }
2242 
2243 
2244 
2245 
2246 
2248 static void get_1d_line_c2e_multipl(int whichquantity, int dir, int interporflux, int bs, int ps, int pe, int be, int i, int j, int k, FTYPE (*p2interp)[NSTORE2][NSTORE3][NPR2INTERP], FTYPE (*yin)[2][NBIGM], struct of_trueijkp *trueijkp)
2249 {
2250  // for NUMPRIMLOOP:
2251  int nprlocalstart,nprlocalend;
2252  int nprlocallist[MAXNPR];
2253  int pllocal;
2254  int numprims;
2255  int plstart;
2256  int pl,pliter;
2257  int mypl;
2258  // others:
2259  int yiniter;
2260  int di2,dj2,dk2;
2261  int i2,j2,k2;
2262  int is2,ie2,js2,je2,ks2,ke2;
2263  int dir_exception[NDIM];
2264 
2265 
2266 
2268  //
2269  // Define which quantities (pl) to operate on
2270  //
2272 
2273  setup_nprlocalist(whichquantity,&nprlocalstart,&nprlocalend,nprlocallist,&numprims);
2274 
2275 
2276 
2277  trueijkp->dir=dir;
2278  trueijkp->iter=dir;
2279  trueijkp->interporflux=interporflux;
2280 
2281 
2282  // determine range of outer loop and range to feed to eno scheme
2283  if(dir==1){
2284 
2285  trueijkp->i=0;
2286  trueijkp->p=CENT;
2287  is2=bs;
2288  ie2=be;
2289 
2290  trueijkp->j=js2=je2=j;
2291 
2292  trueijkp->k=ks2=ke2=k;
2293 
2294  di2=1;
2295  dj2=1;
2296  dk2=1;
2297 
2298  }
2299  else if(dir==2){
2300 
2301  trueijkp->i=is2=ie2=i;
2302 
2303  trueijkp->j=0;
2304  trueijkp->p=CENT;
2305  js2=bs;
2306  je2=be;
2307 
2308  trueijkp->k=ks2=ke2=k;
2309 
2310  di2=1;
2311  dj2=1;
2312  dk2=1;
2313  }
2314  else if(dir==3){
2315 
2316  trueijkp->i=is2=ie2=i;
2317 
2318  trueijkp->j=js2=je2=j;
2319 
2320  trueijkp->k=0;
2321  trueijkp->p=CENT;
2322  ks2=bs;
2323  ke2=be;
2324 
2325  di2=1;
2326  dj2=1;
2327  dk2=1;
2328  }
2329  else{
2330  dualfprintf(fail_file,"No such dir=%d in get_1d_line_c2e_multipl()\n",dir);
2331  myexit(246344576);
2332  }
2333 
2334 
2335  if(DOENOFLUX != NOENOFLUX){
2336  dualfprintf(fail_file,"Reinvestigate why ENO needs this\n");
2337  myexit(9586);
2338  // JCM: expensive and not sure why doing it
2339  // reset to 0 so eno schemes don't care about values there (assume weights set to also 0 there)
2340  pl=0; // fake
2341  for(yiniter=-NBIGBND;yiniter<NBIG+NBIGBND;yiniter++){
2342  yin[pl][0][yiniter] = 0;
2343  }
2344  }
2345 
2346 
2347 
2348  // 1 input, assumed interporflux==ENOINTERPTYPE
2349  yiniter=bs;
2350  SUPERGENLOOP(i2,j2,k2,is2,ie2,js2,je2,ks2,ke2,di2,dj2,dk2){
2351  NUMPRIMLOOP(pliter,pl){ // as with assign_eno_result_c2e_multipl(), more cache friendly to have pl loop inside sinde p2interp has otherwise larger stride to reach other pl's compared to yin
2352  yin[pl][0][yiniter]=MACP0A1(p2interp,i2,j2,k2,pl);
2353  }
2354  yiniter++;
2355  }
2356 
2357 
2358  if(! (interporflux==ENOINTERPTYPE || interporflux==ENOINTERPTYPE4EMF) ){
2359  dualfprintf(fail_file,"get_1d_line_c2e_multipl only handles interporflux==ENOINTERPTYPE or ENOINTERPTYPE4EMF\n");
2360  myexit(25);
2361  }
2362 
2363 
2364 
2365 }
2366 
2367 
2369 static void get_1d_line_c2e(int dir, int interporflux, int pl, int bs, int ps, int pe, int be, int i, int j, int k, FTYPE (*p2interp)[NSTORE2][NSTORE3][NPR2INTERP], FTYPE *yin, struct of_trueijkp *trueijkp)
2370 {
2371  int yiniter;
2372  int di2,dj2,dk2;
2373  int i2,j2,k2;
2374  int is2,ie2,js2,je2,ks2,ke2;
2375  int dir_exception[NDIM];
2376 
2377 
2378  trueijkp->dir=dir;
2379  trueijkp->iter=dir;
2380  trueijkp->interporflux=interporflux;
2381 
2382 
2383  // determine range of outer loop and range to feed to eno scheme
2384  if(dir==1){
2385 
2386  trueijkp->i=0;
2387  trueijkp->p=CENT;
2388  is2=bs;
2389  ie2=be;
2390 
2391  trueijkp->j=js2=je2=j;
2392 
2393  trueijkp->k=ks2=ke2=k;
2394 
2395  di2=1;
2396  dj2=1;
2397  dk2=1;
2398 
2399  }
2400  else if(dir==2){
2401 
2402  trueijkp->i=is2=ie2=i;
2403 
2404  trueijkp->j=0;
2405  trueijkp->p=CENT;
2406  js2=bs;
2407  je2=be;
2408 
2409  trueijkp->k=ks2=ke2=k;
2410 
2411  di2=1;
2412  dj2=1;
2413  dk2=1;
2414  }
2415  else if(dir==3){
2416 
2417  trueijkp->i=is2=ie2=i;
2418 
2419  trueijkp->j=js2=je2=j;
2420 
2421  trueijkp->k=0;
2422  trueijkp->p=CENT;
2423  ks2=bs;
2424  ke2=be;
2425 
2426  di2=1;
2427  dj2=1;
2428  dk2=1;
2429  }
2430  else{
2431  dualfprintf(fail_file,"No such dir=%d in get_1d_line_c2e()\n",dir);
2432  myexit(246344576);
2433  }
2434 
2435 
2436  if(DOENOFLUX != NOENOFLUX){
2437  dualfprintf(fail_file,"Reinvestigate why ENO needs this\n");
2438  myexit(9587);
2439  // JCM: expensive and not sure why doing it
2440  // reset to 0 so eno schemes don't care about values there (assume weights set to also 0 there)
2441  for(yiniter=-NBIGBND;yiniter<NBIG+NBIGBND;yiniter++){
2442  yin[yiniter] = 0;
2443  }
2444  }
2445 
2446 
2447  // 1 input, assumed interporflux==ENOINTERPTYPE
2448  yiniter=bs;
2449  SUPERGENLOOP(i2,j2,k2,is2,ie2,js2,je2,ks2,ke2,di2,dj2,dk2){
2450 
2451  yin[yiniter]=MACP0A1(p2interp,i2,j2,k2,pl);
2452  yiniter++;
2453  }
2454 
2455  if(! (interporflux==ENOINTERPTYPE || interporflux==ENOINTERPTYPE4EMF) ){
2456  dualfprintf(fail_file,"get_1d_line_c2e only handles interporflux==ENOINTERPTYPE or ENOINTERPTYPE4EMF\n");
2457  myexit(25);
2458  }
2459 
2460 
2461 
2462 }
2463 
2465 void get_1d_line_c2e_gammaeffhydro(int dir, int interporflux, int pl, int bs, int ps, int pe, int be, int i, int j, int k, int idel, int jdel, int kdel, FTYPE (*p2interp)[NSTORE2][NSTORE3][NPR2INTERP], FTYPE *yin, struct of_trueijkp *trueijkp)
2466 {
2467  int num,locali,localj,localk,localloc;
2468  int yiniter;
2469  int di2,dj2,dk2;
2470  int i2,j2,k2;
2471  int is2,ie2,js2,je2,ks2,ke2;
2472  int dir_exception[NDIM];
2473  FTYPE rho, u, p, gamma;
2474 
2475 
2476  //SASMARK: not sure if this is req'd
2477  //trueijkp->dir=dir;
2478  //trueijkp->iter=dir;
2479  //trueijkp->interporflux=interporflux;
2480 
2481 
2482  // determine range of outer loop and range to feed to eno scheme
2483  if(dir==1){
2484 
2485  trueijkp->i=0;
2486  trueijkp->p=CENT;
2487  is2=bs;
2488  ie2=be;
2489 
2490  trueijkp->j=js2=je2=j;
2491 
2492  trueijkp->k=ks2=ke2=k;
2493 
2494  di2=1;
2495  dj2=1;
2496  dk2=1;
2497 
2498  }
2499  else if(dir==2){
2500 
2501  trueijkp->i=is2=ie2=i;
2502 
2503  trueijkp->j=0;
2504  trueijkp->p=CENT;
2505  js2=bs;
2506  je2=be;
2507 
2508  trueijkp->k=ks2=ke2=k;
2509 
2510  di2=1;
2511  dj2=1;
2512  dk2=1;
2513  }
2514  else if(dir==3){
2515 
2516  trueijkp->i=is2=ie2=i;
2517 
2518  trueijkp->j=js2=je2=j;
2519 
2520  trueijkp->k=0;
2521  trueijkp->p=CENT;
2522  ks2=bs;
2523  ke2=be;
2524 
2525  di2=1;
2526  dj2=1;
2527  dk2=1;
2528  }
2529 
2530  if(DOENOFLUX != NOENOFLUX){
2531  dualfprintf(fail_file,"Reinvestigate why ENO needs this\n");
2532  myexit(9588);
2533  // JCM: expensive and not sure why doing it
2534  // reset to 0 so eno schemes don't care about values there (assume weights set to also 0 there)
2535  for(yiniter=-NBIGBND;yiniter<NBIG+NBIGBND;yiniter++){
2536  yin[yiniter] = 0;
2537  }
2538  }
2539 
2540 
2541  // 1 input, assumed interporflux==ENOINTERPTYPE
2542  yiniter=bs;
2543  SUPERGENLOOP(i2,j2,k2,is2,ie2,js2,je2,ks2,ke2,di2,dj2,dk2){
2544 
2545  num=i2*idel+j2*jdel+k2*kdel;
2546 
2547  locali=trueijkp->i+di2*num;
2548  localj=trueijkp->j+dj2*num;
2549  localk=trueijkp->k+dk2*num;
2550  localloc=trueijkp->p;
2551 
2552  gamma = MACP0A1(p2interp,i2,j2,k2,VSQ); //Lorentz factor
2553  rho = MACP0A1(p2interp,i2,j2,k2,RHO); //comoving mass density
2554  u = MACP0A1(p2interp,i2,j2,k2,UU); //comoving internal energy
2555  p = pressure_rho0_u_simple(locali,localj,localk,localloc, rho, u ); //comoving plasma pressure
2556 
2557  //an indicator that is giving an effective $\gamma^2$ for a flow
2558  //in a shock where the bulk matter comes to rest, $p = \gamma^2 \rho$, hence $p/\rho$ is an estimate of the square of the $\gamma$-factor
2559  yin[yiniter]= gamma * ( fabs(rho) + fabs(u) + fabs(p) ) / ( fabs(rho) + SQRTMINNUMREPRESENT );
2560  yiniter++;
2561  }
2562 
2563  if(! (interporflux==ENOINTERPTYPE || interporflux==ENOINTERPTYPE4EMF) ){
2564  dualfprintf(fail_file,"get_1d_line_c2e only handles interporflux==ENOINTERPTYPE or ENOINTERPTYPE4EMF\n");
2565  myexit(25);
2566  }
2567 
2568 
2569 
2570 }
2571 
2572 
2574 static void get_1d_line(int dir, int interporflux, int pl, int bs, int ps, int pe, int be, int i, int j, int k, FTYPE (*p2interpm)[NSTORE2][NSTORE3][NPR],FTYPE (*p2interpp)[NSTORE2][NSTORE3][NPR], FTYPE (*yin)[NBIGM], struct of_trueijkp *trueijkp)
2575 {
2576  int yiniter;
2577  int di2,dj2,dk2;
2578  int i2,j2,k2;
2579  int is2,ie2,js2,je2,ks2,ke2;
2580  int dir_exception[NDIM];
2581 
2582 
2583  if(
2584  ((interporflux==ENOFLUXAVG1TYPE)&&(dir==1))||
2585  ((interporflux==ENOFLUXAVG2TYPE)&&(dir==2))||
2586  ((interporflux==ENOFLUXAVG3TYPE)&&(dir==3))
2587  ){
2588  dualfprintf(fail_file,"No such method with interporflux=%d and dir=%d\n",interporflux,dir);
2589  myexit(1);
2590  }
2591 
2592  dir_exception[1] = (interporflux==ENOFLUXAVG2TYPE) || (interporflux==ENOFLUXAVG3TYPE);
2593  dir_exception[2] = (interporflux==ENOFLUXAVG1TYPE) || (interporflux==ENOFLUXAVG3TYPE);
2594  dir_exception[3] = (interporflux==ENOFLUXAVG1TYPE) || (interporflux==ENOFLUXAVG2TYPE);
2595 
2596  trueijkp->dir=dir;
2597  trueijkp->interporflux=interporflux;
2598 
2599  // determine range of outer loop and range to feed to eno scheme
2600  if( ( (!dir_exception[1]) && (dir==1) ) || (interporflux==ENOFLUXAVG1TYPE) ){
2601 
2602  trueijkp->iter=1;
2603 
2604  trueijkp->i=0;
2605  trueijkp->p=CENT; // GODMARK -- probably off a bit
2606  is2=bs;
2607  ie2=be;
2608 
2609  trueijkp->j=js2=je2=j;
2610 
2611  trueijkp->k=ks2=ke2=k;
2612 
2613  di2=1;
2614  dj2=1;
2615  dk2=1;
2616 
2617  }
2618  else if( ( (!dir_exception[2]) && (dir==2) ) || (interporflux==ENOFLUXAVG2TYPE) ){
2619 
2620  trueijkp->iter=2;
2621 
2622  trueijkp->i=is2=ie2=i;
2623 
2624  trueijkp->j=0;
2625  trueijkp->p=CENT;
2626  js2=bs;
2627  je2=be;
2628 
2629  trueijkp->k=ks2=ke2=k;
2630 
2631  di2=1;
2632  dj2=1;
2633  dk2=1;
2634  }
2635  else if( ( (!dir_exception[3]) && (dir==3) ) || (interporflux==ENOFLUXAVG3TYPE) ){
2636 
2637  trueijkp->iter=3;
2638 
2639  trueijkp->i=is2=ie2=i;
2640 
2641  trueijkp->j=js2=je2=j;
2642 
2643  trueijkp->k=0;
2644  trueijkp->p=CENT;
2645  ks2=bs;
2646  ke2=be;
2647 
2648  di2=1;
2649  dj2=1;
2650  dk2=1;
2651  }
2652 
2653 
2654  if(DOENOFLUX != NOENOFLUX){
2655  dualfprintf(fail_file,"Reinvestigate why ENO needs this\n");
2656  myexit(9589);
2657  // JCM: expensive and not sure why doing it
2658  // reset to 0 so eno schemes don't care about values there (assume weights set to also 0 there)
2659  for(yiniter=-NBIGBND;yiniter<NBIG+NBIGBND;yiniter++){
2660  yin[0][yiniter] = yin[1][yiniter] = 0;
2661  }
2662  }
2663 
2664 
2665 
2666  if( (interporflux==ENOINTERPTYPE) || (interporflux==ENOINTERPTYPE4EMF)||(interporflux==ENOFLUXRECONTYPE)||(interporflux==ENOFLUXRECONTYPEGHOSTACTIVE)||(interporflux==ENOQUASIFIELDFLUXRECONTYPE)||(interporflux==ENOFLUXAVG1TYPE)||(interporflux==ENOFLUXAVG2TYPE)||(interporflux==ENOFLUXAVG3TYPE)||(interporflux==ENOAVG2CENTTYPE)||(interporflux==ENOCENT2AVGTYPE) ){// these have only 1 input
2667  yiniter=bs;
2668  SUPERGENLOOP(i2,j2,k2,is2,ie2,js2,je2,ks2,ke2,di2,dj2,dk2){
2669 
2670  yin[0][yiniter]=MACP0A1(p2interpm,i2,j2,k2,pl);
2671  yiniter++;
2672  }
2673  }
2674  else if(interporflux==ENOFLUXSPLITTYPE){ // this method has 2 inputs
2675  yiniter=bs;
2676  SUPERGENLOOP(i2,j2,k2,is2,ie2,js2,je2,ks2,ke2,di2,dj2,dk2){
2677 
2678  yin[0][yiniter]=MACP0A1(p2interpm,i2,j2,k2,pl);
2679  if(p2interpp!=NULL) yin[1][yiniter]=MACP0A1(p2interpp,i2,j2,k2,pl);
2680  yiniter++;
2681  }
2682  }
2683 
2684 
2685 
2686 
2687 
2688 }
2689 
2690 
2691 
2692 
2693 
2694 
2695 
2696 
2698 static void get_1d_line_shockarray(int dir, int interporflux, int bs, int ps, int pe, int be, int i, int j, int k, FTYPE (*shockarray)[NSTORE1][NSTORE2][NSTORE3], FTYPE (*shockindicator)[NBIGM], struct of_trueijkp *trueijkp)
2699 {
2700  int yiniter;
2701  int di2,dj2,dk2;
2702  int i2,j2,k2;
2703  int is2,ie2,js2,je2,ks2,ke2;
2704  int dir_exception[NDIM];
2705 
2706 
2707  if(
2708  ((interporflux==ENOFLUXAVG1TYPE)&&(dir==1))||
2709  ((interporflux==ENOFLUXAVG2TYPE)&&(dir==2))||
2710  ((interporflux==ENOFLUXAVG3TYPE)&&(dir==3))
2711  ){
2712  dualfprintf(fail_file,"No such method with interporflux=%d and dir=%d\n",interporflux,dir);
2713  myexit(1);
2714  }
2715 
2716  dir_exception[1] = (interporflux==ENOFLUXAVG2TYPE) || (interporflux==ENOFLUXAVG3TYPE);
2717  dir_exception[2] = (interporflux==ENOFLUXAVG1TYPE) || (interporflux==ENOFLUXAVG3TYPE);
2718  dir_exception[3] = (interporflux==ENOFLUXAVG1TYPE) || (interporflux==ENOFLUXAVG2TYPE);
2719 
2720  trueijkp->dir=dir;
2721  trueijkp->interporflux=interporflux;
2722 
2723  // determine range of outer loop and range to feed to eno scheme
2724  if( ( (!dir_exception[1]) && (dir==1) ) || (interporflux==ENOFLUXAVG1TYPE) ){
2725 
2726  trueijkp->iter=1;
2727 
2728  trueijkp->i=0;
2729  trueijkp->p=CENT; // GODMARK -- probably off a bit
2730  is2=bs;
2731  ie2=be;
2732 
2733  trueijkp->j=js2=je2=j;
2734 
2735  trueijkp->k=ks2=ke2=k;
2736 
2737  di2=1;
2738  dj2=1;
2739  dk2=1;
2740 
2741  }
2742  else if( ( (!dir_exception[2]) && (dir==2) ) || (interporflux==ENOFLUXAVG2TYPE) ){
2743 
2744  trueijkp->iter=2;
2745 
2746  trueijkp->i=is2=ie2=i;
2747 
2748  trueijkp->j=0;
2749  trueijkp->p=CENT;
2750  js2=bs;
2751  je2=be;
2752 
2753  trueijkp->k=ks2=ke2=k;
2754 
2755  di2=1;
2756  dj2=1;
2757  dk2=1;
2758  }
2759  else if( ( (!dir_exception[3]) && (dir==3) ) || (interporflux==ENOFLUXAVG3TYPE) ){
2760 
2761  trueijkp->iter=3;
2762 
2763  trueijkp->i=is2=ie2=i;
2764 
2765  trueijkp->j=js2=je2=j;
2766 
2767  trueijkp->k=0;
2768  trueijkp->p=CENT;
2769  ks2=bs;
2770  ke2=be;
2771 
2772  di2=1;
2773  dj2=1;
2774  dk2=1;
2775  }
2776 
2777 
2778  yiniter=bs;
2779  SUPERGENLOOP(i2,j2,k2,is2,ie2,js2,je2,ks2,ke2,di2,dj2,dk2){
2781  // ULTRASUPERGODMARK: Really must average input array to correct location. Use of shockindicator in para currently assumes indicator at the cell face in dir-direction after using the result of ficalc() from surrounding center cells
2783  shockindicator[EOMSETMHD][yiniter]=MACP1A0(shockarray,SHOCKPLDIR1+dir-1,i2,j2,k2);
2784 #if(RADSHOCKFLAT&&EOMRADTYPE!=EOMRADNONE)//KORAL
2785  shockindicator[EOMSETRAD][yiniter]=MACP1A0(shockarray,SHOCKRADPLDIR1+dir-1,i2,j2,k2);
2786 #endif
2787  yiniter++;
2788  }
2789 
2790 
2791 
2792 
2793 }
2794 
2795 
2796 
2797 
2798 
2801 static void causal_shift_order(int whichprimtype, int interporflux, int dir, int preforder, int bs, int ps, int pe, int be, int i, int j, int k, int idel, int jdel, int kdel, int *shift, int *minorder, int *maxorder)
2802 {
2803  int i3,j3,k3;
2804  int is3,ie3,js3,je3,ks3,ke3;
2805  int di3, dj3, dk3;
2806  int temporder;
2807  int superdiv;
2808  FTYPE wspeed0l,wspeed0r,wspeed1l,wspeed1r;
2809  FTYPE wspeed0ll,wspeed1ll;
2810  int yiniter;
2811  FTYPE localspeed[2];
2812  int shifttemp;
2813 
2814  int dir_exception[NDIM];
2815 
2816  if(
2817  ((interporflux==ENOFLUXAVG1TYPE)&&(dir==1))||
2818  ((interporflux==ENOFLUXAVG2TYPE)&&(dir==2))||
2819  ((interporflux==ENOFLUXAVG3TYPE)&&(dir==3))
2820  ){
2821  dualfprintf(fail_file,"No such method with interporflux=%d and dir=%d\n",interporflux,dir);
2822  myexit(1);
2823  }
2824 
2825  dir_exception[1] = (interporflux==ENOFLUXAVG2TYPE) || (interporflux==ENOFLUXAVG3TYPE);
2826  dir_exception[2] = (interporflux==ENOFLUXAVG1TYPE) || (interporflux==ENOFLUXAVG3TYPE);
2827  dir_exception[3] = (interporflux==ENOFLUXAVG1TYPE) || (interporflux==ENOFLUXAVG2TYPE);
2828 
2829  if( ( (!dir_exception[1]) && (dir==1) ) || (interporflux==ENOFLUXAVG1TYPE) ){
2830  is3=ps;
2831  ie3=pe;
2832 
2833  js3=je3=j;
2834 
2835  ks3=ke3=k;
2836 
2837  di3=1;
2838  dj3=1;
2839  dk3=1;
2840  }
2841  else if( ( (!dir_exception[2]) && (dir==2) ) || (interporflux==ENOFLUXAVG2TYPE) ){
2842  is3=ie3=i;
2843 
2844  js3=ps;
2845  je3=pe;
2846 
2847  ks3=ke3=k;
2848 
2849  di3=1;
2850  dj3=1;
2851  dk3=1;
2852  }
2853  else if( ( (!dir_exception[3]) && (dir==3) ) || (interporflux==ENOFLUXAVG3TYPE) ){
2854  is3=ie3=i;
2855 
2856  js3=je3=j;
2857 
2858  ks3=ps;
2859  ke3=pe;
2860 
2861  di3=1;
2862  dj3=1;
2863  dk3=1;
2864  }
2865 
2866 
2867  //#if( (STOREWAVESPEEDS)&& ((VCHARTYPE==GLOBALVCHAR)||(VCHARTYPE==LOCALVCHAR)) ) // this procedure makes no sense with GLOBALVCHAR
2868 #if( (STOREWAVESPEEDS==1)&& ((VCHARTYPE==LOCALVCHAR)||(VCHARTYPE==VERYLOCALVCHAR) ) )
2869  // GODMARK: This will use precomputed wave speeds in MACP2A0(wspeed,dir,2,i,j,k)
2870  // wspeed located at cell interface. Take this into account when forming shifter.
2871  // therefore wspeed is wave speeds from interface point of view. For centered quantities should consider average (or max/min) of wspeed.
2872  //
2873  // e.g. for c2e, average wspeed[i] and wspeed[i+1]
2874  // for a2c (as used for ENOFLUXRECONTYPE, not really center, but edge!) then wspeed[i] is correct one
2875  // for a2em/p average wspeed[i] and wspeed[i+1]
2876  //
2877  // shift -order/2-1 (e.g. -2 for WENO5) if flow is superRIGHT (when both left/right chars are +)
2878  // shift order/2-1 (e.g. 2 for WENO5) if flow is superLEFT (when both left/right chars are -)
2879  //
2880  // smoothly vary between and feed integer value (which is interpreted correctly for avg2cent vs. avg2edge vs. cent2edge
2881  yiniter=ps; // note starts at ps not bs since shift only used on points of interest
2882  SUPERGENLOOP(i3,j3,k3,is3,ie3,js3,je3,ks3,ke3,di3,dj3,dk3){ // only over points of interest
2883 
2884  // first determine if should reduce order because of supersonic divergence
2885  // first assume maximum preferred order
2886  temporder=preforder;
2887 
2888  // get standard wavespeed (used by any method)
2889  wspeed0l=GLOBALMACP3A0(wspeed,EOMSETMHD,dir,0,i3,j3,k3);
2890  wspeed1l=GLOBALMACP3A0(wspeed,EOMSETMHD,dir,1,i3,j3,k3);
2891 
2892 
2893  if((interporflux==ENOINTERPTYPE)||(interporflux==ENOINTERPTYPE4EMF)||(interporflux==ENOFLUXSPLITTYPE)||(interporflux==ENOFLUXAVG1TYPE)||(interporflux==ENOFLUXAVG2TYPE)||(interporflux==ENOFLUXAVG3TYPE)){ // quantities are at CENT-dir (assumes ENOFLUXAVG?TYPE is orthogonal to dir)
2894 
2895 
2896  // get grid-on-right wavespeed
2897  wspeed0r=GLOBALMACP3A0(wspeed,EOMSETMHD,dir,0,i3+idel,j3+jdel,k3+kdel);
2898  wspeed1r=GLOBALMACP3A0(wspeed,EOMSETMHD,dir,1,i3+idel,j3+jdel,k3+kdel);
2899 
2900 
2901  // check for superfast divergence
2902  superdiv=0;
2903  if(SUPERFASTDIVREDUCE&&(wspeed1l<0)&&(wspeed0r>0)){
2904  superdiv=1;
2905  temporder=1;
2906  }
2907 
2908  // get correctly positioned wave speed
2909  localspeed[0]=min(wspeed0l,wspeed0r);
2910  localspeed[1]=max(wspeed1l,wspeed1r);
2911 
2912  }
2913  else if((interporflux==ENOFLUXRECONTYPE)||(interporflux==ENOFLUXRECONTYPEGHOSTACTIVE)||(interporflux==ENOQUASIFIELDFLUXRECONTYPE)||(interporflux==ENOAVG2CENTTYPE)||(interporflux==ENOCENT2AVGTYPE)){ // quantities are at FACE-dir
2914 
2915  // get grid-on-right wavespeed
2916  wspeed0r=GLOBALMACP3A0(wspeed,EOMSETMHD,dir,0,i3+idel,j3+jdel,k3+kdel);
2917  wspeed1r=GLOBALMACP3A0(wspeed,EOMSETMHD,dir,1,i3+idel,j3+jdel,k3+kdel);
2918 
2919  // get grid-on-left wavespeed
2920  wspeed0ll=GLOBALMACP3A0(wspeed,EOMSETMHD,dir,0,i3-idel,j3-jdel,k3-kdel);
2921  wspeed1ll=GLOBALMACP3A0(wspeed,EOMSETMHD,dir,1,i3-idel,j3-jdel,k3-kdel);
2922 
2923 
2924  // check for superfast divergence
2925  superdiv=0;
2926  if(SUPERFASTDIVREDUCE&&(wspeed1ll<0)&&(wspeed0r>0)){
2927  superdiv=1;
2928  temporder=1;
2929  }
2930 
2931 
2932  // get correctly positioned wave speed as average of surroundings
2933  // wave speeds at interface
2934  localspeed[0]=wspeed0l;
2935  localspeed[1]=wspeed1l;
2936  }
2937 
2938  // only shift if not diverging superfast
2939  if(superdiv==0){
2940  // linear interpolation between left and right super"sonic" directed shifts back downstream
2941  shifttemp=-((temporder+1)/2)-(temporder+1)*(FTYPE)(localspeed[0]/(localspeed[1]-localspeed[0]+SMALL));
2942  if(shifttemp>(temporder+1)/2) shifttemp=(temporder+1)/2;
2943  if(shifttemp<-((temporder+1)/2)) shifttemp=-((temporder+1)/2);
2944  }
2945  else shifttemp=0; // no shift if super divergence
2946 
2947  // now assign results to arrays of minorder,maxorder, and shifts.
2948  minorder[yiniter]=MIN(MINPREFORDER,temporder);
2949  maxorder[yiniter]=temporder;
2950  shift[yiniter]=shifttemp;
2951  yiniter++;
2952  }
2953 #else
2954  // NO SHIFT OR CHANGE OF ORDER
2955  yiniter=ps; // GODMARK: was bs
2956  SUPERGENLOOP(i3,j3,k3,is3,ie3,js3,je3,ks3,ke3,di3,dj3,dk3){ // only over points of interest //atch correct -- was di3, dj3, dk3
2957  minorder[yiniter]=MIN(preforder,MINPREFORDER); // minimum preferred order
2958  maxorder[yiniter]=preforder;
2959  shift[yiniter]=0;// then no shift since don't have wave speeds (not stored)
2960  yiniter++;
2961  }
2962 #endif
2963 
2964 }
2965 
2966 
2967 
2968 
2970 static int get_V_and_P(int whichprimtype, int interporflux, int dir, int bs, int ps, int pe, int be, int i, int j, int k, int idel, int jdel, int kdel, FTYPE (*yrealin)[2][NBIGM], FTYPE (*Vline)[NBIGM], FTYPE (*Pline)[NBIGM], struct of_trueijkp *trueijkp)
2971 {
2972  int num,locali,localj,localk,localloc;
2973  int di,dj,dk;
2974  int i3,j3,k3;
2975  int is3,ie3,js3,je3,ks3,ke3;
2976  int di3, dj3, dk3;
2977  FTYPE interplistpl[MAXNPR][MAXSPACEORDER];
2978  FTYPE *yrealpl[MAXNPR]; // number of pointers
2979  int plpl;
2980  int startorderi,endorderi;
2981  int yiniter;
2982  int nprlocalstart,nprlocalend;
2983  int nprlocallist[MAXNPR];
2984  int pllocal;
2985  int numprims;
2986  int iii;
2987  FTYPE myprim[MAXNPR],bsq;
2988  struct of_state qdontuse;
2989  struct of_state *qptr=&qdontuse;
2990  int pl,pliter;
2991  struct of_geom geomdontuse;
2992  struct of_geom *ptrgeom=&geomdontuse;
2993 
2994 
2995 
2997  //
2998  // Define which quantities (pl) to operate on
2999  //
3001 
3002  setup_nprlocalist(whichprimtype,&nprlocalstart,&nprlocalend,nprlocallist,&numprims);
3003 
3004 
3005 
3006 
3007 
3008 
3009 
3010 
3011  if(((interporflux!=ENOFLUXAVG1TYPE)&&(dir==1))||(interporflux==ENOFLUXAVG1TYPE)){
3012  is3=bs;
3013  ie3=be;
3014 
3015  js3=je3=j;
3016 
3017  ks3=ke3=k;
3018 
3019  di3=1;
3020  dj3=1;
3021  dk3=1;
3022  }
3023  else if(((interporflux!=ENOFLUXAVG2TYPE)&&(dir==2))||(interporflux==ENOFLUXAVG2TYPE)){
3024  is3=ie3=i;
3025 
3026  js3=bs;
3027  je3=be;
3028 
3029  ks3=ke3=k;
3030 
3031  di3=1;
3032  dj3=1;
3033  dk3=1;
3034  }
3035  else if(((interporflux!=ENOFLUXAVG3TYPE)&&(dir==3))||(interporflux==ENOFLUXAVG3TYPE)){
3036  is3=ie3=i;
3037 
3038  js3=je3=j;
3039 
3040  ks3=bs;
3041  ke3=be;
3042 
3043  di3=1;
3044  dj3=1;
3045  dk3=1;
3046  }
3047  else{
3048  dualfprintf(fail_file,"No such dir=%d in get_V_and_P()\n",dir);
3049  myexit(246344572);
3050  }
3051 
3052  // trueijkp->iter should really be passed
3053  di=(trueijkp->iter==1);
3054  dj=(trueijkp->iter==2);
3055  dk=(trueijkp->iter==3);
3056 
3057 
3058  yiniter=bs; // per-point, get maximum can get
3059  SUPERGENLOOP(i3,j3,k3,is3,ie3,js3,je3,ks3,ke3,di3,dj3,dk3){ // only over points of interest
3060 
3061  num=iii=i3*idel+j3*jdel+k3*kdel;
3062 
3063  // trueijkp->i,trueijkp->j,trueijkp->k,trueijkp->p should really be passed
3064  // needed for EOS in general
3065  locali=trueijkp->i+di*num;
3066  localj=trueijkp->j+dj*num;
3067  localk=trueijkp->k+dk*num;
3068  localloc=trueijkp->p;
3069 
3070 #if(PLINEWITHFIELD || VLINEWITHGDETRHO)
3071  PALLREALLOOP(pl) myprim[pl] = yrealin[pl][0][num];
3072  get_geometry(locali,localj,localk,localloc,ptrgeom);
3073  // need u^\mu and Bsq (currently from b^\mu b_\mu)
3074  // OPTMARK: get_state here is very expensive since only required to have per-centered-point and not per direction, so 3X over computed unless refer to stored data.
3075  // GODMARK: Using stored data presumes yrealin is at CENT, which is true for any use of get_V_and_P() right now
3076  // Actually do have data at left-right for each dir in fluxstate[dir]. So can check where yrealin is located using localglobal.
3077  get_stateforinterpline(myprim,ptrgeom,&qptr); // OPTMARK :Refer to centered state if stored
3078 #endif
3079 
3080 
3081 
3082 #if(VLINEWITHGDETRHO==0)
3083  Vline[EOMSETMHD][yiniter]=yrealin[UU+dir][0][iii];
3084 #else
3085  // \detg rho u^dir (approximately)
3086  // get_geometry_gdetonly(locali,localj,localk,localloc,ptrgeom);
3087  // Vline[yiniter]=(ptrgeom->g)*yrealin[RHO][0][iii]*yrealin[UU+dir][0][iii];
3088  // \detg rho u^dir
3089  Vline[EOMSETMHD][yiniter]=(ptrgeom->gdet)*yrealin[RHO][0][num]*(qptr->ucon[dir]);
3090 #endif
3091 
3092 #if(RADSHOCKFLAT&&EOMRADTYPE!=EOMRADNONE)
3093 #if(VLINEWITHGDETRHO==0)
3094  Vline[EOMSETRAD][yiniter]=yrealin[PRAD0+dir][0][iii];
3095 #else
3096  Vline[EOMSETRAD][yiniter]=(ptrgeom->gdet)*yrealin[PRAD0][0][num]*(qptr->uradcon[dir]); // approximate KORALTODO
3097 #endif
3098 #endif
3099 
3100 
3101 
3102  // OPTMARK: Since assuming at CENT with normal primitive, then already have correct pressure
3103  // P[yiniter]=pressure_rho0_u_simple(locali,localj,localk,localloc,yrealin[RHO][0][iii],yrealin[UU][0][iii]);
3104  Pline[EOMSETMHD][yiniter]=qptr->pressure;
3105 #if(PLINEWITHFIELD==0)
3106  // done
3107 #else
3108  // need total pressure in general
3109  // add magnetic pressure
3110  // OPTMARK: Since assuming at CENT with normal primitive, then already have correct b^2
3111  // bsq = dot(q.bcon, q.bcov); // now store bsq
3112  Pline[EOMSETMHD][yiniter] += 0.5*(qptr->bsq);
3113 #endif
3114 
3115 #if(RADSHOCKFLAT&&EOMRADTYPE!=EOMRADNONE) // KORAL
3116  // add radiation pressure to total pressure if optically thick
3117  FTYPE tautot[NDIM],tautotmax;
3118  //calcfull_tautot(myprim, ptrgeom, tautot, &tautotmax);
3119  calc_tautot(myprim, ptrgeom, NULL, tautot, &tautotmax); // very accurate tautot not necessary, so use Tgas=Trad assumption in opacity
3120  Pline[EOMSETRAD][yiniter]=(4.0/3.0-1.0)*yrealin[PRAD0][0][num]; // radiative pressure
3121  // add radiation pressure to total pressure if optically thick
3122  // KORALNOTE: recall pressure just along diagonal and no velocity in R^\mu_\nu
3123  Pline[EOMSETMHD][yiniter]+=MIN(tautotmax,1.0)*Pline[EOMSETRAD][yiniter];
3124 #endif
3125 
3126 
3127  yiniter++;
3128  }
3129 
3130  return( 0 );
3131 
3132 }
3133 
3134 
3135 
3136 
3138 static int get_shock_indicator(int whichprimtype, int interporflux, int dir, int bs, int ps, int pe, int be, int i, int j, int k, int idel, int jdel, int kdel, FTYPE (*yprim)[2][NBIGM], FTYPE (*Vline)[NBIGM], FTYPE (*Pline)[NBIGM], FTYPE (*shockindicator)[NBIGM], struct of_trueijkp *trueijkp)
3139 {
3140  int i3,j3,k3;
3141  int is3,ie3,js3,je3,ks3,ke3;
3142  int di3, dj3, dk3;
3143 #if(0)
3144  FTYPE yinterplistpl[MAXNPR][MAXSPACEORDER];
3145  FTYPE *ypl[MAXNPR]; // number of pointers
3146  int plpl;
3147  int startorderi,endorderi;
3148  int l;
3149 #endif
3150  int yiniter;
3151  extern FTYPE Ficalc(int dir, FTYPE *V, FTYPE *P);
3152  int nprlocalstart,nprlocalend;
3153  int nprlocallist[MAXNPR];
3154  int pllocal;
3155  int numprims;
3156  int num,locali,localj,localk,localloc;
3157  int di,dj,dk;
3158 
3159 
3160 
3161 
3163  //
3164  // Define which quantities (pl) to operate on
3165  //
3167 
3168  // NO! While Ftilde() doesn't yet use ypl, feed get_shock_indicator() yprim, so need full prim loop
3169  // setup_nprlocalist(whichprimtype,&nprlocalstart,&nprlocalend,nprlocallist,&numprims);
3170 
3171 
3172 
3173 
3174  // bs+2 .. be-2 since Ficalc needs +-2 points
3175 
3176 
3177 
3178  if(((interporflux!=ENOFLUXAVG1TYPE)&&(dir==1))||(interporflux==ENOFLUXAVG1TYPE)){
3179  is3=bs+2;
3180  ie3=be-2;
3181 
3182  js3=je3=j;
3183 
3184  ks3=ke3=k;
3185 
3186  di3=1;
3187  dj3=1;
3188  dk3=1;
3189  }
3190  else if(((interporflux!=ENOFLUXAVG2TYPE)&&(dir==2))||(interporflux==ENOFLUXAVG2TYPE)){
3191  is3=ie3=i;
3192 
3193  js3=bs+2;
3194  je3=be-2;
3195 
3196  ks3=ke3=k;
3197 
3198  di3=1;
3199  dj3=1;
3200  dk3=1;
3201  }
3202  else if(((interporflux!=ENOFLUXAVG3TYPE)&&(dir==3))||(interporflux==ENOFLUXAVG3TYPE)){
3203  is3=ie3=i;
3204 
3205  js3=je3=j;
3206 
3207  ks3=bs+2;
3208  ke3=be-2;
3209 
3210  di3=1;
3211  dj3=1;
3212  dk3=1;
3213  }
3214  else{
3215  dualfprintf(fail_file,"No such interporflux=%d in get_shock_indicator()\n",interporflux);
3216  myexit(9894386);
3217  }
3218 
3219 
3220 
3221  yiniter=bs+2; // starts at bs+2
3222 
3223  if(MAXBND<4){
3224  dualfprintf(fail_file,"MAXBND should be 4 for shockindicator???\n");
3225  myexit(8465684);
3226  }
3227 
3228 
3229 #if(0)
3230  startorderi = - (7)/2; // order=7 fixed for shock detector
3231  endorderi = - startorderi;
3232  // shift pointer
3233  PALLREALLOOP(plpl){
3234  ypl[plpl] = yinterplistpl[plpl] - startorderi;
3235  }
3236 #endif
3237 
3238  // trueijkp->iter should really be passed
3239  di=(trueijkp->iter==1);
3240  dj=(trueijkp->iter==2);
3241  dk=(trueijkp->iter==3);
3242 
3243 
3244  SUPERGENLOOP(i3,j3,k3,is3,ie3,js3,je3,ks3,ke3,di3,dj3,dk3){ // only over points of interest
3245 
3246 #if(0)
3247  // Ficalc doesn't really use ypl right now
3248  PALLREALLOOP(plpl){
3249  // get interpolation points, where y[0] is point of interest for which interpolation is found.
3250  for(l=startorderi;l<=endorderi;l++){
3251  ypl[plpl][l]=yprim[plpl][0][i3*idel+j3*jdel+k3*kdel + l];
3252  }
3253  }
3254 #endif
3255  // dualfprintf(fail_file,"idel=%d jdel=%d kdel=%d :: i3=%d j3=%d k3=%d\n",idel,jdel,kdel,i3,j3,k3);
3256 
3257  //shockindicator[EOMSETMHD][yiniter]=Ficalc(dir,&Vline[EOMSETMHD][yiniter],&Pline[EOMSETMHD][yiniter],ypl);
3258  shockindicator[EOMSETMHD][yiniter]=Ficalc(dir,&Vline[EOMSETMHD][yiniter],&Pline[EOMSETMHD][yiniter]);
3259 
3260 #if(RADSHOCKFLAT&&EOMRADTYPE!=EOMRADNONE)
3261  shockindicator[EOMSETRAD][yiniter]=Ficalc(dir,&Vline[EOMSETRAD][yiniter],&Pline[EOMSETRAD][yiniter]); // KORAL
3262 #endif
3263 
3264 
3265 #if(0 && FLUXDUMP) // DEBUG: (turn off FLUXDUMP in flux.c)
3266  // trueijkp->i,trueijkp->j,trueijkp->k,trueijkp->p should really be passed
3267  // ijkcurr needed for EOS in general
3268  num=i3*idel+j3*jdel+k3*kdel;
3269  locali=trueijkp->i+di*num;
3270  localj=trueijkp->j+dj*num;
3271  localk=trueijkp->k+dk*num;
3272  localloc=trueijkp->p;
3273 
3274  GLOBALMACP0A1(fluxdump,locali,localj,localk,4*NPR + (dir-1)*NPR*5 + NPR*0 + RHO)=shockindicator[EOMSETMHD][yiniter]; // don't look at EOMSETRAD
3275 
3276 #endif
3277 
3278 
3279  yiniter++;
3280  }
3281 
3282  return( 0 );
3283 
3284 }
3285 
3286 
3287 
3289 static int get_contact_indicator(int realisinterp, int whichprimtype, int interporflux, int dir, int bs, int ps, int pe, int be, int i, int j, int k, int idel, int jdel, int kdel, FTYPE (*yin)[2][NBIGM], FTYPE (*Vline)[NBIGM], FTYPE (*Pline)[NBIGM], FTYPE (*etai)[NUMTRUEEOMSETS][NBIGM])
3290 {
3291  int i3,j3,k3;
3292  int is3,ie3,js3,je3,ks3,ke3;
3293  int di3, dj3, dk3;
3294  FTYPE yinterplistpl[MAXNPR][MAXSPACEORDER];
3295  FTYPE *ypl[MAXNPR]; // number of pointers
3296  int plpl, pliter;
3297  int startorderi,endorderi;
3298  int yiniter;
3299  int l;
3300  extern FTYPE etaicalc(int pl, FTYPE *y, FTYPE *V, FTYPE *P);
3301  int nprlocalstart,nprlocalend;
3302  int nprlocallist[MAXNPR];
3303  int pllocal;
3304  int numprims;
3305 
3306 
3307 
3308 
3310  //
3311  // Define which quantities (pl) to operate on
3312  //
3314 
3315  // input is expected to be line of data, so correct to use NUMPRIMLOOP, but in the end only apply if realisinterp==1 and doing pl==RHO
3316  setup_nprlocalist(whichprimtype,&nprlocalstart,&nprlocalend,nprlocallist,&numprims);
3317 
3318 
3319 
3320 
3321  // bs+2 .. be-2 since etaicalc needs +-2 points
3322 
3323 
3324 
3325 
3326  if(((interporflux!=ENOFLUXAVG1TYPE)&&(dir==1))||(interporflux==ENOFLUXAVG1TYPE)){
3327  is3=bs+2;
3328  ie3=be-2;
3329 
3330  js3=je3=j;
3331 
3332  ks3=ke3=k;
3333 
3334  di3=1;
3335  dj3=1;
3336  dk3=1;
3337  }
3338  else if(((interporflux!=ENOFLUXAVG2TYPE)&&(dir==2))||(interporflux==ENOFLUXAVG2TYPE)){
3339  is3=ie3=i;
3340 
3341  js3=bs+2;
3342  je3=be-2;
3343 
3344  ks3=ke3=k;
3345 
3346  di3=1;
3347  dj3=1;
3348  dk3=1;
3349  }
3350  else if(((interporflux!=ENOFLUXAVG3TYPE)&&(dir==3))||(interporflux==ENOFLUXAVG3TYPE)){
3351  is3=ie3=i;
3352 
3353  js3=je3=j;
3354 
3355  ks3=bs+2;
3356  ke3=be-2;
3357 
3358  di3=1;
3359  dj3=1;
3360  dk3=1;
3361  }
3362  else{
3363  dualfprintf(fail_file,"No such interporflux=%d in get_contact_indicator()\n",interporflux);
3364  myexit(9894386);
3365  }
3366 
3367 
3368 
3369  yiniter=bs+2; // bs+2
3370 
3371  if(MAXBND<4){
3372  dualfprintf(fail_file,"MAXBND should be 4 for contactindicator???\n");
3373  myexit(8465684);
3374  }
3375 
3376  startorderi = - (7)/2; // order=7 fixed for shock detector (must make sure MAXBND==5)
3377  endorderi = - startorderi;
3378 
3379  // shift pointer
3380  NUMPRIMLOOP(pliter,plpl){
3381  ypl[plpl] = yinterplistpl[plpl] - startorderi;
3382  }
3383 
3384 
3385  SUPERGENLOOP(i3,j3,k3,is3,ie3,js3,je3,ks3,ke3,di3,dj3,dk3){ // only over points of interest
3386 
3387  NUMPRIMLOOP(pliter,plpl){
3388 
3389 
3390  if(plpl==RHO && realisinterp){
3391  // get interpolation points, where y[0] is point of interest for which interpolation is found.
3392  for(l=startorderi;l<=endorderi;l++){
3393  ypl[plpl][l]=yin[plpl][0][i3*idel+j3*jdel+k3*kdel + l];
3394  }
3395 
3396  etai[plpl][EOMSETMHD][yiniter]=etaicalc(plpl, ypl[plpl], &Vline[EOMSETMHD][yiniter], &Pline[EOMSETMHD][yiniter]);
3397 #if(RADSHOCKFLAT&&EOMRADTYPE!=EOMRADNONE)
3398  etai[plpl][EOMSETRAD][yiniter]=etaicalc(plpl, ypl[plpl], &Vline[EOMSETRAD][yiniter], &Pline[EOMSETRAD][yiniter]);
3399 #endif
3400  }
3401  else{
3402  etai[plpl][EOMSETMHD][yiniter]=0.0;
3403 #if(RADSHOCKFLAT&&EOMRADTYPE!=EOMRADNONE)
3404  etai[plpl][EOMSETRAD][yiniter]=0.0;
3405 #endif
3406  }
3407 
3408  }
3409 
3410  yiniter++;
3411  }
3412 
3413  return( 0 );
3414 
3415 
3416 
3417 }
3418 
3419 
3420 
3421 
3423 static int compute_df_line(int doingweno,int interporflux, int recontype, int whichreduce, int preforder, int pl, int bs, int ps, int pe, int be, int *minorder, int *maxorder, int *shift, FTYPE *yin, FTYPE (*df)[NBIGM])
3424 {
3425  int i;
3426  extern void get_limit_slopes_paraline(FTYPE *dq1l, FTYPE *dq1r, FTYPE *dq2, FTYPE *dq);
3427 
3428 
3429 
3431  //
3432  // check that NUMDFS is sufficient
3433  //
3435 
3436  // if doingweno==0 then assume things needed are like what's needed for paraline
3437  if(NUMDFS<4){
3438  dualfprintf(fail_file,"paraline requires NUMDFS=4 while set to %d\n",NUMDFS);
3439  myexit(8923658);
3440  }
3441  if(DOMONOINTERP!= NOMONOINTERP && NUMDFS<5){
3442  dualfprintf(fail_file,"MONO requires NUMDFS=5 while set to %d\n",NUMDFS);
3443  myexit(8923659);
3444  }
3445 
3446 
3447  // for SWENO: df[0] was one-sided :: df[1] was d^2f :: df[3] was centered diff (but no factor of 0.5!!! GODMARK) df[4] was df centered 2 apart
3448 
3449 #if(NUMDFS>DFONESIDED)
3450  // df
3451  // df (one-sided difference : used to get Dqp,Dqm)
3452  // truly centered on face if y on CENT, but effectively used for CENT
3453  for(i=bs+1;i<=be;i++) df[DFONESIDED][i] = yin[i]-yin[i-1];
3454 #endif
3455 
3456  // df (centered difference : used to get Dqc)
3457  // On CENT if y CENT
3458  // df centered, for Sasha's smoothness indicators
3459 #if(NUMDFS>DFCENT)
3460  for(i=bs+1;i<=be-1;i++) df[DFCENT][i] = 0.5*(yin[i+1]-yin[i-1]);
3461 #endif
3462 
3463 
3464 #if(NUMDFS>DFCENT2APART)
3465  // df centered 2 apart, for Sasha's smoothness indicators
3466  for(i=bs+2;i<=be-2;i++) df[DFCENT2APART][i] = 0.25*(yin[i+2]-yin[i-2]);
3467 #endif
3468 
3469 #if(NUMDFS>DF2OFONESIDED)
3470  // d^2f
3471  for(i=bs+1;i<=be-1;i++) df[DF2OFONESIDED][i] = df[DFONESIDED][i+1]-df[DFONESIDED][i];
3472 #endif
3473 
3474 
3475 #if(NUMDFS>DFMONO)
3476  // PARALIM 3-point MONOTONIZED slope
3477  // on face if y CENT and so df[DFCENT] CENT
3478  for(i=bs+1;i<=be-1;i++){
3479  get_limit_slopes_paraline(&df[DFONESIDED][i], &df[DFONESIDED][i+1], &df[DFCENT][i], &df[DFMONO][i]);
3480  }
3481 #endif
3482 
3483 #if(NUMDFS>DF2OFMONO)
3484  // d^2f (compute second derivative on *monotonized* slopes used for para)
3485  // at CENT if y at CENT
3486  for(i=bs+2;i<=be-1;i++) df[DF2OFMONO][i] = df[DFMONO][i]-df[DFMONO][i-1];
3487 #endif
3488 
3489  // assume all other things computed within PARA-based routines
3490 
3491 
3492 
3493  return(0);
3494 }
3495 
3496 
3499 static int compute_df_line_paraline(int doingweno,int interporflux, int recontype, int whichreduce, int preforder, int pl, int bs, int ps, int pe, int be, int *minorder, int *maxorder, int *shift, FTYPE *yin, FTYPE (*df)[NBIGM])
3500 {
3501  int i;
3502  extern void get_limit_slopes_paraline(FTYPE *dq1l, FTYPE *dq1r, FTYPE *dq2, FTYPE *dq);
3503 
3504 
3505 
3507  //
3508  // check that NUMDFS is sufficient
3509  //
3511 
3512  // if doingweno==0 then assume things needed are like what's needed for paraline
3513 #if(NUMDFS<4)
3514  dualfprintf(fail_file,"paraline requires NUMDFS=4 while set to %d\n",NUMDFS);
3515  myexit(8923658);
3516 #endif
3517 
3518 
3519 #if(NUMDFS>DFONESIDED)
3520  // df
3521  // df (one-sided difference : used to get Dqp,Dqm)
3522  // truly centered on face if y on CENT, but effectively used for CENT
3523  for(i=bs+1;i<=be;i++) df[DFONESIDED][i] = yin[i]-yin[i-1];
3524 #endif
3525 
3526  // df (centered difference : used to get Dqc)
3527  // On CENT if y CENT
3528  // df centered, for Sasha's smoothness indicators
3529 #if(NUMDFS>DFCENT)
3530  for(i=bs+1;i<=be-1;i++) df[DFCENT][i] = 0.5*(yin[i+1]-yin[i-1]);
3531 #endif
3532 
3533 
3534 
3535 #if(NUMDFS>DFMONO)
3536  // PARALIM 3-point MONOTONIZED slope
3537  // on face if y CENT and so df[DFCENT] CENT
3538  for(i=bs+1;i<=be-1;i++){
3539  get_limit_slopes_paraline(&df[DFONESIDED][i], &df[DFONESIDED][i+1], &df[DFCENT][i], &df[DFMONO][i]);
3540  }
3541 #endif
3542 
3543 #if(NUMDFS>DF2OFMONO)
3544  // d^2f (compute second derivative on *monotonized* slopes used for para)
3545  // at CENT if y at CENT
3546  for(i=bs+2;i<=be-1;i++) df[DF2OFMONO][i] = df[DFMONO][i]-df[DFMONO][i-1];
3547 #endif
3548 
3549  // assume all other things computed within PARA-based routines
3550 
3551 
3552 
3553  return(0);
3554 }
3555 
3556 
3558 static int compute_df_line_formono(int doingweno,int interporflux, int recontype, int whichreduce, int preforder, int pl, int bs, int ps, int pe, int be, int *minorder, int *maxorder, int *shift, FTYPE *yin, FTYPE (*df)[NBIGM])
3559 {
3560  int i;
3561 
3562 
3563 
3564 #if(NUMDFS>DFONESIDED)
3565  // df
3566  // df (one-sided difference : used to get Dqp,Dqm)
3567  // truly centered on face if y on CENT, but effectively used for CENT
3568  for(i=bs+1;i<=be;i++) df[DFONESIDED][i] = yin[i]-yin[i-1];
3569 #else
3570  if(DOMONOINTERP==SMONOINTERP){
3571  dualfprintf(fail_file,"Need NUMDFS=%d bigger for SMONO to be used1\n",NUMDFS);
3572  }
3573 #endif
3574 
3575 #if(NUMDFS>DF2OFONESIDED)
3576  // d^2f
3577  for(i=bs+1;i<=be-1;i++) df[DF2OFONESIDED][i] = df[DFONESIDED][i+1]-df[DFONESIDED][i];
3578 #else
3579  if(DOMONOINTERP==SMONOINTERP){
3580  dualfprintf(fail_file,"Need NUMDFS=%d bigger for SMONO to be used2\n",NUMDFS);
3581  }
3582 #endif
3583 
3584 
3585 
3586 
3587 
3588  return(0);
3589 }
3590 
3591 
3593 int compute_df_line_new(int doingweno, int whichprimtype, int interporflux, int recontype, int dir, int whichreduce, int preforder, int bs, int ps, int pe, int be, int *minorder, int *maxorder, int *shift, FTYPE (*yprim)[2][NBIGM], FTYPE *df, FTYPE *stiffindicator, FTYPE (*Vline)[NBIGM], FTYPE (*Pline)[NBIGM], struct of_trueijkp *trueijkp)
3594 {
3595  int i;
3596  int num,locali,localj,localk,localloc;
3597  int pl,pliter;
3598  int di,dj,dk;
3599  FTYPE myprim[MAXNPR];
3600  FTYPE mypriml[MAXNPR];
3601  FTYPE myprimr[MAXNPR];
3602  FTYPE myPgas,myPtot;
3603  // FTYPE U[MAXNPR];
3604  struct of_state qdontuse;
3605  struct of_state *qptr=&qdontuse;
3606  FTYPE bsq;
3607  FTYPE btt;
3608  FTYPE par,parl,parr;
3609  FTYPE uparl,uparr,upar;
3610  FTYPE veffl,veffr,veff;
3611  int nprlocalstart,nprlocalend;
3612  int nprlocallist[MAXNPR];
3613  int pllocal;
3614  int numprims;
3615  FTYPE stifffactor,stifffactor1,stifffactor2;
3616  FTYPE a_Pgasline[NBIGM];
3617  FTYPE (*Pgasline);
3618  FTYPE a_bsqline[NBIGM];
3619  FTYPE (*bsqline);
3620  FTYPE a_parline[NBIGM];
3621  FTYPE (*parline);
3622  FTYPE a_uparline[NBIGM];
3623  FTYPE (*uparline);
3624  FTYPE a_veffline[NBIGM];
3625  FTYPE (*veffline);
3626  struct of_geom geomdontuse;
3627  struct of_geom *ptrgeom=&geomdontuse;
3628 
3629 
3630 
3631 
3632 
3633  if(doingweno==0){
3634  dualfprintf(fail_file,"Shouldn't be in compute_df_line_new if not doing WENO\n");
3635  myexit(7671515);
3636  }
3637 
3639  //
3640  // Define which quantities (pl) to operate on
3641  //
3643 
3644  // we don't deal with interpolated quantities here, only real primitives
3645  // setup_nprlocalist(whichprimtype,&nprlocalstart,&nprlocalend,nprlocallist,&numprims);
3646 
3647 
3648  // shift 1-D arrays
3649  Pgasline=(FTYPE (*)) (&(a_Pgasline[NBIGBND]));
3650  bsqline=(FTYPE (*)) (&(a_bsqline[NBIGBND]));
3651  parline=(FTYPE (*)) (&(a_parline[NBIGBND]));
3652  uparline=(FTYPE (*)) (&(a_uparline[NBIGBND]));
3653  veffline=(FTYPE (*)) (&(a_veffline[NBIGBND]));
3654 
3655 
3656  // trueijkp->iter should really be passed
3657  di=(trueijkp->iter==1);
3658  dj=(trueijkp->iter==2);
3659  dk=(trueijkp->iter==3);
3660 
3661 
3662  for(num=bs;num<=be;num++){
3663 
3664  // trueijkp->i,trueijkp->j,trueijkp->k,trueijkp->p should really be passed
3665  // needed for EOS in general
3666  locali=trueijkp->i+di*num;
3667  localj=trueijkp->j+dj*num;
3668  localk=trueijkp->k+dk*num;
3669  localloc=trueijkp->p;
3670 
3671 
3672  PALLREALLOOP(pl) myprim[pl] = yprim[pl][0][num];
3673  get_geometry(locali,localj,localk,localloc,ptrgeom);
3674 
3675  // SUPERGODMARK: OPTMARK: Assume prim at CENT for now, which is consistent with other places but not really true in general when averaging primitive locations for old a2c method
3676  // if(interporflux==ENOINTERPTYPE){
3677  get_stateforinterpline(myprim,ptrgeom,&qptr);
3678  // }
3679  // else{
3680  // // then can't be sure quantity is at CENT
3681  // get_state(myprim,ptrgeom,qptr);
3682  // }
3683 
3684 
3685  Pgasline[num]=pressure_rho0_u_simple(locali,localj,localk,localloc,myprim[RHO],myprim[UU]); // just gas contribution, add magnetic below
3686 
3687 #if(VLINEWITHGDETRHO==0)
3688  // Vline[EOMSETMHD][num]=myprim[UU+dir];
3689  // Vline[EOMSETMHD][num]=(qptr->ucon[dir])/(qptr->ucon[TT]); // more correct generally since wanted w.r.t. grid
3690  Vline[EOMSETMHD][num]=qptr->ucon[dir]; // more correct generally since wanted w.r.t. grid (need 4-velocity since consistent with equations of motion)
3691 #else
3692  Vline[EOMSETMHD][num]=(ptrgeom->gdet)*yprim[RHO][0][num]*(qptr->ucon[dir]);
3693 #endif
3694 
3695 #if(RADSHOCKFLAT&&EOMRADTYPE!=EOMRADNONE)
3696 #if(VLINEWITHGDETRHO==0)
3697  Vline[EOMSETRAD][num]=qptr->uradcon[dir];
3698 #else
3699  Vline[EOMSETRAD][num]=(ptrgeom->gdet)*yprim[PRAD0][0][num]*(qptr->uradcon[dir]); // approximate KORALTODO
3700 #endif
3701 #endif
3702 
3703 
3704  // bsq = dot(q.bcon, q.bcov);
3705  bsq = qptr->bsq;
3706  btt = -qptr->bcon[dir] * qptr->bcov[TT];
3707  myPgas = Pgasline[num];
3708  myPtot = myPgas + 0.5*bsq;
3709 
3710  // get total pressure
3711  Pline[EOMSETMHD][num]=myPtot;
3712 #if(RADSHOCKFLAT&&EOMRADTYPE!=EOMRADNONE) // KORAL
3713  // add radiation pressure to total pressure if optically thick
3714  FTYPE tautot[NDIM],tautotmax;
3715  // calcfull_tautot(myprim, ptrgeom, tautot, &tautotmax);
3716  calc_tautot(myprim, ptrgeom, NULL, tautot, &tautotmax); // very accurate tautot not necessary, so use Tgas=Trad assumption in opacity
3717  Pline[EOMSETRAD][num]=(4.0/3.0-1.0)*yprim[PRAD0][0][num]; // radiative pressure
3718  // add radiation pressure to total pressure if optically thick
3719  // KORALNOTE: recall pressure just along diagonal and no velocity in R^\mu_\nu
3720  Pline[EOMSETMHD][num]+=MIN(tautotmax,1.0)*Pline[EOMSETRAD][num];
3721 #endif
3722 
3723  bsqline[num]=fabs(bsq);
3724 
3725  parline[num]=fabs(qptr->ifremoverestplus1ud0elseud0*(qptr->ucon[dir]));
3726 
3727  // note sign comes from u^{dir}, but otherwise correct velocity to go into -T^t_t
3728  veffline[num] = sign(Vline[EOMSETMHD][num])*parline[num]/(sqrt(parline[num])+SMALL); // GODMARK: sqrt expensive, not sure how to avoid
3729 
3730  // upar is really just the non-kinetic part of -T^t_t
3731  // compute it like below so don't repeat expensive plus1ud0 calculation
3732  uparline[num] = fabs( (myprim[UU] + myPgas + bsq )*( fabs(qptr->ucon[dir]) *qptr->ucov[TT]) + myPtot + btt );
3733  }
3734 
3735 
3736 
3737 
3739  //
3740  // can go all the way to edge of grid for C2E for stiffness indicator
3741  //
3743  if(recontype==CVT_C2E){
3744  for(num=bs;num<=be;num++){
3745 
3746  // assignments
3747  PALLREALLOOP(pl) myprim[pl] = yprim[pl][0][num];
3748  bsq=bsqline[num];
3749  par=parline[num];
3750  upar=uparline[num];
3751  veff=veffline[num];
3752 
3753  stifffactor = max(max(fabs(bsq/myprim[RHO]),fabs(bsq/myprim[UU])),fabs(par));
3754  // now turn into something that goes from 0 to 1
3755  // assume moderatively relativistic is ok (stifffactor~3, but higher is sufficiently stiff, so truncate at 3)
3756  stiffindicator[num] = max(min(stifffactor-2.0,1.0),0.0);
3757  }
3758  }
3759 
3760 
3761 
3763  //
3764  // Get df type shock indicator (and also stiffindicator if not C2E)
3765  //
3767  for(num=bs+1;num<=be-1;num++){
3768 
3769  // assignments
3770  PALLREALLOOP(pl) myprim[pl] = yprim[pl][0][num];
3771  bsq=bsqline[num];
3772  par=parline[num];
3773  parl=parline[num-1];
3774  parr=parline[num+1];
3775  upar=uparline[num];
3776  uparl=uparline[num-1];
3777  uparr=uparline[num+1];
3778  veff=veffline[num];
3779  veffl=veffline[num-1];
3780  veffr=veffline[num+1];
3781 
3782 
3783  //The expression is:
3784  // -(rho v^2/2 + u) = rho * ucon * (1+ucov) + gam * u * (ucon*ucov) + (gam-1)*u
3785 
3786  // Original dP/P prescription
3787  //df[num] = ( fabs(uparr-uparl) ) / ( fabs(upar) + DBL_MIN );
3788 
3789  //normalized difference using the central value -- think may lead to problems
3790  //df[num] = (fabs(myprim[RHO]*0.5*(parr-parl)) + fabs(0.5*(uparr-uparl)))/(myprim[RHO]*par + upar);
3791 
3792  //normalized difference without using the central value; normalize the difference by the values themselves
3793  //df[num] = ( fabs(myprim[RHO]*(parr-parl)) + fabs(uparr-uparl) ) / ( fabs(myprim[RHO])*(fabs(parr)+fabs(parl)) + fabs(uparr)+fabs(uparl) + DBL_MIN );
3794 
3795  //divide by one point value, not by the sum of the values themselves
3796  // below is like: rho delta(v^2)/u
3797  // df[num] = ( fabs(myprim[RHO]*(parr-parl)) + fabs(uparr-uparl) ) / ( fabs(myprim[RHO]*par) + fabs(upar) + DBL_MIN );
3798 
3799  // below is now like: rho (delta v)^2/u -- consistent with (\delta Mach)^2 that goes into internal energy error estimate
3800  df[num] = ( fabs(myprim[RHO]*(veffr-veffl)*(veffr-veffl)) + fabs(uparr-uparl) ) / ( fabs(myprim[RHO]*par) + fabs(upar) + SMALL );
3801 
3802 
3803 
3804  // equations become stiff when b^2/rho_0, b^2/u, or (\gamma-1) become sufficiently relativistic
3805  // now also include above-computed df that is for supersonic motion that is stiff
3806  // stifffactor = max(max(max(fabs(bsq/myprim[RHO]),fabs(bsq/myprim[UU])),fabs(par)),fabs(df[num]));
3807 
3808  // JCM: for now don't include df[num] since often true and shouldn't really affect C2E that only uses this right now
3809  // JCM: should directly use df[num] in A2C/C2A to determine if to do those operations or not
3810 
3811  if(recontype==CVT_C2E){
3812  // then already done above
3813  }
3814  else if(recontype==CVT_C2A || recontype==CVT_A2C){
3815  // CHANGINGMARK: can use above for c2e but need df[num] version for a2c/c2a/etc. to avoid ....?
3816  // stifffactor = max(max(fabs(bsq/myprim[RHO]),fabs(bsq/myprim[UU])),fabs(par));
3817  stifffactor = fabs(par);
3818  // now turn into something that goes from 0 to 1
3819  stifffactor1 = max(min(stifffactor-2.0,1.0),0.0);
3820  // stifffactor = max(df[num],stifffactor);
3821  // now turn into something that goes from 0 to 1
3822  stifffactor2 = max(min(fabs(df[num])*30.0,1.0),0.0);
3823 
3824  // now turn into something that goes from 0 to 1
3825  stiffindicator[num]=max(min(max(stifffactor1,stifffactor2),1.0),0.0);
3826 
3827  // stiffindicator[num]=1.0;
3828  }
3829 
3830 
3831 
3832 
3833 #if(0 && FLUXDUMP) // DEBUG: (turn off FLUXDUMP in flux.c)
3834  // trueijkp->i,trueijkp->j,trueijkp->k,trueijkp->p should really be passed
3835  // ijkcurr needed for EOS in general
3836  locali=trueijkp->i+di*num;
3837  localj=trueijkp->j+dj*num;
3838  localk=trueijkp->k+dk*num;
3839  localloc=trueijkp->p;
3840 
3841  GLOBALMACP0A1(fluxdump,locali,localj,localk,4*NPR + (dir-1)*NPR*5 + NPR*0 + RHO)=stifffactor;
3842  GLOBALMACP0A1(fluxdump,locali,localj,localk,4*NPR + (dir-1)*NPR*5 + NPR*0 + UU)=stiffindicator[num];
3843  GLOBALMACP0A1(fluxdump,locali,localj,localk,4*NPR + (dir-1)*NPR*5 + NPR*0 + U1)=veff;
3844 
3845  // dualfprintf(fail_file,"%ld %d :: dir=%d :: ijkcurr=%d %d %d :: %21.15g %21.15g %21.15g\n",nstep,steppart,dir,locali,localj,localk,stifffactor,stiffindicator[num],veff);
3846 
3847 #endif
3848 
3849 
3850 
3851 
3852 
3853  }
3854 
3855  return(0);
3856 }
3857 
3858 
3859 
3860 
3861 
3864 static void assign_eno_result_c2e_multipl(int whichquantity, int recontype, int bs, int ps, int pe, int be, int i, int j, int k, int idel, int jdel, int kdel, FTYPE (*yout)[2][NBIGM], FTYPE (*result0)[NSTORE2][NSTORE3][NPR2INTERP], FTYPE (*result1)[NSTORE2][NSTORE3][NPR2INTERP])
3865 {
3866  // for NUMPRIMLOOP:
3867  int nprlocalstart,nprlocalend;
3868  int nprlocallist[MAXNPR];
3869  int pllocal;
3870  int numprims;
3871  int plstart;
3872  int pl,pliter;
3873  int mypl;
3874  // others:
3875  int l;
3876  int ii,jj,kk;
3877 
3878 
3880  //
3881  // Define which quantities (pl) to operate on
3882  //
3884 
3885  setup_nprlocalist(whichquantity,&nprlocalstart,&nprlocalend,nprlocallist,&numprims);
3886 
3887 
3888 
3889 
3890  // where to place in real final array
3891  // originally had idel*l, jdel*l, kdel*l inside loop. Was more expensive and was one of most expensive lines in code.
3892 
3893  // note that NUMPRIMLOOP(pliter,pl) is inside for(l) loop since (*result)[][][] has larger memory jumps for l than (*yout)[][]
3894  // Small difference, since apparently sometimes run and having pl on outside is faster overall, but % usage is better with pl deep inside.
3895  if(idel==1){
3896  jj=j;
3897  kk=k;
3898  for(l=ps;l<=pe;l++){ // inclusive loop
3899  ii=i+l;
3900  NUMPRIMLOOP(pliter,pl) MACP0A1(result0,ii,jj,kk,pl) = yout[pl][0][l];
3901  NUMPRIMLOOP(pliter,pl) MACP0A1(result1,ii,jj,kk,pl) = yout[pl][1][l];
3902  }
3903  }
3904  else if(jdel==1){
3905  ii=i;
3906  kk=k;
3907  for(l=ps;l<=pe;l++){ // inclusive loop
3908  jj=j+l;
3909  NUMPRIMLOOP(pliter,pl) MACP0A1(result0,ii,jj,kk,pl) = yout[pl][0][l];
3910  NUMPRIMLOOP(pliter,pl) MACP0A1(result1,ii,jj,kk,pl) = yout[pl][1][l];
3911  }
3912  }
3913  else if(kdel==1){
3914  ii=i;
3915  jj=j;
3916  for(l=ps;l<=pe;l++){ // inclusive loop
3917  kk=k+l;
3918  NUMPRIMLOOP(pliter,pl) MACP0A1(result0,ii,jj,kk,pl) = yout[pl][0][l];
3919  NUMPRIMLOOP(pliter,pl) MACP0A1(result1,ii,jj,kk,pl) = yout[pl][1][l];
3920  }
3921  }
3922  else{
3923  dualfprintf(fail_file,"No such idel=%d jdel=%d kdel=%d in assign_eno_result_c2e()\n",idel,jdel,kdel);
3924  myexit(3469836);
3925  }
3926 
3927 
3928  if(recontype!=CVT_C2E){
3929  dualfprintf(fail_file,"assign_eno_result_c2e only handles recontype==CVT_C2E\n");
3930  myexit(26);
3931  }
3932 
3933 }
3934 
3935 
3936 
3939 static void assign_eno_result_c2e(int recontype, int pl, int bs, int ps, int pe, int be, int i, int j, int k, int idel, int jdel, int kdel, FTYPE (*yout)[NBIGM], FTYPE (*result0)[NSTORE2][NSTORE3][NPR2INTERP], FTYPE (*result1)[NSTORE2][NSTORE3][NPR2INTERP])
3940 {
3941  int l;
3942  int ii,jj,kk;
3943 
3944  // where to place in real final array
3945  // originally had idel*l, jdel*l, kdel*l inside loop. Was more expensive and was one of most expensive lines in code.
3946 
3947  if(idel==1){
3948  jj=j;
3949  kk=k;
3950  for(l=ps;l<=pe;l++){ // inclusive loop
3951  ii=i+l;
3952  MACP0A1(result0,ii,jj,kk,pl) = yout[0][l];
3953  MACP0A1(result1,ii,jj,kk,pl) = yout[1][l];
3954  }
3955  }
3956  else if(jdel==1){
3957  ii=i;
3958  kk=k;
3959  for(l=ps;l<=pe;l++){ // inclusive loop
3960  jj=j+l;
3961  MACP0A1(result0,ii,jj,kk,pl) = yout[0][l];
3962  MACP0A1(result1,ii,jj,kk,pl) = yout[1][l];
3963  }
3964  }
3965  else if(kdel==1){
3966  ii=i;
3967  jj=j;
3968  for(l=ps;l<=pe;l++){ // inclusive loop
3969  kk=k+l;
3970  MACP0A1(result0,ii,jj,kk,pl) = yout[0][l];
3971  MACP0A1(result1,ii,jj,kk,pl) = yout[1][l];
3972  }
3973  }
3974  else{
3975  dualfprintf(fail_file,"No such idel=%d jdel=%d kdel=%d in assign_eno_result_c2e()\n",idel,jdel,kdel);
3976  myexit(3469836);
3977  }
3978 
3979 
3980  if(recontype!=CVT_C2E){
3981  dualfprintf(fail_file,"assign_eno_result_c2e only handles recontype==CVT_C2E\n");
3982  myexit(26);
3983  }
3984 
3985 }
3986 
3987 
3990 static void assign_eno_result(int recontype, int pl, int bs, int ps, int pe, int be, int i, int j, int k, int idel, int jdel, int kdel, FTYPE (*yout)[NBIGM], FTYPE (*result0)[NSTORE2][NSTORE3][NPR], FTYPE (*result1)[NSTORE2][NSTORE3][NPR])
3991 {
3992  int l;
3993  int ti,tj,tk;
3994 
3995 
3996  if(recontype==CVT_C2E){// these methods generated 2 output values
3997  for(l=ps;l<=pe;l++){ // inclusive loop
3998  MACP0A1(result0,i+l*idel,j+l*jdel,k+l*kdel,pl) = yout[0][l];
3999  MACP0A1(result1,i+l*idel,j+l*jdel,k+l*kdel,pl) = yout[1][l];
4000  }
4001  }
4002  else if( (recontype==CVT_C2A)||(recontype==CVT_A2C) ){// these methods generated 1 output value
4003  for(l=ps;l<=pe;l++){ // inclusive loop
4004  // result0 can be equal to input p2interpm array since now done with that line and dimensionally split, so never operate in another direction on this quantity
4005 
4006  ti=i+l*idel;
4007  tj=j+l*jdel;
4008  tk=k+l*kdel;
4009 
4010  if(crapdebug==0) MACP0A1(result0,i+l*idel,j+l*jdel,k+l*kdel,pl) = yout[0][l];
4011  else{
4012  if(pl!=3) MACP0A1(result0,i+l*idel,j+l*jdel,k+l*kdel,pl) = yout[0][l];
4013  // MACP0A1(result0,i+l*idel,j+l*jdel,k+l*kdel,pl) = yout[0][l];
4014 
4015 
4016  if(ti<-N1BND || ti>N1+N1BND-1 || tj<-N2BND || tj>N2+N2BND-1 || tk<-N3BND || tk>N3+N3BND-1 || pl<0 || pl>U3){
4017  dualfprintf(fail_file,"OUT OF BOUNDS\n");
4018  }
4019  }
4020 
4021 #if(0)
4022  if(crapdebug && pl==3){
4023  dualfprintf(fail_file,"ti=%d tj=%d tk=%d :: i=%d j=%d k=%d l=%d idel=%d jdel=%d kdel=%d pl=%d\n",i+l*idel,j+l*jdel,k+l*kdel,i,j,k,l,idel,jdel,kdel,pl);
4024  }
4025 #endif
4026  }
4027  }
4028 
4029 }
4030 
4031 
4032 
4033 
4034 
4035 
4036 
4037 
Generated on Fri May 20 2016 15:52:34 for HARM by   doxygen 1.8.3.1