HARM
harm and utilities
 All Data Structures Files Functions Variables Typedefs Macros Pages
step_ch.c
Go to the documentation of this file.
1 
7 #include "decs.h"
8 
9 
10 static void setup_rktimestep(int truestep, int *numtimeorders,
11  FTYPE (*p)[NSTORE2][NSTORE3][NPR], FTYPE (*pstag)[NSTORE2][NSTORE3][NPR], FTYPE (*ucons)[NSTORE2][NSTORE3][NPR], FTYPE (*vpot)[NSTORE1+SHIFTSTORE1][NSTORE2+SHIFTSTORE2][NSTORE3+SHIFTSTORE3], FTYPE (*Bhat)[NSTORE2][NSTORE3][NPR],
12  FTYPE (*pk)[NSTORE1][NSTORE2][NSTORE3][NPR],
13  FTYPE (*pii[MAXTIMEORDER])[NSTORE2][NSTORE3][NPR],FTYPE (*pbb[MAXTIMEORDER])[NSTORE2][NSTORE3][NPR],FTYPE (*pff[MAXTIMEORDER])[NSTORE2][NSTORE3][NPR],
14  FTYPE (*uii[MAXTIMEORDER])[NSTORE2][NSTORE3][NPR],FTYPE (*uff[MAXTIMEORDER])[NSTORE2][NSTORE3][NPR],FTYPE (*ucum[MAXTIMEORDER])[NSTORE2][NSTORE3][NPR],
15  FTYPE (*CUf)[NUMDTCUFS],FTYPE (*CUnew)[NUMDTCUFS]);
16 
17 
18 
19 
20 static int pre_stepch(int *dumpingnext, FTYPE (*prim)[NSTORE2][NSTORE3][NPR]);
21 static int post_stepch(int *dumpingnext, FTYPE fullndt, FTYPE (*prim)[NSTORE2][NSTORE3][NPR], FTYPE (*ucons)[NSTORE2][NSTORE3][NPR]);
22 static int step_ch(int truestep, int *dumpingnext, FTYPE *fullndt,FTYPE (*prim)[NSTORE2][NSTORE3][NPR], FTYPE (*pstag)[NSTORE2][NSTORE3][NPR], FTYPE (*ucons)[NSTORE2][NSTORE3][NPR], FTYPE (*vpot)[NSTORE1+SHIFTSTORE1][NSTORE2+SHIFTSTORE2][NSTORE3+SHIFTSTORE3], FTYPE (*Bhat)[NSTORE2][NSTORE3][NPR], FTYPE (*pl_ct)[NSTORE1][NSTORE2][NSTORE3][NPR2INTERP], FTYPE (*pr_ct)[NSTORE1][NSTORE2][NSTORE3][NPR2INTERP],FTYPE (*F1)[NSTORE2][NSTORE3][NPR+NSPECIAL],FTYPE (*F2)[NSTORE2][NSTORE3][NPR+NSPECIAL],FTYPE (*F3)[NSTORE2][NSTORE3][NPR+NSPECIAL],FTYPE (*Atemp)[NSTORE1+SHIFTSTORE1][NSTORE2+SHIFTSTORE2][NSTORE3+SHIFTSTORE3],FTYPE (*uconstemp)[NSTORE2][NSTORE3][NPR]);
23 static int post_advance(int truestep, int *dumpingnext, int timeorder, int numtimeorders, int finalstep, SFTYPE boundtime, SFTYPE fluxtime, FTYPE (*pi)[NSTORE2][NSTORE3][NPR],FTYPE (*pb)[NSTORE2][NSTORE3][NPR],FTYPE (*pf)[NSTORE2][NSTORE3][NPR],FTYPE (*pstag)[NSTORE2][NSTORE3][NPR],FTYPE (*ucons)[NSTORE2][NSTORE3][NPR], FTYPE (*vpot)[NSTORE1+SHIFTSTORE1][NSTORE2+SHIFTSTORE2][NSTORE3+SHIFTSTORE3],FTYPE (*Bhat)[NSTORE2][NSTORE3][NPR], FTYPE (*F1)[NSTORE2][NSTORE3][NPR+NSPECIAL], FTYPE (*F2)[NSTORE2][NSTORE3][NPR+NSPECIAL], FTYPE (*F3)[NSTORE2][NSTORE3][NPR+NSPECIAL], FTYPE (*Atemp)[NSTORE1+SHIFTSTORE1][NSTORE2+SHIFTSTORE2][NSTORE3+SHIFTSTORE3], FTYPE (*uconstemp)[NSTORE2][NSTORE3][NPR]);
24 
25 
26 
27 
28 
29 
30 
31 
32 
34 int step_ch_full(int truestep, FTYPE (*prim)[NSTORE2][NSTORE3][NPR], FTYPE (*pstag)[NSTORE2][NSTORE3][NPR], FTYPE (*ucons)[NSTORE2][NSTORE3][NPR], FTYPE (*vpot)[NSTORE1+SHIFTSTORE1][NSTORE2+SHIFTSTORE2][NSTORE3+SHIFTSTORE3], FTYPE (*Bhat)[NSTORE2][NSTORE3][NPR], FTYPE (*pl_ct)[NSTORE1][NSTORE2][NSTORE3][NPR2INTERP], FTYPE (*pr_ct)[NSTORE1][NSTORE2][NSTORE3][NPR2INTERP],FTYPE (*F1)[NSTORE2][NSTORE3][NPR+NSPECIAL],FTYPE (*F2)[NSTORE2][NSTORE3][NPR+NSPECIAL],FTYPE (*F3)[NSTORE2][NSTORE3][NPR+NSPECIAL],FTYPE (*Atemp)[NSTORE1+SHIFTSTORE1][NSTORE2+SHIFTSTORE2][NSTORE3+SHIFTSTORE3],FTYPE (*uconstemp)[NSTORE2][NSTORE3][NPR])
35 {
36  FTYPE fullndt;
37  // dumpingnext[0] = dumping just after this step
38  // dumpingnext[1] = dumping just after the step after this step
39  int dumpingnext[2];
40 
41 
42 
43  // things to do before taking full step
44  if(truestep) pre_stepch(dumpingnext,prim);
45 
46  // take full step
47  step_ch(truestep,dumpingnext, &fullndt,prim,pstag,ucons,vpot,Bhat,pl_ct, pr_ct, F1, F2, F3,Atemp,uconstemp);
48 
49  // things to do after taking full step
50  if(truestep) post_stepch(dumpingnext, fullndt, prim, ucons);
51 
52 
53  if(truestep){ // don't do if just passing through -- otherwise would end up looping endlessly!
54  // add up contributions to flux through horizon (really inner boundary)
55  if(DOEVOLVEMETRIC && (EVOLVEMETRICSUBSTEP==0 || EVOLVEMETRICSUBSTEP==2) ){
56  compute_new_metric_longsteps(prim,pstag,ucons,vpot,Bhat,pl_ct, pr_ct, F1, F2, F3,Atemp,uconstemp);
57  }
58 
59  // below must come after longstep update so that don't change before metric step taken!
60  control_metric_update();
61 
62  postdt(); // here one can alter variables and try to restart, or implement any post step operations
63 
64 
65  }
66 
67 
68  // must check before MPI operation (since asymmetries would desynchronize cpus)
69  MYFUN(error_check(2),"main.c","error_check",1);
70  //^^ otherwise ok
71  // eventually all cpus come here, either in failure mode or not, and cleanly tell others if failed/exit/dump/etc.
72 
73 
74 
75 
76 
77  return(0);
78 }
79 
80 
81 
83 int step_ch(int truestep, int *dumpingnext, FTYPE *fullndt,FTYPE (*prim)[NSTORE2][NSTORE3][NPR], FTYPE (*pstag)[NSTORE2][NSTORE3][NPR], FTYPE (*ucons)[NSTORE2][NSTORE3][NPR], FTYPE (*vpot)[NSTORE1+SHIFTSTORE1][NSTORE2+SHIFTSTORE2][NSTORE3+SHIFTSTORE3], FTYPE (*Bhat)[NSTORE2][NSTORE3][NPR], FTYPE (*pl_ct)[NSTORE1][NSTORE2][NSTORE3][NPR2INTERP], FTYPE (*pr_ct)[NSTORE1][NSTORE2][NSTORE3][NPR2INTERP],FTYPE (*F1)[NSTORE2][NSTORE3][NPR+NSPECIAL],FTYPE (*F2)[NSTORE2][NSTORE3][NPR+NSPECIAL],FTYPE (*F3)[NSTORE2][NSTORE3][NPR+NSPECIAL],FTYPE (*Atemp)[NSTORE1+SHIFTSTORE1][NSTORE2+SHIFTSTORE2][NSTORE3+SHIFTSTORE3],FTYPE (*uconstemp)[NSTORE2][NSTORE3][NPR])
84 {
85  int step_ch_simplempi(int truestep, int *dumpingnext, FTYPE *fullndt,FTYPE (*prim)[NSTORE2][NSTORE3][NPR], FTYPE (*pstag)[NSTORE2][NSTORE3][NPR], FTYPE (*ucons)[NSTORE2][NSTORE3][NPR], FTYPE (*vpot)[NSTORE1+SHIFTSTORE1][NSTORE2+SHIFTSTORE2][NSTORE3+SHIFTSTORE3], FTYPE (*Bhat)[NSTORE2][NSTORE3][NPR], FTYPE (*pl_ct)[NSTORE1][NSTORE2][NSTORE3][NPR2INTERP], FTYPE (*pr_ct)[NSTORE1][NSTORE2][NSTORE3][NPR2INTERP],FTYPE (*F1)[NSTORE2][NSTORE3][NPR+NSPECIAL],FTYPE (*F2)[NSTORE2][NSTORE3][NPR+NSPECIAL],FTYPE (*F3)[NSTORE2][NSTORE3][NPR+NSPECIAL],FTYPE (*Atemp)[NSTORE1+SHIFTSTORE1][NSTORE2+SHIFTSTORE2][NSTORE3+SHIFTSTORE3],FTYPE (*uconstemp)[NSTORE2][NSTORE3][NPR]);
86  int step_ch_supermpi(int truestep, int *dumpingnext, FTYPE *fullndt,FTYPE (*prim)[NSTORE2][NSTORE3][NPR], FTYPE (*pstag)[NSTORE2][NSTORE3][NPR], FTYPE (*ucons)[NSTORE2][NSTORE3][NPR], FTYPE (*vpot)[NSTORE1+SHIFTSTORE1][NSTORE2+SHIFTSTORE2][NSTORE3+SHIFTSTORE3], FTYPE (*Bhat)[NSTORE2][NSTORE3][NPR], FTYPE (*pl_ct)[NSTORE1][NSTORE2][NSTORE3][NPR2INTERP], FTYPE (*pr_ct)[NSTORE1][NSTORE2][NSTORE3][NPR2INTERP],FTYPE (*F1)[NSTORE2][NSTORE3][NPR+NSPECIAL],FTYPE (*F2)[NSTORE2][NSTORE3][NPR+NSPECIAL],FTYPE (*F3)[NSTORE2][NSTORE3][NPR+NSPECIAL],FTYPE (*Atemp)[NSTORE1+SHIFTSTORE1][NSTORE2+SHIFTSTORE2][NSTORE3+SHIFTSTORE3],FTYPE (*uconstemp)[NSTORE2][NSTORE3][NPR]);
87 
88 
89 
90 #if(SIMULBCCALC==-1)
91  MYFUN(step_ch_simplempi(truestep, dumpingnext, fullndt,prim,pstag,ucons,vpot,Bhat,pl_ct, pr_ct,F1,F2,F3,Atemp,uconstemp),"step_ch.c:step_ch()", "step_ch_simplempi()", 1);
92 #else
93  MYFUN(step_ch_supermpi(truestep, dumpingnext, fullndt,prim,pstag,ucons,vpot,Bhat,pl_ct, pr_ct,F1,F2,F3,Atemp,uconstemp),"step_ch.c:step_ch()", "step_ch_supermpi()", 1);
94 #endif
95 
96 
97  /* done! */
98  return (0);
99 }
100 
101 
102 
105 int pre_stepch(int *dumpingnext, FTYPE (*prim)[NSTORE2][NSTORE3][NPR])
106 {
107 
108 
109 
110 #if(PRODUCTION==0)
111  trifprintf( "\n#step=%ld",realnstep);
112 #endif
113 
114 
115 
116  // if not doing diagnostics, then dumpingnext will be 0
117  // first check if dumping next step
118  // note dt is already set to correct value here
119  dumpingnext[0]=(diag(FUTURE_OUT,t+dt,nstep+1,realnstep+1)==DOINGFUTUREOUT);
120  dumpingnext[1]=(diag(FUTURE_OUT,t+dt+dt/2,nstep+2,realnstep+2)==DOINGFUTUREOUT); // estimate of t+dt+dt/2 good with 1<SAFEFACTOR<2
121 
122 
123 
124 
125  if(dumpingnext[0] || dumpingnext[1]){
126 #if(CALCFARADAYANDCURRENTS)
127  if((WHICHCURRENTCALC==CURRENTCALC0)||(WHICHCURRENTCALC==CURRENTCALC2)){
128  // for this method, all J components are located back in time
129  current_doprecalc(CURTYPEX,prim);
130  current_doprecalc(CURTYPEY,prim);
131  current_doprecalc(CURTYPEZ,prim);
132  // get faraday in case used for time averaging or dumping it
133  current_doprecalc(CURTYPEFARADAY,prim);
134  }
135 #endif
136  }
137 
138 
139  return(0);
140 
141 
142 }
143 
144 
145 
146 
149 int post_stepch(int *dumpingnext, FTYPE fullndt, FTYPE (*prim)[NSTORE2][NSTORE3][NPR], FTYPE (*ucons)[NSTORE2][NSTORE3][NPR])
150 {
151 
152 
153 
154 
155  if(dumpingnext[0]){
156 #if(CALCFARADAYANDCURRENTS)
157  // compute final faradays
158  if(WHICHCURRENTCALC==CURRENTCALC1){
159  // compute faraday components needed for time centering of J
160  current_doprecalc(CURTYPET,prim);
161  // J is located at this time
162  current_doprecalc(CURTYPEX,prim);
163  current_doprecalc(CURTYPEY,prim);
164  current_doprecalc(CURTYPEZ,prim);
165  current_doprecalc(CURTYPEFARADAY,prim); // for diagnostics
166  }
167  // compute current
168  current_calc(GLOBALPOINT(cfaraday));
169 #endif
170  }
171 
172 
173 #if(ACCURATEDIAG==0)
174  // compute flux diagnostics
175  // this doesn't exactly make conservation work -- should have in middle step point using full step. But if PARA, no middle point that's exact.
176  // think about where to put this
177  // GODMARK : use of globals
178  diag_flux(prim,F1, F2, F3, dt); // should use REAL dt, not within a timeorderd RK integration step
179  // horizon_flux(F1,dt); // subsumed
180 #endif
181 
182  // general flux only done on full steps since no requirement for accuracy and code can be expensive computationally
183  diag_flux_general(prim,dt);// Should be full dt, not substep dt.
184 
185 
186 
187 
188  if(DOONESTEPDUACCOUNTING){
189  // do one-step accounting for fixup related things (failures, floors, fixups, etc.)
190  // ucons is in UEVOLVE form (originates from unewglobal in step_ch_full() called in main.c)
191  diag_fixup_allzones(prim, ucons);
192  }
193 
194 
195 
196 #if(PRODUCTION==0)
197  trifprintf( "[d]");
198 #endif
199 
200 
201 
203  //
204  // increment time
205  //
207  t += dt;
208  tstepparti = tsteppartf = t;
209 
210  realnstep++;
211 
213  //
214  // set next timestep
215  //
216  //
217  // Notes:
218  //
219  // dt is computed in a few different places for different steps of the calculation:
220  // 1) flux.c (per dimension over all cells):
221  // A) dtij = cour * dx[dir] / ctop;
222  // B) if (dtij < *ndt){ *ndt = dtij; }
223  //
224  // 2) advance.c: prepare_globaldt() ("min" over each dimension and operator):
225  // A) wavedt = 1. / (1. / ndt1 + 1. / ndt2 + 1. / ndt3);
226  // B) *ndt = defcon * MIN(wavedt , accdt);
227  // C) etc...
228  //
229  // 3) step_ch.c: ("min" over each substep of RK steps):
230  // A) step_ch_simplempi(): if(ndt<lastndt) lastndt=ndt;
231  // B) HERE: to find global minimum over all CPUs
232  // C) HERE: to use safety factor
233  // D) HERE: constrain new dt if near end of calculation
234  //
235  // But, should do #1's MIN over all cells per CPU and #3's global MIN over all CPUs together for perfect MPI consistency.
236  // Further, should compute "wavedt" per cell since otherwise non-local dt_i can impact minimum dt. Can lead to dt smaller by as much as 3X (with no correct numerical method effect) !
237  //
239  // find global minimum value of ndt over all cpus
240  mpifmin(&fullndt);
241  // note that this assignment to fullndt doesn't get returned outside function
242  if (fullndt > SAFE * dt) fullndt = SAFE * dt;
243  dt = fullndt;
244 
246  //
247  // don't step beyond end of run
248  //
250  if(onemorestep){
251  dualfprintf(fail_file,"Got onemorestep dt=%g\n",dt);
252  // check if previous step was onemorestep==1
253  reallaststep=1;
254  dt=SMALL;
255  }
256  else{
257  if (t + dt > tf){
258  dualfprintf(fail_file,"Got t+dt>tf : %g %g %g\n",t,dt,tf);
259  dt = tf - t;
260  onemorestep=1;
261  }
262  else if (t + dt == tf){
263  dualfprintf(fail_file,"Got t+dt==tf : %g %g %g\n",t,dt,tf);
264  reallaststep=1;
265  }
266  // make sure don't get pathological case of dt=0 on last step
267  if(dt<SMALL){
268  dualfprintf(fail_file,"Got dt<SMALL : %g\n",dt);
269  reallaststep=1;
270  dt=SMALL;
271  }
272  }
273 
274 
275 
276 
277  if(dt==0.0 && t>=tf){
278  // somehow got here, then finish
279  dualfprintf(fail_file,"SOMEHOW GOT HERE\n");
280  reallaststep=1;
281  }
282 
283 
284 
285  return(0);
286 
287 
288 
289 }
290 
291 
292 
293 
294 
296 int pre_advance(int timeorder, int numtimeorders, int finalstep, FTYPE (*pi)[NSTORE2][NSTORE3][NPR],FTYPE (*pb)[NSTORE2][NSTORE3][NPR],FTYPE (*pf)[NSTORE2][NSTORE3][NPR])
297 {
298 
299 
301  //
302  // prefixup
303  //
305 #if(EOMTYPE!=EOMFFDE && EOMTYPE!=EOMFFDE2)
306  // force-free and cold GRMHD don't use pre_fixup, currently, even if could
307  MYFUN(pre_fixup(STAGEM1, pi),"step_ch.c:step_ch_simple()", "pre_fixup()", 1);
308 #endif
309 
310 
311  return(0);
312 
313 }
314 
315 
316 
317 
318 
319 
321 int post_advance(int truestep, int *dumpingnext, int timeorder, int numtimeorders, int finalstep, SFTYPE boundtime, SFTYPE fluxtime, FTYPE (*pi)[NSTORE2][NSTORE3][NPR],FTYPE (*pb)[NSTORE2][NSTORE3][NPR],FTYPE (*pf)[NSTORE2][NSTORE3][NPR],FTYPE (*pstag)[NSTORE2][NSTORE3][NPR],FTYPE (*ucons)[NSTORE2][NSTORE3][NPR], FTYPE (*vpot)[NSTORE1+SHIFTSTORE1][NSTORE2+SHIFTSTORE2][NSTORE3+SHIFTSTORE3],FTYPE (*Bhat)[NSTORE2][NSTORE3][NPR], FTYPE (*F1)[NSTORE2][NSTORE3][NPR+NSPECIAL], FTYPE (*F2)[NSTORE2][NSTORE3][NPR+NSPECIAL], FTYPE (*F3)[NSTORE2][NSTORE3][NPR+NSPECIAL], FTYPE (*Atemp)[NSTORE1+SHIFTSTORE1][NSTORE2+SHIFTSTORE2][NSTORE3+SHIFTSTORE3], FTYPE (*uconstemp)[NSTORE2][NSTORE3][NPR])
322 {
323 
324 
325 
326 
327 
329  //
330  // Use tracked A_i to update magnetic field
331  // Required to keep A_i in synch with field and only is different at the machine error level (which grows over long times, so why this is required and hardless)
332  //
334  if(EVOLVEWITHVPOT && TRACKVPOT && (finalstep)){
335  // if evolving with vpot, then assume good enough to use final full timestep A_i to obtain new point fields
336  // less expensive than doing every substep and only wrong at machine error with factors extra for number of substeps
337  // SUPERGODMARK: Check convergence rate and check errors!! SUPERCHANGINGMARK
338 
339  // only do this on the final step where A_i has been set as the cumulative Acum (like ucum) so that not just an arbitrary intermediate step that redefines B's.
340 
341  // NOTEMARK: No, actually want to do on every substep that defines vpot (vpot holds substep (Uf-like) value on substeps except final step where it holds cumulative final value (Ucum-like)) since need staggered field (say inside NS for init.nsbh.c) to be consistent with A_i so that the smooth extrapolation of A_i translates into smooth Bstag^i so that interpolation of Bstag->Bcent gives smooth results.
342  // NOTEMARK: If A_i is set inconsistent with EMF_i, then this step will reveal A_i overridding EMF_i results and probably producing undesirable results.
343  // NOTEMARK: No again, if vpot update is consistent with EMFs (as it should be as code is setup in fluxvpot.c), then Bstag only off by machine error from A_i, so only need to do this infrequently.
344  evolve_withvpot(boundtime, pf, pstag, ucons, vpot, Bhat, F1, F2, F3, Atemp,uconstemp); // boundtime since this time is used for a boundary call inside this function
345  }
346 
347 
348 
349 
351  //
352  // need to compute various things for EOS
353  // use primitive located where EOS will be used (pb)
354  //
356  compute_EOS_parms(WHICHEOS,GLOBALPOINT(EOSextraglobal),pb);
357 
358 
360  //
361  // post-advance bound and fixup
362  //
364 
365 #if( (1 == CHECKSOLUTION || UTOPRIMADJUST == UTOPRIMAVG) )
366  // if CHECKSOLUTION==1, then need values to be bounded right now, since use them to check whether even good solutions are really good.
367  // post_fixup() will use previous time step pff boundary values to fixup_utoprim() if this is not called.
368  // bound advanced values before post_fixup() so fixup_utoprim() has updated boundary values to base fixup on.
369 
370 
371  if(failed>0) dualfprintf(fail_file,"1failed=%d\n",failed);
372 
374  //
375  // bounding 1 layer so fixups have that layer available for fixing up inversion failures
376  // not done if don't care about perfect MPI consistency with failures -- then fixups don't use values across MPI boundaries. Will be less robust in general, but still consistent behavior across MPI boundaries.
377  //
379 
380 #if(PRODUCTION==0)
381  trifprintf("[b1");
382 #endif
383 
384  MYFUN(bound_beforeevolveprim(STAGEM1,finalstep,boundtime,pf,pstag,ucons,UTOPRIMFIXMPICONSISTENT>0),"step_ch.c:post_advance()", "bound_beforeevolveprim()", 1);
385 
386 #if(PRODUCTION==0)
387  trifprintf("]");
388 #endif
389 
390 
391 
392 
393  if(failed>0) dualfprintf(fail_file,"2failed=%d\n",failed);
394 
395 
396 #if(PRODUCTION==0)
397  trifprintf("[x");
398 #endif
399 
400 
402  // done when all timeorders are completed, so stencil used doesn't matter
403  //
404  // postfixup
405  // post_fixup: If this modifies values and then uses the modified values for other modified values, then must bound_prim() after this
406  // if one doesn't care about MPI being same as non-MPI, then can move bound_prim() above to below error_check() and then remove prc<-pv in fixup_utoprim()
407  //
409  MYFUN(post_fixup(STAGEM1,finalstep,boundtime, pf,pb,ucons),"step_ch.c:advance()", "post_fixup()", 1);
410 
411 #if(PRODUCTION==0)
412  trifprintf( "]");
413 #endif
414 
415 
416 
417 
418 #else
419 
421  //
422  // just report problems, don't fix them
423  // Then no need for boundary calls
424  //
426  MYFUN(post_fixup_nofixup(STAGEM1,finalstep,boundtime, pf,pb,ucons),"step_ch.c:advance()", "post_fixup_nofixup()", 1);
427 
428 #endif
429 
430 
432  //
433  // Synch CPUs
434  //
436 
437  // must check before MPI operation (since asymmetries would
438  // desynchronize cpus
439  MYFUN(error_check(1),"step_ch.c", "error_check", 1);
440 
441  // bound final values (comes after post_fixup() since changes made by post_fixup)
442  //#if(MPIEQUALNONMPI)
443 
444 
445 
446 #if(ASYMDIAGCHECK)
447  dualfprintf(fail_file,"1before bound\n");
448  asym_compute_2(pf);
449  dualfprintf(fail_file,"2before bound\n");
450 #endif
451 
452 
453 
454 
455 #if(PRODUCTION==0)
456  trifprintf("[rf");
457 #endif
458 
459 
460 
461  if( DOGRIDSECTIONING){
462  // this must come before last bound() call so boundary conditions set consistently with new section so next step has all values needed in ghost cells
463  // redo all such enerregions since may be time-dependent
464  recompute_fluxpositions(0,timeorder, numtimeorders, nstep,t);
465  }
466 
467 #if(PRODUCTION==0)
468  trifprintf("]");
469 #endif
470 
471 
472 
473  if(timeorder==numtimeorders-1){
474  // ensure total e-m conservation for final solution after all sub-steps
475  int finaluu=1;
476  consfixup_allzones(finaluu,pf,ucons);
477  }
478 
479 
480 
481 
482 #if(PRODUCTION==0)
483  trifprintf( "[b2");
484 #endif
485 
486 
488  //
489  // normal bondary call (only bounds p (centered) and pstag (staggered field) and not A_i or F or other things)
490  // required in general
491  //
493  MYFUN(bound_evolveprim(STAGEM1,finalstep,boundtime,pf,pstag,ucons,USEMPI),"step_ch.c:post_advance()", "bound_evolveprim()", 2);
494 
495  //#endif
496 
497 #if(PRODUCTION==0)
498  trifprintf( "]");
499 #endif
500 
501 
502 
503 
504 
505 
507  //
508  // Compute divb for point-field method
509  //
511  if(truestep){
512  // update Bhat so later can compute divb
513  if(FLUXB==FLUXCTSTAG && extrazones4emf==0 && dofluxreconevolvepointfield==1 && DOENOFLUX == ENOFLUXRECON){
514  //bound_uavg(STAGEM1,unew); // DEBUG
515  // here utoinvert is point conserved field at staggered FACE position
516  field_Bhat_fluxrecon(pf,ucons,Bhat);
517  }
518  }
519 
520 
522  //
523  // Compute current update
524  //
526  if(truestep){ // don't do if just passing through
527 
528  if(dumpingnext[0] || dumpingnext[1]){
529 #if(CALCFARADAYANDCURRENTS)
530 
531 #if(PRODUCTION==0)
532  trifprintf( "[cu");
533 #endif
534 
535  if((WHICHCURRENTCALC==CURRENTCALC0)||(WHICHCURRENTCALC==CURRENTCALC2)){
536  // puts J at the time center, but hard to know if RK is at mid point in time except for midpoint method
537  // compute current_doprecalc if near half-step in time
538  if(
539  ((numtimeorders>=3)&&(timeorder==1))
540  ||((numtimeorders<=2)&&(timeorder==0))
541  )
542  current_doprecalc(CURTYPET,pf); // should be called using half-time step data
543  }
544 
545 #if(PRODUCTION==0)
546  trifprintf( "]");
547 #endif
548 
549 #endif
550  }// end if dumping next
551 
552 
553 
554 
555  }// end if truestep
556 
557 
558 
559 
560 
561 
563  //
564  // Diagnostics
565  //
567  if(truestep){ // don't do if just passing through
568  /* perform diagnostics */
569  // no error check since assume if step_ch passed, diag(1) will pass
570  if (DODIAGS && DODIAGEVERYSUBSTEP ){ //SASMARK -- moved the diags calls here
571  GLOBALPOINT(pdump) = pf;
572  diag(DUMP_OUT,boundtime,nstep,realnstep); // boundtime corresponds to pf for this timeorder
573 #if(PRODUCTION==0)
574  trifprintf( "D");
575 #endif
576  }
577  }
578 
579  return(0);
580 }
581 
582 
583 
584 
585 
586 
589 int step_ch_simplempi(int truestep, int *dumpingnext, FTYPE *fullndt, FTYPE (*prim)[NSTORE2][NSTORE3][NPR], FTYPE (*pstag)[NSTORE2][NSTORE3][NPR], FTYPE (*ucons)[NSTORE2][NSTORE3][NPR], FTYPE (*vpot)[NSTORE1+SHIFTSTORE1][NSTORE2+SHIFTSTORE2][NSTORE3+SHIFTSTORE3], FTYPE (*Bhat)[NSTORE2][NSTORE3][NPR], FTYPE (*pl_ct)[NSTORE1][NSTORE2][NSTORE3][NPR2INTERP], FTYPE (*pr_ct)[NSTORE1][NSTORE2][NSTORE3][NPR2INTERP],FTYPE (*F1)[NSTORE2][NSTORE3][NPR+NSPECIAL],FTYPE (*F2)[NSTORE2][NSTORE3][NPR+NSPECIAL],FTYPE (*F3)[NSTORE2][NSTORE3][NPR+NSPECIAL],FTYPE (*Atemp)[NSTORE1+SHIFTSTORE1][NSTORE2+SHIFTSTORE2][NSTORE3+SHIFTSTORE3],FTYPE (*uconstemp)[NSTORE2][NSTORE3][NPR])
590 {
591  int advance(int truestep, int stage, FTYPE (*pi)[NSTORE2][NSTORE3][NPR],FTYPE (*pb)[NSTORE2][NSTORE3][NPR], FTYPE (*pf)[NSTORE2][NSTORE3][NPR],
592  FTYPE (*pstag)[NSTORE2][NSTORE3][NPR],
593  FTYPE (*pl_ct)[NSTORE1][NSTORE2][NSTORE3][NPR2INTERP], FTYPE (*pr_ct)[NSTORE1][NSTORE2][NSTORE3][NPR2INTERP],
594  FTYPE (*F1)[NSTORE2][NSTORE3][NPR+NSPECIAL],FTYPE (*F2)[NSTORE2][NSTORE3][NPR+NSPECIAL],FTYPE (*F3)[NSTORE2][NSTORE3][NPR+NSPECIAL],
595  FTYPE (*vpot)[NSTORE1+SHIFTSTORE1][NSTORE2+SHIFTSTORE2][NSTORE3+SHIFTSTORE3],
596  FTYPE (*ui)[NSTORE2][NSTORE3][NPR],FTYPE (*uf)[NSTORE2][NSTORE3][NPR], FTYPE (*ucum)[NSTORE2][NSTORE3][NPR],
597  FTYPE *CUf,FTYPE *CUnew,SFTYPE fluxdt, SFTYPE boundtime, SFTYPE fluxtime, int timeorder, int numtimeorders, FTYPE *ndt);
598 
599  int pre_advance(int timeorder, int numtimeorders, int finalstep, FTYPE (*pi)[NSTORE2][NSTORE3][NPR],FTYPE (*pb)[NSTORE2][NSTORE3][NPR],FTYPE (*pf)[NSTORE2][NSTORE3][NPR]);
600  int asym_compute_2(FTYPE (*prim)[NSTORE2][NSTORE3][NPR]);
601  int pl,pliter;
602  //
603  int boundstage;
604  SFTYPE mydt;
605  int stage, stagei,stagef;
606  FTYPE ndt,lastndt;
607  FTYPE (*pi)[NSTORE2][NSTORE3][NPR];
608  FTYPE (*pb)[NSTORE2][NSTORE3][NPR];
609  FTYPE (*pf)[NSTORE2][NSTORE3][NPR];
610  FTYPE (*prevpf)[NSTORE2][NSTORE3][NPR];
611  FTYPE (*pii[MAXTIMEORDER])[NSTORE2][NSTORE3][NPR];
612  FTYPE (*pbb[MAXTIMEORDER])[NSTORE2][NSTORE3][NPR];
613  FTYPE (*pff[MAXTIMEORDER])[NSTORE2][NSTORE3][NPR];
614  FTYPE (*uii[MAXTIMEORDER])[NSTORE2][NSTORE3][NPR];
615  FTYPE (*uff[MAXTIMEORDER])[NSTORE2][NSTORE3][NPR];
616  FTYPE (*ucum[MAXTIMEORDER])[NSTORE2][NSTORE3][NPR];
617  // FTYPE alphaik[MAXSTAGES][MAXSTAGES],betaik[MAXSTAGES];
618  FTYPE CUf[MAXTIMEORDER][NUMDTCUFS],CUnew[MAXTIMEORDER][NUMDTCUFS];
619 
620 
621  int i, j, k;
622  int numtimeorders;
623  int timeorder;
624  int finalstep,initialstep;
625  SFTYPE fluxdt[MAXTIMEORDER], boundtime[MAXTIMEORDER], fluxtime[MAXTIMEORDER];
626 
627 
628 
629 
630 
632  //
633  // setup time-stepping
634  //
636  setup_rktimestep(truestep, &numtimeorders,prim,pstag,ucons,vpot,Bhat,GLOBALPOINT(pk),pii,pbb,pff,uii,uff,ucum,CUf,CUnew);
637 
638 
640  //
641  // Obtain initial time of substep, final time of substep, and true dt used for flux conservation that is used to iterate ucum in advance.c
642  //
644  get_truetime_fluxdt(numtimeorders, dt, CUf, CUnew, fluxdt, boundtime, fluxtime, NULL,NULL);
645 
646 
647  // global debug tracking var
648  numstepparts=numtimeorders;
649  // Initialize old and new dt to be very large so minimization procedure always obtains new dt
650  ndt=lastndt=BIG;
651 
652 
653 
654 
656  //
657  // general temporal loop
658  //
660  for(timeorder=0;timeorder<numtimeorders;timeorder++){
661 
662 
664  //
665  //the starting and the ending times of the current substep
666  //
668  if(timeorder!=numtimeorders-1){
669  tstepparti = t + CUf[timeorder][3] * dt;
670  tsteppartf = t + CUf[timeorder][2] * dt + CUf[timeorder][3] * dt;
671  }
672  else{
673  tstepparti=t;
674  tsteppartf=t+dt;
675  }
676 
677 
679  //
680  // global debug tracking var
681  //
683  steppart=timeorder;
684  if(timeorder==numtimeorders-1) finalstep=1; else finalstep=0;
685  if(timeorder==0) initialstep=1; else initialstep=0;
686 
687 #if(PRODUCTION==0)
688  trifprintf("|%ds",timeorder);
689 #endif
690 
691 
692 
694  //
695  // BEFORE ADVANCE
696  //
698  pre_advance(timeorder, numtimeorders, finalstep, pii[timeorder],pbb[timeorder],pff[timeorder]);
699 
700 
701 
703  //
704  // ADVANCE
705  //
707  if(SPLITNPR){
708 
709 
711  //
712  // SPLITNPR ADVANCE
713  //
715 
716 
717 #if(SPLITNPR&&(USESTOREDSPEEDSFORFLUX==0 || STOREWAVESPEEDS==0))
718 #error Not correct use of SPLITNPR and wavespeeds
719 #endif
720  // Note that even though PLOOP is split-up, some things are still repeated for these 2 advance() calls:
721  // 1) vchar for Riemann solver still needed -- could store vchar on first pass and use on second pass, but then using old B for vchar() on second pass and want to use new B everywhere for second pass
722  // 2) dt is minimized twice -- seems second pass will overwrite first pass
723  // 3) getting of geometry and state everywhere is duplicated
724  // NOTE OF IMPORTANCE: get_state needs all primitives, so must feed non-PLOOP quantities
725  // 4) get_wavespeeds requires ALL quantities during interpolation -- must set USESTOREDSPEEDSFORFLUX==1 and STOREWAVESPEEDS>0 so computes wave speeds across 2 advance() calls and uses centered speed for interface to avoid needing during interpolation for p_l,p_r->wavespeed@interface
726  // Basically flux_compute_general(p_l,p_r needs velocities too) for first pass
727  // 5) p2SFUevolve's get_state needs velocity -- that is, flux_B(B,v) so have to still interpolate v with B -- save v-interpolation (all interpolations) and avoid interpolations on second pass.
728  // Made pleft and pright [3]X in size to store across advance(), dq is already this way
729  // 6) checked NPR and NPR2INTERP uses. All such explicit uses are ok
730  // e.g. grep -e "[^\[]NPR" *.c
731 
732  // Things NOT done on first pass:
733  // 1) source term if NOGDETBi=0
734  // 2) metric update if evolving metric
735  // 3) all non-B quantities if doing PLOOP or PLOOPINTERP
736  // 4) diag_flux, diag_source, diag_source_all -- all use PDIAGLOOP
737  // 5) fixup1zone() or fixup()
738 
739  // Things to fix:
740  // 3) compute_df_line_new during interpolation needs get_state for u and b -- only for pressure indicator for WENO
741  // 4) calculation of bsq requires b^\mu requires B and v/u/urel -- called PLOOP or what?
742  // 5) Check all PLOOP's everywhere to ensure doing ALL NPR if required
743  // 6) could bound only field when bounding between
744 
745 
746 
747  // choice for range of PLOOP
748  // choose to only update magnetic field
749  nprstart=0;
750  nprend=2;
751  nprlist[0]=B1;
752  nprlist[1]=B2;
753  nprlist[2]=B3;
754  // choice for range of PLOOPINTERP
755  // but interpolate everything on first pass (GODMARK: in reality only need to do field on timeorder=0 since last field update is next field pbb??? -- probably only true if doing simple RK methods (1,2))
756 #pragma omp parallel private(pl)
757  {
758  npr2interpstart=0;
759  npr2interpend=NPR2INTERP-1;
760  for(pl=npr2interpstart;pl<=npr2interpend;pl++) npr2interplist[pl]=pl;
761  // choice for range of PLOOPNOTINTERP
762  npr2notinterpstart=0;
763  npr2notinterpend=-1;
764  npr2notinterplist[0]=0;
765  }
766  advancepassnumber=0;
767 
768 
769  // advance (field only)
770  // Only field parts of pff, uff, and ucum updated
771  // on final timeorder, ucum used to get final B that will be used as B for final timeorder for non-field quantities
772  MYFUN(advance(truestep, STAGEM1, pii[timeorder], pbb[timeorder], pff[timeorder], pstag, pl_ct, pr_ct, F1, F2, F3, vpot, uii[timeorder], uff[timeorder], ucum[timeorder], CUf[timeorder], CUnew[timeorder], fluxdt[timeorder], boundtime[timeorder], fluxtime[timeorder], timeorder,numtimeorders,&ndt),"step_ch.c:step_ch_simplempi()", "advance()", 1);
773 
774  // only need to bound field, so control PLOOPMPI
775  nprboundstart=0;
776  nprboundend=2;
777  nprboundlist[0]=B1;
778  nprboundlist[1]=B2;
779  nprboundlist[2]=B3;
780  // with magnetic field, only need to bound
781  MYFUN(bound_evolveprim(STAGEM1,finalstep,boundtime[timeorder],pff[timeorder],pstag,uff[timeorder],USEMPI),"step_ch.c:step_ch_simplempi()", "bound_evolveprim()", 1);
782 
783 
784 
785  // now copy over pff to pbb so Bnew is used in flux calculation of non-field quantities
786  COMPFULLLOOP{
787  PLOOPBONLY(pl) MACP1A1(pbb,timeorder,i,j,k,pl)=MACP1A1(pff,timeorder,i,j,k,pl);
788  }
789 
790 
791  // choice for range of PLOOP
792  // don't update field on second pass (result overwritten in utoprimgen.c even if done)
793  nprstart=0;
794  nprend=NPR-1-3; // no field
795  for(pl=nprstart;pl<=nprend;pl++){
796  if(pl>=B1) nprlist[pl]=pl+3; // skip field
797  else nprlist[pl]=pl;
798  }
799  // choice for range of PLOOPINTERP
800  // Interpolation of only fields on second pass since want to use updated field to compute flux and so at end of both steps first updated field and non-field quantites are fed to inversion for consistent P(Bnew) and Bnew is used
801 #pragma omp parallel private(pl)
802  {
803  npr2interpstart=0;
804  npr2interpend=2;
805  npr2interplist[0]=B1;
806  npr2interplist[1]=B2;
807  npr2interplist[2]=B3;
808  // choice for range of PLOOPNOTINTERP
809  // what not interpolating (all non-field):
810  npr2notinterpstart=0;
811  npr2notinterpend=NPR-1-3; // no field
812  for(pl=npr2notinterpstart;pl<=npr2notinterpend;pl++){
813  if(pl>=B1) npr2notinterplist[pl]=pl+3; // skip field
814  else npr2notinterplist[pl]=pl;
815  }
816  }
817  advancepassnumber=1;
818 
819 
820  // advance (non-field quantities)
821  // only non-field parts of pff, uff, ucum updated
822  MYFUN(advance(truestep, STAGEM1, pii[timeorder], pbb[timeorder], pff[timeorder], pstag, pl_ct, pr_ct, F1, F2, F3, vpot, uii[timeorder], uff[timeorder], ucum[timeorder], CUf[timeorder], CUnew[timeorder], fluxdt[timeorder], boundtime[timeorder], fluxtime[timeorder], timeorder,numtimeorders,&ndt),"step_ch.c:step_ch_simplempi()", "advance()", 1);
823 
824 
825 
827  // go back to defaults
828 
829  // choice for range of PLOOP
830  nprstart=0;
831  nprend=NPR-1;
832  for(pl=nprstart;pl<=nprend;pl++) nprlist[pl]=pl;
833 #pragma omp parallel private(pl)
834  {
835  // choice for range of PLOOPINTERP
836  npr2interpstart=0;
837  npr2interpend=NPR2INTERP-1;
838  for(pl=npr2interpstart;pl<=npr2interpend;pl++) npr2interplist[pl]=pl;
839  // choice for range of PLOOPNOTINTERP
840  npr2notinterpstart=0;
841  npr2notinterpend=-1;
842  npr2notinterplist[0]=0;
843  }
844  advancepassnumber=-1;
845 
846  // default choice for range of PLOOPMPI
847  nprboundstart=0;
848  nprboundend=NPRBOUND-1;
849  for(pl=nprboundstart;pl<=nprboundend;pl++) nprboundlist[pl]=pl;
850 
851 
852 
853 
854  }
855  else{
856 
857 
859  //
860  // NORMAL (not SPLITNPR) ADVANCE
861  //
863 
864  advancepassnumber=-1; // implies do all things, no split
865  // advance (all)
866  MYFUN(advance(truestep, STAGEM1, pii[timeorder], pbb[timeorder], pff[timeorder], pstag, pl_ct, pr_ct, F1, F2, F3, vpot, uii[timeorder], uff[timeorder], ucum[timeorder], CUf[timeorder], CUnew[timeorder], fluxdt[timeorder], boundtime[timeorder], fluxtime[timeorder], timeorder,numtimeorders,&ndt),"step_ch.c:step_ch_simplempi()", "advance()", 1);
867 
868 
869 
870  } // end else if normal advance()
871 
872 
873 
874 
875 
876 
877 
879  //
880  // need to minimize dt over all substeps rather than just last step
881  // lastndt holds cumulative-all-substep ndt
882  //
884  if(ndt<lastndt) lastndt=ndt;
885 
886 
887 
889  //
890  // POST ADVANCE
891  //
893  post_advance(truestep, dumpingnext, timeorder, numtimeorders, finalstep, boundtime[timeorder], fluxtime[timeorder], pii[timeorder],pbb[timeorder],pff[timeorder],pstag,ucons,vpot,Bhat, F1, F2, F3, Atemp, uconstemp);
894 
895 
896 
897  }// end timeorder
898 
899 
900 
902  //
903  // pass back the final minimal dt over all substeps
904  //
906  *fullndt = lastndt;
907 
908 
909 
910  /* done! */
911  return (0);
912 }
913 
914 
915 
917 void get_truetime_fluxdt(int numtimeorders, SFTYPE localdt, FTYPE (*CUf)[NUMDTCUFS], FTYPE (*CUnew)[NUMDTCUFS], SFTYPE *fluxdt, SFTYPE *boundtime, SFTYPE *fluxtime, SFTYPE *tstepparti, SFTYPE *tsteppartf)
918 {
919  int timeorder;
920  SFTYPE ufdt[MAXTIMEORDER],ucumdt[MAXTIMEORDER];
921  SFTYPE oldufdt,olducumdt;
922  FTYPE Ui, dUriemann, dUgeom ;
923 
924 
925  // NOT YET:
927  //
928  //the starting and the ending times of the current substep
929  //
931  // if(timeorder!=numtimeorders-1){
932  // *tstepparti = t + CUf[timeorder][3] * localdt;
933  // *tsteppartf = t + CUf[timeorder][2] * localdt + CUf[timeorder][3] * localdt;
934  // }
935  // else{
936  // *tstepparti=t;
937  // *tsteppartf=t+localdt;
938  // }
939 
940 
941 
942 
943 
944  // initialize
945  for(timeorder=0;timeorder<numtimeorders;timeorder++){
946  fluxdt[timeorder] = 0.0;
947  boundtime[timeorder] = 0.0;
948  fluxtime[timeorder] = 0.0;
949  }
951  //
952  // Get fluxdt (for up to explicit RK5)
953  // See setup_rktimestep() for how Uf^i and unew^i are defined, such that fluxdt[] below just adds all dUe terms appropriately.
954  //
956  for(timeorder=4;timeorder<numtimeorders;timeorder++){ // for >=RK5
957  fluxdt[timeorder-4] += CUnew[timeorder][2]*CUf[timeorder][1]*CUf[timeorder-1][1]*CUf[timeorder-2][1]*CUf[timeorder-3][1]*CUf[timeorder-4][2]*localdt;
958  }
959 
960  for(timeorder=3;timeorder<numtimeorders;timeorder++){ // for >=RK4
961  fluxdt[timeorder-3] += CUnew[timeorder][2]*CUf[timeorder][1]*CUf[timeorder-1][1]*CUf[timeorder-2][1]*CUf[timeorder-3][2]*localdt;
962  }
963 
964  for(timeorder=2;timeorder<numtimeorders;timeorder++){ // for >=RK3
965  fluxdt[timeorder-2] += CUnew[timeorder][2]*CUf[timeorder][1]*CUf[timeorder-1][1]*CUf[timeorder-2][2]*localdt;
966  }
967 
968  for(timeorder=1;timeorder<numtimeorders;timeorder++){ // for >=RK2
969  fluxdt[timeorder-1] += CUnew[timeorder][2]*CUf[timeorder][1]*CUf[timeorder-1][2]*localdt;
970  }
971 
972  for(timeorder=0;timeorder<numtimeorders;timeorder++){ // for >=RK1
973  fluxdt[timeorder] += (CUnew[timeorder][1] + CUnew[timeorder][2]*CUf[timeorder][2])*localdt;
974  }
975 
976  if(nstep==0){
977  for(timeorder=0;timeorder<numtimeorders;timeorder++){
978  dualfprintf(fail_file,"timeorder=%d fluxdt/dt=%g\n",timeorder,fluxdt[timeorder]/localdt);
979  }
980  }
981 
982 
983  // assuming fluxdt used before any calls in advance.c
984  // then represents *amount* of flux'ed stuff in time dt
985  // This is NOT the time of the flux evaluation
986  // This is to be used with diag_flux and other things multiplied by dt for diagnostics only
987  // Other things (e.g. sources()) use same dt as dUriemann automatically, so fluxdt is not for anything but diagnostics
988  // Only exception is updating vpot that is out of place and uses ucum-type updates
989 
990 
991  // loop up to current substep to get current fluxdt
992  Ui=dUgeom=0.0; // don't care about update from non-flux terms
993  FTYPE uf;
994  FTYPE dUnongeomall[MAXTIMEORDER]={0.0};
995  for(timeorder=0;timeorder<numtimeorders;timeorder++){
996  if(timeorder==0){
997  oldufdt=0.0;
998  }
999  else{
1000  oldufdt=ufdt[timeorder-1];
1001  }
1002 
1003  // follows dUtoU() in advance.c
1004  ufdt[timeorder] = UFSET(CUf[timeorder],localdt,Ui,oldufdt,dUriemann,dUgeom,dUnongeomall);
1005  // below is NOT += since just want current change, not all prior changes
1006  // if did +=, then get 1 for timeorder==numtimeorders-1 as required!
1007  ucumdt[timeorder] = UCUMUPDATE(CUnew[timeorder],localdt,Ui,ufdt[timeorder],dUriemann,dUgeom,dUnongeomall);
1008 
1009  // time of pf at end of substep
1010  boundtime[timeorder] = t + ufdt[timeorder];
1011 
1012  // time of pb when used to compute flux
1013  // and time of pb corresponds to time of previous pf (except for first substep where it's just t)
1014  // tstepparti = t + CUf[timeorder][3] * dt;
1015  // tsteppartf = t + CUf[timeorder][2] * dt + CUf[timeorder][3] * dt;
1016  if(timeorder>0) fluxtime[timeorder] = t + CUf[timeorder-1][2] * dt + CUf[timeorder-1][3] * dt;
1017  else fluxtime[timeorder] = t;
1018 
1019  }
1020 
1021 
1022 
1023 
1024 #if(0)
1025 
1026  //
1027  // Get boundtime
1028  //
1029  // Get current time for pf
1030  // The below works because pi=p always and so always just use
1031  // assuming bound called after advance(), then need to get time of next flux evaluation
1032  // This corresponds to time where pb *will be* located
1033  //
1034  //
1036  for(timeorder=0;timeorder<numtimeorders;timeorder++){
1037  if(timeorder<numtimeorders-1){
1038  boundtime[timeorder] = t + localdt*CUf[timeorder][3]+ localdt*CUf[timeorder][2];
1039  }
1040  else{
1041  // if timeorder==numtimeorders-1, then final step and bound should be set for time of flux of next step, which is full t+dt
1042  boundtime[timeorder] = t + localdt;
1043  }
1044  }
1045 #endif
1046 
1047 
1048 #if(0)
1049  // DEBUG:
1050  dualfprintf(fail_file,"FLUXDT/BOUNDTIME: nstep=%ld ",nstep);
1051  for(timeorder=0;timeorder<numtimeorders;timeorder++){
1052  dualfprintf(fail_file,"to=%d fluxdt=%21.15g fluxdtperdt=%21.15g boundtime=%21.15g fluxtime=%21.15g\n",timeorder,fluxdt[timeorder],fluxdt[timeorder]/dt,boundtime[timeorder],fluxtime[timeorder]);
1053  }
1054  dualfprintf(fail_file,"\n");
1055 #endif
1056 
1057 
1058 }
1059 
1060 
1061 
1062 
1063 
1064 
1065 
1066 
1100 int step_ch_supermpi(int truestep, int *dumpingnext, FTYPE *fullndt, FTYPE (*prim)[NSTORE2][NSTORE3][NPR], FTYPE (*pstag)[NSTORE2][NSTORE3][NPR], FTYPE (*ucons)[NSTORE2][NSTORE3][NPR], FTYPE (*vpot)[NSTORE1+SHIFTSTORE1][NSTORE2+SHIFTSTORE2][NSTORE3+SHIFTSTORE3], FTYPE (*Bhat)[NSTORE2][NSTORE3][NPR], FTYPE (*pl_ct)[NSTORE1][NSTORE2][NSTORE3][NPR2INTERP], FTYPE (*pr_ct)[NSTORE1][NSTORE2][NSTORE3][NPR2INTERP],FTYPE (*F1)[NSTORE2][NSTORE3][NPR+NSPECIAL],FTYPE (*F2)[NSTORE2][NSTORE3][NPR+NSPECIAL],FTYPE (*F3)[NSTORE2][NSTORE3][NPR+NSPECIAL],FTYPE (*Atemp)[NSTORE1+SHIFTSTORE1][NSTORE2+SHIFTSTORE2][NSTORE3+SHIFTSTORE3],FTYPE (*uconstemp)[NSTORE2][NSTORE3][NPR])
1101 {
1102  int advance(int truestep, int stage, FTYPE (*pi)[NSTORE2][NSTORE3][NPR],FTYPE (*pb)[NSTORE2][NSTORE3][NPR], FTYPE (*pf)[NSTORE2][NSTORE3][NPR],
1103  FTYPE (*pstag)[NSTORE2][NSTORE3][NPR],
1104  FTYPE (*pl_ct)[NSTORE1][NSTORE2][NSTORE3][NPR2INTERP], FTYPE (*pr_ct)[NSTORE1][NSTORE2][NSTORE3][NPR2INTERP],
1105  FTYPE (*F1)[NSTORE2][NSTORE3][NPR+NSPECIAL],FTYPE (*F2)[NSTORE2][NSTORE3][NPR+NSPECIAL],FTYPE (*F3)[NSTORE2][NSTORE3][NPR+NSPECIAL],
1106  FTYPE (*vpot)[NSTORE1+SHIFTSTORE1][NSTORE2+SHIFTSTORE2][NSTORE3+SHIFTSTORE3],
1107  FTYPE (*ui)[NSTORE2][NSTORE3][NPR],FTYPE (*uf)[NSTORE2][NSTORE3][NPR], FTYPE (*ucum)[NSTORE2][NSTORE3][NPR],
1108  FTYPE *CUf,FTYPE *CUnew,SFTYPE fluxdt, SFTYPE boundtime, SFTYPE fluxtime, int timeorder, int numtimeorders, FTYPE *ndt);
1109  int boundstage;
1110  SFTYPE mydt;
1111  int stage, stagei,stagef;
1112  int timeorder;
1113  FTYPE ndt,lastndt;
1114  FTYPE (*pi)[NSTORE2][NSTORE3][NPR];
1115  FTYPE (*pb)[NSTORE2][NSTORE3][NPR];
1116  FTYPE (*pf)[NSTORE2][NSTORE3][NPR];
1117  FTYPE (*prevpf)[NSTORE2][NSTORE3][NPR];
1118  FTYPE (*pii[MAXTIMEORDER])[NSTORE2][NSTORE3][NPR];
1119  FTYPE (*pbb[MAXTIMEORDER])[NSTORE2][NSTORE3][NPR];
1120  FTYPE (*pff[MAXTIMEORDER])[NSTORE2][NSTORE3][NPR];
1121  FTYPE (*uii[MAXTIMEORDER])[NSTORE2][NSTORE3][NPR];
1122  FTYPE (*uff[MAXTIMEORDER])[NSTORE2][NSTORE3][NPR];
1123  FTYPE (*ucum[MAXTIMEORDER])[NSTORE2][NSTORE3][NPR];
1124  // FTYPE alphaik[MAXSTAGES][MAXSTAGES],betaik[MAXSTAGES];
1125  FTYPE CUf[MAXTIMEORDER][NUMDTCUFS],CUnew[MAXTIMEORDER][NUMDTCUFS];
1126  int i, j, k, pl, pliter;
1127  int numtimeorders;
1128  int finalstep;
1129  // extern int horizon_flux(FTYPE (*F1)[NSTORE2][NSTORE3][NPR+NSPECIAL], SFTYPE Dt);
1130  SFTYPE fluxdt[MAXTIMEORDER], boundtime[MAXTIMEORDER], fluxtime[MAXTIMEORDER];
1131 
1132 
1133 
1134 
1135 
1137  //
1138  // setup time-stepping
1139  //
1141  setup_rktimestep(truestep, &numtimeorders,prim,pstag,ucons,vpot,Bhat,GLOBALPOINT(pk),pii,pbb,pff,uii,uff,ucum,CUf,CUnew);
1142 
1143 
1145  //
1146  // Obtain initial time of substep, final time of substep, and true dt used for flux conservation that is used to iterate ucum in advance.c
1147  //
1149  get_truetime_fluxdt(numtimeorders, dt, CUf, CUnew, fluxdt, boundtime, fluxtime, NULL, NULL);
1150 
1151 
1152 
1153  // SPECIAL BOUNDARY/COMPUTATION MPI METHOD (out of date, and doesn't yet work right even if essentially complete code)
1154  /* check timestep */
1155  // if (dt < MINDT) {
1156  // trifprintf( "timestep too small\n");
1157  // myexit(11);
1158  // }
1159 
1160  lastndt=1E30; // initialize lastndt
1161  for(timeorder=1;timeorder<=numtimeorders;timeorder++){
1162 
1163 
1165  //
1166  //the starting and the ending times of the current substep
1167  //
1169  if(timeorder!=numtimeorders-1){
1170  tstepparti = t + CUf[timeorder][3] * dt;
1171  tsteppartf = t + CUf[timeorder][2] * dt + CUf[timeorder][3] * dt;
1172  }
1173  else{
1174  tstepparti=t;
1175  tsteppartf=t+dt;
1176  }
1177 
1178  if(timeorder==numtimeorders-1) finalstep=1; else finalstep=0;
1179 
1180 
1181 
1182 #if(PRODUCTION==0)
1183  trifprintf("-to%d/%d-",timeorder,numtimeorders);
1184 #endif
1185  if(numtimeorders==2){
1186  // note that pb (used for flux calc which has a stencil
1187  // calculation) must be different from pf so new stencils in
1188  // different stages won't affect stencil calculations -- must
1189  // use old values, not new from most previous temporary stage
1190  //
1191  // pi however can be the same as pf since each pi is replaced 1
1192  // zone at a time with a 0 stencil.
1193  if(timeorder==1){
1194  pi=prim;
1195  pb=prim;
1196  pf=GLOBALPOINT(pk)[0]; // different already, so good for simulbccalc
1197  prevpf=prim; // previous final true array
1198  mydt=0.5*dt;
1199  }
1200  else if(timeorder==2){
1201  pi=prim;
1202  pb=GLOBALPOINT(pk)[0];
1203  pf=prim;
1204  prevpf=GLOBALPOINT(pk)[0];
1205  mydt=dt;
1206  }
1207  }
1208  else if(numtimeorders==1){
1209  pi=prim;
1210  pb=prim;
1211  if(SIMULBCCALC<=0) pf=prim; else pf=GLOBALPOINT(pk)[0]; // need to be different if doing simulbccalc
1212  prevpf=prim;
1213  mydt=dt;
1214  }
1215  if(SIMULBCCALC<=0){ stagei=STAGEM1; stagef=STAGEM1; }
1216  else if(SIMULBCCALC==1) { stagei=STAGE0; stagef=STAGE2;}
1217  else if(SIMULBCCALC==2) { stagei=STAGE0; stagef=STAGE5;}
1218 
1219 
1220  // initialize bound stage
1221  if(SIMULBCCALC>=1) boundstage=STAGE0;
1222  else boundstage=STAGEM1;
1223  for(stage=stagei;stage<=stagef;stage++){
1224 #if(PRODUCTION==0)
1225  if(SIMULBCCALC>=1) trifprintf("!s%d!",stage);
1226 #endif
1227  // setup stage loop
1228 #if(SIMULBCCALC==2)
1229 #if(TYPE2==1)
1230  // GODMARK: isf1, etc. are NOT defined?!
1231  STAGECONDITION(0,N1-1,0,N2-1,isc,iec,jsc,jec);
1232  STAGECONDITION(0,N1,-1,N2,isf1,ief1,jsf1,jef1);
1233  STAGECONDITION(-1,N1,0,N2,isf2,ief2,jsf2,jef2);
1234  STAGECONDITION(0,N1,0,N2,ise,iee,jse,jee);
1235  STAGECONDITION(0,N1,0,N2-1,isf1ct,ief1ct,jsf1ct,jef1ct);
1236  STAGECONDITION(0,N1-1,0,N2,isf2ct,ief2ct,jsf2ct,jef2ct);
1237  STAGECONDITION(-1,N1,-1,N2,isdq,iedq,jsdq,jedq);
1238  STAGECONDITION(-2,N1+1,-2,N2+1,ispdq,iepdq,jspdq,jepdq);
1239  // GODMARK : probably not right for general boundary condition size
1240 #endif
1241 #endif
1242 
1243  // only bounding if safe zones, unsafe needs bz complete
1244  if(stage<STAGE2){
1245  bound_beforeevolveprim(boundstage, finalstep, boundtime[timeorder], prevpf,pstag,ucons,USEMPI); // pstag,ucons not right for supermpi
1246  if(stage!=STAGEM1) boundstage++;
1247  }
1248 
1249  // done here instead of local since pseudo-complicated calculation that might slow the dq calculation if done locally per zone
1250  MYFUN(pre_fixup(stage, prevpf),"step_ch.c:advance()", "pre_fixup()", 1);
1251 
1252  // go from previous solution to new solution
1253  partialstep=timeorder;
1254  // not right for numtimeorders==4 // GODMARK
1255  // advance
1256  MYFUN(advance(truestep,-1, pii[timeorder], pbb[timeorder], pff[timeorder], pstag, pl_ct, pr_ct, F1, F2, F3, vpot, uii[timeorder], uff[timeorder], ucum[timeorder],CUf[timeorder], CUnew[timeorder], fluxdt[timeorder], boundtime[timeorder], fluxtime[timeorder], timeorder,numtimeorders,&ndt),"step_ch.c:step_ch_supermpi()", "advance()", 1);
1257  // must check before MPI operation (since asymmetries would desynchronize cpus)
1258  if(stage<STAGE2){
1259  MYFUN(error_check(1),"step_ch.c", "error_check", 1);
1260  }
1261  if(stage!=STAGEM1){
1262  if(stage<STAGE2){
1263  bound_evolveprim(boundstage,finalstep, boundtime[timeorder], prevpf,pstag,ucons,USEMPI);
1264  boundstage++;
1265  }
1266  }
1267  if(timeorder==numtimeorders){
1268  if(ndt>lastndt) ndt=lastndt; // don't change if last was lower
1269  else lastndt=ndt; // new is lower, keep it
1270  }
1271  }
1272  if(timeorder==numtimeorders){// only do on full step
1273  // find global minimum value of ndt over all cpus
1274  mpifmin(&ndt);
1275  }
1276  // done when all stages are completed, so stencil used doesn't matter
1277 
1278  MYFUN(post_fixup(STAGEM1,finalstep, boundtime[timeorder], pf,pb,ucons),"step_ch.c:advance()", "post_fixup()", 1);
1279  }
1280 
1281 
1282  // pass back the final minimal dt over all substeps
1283  *fullndt = lastndt;
1284 
1285 
1286  // copy the contents to the final working array
1287  if((numtimeorders==1)&&(SIMULBCCALC>=1)) COMPFULLLOOP PLOOP(pliter,pl) MACP0A1(prim,i,j,k,pl)=MACP0A1(pf,i,j,k,pl);
1288 
1289 
1290  /* done! */
1291  return (0);
1292 }
1293 
1294 
1295 
1296 
1297 
1298 
1299 
1300 
1317 void setup_rktimestep(int truestep, int *numtimeorders,
1318  FTYPE (*p)[NSTORE2][NSTORE3][NPR], FTYPE (*pstag)[NSTORE2][NSTORE3][NPR], FTYPE (*ucons)[NSTORE2][NSTORE3][NPR], FTYPE (*vpot)[NSTORE1+SHIFTSTORE1][NSTORE2+SHIFTSTORE2][NSTORE3+SHIFTSTORE3], FTYPE (*Bhat)[NSTORE2][NSTORE3][NPR],
1319  FTYPE (*pk)[NSTORE1][NSTORE2][NSTORE3][NPR],
1320  FTYPE (*pii[MAXTIMEORDER])[NSTORE2][NSTORE3][NPR],FTYPE (*pbb[MAXTIMEORDER])[NSTORE2][NSTORE3][NPR],FTYPE (*pff[MAXTIMEORDER])[NSTORE2][NSTORE3][NPR],
1321  FTYPE (*uii[MAXTIMEORDER])[NSTORE2][NSTORE3][NPR],FTYPE (*uff[MAXTIMEORDER])[NSTORE2][NSTORE3][NPR],FTYPE (*ucum[MAXTIMEORDER])[NSTORE2][NSTORE3][NPR],
1322  FTYPE (*CUf)[NUMDTCUFS],FTYPE (*CUnew)[NUMDTCUFS])
1323 {
1324 
1325 
1326  // initialize CUf and CUnew to be zero
1327  int ii,jj;
1328  for(ii=0;ii<MAXTIMEORDER;ii++){
1329  for(jj=0;jj<NUMDTCUFS;jj++){
1330  CUf[ii][jj]=CUnew[ii][jj]=0.0;
1331  }
1332  }
1333 
1334 
1335 
1337  //
1338  // EXPLICIT TIME STEPPING
1339  //
1341 
1342 
1343  if(TIMETYPE==TIMEEXPLICIT){
1344 
1345 
1346 
1347  // to avoid special copying of final pff->p, always use p as final pff
1348  if(TIMEORDER==4 && truestep){
1349  // RK4 stepping
1350  *numtimeorders=4;
1351 
1352  // Ui ulast dUe(pb) <timestuff> dUi(Uf)
1353  CUf[0][0]=1.0; CUf[0][1]=0.0; CUf[0][2]=0.5; CUf[0][3] = 0.0; CUf[0][4+0] = CUf[0][2];
1354  CUf[1][0]=1.0; CUf[1][1]=0.0; CUf[1][2]=0.5; CUf[1][3] = 0.0; CUf[1][4+1] = CUf[1][2];
1355  CUf[2][0]=1.0; CUf[2][1]=0.0; CUf[2][2]=1.0; CUf[2][3] = 0.0; CUf[2][4+2] = CUf[2][2];
1356  CUf[3][0]=1.0; CUf[3][1]=0.0; CUf[3][2]=1.0; CUf[3][3] = 0.0; CUf[3][4+3] = CUf[3][2];
1357 
1358  // Ui dUe(Ub) Uf <timestuff> dUi(Uf)
1359  CUnew[0][0]=1.0; CUnew[0][1]=1.0/6.0; CUnew[0][2]=0.0; CUnew[0][3] = 0.0; CUnew[0][4+0] = CUnew[0][1];
1360  CUnew[1][0]=0.0; CUnew[1][1]=1.0/3.0; CUnew[1][2]=0.0; CUnew[1][3] = 0.5; CUnew[1][4+1] = CUnew[1][1];
1361  CUnew[2][0]=0.0; CUnew[2][1]=1.0/3.0; CUnew[2][2]=0.0; CUnew[2][3] = 0.5; CUnew[2][4+2] = CUnew[2][1];
1362  CUnew[3][0]=0.0; CUnew[3][1]=1.0/6.0; CUnew[3][2]=0.0; CUnew[3][3] = 1.0; CUnew[3][4+3] = CUnew[3][1];
1363 
1364  //primitive values used for initial state, fluxes, final state (where you output)
1365  pii[0]=p; pbb[0]=p; pff[0]=pk[0]; // produces U1
1366  pii[1]=p; pbb[1]=pk[0]; pff[1]=pk[1]; // produces U2
1367  pii[2]=p; pbb[2]=pk[1]; pff[2]=pk[0]; // produces U3
1368  pii[3]=p; pbb[3]=pk[0]; pff[3]=p; // produces U4 (only dUe part used)
1369 
1370  // GODMARK: use of globals: just scratch anyways
1371  uii[0]=GLOBALPOINT(uinitialglobal); uff[0]=GLOBALPOINT(ulastglobal); ucum[0]=ucons;
1372  uii[1]=GLOBALPOINT(uinitialglobal); uff[1]=GLOBALPOINT(ulastglobal); ucum[1]=ucons;
1373  uii[2]=GLOBALPOINT(uinitialglobal); uff[2]=GLOBALPOINT(ulastglobal); ucum[2]=ucons;
1374  uii[3]=GLOBALPOINT(uinitialglobal); uff[3]=GLOBALPOINT(ulastglobal); ucum[3]=ucons;
1375 
1376  // GODMARK: note that pbstag (staggered field from conserved de-averaging and inversion to primitive no geometry version) is always same memory space and comes from operating on final inverted quantity (ulastglobal or ucons), so just use same quantity for now and avoid adding extra code for pbstag[]
1377  // pbstag[0]=pstagglobal;
1378  // pbstag[1]=pstagglobal;
1379  // pbstag[2]=pstagglobal;
1380  // pbstag[3]=pstagglobal;
1381 
1382  }
1383  else if(TIMEORDER==3 && truestep){
1384  // TVD optimal RK3 method as in Shu's report
1385  *numtimeorders=3;
1386 
1387  // Ui ulastglobal dUe(pb) <timestuff> dUi(Uf)
1388  CUf[0][0]=1.0; CUf[0][1]=0.0; CUf[0][2]=1.0; CUf[0][3] = 0.0; CUf[0][4+0] = CUf[0][2];
1389  CUf[1][0]=3.0/4.0; CUf[1][1]=1.0/4.0; CUf[1][2]=1.0/4.0; CUf[1][3] = 0.0; CUf[1][4+1] = CUf[1][2];
1390  CUf[2][0]=1.0/3.0; CUf[2][1]=2.0/3.0; CUf[2][2]=2.0/3.0; CUf[2][3] = 0.0; CUf[2][4+2] = CUf[2][2];
1391 
1392  // Ui dUe(Ub) Uf <timestuff> dUi(Uf)
1393  // ucons=U3
1394  CUnew[0][0]=0.0; CUnew[0][1]=0.0; CUnew[0][2]=0.0; CUnew[0][3] = 0.0; CUnew[0][4+0] = CUnew[0][1];
1395  CUnew[1][0]=0.0; CUnew[1][1]=0.0; CUnew[1][2]=0.0; CUnew[1][3] = 1.0; CUnew[1][4+1] = CUnew[1][1];
1396  CUnew[2][0]=0.0; CUnew[2][1]=0.0; CUnew[2][2]=1.0; CUnew[2][3] = 1.0/4.0; CUnew[2][4+2] = CUnew[2][1];
1397 
1398  //always starting the substeps from the initial time
1399  pii[0]=p; pbb[0]=p; pff[0]=pk[0]; // produces U1
1400  pii[1]=p; pbb[1]=pk[0]; pff[1]=pk[1]; // produces U2
1401  pii[2]=p; pbb[2]=pk[1]; pff[2]=p; // produces U3
1402 
1403  // GODMARK: use of globals: just scratch anyways
1404  uii[0]=GLOBALPOINT(uinitialglobal); uff[0]=GLOBALPOINT(ulastglobal); ucum[0]=ucons;
1405  uii[1]=GLOBALPOINT(uinitialglobal); uff[1]=GLOBALPOINT(ulastglobal); ucum[1]=ucons;
1406  uii[2]=GLOBALPOINT(uinitialglobal); uff[2]=GLOBALPOINT(ulastglobal); ucum[2]=ucons;
1407 
1408  }
1409  else if(TIMEORDER==2 && truestep){
1410 #if(1)
1411  // midpoint method
1412 
1413  *numtimeorders=2;
1414 
1415  // old ucons not used for this method (i.e. [?][1]=0)
1416  // Ui ulast dUe(pb) <timestuff> dUi(Uf)
1417  CUf[0][0]=1.0; CUf[0][1]=0.0; CUf[0][2]=0.5; CUf[0][3] = 0.0; CUf[0][4+0] = CUf[0][2];
1418  CUf[1][0]=1.0; CUf[1][1]=0.0; CUf[1][2]=1.0; CUf[1][3] = 0.0; CUf[1][4+1] = CUf[1][2];
1419 
1420  // Ui dUe(Ub) Uf <timestuff> dUi(Uf)
1421  // ucons=U2
1422  CUnew[0][0]=0.0; CUnew[0][1]=0.0; CUnew[0][2]=0.0; CUnew[0][3] = 0.0; CUnew[0][4+0] = CUnew[0][1];
1423  CUnew[1][0]=0.0; CUnew[1][1]=0.0; CUnew[1][2]=1.0; CUnew[1][3] = 0.5; CUnew[1][4+1] = CUnew[1][1];
1424 
1425  pii[0]=p; pbb[0]=p; pff[0]=pk[0];
1426  pii[1]=p; pbb[1]=pk[0]; pff[1]=p;
1427 
1428  // GODMARK: use of globals: just scratch anyways
1429  uii[0]=GLOBALPOINT(uinitialglobal); uff[0]=GLOBALPOINT(ulastglobal); ucum[0]=ucons;
1430  uii[1]=GLOBALPOINT(uinitialglobal); uff[1]=GLOBALPOINT(ulastglobal); ucum[1]=ucons;
1431 
1432 #else
1433  *numtimeorders=2;
1434  // TVD RK2 (Chi-Wang Shu 1997 - eq 4.10)
1435  // actually less robust than generic midpoint method above
1436 
1437  // Ui ulast dUe(pb) <timestuff> dUi(Uf)
1438  CUf[0][0]=1.0; CUf[0][1]=0.0; CUf[0][2]=1.0; CUf[0][3] = 0.0; CUf[0][4+0] = CUf[0][2];
1439  CUf[1][0]=0.5; CUf[1][1]=0.5; CUf[1][2]=0.5; CUf[1][3] = 0.0; CUf[1][4+1] = CUf[1][2];
1440 
1441  // Ui dUe(Ub) Uf <timestuff> dUi(Uf)
1442  // ucons=U2
1443  CUnew[0][0]=0.0; CUnew[0][1]=0.0; CUnew[0][2]=0.0; CUnew[0][3] = 0.0; CUnew[0][4+0] = CUnew[0][1];
1444  CUnew[1][0]=0.0; CUnew[1][1]=0.0; CUnew[1][2]=1.0; CUnew[1][3] = 1.0; CUnew[1][4+1] = CUnew[1][1];
1445 
1446  pii[0]=p; pbb[0]=p; pff[0]=pk[0];
1447  pii[1]=p; pbb[1]=pk[0]; pff[1]=p;
1448 
1449  // GODMARK: use of globals: just scratch anyways
1450  uii[0]=GLOBALPOINT(uinitialglobal); uff[0]=GLOBALPOINT(ulastglobal); ucum[0]=ucons;
1451  uii[1]=GLOBALPOINT(uinitialglobal); uff[1]=GLOBALPOINT(ulastglobal); ucum[1]=ucons;
1452 
1453 #endif
1454  }
1455  else if(TIMEORDER==1 || dt==0.0){ // dt==0.0 case is case when just passing through
1456  // Euler method
1457  *numtimeorders=1;
1458 
1459  // Ui ulast dUe(pb) <timestuff> dUi(Uf)
1460  CUf[0][0]=1.0; CUf[0][1]=0.0; CUf[0][2]=1.0; CUf[0][3] = 0.0; CUf[0][4] = CUf[0][2];
1461 
1462  // Ui dUe(Ub) Uf <timestuff> dUi(Uf)
1463  // ucons=U1
1464  CUnew[0][0]=0.0; CUnew[0][1]=0.0; CUnew[0][2]=1.0; CUnew[1][3] = 0.0; CUnew[1][4] = CUnew[0][1];
1465 
1466  pii[0]=p; pbb[0]=p; pff[0]=p;
1467 
1468  // GODMARK: use of globals: just scratch anyways
1469  uii[0]=GLOBALPOINT(uinitialglobal); uff[0]=GLOBALPOINT(ulastglobal); ucum[0]=ucons;
1470 
1471  }
1472  }
1473 
1474 
1475 
1477  //
1478  // IMPLICIT TIME STEPPING
1479  //
1480  // First stage is just implicit solution (no flux needed)
1481  // Last stage is just flux (no implicit source term needed (i.e. no M[timeorder-1] needed), but geometry or any other explicits still needed)
1482  //
1483  // to avoid special copying of final pff->p, always use p as final pff
1484  //
1486  if(TIMETYPE==TIMEIMPLICIT){
1487 
1488  if(TIMEORDER==5 && truestep){
1489  // IMEX3
1490  *numtimeorders=5;
1491 
1492  static FTYPE imp_alpha=0.24169426078821;
1493  static FTYPE imp_beta=0.06042356519705;
1494  static FTYPE imp_eta=0.12915286960590;
1495 
1496  // Ui ulast dUe(pb) <timestuff> dUi(Uf)
1497  CUf[0][0]=1.0; CUf[0][1]=0.0; CUf[0][2]=0.0; CUf[0][3] = 0.0; CUf[0][4+0] = imp_alpha; // M0 no flux step : U0 = [Un] + \alpha dt M0
1498  CUf[1][0]=2.0; CUf[1][1]=-1.0; CUf[1][2]=0.0; CUf[1][3] = 0.0; CUf[1][4+1] = imp_alpha; // M1 no flux step : U1 = [2Un - U0] + \alpha dt M1
1499  CUf[2][0]=1.0; CUf[2][1]=0.0; CUf[2][2]=1.0; CUf[2][3] = 0.0; CUf[2][4+1] = (1.0-imp_alpha); CUf[2][4+2] = imp_alpha; // M2 and F1 : U2 = [Un + dt F1 + (1-\alpha) dt M1] + \alpha dt M2
1500  CUf[3][0]=3.0/4.0; CUf[3][1]=1.0/4.0; CUf[3][2]=1.0/4.0; CUf[3][3] = 0.0; CUf[3][4+0] = imp_beta; CUf[3][4+1] = (-1.0+imp_alpha+4.0*imp_eta)/4.0; CUf[3][4+2] = (2.0-5.0*imp_alpha-4.0*(imp_beta+imp_eta))/4.0; CUf[3][4+3] = imp_alpha; // M3 and F2 : U3 = [(3/4)Un + (1/4)U2 + (dt/4)F2 + \beta dt M0 + (1/4)(-1+\alpha+4\eta) dt M1 + (1/4)(2-5\alpha-4(\beta+\eta)) dt M2] + \alpha dt M3
1501  CUf[4][0]=1.0/3.0; CUf[4][1]=2.0/3.0; CUf[4][2]=2.0/3.0; CUf[4][3] = 0.0; CUf[4][4+0] = (-2.0*imp_beta/3.0); CUf[4][4+1] = (1.0-4.0*imp_eta)/6.0; CUf[4][4+2] = (-1.0+4.0*imp_alpha+4.0*(imp_beta+imp_eta))/6.0; CUf[4][4+3] = 4.0*(1.0-imp_alpha)/6.0; // F3 : U4 = [(1/3)Un + (2/3)U3 + (2dt/3) F3 + (-2beta dt/3) M0 + (1/6)(1-4\eta) dt M1 + (1/6)(-1 + 4\alpha + 4(\beta+\eta)) dt M2 + (1/6)(4(1-\alpha))dt M3]
1502 
1503  // Ui dUe(Ub) Uf <timestuff> dUi(Uf)
1504  // NOTE: Could instead only use final U^{n+1} = U4 as possible for this case and setup above for final step, but current way is more local and really would require less memory (but not implemented to do that)
1505  CUnew[0][0]=1.0/3.0; CUnew[0][1]=0.0; CUnew[0][2]=0.0; CUnew[0][3] = 0.0; CUnew[0][4+0] = (-2.0*imp_beta/3.0); // + (1/3)Un + (-2\beta/3) dt M0
1506  CUnew[1][0]=0.0; CUnew[1][1]=0.0; CUnew[1][2]=0.0; CUnew[1][3] = 0.0; CUnew[1][4+1] = (1.0-4.0*imp_eta)/6.0; // + ((1-4\eta)/6) dt M1
1507  CUnew[2][0]=0.0; CUnew[2][1]=0.0; CUnew[2][2]=0.0; CUnew[2][3] = 0.0; CUnew[2][4+2] = (-1.0+4.0*imp_alpha+4.0*(imp_beta+imp_eta))/6.0; // + ( (-1 + 4\alpha + 4(\beta+\eta))/6) dt M2
1508  CUnew[3][0]=0.0; CUnew[3][1]=0.0; CUnew[3][2]=2.0/3.0; CUnew[3][3] = 0.0; CUnew[3][4+3] = 4.0*(1.0-imp_alpha)/6.0; // + (2/3)U3 + (4(1-\alpha)/6) dt M3
1509  CUnew[4][0]=0.0; CUnew[4][1]=2.0/3.0; CUnew[4][2]=0.0; CUnew[4][3] = 0.0; CUnew[4][4+4] = 0.0; // + (2/3)dt F3
1510 
1511  //primitive values used for initial state, fluxes, final state (where you output)
1512  pii[0]=p; pbb[0]=p; pff[0]=pk[0]; // produces U1
1513  pii[1]=p; pbb[1]=pk[0]; pff[1]=pk[1]; // produces U2
1514  pii[2]=p; pbb[2]=pk[1]; pff[2]=pk[0]; // produces U3
1515  pii[3]=p; pbb[3]=pk[0]; pff[3]=pk[1]; // produces U4 (only dUe part used)
1516  pii[4]=p; pbb[4]=pk[1]; pff[4]=p; // produces U5 (only dUe part used)
1517 
1518  // GODMARK: use of globals: just scratch anyways
1519  uii[0]=GLOBALPOINT(uinitialglobal); uff[0]=GLOBALPOINT(ulastglobal); ucum[0]=ucons;
1520  uii[1]=GLOBALPOINT(uinitialglobal); uff[1]=GLOBALPOINT(ulastglobal); ucum[1]=ucons;
1521  uii[2]=GLOBALPOINT(uinitialglobal); uff[2]=GLOBALPOINT(ulastglobal); ucum[2]=ucons;
1522  uii[3]=GLOBALPOINT(uinitialglobal); uff[3]=GLOBALPOINT(ulastglobal); ucum[3]=ucons;
1523  uii[4]=GLOBALPOINT(uinitialglobal); uff[4]=GLOBALPOINT(ulastglobal); ucum[4]=ucons;
1524 
1525  }
1526  else if(TIMEORDER==4 && truestep){
1527  // IMEX2B
1528  *numtimeorders=4;
1529 
1530  // Ui ulast dUe(pb) <timestuff> dUi(Uf)
1531  CUf[0][0]=1.0; CUf[0][1]=0.0; CUf[0][2]=0.0; CUf[0][3] = 0.0; CUf[0][4+0] = 1.0/4.0; // M0 no flux step : U0 = [Un] + (dt/4)M0
1532  CUf[1][0]=1.0; CUf[1][1]=0.0; CUf[1][2]=0.5; CUf[1][3] = 0.0; CUf[1][4+1] = 1.0/4.0; // F0 and M1 : U1 = [Un + (dt/2)F0] + (dt/4)M1
1533  CUf[2][0]=0.0; CUf[2][1]=1.0; CUf[2][2]=1.0/2.0; CUf[2][3] = 0.0; CUf[2][4+0] = 1.0/3.0; CUf[2][4+1] = 1.0/12.0; CUf[2][4+2] = 1.0/3.0; // F1 and M2 : U2 = [U1 + (dt/2)F1 + (dt/3)M0 + (dt/12)M1] + (dt/3)M2
1534  CUf[3][0]=1.0/3.0; CUf[3][1]=2.0/3.0; CUf[3][2]=1.0/3.0; CUf[3][3] = 0.0; CUf[3][4+0] = 1.0/9.0; CUf[3][4+1] = 1.0/9.0; CUf[3][4+2] = 1.0/9.0; // F2 no new implicit step // U3 = [(1/3)Un + (2/3)U2 + (dt/3)F2 + (dt/9)(M0 + M1 + M2)]
1535 
1536  // Ui dUe(Ub) Uf <timestuff> dUi(Uf)
1537  // NOTE: Could instead only use final U^{n+1} = U3 as possible for this case and setup above for final step, but current way is more local and really would require less memory (but not implemented to do that)
1538  CUnew[0][0]=1.0/3.0; CUnew[0][1]=0.0; CUnew[0][2]=0.0; CUnew[0][3] = 0.0; CUnew[0][4+0] = (1.0/9.0); // + (1/3)Un + (1/9)M0
1539  CUnew[1][0]=0.0; CUnew[1][1]=0.0; CUnew[1][2]=0.0; CUnew[1][3] = 0.0; CUnew[1][4+1] = (1.0/9.0); // + (1/9)dt M1
1540  CUnew[2][0]=0.0; CUnew[2][1]=0.0; CUnew[2][2]=2.0/3.0; CUnew[2][3] = 0.0; CUnew[2][4+2] = (1.0/9.0); // + (2/3)U2 + (1/9)dt M2
1541  CUnew[3][0]=0.0; CUnew[3][1]=1.0/3.0; CUnew[3][2]=0.0; CUnew[3][3] = 0.0; CUnew[3][4+3] = 0.0; // + (1/3)dt F2
1542 
1543  //primitive values used for initial state, fluxes, final state (where you output)
1544  pii[0]=p; pbb[0]=p; pff[0]=pk[0]; // produces U1
1545  pii[1]=p; pbb[1]=pk[0]; pff[1]=pk[1]; // produces U2
1546  pii[2]=p; pbb[2]=pk[1]; pff[2]=pk[0]; // produces U3
1547  pii[3]=p; pbb[3]=pk[0]; pff[3]=p; // produces U4 (only dUe part used)
1548 
1549  // GODMARK: use of globals: just scratch anyways
1550  uii[0]=GLOBALPOINT(uinitialglobal); uff[0]=GLOBALPOINT(ulastglobal); ucum[0]=ucons;
1551  uii[1]=GLOBALPOINT(uinitialglobal); uff[1]=GLOBALPOINT(ulastglobal); ucum[1]=ucons;
1552  uii[2]=GLOBALPOINT(uinitialglobal); uff[2]=GLOBALPOINT(ulastglobal); ucum[2]=ucons;
1553  uii[3]=GLOBALPOINT(uinitialglobal); uff[3]=GLOBALPOINT(ulastglobal); ucum[3]=ucons;
1554 
1555  }
1556  else if(TIMEORDER==3 && truestep){
1557  // IMEX2
1558  *numtimeorders=3;
1559  // static FTYPE sqrt2=1.414213562373095048802;
1560  static FTYPE imp_gamma = 0.2928932188134524755992; // 1.0 - 1.0/sqrt2;
1561 
1562  // Ui ulastglobal dUe(pb) <timestuff> dUi(Uf)
1563  CUf[0][0]=1.0; CUf[0][1]=0.0; CUf[0][2]=0.0; CUf[0][3] = 0.0; CUf[0][4+0] = imp_gamma; // M0 no flux step : U0 = [Un] + \gamma dt M0
1564  CUf[1][0]=(3.0*imp_gamma-1.0)/imp_gamma; CUf[1][1]=(1.0-2.0*imp_gamma)/imp_gamma; CUf[1][2]=1.0; CUf[1][3] = 0.0; CUf[1][4+1] = imp_gamma; // computes F0 and M1 : U1 = [ (1/gamma)(3\gamma-1) Un + (1/\gamma)(1-2\gamma) U0 + dt F0] + \gamma dt M1
1565  CUf[2][0]=1.0/2.0; CUf[2][1]=1.0/2.0; CUf[2][2]=1.0/2.0; CUf[2][3] = 0.0; CUf[2][4+0] = imp_gamma; CUf[2][4+1]=(1.0-imp_gamma)/2.0; // F1 no new implicit step : U2 = [(1/2)Un + (1/2)U1 + (dt/2)F1 + \gamma dt M0 + (1/2)(1-\gamma) dt M1]
1566 
1567  // Ui dUe(Ub) Uf <timestuff> dUi(Uf)
1568  // ucons=U3
1569  // NOTE: Could instead only use final U^{n+1} = U2 as possible for this case and setup above for final step, but current way is more local and really would require less memory (but not implemented to do that)
1570  CUnew[0][0]=0.5; CUnew[0][1]=0.0; CUnew[0][2]=0.0; CUnew[0][3] = 0.0; CUnew[0][4+0]=imp_gamma; // + (1/2)Un + dt\gamma M0
1571  CUnew[1][0]=0.0; CUnew[1][1]=0.0; CUnew[1][2]=0.5; CUnew[1][3] = 0.0; CUnew[1][4+1]=(1.0-imp_gamma)/2.0; // + (1/2)U1 + (1/2)(1-\gamma)M1
1572  CUnew[2][0]=0.0; CUnew[2][1]=1.0/2.0; CUnew[2][2]=0.0; CUnew[2][3] = 0.0; // + (1/2)dt F1
1573 
1574  //always starting the substeps from the initial time
1575  pii[0]=p; pbb[0]=p; pff[0]=pk[0]; // produces U1
1576  pii[1]=p; pbb[1]=pk[0]; pff[1]=pk[1]; // produces U2
1577  pii[2]=p; pbb[2]=pk[1]; pff[2]=p; // produces U3
1578 
1579  // GODMARK: use of globals: just scratch anyways
1580  uii[0]=GLOBALPOINT(uinitialglobal); uff[0]=GLOBALPOINT(ulastglobal); ucum[0]=ucons;
1581  uii[1]=GLOBALPOINT(uinitialglobal); uff[1]=GLOBALPOINT(ulastglobal); ucum[1]=ucons;
1582  uii[2]=GLOBALPOINT(uinitialglobal); uff[2]=GLOBALPOINT(ulastglobal); ucum[2]=ucons;
1583 
1584  }
1585  else if(TIMEORDER==2 && truestep){
1586  // IMEX1B = Backward Euler (1st order implicit, 2nd order explicit)
1587  *numtimeorders=2;
1588 
1589  // old ucons not used for this method (i.e. [?][1]=0)
1590  // Ui ulast dUe(pb) <timestuff> dUi(Uf)
1591  CUf[0][0]=1.0; CUf[0][1]=0.0; CUf[0][2]=0.5; CUf[0][3] = 0.0; CUf[0][4+0] = CUf[0][2];
1592  CUf[1][0]=1.0; CUf[1][1]=0.0; CUf[1][2]=1.0; CUf[1][3] = 0.0; CUf[1][4+1] = CUf[1][2];
1593 
1594  // Ui dUe(Ub) Uf <timestuff> dUi(Uf)
1595  // ucons=U2
1596  CUnew[0][0]=0.0; CUnew[0][1]=0.0; CUnew[0][2]=0.0; CUnew[0][3] = 0.0; CUnew[0][4+0] = CUnew[0][1];
1597  CUnew[1][0]=0.0; CUnew[1][1]=0.0; CUnew[1][2]=1.0; CUnew[1][3] = 0.0; CUnew[1][4+1] = CUnew[1][1];
1598 
1599  pii[0]=p; pbb[0]=p; pff[0]=pk[0];
1600  pii[1]=p; pbb[1]=pk[0]; pff[1]=p;
1601 
1602  // GODMARK: use of globals: just scratch anyways
1603  uii[0]=GLOBALPOINT(uinitialglobal); uff[0]=GLOBALPOINT(ulastglobal); ucum[0]=ucons;
1604  uii[1]=GLOBALPOINT(uinitialglobal); uff[1]=GLOBALPOINT(ulastglobal); ucum[1]=ucons;
1605  }
1606  else if(TIMEORDER==1 || dt==0.0){ // dt==0.0 case is case when just passing through
1607  // Euler method = IMEX1 (1st order implicit, 1st order explicit)
1608  *numtimeorders=1;
1609 
1610  // Ui ulast dUe(pb) <timestuff> dUi(Uf)
1611  CUf[0][0]=1.0; CUf[0][1]=0.0; CUf[0][2]=1.0; CUf[0][3] = 0.0; CUf[0][4] = CUf[0][2];
1612 
1613  // Ui dUe(Ub) Uf <timestuff> dUi(Uf)
1614  // ucons=U1
1615  CUnew[0][0]=0.0; CUnew[0][1]=0.0; CUnew[0][2]=1.0; CUnew[1][3] = 0.0; CUnew[1][4] = CUnew[0][1];
1616 
1617  pii[0]=p; pbb[0]=p; pff[0]=p;
1618 
1619  // GODMARK: use of globals: just scratch anyways
1620  uii[0]=GLOBALPOINT(uinitialglobal); uff[0]=GLOBALPOINT(ulastglobal); ucum[0]=ucons;
1621 
1622  }
1623 
1624 
1625  }
1626 
1627 
1628 
1629 
1630 }
1631 
1632 
1633 
1634 
1635 
1636 
1637 
1638 
1639 
1641 //
1642 // System boundary routines that calls other system routines and user routines
1643 //
1645 
1646 
1647 
1648 // Who calls bound_anypstag():
1649 // v1)
1650 // step_ch.c: bound_evolveprim calls bound_allprim calls bound_anyallprim calls bound_anypstag
1651 // step_ch.c: bound_beforeevolveprim calls bound_anyallprim calls bound_anypstag
1652 // fluxctstag.c : bound_pstag calls bound_anypstag
1653 //
1654 // v2)
1655 // But if bound pstag in fluxctstag.c so Bcent can be interpolated from Bstag, then don't need to do it again!
1656 // fluxctstag.c : bound_pstag calls bound_anypstag
1657 
1658 
1659 
1660 
1663 int bound_evolveprim(int boundstage, int finalstep, SFTYPE boundtime, FTYPE (*prim)[NSTORE2][NSTORE3][NPR], FTYPE (*pstag)[NSTORE2][NSTORE3][NPR], FTYPE (*ucons)[NSTORE2][NSTORE3][NPR], int doboundmpi)
1664 {
1665  int boundvartype=BOUNDPRIMTYPE;
1666 
1667 
1668 
1669  bound_anyprim(boundstage, finalstep, boundtime, boundvartype, prim, pstag, ucons,doboundmpi);
1670  if(unewisavg && BOUNDUNEW){
1671  // SUPERGODMARK:
1672  // if not outflow boundary, then can bound u as p
1673  // desirable since machine errors can be different and lead to secular changes
1674  // esp. in MPI
1675  bound_uavg(boundstage, finalstep, boundtime, boundvartype, ucons,pstag, ucons,doboundmpi);
1676  }
1677 
1678 
1679  // bound_allprim(boundstage,finalstep, boundtime, prim,pstag,ucons,doboundmpi);
1680 
1681  return(0);
1682 
1683 }
1684 
1685 
1688 int bound_beforeevolveprim(int boundstage,int finalstep, SFTYPE boundtime, FTYPE (*prim)[NSTORE2][NSTORE3][NPR], FTYPE (*pstag)[NSTORE2][NSTORE3][NPR], FTYPE (*ucons)[NSTORE2][NSTORE3][NPR], int doboundmpi)
1689 {
1690  int boundvartype=BOUNDPRIMSIMPLETYPE; // tells boundary routines that only bounding 1 layer deep
1691 
1692 
1693  bound_anyprim(boundstage, finalstep, boundtime, boundvartype, prim, pstag, ucons,doboundmpi);
1694  if(unewisavg && BOUNDUNEW){
1695  // SUPERGODMARK:
1696  // if not outflow boundary, then can bound u as p
1697  // desirable since machine errors can be different and lead to secular changes
1698  // esp. in MPI
1699  bound_uavg(boundstage, finalstep, boundtime, boundvartype, ucons,pstag, ucons,doboundmpi);
1700  }
1701 
1702 
1703 
1704  return(0);
1705 
1706 }
1707 
1709 int bound_prim(int boundstage, int finalstep, SFTYPE boundtime, FTYPE (*prim)[NSTORE2][NSTORE3][NPR], FTYPE (*pstag)[NSTORE2][NSTORE3][NPR], FTYPE (*ucons)[NSTORE2][NSTORE3][NPR], int doboundmpi)
1710 {
1711  int boundvartype=BOUNDPRIMTYPE;
1712 
1713 
1714  bound_anyprim(boundstage, finalstep, boundtime, boundvartype, prim, pstag, ucons,doboundmpi);
1715 
1716  return(0);
1717 
1718 }
1719 
1721 int bound_pstag(int boundstage, int finalstep, SFTYPE boundtime, FTYPE (*prim)[NSTORE2][NSTORE3][NPR], FTYPE (*pstag)[NSTORE2][NSTORE3][NPR], FTYPE (*ucons)[NSTORE2][NSTORE3][NPR], int doboundmpi)
1722 {
1723  int boundvartype=BOUNDPSTAGTYPE;
1724 
1725  bound_anypstag(boundstage, finalstep, boundtime, boundvartype, prim, pstag, ucons,doboundmpi);
1726 
1727  return(0);
1728 
1729 }
1730 
1731 
1733 int bound_allprim(int boundstage, int finalstep, SFTYPE boundtime, FTYPE (*prim)[NSTORE2][NSTORE3][NPR], FTYPE (*pstag)[NSTORE2][NSTORE3][NPR], FTYPE (*ucons)[NSTORE2][NSTORE3][NPR], int doboundmpi)
1734 {
1735  int boundvartype=BOUNDPRIMTYPE;
1736 
1737 
1738  bound_anyallprim(boundstage, finalstep, boundtime, boundvartype, prim, pstag,ucons,doboundmpi);
1739 
1740  return(0);
1741 
1742 }
1743 
1744 
1745 
1746 
1748 int bound_anyallprim(int boundstage, int finalstep, SFTYPE boundtime, int boundvartype, FTYPE (*prim)[NSTORE2][NSTORE3][NPR], FTYPE (*pstag)[NSTORE2][NSTORE3][NPR], FTYPE (*ucons)[NSTORE2][NSTORE3][NPR],int doboundmpi)
1749 {
1750  int mystagboundvar;
1751 
1752 
1753 
1754  if(FLUXB==FLUXCTSTAG){
1755  bound_anyprim(boundstage, finalstep, boundtime, boundvartype, prim, pstag, ucons,doboundmpi);
1756  if(boundvartype==BOUNDPRIMTYPE) mystagboundvar=BOUNDPSTAGTYPE;
1757  else if(boundvartype==BOUNDPRIMSIMPLETYPE) mystagboundvar=BOUNDPSTAGSIMPLETYPE;
1758  else mystagboundvar=boundvartype;
1759  bound_anypstag(boundstage, finalstep, boundtime, mystagboundvar, prim, pstag, ucons,doboundmpi);
1760  }
1761  else{
1762  bound_anyprim(boundstage, finalstep, boundtime, boundvartype, prim, pstag, ucons,doboundmpi);
1763  }
1764 
1765 
1766  if(unewisavg && BOUNDUNEW){
1767  // SUPERGODMARK:
1768  // if not outflow boundary, then can bound u as p
1769  // desirable since machine errors can be different and lead to secular changes
1770  // esp. in MPI
1771  bound_uavg(boundstage, finalstep, boundtime, boundvartype, ucons,pstag, ucons,doboundmpi);
1772  }
1773 
1774  return(0);
1775 
1776 }
1777 
1778 
1779 
1781 int bound_uavg(int boundstage, int finalstep, SFTYPE boundtime, int boundvartype, FTYPE (*prim)[NSTORE2][NSTORE3][NPR], FTYPE (*pstag)[NSTORE2][NSTORE3][NPR], FTYPE (*ucons)[NSTORE2][NSTORE3][NPR], int doboundmpi)
1782 {
1783  int mystagboundvar;
1784  FTYPE (*uavg)[NSTORE2][NSTORE3][NPR];
1785 
1786 
1787 
1788 
1789  // assign uavg
1790  uavg=ucons;
1791 
1792 
1793  // for now only worry about bounding staggered fields since that's what can do given "general" boundary condition code
1794  if(1||DOENOFLUX != NOENOFLUX){
1795 
1796  // CHANGINGMARK: DEBUG:
1797  // can modify for diag call if choose to avoid if outflow condition
1798  // or do something simple for outflow just for diagnostics
1799 #if(FULLOUTPUT!=0 && PRODUCTION==0)
1800  bound_anyprim(boundstage, finalstep, boundtime, boundvartype, uavg,pstag, uavg,doboundmpi);
1801 #endif
1802 
1803 
1804  if(DOENOFLUX == ENOFINITEVOLUME){
1805  // then need to bound all conservative quantities
1806  // SUPERGODMARK -- CHANGINGMARK -- need to generalize so bound knows if p or U quantity
1807  // other methods have only field "averaged" so only need to modify it
1808  bound_anyprim(boundstage, finalstep, boundtime, boundvartype, uavg, pstag, uavg,doboundmpi);
1809  }
1810 
1811  if(FLUXB==FLUXCTSTAG){
1812  // bound unew for staggered fields
1813  // SUPERGODMARK -- CHANGINGMARK: need to tell bound if p or U
1814  if(boundvartype==BOUNDPRIMTYPE) mystagboundvar=BOUNDPSTAGTYPE;
1815  else if(boundvartype==BOUNDPRIMSIMPLETYPE) mystagboundvar=BOUNDPSTAGSIMPLETYPE;
1816  else mystagboundvar=boundvartype;
1817 
1818  bound_anypstag(boundstage, finalstep, boundtime, mystagboundvar, uavg, pstag, uavg,doboundmpi);
1819  }
1820  }
1821 
1822  return(0);
1823 
1824 }
1825 
1826 
1827 
1828 
1829 
1834 int bound_anyprim(int boundstage, int finalstep, SFTYPE boundtime, int boundvartype, FTYPE (*prim)[NSTORE2][NSTORE3][NPR], FTYPE (*pstag)[NSTORE2][NSTORE3][NPR], FTYPE (*ucons)[NSTORE2][NSTORE3][NPR], int doboundmpi)
1835 {
1836  int whichpltoavg[NPR];
1837  int ifnotavgthencopy[NPR];
1838  int nprlocalstart,nprlocalend;
1839  int nprlocallist[MAXNPR];
1840  int pl,pliter;
1841  int dir;
1842 
1843 
1844 
1845 
1846  if(DOGRIDSECTIONING){
1847  if((boundstage==STAGE0)||(boundstage==STAGEM1)){
1848  MYFUN(bound_gridsectioning(CENTEREDPRIM,prim,pstag,ucons,finalstep),"step_ch.c:bound_pstag()", "bound_pstag_user()", 1);
1849  }
1850  }
1851 
1852 
1853 
1854  DIMENLOOP(dir){
1855 
1856  if(0 && FLUXB==FLUXCTSTAG){ // apparently need to bound both fields SUPERGODMARK
1857 
1858  PALLLOOP(pl) whichpltoavg[pl]=1;
1859  PLOOPBONLY(pl) whichpltoavg[pl]=0;
1860 
1861  PALLLOOP(pl) ifnotavgthencopy[pl]=0;
1862  PLOOPBONLY(pl) ifnotavgthencopy[pl]=0;
1863 
1864  addremovefromanynpr(REMOVEFROMNPR, whichpltoavg, ifnotavgthencopy, &nprboundstart, &nprboundend, nprboundlist, &nprlocalstart, &nprlocalend, nprlocallist, NULL, NULL);
1865  }
1866  else{
1867  // no change then, doing all variables
1868  }
1869 
1870 
1871 
1872  // pre-post MPI recv's
1873  if(doboundmpi){
1874  MYFUN(bound_mpi_dir(boundstage, finalstep, -dir, boundvartype, prim, NULL, NULL, NULL,NULL),"step_ch.c:bound_prim()", "bound_mpi()", 1);
1875  }
1876 
1877  // real boundary zones
1878  if((boundstage==STAGE0)||(boundstage==STAGEM1)){
1879  MYFUN(bound_prim_user_dir(boundstage, finalstep, boundtime, dir, boundvartype, prim),"step_ch.c:bound_prim()", "bound_prim_user()", 1);
1880  }// end if stage0 or stagem1
1881 
1882  if(doboundmpi){
1883  MYFUN(bound_mpi_dir(boundstage, finalstep, +dir, boundvartype, prim, NULL, NULL, NULL,NULL),"step_ch.c:bound_prim()", "bound_mpi()", 1);
1884  }
1885 
1886  // real boundary zones
1887  if((boundstage==STAGE0)||(boundstage==STAGEM1)){
1888  int ispstag=BOUNDPRIMLOC;
1889  MYFUN(bound_prim_user_after_mpi_dir(boundstage, finalstep, boundtime, dir, boundvartype, ispstag, prim),"step_ch.c:bound_prim()", "bound_prim_user_after_mpi()", 1);
1890  }// end if stage0 or stagem1
1891 
1892 
1894  //
1895  // restore choice for bounding
1896 
1897  if(0 && FLUXB==FLUXCTSTAG){ // apparently need to bound both fields SUPERGODMARK
1898  addremovefromanynpr(RESTORENPR, whichpltoavg, ifnotavgthencopy, &nprboundstart, &nprboundend, nprboundlist, &nprlocalstart, &nprlocalend, nprlocallist, NULL, NULL);
1899  }
1900  }
1901 
1902 
1903  return(0);
1904 }
1905 
1906 
1912 int bound_anypstag(int boundstage, int finalstep, SFTYPE boundtime, int boundvartype, FTYPE (*prim)[NSTORE2][NSTORE3][NPR], FTYPE (*pstag)[NSTORE2][NSTORE3][NPR], FTYPE (*ucons)[NSTORE2][NSTORE3][NPR], int doboundmpi)
1913 {
1914  int pl,pliter;
1915  int nprlocalstart,nprlocalend;
1916  int nprlocallist[MAXNPR];
1917  int dir;
1918  int mystagboundvar;
1919 
1920 
1921  if(boundvartype==BOUNDPRIMTYPE) mystagboundvar=BOUNDPSTAGTYPE;
1922  else if(boundvartype==BOUNDPRIMSIMPLETYPE) mystagboundvar=BOUNDPSTAGSIMPLETYPE;
1923  else mystagboundvar=boundvartype;
1924 
1925 
1926 
1927 
1928 
1929  if(DOGRIDSECTIONING){
1930  if((boundstage==STAGE0)||(boundstage==STAGEM1)){
1931  MYFUN(bound_gridsectioning(STAGGEREDPRIM,prim,pstag,ucons,finalstep),"step_ch.c:bound_pstag()", "bound_pstag_user()", 1);
1932  }
1933  }
1934 
1935 
1936 
1937 
1938  DIMENLOOP(dir){
1939 
1941  //
1942  // save choice for bounding
1943  nprlocalstart=nprboundstart;
1944  nprlocalend=nprboundend;
1945  PMAXNPRLOOP(pl) nprlocallist[pl]=nprboundlist[pl];
1946 
1947  // choose range of PBOUNDLOOP and PLOOPMPI
1948  nprboundstart=0;
1949  nprboundend=2;
1950  nprboundlist[0]=B1;
1951  nprboundlist[1]=B2;
1952  nprboundlist[2]=B3;
1953 
1954 
1955  // pre-post MPI recvs
1956  if(doboundmpi){
1957  MYFUN(bound_mpi_dir(boundstage, finalstep, -dir, mystagboundvar, pstag, NULL, NULL, NULL,NULL),"step_ch.c:bound_pstag()", "bound_mpi()", 1);
1958  }
1959 
1960  // real boundary zones
1961  if((boundstage==STAGE0)||(boundstage==STAGEM1)){
1962  MYFUN(bound_pstag_user_dir(boundstage, finalstep, boundtime, dir,mystagboundvar,pstag),"step_ch.c:bound_pstag()", "bound_pstag_user()", 1);
1963  }// end if stage0 or stagem1
1964 
1965 
1966  if(doboundmpi){
1967  MYFUN(bound_mpi_dir(boundstage, finalstep, +dir, mystagboundvar, pstag, NULL, NULL, NULL,NULL),"step_ch.c:bound_pstag()", "bound_mpi()", 1);
1968  }
1969 
1970  // real boundary zones
1971  if((boundstage==STAGE0)||(boundstage==STAGEM1)){
1972  int ispstag=BOUNDPSTAGLOC;
1973  MYFUN(bound_prim_user_after_mpi_dir(boundstage, finalstep, boundtime, dir, mystagboundvar, ispstag, pstag),"step_ch.c:bound_pstag()", "bound_prim_user_after_mpi()", 1);
1974  }// end if stage0 or stagem1
1975 
1976 
1978  //
1979  // restore choice for bounding
1980  nprboundstart=nprlocalstart;
1981  nprboundend=nprlocalend;
1982  PMAXNPRLOOP(pl) nprboundlist[pl]=nprlocallist[pl];
1983  }
1984 
1985 
1986  return(0);
1987 }
1988 
1989 
1994 int bound_flux(int boundstage, int finalstep, SFTYPE boundtime, int boundvartype, FTYPE (*F1)[NSTORE2][NSTORE3][NPR+NSPECIAL], FTYPE (*F2)[NSTORE2][NSTORE3][NPR+NSPECIAL], FTYPE (*F3)[NSTORE2][NSTORE3][NPR+NSPECIAL], int doboundmpi)
1995 {
1996  int dir;
1997 
1998 
1999 
2000  if(boundvartype==BOUNDFLUXTYPE || boundvartype==BOUNDFLUXSIMPLETYPE){
2001  // then can stay
2002  }
2003  else{
2004  dualfprintf(fail_file,"Shouldn't be in bound_flux() with boundvartype=%d\n",boundvartype);
2005  myexit(2467367);
2006  }
2007 
2008 
2009  if(DOENOFLUX==ENOFLUXRECON && BOUNDFLUXRECON==0){
2010  dualfprintf(fail_file,"DEBUG: Assuming bound_flux() called only for debugging purposes\n");
2011  }
2012 
2013 
2014  // pre-post recv's
2015  if(doboundmpi){
2016  MYFUN(bound_mpi(boundstage, finalstep, -1, boundvartype, NULL, F1, F2, F3, NULL),"step_ch.c:bound_flux()", "bound_mpi()", 1);
2017  }
2018 
2019 
2020  // real boundary zones
2021  if((boundstage==STAGE0)||(boundstage==STAGEM1)){
2022  MYFUN(bound_flux_user(boundstage, finalstep, boundtime, boundvartype,F1,F2,F3),"step_ch.c:bound_flux()", "bound_flux_user()", 1);
2023  }// end if stage0 or stagem1
2024 
2025 
2026  if(doboundmpi){
2027  MYFUN(bound_mpi(boundstage, finalstep, +1, boundvartype, NULL, F1, F2, F3, NULL),"step_ch.c:bound_flux()", "bound_mpi()", 1);
2028  }
2029 
2030 
2031  return(0);
2032 }
2033 
2034 
2036 int bound_vpot(int boundstage, int finalstep, SFTYPE boundtime, int boundvartype, FTYPE (*vpot)[NSTORE1+SHIFTSTORE1][NSTORE2+SHIFTSTORE2][NSTORE3+SHIFTSTORE3], int doboundmpi, int doboundnonmpi)
2037 {
2038  int dir;
2039 
2040 
2041 
2042  if(boundvartype==BOUNDVPOTTYPE || boundvartype==BOUNDVPOTSIMPLETYPE){
2043  // then can stay
2044  }
2045  else{
2046  dualfprintf(fail_file,"Shouldn't be in bound_vpot() with boundvartype=%d\n",boundvartype);
2047  myexit(3483466);
2048  }
2049 
2050 
2051  // pre-post MPI recv's
2052  if(doboundmpi){
2053  MYFUN(bound_mpi(boundstage, finalstep, -1, boundvartype, NULL, NULL, NULL, NULL, vpot),"step_ch.c:bound_flux()", "bound_mpi()", 1);
2054  }
2055 
2056 
2057  // real boundary zones
2058  if(doboundnonmpi){
2059  if((boundstage==STAGE0)||(boundstage==STAGEM1)){
2060  MYFUN(bound_vpot_user(boundstage, finalstep, boundtime, boundvartype,vpot),"step_ch.c:bound_vpot()", "bound_vpot_user()", 1);
2061  }// end if stage0 or stagem1
2062  }
2063 
2064  if(doboundmpi){
2065  MYFUN(bound_mpi(boundstage, finalstep, +1, boundvartype, NULL, NULL, NULL, NULL, vpot),"step_ch.c:bound_flux()", "bound_mpi()", 1);
2066  }
2067 
2068 
2069  return(0);
2070 }
2071 
2072 
2073 
2074 
2075 
2076 
2077 
2080 int bound_pflag(int boundstage, int finalstep, SFTYPE boundtime, PFTYPE (*prim)[NSTORE2][NSTORE3][NUMPFLAGS], int doboundmpi)
2081 {
2082  int boundvartype=BOUNDINTTYPE;
2083 
2084 
2085 
2086 
2087  if(doboundmpi){
2088 
2089  if(UTOPRIMFIXMPICONSISTENT==1){
2090  MYFUN(bound_mpi_int(boundstage, finalstep, -1, boundvartype, prim),"step_ch.c:bound_pflag()", "bound_mpi_int()", 1);
2091  }
2092  else{
2093  // need to fill boundary cells with failure
2094  // otherwise, would have to go deeper into fixups and dependency chain for the UTOPRIMFIXMPICONSISTENT chocie would become too deep
2095  MYFUN(bound_mpi_int_fakeutoprimmpiinconsisent(boundstage, finalstep, -1, boundvartype, prim,UTOPRIMFAILFAKEVALUE),"step_ch.c:bound_pflag()", "bound_mpi_int()", 1);
2096  }
2097 
2098  }
2099 
2100 
2101 
2102 
2103  if((boundstage==STAGE0)||(boundstage==STAGEM1)){
2104 
2105  MYFUN(bound_pflag_user(boundstage, finalstep, boundtime, boundvartype, prim),"step_ch.c:bound_pflag()", "bound_pflag_user()", 1);
2106 
2107  }// end bound stage
2108 
2109 
2110  if(doboundmpi){
2111 
2112  if(UTOPRIMFIXMPICONSISTENT==1){
2113  MYFUN(bound_mpi_int(boundstage, finalstep, +1, boundvartype, prim),"step_ch.c:bound_pflag()", "bound_mpi_int()", 1);
2114  }
2115  else{
2116  // need to fill boundary cells with failure
2117  // otherwise, would have to go deeper into fixups and dependency chain for the UTOPRIMFIXMPICONSISTENT chocie would become too deep
2118  MYFUN(bound_mpi_int_fakeutoprimmpiinconsisent(boundstage, finalstep, +1, boundvartype, prim,UTOPRIMFAILFAKEVALUE),"step_ch.c:bound_pflag()", "bound_mpi_int()", 1);
2119  }
2120 
2121  }
2122 
2123  return(0);
2124 
2125 }
2126 
2127 
2128 
2129 
2130 
2131 
2132 
2133 
2134 
Generated on Fri May 20 2016 15:52:35 for HARM by   doxygen 1.8.3.1