HARM
harm and utilities
 All Data Structures Files Functions Variables Typedefs Macros Pages
set_arrays_multidimen.c
Go to the documentation of this file.
1 // __WORKINGONIT__
2 
14 #include "decs.h"
15 
16 
19 void set_arrays_multidimen()
20 {
21  int i, j, k, pl, pliter, l, m;
22  int pl2;
23  int ii;
24  int jj;
25  int qq;
26  int indexfinalstep,floor,pf, tscale,dtstage;
27  FTYPE valueinit;
28  int dir,interpi,enodebugi;
29  int dimen;
30  int isleftright;
31  struct of_state *myfluxstatetempptr;
32  int firsttimeinloop;
33  extern void set_arrays_multidimen_rad(void);
34 
35 
36 
37 #if(PRODUCTION==0)
38  // initialize things to NAN in order to (hopefully) trigger memory leaks to be noticed
39  valueinit=sqrt(-1.0);
40 #else
41  // avoid this practice for production run since processing nan's slows code and do process some never-initialized/never-used cells for simplicity of code loops
42  valueinit=0.0;
43 #endif
44 
45 
47  //
48  // Basic primitive quantity
49  //
51 
52  GLOBALPOINT(pglobal) = (FTYPE PTRMACP0A1(pglobal,N1M,N2M,N3M,NPR)) (&(BASEMACP0A1(pglobal,N1BND,N2BND,N3BND,0)));
53  FULLLOOP PLOOP(pliter,pl){
54  GLOBALMACP0A1(pglobal,i,j,k,pl) = valueinit;
55  }
56 
57 #if(ANALYTICMEMORY)
58  GLOBALPOINT(panalytic) = (FTYPE PTRMACP0A1(panalytic,N1M,N2M,N3M,NPR)) (&(BASEMACP0A1(panalytic,N1BND,N2BND,N3BND,0)));
59 #if(FIELDSTAGMEM)
60  GLOBALPOINT(pstaganalytic) = (FTYPE PTRMACP0A1(pstaganalytic,N1M,N2M,N3M,NPR)) (&(BASEMACP0A1(pstaganalytic,N1BND,N2BND,N3BND,0)));
61 #endif
62 
63  FULLLOOP PLOOP(pliter,pl){
64  GLOBALMACP0A1(panalytic,i,j,k,pl) = valueinit;
65 #if(FIELDSTAGMEM)
66  GLOBALMACP0A1(pstaganalytic,i,j,k,pl) = valueinit;
67 #endif
68  }
69 #endif
70 
71 
72 #if(NUMPOTHER>0)
73  GLOBALPOINT(pother) = (FTYPE PTRMACP1A0(pother,FILL,N1M,N2M,N3M)) (&(BASEMACP1A0(pother,0,N1BND,N2BND,N3BND)));
74  FULLLOOP for(pl=0;pl<NUMPOTHER;pl++) GLOBALMACP1A0(pother,pl,i,j,k) = -1;
75 #endif
76 
77 
78  // used in fixup.c, higherorder_pointavg.c and kazfulleos.c for compute_Hglobal()
79  GLOBALPOINT(ptemparray) = (FTYPE PTRMACP0A1(ptemparray,N1M,N2M,N3M,NPR)) (&(BASEMACP0A1(ptemparray,N1BND,N2BND,N3BND,0)));
80  FULLLOOP PLOOP(pliter,pl) GLOBALMACP0A1(ptemparray,i,j,k,pl) = valueinit;
81 
82  // used in advance.c to carefully assign conserved quantity per pl
83  GLOBALPOINT(utemparray) = (FTYPE PTRMACP0A1(utemparray,N1M,N2M,N3M,NPR)) (&(BASEMACP0A1(utemparray,N1BND,N2BND,N3BND,0)));
84  FULLLOOP PLOOP(pliter,pl) GLOBALMACP0A1(utemparray,i,j,k,pl) = valueinit;
85 
86 
87 #if(DOEVOLVEMETRIC)
88  GLOBALPOINT(ucumformetric) = (FTYPE PTRMACP0A1(ucumformetric,N1M,N2M,N3M,NPR)) (&(BASEMACP0A1(ucumformetric,N1BND,N2BND,N3BND,0)));
89  FULLLOOP PLOOP(pliter,pl) GLOBALMACP0A1(ucumformetric,i,j,k,pl) = valueinit;
90 #endif
91 
92 
93 
94 #if(STOREFLUXSTATE)
95  GLOBALPOINT(fluxstate) = (struct of_state PTRMACP1A1(fluxstate,FILL,N1M,N2M,N3M,NUMLEFTRIGHT)) (&(BASEMACP1A1(fluxstate,-1,N1BND,N2BND,N3BND,0)));
96  GLOBALPOINT(fluxstatecent) = (struct of_state PTRMACP0A0(fluxstatecent,N1M,N2M,N3M)) (&(BASEMACP0A0(fluxstatecent,N1BND,N2BND,N3BND)));
97 
98  for(isleftright=0;isleftright<NUMLEFTRIGHT+1;isleftright++){
99 
100  firsttimeinloop=1;
101  DIMENLOOP(dimen){
102 
103  if(isleftright==NUMLEFTRIGHT){ // last entry (related to +1 in above isleftright loop) treated as if for centered fluxstatecent
104  if(firsttimeinloop==0) break; // break if dimen iterated since only have fluxstatecent[] that doesn't depend upon dimen
105  else firsttimeinloop=0;
106  }
107  else{
108  // then correct to go over all dimen for fluxstate[dimen]
109  }
110 
111  // loop over i,j,k
112  FULLLOOP{
113 
114  if(isleftright==NUMLEFTRIGHT){ // last entry treated as if for centered fluxstatecent
115  // just way to select centered case so below assignments are in a single code
116  myfluxstatetempptr = &GLOBALMACP0A0(fluxstatecent,i,j,k);
117  }
118  else{
119  myfluxstatetempptr = &GLOBALMACP1A1(fluxstate,dimen,i,j,k,isleftright);
120  }
121 
122 
123 #if(MERGEDC2EA2CMETHOD)
124  PALLLOOP(pl){
125  myfluxstatetempptr->prim[pl]=valueinit;
126  myfluxstatetempptr->EOMFUNCMAC(pl)=valueinit;
127  }
128  DLOOPA(jj){
129  myfluxstatetempptr->Blower[jj]=valueinit;
130  myfluxstatetempptr->vcon[jj]=valueinit;
131  myfluxstatetempptr->gdetBcon[jj]=valueinit;
132  }
133  myfluxstatetempptr->gdet=valueinit;
134  myfluxstatetempptr->overut=valueinit;
135 #endif
136  // always done
137  myfluxstatetempptr->pressure=valueinit;
138  myfluxstatetempptr->entropy=valueinit;
139  myfluxstatetempptr->ifremoverestplus1ud0elseud0=valueinit;
140  DLOOPA(jj){
141  myfluxstatetempptr->ucon[jj]=valueinit;
142  myfluxstatetempptr->ucov[jj]=valueinit;
143  myfluxstatetempptr->bcon[jj]=valueinit;
144  myfluxstatetempptr->bcov[jj]=valueinit;
145  }
146 
147  }
148  }
149  }
150 
151 #endif
152 
153 
154  // below used in fluct.c and as vector potential storage
155  GLOBALPOINT(emf) = (FTYPE PTRMACP1A0(emf,FILL,N1M+SHIFT1,N2M+SHIFT2,N3M+SHIFT3)) (&(BASEMACP1A0(emf,0,N1BND,N2BND,N3BND)));// inner shift still same
156  for(l=0;l<NDIM;l++) FULLLOOP GLOBALMACP1A0(emf,l,i,j,k) = valueinit;
157 
158 
159 #if(FIELDTOTHMEM)
160  // below only used in fluxct.c
161  GLOBALPOINT(vconemf) = (FTYPE PTRMACP0A1(vconemf,N1M,N2M,N3M,NDIM-1)) (&(BASEMACP0A1(vconemf,N1BND,N2BND,N3BND,-U1)));
162  FULLLOOP for(l=U1;l<U1+COMPDIM;l++) GLOBALMACP0A1(vconemf,i,j,k,l) = valueinit;
163 #endif
164 
165 #if(MODIFYEMFORVPOT==MODIFYVPOT || TRACKVPOT>0 || EVOLVEWITHVPOT>0)
166  GLOBALPOINT(vpotarrayglobal) = (FTYPE PTRMACP1A0(vpotarrayglobal,FILL,N1M+SHIFT1,N2M+SHIFT2,N3M+SHIFT3)) (&(BASEMACP1A0(vpotarrayglobal,0,N1BND,N2BND,N3BND)));// inner shift still same
167  for(l=0;l<NUMVPOT;l++) FULLLOOP{
168  GLOBALMACP1A0(vpotarrayglobal,l,i,j,k) = valueinit;
169  }
170 
171 #if(ANALYTICMEMORY)
172  GLOBALPOINT(vpotanalytic) = (FTYPE PTRMACP1A0(vpotanalytic,FILL,N1M+SHIFT1,N2M+SHIFT2,N3M+SHIFT3)) (&(BASEMACP1A0(vpotanalytic,0,N1BND,N2BND,N3BND)));// inner shift still same
173 
174  for(l=0;l<NDIM;l++) FULLLOOP{
175  GLOBALMACP1A0(vpotanalytic,l,i,j,k) = valueinit;
176  }
177 #endif
178 
179 #endif
180 
181 
182 #if(PERCELLDT)
183  GLOBALPOINT(dtijk) = (FTYPE PTRMACP0A1(dtijk,N1M,N2M,N3M,COMPDIM)) (&(BASEMACP0A1(dtijk,N1BND,N2BND,N3BND,-1))); // so access like dtijk[1,2,3]
184  FULLLOOP for(l=1;l<=COMPDIM;l++){
185  GLOBALMACP0A1(dtijk,i,j,k,l) = -1; //valueinit; // SUPERGODMARK: Can improve
186  }
187 #endif
188 
189 
190 #if(STOREWAVESPEEDS>0)
191  GLOBALPOINT(wspeedtemp) = (FTYPE PTRMACP1A1(wspeedtemp,NUMEOMSETS,N1M,N2M,N3M,NUMCS)) (&(BASEMACP1A1(wspeedtemp,0,N1BND,N2BND,N3BND,0)));
192  GLOBALPOINT(wspeed) = (FTYPE PTRMACP3A0(wspeed,FILL,COMPDIM,NUMCS,N1M,N2M,N3M)) (&(BASEMACP3A0(wspeed,0,-1,0,N1BND,N2BND,N3BND))); // shifted so wspeed[1,2,3] accesses the memory
193  FULLLOOP for(qq=0;qq<NUMEOMSETS;qq++) for(l=1;l<=COMPDIM;l++) for(m=0;m<NUMCS;m++) GLOBALMACP3A0(wspeed,qq,l,m,i,j,k) = valueinit;
194 #endif
195 
196 #if(STORESHOCKINDICATOR)
197  GLOBALPOINT(shockindicatorarray) = (FTYPE PTRMACP1A0(shockindicatorarray,NUMSHOCKPLS,N1M,N2M,N3M)) (&(BASEMACP1A0(shockindicatorarray,-SHOCKPLDIR1,N1BND,N2BND,N3BND)));
198  FULLLOOP for(l=SHOCKPLDIR1;l<SHOCKPLDIR1+NUMSHOCKPLS;l++) GLOBALMACP1A0(shockindicatorarray,l,i,j,k) = valueinit;
199 #endif
200 
201 
202 
203 
205  //
206  // TIME-STEPPING
207  //
209 
210  GLOBALPOINT(pk) = (FTYPE PTRMACP1A1(pk,FILL,N1M,N2M,N3M,NPR)) (&(BASEMACP1A1(pk,0,N1BND,N2BND,N3BND,0)));
211  GLOBALPOINT(uinitialglobal) = (FTYPE PTRMACP0A1(uinitialglobal,N1M,N2M,N3M,NPR)) (&(BASEMACP0A1(uinitialglobal,N1BND,N2BND,N3BND,0)));
212  GLOBALPOINT(ulastglobal) = (FTYPE PTRMACP0A1(ulastglobal,N1M,N2M,N3M,NPR)) (&(BASEMACP0A1(ulastglobal,N1BND,N2BND,N3BND,0)));
213  GLOBALPOINT(unewglobal) = (FTYPE PTRMACP0A1(unewglobal,N1M,N2M,N3M,NPR)) (&(BASEMACP0A1(unewglobal,N1BND,N2BND,N3BND,0)));
214  GLOBALPOINT(dUgeomarray) = (FTYPE PTRMACP0A1(dUgeomarray,N1M,N2M,N3M,NPR)) (&(BASEMACP0A1(dUgeomarray,N1BND,N2BND,N3BND,0)));
215  FULLLOOP PLOOP(pliter,pl){
216  DTSTAGELOOP(dtstage) GLOBALMACP1A1(pk,dtstage,i,j,k,pl) = valueinit;
217  GLOBALMACP0A1(uinitialglobal,i,j,k,pl) = valueinit;
218  GLOBALMACP0A1(ulastglobal,i,j,k,pl) = valueinit;
219  if(FULLOUTPUT==0) GLOBALMACP0A1(unewglobal,i,j,k,pl) = valueinit;
220  else GLOBALMACP0A1(unewglobal,i,j,k,pl) = 0;
221  GLOBALMACP0A1(dUgeomarray,i,j,k,pl) = valueinit;
222  }
223 
224 
225 #if(HIGHERORDERMEM||FIELDSTAGMEM) // upoint needed for FV method and STAG for all methods
226  GLOBALPOINT(upointglobal) = (FTYPE PTRMACP0A1(upointglobal,N1M,N2M,N3M,NPR)) (&(BASEMACP0A1(upointglobal,N1BND,N2BND,N3BND,0)));
227  FULLLOOP PLOOP(pliter,pl) GLOBALMACP0A1(upointglobal,i,j,k,pl) = valueinit;
228 
229  GLOBALPOINT(upointglobaluf) = (FTYPE PTRMACP0A1(upointglobaluf,N1M,N2M,N3M,NPR)) (&(BASEMACP0A1(upointglobaluf,N1BND,N2BND,N3BND,0)));
230  FULLLOOP PLOOP(pliter,pl) GLOBALMACP0A1(upointglobaluf,i,j,k,pl) = valueinit;
231 
232  GLOBALPOINT(oldufstore) = (FTYPE PTRMACP0A1(oldufstore,N1M,N2M,N3M,NPR)) (&(BASEMACP0A1(oldufstore,N1BND,N2BND,N3BND,0)));
233  FULLLOOP PLOOP(pliter,pl) GLOBALMACP0A1(oldufstore,i,j,k,pl) = valueinit;
234 #endif
235 
236 
237 
239  //
240  // SPATIAL INTERPOLATION
241  //
243 
244  // Note that PTR does not have N?M+SHIFT? unlike BASE pointer. This is so flux_ct() can avoid extra loops and still not seg fault when accessing beyind F?'s normal range.
245  // Since we want extra space at *bottom* (So can access (e.g.) F1[-N1BND-1] without seg faulting)
246  // For FLUXCD, want extra space at top, so in general add 2 extra spaces (one at bottom and one at top)
247  // Thus, only do SHIFT+NBND for each since that shifts bottom enough, and top space will be there.
248 #if(N1>1)
249  GLOBALPOINT(F1) = (FTYPE PTRMACP0A1(F1,N1M,N2M,N3M,NPR+NSPECIAL)) (&(BASEMACP0A1(F1,SHIFT1+N1BND,SHIFT2+N2BND,SHIFT3+N3BND,0)));
250  FULLLOOP PLOOP(pliter,pl){
251  GLOBALMACP0A1(F1,i,j,k,pl) = valueinit;
252  }
253 #endif
254 #if(N2>1)
255  GLOBALPOINT(F2) = (FTYPE PTRMACP0A1(F2,N1M,N2M,N3M,NPR+NSPECIAL)) (&(BASEMACP0A1(F2,SHIFT1+N1BND,SHIFT2+N2BND,SHIFT3+N3BND,0)));
256  FULLLOOP PLOOP(pliter,pl){
257  GLOBALMACP0A1(F2,i,j,k,pl) = valueinit;
258  }
259 #endif
260 #if(N3>1)
261  GLOBALPOINT(F3) = (FTYPE PTRMACP0A1(F3,N1M,N2M,N3M,NPR+NSPECIAL)) (&(BASEMACP0A1(F3,SHIFT1+N1BND,SHIFT2+N2BND,SHIFT3+N3BND,0)));
262  FULLLOOP PLOOP(pliter,pl){
263  GLOBALMACP0A1(F3,i,j,k,pl) = valueinit;
264  }
265 #endif
266 
267 #if(SPLITMAEMMEM)
268 #if(N1>1)
269  GLOBALPOINT(F1EM) = (FTYPE PTRMACP0A1(F1EM,N1M,N2M,N3M,NPR+NSPECIAL)) (&(BASEMACP0A1(F1EM,SHIFT1+N1BND,SHIFT2+N2BND,SHIFT3+N3BND,0)));
270  FULLLOOP PLOOP(pliter,pl){
271  GLOBALMACP0A1(F1EM,i,j,k,pl) = valueinit;
272  }
273 #endif
274 #if(N2>1)
275  GLOBALPOINT(F2EM) = (FTYPE PTRMACP0A1(F2EM,N1M,N2M,N3M,NPR+NSPECIAL)) (&(BASEMACP0A1(F2EM,SHIFT1+N1BND,SHIFT2+N2BND,SHIFT3+N3BND,0)));
276  FULLLOOP PLOOP(pliter,pl){
277  GLOBALMACP0A1(F2EM,i,j,k,pl) = valueinit;
278  }
279 #endif
280 #if(N3>1)
281  GLOBALPOINT(F3EM) = (FTYPE PTRMACP0A1(F3EM,N1M,N2M,N3M,NPR+NSPECIAL)) (&(BASEMACP0A1(F3EM,SHIFT1+N1BND,SHIFT2+N2BND,SHIFT3+N3BND,0)));
282  FULLLOOP PLOOP(pliter,pl){
283  GLOBALMACP0A1(F3EM,i,j,k,pl) = valueinit;
284  }
285 #endif
286 #endif
287 
288 
289 
290 #if(SPLITNPR||FIELDSTAGMEM)
291  GLOBALPOINT(gp_l) = (FTYPE PTRMACP1A1(gp_l,FILL,N1M,N2M,N3M,NPR2INTERP)) (&(BASEMACP1A1(gp_l,-1,N1BND,N2BND,N3BND,0)));
292  GLOBALPOINT(gp_r) = (FTYPE PTRMACP1A1(gp_r,FILL,N1M,N2M,N3M,NPR2INTERP)) (&(BASEMACP1A1(gp_r,-1,N1BND,N2BND,N3BND,0)));
293 
294  FULLLOOP PINTERPLOOP(pliter,pl){
295  for(l=1;l<=3;l++){
296  GLOBALMACP1A1(gp_l,l,i,j,k,pl) = valueinit;
297  GLOBALMACP1A1(gp_r,l,i,j,k,pl) = valueinit;
298  }
299  }
300 
301 #endif
302 
303 
304  GLOBALPOINT(pleft) = (FTYPE PTRMACP0A1(pleft,N1M,N2M,N3M,NPR2INTERP)) (&(BASEMACP0A1(pleft,N1BND,N2BND,N3BND,0)));
305  GLOBALPOINT(pright) = (FTYPE PTRMACP0A1(pright,N1M,N2M,N3M,NPR2INTERP)) (&(BASEMACP0A1(pright,N1BND,N2BND,N3BND,0)));
306  GLOBALPOINT(prc) = (FTYPE PTRMACP0A1(prc,N1M,N2M,N3M,NPR2INTERP)) (&(BASEMACP0A1(prc,N1BND,N2BND,N3BND,0)));
307 
308  FULLLOOP PINTERPLOOP(pliter,pl){
309  GLOBALMACP0A1(pleft,i,j,k,pl) = valueinit;
310  GLOBALMACP0A1(pright,i,j,k,pl) = valueinit;
311  GLOBALMACP0A1(prc,i,j,k,pl) = valueinit;
312  }
313 
314 
315 
316 
317 #if(FIELDSTAGMEM)
318  // GLOBALPOINT(wspeedcorn) = (FTYPE PTRMACP2A0(wspeedcorn,FILL,NUMCS,N1M,N2M,N3M)) (&(BASEMACP2A0(wspeedcorn,-1,0,N1BND,N2BND,N3BND))); // shifted so wspeedcorn[1,2,3] accesses the memory
319  GLOBALPOINT(pstagglobal) = (FTYPE PTRMACP0A1(pstagglobal,N1M,N2M,N3M,NPR)) (&(BASEMACP0A1(pstagglobal,N1BND,N2BND,N3BND,0)));
320 
321  // -B1 and -U1 are so pbcorninterp[][B1] accesses 0th element of original memory (same for pvcorninterp)
322  // GLOBALPOINT(pbcorninterp) = (FTYPE PTRMACP3A0(pbcorninterp,FILL,COMPDIM,NUMCS,N1M+SHIFT1,N2M+SHIFT2,N3M+SHIFT3)) (&(BASEMACP3A0(pbcorninterp,-1,-B1,0,N1BND,N2BND,N3BND)));
323  GLOBALPOINT(pvbcorninterp) = (FTYPE PTRMACP1A3(pvbcorninterp,COMPDIM,N1M+SHIFT1,N2M+SHIFT2,N3M+SHIFT3,COMPDIM,NUMCS+1,NUMCS)) (&(BASEMACP1A3(pvbcorninterp,-1,N1BND,N2BND,N3BND,-1,0,0))); // only shift by -1 since holds both U1+dir-1 and B1+dir-1 for dir=1,2,3 that will now just be accessed via "dir" alone
324 
325  GLOBALPOINT(geomcornglobal) = (FTYPE PTRMACP1A0(geomcornglobal,FILL,N1M+SHIFT1,N2M+SHIFT2,N3M+SHIFT3)) (&(BASEMACP1A0(geomcornglobal,-1,N1BND,N2BND,N3BND))); // geomcornglobal[1,2,3 for CORN1,2,3]
326 
327 #if(HIGHERORDERMEM)
328  // using NPR so can use normal functions
329  GLOBALPOINT(Bhatglobal) = (FTYPE PTRMACP0A1(Bhatglobal,N1M,N2M,N3M,NPR)) (&(BASEMACP0A1(Bhatglobal,N1BND,N2BND,N3BND,0)));
330  FULLLOOP PLOOP(pliter,pl){
331  GLOBALMACP0A1(Bhatglobal,i,j,k,pl)=valueinit;
332  }
333 
334 #if(ANALYTICMEMORY)
335  GLOBALPOINT(Bhatanalytic) = (FTYPE PTRMACP0A1(Bhatanalytic,N1M,N2M,N3M,NPR)) (&(BASEMACP0A1(Bhatanalytic,N1BND,N2BND,N3BND,0)));
336  FULLLOOP PLOOP(pliter,pl){
337  GLOBALMACP0A1(Bhatanalytic,i,j,k,pl)=valueinit;
338  }
339 #endif
340 
341 #endif
342 
343 
344  FULLLOOP PALLLOOP(pl) GLOBALMACP0A1(pstagglobal,i,j,k,pl) = valueinit;
345 
346  // corner things that have extra at top
347  FULLLOOPP1{
348  for(pl2=1;pl2<=COMPDIM;pl2++) for(pl=1;pl<=COMPDIM;pl++) for(m=0;m<NUMCS+1;m++) for(l=0;l<NUMCS;l++) GLOBALMACP1A3(pvbcorninterp,pl2,i,j,k,pl,m,l)=valueinit;
349  // for(pl2=1;pl2<=COMPDIM;pl2++) for(pl=B1;pl<=B3;pl++) for(l=0;l<NUMCS;l++) GLOBALMACP3A0(pbcorninterp,pl2,pl,l,i,j,k)=valueinit;
350 
351  for(pl2=1;pl2<=COMPDIM;pl2++) GLOBALMACP1A0(geomcornglobal,pl2,i,j,k)=valueinit;
352  }
353 
354 #endif
355 
356 
357 
358 
360  //
361  // OLD SPATIAL INTERPOLATION (and new way to store shock indicator for paraflag or paraline)
362  //
364 
365 #if(DODQMEMORY||STORESHOCKINDICATOR)
366 #if(N1>1)
367  GLOBALPOINT(dq1) = (FTYPE PTRMACP0A1(dq1,N1M,N2M,N3M,NPR2INTERP)) (&(BASEMACP0A1(dq1,N1BND,N2BND,N3BND,0)));
368  FULLLOOP PINTERPLOOP(pliter,pl) GLOBALMACP0A1(dq1,i,j,k,pl) = valueinit;
369 #endif
370 #if(N2>1)
371  GLOBALPOINT(dq2) = (FTYPE PTRMACP0A1(dq2,N1M,N2M,N3M,NPR2INTERP)) (&(BASEMACP0A1(dq2,N1BND,N2BND,N3BND,0)));
372  FULLLOOP PINTERPLOOP(pliter,pl) GLOBALMACP0A1(dq2,i,j,k,pl) = valueinit;
373 #endif
374 #if(N3>1)
375  GLOBALPOINT(dq3) = (FTYPE PTRMACP0A1(dq3,N1M,N2M,N3M,NPR2INTERP)) (&(BASEMACP0A1(dq3,N1BND,N2BND,N3BND,0)));
376  FULLLOOP PINTERPLOOP(pliter,pl) GLOBALMACP0A1(dq3,i,j,k,pl) = valueinit;
377 #endif
378 #endif
379 
380 
381 
382 
383 
384 
386  //
387  // HIGHER ORDER STUFF
388  //
390 
391 #if( HIGHERORDERMEM )
392  GLOBALPOINT(fluxvectemp) = (FTYPE PTRMACP0A1(fluxvectemp,N1M,N2M,N3M,NPR+NSPECIAL)) (&(BASEMACP0A1(fluxvectemp,N1BND,N2BND,N3BND,0)));
393  GLOBALPOINT(Fa) = (FTYPE PTRMACP0A1(Fa,N1M,N2M,N3M,NPR+NSPECIAL)) (&(BASEMACP0A1(Fa,N1BND,N2BND,N3BND,0)));
394  GLOBALPOINT(Fb) = (FTYPE PTRMACP0A1(Fb,N1M,N2M,N3M,NPR+NSPECIAL)) (&(BASEMACP0A1(Fb,N1BND,N2BND,N3BND,0)));
395  trifprintf("Allocating Fa/Fb for UNSPLIT (ALL NOW) flux method\n");
396  GLOBALPOINT(stencilvartemp) = (FTYPE PTRMACP0A1(stencilvartemp,N1M,N2M,N3M,NPR)) (&(BASEMACP0A1(stencilvartemp,N1BND,N2BND,N3BND,0)));
397 
398  // Note that this still references Fa,Fb to save memory
399  GLOBALPOINT(a2cin) = (FTYPE PTRMACP0A1(Fa,N1M,N2M,N3M,NPR+NSPECIAL)) (&(BASEMACP0A1(Fa,N1BND,N2BND,N3BND,0)));
400  GLOBALPOINT(a2cout) = (FTYPE PTRMACP0A1(Fb,N1M,N2M,N3M,NPR+NSPECIAL)) (&(BASEMACP0A1(Fb,N1BND,N2BND,N3BND,0)));
401  trifprintf("Allocating a2cin/a2cout for UNSPLIT (ALL NOW) a2c method\n");
402 
403 #endif
404 
405 
406 
407 
408 
409 
410 #if(FLUXDUMP)
411  GLOBALPOINT(fluxdump)=(FTYPE PTRMACP0A1(fluxdump,N1M,N2M,N3M,NUMFLUXDUMP)) (&(BASEMACP0A1(fluxdump,N1BND,N2BND,N3BND,0)));
412 
413  // normal things
414  FULLLOOP for(pl=0;pl<NUMFLUXDUMP;pl++) GLOBALMACP0A1(fluxdump,i,j,k,pl)=0.0; // not used for evolution -- just dumping -- so ok to ignore if leaking(?)
415 
416 #endif
417 
418 
419 
421  //
422  // DEBUG STUFF USUALLY ON
423  //
425 
426  GLOBALPOINT(pflag) = (PFTYPE PTRMACP0A1(pflag,N1M,N2M,N3M,NUMPFLAGS)) (&(BASEMACP0A1(pflag,N1BND,N2BND,N3BND,0)));
427  FULLLOOP PFLAGLOOP(pf) GLOBALMACP0A1(pflag,i,j,k,pf) = NANPFLAG;
428 
429  GLOBALPOINT(pflagfailorig) = (PFTYPE PTRMACP0A1(pflagfailorig,N1M,N2M,N3M,NUMFAILPFLAGS)) (&(BASEMACP0A1(pflagfailorig,N1BND,N2BND,N3BND,0)));
430  FULLLOOP FAILPFLAGLOOP(pf) GLOBALMACP0A1(pflagfailorig,i,j,k,pf) = NANPFLAG;
431 
432 #if(DODEBUG)
433  GLOBALPOINT(failfloorcount) = (CTYPE PTRMACP0A3(failfloorcount,N1M,N2M,N3M,2,NUMTSCALES,NUMFAILFLOORFLAGS)) (&(BASEMACP0A3(failfloorcount,N1BND,N2BND,N3BND,0,0,0)));
434  FULLLOOP FAILFLOORLOOP(indexfinalstep,tscale,floor) GLOBALMACP0A3(failfloorcount,i,j,k,indexfinalstep,tscale,floor)=valueinit;
435 #endif
436 
437 #if(DOFLOORDIAG)
438  GLOBALPOINT(failfloordu) = (FTYPE PTRMACP0A1(failfloordu,N1M,N2M,N3M,NPR)) (&(BASEMACP0A1(failfloordu,N1BND,N2BND,N3BND,0)));
439  FULLLOOP{
440  PALLLOOP(pl) GLOBALMACP0A1(failfloordu,i,j,k,pl)=valueinit;
441  }
442 #endif
443 #if(DODISSMEASURE)
444  GLOBALPOINT(dissmeasurearray) = (FTYPE PTRMACP0A1(dissmeasurearray,N1M,N2M,N3M,NSPECIAL+1+3*2)) (&(BASEMACP0A1(dissmeasurearray,N1BND,N2BND,N3BND,0)));
445  FULLLOOP{
446  // for(pl=0;pl<NSPECIAL+1+3*2;pl++) GLOBALMACP0A1(dissmeasurearray,i,j,k,pl)=valueinit;
447  for(pl=0;pl<NSPECIAL+1+3*2;pl++) GLOBALMACP0A1(dissmeasurearray,i,j,k,pl)=0.0; // __WORKINGONIT__ For now so can use past measure to help set Fi and Firad
448  // for dump to be clean for unused things
449  if(N1==1){
450  dir=1; pl=NSPECIAL+1+dir-1; GLOBALMACP0A1(dissmeasurearray,i,j,k,pl)=0.0;
451  dir=1; pl=NSPECIAL+1+3+dir-1; GLOBALMACP0A1(dissmeasurearray,i,j,k,pl)=0.0;
452  }
453  if(N2==1){
454  dir=2; pl=NSPECIAL+1+dir-1; GLOBALMACP0A1(dissmeasurearray,i,j,k,pl)=0.0;
455  dir=2; pl=NSPECIAL+1+3+dir-1; GLOBALMACP0A1(dissmeasurearray,i,j,k,pl)=0.0;
456  }
457  if(N3==1){
458  dir=3; pl=NSPECIAL+1+dir-1; GLOBALMACP0A1(dissmeasurearray,i,j,k,pl)=0.0;
459  dir=3; pl=NSPECIAL+1+3+dir-1; GLOBALMACP0A1(dissmeasurearray,i,j,k,pl)=0.0;
460  }
461  }
462 #endif
463 
465  //
466  // other diagnostics
467  //
469 
470 
471 #if(DODISS)
472  GLOBALPOINT(dissfunpos) = (FTYPE PTRMACP0A1(dissfunpos,N1M,N2M,N3M,NUMDISSFUNPOS)) (&(BASEMACP0A1(dissfunpos,N1BND,N2BND,N3BND,0)));
473 #endif
474 
475 
476 #if(CALCFARADAYANDCURRENTS)
477  // this faraday needed for current calculation
478  GLOBALPOINT(cfaraday) = (FTYPE PTRMACP0A2(cfaraday,N1M,N2M,N3M,NUMCURRENTSLOTS,3)) (&(BASEMACP0A2(cfaraday,N1BND,N2BND,N3BND,0,0)));
479  GLOBALPOINT(fcon) = (FTYPE PTRMACP0A1(fcon,N1M,N2M,N3M,NUMFARADAY)) (&(BASEMACP0A1(fcon,N1BND,N2BND,N3BND,0)));
480  GLOBALPOINT(jcon) = (FTYPE PTRMACP0A1(jcon,N1M,N2M,N3M,NDIM)) (&(BASEMACP0A1(jcon,N1BND,N2BND,N3BND,0)));
481 
482  FULLLOOP{
483  for(pl=0;pl<NUMCURRENTSLOTS;pl++) for(l=0;l<3;l++){
484 #if(DOGRIDSECTIONING)
485  // GODMARK: if grid moves while computing current, then near boundary current will be undefined for time derivative.
486  // Could compute J^i using substeps, but harder to do.
487  // So for now just assume border-region of current is poorly computed and just avoid nan's
488  GLOBALMACP0A2(cfaraday,i,j,k,pl,l)=0.0;
489 #else
490  GLOBALMACP0A2(cfaraday,i,j,k,pl,l)=valueinit;
491 #endif
492  }
493  for(pl=0;pl<NUMFARADAY;pl++){
494  GLOBALMACP0A1(fcon,i,j,k,pl)=valueinit;
495  }
496  for(pl=0;pl<NDIM;pl++){
497  GLOBALMACP0A1(jcon,i,j,k,pl)=valueinit;
498  }
499  }
500 
501 #endif
502 
503 
505  //
506  // AVG diagnostics
507  //
509 
510  // assume time average stuff gets zeroed in avg routine
511 #if(DOAVG)
512  GLOBALPOINT(normalvarstavg) = (FTYPE PTRMACP0A1(normalvarstavg,N1M,N2M,N3M,NUMNORMDUMP)) (&(BASEMACP0A1(normalvarstavg,N1BND,N2BND,N3BND,0)));
513  GLOBALPOINT(anormalvarstavg) = (FTYPE PTRMACP0A1(anormalvarstavg,N1M,N2M,N3M,NUMNORMDUMP)) (&(BASEMACP0A1(anormalvarstavg,N1BND,N2BND,N3BND,0)));
514 
515  GLOBALPOINT(fcontavg) = (FTYPE PTRMACP0A1(fcontavg,N1M,N2M,N3M,NUMFARADAY)) (&(BASEMACP0A1(fcontavg,N1BND,N2BND,N3BND,0)));
516  GLOBALPOINT(fcovtavg) = (FTYPE PTRMACP0A1(fcovtavg,N1M,N2M,N3M,NUMFARADAY)) (&(BASEMACP0A1(fcovtavg,N1BND,N2BND,N3BND,0)));
517 
518  GLOBALPOINT(afcontavg) = (FTYPE PTRMACP0A1(afcontavg,N1M,N2M,N3M,NUMFARADAY)) (&(BASEMACP0A1(afcontavg,N1BND,N2BND,N3BND,0)));
519  GLOBALPOINT(afcovtavg) = (FTYPE PTRMACP0A1(afcovtavg,N1M,N2M,N3M,NUMFARADAY)) (&(BASEMACP0A1(afcovtavg,N1BND,N2BND,N3BND,0)));
520 
521  GLOBALPOINT(massfluxtavg) = (FTYPE PTRMACP0A1(massfluxtavg,N1M,N2M,N3M,NDIM)) (&(BASEMACP0A1(massfluxtavg,N1BND,N2BND,N3BND,0)));
522  GLOBALPOINT(amassfluxtavg) = (FTYPE PTRMACP0A1(amassfluxtavg,N1M,N2M,N3M,NDIM)) (&(BASEMACP0A1(amassfluxtavg,N1BND,N2BND,N3BND,0)));
523 
524  GLOBALPOINT(othertavg) = (FTYPE PTRMACP0A1(othertavg,N1M,N2M,N3M,NUMOTHER)) (&(BASEMACP0A1(othertavg,N1BND,N2BND,N3BND,0)));
525  GLOBALPOINT(aothertavg) = (FTYPE PTRMACP0A1(aothertavg,N1M,N2M,N3M,NUMOTHER)) (&(BASEMACP0A1(aothertavg,N1BND,N2BND,N3BND,0)));
526 
527 #if(CALCFARADAYANDCURRENTS)
528  GLOBALPOINT(jcontavg) = (FTYPE PTRMACP0A1(jcontavg,N1M,N2M,N3M,NDIM)) (&(BASEMACP0A1(jcontavg,N1BND,N2BND,N3BND,0)));
529  GLOBALPOINT(jcovtavg) = (FTYPE PTRMACP0A1(jcovtavg,N1M,N2M,N3M,NDIM)) (&(BASEMACP0A1(jcovtavg,N1BND,N2BND,N3BND,0)));
530 
531  GLOBALPOINT(ajcontavg) = (FTYPE PTRMACP0A1(ajcontavg,N1M,N2M,N3M,NDIM)) (&(BASEMACP0A1(ajcontavg,N1BND,N2BND,N3BND,0)));
532  GLOBALPOINT(ajcovtavg) = (FTYPE PTRMACP0A1(ajcovtavg,N1M,N2M,N3M,NDIM)) (&(BASEMACP0A1(ajcovtavg,N1BND,N2BND,N3BND,0)));
533 #endif
534 
535  GLOBALPOINT(tudtavg) = (FTYPE PTRMACP0A1(tudtavg,N1M,N2M,N3M,NUMSTRESSTERMS)) (&(BASEMACP0A1(tudtavg,N1BND,N2BND,N3BND,0)));
536  GLOBALPOINT(atudtavg) = (FTYPE PTRMACP0A1(atudtavg,N1M,N2M,N3M,NUMSTRESSTERMS)) (&(BASEMACP0A1(atudtavg,N1BND,N2BND,N3BND,0)));
537 #endif
538 
539 
540 
541 
542 
543  /* grid functions */
544  // GODMARK: for axisymmetric space-times, may want to keep metric functions 2D to save memory
545 
546  // these have 1 extra value on outer edge. Shift for real pointer no different
547 #if(NEWMETRICSTORAGE)
548  // new way
549  GLOBALPOINT(gdetgeom) = (struct of_gdetgeom PTRMETMACP0A1(gdetgeom,N1M+SHIFT1,N2M+SHIFT2,N3M+SHIFT3,NPG)) (&(BASEMETMACP0A1(gdetgeom,N1BND,N2BND,N3BND,0)));
550 
551  GLOBALPOINT(gdetgeomnormal) = (struct of_gdetgeom PTRMETMACP1A0(gdetgeomnormal,NPG,N1M+SHIFT1,N2M+SHIFT2,N3M+SHIFT3)) (&(BASEMETMACP1A0(gdetgeomnormal,0,N1BND,N2BND,N3BND)));
552 
553  GLOBALPOINT(compgeom) = (struct of_compgeom PTRMETMACP1A0(compgeom,FILL,N1M+SHIFT1,N2M+SHIFT2,N3M+SHIFT3)) (&(BASEMETMACP1A0(compgeom,0,N1BND,N2BND,N3BND)));
554 #if(DOEVOLVEMETRIC)
555  GLOBALPOINT(compgeomlast) = (struct of_compgeom PTRMETMACP1A0(compgeomlast,FILL,N1M+SHIFT1,N2M+SHIFT2,N3M+SHIFT3)) (&(BASEMETMACP1A0(compgeomlast,0,N1BND,N2BND,N3BND)));
556 #endif
557 
558 #else //else if old way
559  GLOBALPOINT(gcon) = (FTYPE PTRMETMACP1A1(gcon,FILL,N1M+SHIFT1,N2M+SHIFT2,N3M+SHIFT3,SYMMATRIXNDIM)) (&(BASEMETMACP1A1(gcon,0,N1BND,N2BND,N3BND,0)));
560  GLOBALPOINT(gcov) = (FTYPE PTRMETMACP1A1(gcov,FILL,N1M+SHIFT1,N2M+SHIFT2,N3M+SHIFT3,SYMMATRIXNDIM)) (&(BASEMETMACP1A1(gcov,0,N1BND,N2BND,N3BND,0)));
561  GLOBALPOINT(gcovpert) = (FTYPE PTRMETMACP1A1(gcovpert,FILL,N1M+SHIFT1,N2M+SHIFT2,N3M+SHIFT3,NDIM)) (&(BASEMETMACP1A1(gcovpert,0,N1BND,N2BND,N3BND,0)));
562  GLOBALPOINT(gdet) = (FTYPE PTRMETMACP1A0(gdet,FILL,N1M+SHIFT1,N2M+SHIFT2,N3M+SHIFT3)) (&(BASEMETMACP1A0(gdet,0,N1BND,N2BND,N3BND)));
563 
564 #if(WHICHEOM!=WITHGDET)
565  GLOBALPOINT(eomfunc) = (FTYPE PTRMETMACP1A1(eomfunc,FILL,N1M+SHIFT1,N2M+SHIFT2,N3M+SHIFT3,NPR)) (&(BASEMETMACP1A1(eomfunc,0,N1BND,N2BND,N3BND,0)));
566 #endif
567 #if(GDETVOLDIFF)
568  GLOBALPOINT(gdetvol) = (FTYPE PTRMETMACP1A0(gdetvol,FILL,N1M+SHIFT1,N2M+SHIFT2,N3M+SHIFT3)) (&(BASEMETMACP1A0(gdetvol,0,N1BND,N2BND,N3BND)));
569 #endif
570  GLOBALPOINT(alphalapse) = (FTYPE PTRMETMACP1A0(alphalapse,FILL,N1M+SHIFT1,N2M+SHIFT2,N3M+SHIFT3)) (&(BASEMETMACP1A0(alphalapse,0,N1BND,N2BND,N3BND)));
571  GLOBALPOINT(betasqoalphasq) = (FTYPE PTRMETMACP1A0(betasqoalphasq,FILL,N1M+SHIFT1,N2M+SHIFT2,N3M+SHIFT3)) (&(BASEMETMACP1A0(betasqoalphasq,0,N1BND,N2BND,N3BND)));
572  GLOBALPOINT(beta) = (FTYPE PTRMETMACP1A1(beta,FILL,N1M+SHIFT1,N2M+SHIFT2,N3M+SHIFT3,NDIM)) (&(BASEMETMACP1A1(beta,0,N1BND,N2BND,N3BND,0)));
573 
574 #if(DOEVOLVEMETRIC)
575  GLOBALPOINT(gcovlast) = (FTYPE PTRMETMACP1A1(gcovlast,FILL,N1M+SHIFT1,N2M+SHIFT2,N3M+SHIFT3,SYMMATRIXNDIM)) (&(BASEMETMACP1A1(gcovlast,0,N1BND,N2BND,N3BND,0)));
576  GLOBALPOINT(gcovpertlast) = (FTYPE PTRMETMACP1A1(gcovpertlast,FILL,N1M+SHIFT1,N2M+SHIFT2,N3M+SHIFT3,NDIM)) (&(BASEMETMACP1A1(gcovpertlast,0,N1BND,N2BND,N3BND,0)));
577  GLOBALPOINT(alphalapselast) = (FTYPE PTRMETMACP1A0(alphalapselast,FILL,N1M+SHIFT1,N2M+SHIFT2,N3M+SHIFT3)) (&(BASEMETMACP1A0(alphalapselast,0,N1BND,N2BND,N3BND)));
578 #endif
579 
580 #endif // end if old way
581 
582 
583 #if(DOSTOREPOSITIONDATA)
584  GLOBALPOINT(dxdxpstore) = (FTYPE PTRMETMACP1A2(dxdxpstore,FILL,N1M+SHIFT1,N2M+SHIFT2,N3M+SHIFT3,NDIM,NDIM)) (&(BASEMETMACP1A2(dxdxpstore,0,N1BND,N2BND,N3BND,0,0)));
585  GLOBALPOINT(idxdxpstore) = (FTYPE PTRMETMACP1A2(idxdxpstore,FILL,N1M+SHIFT1,N2M+SHIFT2,N3M+SHIFT3,NDIM,NDIM)) (&(BASEMETMACP1A2(idxdxpstore,0,N1BND,N2BND,N3BND,0,0)));
586 
587  GLOBALPOINT(Xstore) = (FTYPE PTRMACP1A1(Xstore,FILL,N1M+SHIFT1*3,N2M+SHIFT2*3,N3M+SHIFT3*3,NDIM)) (&(BASEMACP1A1(Xstore,0,N1BND+SHIFT1,N2BND+SHIFT2,N3BND+SHIFT3,0))); // Xstore is not reduced if CARTMINKMETRIC
588  GLOBALPOINT(Vstore) = (FTYPE PTRMACP1A1(Vstore,FILL,N1M+SHIFT1*3,N2M+SHIFT2*3,N3M+SHIFT3*3,NDIM)) (&(BASEMACP1A1(Vstore,0,N1BND+SHIFT1,N2BND+SHIFT2,N3BND+SHIFT3,0))); // Vstore is not reduced if CARTMINKMETRIC
589 #endif
590 
591 
592  // rest are always located at CENT
593  GLOBALPOINT(conn) = (FTYPE PTRMETMACP0A3(conn,N1M,N2M,N3M,NDIM,NDIM,NDIM)) (&(BASEMETMACP0A3(conn,N1BND,N2BND,N3BND,0,0,0)));
594  GLOBALPOINT(conn2) = (FTYPE PTRMETMACP0A1(conn2,N1M,N2M,N3M,NDIM)) (&(BASEMETMACP0A1(conn2,N1BND,N2BND,N3BND,0)));
595 
596 #if(VOLUMEDIFF)
597  GLOBALPOINT(idxvol) = (FTYPE PTRMETMACP0A1(idxvol,N1M,N2M,N3M,NDIM)) (&(BASEMETMACP0A1(idxvol,N1BND,N2BND,N3BND,0)));
598 #endif
599 
600 
601  // initialize global pointers for multi-dimensional arrays for radiation
602  // KORAL
603  set_arrays_multidimen_rad();
604 
605 
606 }
Generated on Fri May 20 2016 15:52:35 for HARM by   doxygen 1.8.3.1