HARM
harm and utilities
 All Data Structures Files Functions Variables Typedefs Macros Pages
bounds.tools.c
Go to the documentation of this file.
1 
13 #include "decs.h"
14 
15 // see definit.h or init.h for some defines here.
16 
17 
24 int bound_x1dn_analytic(int boundstage, int finalstep, SFTYPE boundtime, int whichdir, int boundvartype, int *dirprim, int ispstag, FTYPE (*prim)[NSTORE2][NSTORE3][NPR])
25 {
26  int i,j,k;
27  int pl,pliter;
28 
29 #if(ANALYTICMEMORY==0)
30  dualfprintf(fail_file,"ANALYTICMEMORY==0 but called bound_x1dn_analytic\n");
31  myexit(786752);
32 #endif
33 
34 
35  // then no need to set pglobal or pstagglobal since evolution sets appropriately in constrained way on a well-defined computational box
36  // dualfprintf(fail_file,"No use for this currently\n");
37  // myexit(34269834);
38 
39 
40 
41  if(ispstag){
43  // note we assume "i=0" is boundary cell to be fixed
44  // This ensures divb=0, but may be inconsistent with code's treatement of true i=0 if doing OUTFLOW
46  pl=B1; MACP0A1(prim,i,j,k,pl) = GLOBALMACP0A1(pstaganalytic,i,j,k,pl);
47  }
48  if(WITHINACTIVEBNDSECTIONX1DN(i,j,k)){
49  PBOUNDLOOP(pliter,pl) if(pl!=B1) MACP0A1(prim,i,j,k,pl) = GLOBALMACP0A1(pstaganalytic,i,j,k,pl);
50  }
51  }
52  }
53  else{
55  if(WITHINACTIVEBNDSECTIONX1DN(i,j,k)){
56  PBOUNDLOOP(pliter,pl) MACP0A1(prim,i,j,k,pl) = GLOBALMACP0A1(panalytic,i,j,k,pl);
57  }
58  }
59 
60  }
61 
62 
63  return(0);
64 }
65 
66 
73 int bound_x1up_analytic(int boundstage, int finalstep, SFTYPE boundtime, int whichdir, int boundvartype, int *dirprim, int ispstag, FTYPE (*prim)[NSTORE2][NSTORE3][NPR])
74 {
75  int i,j,k;
76  int pl,pliter;
77 
78 
79 #if(ANALYTICMEMORY==0)
80  dualfprintf(fail_file,"ANALYTICMEMORY==0 but called bound_x1up_analytic\n");
81  myexit(786753);
82 #endif
83 
84 
85  // KORALTODO: Check below
86  // then no need to set pglobal or pstagglobal since evolution sets appropriately in constrained way on a well-defined computational box
87  // dualfprintf(fail_file,"No use for this currently\n");
88  // myexit(34269834);
89 
90 
91  if(ispstag){
94  pl=B1; MACP0A1(prim,i,j,k,pl) = GLOBALMACP0A1(pstaganalytic,i,j,k,pl);
95  }
96  if(WITHINACTIVEBNDSECTIONX1UP(i,j,k)){
97  PBOUNDLOOP(pliter,pl) if(pl!=B1) MACP0A1(prim,i,j,k,pl) = GLOBALMACP0A1(pstaganalytic,i,j,k,pl);
98  }
99  }
100  }
101  else{
102  COMPFULLLOOP{
103  if(WITHINACTIVEBNDSECTIONX1UP(i,j,k)){
104  PBOUNDLOOP(pliter,pl) MACP0A1(prim,i,j,k,pl) = GLOBALMACP0A1(panalytic,i,j,k,pl);
105  }
106  }
107 
108  }
109 
110 
111  return(0);
112 }
113 
114 
121 int bound_x2dn_analytic(int boundstage, int finalstep, SFTYPE boundtime, int whichdir, int boundvartype, int *dirprim, int ispstag, FTYPE (*prim)[NSTORE2][NSTORE3][NPR])
122 {
123  int i,j,k;
124  int pl,pliter;
125 
126 #if(ANALYTICMEMORY==0)
127  dualfprintf(fail_file,"ANALYTICMEMORY==0 but called bound_x2dn_analytic\n");
128  myexit(786752);
129 #endif
130 
131 
132  // KORALTODO: Check if below true.
133  // then no need to set pglobal or pstagglobal since evolution sets appropriately in constrained way on a well-defined computational box
134  // dualfprintf(fail_file,"No use for this currently\n");
135  // myexit(34269834);
136 
137 
138 
139  if(ispstag){
140  COMPFULLLOOP{
141  // note we assume "j=0" is boundary cell to be fixed
142  // This ensures divb=0, but may be inconsistent with code's treatement of true j=0 if doing OUTFLOW
144  pl=B2; MACP0A1(prim,i,j,k,pl) = GLOBALMACP0A1(pstaganalytic,i,j,k,pl);
145  }
146  if(WITHINACTIVEBNDSECTIONX2DN(i,j,k)){
147  PBOUNDLOOP(pliter,pl) if(pl!=B2) MACP0A1(prim,i,j,k,pl) = GLOBALMACP0A1(pstaganalytic,i,j,k,pl);
148  }
149  }
150  }
151  else{
152  COMPFULLLOOP{
153  if(WITHINACTIVEBNDSECTIONX2DN(i,j,k)){
154  PBOUNDLOOP(pliter,pl) MACP0A1(prim,i,j,k,pl) = GLOBALMACP0A1(panalytic,i,j,k,pl);
155  }
156  }
157 
158  }
159 
160 
161  return(0);
162 }
163 
164 
171 int bound_x2up_analytic(int boundstage, int finalstep, SFTYPE boundtime, int whichdir, int boundvartype, int *dirprim, int ispstag, FTYPE (*prim)[NSTORE2][NSTORE3][NPR])
172 {
173  int i,j,k;
174  int pl,pliter;
175 
176 
177 #if(ANALYTICMEMORY==0)
178  dualfprintf(fail_file,"ANALYTICMEMORY==0 but called bound_x2up_analytic\n");
179  myexit(786753);
180 #endif
181 
182 
183 
184  // then no need to set pglobal or pstagglobal since evolution sets appropriately in constrained way on a well-defined computational box
185  // dualfprintf(fail_file,"No use for this currently\n");
186  // myexit(34269834);
187 
188 
189  if(ispstag){
190  COMPFULLLOOP{
192  pl=B2; MACP0A1(prim,i,j,k,pl) = GLOBALMACP0A1(pstaganalytic,i,j,k,pl);
193  }
194  if(WITHINACTIVEBNDSECTIONX2UP(i,j,k)){
195  PBOUNDLOOP(pliter,pl) if(pl!=B2) MACP0A1(prim,i,j,k,pl) = GLOBALMACP0A1(pstaganalytic,i,j,k,pl);
196  }
197  }
198  }
199  else{
200  COMPFULLLOOP{
201  if(WITHINACTIVEBNDSECTIONX2UP(i,j,k)){
202  PBOUNDLOOP(pliter,pl) MACP0A1(prim,i,j,k,pl) = GLOBALMACP0A1(panalytic,i,j,k,pl);
203  }
204  }
205 
206  }
207 
208 
209  return(0);
210 }
211 
212 
213 
220 int bound_x3dn_analytic(int boundstage, int finalstep, SFTYPE boundtime, int whichdir, int boundvartype, int *dirprim, int ispstag, FTYPE (*prim)[NSTORE2][NSTORE3][NPR])
221 {
222  int i,j,k;
223  int pl,pliter;
224 
225 #if(ANALYTICMEMORY==0)
226  dualfprintf(fail_file,"ANALYTICMEMORY==0 but called bound_x3dn_analytic\n");
227  myexit(786752);
228 #endif
229 
230 
231  // KORALTODO: Check if below true.
232  // then no need to set pglobal or pstagglobal since evolution sets appropriately in constrained way on a well-defined computational box
233  // dualfprintf(fail_file,"No use for this currently\n");
234  // myexit(34269834);
235 
236 
237 
238  if(ispstag){
239  COMPFULLLOOP{
240  // note we assume "k=0" is boundary cell to be fixed
241  // This ensures divb=0, but may be inconsistent with code's treatement of true k=0 if doing OUTFLOW
243  pl=B3; MACP0A1(prim,i,j,k,pl) = GLOBALMACP0A1(pstaganalytic,i,j,k,pl);
244  }
245  if(WITHINACTIVEBNDSECTIONX3DN(i,j,k)){
246  PBOUNDLOOP(pliter,pl) if(pl!=B3) MACP0A1(prim,i,j,k,pl) = GLOBALMACP0A1(pstaganalytic,i,j,k,pl);
247  }
248  }
249  }
250  else{
251  COMPFULLLOOP{
252  if(WITHINACTIVEBNDSECTIONX3DN(i,j,k)){
253  PBOUNDLOOP(pliter,pl) MACP0A1(prim,i,j,k,pl) = GLOBALMACP0A1(panalytic,i,j,k,pl);
254  }
255  }
256 
257  }
258 
259 
260  return(0);
261 }
262 
263 
270 int bound_x3up_analytic(int boundstage, int finalstep, SFTYPE boundtime, int whichdir, int boundvartype, int *dirprim, int ispstag, FTYPE (*prim)[NSTORE2][NSTORE3][NPR])
271 {
272  int i,j,k;
273  int pl,pliter;
274 
275 
276 #if(ANALYTICMEMORY==0)
277  dualfprintf(fail_file,"ANALYTICMEMORY==0 but called bound_x3up_analytic\n");
278  myexit(786753);
279 #endif
280 
281 
282 
283  // then no need to set pglobal or pstagglobal since evolution sets appropriately in constrained way on a well-defined computational box
284  // dualfprintf(fail_file,"No use for this currently\n");
285  // myexit(34269834);
286 
287 
288  if(ispstag){
289  COMPFULLLOOP{
291  pl=B3; MACP0A1(prim,i,j,k,pl) = GLOBALMACP0A1(pstaganalytic,i,j,k,pl);
292  }
293  if(WITHINACTIVEBNDSECTIONX3UP(i,j,k)){
294  PBOUNDLOOP(pliter,pl) if(pl!=B3) MACP0A1(prim,i,j,k,pl) = GLOBALMACP0A1(pstaganalytic,i,j,k,pl);
295  }
296  }
297  }
298  else{
299  COMPFULLLOOP{
300  if(WITHINACTIVEBNDSECTIONX3UP(i,j,k)){
301  PBOUNDLOOP(pliter,pl) MACP0A1(prim,i,j,k,pl) = GLOBALMACP0A1(panalytic,i,j,k,pl);
302  }
303  }
304 
305  }
306 
307 
308  return(0);
309 }
310 
311 
312 
314 int bound_x1dn_outflow_simple(
315  int boundstage, int finalstep, SFTYPE boundtime, int whichdir, int boundvartype, int *dirprim, int ispstag, FTYPE (*prim)[NSTORE2][NSTORE3][NPR],
316  int *inboundloop,
317  int *outboundloop,
318  int *innormalloop,
319  int *outnormalloop,
320  int (*inoutlohi)[NUMUPDOWN][NDIM],
321  int riin, int riout, int rjin, int rjout, int rkin, int rkout,
322  int *dosetbc,
323  int enerregion,
324  int *localenerpos
325  )
326 {
327 
328 
329 
330 #pragma omp parallel // assume don't require EOS
331  {
332  int i,j,k,pl,pliter;
333  struct of_geom geom[NPR],rgeom[NPR];
334  FTYPE vcon[NDIM]; // coordinate basis vcon
335 #if(WHICHVEL==VEL3)
336  int failreturn;
337 #endif
338  int ri, rj, rk; // reference i,j,k
339  FTYPE prescale[NPR];
340  int jj,kk;
341  struct of_geom geomdontuse[NPR];
342  struct of_geom *ptrgeom[NPR];
343  struct of_geom rgeomdontuse[NPR];
344  struct of_geom *ptrrgeom[NPR];
345 
346 
347  // assign memory
348  PALLLOOP(pl){
349  ptrgeom[pl]=&(geomdontuse[pl]);
350  ptrrgeom[pl]=&(rgeomdontuse[pl]);
351  }
352 
353 
354  if(((BCtype[X1DN]==OUTFLOW || BCtype[X1DN]==HORIZONOUTFLOW))||(BCtype[X1DN]==FIXEDOUTFLOW)||(BCtype[X1DN]==FREEOUTFLOW)){
355 
356 
357  if ( (totalsize[1]>1) && (mycpupos[1] <= 0)) {
358  /* inner r boundary condition: u, just copy */
359 
360  // PBOUNDLOOP(pliter,pl) dualfprintf(fail_file,"WTF1A: pl=%d pliter=%d : %d %d %d : nstep=%ld %d %g\n",pl,pliter,nprboundstart,nprboundend,nprboundlist[pliter],nstep,steppart,t);
361  // dualfprintf(fail_file,"WTF1A: ri=%d loopin=%d %d\n",riin,inoutlohi[POINTDOWN][POINTDOWN][1],inoutlohi[POINTDOWN][POINTUP][1]);
362 
363 
366 #pragma omp for schedule(OPENMPSCHEDULE(),OPENMPCHUNKSIZE(blocksize))
369 
370 
371 
372  ri=riin;
373  rj=j;
374  rk=k;
375 
376 
377  PALLLOOP(pl) get_geometry(ri, rj, rk, dirprim[pl], ptrrgeom[pl]);
378 
379  LOOPBOUND1IN{ // bound entire region inside non-evolved portion of grid
380  PBOUNDLOOP(pliter,pl) MACP0A1(prim,i,j,k,pl) = MACP0A1(prim,ri,rj,rk,pl);
381  //if(MACP0A1(prim,i,j,k,URAD1)>=0.0) MACP0A1(prim,i,j,k,URAD1)=0.0;
382  }
383 
384 
385 
386 
387  if(ispstag==0){
388 
389  if(((BCtype[X1DN]==OUTFLOW || BCtype[X1DN]==HORIZONOUTFLOW))||(BCtype[X1DN]==FIXEDOUTFLOW)){
390  // GODMARK: assume all velocities at same location when doing inflow check
392 #if(WHICHVEL==VEL4)
393  get_geometry(i, j, k, dirprim[U1], ptrgeom[U1]);
394  inflow_check_4vel(1,MAC(prim,i,j,k),NULL,ptrgeom[U1], 0) ;
395 #elif(WHICHVEL==VEL3)
396  get_geometry(i, j, k, dirprim[U1], ptrgeom[U1]);
397  inflow_check_3vel(1,MAC(prim,i,j,k),NULL,ptrgeom[U1], 0) ;
398  // projection may not preserve u^t to be real and rho>rhoscal u>uuscal
399 #if(JONCHECKS)
400  if(jonchecks){
401  //fixup1zone(0,MAC(prim,i,j,k),ptrgeom[U1],0);
402  failreturn=check_pr(MAC(prim,i,j,k),MAC(prim,i,j,k),ptrgeom[U1],-3);
403  if(failreturn){
404  dualfprintf(fail_file,"Bad boundary zone, couldn't fix: i=%d j=%d k=%d\n",startpos[1]+i,startpos[2]+j,startpos[3]+k);
405  if (fail(i,j,k,FAIL_BCFIX) >= 1) return (1);
406  }
407  }
408 #endif
409 #elif(WHICHVEL==VELREL4)
410  get_geometry(i,j,k,dirprim[U1],ptrgeom[U1]) ;
411  inflow_check_rel4vel(1,MAC(prim,i,j,k),NULL,ptrgeom[U1],0) ;
412  if(limit_gamma(0,GAMMAMAX,GAMMAMAXRAD,MAC(prim,i,j,k),NULL,ptrgeom[U1],0)>=1)
413  FAILSTATEMENT("bounds.c:bound_prim()", "limit_gamma()", 1);
414 #endif
415  }
416  } // end if not allowing inflow
417  }
418 
419 
420  }// end 2 3
421 
422  }// end if mycpupos[1]==0
423 
424 
425 
426  }
427  else{
428  dualfprintf(fail_file,"Shouldn't be here in bounds: dir=%d\n",X1DN);
429  myexit(3946836);
430  }
431 
432 
433  }// end parallel region
434 
435 
436  return(0);
437 }
438 
440 int bound_x1up_outflow_simple(
441  int boundstage, int finalstep, SFTYPE boundtime, int whichdir, int boundvartype, int *dirprim, int ispstag, FTYPE (*prim)[NSTORE2][NSTORE3][NPR],
442  int *inboundloop,
443  int *outboundloop,
444  int *innormalloop,
445  int *outnormalloop,
446  int (*inoutlohi)[NUMUPDOWN][NDIM],
447  int riin, int riout, int rjin, int rjout, int rkin, int rkout,
448  int *dosetbc,
449  int enerregion,
450  int *localenerpos
451  )
452 {
453 
454 
455 
456 
457 #pragma omp parallel // assume don't require EOS
458  {
459 
460  int i,j,k,pl,pliter;
461  FTYPE vcon[NDIM]; // coordinate basis vcon
462 #if(WHICHVEL==VEL3)
463  int failreturn;
464 #endif
465  int ri, rj, rk; // reference i,j,k
466  FTYPE prescale[NPR];
467  int jj,kk;
468  struct of_geom geomdontuse[NPR];
469  struct of_geom *ptrgeom[NPR];
470  struct of_geom rgeomdontuse[NPR];
471  struct of_geom *ptrrgeom[NPR];
472 
473  // assign memory
474  PALLLOOP(pl){
475  ptrgeom[pl]=&(geomdontuse[pl]);
476  ptrrgeom[pl]=&(rgeomdontuse[pl]);
477  }
478 
479 
480 
481  if(((BCtype[X1UP]==OUTFLOW || BCtype[X1UP]==HORIZONOUTFLOW))||(BCtype[X1UP]==FIXEDOUTFLOW)||(BCtype[X1UP]==FREEOUTFLOW)){
482 
483  // outer r BC:
484  if ( (totalsize[1]>1) && (mycpupos[1] == ncpux1 - 1) ) {
485 
486 
487  // PBOUNDLOOP(pliter,pl) dualfprintf(fail_file,"WTF2A: pl=%d pliter=%d : %d %d %d : nstep=%ld %d %g\n",pl,pliter,nprboundstart,nprboundend,nprboundlist[pliter],nstep,steppart,t);
488  // dualfprintf(fail_file,"WTF2A: ri=%d loop=%d %d\n",riout,inoutlohi[POINTUP][POINTDOWN][1],inoutlohi[POINTUP][POINTUP][1]);
489 
490 
493 #pragma omp for schedule(OPENMPSCHEDULE(),OPENMPCHUNKSIZE(blocksize))
496 
497 
498  ri=riout;
499  rj=j;
500  rk=k;
501 
502  PALLLOOP(pl) get_geometry(ri, rj, rk, dirprim[pl], ptrrgeom[pl]);
503 
504  LOOPBOUND1OUT PBOUNDLOOP(pliter,pl){
505  MACP0A1(prim,i,j,k,pl) = MACP0A1(prim,ri,rj,rk,pl);
506  }
507 
508 
509 
510 
511 
512  if(ispstag==0){
513 
514  if(((BCtype[X1UP]==OUTFLOW || BCtype[X1UP]==HORIZONOUTFLOW))||(BCtype[X1UP]==FIXEDOUTFLOW)){
515 
517 #if(WHICHVEL==VEL4)
518  get_geometry(i, j, k, dirprim[U1], ptrgeom[U1]);
519  inflow_check_4vel(1,MAC(prim,i,j,k),NULL,ptrgeom[U1],0) ;
520 #elif(WHICHVEL==VEL3)
521  get_geometry(i, j, k, dirprim[U1], ptrgeom[U1]);
522  inflow_check_3vel(1,MAC(prim,i,j,k),NULL,ptrgeom[U1],0) ;
523  // projection may not preserve u^t to be real and rho>rhoscal u>uuscal
524 #if(JONCHECKS)
525  if(jonchecks){
526  //fixup1zone(0,MAC(prim,i,j,k),ptrgeom[U1],0);
527  failreturn=check_pr(MAC(prim,i,j,k),MAC(prim,i,j,k),ptrgeom[U1],-3);
528  if(failreturn){
529  dualfprintf(fail_file,"Bad boundary zone, couldn't fix: i=%d j=%d k=%d\n",startpos[1]+i,startpos[2]+j,startpos[3]+k);
530  if (fail(i,j,k,FAIL_BCFIX) >= 1) return (1);
531  }
532  }
533 #endif
534 #elif(WHICHVEL==VELREL4)
535  get_geometry(i,j,k,dirprim[U1],ptrgeom[U1]) ;
536  inflow_check_rel4vel(1,MAC(prim,i,j,k),NULL,ptrgeom[U1],0) ;
537  if(limit_gamma(0,GAMMAMAX,GAMMAMAXRAD,MAC(prim,i,j,k),NULL,ptrgeom[U1], 0)>=1)
538  FAILSTATEMENT("bounds.c:bound_prim()", "limit_gamma()", 2);
539 #endif
540  }// end over i
541  }// end if not allowing inflow
542 
543  }
544 
545 
546 
547  }// end 2 3
548  }// end if mycpu is correct
549 
550  }
551  else{
552  dualfprintf(fail_file,"Shouldn't be here in bounds: dir=%d\n",X1UP);
553  myexit(3946837);
554  }
555 
556  }// end parallel region
557 
558 
559  return(0);
560 }
561 
562 
563 
564 
565 
566 
567 
568 
569 
570 
571 
572 
573 
574 
575 
576 
579  int boundstage, int finalstep, SFTYPE boundtime, int whichdir, int boundvartype, int *dirprim, int ispstag, FTYPE (*prim)[NSTORE2][NSTORE3][NPR],
580  int *inboundloop,
581  int *outboundloop,
582  int *innormalloop,
583  int *outnormalloop,
584  int (*inoutlohi)[NUMUPDOWN][NDIM],
585  int riin, int riout, int rjin, int rjout, int rkin, int rkout,
586  int *dosetbc,
587  int enerregion,
588  int *localenerpos
589  )
590 {
591 
592 
593 
594 #pragma omp parallel // assume don't require EOS
595  {
596  int i,j,k,pl,pliter;
597  struct of_geom geom[NPR],rgeom[NPR];
598  FTYPE vcon[NDIM]; // coordinate basis vcon
599 #if(WHICHVEL==VEL3)
600  int failreturn;
601 #endif
602  int ri, rj, rk; // reference i,j,k
603  FTYPE prescale[NPR];
604  int horizonti;
605  int jj,kk;
606  struct of_geom geomdontuse[NPR];
607  struct of_geom *ptrgeom[NPR];
608  struct of_geom rgeomdontuse[NPR];
609  struct of_geom *ptrrgeom[NPR];
610 
611 
612  // assign memory
613  PALLLOOP(pl){
614  ptrgeom[pl]=&(geomdontuse[pl]);
615  ptrrgeom[pl]=&(rgeomdontuse[pl]);
616  }
617 
618 
620 
621 
622  if ( (totalsize[1]>1) && (mycpupos[1] <= horizoncpupos1)) { // now all CPUs inside CPU with horizon will be using this (GODMARK: reference value needs to be chosen somehow for CPUs not on active grid)
623  /* inner r boundary condition: u, just copy */
624 
627 #pragma omp for schedule(OPENMPSCHEDULE(),OPENMPCHUNKSIZE(blocksize))
630 
631 
632 
633  ri=riin;
634  rj=j;
635  rk=k;
636 
637 
638 #if(HORIZONEXTRAP==0)
639  PALLLOOP(pl) get_geometry(ri, rj, rk, dirprim[pl], ptrrgeom[pl]);
640 
641  LOOPBOUND1INSPECIAL{ // bound entire region inside non-evolved portion of grid
642  PBOUNDLOOP(pliter,pl) MACP0A1(prim,i,j,k,pl) = MACP0A1(prim,ri,rj,rk,pl);
643  }
644 #elif(HORIZONEXTRAP==1)
645  PALLLOOP(pl) get_geometry(ri, rj, rk, dirprim[pl], ptrrgeom[pl]);
646 
647  LOOPBOUND1INSPECIAL{ // bound entire region inside non-evolved portion of grid
648 
649  PALLLOOP(pl) get_geometry(i, j, k, dirprim[pl], ptrgeom[pl]);
650 
651  PBOUNDLOOP(pliter,pl){
652  if(pl>=RHO && pl<=UU){
653  MACP0A1(prim,i,j,k,pl) = MACP0A1(prim,ri,rj,rk,pl) * (ptrrgeom[pl]->gdet/ptrgeom[pl]) ;
654  }
655  }
656  PBOUNDLOOP(pliter,pl){
657  if(pl==U1){ // treat U1 as special
658  MACP0A1(prim,i,j,k,pl) = MACP0A1(prim,ri,rj,rk,pl) * (1. - (i-ri)*dx[1]) ;
659  for(pl=U2;pl<=U3;pl++){
660  MACP0A1(prim,i,j,k,pl) = MACP0A1(prim,ri,rj,rk,pl) * (1. + (i-ri)*dx[1]) ;
661  }
662  }
663  }
664  PBOUNDLOOP(pliter,pl){
665  if(pl==B1){ // treat B1 special
666  MACP0A1(prim,i,j,k,pl) = MACP0A1(prim,ri,rj,rk,pl) * (ptrrgeom[pl]->gdet/ptrgeom[pl]);
667  for(pl=B2;pl<=B3;pl++){
668  MACP0A1(prim,i,j,k,pl) = MACP0A1(prim,ri,rj,rk,pl) * (1. + (i-ri)*dx[1]) ;
669  }
670  }
671  }
672  PBOUNDLOOP(pliter,pl){
673  if(pl>=B3+1 && pl<NPRBOUND){
674  MACP0A1(prim,i,j,k,pl) = MACP0A1(prim,ri,rj,rk,pl) * (ptrrgeom[pl]->gdet/ptrgeom[pl]);
675  }
676  }
677  }
678 #elif(HORIZONEXTRAP==2)
679  get_geometry(ri, rj, rk, dirprim[0], ptrrgeom[0]);
680  rescale(DORESCALE,1,MAC(prim,ri,rj,rk),ptrrgeom[0],prescale);
682  // set guess
683  PBOUNDLOOP(pliter,pl) MACP0A1(prim,i,j,k,pl)=MACP0A1(prim,ri,rj,k,pl);
684  get_geometry(i, j, k, dirprim[0], ptrgeom[0]);
685  rescale(UNRESCALE,1,MAC(prim,i,j,k),ptrgeom[0],prescale);
686  }
687 #elif(HORIZONEXTRAP==3)
688  extrapfunc(X1DN,j,k,boundstage,finalstep,boundtime,whichdir,boundvartype,dirprim,ispstag,prim,inboundloop,outboundloop,innormalloop,outnormalloop,inoutlohi,riin,riout,rjin,rjout,rkin,rkout,dosetbc,enerregion,localenerpos);
689 #endif
690 
691 
692 
693 
694  if(ispstag==0){
695 
697  // GODMARK: assume all velocities at same location when doing inflow check
699 #if(WHICHVEL==VEL4)
700  get_geometry(i, j, k, dirprim[U1], ptrgeom[U1]);
701  inflow_check_4vel(1,MAC(prim,i,j,k),NULL,ptrgeom[U1], 0) ;
702 #elif(WHICHVEL==VEL3)
703  get_geometry(i, j, k, dirprim[U1], ptrgeom[U1]);
704  inflow_check_3vel(1,MAC(prim,i,j,k),NULL,ptrgeom[U1], 0) ;
705  // projection may not preserve u^t to be real and rho>rhoscal u>uuscal
706 #if(JONCHECKS)
707  if(jonchecks){
708  //fixup1zone(0,MAC(prim,i,j,k),ptrgeom[U1],0);
709  failreturn=check_pr(MAC(prim,i,j,k),MAC(prim,i,j,k),ptrgeom[U1],-3);
710  if(failreturn){
711  dualfprintf(fail_file,"Bad boundary zone, couldn't fix: i=%d j=%d k=%d\n",startpos[1]+i,startpos[2]+j,startpos[3]+k);
712  if (fail(i,j,k,FAIL_BCFIX) >= 1) return (1);
713  }
714  }
715 #endif
716 #elif(WHICHVEL==VELREL4)
717  get_geometry(i,j,k,dirprim[U1],ptrgeom[U1]) ;
718  inflow_check_rel4vel(1,MAC(prim,i,j,k),NULL,ptrgeom[U1],0) ;
719  if(limit_gamma(0,GAMMAMAX,GAMMAMAXRAD,MAC(prim,i,j,k),NULL,ptrgeom[U1],0)>=1)
720  FAILSTATEMENT("bounds.c:bound_prim()", "limit_gamma()", 1);
721 #endif
722  }
723  } // end if not allowing inflow
724  }
725 
726 
727  }// end 2 3
728 
729  }// end if mycpupos[1]==0
730 
731 
732 
733  }
734  else{
735  dualfprintf(fail_file,"Shouldn't be here in bounds: horizonextrap version: dir=%d\n",X1DN);
736  myexit(3946838);
737  }
738 
739 
740  }// end parallel region
741 
742 
743  return(0);
744 }
745 
746 
747 
748 
749 
750 
751 
752 
753 
754 
755 
756 
759  int boundstage, int finalstep, SFTYPE boundtime, int whichdir, int boundvartype, int *dirprim, int ispstag, FTYPE (*prim)[NSTORE2][NSTORE3][NPR],
760  int *inboundloop,
761  int *outboundloop,
762  int *innormalloop,
763  int *outnormalloop,
764  int (*inoutlohi)[NUMUPDOWN][NDIM],
765  int riin, int riout, int rjin, int rjout, int rkin, int rkout,
766  int *dosetbc,
767  int enerregion,
768  int *localenerpos
769  )
770 {
771 
772 
773 
774 
775 #pragma omp parallel // assume don't require EOS
776  {
777 
778  int i,j,k,pl,pliter;
779  FTYPE vcon[NDIM]; // coordinate basis vcon
780 #if(WHICHVEL==VEL3)
781  int failreturn;
782 #endif
783  int ri, rj, rk; // reference i,j,k
784  FTYPE prescale[NPR];
785  int jj,kk;
786  struct of_geom geomdontuse[NPR];
787  struct of_geom *ptrgeom[NPR];
788  struct of_geom rgeomdontuse[NPR];
789  struct of_geom *ptrrgeom[NPR];
790 
791  // assign memory
792  PALLLOOP(pl){
793  ptrgeom[pl]=&(geomdontuse[pl]);
794  ptrrgeom[pl]=&(rgeomdontuse[pl]);
795  }
796 
797 
798 
800 
801  // outer r BC:
802  if ( (totalsize[1]>1) && (mycpupos[1] == ncpux1 - 1) ) {
803 
804 
807 #pragma omp for schedule(OPENMPSCHEDULE(),OPENMPCHUNKSIZE(blocksize))
810 
811 
812  ri=riout;
813  rj=j;
814  rk=k;
815 
816 
817 #if(OUTEREXTRAP==0)
818  PALLLOOP(pl) get_geometry(ri, rj, rk, dirprim[pl], ptrrgeom[pl]);
819 
820  LOOPBOUND1OUT PBOUNDLOOP(pliter,pl) MACP0A1(prim,i,j,k,pl) = MACP0A1(prim,ri,rj,rk,pl);
821 #elif(OUTEREXTRAP==1)
822  PALLLOOP(pl) get_geometry(ri, rj, rk, dirprim[pl], ptrrgeom[pl]);
823 
825  PALLLOOP(pl) get_geometry(i, j, k, dirprim[pl], ptrgeom[pl]);
826 
827  PBOUNDLOOP(pliter,pl){
828  if(pl>=RHO && pl<=UU){
829  MACP0A1(prim,i,j,k,pl) = MACP0A1(prim,ri,rj,rk,pl) * (ptrrgeom[pl]->gdet/ptrgeom[pl]) ;
830  }
831  }
832  PBOUNDLOOP(pliter,pl){
833  if(pl==U1){ // treat U1 as special
834  MACP0A1(prim,i,j,k,pl) = MACP0A1(prim,ri,rj,rk,pl) * (1. - 2*(i-ri)*dx[1]) ;
835  for(pl=U2;pl<=U3;pl++){
836  MACP0A1(prim,i,j,k,pl) = MACP0A1(prim,ri,rj,rk,pl) * (1. - (i-ri)*dx[1]) ;
837  }
838  }
839  }
840  PBOUNDLOOP(pliter,pl){
841  if(pl==B1){ // treat B1 special
842  MACP0A1(prim,i,j,k,pl) = MACP0A1(prim,ri,rj,rk,pl) * (ptrrgeom[pl]->gdet/ptrgeom[pl]) ;
843  for(pl=B2;pl<=B3;pl++){
844  MACP0A1(prim,i,j,k,pl) = MACP0A1(prim,ri,rj,rk,pl) * (1. - (i-ri)*dx[1]) ;
845  }
846  }
847  }
848  PBOUNDLOOP(pliter,pl){
849  if(pl>=B3+1 && pl<NPRBOUND){
850  MACP0A1(prim,i,j,k,pl) = MACP0A1(prim,ri,rj,rk,pl) * (ptrrgeom[pl]->gdet/ptrgeom[pl]) ;
851  }
852  }
853  }
854 #elif(OUTEREXTRAP==2)
855  get_geometry(ri, rj, rk, dirprim[0], ptrrgeom[0]);
856  rescale(DORESCALE,1,MAC(prim,ri,rj,rk),ptrrgeom[0],prescale);
858  // set guess
859  PBOUNDLOOP(pliter,pl) MACP0A1(prim,i,j,k,pl)=MACP0A1(prim,ri,rj,rk,pl);
860  get_geometry(i, j, k, dirprim[0], ptrgeom[0]);
861  rescale(UNRESCALE,1,MAC(prim,i,j,k),ptrgeom[0],prescale);
862  }
863 #elif(OUTEREXTRAP==3)
864  extrapfunc(X1UP,j,k,boundstage,finalstep,boundtime,whichdir,boundvartype,dirprim,ispstag,prim,inboundloop,outboundloop,innormalloop,outnormalloop,inoutlohi,riin,riout,rjin,rjout,rkin,rkout,dosetbc,enerregion,localenerpos);
865 #endif
866 
867 
868 
869 
870 
871  if(ispstag==0){
872 
874 
876 #if(WHICHVEL==VEL4)
877  get_geometry(i, j, k, dirprim[U1], ptrgeom[U1]);
878  inflow_check_4vel(1,MAC(prim,i,j,k),NULL,ptrgeom[U1],0) ;
879 #elif(WHICHVEL==VEL3)
880  get_geometry(i, j, k, dirprim[U1], ptrgeom[U1]);
881  inflow_check_3vel(1,MAC(prim,i,j,k),NULL,ptrgeom[U1],0) ;
882  // projection may not preserve u^t to be real and rho>rhoscal u>uuscal
883 #if(JONCHECKS)
884  if(jonchecks){
885  //fixup1zone(0,MAC(prim,i,j,k),ptrgeom[U1],0);
886  failreturn=check_pr(MAC(prim,i,j,k),MAC(prim,i,j,k),ptrgeom[U1],-3);
887  if(failreturn){
888  dualfprintf(fail_file,"Bad boundary zone, couldn't fix: i=%d j=%d k=%d\n",startpos[1]+i,startpos[2]+j,startpos[3]+k);
889  if (fail(i,j,k,FAIL_BCFIX) >= 1) return (1);
890  }
891  }
892 #endif
893 #elif(WHICHVEL==VELREL4)
894  get_geometry(i,j,k,dirprim[U1],ptrgeom[U1]) ;
895  inflow_check_rel4vel(1,MAC(prim,i,j,k),NULL,ptrgeom[U1],0) ;
896  if(limit_gamma(0,GAMMAMAX,GAMMAMAXRAD,MAC(prim,i,j,k),NULL,ptrgeom[U1], 0)>=1)
897  FAILSTATEMENT("bounds.c:bound_prim()", "limit_gamma()", 2);
898 #endif
899  }// end over i
900  }// end if not allowing inflow
901 
902  }
903 
904 
905 
906  }// end 2 3
907  }// end if mycpu is correct
908 
909  }
910  else{
911  dualfprintf(fail_file,"Shouldn't be here in bounds\n");
912  myexit(3946839);
913  }
914 
915  }// end parallel region
916 
917 
918  return(0);
919 }
920 
921 
922 
923 
924 
925 
928  int boundstage, int finalstep, SFTYPE boundtime, int whichdir, int boundvartype, int *dirprim, int ispstag, FTYPE (*prim)[NSTORE2][NSTORE3][NPR],
929  int *inboundloop,
930  int *outboundloop,
931  int *innormalloop,
932  int *outnormalloop,
933  int (*inoutlohi)[NUMUPDOWN][NDIM],
934  int riin, int riout, int rjin, int rjout, int rkin, int rkout,
935  int *dosetbc,
936  int enerregion,
937  int *localenerpos
938  )
939 {
940 
941 
942 #pragma omp parallel // assume don't require EOS
943  {
944  int i,j,k,pl,pliter;
945  FTYPE vcon[NDIM]; // coordinate basis vcon
946 #if(WHICHVEL==VEL3)
947  int failreturn;
948 #endif
949  int ri, rj, rk; // reference i,j,k
950  FTYPE prescale[NPR];
951  int jj,kk;
952 
953 
954 
955  if((BCtype[X1DN]==R0SING)||(BCtype[X1DN]==SYMM)||(BCtype[X1DN]==ASYMM) ){
956 
957 
958 
959  /* inner radial BC (preserves u^t rho and u) */
960  if ( (totalsize[1]>1) && (mycpupos[1] == 0) ) {
962 
963  { // start block
966 #pragma omp for schedule(OPENMPSCHEDULE(),OPENMPCHUNKSIZE(blocksize))
969 
970  rj=j;
971  rk=k;
972  LOOPBOUND1IN{
973  PBOUNDLOOP(pliter,pl){
974  // SECTIONMARK: assume r=0 singularity can't move
975  if(dirprim[pl]==FACE1 || dirprim[pl]==CORN3 || dirprim[pl]==CORN2 || dirprim[pl]==CORNT ) ri = -i; // FACE1 values
976  else ri=-i-1; // "CENT" values for purposes of this BC
977  MACP0A1(prim,i,j,k,pl) = MACP0A1(prim,ri,rj,rk,pl);
978  }// over pl
979  }// over boundary zones
980  }
981  }// end block
982 
983 
984 
985  if((BCtype[X1DN]==R0SING)||(BCtype[X1DN]==ASYMM) ){
986 
987  /* make sure b and u are antisymmetric at the poles (preserves u^t rho and u) */
989 
992 #pragma omp for schedule(OPENMPSCHEDULE(),OPENMPCHUNKSIZE(blocksize))
995 
996 
997  // SECTIONMARK: assume r=0 singularity can't move
998  i=0;
999  PBOUNDLOOP(pliter,pl){
1000  if(pl==U1 || pl==URAD1 || pl==B1){
1001  if(dirprim[pl]==FACE1 || dirprim[pl]==CORN3 || dirprim[pl]==CORN2 || dirprim[pl]==CORNT ){
1002  MACP0A1(prim,i,j,k,pl) = 0.0;
1003  }
1004  }// else don't do this pl
1005  } // end over pl
1006 
1007  LOOPBOUND1IN {
1008  PBOUNDLOOP(pliter,pl){
1009  if(pl==U1 || pl==URAD1 || pl==B1){
1010  MACP0A1(prim,i,j,k,pl) *= -1.;
1011  }// end if right pl
1012  } // end over pl
1013  } // end over boundary zones
1014  }// end loop 23
1015  }
1016  } //end if inner CPU wall
1017  }
1018  else{
1019  dualfprintf(fail_file,"Shouldn't be here in bounds\n");
1020  myexit(3946840);
1021  }
1022 
1023  } // end parallel region
1024 
1025  if(BCtype[X1DN]==R0SING){
1026  bound_x1dn_r0singfixinterior(boundstage,finalstep, boundtime, whichdir, boundvartype, dirprim,ispstag, prim,inboundloop,outboundloop,innormalloop,outnormalloop,inoutlohi,riin,riout,rjin,rjout,rkin,rkout,dosetbc,enerregion,localenerpos);
1027  }
1028 
1029 
1030 
1031  return(0);
1032 
1033 }
1034 
1035 
1036 
1037 
1039 int bound_x2dn_outflow_simple(
1040  int boundstage, int finalstep, SFTYPE boundtime, int whichdir, int boundvartype, int *dirprim, int ispstag, FTYPE (*prim)[NSTORE2][NSTORE3][NPR],
1041  int *inboundloop,
1042  int *outboundloop,
1043  int *innormalloop,
1044  int *outnormalloop,
1045  int (*inoutlohi)[NUMUPDOWN][NDIM],
1046  int riin, int riout, int rjin, int rjout, int rkin, int rkout,
1047  int *dosetbc,
1048  int enerregion,
1049  int *localenerpos
1050  )
1051 {
1052 
1053 
1054 
1055 #pragma omp parallel // assume don't require EOS
1056  {
1057  int i,j,k,pl,pliter;
1058  int localj1,globalj1;
1059  int localj2,globalj2;
1060  int rj1,rj2;
1061  struct of_geom geom[NPR],rgeom[NPR];
1062  FTYPE vcon[NDIM]; // coordinate basis vcon
1063 #if(WHICHVEL==VEL3)
1064  int failreturn;
1065 #endif
1066  int ri, rj, rk; // reference i,j,k
1067  FTYPE prescale[NPR];
1068  int jj,kk;
1069  struct of_geom geomdontuse[NPR];
1070  struct of_geom *ptrgeom[NPR];
1071  struct of_geom rgeomdontuse[NPR];
1072  struct of_geom *ptrrgeom[NPR];
1073 
1074 
1075  // assign memory
1076  PALLLOOP(pl){
1077  ptrgeom[pl]=&(geomdontuse[pl]);
1078  ptrrgeom[pl]=&(rgeomdontuse[pl]);
1079  }
1080 
1081 
1082  if(((BCtype[X2DN]==OUTFLOW || BCtype[X2DN]==HORIZONOUTFLOW))||(BCtype[X2DN]==FIXEDOUTFLOW)||(BCtype[X1DN]==FREEOUTFLOW)){
1083 
1084 
1085  if ( (totalsize[2]>1) && (mycpupos[2] <= 0)) {
1086  /* inner r boundary condition: u, just copy */
1087 
1090 #pragma omp for schedule(OPENMPSCHEDULE(),OPENMPCHUNKSIZE(blocksize))
1093 
1094 
1095 
1096  ri=i;
1097  rj=rjin;
1098  rk=k;
1099 
1100 
1101  PALLLOOP(pl) get_geometry(ri, rj, rk, dirprim[pl], ptrrgeom[pl]);
1102 
1103  LOOPBOUND2IN{ // bound entire region inside non-evolved portion of grid
1104  PBOUNDLOOP(pliter,pl){
1105  // dualfprintf(fail_file,"pl=%d i=%d j=%d k=%d : ri=%d rj=%d rk=%d\n",pl,i,j,k,ri,rj,rk);
1106  MACP0A1(prim,i,j,k,pl) = MACP0A1(prim,ri,rj,rk,pl);
1107  }
1108  }
1109 
1110 
1111 
1112 
1113  if(ispstag==0){
1114 
1115  if((BCtype[X2DN]==OUTFLOW)||(BCtype[X2DN]==FIXEDOUTFLOW)){
1116  // GODMARK: assume all velocities at same location when doing inflow check
1117  LOOPBOUND2IN{
1118 #if(WHICHVEL==VEL4)
1119  get_geometry(i, j, k, dirprim[U2], ptrgeom[U2]);
1120  inflow_check_4vel(2,MAC(prim,i,j,k),NULL,ptrgeom[U2], 0) ;
1121 #elif(WHICHVEL==VEL3)
1122  get_geometry(i, j, k, dirprim[U2], ptrgeom[U2]);
1123  inflow_check_3vel(2,MAC(prim,i,j,k),NULL,ptrgeom[U2], 0) ;
1124  // projection may not preserve u^t to be real and rho>rhoscal u>uuscal
1125 #if(JONCHECKS)
1126  if(jonchecks){
1127  //fixup1zone(0,MAC(prim,i,j,k),ptrgeom[U2],0);
1128  failreturn=check_pr(MAC(prim,i,j,k),MAC(prim,i,j,k),ptrgeom[U2],-3);
1129  if(failreturn){
1130  dualfprintf(fail_file,"Bad boundary zone, couldn't fix: i=%d j=%d k=%d\n",startpos[1]+i,startpos[2]+j,startpos[3]+k);
1131  if (fail(i,j,k,FAIL_BCFIX) >= 1) return (1);
1132  }
1133  }
1134 #endif
1135 #elif(WHICHVEL==VELREL4)
1136  get_geometry(i,j,k,dirprim[U2],ptrgeom[U2]) ;
1137  inflow_check_rel4vel(2,MAC(prim,i,j,k),NULL,ptrgeom[U2],0) ;
1138  if(limit_gamma(0,GAMMAMAX,GAMMAMAXRAD,MAC(prim,i,j,k),NULL,ptrgeom[U2],0)>=1)
1139  FAILSTATEMENT("bounds.c:bound_prim()", "limit_gamma()", 1);
1140 #endif
1141  }
1142  } // end if not allowing inflow
1143  }
1144 
1145 
1146  }// end 2 3
1147 
1148  }// end if mycpupos[2]==0
1149 
1150 
1151 
1152  }
1153  else{
1154  dualfprintf(fail_file,"Shouldn't be here in bounds\n");
1155  myexit(39377841);
1156  }
1157 
1158 
1159  }// end parallel region
1160 
1161 
1162  return(0);
1163 }
1164 
1166 int bound_x2up_outflow_simple(
1167  int boundstage, int finalstep, SFTYPE boundtime, int whichdir, int boundvartype, int *dirprim, int ispstag, FTYPE (*prim)[NSTORE2][NSTORE3][NPR],
1168  int *inboundloop,
1169  int *outboundloop,
1170  int *innormalloop,
1171  int *outnormalloop,
1172  int (*inoutlohi)[NUMUPDOWN][NDIM],
1173  int riin, int riout, int rjin, int rjout, int rkin, int rkout,
1174  int *dosetbc,
1175  int enerregion,
1176  int *localenerpos
1177  )
1178 {
1179 
1180 
1181 
1182 
1183 #pragma omp parallel // assume don't require EOS
1184  {
1185 
1186  int i,j,k,pl,pliter;
1187  int localj1,globalj1;
1188  int localj2,globalj2;
1189  int rj1,rj2;
1190  FTYPE vcon[NDIM]; // coordinate basis vcon
1191 #if(WHICHVEL==VEL3)
1192  int failreturn;
1193 #endif
1194  int ri, rj, rk; // reference i,j,k
1195  FTYPE prescale[NPR];
1196  int jj,kk;
1197  struct of_geom geomdontuse[NPR];
1198  struct of_geom *ptrgeom[NPR];
1199  struct of_geom rgeomdontuse[NPR];
1200  struct of_geom *ptrrgeom[NPR];
1201 
1202  // assign memory
1203  PALLLOOP(pl){
1204  ptrgeom[pl]=&(geomdontuse[pl]);
1205  ptrrgeom[pl]=&(rgeomdontuse[pl]);
1206  }
1207 
1208 
1209 
1210  if((BCtype[X2UP]==OUTFLOW)||(BCtype[X2UP]==FIXEDOUTFLOW)||(BCtype[X2UP]==FREEOUTFLOW)){
1211 
1212  // outer r BC:
1213  if ( (totalsize[2]>1) && (mycpupos[2] == ncpux2 - 1) ) {
1214 
1215 
1218 #pragma omp for schedule(OPENMPSCHEDULE(),OPENMPCHUNKSIZE(blocksize))
1221 
1222 
1223  ri=i;
1224  rj=rjout;
1225  rk=k;
1226 
1227  PALLLOOP(pl) get_geometry(ri, rj, rk, dirprim[pl], ptrrgeom[pl]);
1228 
1229  LOOPBOUND2OUT PBOUNDLOOP(pliter,pl) MACP0A1(prim,i,j,k,pl) = MACP0A1(prim,ri,rj,rk,pl);
1230 
1231 
1232 
1233 
1234 
1235  if(ispstag==0){
1236 
1237  if((BCtype[X2UP]==OUTFLOW)||(BCtype[X2UP]==FIXEDOUTFLOW)){
1238 
1239  LOOPBOUND2OUT{
1240 #if(WHICHVEL==VEL4)
1241  get_geometry(i, j, k, dirprim[U2], ptrgeom[U2]);
1242  inflow_check_4vel(2,MAC(prim,i,j,k),NULL,ptrgeom[U2],0) ;
1243 #elif(WHICHVEL==VEL3)
1244  get_geometry(i, j, k, dirprim[U2], ptrgeom[U2]);
1245  inflow_check_3vel(2,MAC(prim,i,j,k),NULL,ptrgeom[U2],0) ;
1246  // projection may not preserve u^t to be real and rho>rhoscal u>uuscal
1247 #if(JONCHECKS)
1248  if(jonchecks){
1249  //fixup1zone(0,MAC(prim,i,j,k),ptrgeom[U2],0);
1250  failreturn=check_pr(MAC(prim,i,j,k),MAC(prim,i,j,k),ptrgeom[U2],-3);
1251  if(failreturn){
1252  dualfprintf(fail_file,"Bad boundary zone, couldn't fix: i=%d j=%d k=%d\n",startpos[1]+i,startpos[2]+j,startpos[3]+k);
1253  if (fail(i,j,k,FAIL_BCFIX) >= 1) return (1);
1254  }
1255  }
1256 #endif
1257 #elif(WHICHVEL==VELREL4)
1258  get_geometry(i,j,k,dirprim[U2],ptrgeom[U2]) ;
1259  inflow_check_rel4vel(2,MAC(prim,i,j,k),NULL,ptrgeom[U2],0) ;
1260  if(limit_gamma(0,GAMMAMAX,GAMMAMAXRAD,MAC(prim,i,j,k),NULL,ptrgeom[U2], 0)>=1)
1261  FAILSTATEMENT("bounds.c:bound_prim()", "limit_gamma()", 2);
1262 #endif
1263  }// end over i
1264  }// end if not allowing inflow
1265 
1266  }
1267 
1268 
1269 
1270  }// end 2 3
1271  }// end if mycpu is correct
1272 
1273  }
1274  else{
1275  dualfprintf(fail_file,"Shouldn't be here in bounds\n");
1276  myexit(39324842);
1277  }
1278 
1279  }// end parallel region
1280 
1281 
1282  return(0);
1283 }
1284 
1285 
1286 
1287 
1288 
1289 
1290 
1291 
1292 
1293 
1294 
1295 
1296 
1297 
1298 
1299 
1306  int whichcall,
1307  int boundstage, int finalstep, SFTYPE boundtime, int whichdir, int boundvartype, int *dirprim, int ispstag, FTYPE (*prim)[NSTORE2][NSTORE3][NPR],
1308  int *inboundloop,
1309  int *outboundloop,
1310  int *innormalloop,
1311  int *outnormalloop,
1312  int (*inoutlohi)[NUMUPDOWN][NDIM],
1313  int riin, int riout, int rjin, int rjout, int rkin, int rkout,
1314  int *dosetbc,
1315  int enerregion,
1316  int *localenerpos
1317  )
1318 {
1319 
1320 
1321 
1322 #pragma omp parallel // assume don't require EOS
1323  {
1324  int i,j,k,pl,pliter;
1325  FTYPE vcon[NDIM]; // coordinate basis vcon
1326 #if(WHICHVEL==VEL3)
1327  int failreturn;
1328 #endif
1329  int ri, rj, rk; // reference i,j,k
1330  FTYPE prescale[NPR];
1331  int jj,kk;
1332 
1333 
1334 
1335  if(BCtype[X2DN]==POLARAXIS && special3dspc){
1336 
1337  /* inner polar BC (preserves u^t rho and u) */
1338  if ( (totalsize[2]>1) && (mycpupos[2] == 0) ) {
1339 
1340  if(whichcall==1 && ncpux3==1){
1341  // if ncpux3>1 then this is handled by MPI
1343 
1345 #pragma omp for schedule(OPENMPSCHEDULE(),OPENMPCHUNKSIZE(blocksize))
1348 
1349  ri=i;
1350  rk=(k+(int)N3/2)%N3;
1351  LOOPBOUND2IN{
1352  PBOUNDLOOP(pliter,pl){
1353  // assume polar axis can't move SECTIONMARK
1354  if(dirprim[pl]==FACE2 || dirprim[pl]==CORN3 || dirprim[pl]==CORN1 || dirprim[pl]==CORNT ) rj = -j; // FACE2 values
1355  else rj=-j-1; // "CENT" values for purposes of this BC
1356  MACP0A1(prim,i,j,k,pl) = MACP0A1(prim,ri,rj,rk,pl);
1357 
1358  // flip sign
1359  if(pl==U1) MACP0A1(prim,i,j,k,pl) *= SIGNFLIPU1;
1360  if(pl==B1) MACP0A1(prim,i,j,k,pl) *= SIGNFLIPB1;
1361  if(pl==URAD1) MACP0A1(prim,i,j,k,pl) *= SIGNFLIPU1;
1362  if(pl==U2) MACP0A1(prim,i,j,k,pl) *= SIGNFLIPU2;
1363  if(pl==B2) MACP0A1(prim,i,j,k,pl) *= SIGNFLIPB2;
1364  if(pl==URAD2) MACP0A1(prim,i,j,k,pl) *= SIGNFLIPU2;
1365  if(pl==U3) MACP0A1(prim,i,j,k,pl) *= SIGNFLIPU3;
1366  if(pl==B3) MACP0A1(prim,i,j,k,pl) *= SIGNFLIPB3;
1367  if(pl==URAD3) MACP0A1(prim,i,j,k,pl) *= SIGNFLIPU3;
1368 
1369 #if(DEBUGINOUTLOOPS)
1370  dualfprintf(fail_file,"INNER pole1: ispstag=%d pl=%d :: ri=%d rj=%d rk=%d i=%d j=%d k=%d\n",ispstag,pl,ri,rj,rk,i,j,k);
1371  if(!isfinite(MACP0A1(prim,ri,rj,rk,pl))){
1372  dualfprintf(fail_file,"INNER pole1: caught copying nan ri=%d rj=%d rk=%d pl=%d\n",ri,rj,rk,pl);
1373  }
1374 #endif
1375 
1376 
1377  }// end over pl
1378  }// end over j
1379 
1380  // also do j=0 (this just makes B2 @ FACE2-type location at j=0 at k and rk the same in correct sense)
1381  j=0;
1382  PBOUNDLOOP(pliter,pl){
1383  if(dirprim[pl]==FACE2 || dirprim[pl]==CORN3 || dirprim[pl]==CORN1 || dirprim[pl]==CORNT ){
1384  // only do j=0 if at pole exactly
1385  rj = -j; // FACE2 values
1386  MACP0A1(prim,i,j,k,pl) = MACP0A1(prim,ri,rj,rk,pl);
1387 
1388  if(pl==U1) MACP0A1(prim,i,j,k,pl) *= SIGNFLIPU1;
1389  if(pl==B1) MACP0A1(prim,i,j,k,pl) *= SIGNFLIPB1;
1390  if(pl==URAD1) MACP0A1(prim,i,j,k,pl) *= SIGNFLIPU1;
1391  // flip sign
1392  if(pl==U2) MACP0A1(prim,i,j,k,pl) *= SIGNFLIPU2;
1393  if(pl==B2) MACP0A1(prim,i,j,k,pl) *= SIGNFLIPB2;
1394  if(pl==URAD2) MACP0A1(prim,i,j,k,pl) *= SIGNFLIPU2;
1395  // NOTEMARK: No U3,B3 staggered yet, so below not used
1396  if(pl==U3) MACP0A1(prim,i,j,k,pl) *= SIGNFLIPU3;
1397  if(pl==B3) MACP0A1(prim,i,j,k,pl) *= SIGNFLIPB3;
1398  if(pl==URAD3) MACP0A1(prim,i,j,k,pl) *= SIGNFLIPU3;
1399 
1400 #if(DEBUGINOUTLOOPS)
1401  dualfprintf(fail_file,"INNER pole2: ispstag=%d pl=%d :: ri=%d rj=%d rk=%d i=%d j=%d k=%d\n",ispstag,pl,ri,rj,rk,i,j,k);
1402  if(!isfinite(MACP0A1(prim,ri,rj,rk,pl))){
1403  dualfprintf(fail_file,"INNER pole2: caught copying nan ri=%d rj=%d rk=%d pl=%d\n",ri,rj,rk,pl);
1404  }
1405 #endif
1406 
1407 
1408  }
1409  }
1410 
1411  }// end over i,k
1412  }// end if ncpux3==1
1413 
1414  // SUPERGODMARK: continue to use for now
1415  // only do poledeath() after MPI call (i.e. whichcall==2)
1416  if(BCtype[X2DN]==POLARAXIS && (whichcall==2 && ncpux3>1 || whichcall==1 && ncpux3==1) ){
1417  if(POLEDEATH || POLEGAMMADEATH) poledeath(X2DN,boundstage,finalstep,boundtime,whichdir,boundvartype,dirprim,ispstag,prim,inboundloop,outboundloop,innormalloop,outnormalloop,inoutlohi,riin,riout,rjin,rjout,rkin,rkout,dosetbc,enerregion,localenerpos);
1418  if(POLESMOOTH) polesmooth(X2DN,boundstage,finalstep,boundtime,whichdir,boundvartype,dirprim,ispstag,prim,inboundloop,outboundloop,innormalloop,outnormalloop,inoutlohi,riin,riout,rjin,rjout,rkin,rkout,dosetbc,enerregion,localenerpos);
1419  }
1420 
1421 
1422  }//end if inner CPU wall
1423 
1424 
1425 
1426 
1427 
1428  }
1429  else{
1430  dualfprintf(fail_file,"Shouldn't be here in bounds\n");
1431  myexit(3946843);
1432  }
1433 
1434 
1435 
1436  } // end parallel region
1437 
1438 
1439 
1440 
1441  //end if special transmissive BCs at poles
1442  return(0);
1443 }
1444 
1445 
1446 
1447 
1448 
1451  int boundstage, int finalstep, SFTYPE boundtime, int whichdir, int boundvartype, int *dirprim, int ispstag, FTYPE (*prim)[NSTORE2][NSTORE3][NPR],
1452  int *inboundloop,
1453  int *outboundloop,
1454  int *innormalloop,
1455  int *outnormalloop,
1456  int (*inoutlohi)[NUMUPDOWN][NDIM],
1457  int riin, int riout, int rjin, int rjout, int rkin, int rkout,
1458  int *dosetbc,
1459  int enerregion,
1460  int *localenerpos
1461  )
1462 {
1463 
1464 
1465 
1466 #pragma omp parallel // assume don't require EOS
1467  {
1468  int i,j,k,pl,pliter;
1469  FTYPE vcon[NDIM]; // coordinate basis vcon
1470 #if(WHICHVEL==VEL3)
1471  int failreturn;
1472 #endif
1473  int ri, rj, rk; // reference i,j,k
1474  FTYPE prescale[NPR];
1475  int jj,kk;
1476 
1477 
1478  if((BCtype[X2DN]==POLARAXIS)||(BCtype[X2DN]==SYMM)||(BCtype[X2DN]==ASYMM) ){
1479 
1480 
1481  if ( (totalsize[2]>1) && (mycpupos[2] == 0) ){
1482 
1484 
1485  {// block region
1487 #pragma omp for schedule(OPENMPSCHEDULE(),OPENMPCHUNKSIZE(blocksize))
1490 
1491  ri=i;
1492  rk=k;
1493  LOOPBOUND2IN{
1494  PBOUNDLOOP(pliter,pl){
1495  // assume polar axis can't move SECTIONMARK
1496  if(dirprim[pl]==FACE2 || dirprim[pl]==CORN3 || dirprim[pl]==CORN1 || dirprim[pl]==CORNT ) rj = -j; // FACE2 values
1497  else rj=-j-1; // "CENT" values for purposes of this BC
1498  MACP0A1(prim,i,j,k,pl) = MACP0A1(prim,ri,rj,rk,pl);
1499  }
1500  }
1501  }
1502  }// end blocked region
1503 
1504 
1505  if((BCtype[X2DN]==POLARAXIS)||(BCtype[X2DN]==ASYMM) ){
1506 
1507  /* make sure b and u are antisymmetric at the poles (preserves u^t rho and u) */
1509 
1511 #pragma omp for schedule(OPENMPSCHEDULE(),OPENMPCHUNKSIZE(blocksize))
1514 
1515 
1516  // assume polar axis can't move SECTIONMARK
1517  j=0;
1518  PBOUNDLOOP(pliter,pl){
1519  if(pl==U2 || pl==URAD2 || pl==B2){
1520  if(dirprim[pl]==FACE2 || dirprim[pl]==CORN3 || dirprim[pl]==CORN1 || dirprim[pl]==CORNT ){
1521  MACP0A1(prim,i,j,k,pl) = 0.0;
1522  }
1523  }// else don't do this pl
1524  } // end over pl
1525 
1526 
1527  if(BCtype[X2DN]==POLARAXIS){
1528  LOOPBOUND2IN {
1529  PBOUNDLOOP(pliter,pl){
1530  if(pl==U1) MACP0A1(prim,i,j,k,pl) *= SIGNFLIPU1;
1531  if(pl==B1) MACP0A1(prim,i,j,k,pl) *= SIGNFLIPB1;
1532  if(pl==URAD1) MACP0A1(prim,i,j,k,pl) *= SIGNFLIPU1;
1533  if(pl==U2) MACP0A1(prim,i,j,k,pl) *= SIGNFLIPU2;
1534  if(pl==B2) MACP0A1(prim,i,j,k,pl) *= SIGNFLIPB2;
1535  if(pl==URAD2) MACP0A1(prim,i,j,k,pl) *= SIGNFLIPU2;
1536  if(pl==U3) MACP0A1(prim,i,j,k,pl) *= SIGNFLIPU3;
1537  if(pl==B3) MACP0A1(prim,i,j,k,pl) *= SIGNFLIPB3;
1538  if(pl==URAD3) MACP0A1(prim,i,j,k,pl) *= SIGNFLIPU3;
1539  } // end over pl
1540  } // end over boundary cells
1541  }// end of polar axis type
1542  else if(BCtype[X2DN]==ASYMM){
1543  LOOPBOUND2IN {
1544  PBOUNDLOOP(pliter,pl){
1545  if(pl==U2 || pl==URAD2 || pl==B2) MACP0A1(prim,i,j,k,pl) *= -1.0;
1546  } // end over pl
1547  } // end over boundary cells
1548  }// end of polar axis type
1549  }// end loop 13
1550 
1551 
1552  }// end if polar or asym condition
1553 
1554  if(BCtype[X2DN]==POLARAXIS){
1555  if(POLEDEATH || POLEGAMMADEATH) poledeath(X2DN,boundstage,finalstep,boundtime,whichdir,boundvartype,dirprim,ispstag,prim,inboundloop,outboundloop,innormalloop,outnormalloop,inoutlohi,riin,riout,rjin,rjout,rkin,rkout,dosetbc,enerregion,localenerpos);
1556  if(POLESMOOTH) polesmooth(X2DN,boundstage,finalstep,boundtime,whichdir,boundvartype,dirprim,ispstag,prim,inboundloop,outboundloop,innormalloop,outnormalloop,inoutlohi,riin,riout,rjin,rjout,rkin,rkout,dosetbc,enerregion,localenerpos);
1557  }
1558 
1559  }// end if inner CPU wall
1560 
1561 
1562 
1563  }// end if polar asym or asym
1564  else{
1565  dualfprintf(fail_file,"Shouldn't be here in bounds\n");
1566  myexit(3946844);
1567  }
1568 
1569  }// end parallel region
1570 
1571 
1572 
1573  return(0);
1574 }
1575 
1576 
1577 
1578 
1579 
1582  int whichcall,
1583  int boundstage, int finalstep, SFTYPE boundtime, int whichdir, int boundvartype, int *dirprim, int ispstag, FTYPE (*prim)[NSTORE2][NSTORE3][NPR],
1584  int *inboundloop,
1585  int *outboundloop,
1586  int *innormalloop,
1587  int *outnormalloop,
1588  int (*inoutlohi)[NUMUPDOWN][NDIM],
1589  int riin, int riout, int rjin, int rjout, int rkin, int rkout,
1590  int *dosetbc,
1591  int enerregion,
1592  int *localenerpos
1593  )
1594 {
1595 
1596 
1597 
1598 
1599 
1600 #pragma omp parallel // assume don't require EOSs
1601  {
1602  int i,j,k,pl,pliter;
1603  FTYPE vcon[NDIM]; // coordinate basis vcon
1604 #if(WHICHVEL==VEL3)
1605  int failreturn;
1606 #endif
1607  int ri, rj, rk; // reference i,j,k
1608  FTYPE prescale[NPR];
1609  int jj,kk;
1610 
1611 
1612 
1613  if(BCtype[X2UP]==POLARAXIS && special3dspc){
1614 
1615 
1616  if ( (totalsize[2]>1) && (mycpupos[2] == ncpux2-1) ) {
1617 
1618 
1619  if(whichcall==1 && ncpux3==1){
1620  // if ncpux3>1 then this is handled by MPI
1621 
1624 #pragma omp for schedule(OPENMPSCHEDULE(),OPENMPCHUNKSIZE(blocksize))
1627 
1628 
1629  ri=i;
1630  rk=(k+(int)N3/2)%N3;
1631  LOOPBOUND2OUT{
1632  PBOUNDLOOP(pliter,pl){
1633  // assume polar axis can't move SECTIONMARK
1634  if(dirprim[pl]==FACE2 || dirprim[pl]==CORN3 || dirprim[pl]==CORN1 || dirprim[pl]==CORNT ) rj = 2*N2-j; // FACE2 values
1635  else rj=2*(N2-1)-j+1; // "CENT" values for purposes of this BC
1636  MACP0A1(prim,i,j,k,pl) = MACP0A1(prim,ri,rj,rk,pl);
1637 
1638  // flip sign
1639  if(pl==U1) MACP0A1(prim,i,j,k,pl) *= SIGNFLIPU1;
1640  if(pl==B1) MACP0A1(prim,i,j,k,pl) *= SIGNFLIPB1;
1641  if(pl==URAD1) MACP0A1(prim,i,j,k,pl) *= SIGNFLIPU1;
1642  if(pl==U2) MACP0A1(prim,i,j,k,pl) *= SIGNFLIPU2;
1643  if(pl==B2) MACP0A1(prim,i,j,k,pl) *= SIGNFLIPB2;
1644  if(pl==URAD2) MACP0A1(prim,i,j,k,pl) *= SIGNFLIPU2;
1645  if(pl==U3) MACP0A1(prim,i,j,k,pl) *= SIGNFLIPU3;
1646  if(pl==B3) MACP0A1(prim,i,j,k,pl) *= SIGNFLIPB3;
1647  if(pl==URAD3) MACP0A1(prim,i,j,k,pl) *= SIGNFLIPU3;
1648 
1649 #if(DEBUGINOUTLOOPS)
1650  dualfprintf(fail_file,"OUTER pole1: ispstag=%d pl=%d :: ri=%d rj=%d rk=%d i=%d j=%d k=%d\n",ispstag,pl,ri,rj,rk,i,j,k);
1651  if(!isfinite(MACP0A1(prim,ri,rj,rk,pl))){
1652  dualfprintf(fail_file,"OUTER pole1: caught copying nan ri=%d rj=%d rk=%d pl=%d\n",ri,rj,rk,pl);
1653  }
1654 #endif
1655 
1656 
1657  }// end over pl
1658  }// end over j
1659 
1660  // also do j=N2 (this just makes B2 @ FACE2-type location at j=N2 at k and rk the same in correct sense)
1661  j=N2;
1662  PBOUNDLOOP(pliter,pl){
1663  if(dirprim[pl]==FACE2 || dirprim[pl]==CORN3 || dirprim[pl]==CORN1 || dirprim[pl]==CORNT ){
1664  // only do j=0 if at pole exactly
1665  rj = 2*N2-j; // FACE2 values
1666  MACP0A1(prim,i,j,k,pl) = MACP0A1(prim,ri,rj,rk,pl);
1667 
1668  // flip sign
1669  if(pl==U1) MACP0A1(prim,i,j,k,pl) *= SIGNFLIPU1;
1670  if(pl==B1) MACP0A1(prim,i,j,k,pl) *= SIGNFLIPB1;
1671  if(pl==URAD1) MACP0A1(prim,i,j,k,pl) *= SIGNFLIPU1;
1672  if(pl==U2) MACP0A1(prim,i,j,k,pl) *= SIGNFLIPU2;
1673  if(pl==B2) MACP0A1(prim,i,j,k,pl) *= SIGNFLIPB2;
1674  if(pl==URAD2) MACP0A1(prim,i,j,k,pl) *= SIGNFLIPU2;
1675  if(pl==U3) MACP0A1(prim,i,j,k,pl) *= SIGNFLIPU3;
1676  if(pl==B3) MACP0A1(prim,i,j,k,pl) *= SIGNFLIPB3;
1677  if(pl==URAD3) MACP0A1(prim,i,j,k,pl) *= SIGNFLIPU3;
1678 
1679 #if(DEBUGINOUTLOOPS)
1680  dualfprintf(fail_file,"OUTER pole2: ispstag=%d pl=%d :: ri=%d rj=%d rk=%d i=%d j=%d k=%d\n",ispstag,pl,ri,rj,rk,i,j,k);
1681  if(!isfinite(MACP0A1(prim,ri,rj,rk,pl))){
1682  dualfprintf(fail_file,"OUTER pole2: caught copying nan ri=%d rj=%d rk=%d pl=%d\n",ri,rj,rk,pl);
1683  }
1684 #endif
1685 
1686  }
1687  }
1688 
1689  }// end over i,k
1690  }// end if ncpux3==1
1691 
1692  // SUPERGODMARK: continue to use for now
1693  // only do poledeath() after MPI call (i.e. whichcall==2)
1694  if(BCtype[X2UP]==POLARAXIS && (whichcall==2 && ncpux3>1 || whichcall==1 && ncpux3==1) ){
1695  if(POLEDEATH || POLEGAMMADEATH) poledeath(X2UP,boundstage,finalstep,boundtime,whichdir,boundvartype,dirprim,ispstag,prim,inboundloop,outboundloop,innormalloop,outnormalloop,inoutlohi,riin,riout,rjin,rjout,rkin,rkout,dosetbc,enerregion,localenerpos);
1696  if(POLESMOOTH) polesmooth(X2UP,boundstage,finalstep,boundtime,whichdir,boundvartype,dirprim,ispstag,prim,inboundloop,outboundloop,innormalloop,outnormalloop,inoutlohi,riin,riout,rjin,rjout,rkin,rkout,dosetbc,enerregion,localenerpos);
1697  }
1698 
1699  }// end if outer CPU wall
1700 
1701  }
1702  else{
1703  dualfprintf(fail_file,"Shouldn't be here in bounds\n");
1704  myexit(3946845);
1705  }
1706 
1707  }// end parallel region
1708 
1709 
1710  return(0);
1711 }
1712 
1713 
1714 
1715 
1716 
1719  int boundstage, int finalstep, SFTYPE boundtime, int whichdir, int boundvartype, int *dirprim, int ispstag, FTYPE (*prim)[NSTORE2][NSTORE3][NPR],
1720  int *inboundloop,
1721  int *outboundloop,
1722  int *innormalloop,
1723  int *outnormalloop,
1724  int (*inoutlohi)[NUMUPDOWN][NDIM],
1725  int riin, int riout, int rjin, int rjout, int rkin, int rkout,
1726  int *dosetbc,
1727  int enerregion,
1728  int *localenerpos
1729  )
1730 {
1731 
1732 
1733 
1734 
1735 
1736 #pragma omp parallel // assume don't require EOS
1737  {
1738  int i,j,k,pl,pliter;
1739  FTYPE vcon[NDIM]; // coordinate basis vcon
1740 #if(WHICHVEL==VEL3)
1741  int failreturn;
1742 #endif
1743  int ri, rj, rk; // reference i,j,k
1744  FTYPE prescale[NPR];
1745  int jj,kk;
1746 
1747 
1748 
1749 
1750 
1751  if((BCtype[X2UP]==POLARAXIS)||(BCtype[X2UP]==SYMM)||(BCtype[X2UP]==ASYMM) ){
1752 
1753 
1754  if ( (totalsize[2]>1) && (mycpupos[2] == ncpux2-1) ) {
1755 
1757  {// block region
1759 #pragma omp for schedule(OPENMPSCHEDULE(),OPENMPCHUNKSIZE(blocksize))
1762 
1763 
1764  ri=i;
1765  rk=k;
1766  LOOPBOUND2OUT{
1767  PBOUNDLOOP(pliter,pl){
1768  // assume polar axis can't move SECTIONMARK
1769  if(dirprim[pl]==FACE2 || dirprim[pl]==CORN3 || dirprim[pl]==CORN1 || dirprim[pl]==CORNT ) rj = 2*N2-j; // FACE2 values
1770  else rj=2*(N2-1)-j+1; // "CENT" values for purposes of this BC
1771  MACP0A1(prim,i,j,k,pl) = MACP0A1(prim,ri,rj,rk,pl);
1772  }
1773  }
1774  }
1775  }
1776 
1777  if((BCtype[X2UP]==POLARAXIS)||(BCtype[X2UP]==ASYMM) ){
1778 
1779  /* make sure b and u are antisymmetric at the poles (preserves u^t rho and u) */
1782 #pragma omp for schedule(OPENMPSCHEDULE(),OPENMPCHUNKSIZE(blocksize))
1785 
1786 
1787  j=N2;
1788  PBOUNDLOOP(pliter,pl){
1789  if(pl==U2 || pl==URAD2 || pl==B2){
1790  if(dirprim[pl]==FACE2 || dirprim[pl]==CORN3 || dirprim[pl]==CORN1 || dirprim[pl]==CORNT ){
1791  MACP0A1(prim,i,j,k,pl) = 0.0;
1792  }
1793  }// else don't do this pl
1794  } // end over pl
1795 
1796 
1797  if(BCtype[X2UP]==POLARAXIS){
1798  LOOPBOUND2OUT {
1799  PBOUNDLOOP(pliter,pl){
1800  if(pl==U1) MACP0A1(prim,i,j,k,pl) *= SIGNFLIPU1;
1801  if(pl==B1) MACP0A1(prim,i,j,k,pl) *= SIGNFLIPB1;
1802  if(pl==URAD1) MACP0A1(prim,i,j,k,pl) *= SIGNFLIPU1;
1803  if(pl==U2) MACP0A1(prim,i,j,k,pl) *= SIGNFLIPU2;
1804  if(pl==B2) MACP0A1(prim,i,j,k,pl) *= SIGNFLIPB2;
1805  if(pl==URAD2) MACP0A1(prim,i,j,k,pl) *= SIGNFLIPU2;
1806  if(pl==U3) MACP0A1(prim,i,j,k,pl) *= SIGNFLIPU3;
1807  if(pl==B3) MACP0A1(prim,i,j,k,pl) *= SIGNFLIPB3;
1808  if(pl==URAD3) MACP0A1(prim,i,j,k,pl) *= SIGNFLIPU3;
1809  } // end over pl
1810  } // end over bondary cells
1811  }// end if polar axis
1812  else if(BCtype[X2UP]==ASYMM){
1813  LOOPBOUND2OUT {
1814  PBOUNDLOOP(pliter,pl){
1815  if(pl==U2 || pl==URAD2 || pl==B2) MACP0A1(prim,i,j,k,pl) *= -1.0;
1816  } // end over pl
1817  } // end over bondary cells
1818  }
1819  }// end loop 13
1820 
1821  }// end if reflecting type conditions
1822 
1823 
1824  if(BCtype[X2UP]==POLARAXIS){
1825  if(POLEDEATH || POLEGAMMADEATH) poledeath(X2UP,boundstage,finalstep,boundtime,whichdir,boundvartype,dirprim,ispstag,prim,inboundloop,outboundloop,innormalloop,outnormalloop,inoutlohi,riin,riout,rjin,rjout,rkin,rkout,dosetbc,enerregion,localenerpos);
1826  if(POLESMOOTH) polesmooth(X2UP,boundstage,finalstep,boundtime,whichdir,boundvartype,dirprim,ispstag,prim,inboundloop,outboundloop,innormalloop,outnormalloop,inoutlohi,riin,riout,rjin,rjout,rkin,rkout,dosetbc,enerregion,localenerpos);
1827  }
1828 
1829  }// end if mycpupos[2]==ncpux2-1
1830 
1831 
1832 
1833  }
1834  else{
1835  dualfprintf(fail_file,"Shouldn't be here in bounds\n");
1836  myexit(3946846);
1837  }
1838 
1839  }// end parallel region
1840 
1841 
1842  return(0);
1843 }
1844 
1845 
1846 
1847 
1850  int boundstage, int finalstep, SFTYPE boundtime, int whichdir, int boundvartype, int *dirprim, int ispstag, FTYPE (*prim)[NSTORE2][NSTORE3][NPR],
1851  int *inboundloop,
1852  int *outboundloop,
1853  int *innormalloop,
1854  int *outnormalloop,
1855  int (*inoutlohi)[NUMUPDOWN][NDIM],
1856  int riin, int riout, int rjin, int rjout, int rkin, int rkout,
1857  int *dosetbc,
1858  int enerregion,
1859  int *localenerpos
1860  )
1861 {
1862 
1863 
1864 
1865 #pragma omp parallel // assume don't require EOS
1866  {
1867  int i,j,k,pl,pliter;
1868  FTYPE vcon[NDIM]; // coordinate basis vcon
1869 #if(WHICHVEL==VEL3)
1870  int failreturn;
1871 #endif
1872  int ri, rj, rk; // reference i,j,k
1873  FTYPE prescale[NPR];
1874  int jj,kk;
1875 
1876 
1877  if( (BCtype[X1DN]==PERIODIC)&&(BCtype[X1UP]==PERIODIC) ){
1878 
1879 
1880  // periodic x1
1881  if ( (totalsize[1]>1) && (mycpupos[1] == 0)&&(mycpupos[1] == ncpux1 - 1) ) {
1882  // just copy from one side to another
1883 
1885 
1887 #pragma omp for schedule(OPENMPSCHEDULE(),OPENMPCHUNKSIZE(blocksize))
1890 
1891 
1892  // copy from upper side to lower boundary zones
1893  ri=riout;
1894  rj=j;
1895  rk=k;
1896  LOOPBOUND1IN PBOUNDLOOP(pliter,pl){
1897  MACP0A1(prim,i,j,k,pl) = MACP0A1(prim,ri+1+i,rj,rk,pl);
1898 
1899 
1900 #if(DEBUGINOUTLOOPS)
1901  dualfprintf(fail_file,"INNER X1: ispstag=%d pl=%d :: ri=%d rj=%d rk=%d (otherri=%d) i=%d j=%d k=%d\n",ispstag,pl,ri,rj,rk,ri+SHIFT1+i,i,j,k);
1902  if(!isfinite(MACP0A1(prim,ri+SHIFT1+i,rj,rk,pl))){
1903  dualfprintf(fail_file,"INNER X1: caught copying nan ri=%d rj=%d rk=%d pl=%d\n",ri,rj,rk,pl);
1904  }
1905 #endif
1906 
1907  }
1908 
1909  // copy from lower side to upper boundary zones
1910  ri=riin;
1911  rj=j;
1912  rk=k;
1913  LOOPBOUND1OUT PBOUNDLOOP(pliter,pl){
1914  MACP0A1(prim,i,j,k,pl) = MACP0A1(prim,ri+(i-N1),rj,rk,pl);
1915 
1916 #if(DEBUGINOUTLOOPS)
1917  dualfprintf(fail_file,"OUTER X1: ispstag=%d pl=%d :: ri=%d rj=%d rk=%d (otherri=%d) i=%d j=%d k=%d\n",ispstag,pl,ri,rj,ri+(i-N1),rk,i,j,k);
1918  if(!isfinite(MACP0A1(prim,ri+(i-N1),rj,rk,pl))){
1919  dualfprintf(fail_file,"INNER X1: caught copying nan ri=%d rj=%d rk=%d pl=%d\n",ri,rj,rk,pl);
1920  }
1921 #endif
1922 
1923  }
1924  }
1925  }
1926 
1927 
1928  }
1929  else{
1930  dualfprintf(fail_file,"Shouldn't be here in bounds\n");
1931  myexit(1946847);
1932  }
1933 
1934 
1935  }// end parallel region
1936 
1937  return(0);
1938 }
1939 
1942  int boundstage, int finalstep, SFTYPE boundtime, int whichdir, int boundvartype, int *dirprim, int ispstag, FTYPE (*prim)[NSTORE2][NSTORE3][NPR],
1943  int *inboundloop,
1944  int *outboundloop,
1945  int *innormalloop,
1946  int *outnormalloop,
1947  int (*inoutlohi)[NUMUPDOWN][NDIM],
1948  int riin, int riout, int rjin, int rjout, int rkin, int rkout,
1949  int *dosetbc,
1950  int enerregion,
1951  int *localenerpos
1952  )
1953 {
1954 
1955 
1956 
1957 #pragma omp parallel // assume don't require EOS
1958  {
1959  int i,j,k,pl,pliter;
1960  FTYPE vcon[NDIM]; // coordinate basis vcon
1961 #if(WHICHVEL==VEL3)
1962  int failreturn;
1963 #endif
1964  int ri, rj, rk; // reference i,j,k
1965  FTYPE prescale[NPR];
1966  int jj,kk;
1967 
1968 
1969  if( (BCtype[X2DN]==PERIODIC)&&(BCtype[X2UP]==PERIODIC) ){
1970 
1971 
1972  // periodic x2
1973  if ( (totalsize[2]>1) && (mycpupos[2] == 0)&&(mycpupos[2] == ncpux2 - 1) ) {
1974  // just copy from one side to another
1975 
1977 
1979 #pragma omp for schedule(OPENMPSCHEDULE(),OPENMPCHUNKSIZE(blocksize))
1982 
1983 
1984  // copy from upper side to lower boundary zones
1985  ri=i;
1986  rj=rjout;
1987  rk=k;
1988  LOOPBOUND2IN PBOUNDLOOP(pliter,pl){
1989  MACP0A1(prim,i,j,k,pl) = MACP0A1(prim,ri,rj+SHIFT2+j,rk,pl);
1990 
1991 
1992 #if(DEBUGINOUTLOOPS)
1993  dualfprintf(fail_file,"INNER X2: ispstag=%d pl=%d :: ri=%d rj=%d rk=%d (otherrj=%d) i=%d j=%d k=%d\n",ispstag,pl,ri,rj,rk,rj+SHIFT2+j,i,j,k);
1994  dualfprintf(fail_file,"ref: %g\n",MACP0A1(prim,ri,rj+SHIFT2+j,rk,pl));
1995  if(!isfinite(MACP0A1(prim,ri,rj+SHIFT2+j,rk,pl))){
1996  dualfprintf(fail_file,"INNER X2: caught copying nan ri=%d rj=%d rk=%d pl=%d\n",ri,rj,rk,pl);
1997  myexit(29533246346);
1998  }
1999 #endif
2000 
2001  }
2002 
2003  // copy from lower side to upper boundary zones
2004  ri=i;
2005  rj=rjin;
2006  rk=k;
2007  LOOPBOUND2OUT PBOUNDLOOP(pliter,pl){
2008  MACP0A1(prim,i,j,k,pl) = MACP0A1(prim,ri,rj+(j-N2),rk,pl);
2009 
2010 #if(DEBUGINOUTLOOPS)
2011  dualfprintf(fail_file,"OUTER X2: ispstag=%d pl=%d :: ri=%d rj=%d rk=%d (otherrj=%d) i=%d j=%d k=%d\n",ispstag,pl,ri,rj,rk,rj+(j-N2),i,j,k);
2012  dualfprintf(fail_file,"ref: %g\n",MACP0A1(prim,ri,rj+(j-N2),rk,pl));
2013  if(!isfinite(MACP0A1(prim,ri,rj+(j-N2),rk,pl))){
2014  dualfprintf(fail_file,"INNER X2: caught copying nan ri=%d rj=%d rk=%d pl=%d\n",ri,rj,rk,pl);
2015  myexit(3294863498634);
2016  }
2017 #endif
2018 
2019  }
2020  }
2021  }
2022 
2023 
2024  }
2025  else{
2026  dualfprintf(fail_file,"Shouldn't be here in bounds\n");
2027  myexit(2946848);
2028  }
2029 
2030 
2031  }// end parallel region
2032 
2033  return(0);
2034 }
2035 
2036 
2037 
2038 
2039 
2040 
2042 int bound_x3dn_outflow_simple(
2043  int boundstage, int finalstep, SFTYPE boundtime, int whichdir, int boundvartype, int *dirprim, int ispstag, FTYPE (*prim)[NSTORE2][NSTORE3][NPR],
2044  int *inboundloop,
2045  int *outboundloop,
2046  int *innormalloop,
2047  int *outnormalloop,
2048  int (*inoutlohi)[NUMUPDOWN][NDIM],
2049  int riin, int riout, int rjin, int rjout, int rkin, int rkout,
2050  int *dosetbc,
2051  int enerregion,
2052  int *localenerpos
2053  )
2054 {
2055 
2056 
2057 
2058 #pragma omp parallel // assume don't require EOS
2059  {
2060  int i,j,k,pl,pliter;
2061  int localj1,globalj1;
2062  int localk2,globalk2;
2063  int rk1,rk2;
2064  struct of_geom geom[NPR],rgeom[NPR];
2065  FTYPE vcon[NDIM]; // coordinate basis vcon
2066 #if(WHICHVEL==VEL3)
2067  int failreturn;
2068 #endif
2069  int ri, rj, rk; // reference i,j,k
2070  FTYPE prescale[NPR];
2071  int jj,kk;
2072  struct of_geom geomdontuse[NPR];
2073  struct of_geom *ptrgeom[NPR];
2074  struct of_geom rgeomdontuse[NPR];
2075  struct of_geom *ptrrgeom[NPR];
2076 
2077 
2078  // assign memory
2079  PALLLOOP(pl){
2080  ptrgeom[pl]=&(geomdontuse[pl]);
2081  ptrrgeom[pl]=&(rgeomdontuse[pl]);
2082  }
2083 
2084 
2085  if((BCtype[X3DN]==OUTFLOW)||(BCtype[X3DN]==FIXEDOUTFLOW)||(BCtype[X1DN]==FREEOUTFLOW)){
2086 
2087 
2088  if ( (totalsize[3]>1) && (mycpupos[3] <= 0)) {
2089  /* inner r boundary condition: u, just copy */
2090 
2093 #pragma omp for schedule(OPENMPSCHEDULE(),OPENMPCHUNKSIZE(blocksize))
2096 
2097 
2098 
2099  ri=i;
2100  rj=j;
2101  rk=rkin;
2102 
2103 
2104  PALLLOOP(pl) get_geometry(ri, rj, rk, dirprim[pl], ptrrgeom[pl]);
2105 
2106  LOOPBOUND3IN{ // bound entire region inside non-evolved portion of grid
2107  PBOUNDLOOP(pliter,pl) MACP0A1(prim,i,j,k,pl) = MACP0A1(prim,ri,rj,rk,pl);
2108  }
2109 
2110 
2111 
2112 
2113  if(ispstag==0){
2114 
2115  if((BCtype[X3DN]==OUTFLOW)||(BCtype[X3DN]==FIXEDOUTFLOW)){
2116  // GODMARK: assume all velocities at same location when doing inflow check
2117  LOOPBOUND3IN{
2118 #if(WHICHVEL==VEL4)
2119  get_geometry(i, j, k, dirprim[U3], ptrgeom[U3]);
2120  inflow_check_4vel(3,MAC(prim,i,j,k),NULL,ptrgeom[U3], 0) ;
2121 #elif(WHICHVEL==VEL3)
2122  get_geometry(i, j, k, dirprim[U3], ptrgeom[U3]);
2123  inflow_check_3vel(3,MAC(prim,i,j,k),NULL,ptrgeom[U3], 0) ;
2124  // projection may not preserve u^t to be real and rho>rhoscal u>uuscal
2125 #if(JONCHECKS)
2126  if(jonchecks){
2127  //fixup1zone(0,MAC(prim,i,j,k),ptrgeom[U3],0);
2128  failreturn=check_pr(MAC(prim,i,j,k),MAC(prim,i,j,k),ptrgeom[U3],-3);
2129  if(failreturn){
2130  dualfprintf(fail_file,"Bad boundary zone, couldn't fix: i=%d j=%d k=%d\n",startpos[1]+i,startpos[2]+j,startpos[3]+k);
2131  if (fail(i,j,k,FAIL_BCFIX) >= 1) return (1);
2132  }
2133  }
2134 #endif
2135 #elif(WHICHVEL==VELREL4)
2136  get_geometry(i,j,k,dirprim[U3],ptrgeom[U3]) ;
2137  inflow_check_rel4vel(3,MAC(prim,i,j,k),NULL,ptrgeom[U3],0) ;
2138  if(limit_gamma(0,GAMMAMAX,GAMMAMAXRAD,MAC(prim,i,j,k),NULL,ptrgeom[U3],0)>=1)
2139  FAILSTATEMENT("bounds.c:bound_prim()", "limit_gamma()", 1);
2140 #endif
2141  }
2142  } // end if not allowing inflow
2143  }
2144 
2145 
2146  }// end 1 2
2147 
2148  }// end if mycpupos[2]==0
2149 
2150 
2151 
2152  }
2153  else{
2154  dualfprintf(fail_file,"Shouldn't be here in bounds\n");
2155  myexit(312475849);
2156  }
2157 
2158 
2159  }// end parallel region
2160 
2161 
2162  return(0);
2163 }
2164 
2166 int bound_x3up_outflow_simple(
2167  int boundstage, int finalstep, SFTYPE boundtime, int whichdir, int boundvartype, int *dirprim, int ispstag, FTYPE (*prim)[NSTORE2][NSTORE3][NPR],
2168  int *inboundloop,
2169  int *outboundloop,
2170  int *innormalloop,
2171  int *outnormalloop,
2172  int (*inoutlohi)[NUMUPDOWN][NDIM],
2173  int riin, int riout, int rjin, int rjout, int rkin, int rkout,
2174  int *dosetbc,
2175  int enerregion,
2176  int *localenerpos
2177  )
2178 {
2179 
2180 
2181 
2182 
2183 #pragma omp parallel // assume don't require EOS
2184  {
2185 
2186  int i,j,k,pl,pliter;
2187  int localk1,globalk1;
2188  int localk2,globalk2;
2189  int rk1,rk2;
2190  FTYPE vcon[NDIM]; // coordinate basis vcon
2191 #if(WHICHVEL==VEL3)
2192  int failreturn;
2193 #endif
2194  int ri, rj, rk; // reference i,j,k
2195  FTYPE prescale[NPR];
2196  int jj,kk;
2197  struct of_geom geomdontuse[NPR];
2198  struct of_geom *ptrgeom[NPR];
2199  struct of_geom rgeomdontuse[NPR];
2200  struct of_geom *ptrrgeom[NPR];
2201 
2202  // assign memory
2203  PALLLOOP(pl){
2204  ptrgeom[pl]=&(geomdontuse[pl]);
2205  ptrrgeom[pl]=&(rgeomdontuse[pl]);
2206  }
2207 
2208 
2209 
2210  if((BCtype[X3UP]==OUTFLOW)||(BCtype[X3UP]==FIXEDOUTFLOW)||(BCtype[X3UP]==FREEOUTFLOW)){
2211 
2212  // outer r BC:
2213  if ( (totalsize[3]>1) && (mycpupos[3] == ncpux3 - 1) ) {
2214 
2215 
2218 #pragma omp for schedule(OPENMPSCHEDULE(),OPENMPCHUNKSIZE(blocksize))
2221 
2222 
2223  ri=i;
2224  rj=j;
2225  rk=rkout;
2226 
2227  PALLLOOP(pl) get_geometry(ri, rj, rk, dirprim[pl], ptrrgeom[pl]);
2228 
2229  LOOPBOUND3OUT PBOUNDLOOP(pliter,pl) MACP0A1(prim,i,j,k,pl) = MACP0A1(prim,ri,rj,rk,pl);
2230 
2231 
2232 
2233 
2234 
2235  if(ispstag==0){
2236 
2237  if((BCtype[X3UP]==OUTFLOW)||(BCtype[X3UP]==FIXEDOUTFLOW)){
2238 
2239  LOOPBOUND3OUT{
2240 #if(WHICHVEL==VEL4)
2241  get_geometry(i, j, k, dirprim[U3], ptrgeom[U3]);
2242  inflow_check_4vel(3,MAC(prim,i,j,k),NULL,ptrgeom[U3],0) ;
2243 #elif(WHICHVEL==VEL3)
2244  get_geometry(i, j, k, dirprim[U3], ptrgeom[U3]);
2245  inflow_check_3vel(3,MAC(prim,i,j,k),NULL,ptrgeom[U3],0) ;
2246  // projection may not preserve u^t to be real and rho>rhoscal u>uuscal
2247 #if(JONCHECKS)
2248  if(jonchecks){
2249  //fixup1zone(0,MAC(prim,i,j,k),ptrgeom[U3],0);
2250  failreturn=check_pr(MAC(prim,i,j,k),MAC(prim,i,j,k),ptrgeom[U3],-3);
2251  if(failreturn){
2252  dualfprintf(fail_file,"Bad boundary zone, couldn't fix: i=%d j=%d k=%d\n",startpos[1]+i,startpos[2]+j,startpos[3]+k);
2253  if (fail(i,j,k,FAIL_BCFIX) >= 1) return (1);
2254  }
2255  }
2256 #endif
2257 #elif(WHICHVEL==VELREL4)
2258  get_geometry(i,j,k,dirprim[U3],ptrgeom[U3]) ;
2259  inflow_check_rel4vel(3,MAC(prim,i,j,k),NULL,ptrgeom[U3],0) ;
2260  if(limit_gamma(0,GAMMAMAX,GAMMAMAXRAD,MAC(prim,i,j,k),NULL,ptrgeom[U3], 0)>=1)
2261  FAILSTATEMENT("bounds.c:bound_prim()", "limit_gamma()", 2);
2262 #endif
2263  }// end over i
2264  }// end if not allowing inflow
2265 
2266  }
2267 
2268 
2269 
2270  }// end 1 2
2271  }// end if mycpu is correct
2272 
2273  }
2274  else{
2275  dualfprintf(fail_file,"Shouldn't be here in bounds\n");
2276  myexit(315875850);
2277  }
2278 
2279  }// end parallel region
2280 
2281 
2282  return(0);
2283 }
2284 
2285 
2286 
2287 
2288 
2289 
2290 
2291 
2292 
2293 
2294 
2295 
2296 
2297 
2298 
2299 
2300 
2301 
2302 
2303 
2306  int boundstage, int finalstep, SFTYPE boundtime, int whichdir, int boundvartype, int *dirprim, int ispstag, FTYPE (*prim)[NSTORE2][NSTORE3][NPR],
2307  int *inboundloop,
2308  int *outboundloop,
2309  int *innormalloop,
2310  int *outnormalloop,
2311  int (*inoutlohi)[NUMUPDOWN][NDIM],
2312  int riin, int riout, int rjin, int rjout, int rkin, int rkout,
2313  int *dosetbc,
2314  int enerregion,
2315  int *localenerpos
2316  )
2317 {
2318 
2319 
2320 
2321 #pragma omp parallel // assume don't require EOS
2322  {
2323  int i,j,k,pl,pliter;
2324  FTYPE vcon[NDIM]; // coordinate basis vcon
2325 #if(WHICHVEL==VEL3)
2326  int failreturn;
2327 #endif
2328  int ri, rj, rk; // reference i,j,k
2329  FTYPE prescale[NPR];
2330  int jj,kk;
2331 
2332 
2333  if( (BCtype[X3DN]==PERIODIC)&&(BCtype[X3UP]==PERIODIC) ){
2334 
2335 
2336  // periodic x3
2337  if ( (totalsize[3]>1) && (mycpupos[3] == 0)&&(mycpupos[3] == ncpux3 - 1) ) {
2338  // just copy from one side to another
2339 
2341 
2343 #pragma omp for schedule(OPENMPSCHEDULE(),OPENMPCHUNKSIZE(blocksize))
2346 
2347 
2348  // copy from upper side to lower boundary zones
2349  ri=i;
2350  rj=j;
2351  rk=rkout;
2352  LOOPBOUND3IN PBOUNDLOOP(pliter,pl){
2353  MACP0A1(prim,i,j,k,pl) = MACP0A1(prim,ri,rj,rk+1+k,pl);
2354 
2355 
2356 #if(DEBUGINOUTLOOPS)
2357  dualfprintf(fail_file,"INNER X3: ispstag=%d pl=%d :: ri=%d rj=%d rk=%d (otherrk=%d) i=%d j=%d k=%d\n",ispstag,pl,ri,rj,rk,rk+SHIFT3+k,i,j,k);
2358  if(!isfinite(MACP0A1(prim,ri,rj,rk+SHIFT3+k,pl))){
2359  dualfprintf(fail_file,"INNER X3: caught copying nan ri=%d rj=%d rk=%d pl=%d\n",ri,rj,rk,pl);
2360  }
2361 #endif
2362 
2363  }
2364 
2365  // copy from lower side to upper boundary zones
2366  ri=i;
2367  rj=j;
2368  rk=rkin;
2369  LOOPBOUND3OUT PBOUNDLOOP(pliter,pl){
2370  MACP0A1(prim,i,j,k,pl) = MACP0A1(prim,ri,rj,rk+(k-N3),pl);
2371 
2372 #if(DEBUGINOUTLOOPS)
2373  dualfprintf(fail_file,"OUTER X3: ispstag=%d pl=%d :: ri=%d rj=%d rk=%d (otherrk=%d) i=%d j=%d k=%d\n",ispstag,pl,ri,rj,rk,rk+(k-N3),i,j,k);
2374  if(!isfinite(MACP0A1(prim,ri,rj,rk+(k-N3),pl))){
2375  dualfprintf(fail_file,"INNER X3: caught copying nan ri=%d rj=%d rk=%d pl=%d\n",ri,rj,rk,pl);
2376  }
2377 #endif
2378 
2379  }
2380  }
2381  }
2382 
2383 
2384  }
2385  else{
2386  dualfprintf(fail_file,"Shouldn't be here in bounds\n");
2387  myexit(3946851);
2388  }
2389 
2390 
2391  }// end parallel region
2392 
2393  return(0);
2394 }
2395 
2396 
2397 
2398 
2401  int boundstage, int finalstep, SFTYPE boundtime, int whichdir, int boundvartype, int *dirprim, int ispstag, FTYPE (*prim)[NSTORE2][NSTORE3][NPR],
2402  int *inboundloop,
2403  int *outboundloop,
2404  int *innormalloop,
2405  int *outnormalloop,
2406  int (*inoutlohi)[NUMUPDOWN][NDIM],
2407  int riin, int riout, int rjin, int rjout, int rkin, int rkout,
2408  int *dosetbc,
2409  int enerregion,
2410  int *localenerpos
2411  )
2412 {
2413  int i,j,k,pl,pliter;
2414  int locali1,globali1;
2415  int locali2,globali2;
2416  int ri1,ri2;
2417  FTYPE vcon[NDIM]; // coordinate basis vcon
2418 #if(WHICHVEL==VEL3)
2419  int failreturn;
2420 #endif
2421  int ri, rj, rk; // reference i,j,k
2422  FTYPE prescale[NPR];
2423  int horizonti;
2424  int jj,kk;
2425 
2426 
2427 
2428  if( (totalsize[1]>1) && BCtype[X1DN]==R0SING){ // global condition that all CPUs know about
2429 
2431  //
2432  // now deal with region interior to the horizon
2433  // This just makes values reasonable for diagnostics to be ignorant of loop ranges -- results are NOT used! -- no FLUXCTSTAG issues
2434  // similar to bound_spacetime_inside_horizon() in metric.c
2435 
2436 
2437  horizonti=N1*horizoncpupos1+horizoni;
2438 
2439  if(horizonti>0){ // global condition that all CPUs know about
2440  // for self-gravity, only loop over region outside horizon
2441  enerregion=OUTSIDEHORIZONENERREGION;
2442  enerpos=enerposreg[enerregion];
2443 
2444 
2445  // copy from very outer boundary or from location on grid defined by where last ghost cell is inside horizon
2446 
2447  // this doesn't do anything to cells inside horizon but outside fake region
2448  // apparently in this case uu1>0 even though inside the horizon
2449  globali1=N1*horizoncpupos1+horizoni-N1BND;
2450  if(globali1<startpos[1]-N1BND) locali1=-N1BND;
2451  else if(globali1>endpos[1]+N1BND-1) locali1=N1+N1BND-1;
2452  else locali1=horizoni-N1BND;
2453  ri1=locali1;
2454 
2455 
2456 
2457  // include all zones inside horizon as outflow so remove any peak that was accreted as black hole
2458  // flux through horizon of mass should be consistent with either method, but seems more reasonable to do below
2459  globali2=N1*horizoncpupos1+horizoni;
2460  if(globali2<startpos[1]) locali2=0;
2461  else if(globali2>endpos[1]-1) locali2=N1-1;
2462  else locali2=horizoni;
2463  ri2=locali2;
2464  // GODMARK: somehow leads to density sticking to value despite much lower active density
2465 
2466 
2467  LOOPF{
2468 
2469  rj=j;
2470  rk=k;
2471 
2472  // if(WITHINENERREGION(enerpos,i,j,k)){
2473  // then do nothing
2474  // }
2475  //else{
2476  // assume horizon on negative side of "i" so don't modify right side of "i" that happens to be in the outer radial boundary
2477 
2478  if(i < globali1-startpos[1]){
2479  // then inside horizon
2480 
2481  // set primitives to be trivialized
2482 
2483  PBOUNDLOOP(pliter,pl) MACP0A1(prim,i,j,k,pl) = MACP0A1(prim,ri1,rj,rk,pl);
2484 
2485  // reset pflag in case old failure inside horizon
2486 #if(UTOPRIMADJUST==UTOPRIMAVG)
2487  GLOBALMACP0A1(pflag,i,j,k,FLAGUTOPRIMFAIL)=UTOPRIMNOFAIL;
2488 #endif
2489 
2490 
2491  if(ispstag==0){
2492  // making it small artificially forces timestep to be small since timestep still limited within fake boundary region
2493  // MACP0A1(prim,i,j,k,RHO) = SMALL;
2494  // assumed acting on relative 4-velocity that stays physical no matter what
2495  if(MACP0A1(prim,i,j,k,U1)>0) MACP0A1(prim,i,j,k,U1)=-0.5*(fabs(MACP0A1(prim,ri1,rj,rk,U1))+fabs(MACP0A1(prim,ri1+1,rj,rk,U1)));
2496  }
2497  }
2498 
2499  if(0 && i < globali2-startpos[1]){
2500  // then on or inside horizon and inside N1BND more grid cells for interpolation
2501 
2502  if(ispstag==0){
2503  // assumed acting on relative 4-velocity that stays physical no matter what
2504  if(MACP0A1(prim,i,j,k,U1)>0) MACP0A1(prim,i,j,k,U1)=-0.5*(fabs(MACP0A1(prim,ri1,rj,rk,U1))+fabs(MACP0A1(prim,ri1+1,rj,rk,U1)));
2505  }
2506 
2507  }
2508  //}
2509  }
2510 
2511 
2512  }
2513  }
2514  else{
2515  dualfprintf(fail_file,"Shouldn't be here in bounds\n");
2516  myexit(3946852);
2517  }
2518 
2519 
2520 
2521  return(0);
2522 }
2523 
2524 
2525 
2526 
2527 
2528 
2529 
2532  int boundstage, int finalstep, SFTYPE boundtime, int whichdir, int boundvartype, int *dirprim, int ispstag, FTYPE (*prim)[NSTORE2][NSTORE3][NPR],
2533  int *inboundloop,
2534  int *outboundloop,
2535  int *innormalloop,
2536  int *outnormalloop,
2537  int (*inoutlohi)[NUMUPDOWN][NDIM],
2538  int riin, int riout, int rjin, int rjout, int rkin, int rkout,
2539  int *dosetbc,
2540  int enerregion,
2541  int *localenerpos
2542  )
2543 {
2544  int i,j,k,pl,pliter;
2545  int locali1,globali1;
2546  int locali2,globali2;
2547  int ri1,ri2;
2548  FTYPE vcon[NDIM]; // coordinate basis vcon
2549 #if(WHICHVEL==VEL3)
2550  int failreturn;
2551 #endif
2552  int ri, rj, rk; // reference i,j,k
2553  FTYPE prescale[NPR];
2554  int horizonti;
2555  int jj,kk;
2556 
2557 
2558 
2559 #if(PRODUCTION==0)
2560  // can only check after all directions done, and we assume whichdir==3 last to be done
2561  int trigger=0;
2562  if(ispstag){
2563  ZSLOOP(inoutlohi[POINTDOWN][POINTDOWN][1], inoutlohi[POINTUP][POINTUP][1],inoutlohi[POINTDOWN][POINTDOWN][2], inoutlohi[POINTUP][POINTUP][2],inoutlohi[POINTDOWN][POINTDOWN][3], inoutlohi[POINTUP][POINTUP][3]){
2564  // Can't use FULLLOOP since different boundary types go different depths into boundary cells
2565  // PALLLOOP(pl){
2566  PBOUNDLOOP(pliter,pl){
2567  // for staggered field, avoid inner-most cell since isn't needed or computed and so can be NaN
2568  if(pl==B1 && i>=-N1BND+SHIFT1 || pl==B2 && j>=-N2BND+SHIFT2 || pl==B3 && k>=-N3BND+SHIFT3){
2569  if(!finite(MACP0A1(prim,i,j,k,pl))){
2570  trigger++;
2571  dualfprintf(fail_file,"whichdir=%d ispstag=%d trigger=%d :: BC didn't set properly: #1: i=%d j=%d k=%d ti=%d tj=%d tk=%d pl=%d\n",whichdir,ispstag,trigger,i,j,k,startpos[1]+i,startpos[2]+j,startpos[3]+k,pl);
2572  }
2573  }
2574  }
2575  }
2576  }
2577  else{
2578  ZSLOOP(inoutlohi[POINTDOWN][POINTDOWN][1], inoutlohi[POINTUP][POINTUP][1],inoutlohi[POINTDOWN][POINTDOWN][2], inoutlohi[POINTUP][POINTUP][2],inoutlohi[POINTDOWN][POINTDOWN][3], inoutlohi[POINTUP][POINTUP][3]){
2579  // PALLLOOP(pl){
2580  PBOUNDLOOP(pliter,pl){
2581  if(!finite(MACP0A1(prim,i,j,k,pl))){
2582  trigger++;
2583  dualfprintf(fail_file,"whichdir=%d ispstag=%d trigger=%d :: BC didn't set properly: #1: i=%d j=%d k=%d ti=%d tj=%d tk=%d pl=%d\n",whichdir,ispstag,trigger,i,j,k,startpos[1]+i,startpos[2]+j,startpos[3]+k,pl);
2584  }
2585  }
2586  }
2587  }
2588  // if(trigger>0) myexit(92489346); // first calls to bound will have no field, for example, so can't exit
2589 #endif
2590 
2591 
2592 
2593  return (0);
2594 }
2595 
2596 
2597 
2598 
2599 
2602 //
2605 int extrapfunc(int boundary, int j,int k,
2606  int boundstage, int finalstep, SFTYPE boundtime, int whichdir, int boundvartype, int *dirprim, int ispstag, FTYPE (*prim)[NSTORE2][NSTORE3][NPR],
2607  int *inboundloop,
2608  int *outboundloop,
2609  int *innormalloop,
2610  int *outnormalloop,
2611  int (*inoutlohi)[NUMUPDOWN][NDIM],
2612  int riin, int riout, int rjin, int rjout, int rkin, int rkout,
2613  int *dosetbc,
2614  int enerregion,
2615  int *localenerpos
2616  )
2617 {
2618  int i,pl,pliter;
2619  int ri,rj,rk;
2620  int iloopstart,iloopend,iloopstep;
2621  int ri2,ri3;
2622  // for Bd3 copy
2623  FTYPE Bd3,Bd3ri2,Bd3ri3;
2624  int jj;
2625  FTYPE Bu1,Bu2,gcon03,gcon13,gcon23,gcon33;
2626  FTYPE gcov01,gcov02,gcov11,gcov12,gcov21,gcov22,gcov03,gcov13,gcov23;
2627  // Mdot copy
2628  FTYPE Mdot,ucon[NDIM],uconref[NDIM],uradconref[NDIM];
2629  FTYPE Mraddot,uradcon[NDIM];
2630  FTYPE Xr[NPR][NDIM],Vr[NPR][NDIM],X[NPR][NDIM],V[NPR][NDIM];
2631  FTYPE primtemp[NPR];
2632  FTYPE *ucontouse, ucon2[NDIM];
2633  int itervphi;
2634  FTYPE uconrefri2[NDIM],uradconrefri2[NDIM];
2635  FTYPE Mdotri2;
2636  FTYPE uconrefri3[NDIM],uradconrefri3[NDIM];
2637  FTYPE othersref[NUMOTHERSTATERESULTS],othersrefri2[NUMOTHERSTATERESULTS],othersrefri3[NUMOTHERSTATERESULTS],others[NUMOTHERSTATERESULTS],othersrefnew[NUMOTHERSTATERESULTS];
2638  FTYPE othersradref[NUMOTHERSTATERESULTS],othersradrefri2[NUMOTHERSTATERESULTS],othersradrefri3[NUMOTHERSTATERESULTS],othersrad[NUMOTHERSTATERESULTS],othersradrefnew[NUMOTHERSTATERESULTS];
2639  FTYPE Mdotri3;
2640  FTYPE Dqp,Dqm,Dqc,dqMdot;
2641  FTYPE myMdot;
2642  FTYPE dq[NPR];
2643  FTYPE dqucon[NDIM];
2644  FTYPE dqBd3,myBd3;
2645  FTYPE dqgdetpl[NPR];
2646  FTYPE dqorthopl[NPR];
2647  FTYPE dqlogdensity[NPR];
2648  FTYPE uconreftouse[NDIM],uradconreftouse[NDIM];
2649  FTYPE xfrac,yfrac;
2650  FTYPE ftemp;
2651  FTYPE signdq;
2652  FTYPE xdqfrac,ydqfrac;
2653  FTYPE ftemp2,linearvalue,expvalue;
2654  FTYPE D0,uconrefnew[NDIM],uradconrefnew[NDIM];
2655  FTYPE mydqlog;
2656  FTYPE igdetnosing;
2657  struct of_geom geomdontuse[NPR];
2658  struct of_geom *ptrgeom[NPR];
2659  struct of_geom rgeomdontuse[NPR];
2660  struct of_geom *ptrrgeom[NPR];
2661  struct of_geom ri2geomdontuse[NPR];
2662  struct of_geom *ptrri2geom[NPR];
2663  struct of_geom ri3geomdontuse[NPR];
2664  struct of_geom *ptrri3geom[NPR];
2665 
2666 
2667  // assign memory
2668  PALLLOOP(pl){
2669  ptrgeom[pl]=&(geomdontuse[pl]);
2670  ptrrgeom[pl]=&(rgeomdontuse[pl]);
2671  ptrri2geom[pl]=&(ri2geomdontuse[pl]);
2672  ptrri3geom[pl]=&(ri3geomdontuse[pl]);
2673  }
2674 
2675 
2676  rj=j;
2677  rk=k;
2678 
2679 
2680  // setup multiple reference points
2681  if(boundary==X1DN){
2682  ri = riin;
2683  ri2=ri+1;
2684  ri3=ri+2;
2685  // iterate up (due to use of single loop form)
2686  iloopstart=LOWERBOUND1;
2687  iloopend=ri-1;
2688  iloopstep=1;
2689  signdq=1.0;
2690  }
2691  else if(boundary==X1UP){
2692  ri = riout;
2693  ri2=ri-1;
2694  ri3=ri-2;
2695  // iterate up (due to use of single loop form)
2696  iloopstart=ri+1;
2697  iloopend=UPPERBOUND1;
2698  iloopstep=1;
2699  signdq=-1.0;
2700  }
2701  else{
2702  dualfprintf(fail_file,"extrapfunc not setup for this boundary = %d\n",boundary);
2703  myexit(811751);
2704  }
2705 
2706 
2707 
2709  //
2710  // get coordinates of reference location
2711  //
2713  PBOUNDLOOP(pliter,pl){
2714  get_geometry(ri, rj, rk, dirprim[pl], ptrrgeom[pl]);
2715  get_geometry(ri2, rj, rk, dirprim[pl], ptrri2geom[pl]);
2716  get_geometry(ri3, rj, rk, dirprim[pl], ptrri3geom[pl]);
2717 
2718  bl_coord_ijk_2(ri,rj,rk,dirprim[pl],Xr[pl],Vr[pl]); // reference locations for B2/U2
2719  }
2720 
2721 
2722 
2724  //
2725  // obtain reference B_\phi
2726  //
2728  Bd3=0.0; SLOOPA(jj) Bd3 += MACP0A1(prim,ri,rj,rk,B1+jj-1)*(ptrrgeom[B1+jj-1]->gcov[GIND(3,jj)]);
2729  Bd3ri2=0.0; SLOOPA(jj) Bd3ri2 += MACP0A1(prim,ri2,rj,rk,B1+jj-1)*(ptrri2geom[B1+jj-1]->gcov[GIND(3,jj)]);
2730  Bd3ri3=0.0; SLOOPA(jj) Bd3ri3 += MACP0A1(prim,ri3,rj,rk,B1+jj-1)*(ptrri3geom[B1+jj-1]->gcov[GIND(3,jj)]);
2731  // B_\phi slope
2732  Dqp=Bd3ri3-Bd3ri2;
2733  Dqm=Bd3ri2-Bd3;
2734  Dqc=Bd3ri3-Bd3;
2735  Dqc*=0.5;
2736  dqBd3 = signdq*MINMOD(MINMOD(Dqp,Dqm),Dqc);
2737 
2738 
2739 
2741  //
2742  // reference value for \detg rho u^{x1}
2743  //
2745 
2746 
2747  if(ispstag==0){
2748  ucon_calc(MAC(prim,ri,rj,rk),ptrrgeom[U1],uconref,othersref);
2749  Mdot = (ptrrgeom[U1]->gdet)*MACP0A1(prim,ri,rj,rk,RHO)*uconref[1];
2750 
2751  // ri2
2752  ucon_calc(MAC(prim,ri2,rj,rk),ptrri2geom[U1],uconrefri2,othersrefri2);
2753  Mdotri2 = (ptrri2geom[U1]->gdet)*MACP0A1(prim,ri2,rj,rk,RHO)*uconrefri2[1];
2754 
2755  // ri3
2756  ucon_calc(MAC(prim,ri3,rj,rk),ptrri3geom[U1],uconrefri3,othersrefri3);
2757  Mdotri3 = (ptrri3geom[U1]->gdet)*MACP0A1(prim,ri3,rj,rk,RHO)*uconrefri3[1];
2758 
2759  // Mdot slope
2760  Dqp=Mdotri3-Mdotri2;
2761  Dqm=Mdotri2-Mdot;
2762  Dqc=Mdotri3-Mdot;
2763  Dqc*=0.5;
2764  dqMdot = signdq*MINMOD(MINMOD(Dqp,Dqm),Dqc);
2765 
2766  if(EOMRADTYPE!=EOMRADNONE){
2767  // KORALTODO: Could treat radiation stuff equal to MHD, but here MHD is just getting Mdot, that has no real counterpart for radiation
2768  ucon_calc(&MACP0A1(prim,ri,rj,rk,URAD1-U1),ptrrgeom[URAD1],uradconref,othersradref);
2769  ucon_calc(&MACP0A1(prim,ri2,rj,rk,URAD1-U1),ptrri2geom[URAD1],uradconrefri2,othersradrefri2);
2770  ucon_calc(&MACP0A1(prim,ri3,rj,rk,URAD1-U1),ptrri3geom[URAD1],uradconrefri3,othersradrefri3);
2771  }
2772 
2773  }
2774 
2775 
2776 
2778  //
2779  // prim[pl] slope
2780  //
2782  PBOUNDLOOP(pliter,pl){
2783  // determine MINM slope for extrapolation
2784  Dqp=MACP0A1(prim,ri3,rj,rk,pl)-MACP0A1(prim,ri2,rj,rk,pl);
2785  Dqm=MACP0A1(prim,ri2,rj,rk,pl)-MACP0A1(prim,ri,rj,rk,pl);
2786  Dqc=MACP0A1(prim,ri3,rj,rk,pl)-MACP0A1(prim,ri,rj,rk,pl);
2787  Dqc*=0.5;
2788  dq[pl] = signdq*MINMOD(MINMOD(Dqp,Dqm),Dqc);
2789  }
2790 
2791 
2793  //
2794  // log(prim[pl]) slope
2795  //
2797  PBOUNDLOOP(pliter,pl){
2798  // determine MINM slope for extrapolation
2799  Dqp=log(SMALL+fabs(MACP0A1(prim,ri3,rj,rk,pl)))-log(SMALL+fabs(MACP0A1(prim,ri2,rj,rk,pl)));
2800  Dqm=log(SMALL+fabs(MACP0A1(prim,ri2,rj,rk,pl)))-log(SMALL+fabs(MACP0A1(prim,ri,rj,rk,pl)));
2801  Dqc=log(SMALL+fabs(MACP0A1(prim,ri3,rj,rk,pl)))-log(SMALL+fabs(MACP0A1(prim,ri,rj,rk,pl)));
2802  Dqc*=0.5;
2803  dqlogdensity[pl] = signdq*MINMOD(MINMOD(Dqp,Dqm),Dqc);
2804  }
2805 
2807  //
2808  // gdet * prim[pl] slope
2809  //
2811  PBOUNDLOOP(pliter,pl){
2812  // determine MINM slope for extrapolation
2813  Dqp=(ptrri3geom[pl]->gdet)*MACP0A1(prim,ri3,rj,rk,pl)-(ptrri2geom[pl]->gdet)*MACP0A1(prim,ri2,rj,rk,pl);
2814  Dqm=(ptrri2geom[pl]->gdet)*MACP0A1(prim,ri2,rj,rk,pl)-(ptrrgeom[pl]->gdet)*MACP0A1(prim,ri,rj,rk,pl);
2815  Dqc=(ptrri3geom[pl]->gdet)*MACP0A1(prim,ri3,rj,rk,pl)-(ptrrgeom[pl]->gdet)*MACP0A1(prim,ri,rj,rk,pl);
2816  Dqc*=0.5;
2817  dqgdetpl[pl] = signdq*MINMOD(MINMOD(Dqp,Dqm),Dqc);
2818  }
2819 
2821  //
2822  // sqrt(gv3ii) * prim[pl] slope
2823  //
2825  int ddir;
2826  PBOUNDLOOP(pliter,pl){
2827  // determine MINM slope for extrapolation
2828  ddir=0; // default
2829  if(pl>=U1 || pl<=U3) ddir=pl-U1+1;
2830  if(pl>=URAD1 || pl<=URAD3) ddir=pl-URAD1+1;
2831  if(pl>=B1 || pl<=B3) ddir=pl-B1+1;
2832  Dqp=sqrt(fabs(ptrri3geom[pl]->gcov[GIND(ddir,ddir)]))*MACP0A1(prim,ri3,rj,rk,pl)-sqrt(fabs(ptrri2geom[pl]->gcov[GIND(ddir,ddir)]))*MACP0A1(prim,ri2,rj,rk,pl);
2833  Dqm=sqrt(fabs(ptrri2geom[pl]->gcov[GIND(ddir,ddir)]))*MACP0A1(prim,ri2,rj,rk,pl)-sqrt(fabs(ptrrgeom[pl]->gcov[GIND(ddir,ddir)]))*MACP0A1(prim,ri,rj,rk,pl);
2834  Dqc=sqrt(fabs(ptrri3geom[pl]->gcov[GIND(ddir,ddir)]))*MACP0A1(prim,ri3,rj,rk,pl)-sqrt(fabs(ptrrgeom[pl]->gcov[GIND(ddir,ddir)]))*MACP0A1(prim,ri,rj,rk,pl);
2835  Dqc*=0.5;
2836  dqorthopl[pl] = signdq*MINMOD(MINMOD(Dqp,Dqm),Dqc);
2837  }
2838 
2839 
2841  //
2842  // u^i slope
2843  //
2845  if(ispstag==0){
2846  DLOOPA(jj){
2847  // determine MINM slope for extrapolation
2848  Dqp=uconrefri3[jj]-uconrefri2[jj];
2849  Dqm=uconrefri2[jj]-uconref[jj];
2850  Dqc=uconrefri3[jj]-uconref[jj];
2851  Dqc*=0.5;
2852  dqucon[jj] = signdq*MINMOD(MINMOD(Dqp,Dqm),Dqc);
2853  }
2854 
2855  // get fraction of farther away reference value to use
2856  // Fraction of normal reference value is determined by how similar u^t is for (e.g. for inner boundary) ri and ri+1
2857  xfrac = fabs(uconref[TT]-uconrefri2[TT])/uconref[TT];
2858 #define XFRAC1 (0.1) // 0.1 was chosen from experience with BH+torus and what the u^t value was doing at the pole near the horizon -- JCM found once d(u^t)>0.1 then started to grow unbound. Certainly related to resolution so assume resolving of horizon is generally similar to the 64^2 case with Rout=40 and exp grid
2859 #define YFRAC1 (0.0)
2860 #define XFRAC2 (0.2)
2861 #define YFRAC2 (1.0)
2862 #define YFRAC12(frac) (YFRAC1 + (YFRAC2-YFRAC1)/(XFRAC2-XFRAC1)*(frac-XFRAC1))
2863 
2864  if(xfrac<XFRAC1) yfrac=YFRAC1;
2865  else if(xfrac<XFRAC2) yfrac=YFRAC12(xfrac);
2866  else yfrac=YFRAC2;
2867 
2868 
2869 #define XDQFRAC1 (0.1)
2870 #define YDQFRAC1 (0.0)
2871 #define XDQFRAC2 (0.2)
2872 #define YDQFRAC2 (1.0)
2873 #define YDQFRAC12(frac) (YDQFRAC1 + (YDQFRAC2-YDQFRAC1)/(XDQFRAC2-XDQFRAC1)*(frac-XDQFRAC1))
2874 
2875  // limit slope in extrapolation if excessive slope
2876  xdqfrac = fabs(uconref[TT]-uconrefri2[TT])/uconref[TT];
2877  if(xdqfrac<XDQFRAC1) ydqfrac=YDQFRAC1;
2878  else if(xdqfrac<XDQFRAC2) ydqfrac= YDQFRAC12(xdqfrac);
2879  else ydqfrac=YDQFRAC2;
2880 
2881  // modify dqucon based upon ydqfrac
2882  // Uses original MINM slope unless slope is too large, then start to decrease to no slope
2883  // only apply to velocity related slopes
2884  DLOOPA(jj){
2885  dqucon[jj] = dqucon[jj]*(1.0-ydqfrac);
2886  dq[UU+jj] = dq[UU+jj]*(1.0-ydqfrac);
2887  dqlogdensity[UU+jj] = dqlogdensity[UU+jj]*(1.0-ydqfrac);
2888  dqMdot = dqMdot*(1.0-ydqfrac);
2889  dqgdetpl[UU+jj] = dqgdetpl[UU+jj]*(1.0-ydqfrac);
2890  }
2891 
2892 
2893 
2894  }// end if ispstag==0
2895  else{
2896  yfrac=ydqfrac=0.0;
2897  }
2898 
2899 
2900 
2901 
2902 
2903 
2904 
2906  //
2907  // LOOP over i and pl
2908  //
2910 
2911  LOOPBOUNDGENMORE(i,iloopstart,iloopend,iloopstep){
2912 
2913 
2914 
2916  //
2917  // default copy
2918  //
2920  PBOUNDLOOP(pliter,pl) MACP0A1(prim,i,j,k,pl)=MACP0A1(prim,ri,rj,rk,pl);
2921 
2923  //
2924  // get coordinates of local location for all quantities
2925  //
2927  PBOUNDLOOP(pliter,pl){
2928  get_geometry(i, j, k, dirprim[pl], ptrgeom[pl]);
2929  bl_coord_ijk_2(i,j,k,dirprim[pl],X[pl],V[pl]); // reference locations for B2/U2
2930  }
2931 
2932 
2933 
2934 
2935 
2937  //
2938  // NON-magnetic field parts (\rho,u,u^i, etc.)
2939  //
2941  if(ispstag==0){
2942 
2943 
2945  //
2946  // Other Densities
2947  //
2949 
2951  // now focus on accuracy for polar regions to maintain stability (avoid dissipation terms doing funny things)
2952 
2953  // using below Bondi-like scaling with Mdot scaling for ucon leads to space-like answers Near BH densities should become roughly constant, and this doesn't lead to problems, so use instead
2954  // pl=RHO; MACP0A1(prim,i,j,k,pl) = MACP0A1(prim,ri,rj,rk,pl)*pow(V[pl][1]/Vr[pl][1],-3.0/2.0);
2955  // pl=UU; MACP0A1(prim,i,j,k,pl) = MACP0A1(prim,ri,rj,rk,pl)*pow(V[pl][1]/Vr[pl][1],-5.0/2.0);
2956 
2957 #define POWEREXTRAFRAC (0.1)
2958 #define POWERRATIO (1.0+POWEREXTRAFRAC) // ratio of densities not allow to go bigger than this when using dqlogdensity
2959 
2960 
2961 #if(1)
2962 
2963  PBOUNDLOOP(pliter,pl){
2964  // MUST BE POSITIVE DEFINITE!!!
2965  if(pl!=B1 && pl!=B2 && pl!=B3 && pl!=URAD1 && pl!=URAD2 && pl!=URAD3 && SCALARPL(pl)==0){
2966 
2967  // only use linear if exponentiation causes growth of value, not decreasion
2968  ftemp = exp(-signdq*dqlogdensity[pl]) - POWERRATIO;
2969 
2970  // limit interpolation strength to fixed log-slope
2971  if(ftemp>=0.0) mydqlog=log(POWERRATIO)/(-signdq);
2972  else mydqlog = dqlogdensity[pl];
2973 
2974  // log extrap (very speculative and can cause problems if used alone when (say) density is super low on ri+1 and relatively high on ri, then i will be super huge
2975  // should use this when values that go into slope are much different, or equally when dqlogdensity is large
2976  MACP0A1(prim,i,j,k,pl) = exp(log(SMALL+fabs(MACP0A1(prim,ri,rj,rk,pl))) + mydqlog*(i-ri));
2977 
2978  // DEBUG:
2979  // dualfprintf(fail_file,"i=%d j=%d pl=%d ftemp=%21.15g linearvalue=%21.15g expvalue=%21.15g final=%21.15g\n",i,j,pl,ftemp,linearvalue,expvalue,MACP0A1(prim,i,j,k,pl));
2980 
2981  }
2982  }
2983 
2984 
2985 
2986 #else
2987 
2988 
2989 
2990  // override using OLD WAY for densities
2991 #if(0)
2992  // linear extrap (leads to negative values, which can cause problems with flux inversions)
2993  PBOUNDLOOP(pliter,pl) if(pl==RHO) MACP0A1(prim,i,j,k,pl) = MACP0A1(prim,ri,rj,rk,pl) + dq[pl]*(i-ri);
2994  PBOUNDLOOP(pliter,pl) if(pl==UU) MACP0A1(prim,i,j,k,pl) = MACP0A1(prim,ri,rj,rk,pl) + dq[pl]*(i-ri);
2995  PBOUNDLOOP(pliter,pl) if(pl==URAD0) MACP0A1(prim,i,j,k,pl) = MACP0A1(prim,ri,rj,rk,pl) + dq[pl]*(i-ri);
2996 #else
2997  // log extrap
2998  PBOUNDLOOP(pliter,pl) if(pl==RHO) MACP0A1(prim,i,j,k,pl) = exp(log(SMALL+fabs(MACP0A1(prim,ri,rj,rk,pl))) + dqlogdensity[pl]*(i-ri));
2999  PBOUNDLOOP(pliter,pl) if(pl==UU) MACP0A1(prim,i,j,k,pl) = exp(log(SMALL+fabs(MACP0A1(prim,ri,rj,rk,pl))) + dqlogdensity[pl]*(i-ri));
3000  PBOUNDLOOP(pliter,pl) if(pl==URAD0) MACP0A1(prim,i,j,k,pl) = exp(log(SMALL+fabs(MACP0A1(prim,ri,rj,rk,pl))) + dqlogdensity[pl]*(i-ri));
3001 #endif
3002 
3003 
3004 #endif
3005 
3006 
3007 
3008 
3010  //
3011  // Velocity
3012  //
3014 
3015 #if(0)
3016 
3017 
3018  // appears too speculative and use of v^3 causes ucon2pr to lead to no solution within boundary conditions -- suggestive of some inconsistency
3019  // (primitive can be anything)
3020 
3021  // combined with \rho_0=constant, this implies \detg \rho_0 u^x1 == constant as required for Mdot=constant for radial flow
3022  // pl=U1; MACP0A1(prim,i,j,k,pl) = MACP0A1(prim,ri,rj,rk,pl)*fabs((ptrrgeom[pl]->gdet)/(ptrgeom[pl]->gdet));
3023  myMdot = Mdot + dqMdot*(i-ri);
3024  ucon[1] = myMdot/((ptrgeom[U1]->gdet)*MACP0A1(prim,i,j,k,RHO));
3025  // assume ucon[2] and ucon[3] are just copied directly
3026  // assume mostly radial flow, but visibly appears like B2 so treat similarly so no large kink in behavior
3027  ucon[2] = uconref[2] + dqucon[2]*(i-ri);
3028  // assume v^\phi \sim \Omega doesn't change much
3029  ucon[3] = uconref[3] + dqucon[3]*(i-ri);
3030 
3031  // using non-computed u^t in ucon2pr means slightly inconsistent with desired u^t[u^i], but assume approximately correct and is sufficient for extrapolation to get primitive
3032  // ucon[TT] = uconref[TT];
3033 
3034  // fill so can use standard functions
3035  primtemp[U1] = ucon[1];
3036  primtemp[U2] = ucon[2];
3037  primtemp[U3] = ucon[3];
3038  // get real ucon with real u^t
3039  ucon_calc_4vel_bothut(primtemp, ptrgeom[U1], ucon, ucon2, others);
3040  // now get primitive velocity (choosing 4-velocity branch closest to reference velocity)
3041  // using the below causes flipping of choice in accretion flow region of boundary conditions due to being too far away from original uconref[TT]
3042  // if(fabs(ucon[TT]-uconref[TT])<fabs(ucon2[TT]-uconref[TT])) ucontouse=ucon;
3043  // else ucontouse=ucon2;
3044  ucontouse=ucon;
3045 
3047  //
3048  // iterate to get closer to v^\phi\sim constant (safer than setting v^i directly and using vcon2pr())
3049  //
3051 
3052  for(itervphi=0;itervphi<=NUMITERVPHI;itervphi++){
3053  // assume v^\phi \sim \Omega doesn't change much
3054  primtemp[U3] = uconref[3]/uconref[TT]*ucontouse[TT];
3055  // get real ucon with real u^t
3056  ucon_calc_4vel_bothut(primtemp, ptrgeom[U1], ucon, ucon2, others);
3057  // now get primitive velocity (choosing 4-velocity branch closest to reference velocity)
3058  // if(fabs(ucon[TT]-uconref[TT])<fabs(ucon2[TT]-uconref[TT])) ucontouse=ucon;
3059  // else ucontouse=ucon2;
3060  ucontouse=ucon;
3061  }
3062 
3064  // get final primitive from ucon
3065  ucon2pr(WHICHVEL,ucontouse,ptrgeom[U1],MAC(prim,i,j,k)); // fills only U1,U2,U3
3066 
3067 
3068 #elif(0)
3069 
3070  // GODMARK: Was causing problems for Rebecca
3071 
3072 
3073  // BH-pole corner angular sucking drive u^t>>0 and then failures occur
3074  // avoid exaggerating via unwanted radial sucking by using reference value with
3075 
3076  // whether to conserve (at least) D=\rho_0 u^t when modifying \gamma
3077  // see notes in fixup.c for limit_gamma() regarding why this is not a good idea
3078 #define DO_CONSERVE_D_INBOUNDS 0 // CHANGINGMARK
3079 
3080 
3081 
3082 #if(DO_CONSERVE_D_INBOUNDS)
3083  // GODMARK: here we actually modify active values of quantities (should we preserve D=\rho_0 u^t ? -- i.e. change \rho_0?)
3084  // This correction comes after any other velocity modifications
3085  // of all conservation laws, at least conserve particle number
3086  D0 = MACP0A1(prim,ri,rj,rk,RHO)*uconref[TT];
3087 #endif
3088 
3089  PBOUNDLOOP(pliter,pl){
3090  if(pl==U1 || pl==U2 || pl==U3 || pl==URAD1 || pl==URAD2 || pl==URAD3){
3091 
3092  // switch away from using near-BC value as reference if going bad since don't want to exaggerate it
3093  ftemp = MACP0A1(prim,ri,rj,rk,pl)*(1.0-yfrac) + MACP0A1(prim,ri2,rj,rk,pl)*yfrac;
3094  // interpolate
3095  MACP0A1(prim,i,j,k,pl) = (ftemp + dq[pl]*(i-ri));
3096 
3097 
3098  // dualfprintf(fail_file,"i=%d j=%d k=%d pl=%d ftemp=%21.15g yfrac=%21.15g dq[pl]=%21.15g i=%d ri=%d ri2=%d ri3=%d rj=%d rk=%d :: prim=%21.15g pri2=%21.15g pri3=%21.15g\n",i,j,k,pl,ftemp,yfrac,dq[pl],i,ri,ri2,ri3,rj,rk,MACP0A1(prim,i,j,k,pl),MACP0A1(prim,ri2,rj,rk,pl),MACP0A1(prim,ri3,rj,rk,pl)); // CHANGINGMARK
3099 
3100 
3101 
3102  if(V[pl][RR]<RADIUSTOAVOIDRADIALSUCK){
3103  // only do this if near a BH
3104  // also use yfrac to reset on-grid value since apparently it's going bad
3105  // this also keeps interpolation scheme passing through to boundary and so avoiding unwanted dissipation that can evolve flow away from sanity
3106  // linear extrap
3107  // GODMARK: Should probably make sure that u^r doesn't change sign?!
3108  MACP0A1(prim,ri,rj,rk,pl) = MACP0A1(prim,ri,rj,rk,pl)*(1.0-yfrac) + (MACP0A1(prim,ri2,rj,rk,pl) + dq[pl]*(ri-(ri2)))*yfrac;
3109  }
3110 
3111  }// end over velocities
3112  }
3113 #if(DO_CONSERVE_D_INBOUNDS)
3114  // get new u^t
3115  ucon_calc(MAC(prim,ri,rj,rk),ptrrgeom[U1],uconrefnew, othersrefnew);
3116  // recompute \rho_0 to conserve particle number (generally will increase \rho_0)
3117  MACP0A1(prim,ri,rj,rk,RHO) = D0/uconrefnew[TT];
3118 #endif
3119 
3120 
3121 
3122 
3123 #elif(0)
3124  // linear extrap
3125  PBOUNDLOOP(pliter,pl){
3126  if(pl==U1 || pl==U2 || pl==U3 || pl==URAD1 || pl==URAD2 || pl==URAD3){
3127 
3128  // switch away from using near-BC value as reference if going bad since don't want to exaggerate it
3129  ftemp = MACP0A1(prim,ri,rj,rk,pl)*(1.0-yfrac) + MACP0A1(prim,ri3,rj,rk,pl)*yfrac;
3130  // interpolate
3131  MACP0A1(prim,i,j,k,pl) = ftemp + dq[pl]*(i-ri);
3132  }
3133  }
3134 
3135 
3136 #elif(1)
3137  // linear extrap for velocities
3138  PBOUNDLOOP(pliter,pl){
3139  if(pl==U1 || pl==U2 || pl==U3 || pl==URAD1 || pl==URAD2 || pl==URAD3){
3140  //if(pl==U1 || pl==U2 || pl==U3){
3141  ftemp = MACP0A1(prim,ri,rj,rk,pl);
3142  // interpolate
3143  MACP0A1(prim,i,j,k,pl) = ftemp + dq[pl]*(i-ri);
3144  }
3145 
3146  if(0){
3147  // SUPERMARK: KORALTODO
3148  if(pl==URAD1 || pl==URAD2 || pl==URAD3){
3149  ftemp = MACP0A1(prim,ri,rj,rk,pl);
3150  // interpolate
3151  MACP0A1(prim,i,j,k,pl) = ftemp;
3152  }
3153  }
3154 
3155  }
3156 
3157 
3158 #if(0)
3159  ucon_calc(&MACP0A1(prim,ri,rj,rk,0),ptrgeom[U1],ucon,others);
3160  if(boundary==X1DN && (uconref[1]>0.0 || ucon[1]>0.0)) MACP0A1(prim,i,j,k,U1)=0.0;
3161  if(boundary==X1UP && (uconref[1]<0.0 || ucon[1]<0.0)) MACP0A1(prim,i,j,k,U1)=0.0;
3162  if(boundary==X2DN && (uconref[2]>0.0 || ucon[2]>0.0)) MACP0A1(prim,i,j,k,U2)=0.0;
3163  if(boundary==X2UP && (uconref[2]<0.0 || ucon[2]<0.0)) MACP0A1(prim,i,j,k,U2)=0.0;
3164  if(boundary==X3DN && (uconref[3]>0.0 || ucon[3]>0.0)) MACP0A1(prim,i,j,k,U3)=0.0;
3165  if(boundary==X3UP && (uconref[3]<0.0 || ucon[3]<0.0)) MACP0A1(prim,i,j,k,U3)=0.0;
3166 
3167  if(EOMRADTYPE!=EOMRADNONE){
3168  ucon_calc(&MACP0A1(prim,ri,rj,rk,URAD1-U1),ptrgeom[URAD1],uradcon,othersrad);
3169  if(boundary==X1DN && (uradconref[1]>0.0 || uradcon[1]>0.0)) MACP0A1(prim,i,j,k,URAD1)=0.0;
3170  if(boundary==X1UP && (uradconref[1]<0.0 || uradcon[1]<0.0)) MACP0A1(prim,i,j,k,URAD1)=0.0;
3171  if(boundary==X2DN && (uradconref[2]>0.0 || uradcon[2]>0.0)) MACP0A1(prim,i,j,k,URAD2)=0.0;
3172  if(boundary==X2UP && (uradconref[2]<0.0 || uradcon[2]<0.0)) MACP0A1(prim,i,j,k,URAD2)=0.0;
3173  if(boundary==X3DN && (uradconref[3]>0.0 || uradcon[3]>0.0)) MACP0A1(prim,i,j,k,URAD3)=0.0;
3174  if(boundary==X3UP && (uradconref[3]<0.0 || uradcon[3]<0.0)) MACP0A1(prim,i,j,k,URAD3)=0.0;
3175  }
3176 #endif
3177 
3178 
3179 #endif
3180 
3181 
3182  } // end if ispstag==0
3183 
3184 
3185 
3187  //
3188  // MAGNETIC FIELD
3189  //
3191 
3193  //
3194  // assume \detg B^x1 == constant (well, extrapolated now)
3195  //
3196  // GODMARK: Should probably make sure that B^r and B^\phi don't change sign! (well, maybe at least B^r)
3197  // Should probably preserve signature of B^\phi/B^r = -\Omega_F/c -- suggests interpolating B^\phi/B^r instead of B_\phi -- well, can't divide by B^r
3198  //
3200  if(VARTOINTERP==PRIMTOINTERP_GDETFULLVERSION_WALD && boundary==X1UP){ // WALD:
3201  int ddirl;
3202  for(pl=B1;pl<=B2;pl++){
3203  if(pl==B1) ddirl=1;
3204  if(pl==B2) ddirl=2;
3205  MACP0A1(prim,i,j,k,pl) = (MACP0A1(prim,ri,rj,rk,pl)*sqrt(fabs(ptrrgeom[pl]->gcov[GIND(ddirl,ddirl)])) + 0.0*dqorthopl[pl]*(i-ri))/sqrt(fabs(ptrgeom[pl]->gcov[GIND(ddirl,ddirl)]));
3206  // MACP0A1(prim,i,j,k,pl) = (MACP0A1(prim,ri,rj,rk,pl)*sqrt(fabs(ptrrgeom[pl]->gcov[GIND(ddirl,ddirl)])) + dqorthopl[pl]*(i-ri))/sqrt(fabs(ptrgeom[pl]->gcov[GIND(ddirl,ddirl)]));
3207  // MACP0A1(prim,i,j,k,pl) = (MACP0A1(prim,ri,rj,rk,pl));
3208 
3209  }
3210 
3211  }
3212  else{
3213  for(pl=B1;pl<=B2;pl++){
3214  igdetnosing=sign(ptrgeom[pl]->gdet)/(fabs(ptrgeom[pl]->gdet)+SMALL); // avoids 0.0 for any sign of ptrgeom->gdet
3215  MACP0A1(prim,i,j,k,pl) = (MACP0A1(prim,ri,rj,rk,pl)*(ptrrgeom[pl]->gdet) + dqgdetpl[pl]*(i-ri))*igdetnosing;
3216  }
3217  }
3218 
3219 
3220 
3221 #define EXTRAPBD3 0
3222 
3223  if(EXTRAPBD3){
3225  //
3226  // obtian generally correct B^\phi[B^r,B^\theta,B_\phi]
3227  //
3229  if(dirprim[B3]==FACE3){
3230  // then assume staggered method
3231  // need to average fields
3232  // Bu1=0.25*(MACP0A1(prim,i,j,k,B1)+MACP0A1(prim,ip1,j,k,B1)+MACP0A1(prim,i,j,km1,B1)+MACP0A1(prim,ip1,j,km1,B1));
3233  //Bu2=0.25*(MACP0A1(prim,i,j,k,B2)+MACP0A1(prim,i,jp1,k,B2)+MACP0A1(prim,i,j,km1,B2)+MACP0A1(prim,i,jp1,km1,B2));
3234  // assume average just results in the same value since can't average non-set values
3235  Bu1=MACP0A1(prim,i,j,k,B1);
3236  Bu2=MACP0A1(prim,i,j,k,B2);
3237  }
3238  else{
3239  // then assume all fields at CENT
3240  Bu1=MACP0A1(prim,i,j,k,B1);
3241  Bu2=MACP0A1(prim,i,j,k,B2);
3242  }
3243  gcon03=ptrgeom[B3]->gcon[GIND(0,3)];
3244  gcon13=ptrgeom[B3]->gcon[GIND(1,3)];
3245  gcon23=ptrgeom[B3]->gcon[GIND(2,3)];
3246  gcon33=ptrgeom[B3]->gcon[GIND(3,3)];
3247 
3248  gcov01=ptrgeom[B3]->gcov[GIND(0,1)];
3249  gcov02=ptrgeom[B3]->gcov[GIND(0,2)];
3250  gcov11=ptrgeom[B3]->gcov[GIND(1,1)];
3251  gcov12=gcov21=ptrgeom[B3]->gcov[GIND(1,2)];
3252  gcov22=ptrgeom[B3]->gcov[GIND(2,2)];
3253  gcov03=ptrgeom[B3]->gcov[GIND(0,3)];
3254  gcov13=ptrgeom[B3]->gcov[GIND(1,3)];
3255  gcov23=ptrgeom[B3]->gcov[GIND(2,3)];
3256 
3257  // Bd3fromBu.nb (just moved signs)
3258  myBd3=Bd3 + dqBd3*(i-ri);
3259  ftemp=(1.0 - gcon03*gcov03 - gcon13*gcov13 - gcon23*gcov23);
3260  igdetnosing=sign(ftemp)/(fabs(ftemp)+SMALL);
3261  pl=B3; MACP0A1(prim,i,j,k,pl) = (myBd3*gcon33 + Bu1*gcon03*gcov01 + Bu2*gcon03*gcov02 + Bu1*gcon13*gcov11 + Bu2*gcon13*gcov12 + Bu1*gcon23*gcov21 + Bu2*gcon23*gcov22)*igdetnosing;
3262 
3263  // old B_\phi imprecise copy (need to avoid singularity anywways)
3264  // pl=B3; MACP0A1(prim,i,j,k,pl) = MACP0A1(prim,ri,rj,rk,pl)*((ptrrgeom[pl]->gcov[GIND(3,3)])/(ptrgeom[pl]->gcov[GIND(3,3)]));
3265  }
3266  else{
3267  // maybe more consistent with interpolation
3268  // also seems to be less large values in general near poles for tilted spins
3269  pl=B3;
3270  igdetnosing=sign(ptrgeom[pl]->gdet)/(fabs(ptrgeom[pl]->gdet)+SMALL); // avoids 0.0 for any sign of ptrgeom->gdet
3271  MACP0A1(prim,i,j,k,pl) = (MACP0A1(prim,ri,rj,rk,pl)*(ptrrgeom[pl]->gdet) + dqgdetpl[pl]*(i-ri))*igdetnosing;
3272  }
3273 
3274 
3275 
3276 
3277  } // end LOOP over i and pl
3278 
3279 
3280 
3281 
3282  return(0);
3283 
3284 
3285 
3286 }
3287 
3288 
3289 
3290 
3291 #define UTHETAPOLEDEATH ((PRIMTOINTERP_3VELREL_GAMMAREL_DXDXP==VARTOINTERP)?(1):(0)) //only interpolate u^\theta instead of u^2 here when doing the same in interp
3292 
3293 #if(0 == UTHETAPOLEDEATH)
3294 # define MACP0A1mod(prim,ri,rj,rk,pl) MACP0A1(prim,ri,rj,rk,pl)
3295 #else
3296 
3297 static FTYPE MACP0A1mod(FTYPE (*prim)[NSTORE2][NSTORE3][NPR], int ri, int rj, int rk, int pl);
3298 
3299 static FTYPE MACP0A1mod(FTYPE (*prim)[NSTORE2][NSTORE3][NPR], int ri, int rj, int rk, int pl)
3300 {
3301  FTYPE dxdxp[NDIM][NDIM];
3302  int dir;
3303 
3304 #if(PRODUCTION==0)
3305  if ( dir < U2 || dir > B3 ) {
3306  dualfprintf( fail_file, "dir cannot be < U2 or > B3 in MACP0A1mod()\n" );
3307  exit(3232);
3308  }
3309 #endif
3310 
3311  dir = 1 + (pl-U1)%3; //direction that corresponds to pl (e.g., 2 for U2 or B2); assumes contiguous ordering of pl: <..>,U1,U2,U3,B1,B2,B3,<..>
3312 
3313  dxdxprim_ijk(ri, rj, rk, CENT, dxdxp);
3314 
3315  return( MACP0A1(prim,ri,rj,rk,pl) * dxdxp[dir][dir] );
3316 }
3317 #endif
3318 
3319 
3320 
3321 
3322 
3323 
3324 
3325 
3326 
3327 
3330 int poledeath(int whichx2,
3331  int boundstage, int finalstep, SFTYPE boundtime, int whichdir, int boundvartype, int *dirprim, int ispstag, FTYPE (*prim)[NSTORE2][NSTORE3][NPR],
3332  int *inboundloop,
3333  int *outboundloop,
3334  int *innormalloop,
3335  int *outnormalloop,
3336  int (*inoutlohi)[NUMUPDOWN][NDIM],
3337  int riin, int riout, int rjin, int rjout, int rkin, int rkout,
3338  int *dosetbc,
3339  int enerregion,
3340  int *localenerpos)
3341 {
3342 
3343 
3344 
3345 
3346  // when poledeath() called, already in parallel region
3347  //#pragma omp parallel // assume don't require EOS
3348  // {
3349 
3350  int jj;
3351  int i,j,k;
3352  int pl,pliter;
3353  int rim1,rip1;
3354  int ri,rj,rk;
3355  int rj0;
3356  int rjstag;
3357  int rjtest,rjstag0,deathstagjs0,deathstagje0,rjstagtest,deathstagjstest,deathstagjetest;
3358  int poleloc;
3359  int poleloccent;
3360  int deathjs0,deathje0;
3361  int deathjs,deathje;
3362  int deathstagjs,deathstagje;
3363  int gammadeathjs,gammadeathje;
3364  FTYPE X[NPR][NDIM];
3365  FTYPE V[NPR][NDIM];
3366  FTYPE Xr[NPR][NDIM];
3367  FTYPE X0[NDIM];
3368  int lowrho,lowuu,highu1;
3369  int deathjstest,deathjetest;
3370  FTYPE gamma,gammarj0,gammarjtest;
3371  FTYPE qsq,qsqrj0,qsqrjtest;
3372  FTYPE ftemp;
3373  FTYPE ucon[NDIM];
3374  FTYPE others[NUMOTHERSTATERESULTS];
3375  FTYPE uconrad[NDIM];
3376  FTYPE othersrad[NUMOTHERSTATERESULTS];
3377  FTYPE gammavaluelimit,gammaradvaluelimit;
3378  int doavginradius[NPR];
3379  int pl2;
3380  struct of_geom geomdontuse[NPR];
3381  struct of_geom *ptrgeom[NPR];
3382  struct of_geom rgeomdontuse[NPR];
3383  struct of_geom *ptrrgeom[NPR];
3384  FTYPE dxdxp[NDIM][NDIM];
3385  FTYPE prdiag[NPR],pr[NPR];
3386  int jstep;
3387  struct of_geom geompoledontuse;
3388  struct of_geom *ptrgeompole=&geompoledontuse;
3389 
3390 
3391  // NOTE: If bounding less than full quantities, this won't be good diag quantity.
3392  int testplexist[MAXNPR]={0};
3393  int candodiag=1;
3394  PBOUNDLOOP(pliter,pl) testplexist[pl]=1;
3395  PLOOP(pliter,pl){ // yes, should be PLOOP(pliter,pl)
3396  if(testplexist[pl]!=1){
3397  candodiag=0;
3398  }
3399  }
3400 
3401 
3402  // OPENMPMARK: Can't do this check because if reduce to 1 thread (even in OpenMP) then omp_in_parallel() is false!
3403  //#if(USEOPENMP)
3404  // // ensure within parallel region
3405  // if(!omp_in_parallel()){
3406  // dualfprintf(fail_file,"poledeath() called outside parallel region\n");
3407  // myexit(84968346);
3408  // }
3409  //#endif
3410 
3411 
3413  //
3414  // assign memory
3415  //
3417  PALLLOOP(pl){
3418  ptrgeom[pl]=&(geomdontuse[pl]);
3419  ptrrgeom[pl]=&(rgeomdontuse[pl]);
3420  }
3421 
3422 
3423 
3424 
3425 
3426 
3427 
3428 
3429 
3430 
3431 
3432 
3433 
3435  //
3436  // Loop over i,k
3437  //
3439 
3440  if(POLEDEATH){
3441 
3444 #pragma omp for schedule(OPENMPSCHEDULE(),OPENMPCHUNKSIZE(blocksize))
3447 
3448 
3449 
3451  //
3452  // setup loop ranges
3453  //
3455 
3456  int poledeathreal,polegammadeathreal;
3457  FTYPE Vtemp[NDIM];
3458  // note that doesn't matter the order of the j-loop since always using reference value (so for loop doesn't need change in <= to >=)
3459  if(whichx2==X2DN){
3460 
3461  //#define RADIUSOUTERDEATHMORE 300.0
3462 #define RADIUSOUTERDEATHMORE (BIG) //(OUTERDEATHRADIUS) // i.e. inactive
3463 
3464  bl_coord_ijk(i,0,k,CENT,Vtemp);
3465  if(Vtemp[1]>RADIUSOUTERDEATHMORE){
3466  poledeathreal=N2BND;
3467  polegammadeathreal=N2BND;
3468  }
3469  else{
3470  poledeathreal=POLEDEATH;
3471  polegammadeathreal=POLEGAMMADEATH;
3472  }
3473 
3474 
3475  jstep=-1; // direction of j loop, so starts with active cells in case modify boundary cells as dependent upon active cells.
3476 
3477  rj0 = poledeathreal;
3478  rjtest = rj0+DEATHEXPANDAMOUNT; // used to ensure near pole the density doesn't drop suddenly
3479  poleloc = 0;
3480  poleloccent = 0;
3481  // for poledeathreal==2, then deathjs,je=-2,-1,0,1 as required for CENT quantities rj=2
3482  deathjs0 = 0-poledeathreal;
3483  deathje0 = 0+poledeathreal-1;
3484 
3485  deathjstest = deathjs0-DEATHEXPANDAMOUNT;
3486  deathjetest = deathje0+DEATHEXPANDAMOUNT;
3487  if(deathjstest<inoutlohi[POINTDOWN][POINTDOWN][2]) deathjstest=inoutlohi[POINTDOWN][POINTDOWN][2];
3488  if(deathjetest<inoutlohi[POINTDOWN][POINTDOWN][2]) deathjetest=inoutlohi[POINTDOWN][POINTDOWN][2];
3489 
3490  // assume for poledeathreal==1 that B2 set correctly as 0 on pole and only do something if poledeathreal>1
3491  // if poledeathreal==2 then B2 set at -1,0,1 and will correctly set B2 to 0 at pole rj=2
3492  rjstag0 = rj0;
3493  deathstagjs0 = 0-poledeathreal+1;
3494  deathstagje0 = 0+poledeathreal-1;
3495 
3496  rjstagtest = rjtest;
3497  deathstagjstest = deathstagjs0-DEATHEXPANDAMOUNT;
3498  deathstagjetest = deathstagje0+DEATHEXPANDAMOUNT;
3499  if(deathstagjstest<inoutlohi[POINTDOWN][POINTDOWN][2]) deathstagjstest=inoutlohi[POINTDOWN][POINTDOWN][2];
3500  if(deathstagjetest<inoutlohi[POINTDOWN][POINTDOWN][2]) deathstagjetest=inoutlohi[POINTDOWN][POINTDOWN][2];
3501 
3502  // assumes velocity is always CENT
3503  gammadeathjs=0-polegammadeathreal;
3504  gammadeathje=0+polegammadeathreal-1;
3505  if(gammadeathjs<inoutlohi[POINTDOWN][POINTDOWN][2]) gammadeathjs=inoutlohi[POINTDOWN][POINTDOWN][2];
3506  if(gammadeathje<inoutlohi[POINTDOWN][POINTDOWN][2]) gammadeathje=inoutlohi[POINTDOWN][POINTDOWN][2];
3507 
3508  // NO, don't do below. Since poledeath called after MPI, need to set ghost cells as consistent with how active cells would have been set by other part of grid or other CPUs
3509  // if(special3dspc){
3510  // // then assume poledeath called *after* MPI (so have full and correct information across pole), so only should modify active cells and not boundary cells
3511  // deathjs0 = 0;
3512  // deathjstest = 0;
3513  // deathstagjs0 = 0;
3514  // deathstagjstest = 0;
3515  // gammadeathjs=0;
3516  // }
3517 
3518  }
3519  else if(whichx2==X2UP){
3520 
3521  bl_coord_ijk(i,N2BND,k,CENT,Vtemp);
3522  if(Vtemp[1]>RADIUSOUTERDEATHMORE){
3523  poledeathreal=N2BND;
3524  polegammadeathreal=N2BND;
3525  }
3526  else{
3527  poledeathreal=POLEDEATH;
3528  polegammadeathreal=POLEGAMMADEATH;
3529  }
3530 
3531  rj0=N2-1-poledeathreal;
3532  rjtest = rj0-DEATHEXPANDAMOUNT;
3533  poleloc=N2;
3534  poleloccent=N2-1;
3535  // if poledeathreal==2 then CENTs set at N2-2,N2-1,N2,N2+1 rj=N2-3
3536  deathjs0 = N2-1+1-poledeathreal;
3537  deathje0 = N2-1+poledeathreal;
3538 
3539  deathjstest = deathjs0-DEATHEXPANDAMOUNT;
3540  deathjetest = deathje0+DEATHEXPANDAMOUNT;
3541  if(deathjstest>inoutlohi[POINTUP][POINTUP][2]) deathjstest=inoutlohi[POINTUP][POINTUP][2];
3542  if(deathjetest>inoutlohi[POINTUP][POINTUP][2]) deathjetest=inoutlohi[POINTUP][POINTUP][2];
3543 
3544  // if poledeathreal==2, then B2 is set at N2-1,N2,N2+1 rj=N2-3
3545  if(dirprim[B2]==FACE2) rjstag0=N2-poledeathreal;
3546  else if(dirprim[B2]==CENT) rjstag0=rj0;
3547  deathstagjs0 = N2+1-poledeathreal;
3548  deathstagje0 = N2-1+poledeathreal;
3549 
3550  rjstagtest = rjtest;
3551  deathstagjstest = deathstagjs0-DEATHEXPANDAMOUNT;
3552  deathstagjetest = deathstagje0+DEATHEXPANDAMOUNT;
3553  if(deathstagjstest>inoutlohi[POINTUP][POINTUP][2]) deathstagjstest=inoutlohi[POINTUP][POINTUP][2];
3554  if(deathstagjetest>inoutlohi[POINTUP][POINTUP][2]) deathstagjetest=inoutlohi[POINTUP][POINTUP][2];
3555 
3556 
3557  // assumes velocity is always CENT . If POLEDEATH==2, N2-2,N2-1,N2,N2+1
3558  gammadeathjs = N2-1+1-polegammadeathreal;
3559  gammadeathje = N2-1+polegammadeathreal;
3560  if(gammadeathjs>inoutlohi[POINTUP][POINTUP][2]) gammadeathjs=inoutlohi[POINTUP][POINTUP][2];
3561  if(gammadeathje>inoutlohi[POINTUP][POINTUP][2]) gammadeathje=inoutlohi[POINTUP][POINTUP][2];
3562 
3563 
3564  // if(special3dspc){
3565  // // then assume poledeath called *after* MPI (so have full and correct information across pole), so only should modify active cells and not boundary cells
3566  // deathje0 = N2-1;
3567  // deathjetest = N2-1;
3568  // deathstagje0 = N2-1;
3569  // deathstagjetest = N2-1;
3570  // gammadeathje=N2-1;
3571  // }
3572 
3573  }
3574 
3576  FTYPE Rhorref=rhor_calc(0);
3577  if(ispstag==0){
3578  // even problems with poledeath for within horizon, can lead to Erf blowing up, so why 1|| below
3579  if( (1||startpos[1]+i>N1BND || Vtemp[1]>0.9*Rhorref) && (Vtemp[1]<OUTERDEATHRADIUS && OUTERDEATH==1 || OUTERDEATH==0)) continue;
3580  }
3581  }
3582 
3584  //
3585  // Skip this i,k if doing limited poledeath()
3586  //
3588 
3589  int dontskippoledeath;
3590  dontskippoledeath=0;
3591  if(IFLIMITPOLEDEATHIOUT>0){
3592  if(mycpupos[1]==ncpux1-1){
3593  if(i>N1-1-IFLIMITPOLEDEATHIOUT){
3594  // then force use of poledeath (no skip allowed)
3595  dontskippoledeath=1;
3596  //dualfprintf(fail_file,"OUTERDONTSKIP: %d %d %d\n",i,j,k);
3597  }
3598  }
3599  }
3600 
3601 
3602  if(IFLIMITPOLEDEATH>0 && ispstag==0 && dontskippoledeath==0){ // only go here if still possible to skip poledeath
3603  // assume only can skip if centered primitives (staggered are field, and nominally don't change them)
3604 
3605  FTYPE Vpole[NDIM];
3606 
3607  // use centered cell value (i.e. pl=RHO) to determine radius, and skip rest of this "j" if radius outside
3608  bl_coord_ijk(i,poleloc,k,FACE2,Vpole); // FACE2 at CENT for radius, but using FACE2 allows us to use j=poleloc
3609 
3610  // get u^\mu
3611  j=poleloccent;
3612  get_geometry(i, j, k, CENT, ptrgeompole);
3613  ucon_calc(MAC(prim,i,j,k),ptrgeompole,ucon, others);
3614 
3615 
3616  if(IFLIMITPOLEDEATH==1){
3617  if(ucon[RR]>=0.0){
3618  //dualfprintf(fail_file,"SKIP1: %d %d %d : %g %g\n",i,j,k,ucon[RR],Vpole[1]);
3619  continue;
3620  }
3621  }
3622  else if(IFLIMITPOLEDEATH==2){
3623  if(Vpole[1]>=RADIUSLIMITPOLEDEATHIN){
3624  //dualfprintf(fail_file,"SKIP2: %d %d %d : %g %g\n",i,j,k,ucon[RR],Vpole[1]);
3625  continue;
3626  }
3627  }
3628  else if(IFLIMITPOLEDEATH==3){
3629  if(ucon[RR]>=0.0 && Vpole[1]>=RADIUSLIMITPOLEDEATHIN){ // i.e., do poledeath unless *both* u^r>0 *and* beyond given radius
3630  //dualfprintf(fail_file,"SKIP3: %d %d %d : %g %g\n",i,j,k,ucon[RR],Vpole[1]);
3631  continue;
3632  }
3633  }
3634 
3635 
3636  //dualfprintf(fail_file,"NOTSKIP: %d %d %d : %g %g\n",i,j,k,ucon[RR],Vpole[1]);
3637 
3638  }
3639 
3640 
3641 
3642 
3643 
3644  // set reference positions in i,k
3645  ri = i;
3646  rk = k;
3647 
3648 
3650  //
3651  // set radial positions to average over
3652  //
3654  if(ri==inoutlohi[POINTDOWN][POINTDOWN][1]){
3655  // shift up, but not including ri
3656  rim1=ri+1;
3657  rip1=ri+2;
3658  }
3659  else if(ri==inoutlohi[POINTUP][POINTUP][1]){
3660  // shift down, but not including ri
3661  rim1=ri-2;
3662  rip1=ri-1;
3663  }
3664  else{
3665  // around ri
3666  rim1=ri-1;
3667  rip1=ri+1;
3668  }
3669 
3670 
3671 
3672 
3673 
3675  // first check for significant density drops near the axis or Lorentz factor spikes up near axis
3676  lowrho=lowuu=0;
3677  highu1=0;
3678 
3679 
3680  if(0){ // disabled for now
3681 
3682  get_geometry(ri, rj0, rk, dirprim[U1], ptrrgeom[U1]);
3683  gamma_calc(MAC(prim,ri,rj0,rk), ptrrgeom[U1],&gammarj0,&qsqrj0);
3684 
3685  get_geometry(ri, rjtest, rk, dirprim[U1], ptrrgeom[U1]);
3686  gamma_calc(MAC(prim,ri,rjtest,rk), ptrrgeom[U1],&gammarjtest,&qsqrjtest);
3687 
3688 
3689 
3690  for (j = deathjs0; j <= deathje0; j++) { // (no prims modified here, so no need for special j loop or diag_fixup)
3691 
3692  //POLEDENSITYDROPFACTOR
3693  //POLEGAMMAJUMPFACTOR
3694 
3695 
3696  pl=RHO;
3697  if(fabs(MACP0A1(prim,i,j,k,pl))< fabs(MACP0A1(prim,ri,rj0,rk,pl))/POLEDENSITYDROPFACTOR || fabs(MACP0A1(prim,i,j,k,pl))< fabs(MACP0A1(prim,ri,rjtest,rk,pl))/POLEDENSITYDROPFACTOR){
3698  lowrho=1;
3699  }
3700  pl=UU;
3701  if(fabs(MACP0A1(prim,i,j,k,pl))< fabs(MACP0A1(prim,ri,rj0,rk,pl)/POLEDENSITYDROPFACTOR) || fabs(MACP0A1(prim,i,j,k,pl))< fabs(MACP0A1(prim,ri,rjtest,rk,pl))/POLEDENSITYDROPFACTOR){
3702  lowuu=1;
3703  }
3704 
3705 
3706  // get Lorentz factor
3707  get_geometry(i, j, k, dirprim[U1], ptrrgeom[U1]);
3708  gamma_calc(MAC(prim,i,j,k), ptrrgeom[U1],&gamma,&qsq);
3709 
3710 
3711  // pl=U1;
3712  // if(fabs(MACP0A1(prim,i,j,k,pl))> fabs(MACP0A1(prim,ri,rj0,rk,pl))/POLEJUMPFACTOR || fabs(MACP0A1(prim,i,j,k,pl))> fabs(MACP0A1(prim,ri,rjtest,rk,pl))/POLEJUMPFACTOR ){
3713  // highu1=1;
3714  // }
3715 
3716  // check variation of Lorentz factor
3717  if(gamma> gammarj0*POLEGAMMAJUMPFACTOR || gamma> gammarjtest*POLEGAMMAJUMPFACTOR ){
3718  highu1=1;
3719  }
3720 
3721 
3722  }
3723 
3724  if(lowuu || lowrho || highu1){
3725  rj=rjtest; // expand copy procedure to use better reference value
3726  deathjs=deathjstest;
3727  deathje=deathjetest;
3728 
3729  rjstag=rjstagtest; // expand copy procedure to use better reference value
3730  deathstagjs=deathstagjstest;
3731  deathstagje=deathstagjetest;
3732  }
3733  else{
3734  // then use normal reference and range
3735  rj=rj0;
3736  deathjs=deathjs0;
3737  deathje=deathje0;
3738 
3739  rjstag=rjstag0;
3740  deathstagjs=deathstagjs0;
3741  deathstagje=deathstagje0;
3742  }
3743  }
3744  else{
3745  rj=rjtest; // expand copy procedure to use better reference value
3746  deathjs=deathjstest;
3747  deathje=deathjetest;
3748 
3749  rjstag=rjstagtest; // expand copy procedure to use better reference value
3750  deathstagjs=deathstagjstest;
3751  deathstagje=deathstagjetest;
3752  }
3753 
3754 
3755 
3756 
3758  //
3759  // reference location geometry and coordinates
3760  //
3762 
3763  PBOUNDLOOP(pliter,pl) {
3764  coord_ijk(ri,rj,rk,dirprim[pl],Xr[pl]); // reference locations for B2/U2
3765  if(pl==B2 && dirprim[B2]==FACE2){
3766  get_geometry(ri, rjstag, rk, dirprim[pl], ptrrgeom[pl]);
3767  }
3768  else{
3769  get_geometry(ri, rj, rk, dirprim[pl], ptrrgeom[pl]);
3770  }
3771  }
3772 
3773 
3774 
3775 
3776 
3777 
3778  //#define DEATHLOOP2 for (j = MIN(deathstagjs,deathjs); j <= MAX(deathstagje,deathje); j++) {
3779  // for (j = MIN(deathstagjs,deathjs); j <= MAX(deathstagje,deathje); j++) {
3780  // always iterate from active to ghost cell, so can modify active and have ready for ghost
3781  // didn't redefine js and je, they always go from small to large values, which is why ordering below
3782 #define DEATHLOOP2(whichx2) for (j = (whichx2==X2UP ? MIN(deathstagjs,deathjs) : MAX(deathstagje,deathje)) ; (whichx2==X2UP ? (j <= MAX(deathstagje,deathje)) : (j >= MIN(deathstagjs,deathjs)) ) ; (whichx2==X2UP ? j++ : j--) )
3783 
3785  //
3786  // For RHO,UU,U1,U3,B(1,2,3), and other primitives assumed to be like density
3787  //
3788  // Loop over points to be modified (prim's are only modifed HERE!)
3789  //
3791 
3792  DEATHLOOP2(whichx2){
3793 
3794 
3795 
3797  //
3798  // setup initial pr
3799  //
3801  //
3802  PBOUNDLOOP(pliter,pl) prdiag[pl]=MACP0A1(prim,i,j,k,pl);
3803  int madechange=0;
3804 
3805 
3806 
3807 
3808  PBOUNDLOOP(pliter,pl) {
3809  // pole location
3810  coord_ijk(i,poleloc,k,FACE2,X0); // pole locations for B2/U2
3811 
3812  // current location
3813  bl_coord_ijk_2(i,j,k,dirprim[pl],X[pl],V[pl]);
3814  get_geometry(i, j, k, dirprim[pl], ptrgeom[pl]);
3815 
3816  if(AVERAGEINRADIUS && (V[pl][RR]>RADIUSTOSTARTAVERAGING)){
3817  doavginradius[pl]=1;
3818  }
3819  else doavginradius[pl]=0;
3820 
3821  }// end pboundloop
3822 
3823 
3824 
3825 
3827  //
3828  // velocity and DENSITY (rho,u) POLEDATH
3829  //
3831  if(ispstag==0){
3832  // symmetric (if reflecting BC at pole) quantities
3833  if(j>=deathjs && j<=deathje){
3834 
3835 
3837  // u1, and u3, average in radius too!
3838  // copying this means copying \Omega_F in magnetically-dominated regime beyond LC
3839  PBOUNDLOOP(pliter,pl) {
3840  if(!(pl==U1 || pl==U3 || pl==URAD1 || pl==URAD3)) continue;
3841 
3842  FTYPE myvalue;
3843  myvalue=MACP0A1(prim,i,j,k,pl);
3844  FTYPE avgresult;
3845  if(doavginradius[pl]){
3846  avgresult=THIRD*(MACP0A1(prim,rim1,rj,rk,pl)+MACP0A1(prim,ri,rj,rk,pl)+MACP0A1(prim,rip1,rj,rk,pl));
3847  }
3848  else{
3849  avgresult=MACP0A1(prim,ri,rj,rk,pl);
3850  }
3851 
3852  if(1||fabs(avgresult)<myvalue){// only use reference if smaller value of velocity // well, for now, just copy over velocity under no additional conditions
3853  MACP0A1(prim,i,j,k,pl) = avgresult;
3854  madechange++;
3855  }
3856  }
3857 
3858 
3859 
3860  if(EOMTYPE!=EOMFFDE && EOMTYPE!=EOMFFDE2){
3862  // for densities
3863  // this helps remove drop-outs in density at high b^2/\rho_0 and high b^2/u
3864  PBOUNDLOOP(pliter,pl) {
3865  if(!(pl==RHO || pl==UU || pl==ENTROPY || pl==URAD0)) continue;
3866 
3867  FTYPE myvalue;
3868  myvalue=MACP0A1(prim,i,j,k,pl);
3869  FTYPE avgresult;
3870  if(doavginradius[pl]){
3871  avgresult=THIRD*(MACP0A1(prim,rim1,rj,rk,pl)+MACP0A1(prim,ri,rj,rk,pl)+MACP0A1(prim,rip1,rj,rk,pl));
3872  }
3873  else{
3874  avgresult=MACP0A1(prim,ri,rj,rk,pl);
3875  }
3876 
3877  if(1||fabs(avgresult)>myvalue){// only use reference if it is larger than polar value. That is, trying to avoid evacuation, and using smaller value of density only makes things worse, not better.
3878  MACP0A1(prim,i,j,k,pl) = avgresult;
3879  madechange++;
3880  }
3881 
3882  // if(doavginradius[pl]) MACP0A1(prim,i,j,k,pl) = THIRD*(MACP0A1(prim,rim1,rj,rk,pl)+MACP0A1(prim,ri,rj,rk,pl)+MACP0A1(prim,rip1,rj,rk,pl));
3883  // else MACP0A1(prim,i,j,k,pl) = MACP0A1(prim,ri,rj,rk,pl);
3884  // madechange++;
3885  }
3886 
3887 
3888  if(BCSIGMACONSTATPOLE==1){
3889 
3890  pl=RHO;
3891  FTYPE bsqr,sigmar;
3892 
3893  if (bsq_calc(MAC(prim,ri,rj,rk), ptrrgeom[pl], &bsqr) >= 1) FAILSTATEMENT("bounds.tools.c:poledeath()", "bsq_calc()", 1);
3894  sigmar=bsqr/(2.0*fabs(MACP0A1(prim,ri,rj,rk,pl)));
3895 
3896  // ensure constant sigma at pole
3897  FTYPE bsq,sigma;
3898  if (bsq_calc(MAC(prim,i,j,k), ptrgeom[pl], &bsq) >= 1) FAILSTATEMENT("bounds.tools.c:poledeath()", "bsq_calc()", 2);
3899  sigma=bsq/(2.0*fabs(MACP0A1(prim,i,j,k,pl)));
3900 
3901  if(sigma>BSQORHOLIMIT*0.5 && sigmar<BSQORHOLIMIT*0.5){
3902  MACP0A1(prim,i,j,k,pl) = fabs(MACP0A1(prim,i,j,k,pl)*(sigma/sigmar));
3903  madechange++;
3904  }
3905  // do nothing different than simple copy in any other cases. NOTE: If not high sigma, else if near-pole is low sigma, then this feeds in mass crazily
3906  }// end if BCSIGMACONSTATPOLE==1
3907  }// end if EOMTYPE!=EOMFFDE && EOMTYPE!=EOMFFDE2
3908 
3909 
3910  }// end if correct j
3911  }// end if ispstag==0
3912 
3913 
3914 
3915 
3917  //
3918  // symmetric (if reflecting BC at pole) quantities (e.g. B1)
3919  //
3921  if(j>=deathjs && j<=deathje){
3922  // B1 left alone
3923  }
3924 
3925 
3926 
3928  //
3929  // antisymmetric quantities (e.g. B2)
3930  //
3932  if(
3933  dirprim[B2]==FACE2 && j>=deathstagjs && j<=deathstagje ||
3934  dirprim[B2]==CENT && j>=deathjs && j<=deathje
3935  ){
3936 
3937  if(N2==1){
3939  // B2:
3940 #if(POLEINTERPTYPE==0)
3941  // if flow converges toward pole, then this loses information about the velocity and field approaching the pole
3942  // anti-symmetric (if reflecting BC at pole) quantities:
3943  MACP0A1(prim,i,j,k,B2) = 0.;
3944  madechange++;
3945 
3946 #elif(POLEINTERPTYPE==1 || POLEINTERPTYPE==2)
3947  // anti-symmetric (if reflecting BC at pole):
3948  // assume X[2] goes through 0 at the pole and isn't positive definite
3949  pl=B2;
3950  if(doavginradius[pl]) ftemp=THIRD*(MACP0A1(prim,rim1,rjstag,rk,pl) + MACP0A1(prim,ri,rjstag,rk,pl) + MACP0A1(prim,rip1,rjstag,rk,pl));
3951  else ftemp=MACP0A1(prim,ri,rjstag,rk,pl);
3952  MACP0A1(prim,i,j,k,pl) = ftemp + (X[pl][2]-Xr[pl][2])*(ftemp-0.0)/(Xr[pl][2]-X0[2]);
3953  madechange++;
3954 
3955 #elif(POLEINTERPTYPE==3)
3956  // then don't modify B2 -- trying to avoid instabilities related to divb=0 violation. And seems B2 behaves ok
3957 #endif
3958  }
3959  else{
3960  // then don't modify B2 to ensure divB=0
3961  }
3962 
3963  }
3964 
3965 
3966 
3967 
3968 
3970  //
3971  // other symmetric quantities (e.g. B3)
3972  //
3974  if(j>=deathjs && j<=deathje){
3975 
3977  // B3:
3978  pl=B3;
3979 
3980  if(N3==1){
3981 
3982 #if(POLEINTERPTYPE==0 || POLEINTERPTYPE==1)
3983  // symmetric:
3984  if(doavginradius[pl]) ftemp=THIRD*(MACP0A1(prim,rim1,rj,rk,pl) + MACP0A1(prim,ri,rj,rk,pl) + MACP0A1(prim,rip1,rj,rk,pl));
3985  else ftemp=MACP0A1(prim,ri,rj,rk,pl);
3986  MACP0A1(prim,i,j,k,B3) = ftemp;
3987  madechange++;
3988 
3989 #elif(POLEINTERPTYPE==2)
3990  // symmetric:
3991  // approximate B_\phi copy, which (unlike copying B3) can resolve singular currents on axis
3992  if(doavginradius[pl]) ftemp=THIRD*(MACP0A1(prim,rim1,rj,rk,pl) + MACP0A1(prim,ri,rj,rk,pl) + MACP0A1(prim,rip1,rj,rk,pl));
3993  else ftemp=MACP0A1(prim,ri,rj,rk,pl);
3994  MACP0A1(prim,i,j,k,pl) = ftemp*fabs((ptrrgeom[pl]->gcov[GIND(3,3)])/(ptrgeom[pl]->gcov[GIND(3,3)]));
3995  madechange++;
3996 
3997 #elif(POLEINTERPTYPE==3)
3998  // symmetric:
3999  if(doavginradius[pl]) ftemp=THIRD*(MACP0A1(prim,rim1,rj,rk,pl) + MACP0A1(prim,ri,rj,rk,pl) + MACP0A1(prim,rip1,rj,rk,pl));
4000  else ftemp=MACP0A1(prim,ri,rj,rk,pl);
4001  MACP0A1(prim,i,j,k,B3) = ftemp;
4002  madechange++;
4003 #endif
4004  }// end if N3==1
4005  else{
4006  // then do nothing if in 3D
4007  }
4008 
4009 
4010 
4011 
4013  //
4014  // do rest if any -- assumed at CENT
4015  //
4017  if(ispstag==0){
4018  PBOUNDLOOP(pliter,pl){
4019  if((pl==RHO || pl==UU || pl==U1 || pl==U2 || pl==U3 || pl==ENTROPY || pl==B1 || pl==B2 || pl==B3 || pl==URAD0 || pl==URAD1 || pl==URAD2 || pl==URAD3)) continue;
4020 
4021  if(doavginradius[pl]) ftemp=THIRD*(MACP0A1(prim,rim1,rj,rk,pl) + MACP0A1(prim,ri,rj,rk,pl) + MACP0A1(prim,rip1,rj,rk,pl));
4022  else ftemp=MACP0A1(prim,ri,rj,rk,pl);
4023  MACP0A1(prim,i,j,k,pl)=ftemp;
4024  madechange++;
4025  }
4026  }
4027 
4028  }// end over CENT-type quantities
4029 
4030 
4031 
4032  if(madechange && candodiag){
4034  //
4035  // accounting for on-grid changes
4036  //
4038  int modcons=(DOENOFLUX != NOENOFLUX);
4039  FTYPE *ucons;
4040  ucons=GLOBALMAC(unewglobal,i,j,k); // GODMARK -- need to pass ucons to bounds if really want !=NOENOFLUX method to work
4041  // GODMARK: assume all quantities at the same location since only ispstag==0 modifies relevant primitves, so ptrgeom[pl]->ptrgoem
4042  // in general, not sure which pl really exists at this point, so pick first in PBOUNDLOOP loop
4043  struct of_geom *fixupptrgeom;
4044  PBOUNDLOOP(pliter,pl) {
4045  fixupptrgeom=ptrgeom[pl];
4046  break;
4047  }
4048  PBOUNDLOOP(pliter,pl) pr[pl]=MACP0A1(prim,i,j,k,pl);
4049  int doingmhdfixup=1;
4050  diag_fixup(modcons,prdiag,pr,ucons,fixupptrgeom,finalstep,doingmhdfixup,COUNTBOUND1);
4051  PBOUNDLOOP(pliter,pl) prdiag[pl]=pr[pl];
4052  }
4053 
4054 
4055  }// end loop over j
4056 
4057 
4058 
4059 
4060 
4061 
4063  //
4064  // FOR U2 alone!
4065  //
4066  // Loop over points to be modified (prim's are only modifed HERE!)
4067  //
4068  // independent loop over j from above density loop because U2 changes may required RHO and those must already all be changed for all j
4069  //
4071 
4072  DEATHLOOP2(whichx2){
4073 
4074 
4075 
4077  //
4078  // setup initial pr
4079  //
4081  PBOUNDLOOP(pliter,pl) prdiag[pl]=MACP0A1(prim,i,j,k,pl);
4082  int madechange=0;
4083 
4084 
4085 
4086 
4087  PBOUNDLOOP(pliter,pl) {
4088  // pole location
4089  coord_ijk(i,poleloc,k,FACE2,X0); // pole locations for B2/U2
4090 
4091  // current location
4092  bl_coord_ijk_2(i,j,k,dirprim[pl],X[pl],V[pl]);
4093  get_geometry(i, j, k, dirprim[pl], ptrgeom[pl]);
4094 
4095  if(AVERAGEINRADIUS && (V[pl][RR]>RADIUSTOSTARTAVERAGING)){
4096  doavginradius[pl]=1;
4097  }
4098  else doavginradius[pl]=0;
4099 
4100  }
4101 
4102 
4103 
4104  if(ispstag==0){
4105 
4106  if(j>=deathjs && j<=deathje){
4108  int iteru2;
4109  int plrho;
4110  for(iteru2=0;iteru2<=1;iteru2++){
4111  // U2:
4112  if(iteru2==0){
4113  pl=U2;
4114  plrho=RHO;
4115  }
4116  else if(iteru2==1){
4117  pl=URAD2;
4118  plrho=URAD0;
4119  }
4120  if(URAD0<0) continue; // skip of not doing radiation
4121 
4122 
4123 
4124  if(special3dspc){
4125  // U2 not necessarily anti-symmetric in this case.
4126  // This will be kinda odd if POLEDEATH>1 due to comparing non-local regions.
4127  // But, for now, POLEDEATH<=1 is set so make sense.
4128 
4129  int jother;
4130  FTYPE signD;
4131  if(j<N2/2){
4132  jother=-1-j;
4133  }
4134  else{
4135  jother=N2-1+(N2-j);
4136  }
4137 
4138  if(j>=jother) signD=+1.0;
4139  else signD=-1.0;
4140 
4141 
4142  // see if sucking on pole
4143  FTYPE rhovjhere,rhovjother,rhovDiff;
4144  rhovjhere =MACP0A1(prim,i,j,k,plrho)*MACP0A1mod(prim,i,j,k,pl);
4145  rhovjother=MACP0A1(prim,i,jother,k,plrho)*MACP0A1mod(prim,i,jother,k,pl);
4146 
4147  rhovDiff=signD*(rhovjhere-rhovjother);
4148  // same gdet, so no need to multiply both by same factor for below test
4149  // make change to both simultaneously so that when other j is hit, rhovDiff condition is no longer hit and all is consistent as if separate memory field for entire poledeath before final copy-over to primitive memory space.
4150  // But do active grid cells first (determined by DEATHLOOPJ) so diag_fixup() occurs on active region for accounting.
4151  if(rhovDiff>0.0){
4152  // then sucking on pole
4153  // must change active and ghost cells consistently for any number of CPUs
4154  // so average-out the suck to zero suck by changing the values equally (as weighted by mass)
4155  FTYPE dU2 =-signD*rhovDiff*0.5/MACP0A1(prim,i,j,k,plrho);
4156  FTYPE U2jhere = MACP0A1mod(prim,i,j,k,pl);
4157  FTYPE U2jother= MACP0A1mod(prim,i,jother,k,pl);
4158 
4159  if( (fabs(U2jhere)>fabs(dU2))&&(fabs(U2jother)>fabs(dU2)) ){
4160  // only change if in both cases we lower the velocity, not increase.
4161  MACP0A1(prim,i,j,k,pl) += dU2;
4162  MACP0A1(prim,i,jother,k,pl) -= dU2;
4163  madechange++;
4164  }
4165  else{ // then jhere or jother is changed by an absolute magnitude more than its value, which we want to avoid
4166  // crush a bit only as much as leaves smaller value changed as much as possible without increasing its magnitude
4167  if(fabs(U2jhere)<fabs(U2jother)){
4168  // then drop jhere as closest to fixed value (100% change), and reset jother with same value
4169  MACP0A1(prim,i,j,k,pl) *= -1.0;
4170  MACP0A1(prim,i,jother,k,pl) = MACP0A1(prim,i,j,k,pl);
4171  madechange++;
4172  }
4173  else{
4174  // then drops down to value matching other side so D=0 in the end still, so still no sucking.
4175  MACP0A1(prim,i,jother,k,pl) *= -1.0;
4176  MACP0A1(prim,i,j,k,pl) = MACP0A1(prim,i,jother,k,pl);
4177  madechange++;
4178  }
4179  }// end else abs mag of change is larger than 100% for one of the values
4180 
4181  // MACP0A1(prim,i,j,k,pl) =0.0;
4182  // madechange++;
4183  // MACP0A1(prim,i,jother,k,pl) += +rhovDiff*0.5/MACP0A1(prim,i,jother,k,plrho); // this taken care of by other j in ghost region
4184 
4185  // Note that for anti-symmetric U2 (i.e. reflective BCs around pole) this is same as crushing regularization leading to pl->0
4186  }
4187  else{
4188  // then just enforce linear behavior near pole
4189  // NOT YET
4190  }
4191 
4192 
4193  }
4194  else{
4195 
4196 
4197 #if(POLEINTERPTYPE==0)
4198  // if flow converges toward pole, then this loses information about the velocity and field approaching the pole
4199  // anti-symmetric (if reflecting BC at pole) quantities:
4200  MACP0A1(prim,i,j,k,pl) = 0.;
4201  madechange++;
4202 
4203 #elif(POLEINTERPTYPE==1 || POLEINTERPTYPE==2)
4204  // anti-symmetric (if reflecting BC at pole):
4205  // assume X[2] goes through 0 at the pole and isn't positive definite
4206  if(doavginradius[pl]) ftemp=THIRD*(MACP0A1mod(prim,rim1,rj,rk,pl) + MACP0A1mod(prim,ri,rj,rk,pl) + MACP0A1mod(prim,rip1,rj,rk,pl));
4207  else ftemp=MACP0A1mod(prim,ri,rj,rk,pl);
4208  MACP0A1(prim,i,j,k,pl) = ftemp + (X[pl][2]-Xr[pl][2])*(ftemp-0.0)/(Xr[pl][2]-X0[2]);
4209  madechange++;
4210 
4211 #elif(POLEINTERPTYPE==3)
4212 
4213  // anti-symmetric (if reflecting BC at pole):
4214 
4215  // assume X[2] goes through 0 at the pole and isn't positive definite
4216 
4217  // choose reference value
4218  // ftemp=THIRD*(MACP0A1mod(prim,rim1,rj,rk,pl) + MACP0A1mod(prim,ri,rj,rk,pl) + MACP0A1mod(prim,rip1,rj,rk,pl));
4219  ftemp=MACP0A1mod(prim,ri,rj,rk,pl);
4220 
4221  if(whichx2==X2DN && ftemp>0.0){
4222  // then sucking on \theta=0 pole
4223  // try to minimize sucking on pole by finding minimum U2 around
4224  for(jj=0;jj<=rj+DEATHEXPANDAMOUNT;jj++) ftemp=MIN(ftemp,MACP0A1mod(prim,ri,jj,rk,pl));
4225  if(doavginradius[pl]){
4226  for(jj=0;jj<=rj+DEATHEXPANDAMOUNT;jj++) ftemp=MIN(ftemp,MACP0A1mod(prim,rip1,jj,rk,pl));
4227  for(jj=0;jj<=rj+DEATHEXPANDAMOUNT;jj++) ftemp=MIN(ftemp,MACP0A1mod(prim,rim1,jj,rk,pl));
4228  }
4229 
4230  ftemp=0.0; // try crushing sucking GODMARK
4231 
4232  // assume ftemp is at reference location
4233  MACP0A1(prim,i,j,k,pl) = ftemp + (X[pl][2]-Xr[pl][2])*(ftemp-0.0)/(Xr[pl][2]-X0[2]);
4234  madechange++;
4235  }
4236  else if(whichx2==X2UP && ftemp<0.0){
4237  // then sucking on \theta=\pi pole
4238  for(jj=N2-1;jj>=rj-DEATHEXPANDAMOUNT;jj--) ftemp=MAX(ftemp,MACP0A1mod(prim,ri,jj,rk,pl));
4239  if(doavginradius[pl]){
4240  for(jj=N2-1;jj>=rj-DEATHEXPANDAMOUNT;jj--) ftemp=MAX(ftemp,MACP0A1mod(prim,rip1,jj,rk,pl));
4241  for(jj=N2-1;jj>=rj-DEATHEXPANDAMOUNT;jj--) ftemp=MAX(ftemp,MACP0A1mod(prim,rim1,jj,rk,pl));
4242  }
4243 
4244  ftemp=0.0; // try crushing sucking GODMARK
4245 
4246  // assume ftemp is at reference location (same formula for both poles)
4247  MACP0A1(prim,i,j,k,pl) = ftemp + (X[pl][2]-Xr[pl][2])*(ftemp-0.0)/(Xr[pl][2]-X0[2]);
4248  madechange++;
4249  }
4250  else{
4251  // otherwise enforce natural regular linear behavior on U2
4252  if(doavginradius[pl]) ftemp=THIRD*(MACP0A1mod(prim,rim1,rj,rk,pl) + MACP0A1mod(prim,ri,rj,rk,pl) + MACP0A1mod(prim,rip1,rj,rk,pl));
4253  else ftemp=MACP0A1mod(prim,ri,rj,rk,pl);
4254 
4255  MACP0A1(prim,i,j,k,pl) = ftemp + (X[pl][2]-Xr[pl][2])*(ftemp-0.0)/(Xr[pl][2]-X0[2]);
4256  madechange++;
4257  }
4258 
4259 #endif // endif POLEINTERPTYPE==3
4260  }// end if special3dspc==0
4261 
4262 
4263 #if( UTHETAPOLEDEATH )
4264  if(pl==U2){
4265  //if interpolated u^\theta, now convert back to u^2
4266  dxdxprim_ijk(i, j, k, CENT, dxdxp);
4267  MACP0A1(prim,i,j,k,pl) /= dxdxp[1 + (pl-U1)%3][1 + (pl-U1)%3];
4268  madechange++;
4269  }
4270  if(pl==URAD2){
4271  //if interpolated u^\theta, now convert back to u^2
4272  dxdxprim_ijk(i, j, k, CENT, dxdxp);
4273  MACP0A1(prim,i,j,k,pl) /= dxdxp[1 + (pl-URAD1)%3][1 + (pl-URAD1)%3];
4274  madechange++;
4275  }
4276 #endif
4277 
4278  }// over U2 and URAD2
4279  }// end if correct j range
4280  }// end if ispstag==0
4281 
4282 
4283 
4284  if(madechange&&candodiag){ // only need accounting if actually made a change
4286  //
4287  // accounting for on-grid changes
4288  //
4290  int modcons=(DOENOFLUX != NOENOFLUX);
4291  FTYPE *ucons;
4292  ucons=GLOBALMAC(unewglobal,i,j,k); // GODMARK -- need to pass ucons to bounds if really want !=NOENOFLUX method to work
4293  // GODMARK: assume all quantities at the same location since only ispstag==0 modifies relevant primitves, so ptrgeom[pl]->ptrgoem
4294  // in general, not sure which pl really exists at this point, so pick first in PBOUNDLOOP loop
4295  struct of_geom *fixupptrgeom;
4296  PBOUNDLOOP(pliter,pl) {
4297  fixupptrgeom=ptrgeom[pl];
4298  break;
4299  }
4300 
4301  PBOUNDLOOP(pliter,pl) pr[pl]=MACP0A1(prim,i,j,k,pl);
4302  int doingmhdfixup=1;
4303  diag_fixup(modcons,prdiag,pr,ucons,fixupptrgeom,finalstep,doingmhdfixup,COUNTBOUND1);
4304  PBOUNDLOOP(pliter,pl) prdiag[pl]=pr[pl];
4305  }
4306 
4307 
4308  }// end loop over j
4309  }// end loop 13
4310  }// end if POLEDEATH
4311 
4312 
4313 
4314 
4315 
4316 
4317 
4318 
4319 
4320  if(POLEGAMMADEATH){ // should come after POLEDEATH to be useful
4321  if(ispstag==0){
4322  // fixup
4323 
4324 
4325  // assume only dealing with velocities at same location (loop assumes CENT)
4326  pl=U1;
4327 
4328 
4331 #pragma omp for schedule(OPENMPSCHEDULE(),OPENMPCHUNKSIZE(blocksize))
4334 
4335 
4337  //
4338  // setup loop ranges
4339  //
4341 
4342 
4344  // BELOW JUST COPY OF ABOVE
4345  int poledeathreal,polegammadeathreal;
4346  FTYPE Vtemp[NDIM];
4347  // note that doesn't matter the order of the j-loop since always using reference value (so for loop doesn't need change in <= to >=)
4348  if(whichx2==X2DN){
4349 
4350  bl_coord_ijk(i,0,k,CENT,Vtemp);
4351  if(Vtemp[1]>RADIUSOUTERDEATHMORE){
4352  poledeathreal=N2BND;
4353  polegammadeathreal=N2BND;
4354  }
4355  else{
4356  poledeathreal=POLEDEATH;
4357  polegammadeathreal=POLEGAMMADEATH;
4358  }
4359 
4360 
4361  jstep=-1; // direction of j loop, so starts with active cells in case modify boundary cells as dependent upon active cells.
4362 
4363  rj0 = poledeathreal;
4364  rjtest = rj0+DEATHEXPANDAMOUNT; // used to ensure near pole the density doesn't drop suddenly
4365  poleloc = 0;
4366  poleloccent = 0;
4367  // for poledeathreal==2, then deathjs,je=-2,-1,0,1 as required for CENT quantities rj=2
4368  deathjs0 = 0-poledeathreal;
4369  deathje0 = 0+poledeathreal-1;
4370 
4371  deathjstest = deathjs0-DEATHEXPANDAMOUNT;
4372  deathjetest = deathje0+DEATHEXPANDAMOUNT;
4373  if(deathjstest<inoutlohi[POINTDOWN][POINTDOWN][2]) deathjstest=inoutlohi[POINTDOWN][POINTDOWN][2];
4374  if(deathjetest<inoutlohi[POINTDOWN][POINTDOWN][2]) deathjetest=inoutlohi[POINTDOWN][POINTDOWN][2];
4375 
4376  // assume for poledeathreal==1 that B2 set correctly as 0 on pole and only do something if poledeathreal>1
4377  // if poledeathreal==2 then B2 set at -1,0,1 and will correctly set B2 to 0 at pole rj=2
4378  rjstag0 = rj0;
4379  deathstagjs0 = 0-poledeathreal+1;
4380  deathstagje0 = 0+poledeathreal-1;
4381 
4382  rjstagtest = rjtest;
4383  deathstagjstest = deathstagjs0-DEATHEXPANDAMOUNT;
4384  deathstagjetest = deathstagje0+DEATHEXPANDAMOUNT;
4385  if(deathstagjstest<inoutlohi[POINTDOWN][POINTDOWN][2]) deathstagjstest=inoutlohi[POINTDOWN][POINTDOWN][2];
4386  if(deathstagjetest<inoutlohi[POINTDOWN][POINTDOWN][2]) deathstagjetest=inoutlohi[POINTDOWN][POINTDOWN][2];
4387 
4388  // assumes velocity is always CENT
4389  gammadeathjs=0-polegammadeathreal;
4390  gammadeathje=0+polegammadeathreal-1;
4391  if(gammadeathjs<inoutlohi[POINTDOWN][POINTDOWN][2]) gammadeathjs=inoutlohi[POINTDOWN][POINTDOWN][2];
4392  if(gammadeathje<inoutlohi[POINTDOWN][POINTDOWN][2]) gammadeathje=inoutlohi[POINTDOWN][POINTDOWN][2];
4393 
4394  // NO, don't do below. Since poledeath called after MPI, need to set ghost cells as consistent with how active cells would have been set by other part of grid or other CPUs
4395  // if(special3dspc){
4396  // // then assume poledeath called *after* MPI (so have full and correct information across pole), so only should modify active cells and not boundary cells
4397  // deathjs0 = 0;
4398  // deathjstest = 0;
4399  // deathstagjs0 = 0;
4400  // deathstagjstest = 0;
4401  // gammadeathjs=0;
4402  // }
4404 
4405  }
4406  else if(whichx2==X2UP){
4407 
4408 
4409  bl_coord_ijk(i,N2BND,k,CENT,Vtemp);
4410  if(Vtemp[1]>RADIUSMOREDEATH){
4411  poledeathreal=N2BND;
4412  polegammadeathreal=N2BND;
4413  }
4414  else{
4415  poledeathreal=POLEDEATH;
4416  polegammadeathreal=POLEGAMMADEATH;
4417  }
4418 
4419  rj0=N2-1-poledeathreal;
4420  rjtest = rj0-DEATHEXPANDAMOUNT;
4421  poleloc=N2;
4422  poleloccent=N2-1;
4423  // if poledeathreal==2 then CENTs set at N2-2,N2-1,N2,N2+1 rj=N2-3
4424  deathjs0 = N2-1+1-poledeathreal;
4425  deathje0 = N2-1+poledeathreal;
4426 
4427  deathjstest = deathjs0-DEATHEXPANDAMOUNT;
4428  deathjetest = deathje0+DEATHEXPANDAMOUNT;
4429  if(deathjstest>inoutlohi[POINTUP][POINTUP][2]) deathjstest=inoutlohi[POINTUP][POINTUP][2];
4430  if(deathjetest>inoutlohi[POINTUP][POINTUP][2]) deathjetest=inoutlohi[POINTUP][POINTUP][2];
4431 
4432  // if poledeathreal==2, then B2 is set at N2-1,N2,N2+1 rj=N2-3
4433  if(dirprim[B2]==FACE2) rjstag0=N2-poledeathreal;
4434  else if(dirprim[B2]==CENT) rjstag0=rj0;
4435  deathstagjs0 = N2+1-poledeathreal;
4436  deathstagje0 = N2-1+poledeathreal;
4437 
4438  rjstagtest = rjtest;
4439  deathstagjstest = deathstagjs0-DEATHEXPANDAMOUNT;
4440  deathstagjetest = deathstagje0+DEATHEXPANDAMOUNT;
4441  if(deathstagjstest>inoutlohi[POINTUP][POINTUP][2]) deathstagjstest=inoutlohi[POINTUP][POINTUP][2];
4442  if(deathstagjetest>inoutlohi[POINTUP][POINTUP][2]) deathstagjetest=inoutlohi[POINTUP][POINTUP][2];
4443 
4444 
4445  // assumes velocity is always CENT . If POLEDEATH==2, N2-2,N2-1,N2,N2+1
4446  gammadeathjs = N2-1+1-polegammadeathreal;
4447  gammadeathje = N2-1+polegammadeathreal;
4448  if(gammadeathjs>inoutlohi[POINTUP][POINTUP][2]) gammadeathjs=inoutlohi[POINTUP][POINTUP][2];
4449  if(gammadeathje>inoutlohi[POINTUP][POINTUP][2]) gammadeathje=inoutlohi[POINTUP][POINTUP][2];
4450 
4451 
4452  // if(special3dspc){
4453  // // then assume poledeath called *after* MPI (so have full and correct information across pole), so only should modify active cells and not boundary cells
4454  // deathje0 = N2-1;
4455  // deathjetest = N2-1;
4456  // deathstagje0 = N2-1;
4457  // deathstagjetest = N2-1;
4458  // gammadeathje=N2-1;
4459  // }
4460 
4461  }
4463 
4464 
4466  FTYPE Rhorref=rhor_calc(0);
4467  if(ispstag==0){
4468  // even problems with poledeath for within horizon, can lead to Erf blowing up, so why 1|| below
4469  if( (1||startpos[1]+i>N1BND || Vtemp[1]>0.9*Rhorref) && (Vtemp[1]<OUTERDEATHRADIUS && OUTERDEATH==1 || OUTERDEATH==0)) continue;
4470  }
4471  }
4472 
4473 
4474 
4475  for (j = gammadeathjs; j <= gammadeathje; j++) { // currently not multiple-point dependent, so normal j loop is fine
4476 
4477 
4479  //
4480  // setup initial pr (already using pl, so use pl2 for prdiag assignment/setup)
4481  //
4483  PBOUNDLOOP(pliter,pl2) prdiag[pl2]=MACP0A1(prim,i,j,k,pl2);
4484  int madechange=0;
4485 
4486 
4488  //
4489  // forever-kept debug
4490  //
4492  PBOUNDLOOP(pliter,pl2){
4493  if( !isfinite( MACP0A1(prim,i,j,k,pl2))){
4494  dualfprintf(fail_file,"BNDNAN: ispstag=%d t=%21.15g nstep=%ld steppart=%d :: i=%d j=%d k=%d pl2=%d prim=%21.15g\n",ispstag,t,nstep,steppart,i,j,k,pl2,MACP0A1(prim,i,j,k,pl2));
4495  }
4496  }
4497 
4498 
4499 
4501  //
4502  // get u^r and limit if necessary
4503  //
4505  get_geometry(i, j, k, dirprim[pl], ptrgeom[pl]);
4506 
4507  ucon_calc(&MACP0A1(prim,i,j,k,0),ptrgeom[pl],ucon, others);
4508  if(URAD0>=0){
4509  ucon_calc(&MACP0A1(prim,i,j,k,URAD1-U1),ptrgeom[pl],uconrad, othersrad);
4510  }
4511  // only limit velocity if outgoing relative to grid (GODMARK: only valid in KS or BL-like coordinates such that u^r>0 means outgoing w.r.t. an observer at infinity)
4512  // if(ucon[RR]>0.0||uconrad[RR]>0.0){
4513 
4514  // get V[RR]
4515  bl_coord_ijk_2(i,j,k,dirprim[pl],X[pl],V[pl]);
4516 
4517 
4518  // check if did limit
4519  int didlimit=0;
4520 
4522  //
4523  // FLUID
4524  //
4526  gammavaluelimit = GAMMAMAX; // default is no limit
4527 
4528  if(V[pl][RR]>GAMMAPOLEOUTGOINGRADIUSIN){
4529 
4530  if(V[pl][RR]<=GAMMAPOLEOUTGOINGRADIUS){
4531  // flat \gamma limit up to GAMMAPOLEOUTGOINGRADIUS
4532  if(ucon[RR]>0.0){
4533  didlimit=1;
4534  gammavaluelimit = GAMMAPOLEOUTGOING;
4535  }
4536  else if(ucon[RR]<0.0){
4537  didlimit=1;
4538  gammavaluelimit = GAMMAPOLEINGOINGOUT;
4539  }
4540  else gammavaluelimit = GAMMAMAX;
4541 
4542  }
4543  else{
4544  if(1||ucon[RR]>0.0){ // both ingoing and outgoing components limited in same way
4545  didlimit=1;
4547  gammavaluelimit=MIN(gammavaluelimit,GAMMAPOLEOUTGOINGMAX);
4548  }
4549  else gammavaluelimit=GAMMAMAX;
4550 
4551  }
4552 
4553  }// end if u^r>0
4554  else{
4555  didlimit=1;
4556  gammavaluelimit=GAMMAPOLEINGOING;
4557  }
4558 
4559 
4560 
4562  //
4563  // RADIATION
4564  //
4566  gammaradvaluelimit = GAMMAMAXRAD; // default is no limit
4567 
4568  if(URAD0>=0){
4569  if(V[pl][RR]>GAMMARADPOLEOUTGOINGRADIUSIN){
4570 
4571  if(V[pl][RR]<=GAMMARADPOLEOUTGOINGRADIUS){
4572 
4573  if(uconrad[RR]>0.0){
4574  didlimit=1;
4575  gammaradvaluelimit = GAMMARADPOLEOUTGOING;
4576  }
4577  else if(uconrad[RR]<0.0){
4578  didlimit=1;
4579  gammaradvaluelimit = GAMMARADPOLEINGOINGOUT;
4580  }
4581  else gammaradvaluelimit = GAMMAMAXRAD;
4582  }
4583  else{
4584  if(1||uconrad[RR]>0.0){ // both ingoing and outgoing components limited in same way
4585  didlimit=1;
4587  gammaradvaluelimit=MIN(gammaradvaluelimit,GAMMARADPOLEOUTGOINGMAX);
4588  }
4589  else gammaradvaluelimit=GAMMAMAXRAD;
4590  }
4591  }// end if u^r>0
4592  else{
4593  didlimit=1;
4594  gammaradvaluelimit=GAMMARADPOLEINGOING;
4595  }
4596  }
4597 
4598 
4600  //
4601  // Limit gamma's
4602  //
4604  if(didlimit){
4605  //if(j==0 && (ucon[TT]>10 || uconrad[TT]>10)){ // DEBUG
4606  // dualfprintf(fail_file,"GODHERE1: ur=%g %g : ut=%g %g : limit= %g %g\n",ucon[RR]*sqrt(ptrgeom[pl]->gcov[GIND(RR,RR)]),uconrad[RR]*sqrt(ptrgeom[pl]->gcov[GIND(RR,RR)]),ucon[TT],uconrad[TT],gammavaluelimit,gammaradvaluelimit);
4607  //}
4608 
4609  limit_gamma(0,gammavaluelimit,gammaradvaluelimit,MAC(prim,i,j,k),NULL,ptrgeom[pl],-1);
4610 
4611  // if(j==0 && (ucon[TT]>10 || uconrad[TT]>10)){ // DEBUG
4612  // ucon_calc(&MACP0A1(prim,i,j,k,0),ptrgeom[pl],ucon, others);
4613  // if(URAD0>=0){
4614  // ucon_calc(&MACP0A1(prim,i,j,k,URAD1-U1),ptrgeom[pl],uconrad, othersrad);
4615  // }
4616  // dualfprintf(fail_file,"GODHERE2: ur=%g %g : ut=%g %g : limit= %g %g\n",ucon[RR]*sqrt(ptrgeom[pl]->gcov[GIND(RR,RR)]),uconrad[RR]*sqrt(ptrgeom[pl]->gcov[GIND(RR,RR)]),ucon[TT],uconrad[TT],gammavaluelimit,gammaradvaluelimit);
4617  // }
4618  madechange++;
4619  }
4620 
4621 
4622 
4624  //
4625  // Accounting for gamma changes
4626  //
4628  if(madechange&&candodiag){
4630  //
4631  // accounting for on-grid changes
4632  //
4634  int modcons=(DOENOFLUX != NOENOFLUX);
4635  FTYPE *ucons;
4636  ucons=GLOBALMAC(unewglobal,i,j,k); // GODMARK -- need to pass ucons to bounds if really want !=NOENOFLUX method to work
4637  // GODMARK: assume all quantities at the same location since ispstag==0 is assumed in this section, so ptrgeom[pl]->ptrgoem
4638  // in general, not sure which pl2 really exists at this point, so pick first in PBOUNDLOOP loop
4639  struct of_geom *fixupptrgeom;
4640  // PBOUNDLOOP(pliter,pl2) {
4641  pl2=pl; // only 1 ptrgeom defined
4642  fixupptrgeom=ptrgeom[pl2];
4643  // break;
4644  // }
4645 
4646  PBOUNDLOOP(pliter,pl2) pr[pl2]=MACP0A1(prim,i,j,k,pl2);
4647  int doingmhdfixup=1;
4648  diag_fixup(modcons,prdiag,pr,ucons,fixupptrgeom,finalstep,doingmhdfixup,COUNTBOUND2);
4649  PBOUNDLOOP(pliter,pl2) prdiag[pl2]=pr[pl2];
4650 
4651  }
4652 
4653 
4654  }// end over j
4655  }// end loop 13
4656  }// end if ispstag==0
4657  }// end if POLEGAMMADEATH
4658 
4659 
4660 
4661  return(0);
4662 
4663 
4664 }
4665 
4666 
4667 
4668 
4669 
4670 
4671 
4672 #define DEBUGPOLESMOOTH 0
4673 
4677 int polesmooth(int whichx2,
4678  int boundstage, int finalstep, SFTYPE boundtime, int whichdir, int boundvartype, int *dirprim, int ispstag, FTYPE (*prim)[NSTORE2][NSTORE3][NPR],
4679  int *inboundloop,
4680  int *outboundloop,
4681  int *innormalloop,
4682  int *outnormalloop,
4683  int (*inoutlohi)[NUMUPDOWN][NDIM],
4684  int riin, int riout, int rjin, int rjout, int rkin, int rkout,
4685  int *dosetbc,
4686  int enerregion,
4687  int *localenerpos)
4688 {
4689 
4690  int i, j, k;
4691  int j0, dj, stopj, rj;
4692  int fakej;
4693  FTYPE X[NDIM], V[NDIM];
4694  FTYPE r, th, ph;
4695  FTYPE dxdxp[NDIM][NDIM], idxdxp[NDIM][NDIM];
4696  FTYPE *pr;
4697  static FTYPE *fullpr;
4698  FTYPE cartavgpr[NPR], spcavgpr[NPR],spcpr[NPR];
4699  static int firsttime=1;
4700  int pliter,pl;
4701 
4702  // return(0);
4703 
4704 
4705 
4706 
4708  //
4709  // Check/report some conditions to ensure can allow polesmooth()
4710  //
4712 
4713  // no need to process staggered primitive that are currently only field components
4714  if(ispstag) return(0);
4715 
4716  if(firsttime && special3dspc==0 && N3>1){ // not valid use of polesmooth()
4717  dualfprintf(fail_file,"SUPERWARNING: Must enable special3dspc==1 for it to be useful when N3>1\n");
4718  }
4719 
4720  if(firsttime && dofull2pi==0 && N3>1){
4721  dualfprintf(fail_file,"SUPERWARNING: Must enable dofull2pi==1 and have full 2\\pi doain when using polesmooth() for it to be useful when N3>1\n");
4722  }
4723 
4724 
4725 #if(DEBUGPOLESMOOTH)
4726  dualfprintf(fail_file,"Got here t=%g\n",t);
4727 #endif
4728 
4729 
4730 
4731 
4732 
4733 
4734 
4735 
4736 
4738  //
4739  // Allocate full span required for each myid
4740  //
4742  // need k slowest index so can easily do MPI sends/receives
4743 #define MAPFULLPR(i,fakej,k,pl) ( (fakej)*(N3BND+N3+N3BND)*(N1BND+N1+N1BND)*NPR + ((k)+N3BND)*(N1BND+N1+N1BND)*NPR + ((i)+N1BND)*NPR + (pl) )
4744  // so access as: fullpr[MAPFULLPR(i,fakej,k,pl)] . Order of i,fakej,k,pl is just to be similar to normal order of i,j,k,pl
4745 #define N1TOTFULLPR (N1BND+N1+N1BND)
4746 #define N2TOTFULLPR (ncpux3>1&&USEMPI&&N3>1 ? 1 : 2) // need 2 j positions if single cpu
4747 #define N3TOTFULLPR (N3BND+totalsize[3]+N3BND)
4748  if(firsttime){
4749  // not necessary if USEMPI==0 or ncpux3==1, but still do since not expensive
4750  fullpr=(FTYPE *)malloc(sizeof(FTYPE)*(N2TOTFULLPR*N3TOTFULLPR*N1TOTFULLPR*NPR));
4751  if(fullpr==NULL){
4752  dualfprintf(fail_file,"Cannot allocate fullpr\n");
4753  myexit(195367346);
4754  }
4755  // fullpr never freed, but repeatedly used
4756  }
4757 
4758 
4759 
4760 
4762  //
4763  // Setup which j=rj to use as reference and Setup interior j loop
4764  //
4766  if (whichx2==X2DN) {
4767  j0 = -POLESMOOTH; //starting from this j including this j
4768  rj = POLESMOOTH; //until this j but not including this j
4769  dj = 1;
4770  // stopj=rj+1; // can also average-out j=rj
4771  stopj=rj;
4772  fakej=0;
4773  }
4774  else{
4775  j0 = N2-1+POLESMOOTH; //starting from this j including this j
4776  rj = N2-1-POLESMOOTH; //until this j but not including this j
4777  dj = -1;
4778  // stopj=rj-1; // can also average-out j=rj
4779  stopj=rj;
4780  if(N2TOTFULLPR==1) fakej=0;
4781  else fakej=1; // 2 j memory slots since on 1 cpu
4782  }
4783 
4784 
4785 
4786 
4788  //
4789  // if MPI, then setup memory and pre-post recv's
4790  //
4792  static MPI_Request *srequest;
4793  static MPI_Request *rrequest;
4794 
4795  if(USEMPI && ncpux3>1 && N3>1){
4796 
4797  if(firsttime){
4798  srequest=(MPI_Request *)malloc(sizeof(MPI_Request)*ncpux3);
4799  if(srequest==NULL){
4800  dualfprintf(fail_file,"Cannot allocate srequest\n");
4801  myexit(19523766);
4802  }
4803  rrequest=(MPI_Request *)malloc(sizeof(MPI_Request)*ncpux3);
4804  if(rrequest==NULL){
4805  dualfprintf(fail_file,"Cannot allocate rrequest\n");
4806  myexit(346962762);
4807  }
4808  }
4809 
4810 
4811 
4812  if(mycpupos[2]==0 || mycpupos[2]==ncpux2-1){ // only done for j=rj near physical poles
4813  int posk;
4814 
4815  // transfer myid mycpupos[3] pr data to all other relevant cores (all other cores at this mycpupos[1],mycpupos[2] for all mycpupos[3])
4816  // non-blocking (all sends's occur at once)
4817 
4818  // receive all other portions of pr and fill local myid's fullpr data
4819  // non-blocking (all recv's occur at once)
4820 
4821  int count=N3*N1M*NPR;
4822 
4823  for(posk=0;posk<ncpux3;posk++){ // posk is absolute mycpupos[3] value
4824  if(posk!=mycpupos[3]){
4825 
4826  int originmyid=posk*ncpux2*ncpux1 + mycpupos[2]*ncpux1 + mycpupos[1];
4827  int recvtag=TAGSTARTBOUNDMPIPOLESMOOTH + originmyid + numprocs*mycpupos[3]; // recvtag used by other myid. Receiving for my mycpupos[3]
4828 
4829  // do recv
4830  MPI_Irecv(&fullpr[MAPFULLPR(-N1BND,fakej,posk*N3,0)] , count , MPI_FTYPE , MPIid[originmyid] , recvtag , MPI_COMM_GRMHD , &rrequest[posk]);
4831 #if(DEBUGPOLESMOOTH)
4832  dualfprintf(fail_file,"MPI_Irecv: posk=%d : %d %d\n",posk,originmyid,recvtag);
4833 #endif
4834 
4835  }
4836  }
4837  }
4838  }
4839 
4840 
4841 
4842 
4843 
4844 
4845 
4846 
4847 
4848 
4849 
4850 
4852  //
4853  // fullpr contains SPC versions of U1-U3 and PRIMECOORDS versions of other things (including scalars)
4854  // fill-in myid portion of fullpr
4855  // only a single j=rj
4856  //
4858  // do LOOPF1 in case near physical inner-radial or outer-radial boundaries that have already been set otherwise
4859  // Do transformation to SPC and Cart here to keep calculation distributed *and* because otherwise would have to fix "k" since don't have dxdxp or bl_coord stored for other k off this MPI process.
4860  LOOPF1 LOOPN3{
4861 
4863  // TRANSFORM PRIMECOORDS to quasi-SPC
4864  // Do here during fullpr store instead of later, because otherwise have to assume dxdxp[] does not vary in k and have to choose fixed k for dxdxp
4865  dxdxprim_ijk(i, rj, k, CENT, dxdxp);
4866  // No need to get_geometry() (that would be used to get lab-frame quantities) or to get ucon using pr2ucon(). Just process original frame quantities.
4867  // No need to get dxdxp using dxdxpprim_ijk() near pole except to obtain dominant scale. Major issue is twisting of coordinates, and dxdxp[1,2] and dxdxp[2,1] aren't large enough near pole to introduce major twisting.
4868 
4869  // set pr to pull from prim[]
4870  pr = MAC(prim,i,rj,k);
4871 
4873  // TRANSFORM PRIMECOORDS -> quasi-SPC
4874  PBOUNDLOOP(pliter,pl) if(pl<U1 || pl>B3) spcpr[pl] = pr[pl];
4875 
4876  // get quasi-SPC coordinates (avoid use of full dxdxp -- just need to get dimensional scalings correct for most part)
4877  spcpr[U1]=dxdxp[RR][RR]*pr[U1] + dxdxp[RR][TH]*pr[U2];
4878  spcpr[U2]=dxdxp[TH][RR]*pr[U1] + dxdxp[TH][TH]*pr[U2];
4879  spcpr[U3]=dxdxp[PH][PH]*pr[U3];
4880 
4881 
4882  // store in fullpr the SPC versions of U1-U3 and other scalars. Might as well avoid B1-B3 overwritten later.
4883  PBOUNDLOOP(pliter,pl) if(pl<B1 || pl>B3) fullpr[MAPFULLPR(i,fakej,mycpupos[3]*N3+k,pl)] = spcpr[pl];
4884 
4886  // TRANSFORM quasi-SPC to quasi-CART
4887  // quasi-Cart x, y, z (see tiltedAphi.m)
4888  // Do here so bl_coord() doesn't have to use fixed k or recompute for non-stored k values if looping over totalsize[3] per core later
4889  // get SPC V=t,r,\theta,\phi
4890  bl_coord_ijk_2(i,rj,k,CENT,X,V);
4891  r = V[1];
4892  th = V[2];
4893  ph = V[3];
4894  // No need to get_geometry() (that would be used to get lab-frame quantities) or to get ucon using pr2ucon(). Just process original frame quantities.
4895  // No need to get dxdxp using dxdxpprim_ijk() near pole except to obtain dominant scale. Major issue is twisting of coordinates, and dxdxp[1,2] and dxdxp[2,1] aren't large enough near pole to introduce major twisting.
4896 
4897  // put CART[U1-U3] into prfull[B1-B3] for storage
4898  fullpr[MAPFULLPR(i,fakej,mycpupos[3]*N3+k,B1)] = -(r*sin(ph)*sin(th)*spcpr[U3]) + cos(ph)*sin(th)*spcpr[U1] + r*cos(ph)*cos(th)*spcpr[U2];
4899  fullpr[MAPFULLPR(i,fakej,mycpupos[3]*N3+k,B2)] = r*cos(ph)*sin(th)*spcpr[U3] + sin(ph)*sin(th)*spcpr[U1] + r*cos(th)*sin(ph)*spcpr[U2];
4900  fullpr[MAPFULLPR(i,fakej,mycpupos[3]*N3+k,B3)] = cos(th)*spcpr[U1] - r*sin(th)*spcpr[U2];
4901 
4902 #if(DEBUGPOLESMOOTH)
4903  PBOUNDLOOP(pliter,pl) dualfprintf(fail_file,"Got hereINITIAL: t=%g i=%d j=%d k=%d : pl=%d pr=%g spc=%g cart=%g i/dxdxp11=%g %g\n",t,i,j,k,pl,pr[pl],spcpr[pl],fullpr[MAPFULLPR(i,fakej,mycpupos[3]*N3+k,pl)],idxdxp[RR][RR],dxdxp[RR][RR]);
4904 #endif
4905  }
4906 
4907 
4908 
4909 
4910 
4911 
4912 
4913 
4915  //
4916  // if MPI, then send my portion to other posk's and insert portions from all other posk's
4917  //
4919  if(USEMPI && ncpux3>1 && N3>1){
4920  MPI_Status mpistatus;
4921 
4922 
4923 
4924  if(mycpupos[2]==0 || mycpupos[2]==ncpux2-1){ // only done for j=rj near physical poles
4925  int posk;
4926 
4927  // transfer myid mycpupos[3] pr data to all other relevant cores (all other cores at this mycpupos[1],mycpupos[2] for all mycpupos[3])
4928  // non-blocking (all sends's occur at once)
4929 
4930  // receive all other portions of pr and fill local myid's fullpr data
4931  // non-blocking (all recv's occur at once)
4932 
4933  int count=N3*N1M*NPR;
4934 
4935  // do sends, with option to handshake first to ensure recv is ready and posted to avoid unexpected messages filling up buffer.
4936  for(posk=0;posk<ncpux3;posk++){ // posk is absolute mycpupos[3] value
4937  if(posk!=mycpupos[3]){
4938 
4939  int originmyid=posk*ncpux2*ncpux1 + mycpupos[2]*ncpux1 + mycpupos[1];
4940  int recvtag=TAGSTARTBOUNDMPIPOLESMOOTH + originmyid + numprocs*mycpupos[3]; // recvtag used by other myid. Receiving for my mycpupos[3]
4941 
4942  int destmyid=posk*ncpux2*ncpux1 + mycpupos[2]*ncpux1 + mycpupos[1];
4943  int sendtag=TAGSTARTBOUNDMPIPOLESMOOTH + myid + numprocs*posk; // sending to mycpupos[3]=posk (large sendtag space to ensure no duplicates)
4944 
4945  if(MPIFLOWCONTROL==1){
4946  // handshake before each send but after all recv's posted
4947  int nothingsend=0;
4948  int nothingrecv=0;
4949  int maxtag = numprocs*ncpux3; // might be somewhat limiting on numprocs for large number of ncpux3
4950 
4951  MPI_Sendrecv(
4952  &nothingsend,0,MPI_INT,
4953  MPIid[destmyid],
4954  maxtag + sendtag,
4955 
4956  &nothingrecv,0,MPI_INT,
4957  MPIid[originmyid],
4958  maxtag + recvtag,
4959 
4960  MPI_COMM_GRMHD,MPI_STATUS_IGNORE);
4961 
4962  } // end if doing FLOWCONTROL
4963 
4964 
4965  // now send
4966  MPI_Isend(&fullpr[MAPFULLPR(-N1BND,fakej,mycpupos[3]*N3,0)] , count , MPI_FTYPE , MPIid[destmyid] , sendtag , MPI_COMM_GRMHD , &srequest[posk]);
4967 #if(DEBUGPOLESMOOTH)
4968  dualfprintf(fail_file,"MPI_Isend: posk=%d : %d %d\n",posk,destmyid,sendtag);
4969 #endif
4970  }
4971  }
4972 
4973 
4974 
4975 
4976 
4977  // wait for all data to be recv'ed before moving to average that uses fullpr data and requires all data to be present in fullpr
4978  // wait for recv's
4979  for(posk=0;posk<ncpux3;posk++){
4980  if(posk!=mycpupos[3]){
4981  MPI_Wait(&rrequest[posk],&mpistatus); // assume successful so don't check mpistatus
4982 #if(DEBUGPOLESMOOTH)
4983  dualfprintf(fail_file,"MPI_Wait: posk=%d\n",posk);
4984 #endif
4985  }
4986  }
4987 
4988 
4989  // finally wait for sends
4990  for(posk=0;posk<ncpux3;posk++){
4991  if(posk!=mycpupos[3]){
4992  MPI_Wait(&srequest[posk],&mpistatus); // assume successful so don't check mpistatus
4993 #if(DEBUGPOLESMOOTH)
4994  dualfprintf(fail_file,"MPI_Wait: posk=%d\n",posk);
4995 #endif
4996  }
4997  }
4998 
4999 
5000  }// end if physical polar myid
5001  }// end if USEMPI and ncpux3>1 so have to transfer to other procs
5002 
5003 
5004 #if(DEBUGPOLESMOOTH)
5005  LOOPF1 for(k=0;k<totalsize[3];k++){
5006  PBOUNDLOOP(pliter,pl) if(pl==RHO || pl==B3) dualfprintf(fail_file,"FULLPRCHECK: i=%d k=%d pl=%d fullpr=%g\n",i,k,pl,fullpr[MAPFULLPR(i,fakej,k,pl)]);
5007  }
5008 #endif
5009 
5010 
5011 
5012 
5014  //
5015  // Loop over i (per i, no averaging or extrapolation in i)
5016  //
5018  LOOPF1{ // full i to account for already-assigned real boundary cells
5019 
5020 
5021 #if(DEBUGPOLESMOOTH)
5022  dualfprintf(fail_file,"Got here t=%g i=%d\n",t,i);
5023 #endif
5024 
5026  //
5027  // zero-out cumulating primitives
5028  //
5030  PBOUNDLOOP(pliter,pl){
5031  cartavgpr[pl]=0.0;
5032  spcavgpr[pl]=0.0;
5033  }
5034 
5035 
5036 
5038  //
5039  // loop over totalsize[3] for k and obtain average
5040  //
5042  for(k=0;k<totalsize[3];k++){ // only over active domain for averaging. Over full-k for fullpr
5043 
5044 #if(DEBUGPOLESMOOTH)
5045  dualfprintf(fail_file,"Got here t=%g i=%d k=%d\n",t,i,k);
5046 #endif
5047 
5048  // Averaged quasi-SPC versions (note spcavgpr[B1-B3] not used and fullpr[B1-B3] stored quasi-CART U1-U3, so might as well avoid B1-B3)
5049  PBOUNDLOOP(pliter,pl) if(pl<B1 || pl>B3) spcavgpr[pl] += fullpr[MAPFULLPR(i,fakej,k,pl)];
5050 
5051  // Averaged non-velocities for quasi-CART (avoid cartavgpr[pl=U1-U3] that cumulate next. Also go ahead and avoid B1-B3 since fullpr[B1-B3] filled will cart U1-U3
5052  PBOUNDLOOP(pliter,pl) if(pl<U1 || pl>B3) cartavgpr[pl] += fullpr[MAPFULLPR(i,fakej,k,pl)];
5053 
5054  // cumulate quasi-CART U1-U3 (stored in fullpr[B1-B3])
5055  cartavgpr[U1] += fullpr[MAPFULLPR(i,fakej,k,B1)];
5056  cartavgpr[U2] += fullpr[MAPFULLPR(i,fakej,k,B2)];
5057  cartavgpr[U3] += fullpr[MAPFULLPR(i,fakej,k,B3)];
5058 
5059  }// end cumulation over k at j=rj
5060 
5061 
5063  //
5064  // Obtain average by dividing by numer of terms (assume use *all* \phi cells per i,j)
5065  // for USEMPI==1 && ncpux3>1, assume all CPUs got same average so no need to share average. Will be same to machine precision, which is sufficient since no inconsistency across MPI boundaries in the end -- these machine errors just make average different for each k, which is perfectly fine.
5066  //
5068  PBOUNDLOOP(pliter,pl){
5069  if(pl<B1 || pl>B3){ // skip unused B1-B3
5070  cartavgpr[pl] /= (FTYPE)totalsize[3];
5071  spcavgpr[pl] /= (FTYPE)totalsize[3];
5072  }
5073  }
5074 
5075  // account for dofull2pi==0 or 2D axisymmetry
5076  // in such a case, can't have net flow across axis
5077  // (so similar to poledeath() in 2D limit or not-fully-3d dofull2pi=0 limit)
5078  if(dofull2pi==0 || totalsize[3]==1){
5079  cartavgpr[U1] = 0.0; // net vx=0 across axis
5080  cartavgpr[U2] = 0.0; // net vy=0 across axis
5081  }
5082 
5083 #if(DEBUGPOLESMOOTH)
5084  PBOUNDLOOP(pliter,pl) dualfprintf(fail_file,"Got hereCARTSPC: t=%g i=%d pl=%d %g %g\n",t,i,pl,cartavgpr[pl],spcavgpr[pl]);
5085 #endif
5086 
5087 
5088 
5090  //
5091  // Populate the interior (to stopj) j cells with averaged primitives
5092  //
5094  LOOPF3{// over full domain for assignment of the average since boundary call for periodic x3 may already (as is normal) be done.
5095  for (j=j0; j != stopj; j+=dj) { // over interior j to stopj that is (typically) rj
5096 
5097 
5099  // TRANSFORM CART to quasi-SPC
5100  // Current interior-pole location (j, which is inside previous rj location)
5101  bl_coord_ijk_2(i,j,k,CENT,X,V);
5102  r = V[1];
5103  th = V[2];
5104  ph = V[3];
5105 
5106  // Set non-velocity SPC (cartavgpr[B1-B3] has nothing, so avoid)
5107  PBOUNDLOOP(pliter,pl) if(pl<U1 || pl>B3) spcpr[pl] = spcavgpr[pl]; // could also use cartavgpr too. Same results for scalars.
5108 
5109  // Set quasi-SPC for i,j,k from phi-averaged quasi-Cart from i,j=rj,all k
5110  // just inverse transformation of above cartavgpr(pr)
5111  spcpr[U1] = cos(ph)*sin(th)*cartavgpr[U1] + sin(ph)*sin(th)*cartavgpr[U2] + cos(th)*cartavgpr[U3];
5112  spcpr[U2] = pow(r,-1)*(cos(ph)*cos(th)*cartavgpr[U1] + cos(th)*sin(ph)*cartavgpr[U2] - sin(th)*cartavgpr[U3]);
5113  spcpr[U3] = (1./sin(th))*pow(r,-1)*(-sin(ph)*cartavgpr[U1] + cos(ph)*cartavgpr[U2]);
5114 
5115  // The average of v_x and v_y around axis removes rotation around axis, so add back-in the phi-averaged rotational velocity
5116  spcpr[U3] += spcavgpr[U3];
5117 
5118  // TODOMARK GODMARK: Why not also add back-in average compression/expansion around axis in \theta direction by doing:
5119  // Might be better for shocks near axis
5120  // spcpr[U2] += spcavgpr[U2];
5121 
5122  // NOTEMARK: spcpr[U1]=r is just like cartavgpr[U3]=z, so no recovery of net motion required.
5123 
5124 
5126  // TRANSFORM quasi-SPC to PRIMECOORDS
5127 
5128  // set pr to assign
5129  pr = MAC(prim,i,j,k);
5130 
5131  // Set other non-velocity, non-field things (DO NOT OVERWRITE FIELD!)
5132  PBOUNDLOOP(pliter,pl) if(pl<U1 || pl>B3) pr[pl] = spcpr[pl];
5133 
5134  // Since using simplified dxdxp above, must use inverse of that for consistency. Cannot use full idxdxp unless used full dxdxp.
5135  // idxdxprim_ijk(i, j, k, CENT, idxdxp);
5136  dxdxprim_ijk(i, j, k, CENT, dxdxp);
5137  idxdxp[RR][RR]=dxdxp[TH][TH]/(dxdxp[TH][TH]*dxdxp[RR][RR]-dxdxp[TH][RR]*dxdxp[RR][TH]);
5138  idxdxp[RR][TH]=dxdxp[RR][TH]/(dxdxp[TH][RR]*dxdxp[RR][TH]-dxdxp[TH][TH]*dxdxp[RR][RR]);
5139  idxdxp[TH][RR]=dxdxp[TH][RR]/(dxdxp[TH][RR]*dxdxp[RR][TH]-dxdxp[TH][TH]*dxdxp[RR][RR]);
5140  idxdxp[TH][TH]=dxdxp[RR][RR]/(dxdxp[TH][TH]*dxdxp[RR][RR]-dxdxp[TH][RR]*dxdxp[RR][TH]);
5141  idxdxp[PH][PH]=1.0/dxdxp[PH][PH];
5142 
5143  // note that the below idxdp[] is transposed compared to how would act on u_\mu
5144  pr[U1]=idxdxp[RR][RR]*spcpr[U1] + idxdxp[RR][TH]*spcpr[U2];
5145  pr[U2]=idxdxp[TH][RR]*spcpr[U1] + idxdxp[TH][TH]*spcpr[U2];
5146  pr[U3]=idxdxp[PH][PH]*spcpr[U3];
5147 
5148  // DONE! Have full pr=prim[] set now
5149 
5150 
5151 
5152  // NOW set boundary conditions consistent with full3d treatment
5153  // NO! Boundary cell values only change in \theta, not \phi as physical cells would, for given j. So no changes to sign should occur.
5154  // if(j<0 && whichx2==X2DN || j>totalsize[2]-1 && whichx2==X2UP){
5155  // pr[U1]*=SIGNFLIP12;
5156  // pr[U2]*=SIGNFLIPU2;
5157  // pr[U3]*=SIGNFLIPU3;
5158  // // Note that B2,B3 are not ever modified, so no need for sign change
5159  // }
5160 
5161 
5162 #if(DEBUGPOLESMOOTH)
5163  PBOUNDLOOP(pliter,pl) dualfprintf(fail_file,"Got hereFINAL: t=%g i=%d j=%d k=%d : pl=%d pr=%g\n",t,i,j,k,pl,pr[pl]);
5164 #endif
5165 
5166 
5167  }// end over each j
5168  }// end over each k
5169 
5170 
5171  }// end over each i
5172 
5173 
5174  firsttime=0;
5175 
5176 
5177  return(0);
5178 
5179 }
5180 
5181 
5182 
5183 
5186 void user1_adjust_fluxcttoth_emfs(SFTYPE time, FTYPE (*prim)[NSTORE2][NSTORE3][NPR], FTYPE (*emf)[NSTORE1+SHIFTSTORE1][NSTORE2+SHIFTSTORE2][NSTORE3+SHIFTSTORE3] )
5187 {
5188  int i, j, k;
5189  int dir;
5190  int dimen;
5191  int dirsign;
5192 
5193 
5194  if(ADJUSTFLUXCT==0){
5195  // then no need to set EMF's to zero to maintain divb=0, although can choose to still set EMF's to constant along appropriate directions to maintain stationarity all the way up to including the boundary. Not sure why that would be important to maintain.
5196  return;
5197  }
5198 
5199 
5200  if( DOGRIDSECTIONING==0) { //don't do anything if sectioning not done
5201  return;
5202  }
5203 
5204 
5205  if(BCtype[X1DN]==FIXEDUSEPANALYTIC || BCtype[X1UP]==FIXEDUSEPANALYTIC){
5206  // only do if fixing BCs, not (e.g.) outflowing
5207 
5208  //x1-dimension
5209  dimen = 1;
5210  DIRSIGNLOOP(dirsign){
5211  dir = DIRFROMDIMEN( dimen, dirsign );
5212 
5213  if(dir==X1DN && BCtype[X1DN]==FIXEDUSEPANALYTIC || dir==X1UP && BCtype[X1UP]==FIXEDUSEPANALYTIC){ // otherwise don't do
5214 
5215  i = dofluxreg[ACTIVEREGION][dir];
5216 
5217  //if boundary is not on this processor, do not modify emf's
5218  if( i < -MAXBND ) continue;
5219  if( i > N1-1+MAXBND) continue;
5220 
5221  //the boundary is on the processor, so reset emf's to zero at the boundary
5222  if(dir==X1UP && BCtype[dir]==FIXEDUSEPANALYTIC){ // i.e. don't reset EMF2=0 for lower boundary since that corresponds to rotation
5224  MACP1A0(emf,2,i,j,k) = 0.0;
5225  }
5226  }
5228  MACP1A0(emf,3,i,j,k) = 0.0;
5229  }
5230  }
5231  }
5232  }
5233 
5234 
5235 #if(0) // below not currently used
5236 
5237  //x2-dimension
5238  dimen = 2;
5239 
5240 
5241  //SASMARKXXX -- DO demand stationarity at the polar axis
5242  DIRSIGNLOOP(dirsign){
5243  dir = DIRFROMDIMEN( dimen, dirsign );
5244 
5245  j = dofluxreg[ACTIVEREGION][dir];
5246 
5247  //if boundary is not on this processor, do not modify emf's
5248  if( j < -MAXBND ) continue;
5249  if( j > N2-1+MAXBND ) continue;
5250 
5251  //the boundary is on the processor, so reset emf's to zero at the boundary
5253  MACP1A0(emf,1,i,j,k) = 0.0;
5254  MACP1A0(emf,3,i,j,k) = 0.0;
5255  }
5256  }
5257 
5258  //x3-dimension
5259  dimen = 3;
5260 
5261 
5262  DIRSIGNLOOP(dirsign){
5263  dir = DIRFROMDIMEN( dimen, dirsign );
5264 
5265  k = dofluxreg[ACTIVEREGION][dir];
5266 
5267  //if boundary is not on this processor, do not modify emf's
5268  if( k < -MAXBND ) continue;
5269  if( k > N3-1+MAXBND ) continue;
5270 
5271  //the boundary is on the processor, so reset emf's to zero at the boundary
5273  MACP1A0(emf,1,i,j,k) = 0.0;
5274  MACP1A0(emf,2,i,j,k) = 0.0;
5275  }
5276  }
5277 #endif
5278 
5279 
5280 }
5281 
5285 void user1_adjust_fluxctstag_emfs(SFTYPE time, FTYPE (*prim)[NSTORE2][NSTORE3][NPR], int *Nvec, FTYPE (*fluxvec[NDIM])[NSTORE2][NSTORE3][NPR+NSPECIAL])
5286 {
5287  int i, j, k;
5288  int dir;
5289  int dimen;
5290  int dirsign;
5291 
5292  if(ADJUSTFLUXCT==0){
5293  // then no need to set EMF's to zero to maintain divb=0, although can choose to still set EMF's to constant along appropriate directions to maintain stationarity all the way up to including the boundary. Not sure why that would be important to maintain.
5294  return;
5295  }
5296 
5297  if( DOGRIDSECTIONING==0) { //don't do anything if sectioning not done
5298  return;
5299  }
5300 
5301 
5303  // only do if fixing BCs, not (e.g.) outflowing
5304 
5305  //x1-dimension
5306  dimen = 1;
5307  DIRSIGNLOOP(dirsign){
5308  dir = DIRFROMDIMEN( dimen, dirsign );
5309 
5310 
5311  if(dir==X1DN && BCtype[X1DN]==FIXEDUSEPANALYTIC || dir==X1UP && BCtype[X1UP]==FIXEDUSEPANALYTIC){ // otherwise don't do
5312 
5313  //if boundary is not on this processor, do not modify emf's
5314  i = dofluxreg[ACTIVEREGION][dir];
5315  if( i < -MAXBND ) continue;
5316  if( i > N1-1+MAXBND) continue;
5317 
5318 
5319  //the boundary is on the processor, so reset emf's to zero at the boundary
5321  // EMF[3]:
5322 #if(N1>1)
5323  MACP1A1(fluxvec,1,i,j,k,B2) = 0.0;
5324  MACP1A1(fluxvec,1,i,j,k,B1) = 0.0; // always zero
5325 #endif
5326 #if(N2>1)
5327  MACP1A1(fluxvec,2,i,j,k,B1) = 0.0;
5328  MACP1A1(fluxvec,2,i,j,k,B2) = 0.0; // always zero
5329 #endif
5330  }
5331  }
5332 
5333 
5334  if(dir==X1UP && BCtype[dir]==FIXEDUSEPANALYTIC){ // i.e. don't reset EMF2=0 for lower boundary since that corresponds to rotation
5335 
5336  //if boundary is not on this processor, do not modify emf's
5337  i = dofluxreg[ACTIVEREGION][dir];
5338  if( i < -MAXBND ) continue;
5339  if( i > N1-1+MAXBND) continue;
5340 
5341 
5343  // EMF[2]:
5344 #if(N1>1)
5345  MACP1A1(fluxvec,1,i,j,k,B3) = 0.0;
5346  MACP1A1(fluxvec,1,i,j,k,B1) = 0.0; // always zero
5347 #endif
5348 #if(N3>1)
5349  MACP1A1(fluxvec,3,i,j,k,B1) = 0.0;
5350  MACP1A1(fluxvec,3,i,j,k,B3) = 0.0; // always zero
5351 #endif
5352  }
5353  }
5354 
5355 
5356  }
5357  }
5358 
5359 
5360 
5361 }
5362 
5363 
5364 
5365 
5366 
5367 
5370 void user1_adjust_final_fluxes(SFTYPE time, FTYPE (*prim)[NSTORE2][NSTORE3][NPR], int *Nvec, FTYPE (*fluxvec[NDIM])[NSTORE2][NSTORE3][NPR+NSPECIAL])
5371 {
5372  int i, j, k;
5373  int dir;
5374  int dimen;
5375  int dirsign;
5376 
5377  if( DOGRIDSECTIONING==0) { //don't do anything if sectioning not done
5378  return;
5379  }
5380 
5381 
5382  //x1-dimension
5383  dimen = 1;
5384  DIRSIGNLOOP(dirsign){
5385  dir = DIRFROMDIMEN( dimen, dirsign );
5386 
5387 
5388  //if boundary is not on this processor, do not modify emf's
5389  i = dofluxreg[ACTIVEREGION][dir];
5390  if( i < -MAXBND ) continue;
5391  if( i > N1-1+MAXBND) continue;
5392 
5393 
5394  //the boundary is on the processor, so reset emf's to zero at the boundary
5396  // EMF[3]:
5397 #if(N1>1)
5398  MACP1A1(fluxvec,1,i,j,k,B2) = 0.0;
5399  MACP1A1(fluxvec,1,i,j,k,B1) = 0.0; // always zero
5400 #endif
5401 #if(N2>1)
5402  MACP1A1(fluxvec,2,i,j,k,B1) = 0.0;
5403  MACP1A1(fluxvec,2,i,j,k,B2) = 0.0; // always zero
5404 #endif
5405  // EMF[2]:
5406 #if(N1>1)
5407  MACP1A1(fluxvec,1,i,j,k,B3) = 0.0;
5408  MACP1A1(fluxvec,1,i,j,k,B1) = 0.0; // always zero
5409 #endif
5410 #if(N3>1)
5411  MACP1A1(fluxvec,3,i,j,k,B1) = 0.0;
5412  MACP1A1(fluxvec,3,i,j,k,B3) = 0.0; // always zero
5413 #endif
5414  }
5415  }
5416 
5417 
5418 
5419 
5420 }
5421 
5422 
5423 
5424 
5425 
5431 void check_spc_singularities_user(void)
5432 {
5433  int i,j,k;
5434  int poleloop;
5435  int needzero;
5436  int singfound;
5437  int nonsingfound;
5438  int numlocs=NPG,indloc,loc;
5439  int pl,pliter;
5441  int doprintout;
5442 
5443  if(nstep==0) doprintout=1;
5444  else doprintout=0; // avoid print out when evolving metric and recalling this -- assume faily similar situation as at nstep==0
5445 
5447  //
5448  // X1 inner r=0 singularity
5449  //
5451 
5452 
5453  /* inner radial BC (preserves u^t rho and u) */
5454  if(mycpupos[1] == 0 && BCtype[X1DN]==R0SING){
5455  i=0;
5456 
5457  // LOOPX1dir{
5458  // No, should really be full loops
5459  LOOPF_23{
5460 
5461  for(indloc=0;indloc<numlocs;indloc++){
5462 
5463  if(indloc==0) loc=FACE1;
5464  else if(indloc==1) loc=CORN3;
5465  else if(indloc==2) loc=CORN2;
5466  else if(indloc==3) loc=CORNT;
5467 
5468  // get local metric quantities for this loc,i,j,k
5469  GETLOCALMETRIC(loc,i,j,k);
5470 
5471  singfound=0;
5472  singfound+=(localgdet[0]!=0.0);
5473 #if(WHICHEOM!=WITHGDET)
5474  PLOOP(pliter,pl) singfound+=(LOCALEOMFUNCMAC(pl)!=0.0);
5475 #endif
5476 
5477  if(singfound){
5478  if(doprintout) dualfprintf(fail_file,"Detected singularity at i=%d j=%d k=%d loc=%d with gdet=%21.15g so resetting it to 0.0\n",i,j,k,loc,localgdet[0]);
5479  localgdet[0]=0.0;
5480  if(NEWMETRICSTORAGE){
5481  GLOBALMETMACP1A0(compgeom,loc,i,j,k).gdet=0.0;
5482  GLOBALMETMACP1A0(compgeom,loc,i,j,k).igdetnosing=0.0;
5483  GLOBALMETMACP1A0(gdetgeomnormal,loc,i,j,k).gdet=0.0;
5484  GLOBALMETMACP1A0(gdetgeomnormal,loc,i,j,k).igdetnosing=0.0;
5485  GLOBALMETMACP0A1(gdetgeom,i,j,k,loc).gdet=0.0;
5486  GLOBALMETMACP0A1(gdetgeom,i,j,k,loc).igdetnosing=0.0;
5487  }
5488 #if(WHICHEOM!=WITHGDET)
5489  PLOOP(pliter,pl) LOCALEOMFUNCMAC(pl)=0.0;
5490 #endif
5491 
5492 
5493  }
5494  }// end over indloc
5495  }// end loop 23
5496  }
5497 
5498 
5499 
5501  //
5502  // X2 POLARAXES
5503  //
5505 
5506 
5507 
5508  if(special3dspc && COORDSINGFIX>0){
5509  // Then presume really doing 3D problem with \phi direction
5510  // Then presume that (for simplicity of the method) offsetting pole by \epsilon from \theta=0,\pi so that \detg is generally non-zero
5511  // Then flux through pole, but transmissive BCs at pole ensure flux into "\epsilon hole" are same as out of the hole on other side
5512  // In any case that flux is very small. Key point is that for d_t(\detg B2) = -d_1(\detg E3) + d_3(\detg E1), the small angular term in \detg cancels leaving regular B2 on the pole that enters interpolations of B2.
5513 
5514  // Check that \detg is NOT zero now
5515  needzero=0;
5516  }
5517  else{
5518  // Then presume no \epsilon offset at pole
5519  needzero=1;
5520  }
5521 
5522 
5523 
5524  for(poleloop=0;poleloop<=1;poleloop++){
5525  if(poleloop==0 && mycpupos[2] == 0 && BCtype[X2DN]==POLARAXIS){
5526  j=0;
5527  }
5528  else if(poleloop==1 && mycpupos[2] == ncpux2-1 && BCtype[X2UP]==POLARAXIS){
5529  j=N2-1+SHIFT2;
5530  }
5531  else continue;
5532 
5533 
5534  // LOOPX2dir{
5535  // No, should really be full loops:
5536  LOOPF_13{
5537 
5538  for(indloc=0;indloc<numlocs;indloc++){
5539 
5540  if(indloc==0) loc=FACE2;
5541  else if(indloc==1) loc=CORN3;
5542  else if(indloc==2) loc=CORN1;
5543  else if(indloc==3) loc=CORNT;
5544 
5545  // get local metric quantities for this loc,i,j,k
5546  GETLOCALMETRIC(loc,i,j,k);
5547 
5548  singfound=0;
5549  singfound+=(fabs(localgdet[0])>0.0);
5550 #if(WHICHEOM!=WITHGDET)
5551  PLOOP(pliter,pl) singfound+=(fabs(LOCALEOMFUNCMAC(pl))>0.0);
5552 #endif
5553 
5554  nonsingfound=0;
5555  nonsingfound+=(fabs(localgdet[0])<100*NUMEPSILON);
5556 #if(WHICHEOM!=WITHGDET)
5557  PLOOP(pliter,pl) nonsingfound+=(fabs(LOCALEOMFUNCMAC(pl))<100*NUMEPSILON);
5558 #endif
5559 
5560  if(needzero && singfound){
5561  if(doprintout) dualfprintf(fail_file,"Detected singularity at i=%d j=%d k=%d loc=%d with gdet=%21.15g so resetting it to 0.0\n",i,j,k,loc,localgdet[0]);
5562  localgdet[0]=0.0;
5563  if(NEWMETRICSTORAGE){
5564  GLOBALMETMACP1A0(compgeom,loc,i,j,k).gdet=0.0;
5565  GLOBALMETMACP1A0(compgeom,loc,i,j,k).igdetnosing=0.0;
5566  GLOBALMETMACP1A0(gdetgeomnormal,loc,i,j,k).gdet=0.0;
5567  GLOBALMETMACP1A0(gdetgeomnormal,loc,i,j,k).igdetnosing=0.0;
5568  GLOBALMETMACP0A1(gdetgeom,i,j,k,loc).gdet=0.0;
5569  GLOBALMETMACP0A1(gdetgeom,i,j,k,loc).igdetnosing=0.0;
5570  }
5571 #if(WHICHEOM!=WITHGDET)
5572  PLOOP(pliter,pl) LOCALEOMFUNCMAC(pl)=0.0;
5573 #endif
5574 
5575  }
5576  else if(needzero==0 && nonsingfound){
5577  dualfprintf(fail_file,"Detected |\\detg|<=NUMEPSILON at i=%d j=%d k=%d loc=%d with gdet=%21.15g: Must be non-zero\n",i,j,k,loc,localgdet[0]);
5578  myexit(9813523);
5579  }
5580  }// end over indloc
5581  }// end loop 13
5582  }// end over poles
5583 
5584 }
5585 
5586 
5587 
5591 int debugspecial3dspc(int which, int whichdir, int ispstag, FTYPE (*prim)[NSTORE2][NSTORE3][NPR])
5592 {
5593  int i,j,k,pliter,pl;
5594 
5595 
5596  dualfprintf(fail_file,"nstep=%ld steppart=%d which=%d whichdir=%d ispstag=%d\n",nstep,steppart,which,whichdir,ispstag);
5597  for(pl=0;pl<=7;pl++){
5598  if(ispstag==1 && (pl<5 || pl>7)) continue;
5599  dualfprintf(fail_file,"pl=%d\n",pl);
5600  for(k=-N3BND;k<=N3BND;k++){ // was looking at just k=0, but seeing now if all other k are good.
5601  dualfprintf(fail_file,"k=%d\n",k);
5602  dualfprintf(fail_file,"%02s | i=\n ","j");
5603  for(i=-N1BND;i<=N1BND;i++){
5604  if(ispstag==1 && (pl==5 && i==-N1BND)) dualfprintf(fail_file,"%021s ","NA");
5605  else dualfprintf(fail_file,"%19s%02d "," ",i);
5606  }
5607  dualfprintf(fail_file,"\n");
5608  for(j=-N2BND;j<=N2BND;j++){
5609  dualfprintf(fail_file,"%02d ",j);
5610  for(i=-N1BND;i<=N1BND;i++){
5611  if(ispstag==1 && (pl==5 && i==-N1BND || pl==6 && j==-N2BND)) dualfprintf(fail_file, "%21s ","NA");
5612  else dualfprintf(fail_file, "%21.15g ",MACP0A1(prim,i,j,k,pl));
5613  }
5614  dualfprintf(fail_file,"\n");
5615  }
5616  dualfprintf(fail_file,"\n");
5617  }
5618  }// end over pl
5619 
5620 
5621  return(0);
5622 
5623 }
5624 
5625 
5626 
5628 void debugfixupaltdeath_bc_old(FTYPE (*prim)[NSTORE2][NSTORE3][NPR])
5629 {
5630  // hack to get rid of bad region at large distances (e.g. when restarted)
5631  int i,j,k;
5632  struct of_geom geomdontuse;
5633  struct of_geom *ptrgeom=&geomdontuse;
5634  FTYPE X[NDIM],V[NDIM];
5635  FTYPE *prfix,*ufix;
5636  int jjj;
5637 
5638 
5639  FULLLOOP{
5640  prfix=&MACP0A1(prim,i,j,k,0);//&GLOBALMACP0A1(pglobal,i,j,k,0);
5641  ufix=&GLOBALMACP0A1(unewglobal,i,j,k,0);
5642 
5643  // get geometry for center pre-interpolated values
5644  get_geometry(i, j, k, CENT, ptrgeom);
5645  bl_coord_ijk_2(i,j,k,CENT,X, V) ;
5646 
5647  // FTYPE bsq=0.0;
5648  // bsq_calc(prfix,ptrgeom,&bsq);
5649 
5650 
5651  FTYPE Rbadout=OUTERDEATHRADIUS; //MIN(500.0,0.5*Rout);
5652 
5653 
5654  if(V[1]>Rbadout ){
5655 
5656  //prfix[RHO] = 1E-10*pow(V[1]/500.0,-1.5);
5657 
5658  prfix[UU]= MIN(prfix[RHO],prfix[UU]); // no more than u/rho=1
5659  if(DOENOFLUX != NOENOFLUX){
5660  ufix[UU]=MAX(-prfix[UU],ufix[UU]);
5661  ufix[ENTROPY] = ufix[UU];
5662  ufix[ENTROPY] = MAX(0.0001,MIN(ufix[ENTROPY],1.0)); // like u/rho=1
5663  }
5664 
5665  if(URAD0>=0){
5666  //prfix[URAD0] = 2E-10*pow(V[1]/500.0,-1.5);
5667  prfix[URAD0] = MIN(MIN(prfix[RHO],prfix[URAD0]),prfix[UU]); // no more than Erf/rho=1 and Erf/u=1
5668  if(DOENOFLUX != NOENOFLUX){
5669  ufix[URAD0]=MAX(-prfix[URAD0],ufix[URAD0]);
5670  }
5671  }
5672 
5673  // limit_gamma(0,1.5,GAMMAMAXRAD,prfix,NULL,ptrgeom,-1);
5674  limit_gamma(0,OUTERDEATHGAMMAMAX,OUTERDEATHGAMMAMAXRAD,prfix,NULL,ptrgeom,-1);
5675 
5676 
5677 
5678  if(prfix[U1]<0.0){
5679  prfix[U1]=0.0;
5680  if(DOENOFLUX != NOENOFLUX){
5681  ufix[U1]=0.0;
5682  }
5683  }
5684  // ufix[U2]=ufix[U3]=0.0;
5685  // SLOOPA(jjj) ufix[U1+jjj-1]=prfix[U1+jjj-1] = 0.0;
5686 
5687  if(URAD0>=0){
5688  if(prfix[URAD1]<0.0){
5689  prfix[URAD1]=0.0;
5690  if(DOENOFLUX != NOENOFLUX){
5691  ufix[URAD1]=0.0;
5692  }
5693  }
5694  }
5695 
5696  }
5697  }
5698 }
5699 
5700 
5702 void debugfixupaltdeath_bc(FTYPE (*prim)[NSTORE2][NSTORE3][NPR])
5703 {
5704  // hack to get rid of bad region at large distances (e.g. when restarted)
5705  int i,j,k;
5706  struct of_geom geomdontuse;
5707  struct of_geom *ptrgeom=&geomdontuse;
5708  FTYPE X[NDIM],V[NDIM];
5709  FTYPE *prfix,*ufix;
5710  int jjj;
5711 
5712  if(VARTOINTERP==PRIMTOINTERP_GDETFULLVERSION_WALD){
5713 
5714  // 0 : don't do
5715  // 1: original method
5716  // 2: new method1
5717 #define AVOIDCS 0
5718 
5719  if(AVOIDCS==1){
5720  FULLLOOP{
5721  prfix=&MACP0A1(prim,i,j,k,0);//&GLOBALMACP0A1(pglobal,i,j,k,0);
5722  ufix=&GLOBALMACP0A1(unewglobal,i,j,k,0);
5723 
5724  // get geometry for center pre-interpolated values
5725  //get_geometry(i, j, k, CENT, ptrgeom);
5726  //bl_coord_ijk_2(i,j,k,CENT,X, V) ;
5727 
5728  // FTYPE bsq=0.0;
5729  // bsq_calc(prfix,ptrgeom,&bsq);
5730 
5731 
5732  // symmetric around equator for centered zones
5733  int ti=startpos[1]+i;
5734  int tj=startpos[2]+j;
5735  int tk=startpos[3]+k;
5736  if((tj>totalsize[2]/2-3)&&(tj<totalsize[2]/2+2)){
5737  prfix[U2]=0.0;
5738  }
5739 
5740  }
5741  }
5742  else if(AVOIDCS==2){
5743  FULLLOOP{
5744 
5745  // get geometry for center pre-interpolated values
5746  get_geometry(i, j, k, CENT, ptrgeom);
5747  bl_coord_ijk_2(i,j,k,CENT,X, V) ;
5748 
5749 
5750  FTYPE Rhorlocal=1.0+sqrt(1.0-a*a);
5751  // Rhorlocal=2.0; // use ergosphere instead.
5752  if(V[1]<Rhorlocal){
5753 
5754  prfix=&MACP0A1(prim,i,j,k,0);//&GLOBALMACP0A1(pglobal,fii,fjj,fkk,0);
5755  if(DOENOFLUX != NOENOFLUX){
5756  ufix=&GLOBALMACP0A1(unewglobal,i,j,k,0);
5757  }
5758  prfix[U2]=0.0;
5759  }
5760 
5761  }// over FULLLOOP
5762  }
5763  else if(AVOIDCS==3){
5764  FULLLOOP{
5765 
5766  // get geometry for center pre-interpolated values
5767  //get_geometry(i, j, k, CENT, ptrgeom);
5768  //bl_coord_ijk_2(i,j,k,CENT,X, V) ;
5769 
5770 
5771  FTYPE Rhorlocal=1.0+sqrt(1.0-a*a);
5772  // if(V[1]<Rhorlocal){
5773  FTYPE maxbsq=0.0,mybsq=0.0;
5774  int ii,jj,kk;
5775  for(ii=-N1BND;ii<=N1BND;ii++){
5776  for(jj=-N2BND;jj<=N2BND;jj++){
5777  for(kk=-N3BND;kk<=N3BND;kk++){
5778  int fii=MIN(MAX(i+ii,-N1BND),N1-1+N1BND);
5779  int fjj=MIN(MAX(j+jj,-N2BND),N2-1+N2BND);
5780  int fkk=MIN(MAX(k+kk,-N3BND),N3-1+N3BND);
5781 
5782  prfix=&MACP0A1(prim,fii,fjj,fkk,0);//&GLOBALMACP0A1(pglobal,fii,fjj,fkk,0);
5783  ufix=&GLOBALMACP0A1(unewglobal,fii,fjj,fkk,0);
5784 
5785 
5786  // get geometry for center pre-interpolated values
5787  get_geometry(fii,fjj,fkk, CENT, ptrgeom);
5788  //bl_coord_ijk_2(fii,fjj,fkk,CENT,X, V) ;
5789 
5790 
5791  FTYPE bsq=0.0;
5792  bsq_calc(prfix,ptrgeom,&bsq);
5793 
5794  maxbsq=MAX(maxbsq,bsq);
5795 
5796  if(ii==0 && jj==0 && kk==0) mybsq=bsq;
5797 
5798  }
5799  }
5800  }
5801 
5802  if(mybsq*10.0<maxbsq){
5803 
5804  dualfprintf(fail_file,"i=%d j=%d k=%d mybsq=%g maxbsq=%g\n",i,j,k,mybsq,maxbsq);
5805 
5806  // prfix=&MACP0A1(prim,i,j,k,0);//&GLOBALMACP0A1(pglobal,fii,fjj,fkk,0);
5807  // ufix=&GLOBALMACP0A1(unewglobal,i,j,k,0);
5808  // prfix[U2]=0.0;
5809 
5810  for(ii=-N1BND/2;ii<=N1BND/2;ii++){
5811  for(jj=-N2BND/2;jj<=N2BND/2;jj++){
5812  for(kk=-N3BND/2;kk<=N3BND/2;kk++){
5813  int fii=MIN(MAX(i+ii,-N1BND),N1-1+N1BND);
5814  int fjj=MIN(MAX(j+jj,-N2BND),N2-1+N2BND);
5815  int fkk=MIN(MAX(k+kk,-N3BND),N3-1+N3BND);
5816 
5817  prfix=&MACP0A1(prim,fii,fjj,fkk,0);//&GLOBALMACP0A1(pglobal,fii,fjj,fkk,0);
5818  ufix=&GLOBALMACP0A1(unewglobal,fii,fjj,fkk,0);
5819 
5820  prfix[U2]=0.0;
5821  }
5822  }
5823  }
5824 
5825  }
5826 
5827  }// over FULLLOOP
5828  }
5829 
5830  // do nothing, since when including low density and high b^2/rho, has issues.
5831  // return; // flipping back to return because shifting floor in rho,u didn't help MHD heating.
5832  }
5833  else{
5834  // allow outerdeath
5835 
5836  // don't allow in general
5837  //return;
5838  }
5839 
5840 
5841  if(0){
5842 
5843 
5844 #define DAMPECOV 0
5845 #define DAMPGAMMA 0
5846 #define OUTERDEATHDAMP 0
5847 
5848  if(DAMPGAMMA || OUTERDEATHDAMP || DAMPECOV){
5849 
5850  FULLLOOP{
5851  prfix=&MACP0A1(prim,i,j,k,0);//&GLOBALMACP0A1(pglobal,i,j,k,0);
5852  ufix=&GLOBALMACP0A1(unewglobal,i,j,k,0);
5853 
5854  // get geometry for center pre-interpolated values
5855  get_geometry(i, j, k, CENT, ptrgeom);
5856  bl_coord_ijk_2(i,j,k,CENT,X, V) ;
5857 
5858  FTYPE bsq=0.0;
5859  bsq_calc(prfix,ptrgeom,&bsq);
5860 
5861 
5862  if(DAMPGAMMA){
5863  if(V[1]>t-100.0 && steppart==0 && nstep%2==0){
5864  FTYPE gamma,qsq;
5865  gamma_calc(prfix,ptrgeom,&gamma,&qsq);
5866  FTYPE gammanew=((gamma-1.0)*0.1)+1.0;
5867  limit_gamma(0,gammanew,OUTERDEATHGAMMAMAXRAD,prfix,NULL,ptrgeom,-1);
5868  //prfix[U1]=0.0;
5869  //dualfprintf(fail_file,"problem: gamma=%g gammanew=%g\n",gamma,gammanew);
5870  }
5871  }
5872 
5873  if(DAMPECOV){
5874  if(V[1]>2.5 && steppart==0 && nstep%2==0){
5875  extern int init_postvpot(int i, int j, int k, FTYPE *pr, FTYPE *pstag, FTYPE *ucons);
5876  init_postvpot(i, j, k, prfix, NULL, NULL);
5877  }
5878  }
5879 
5880  if(OUTERDEATHDAMP){
5881 
5882  FTYPE Rbadout=OUTERDEATHRADIUS; //MIN(500.0,0.5*Rout);
5883 
5884  Rbadout=100.0 + t;
5885 
5886  if(V[1]>Rbadout || V[1]>t && (mycpupos[2]==0 && (j<=2) || mycpupos[2]==ncpux2-1 && (j>=N2-1-2)) ){
5887 
5888  if(VARTOINTERP==PRIMTOINTERP_GDETFULLVERSION_WALD){
5889 
5890  FTYPE prfloor[NPR],prceiling[NPR];
5891  set_density_floors(ptrgeom,prfix,prfloor,prceiling);
5892 
5893  prfix[RHO]=prfloor[RHO];
5894  prfix[UU]=prfloor[UU];
5895 
5896  limit_gamma(0,OUTERDEATHGAMMAMAX,OUTERDEATHGAMMAMAXRAD,prfix,NULL,ptrgeom,-1);
5897 
5898  }
5899  else{
5900 
5901 
5902  //prfix[RHO] = 1E-10*pow(V[1]/500.0,-1.5);
5903 
5904  prfix[UU]= MIN(50.0*prfix[RHO],prfix[UU]); // no more than u/rho=10
5905  // prfix[UU] = MIN(prfix[UU],
5906  if(DOENOFLUX != NOENOFLUX){
5907  ufix[UU]=MAX(-prfix[UU],ufix[UU]);
5908  ufix[ENTROPY] = ufix[UU];
5909  ufix[ENTROPY] = MAX(0.0001,MIN(ufix[ENTROPY],1.0)); // like u/rho=1
5910  }
5911 
5912  if(URAD0>=0){
5913  //prfix[URAD0] = 2E-10*pow(V[1]/500.0,-1.5);
5914  prfix[URAD0] = MIN(MIN(prfix[RHO],prfix[URAD0]),prfix[UU]); // no more than Erf/rho=1 and Erf/u=1
5915  if(DOENOFLUX != NOENOFLUX){
5916  ufix[URAD0]=MAX(-prfix[URAD0],ufix[URAD0]);
5917  }
5918  }
5919 
5920  // limit_gamma(0,1.5,GAMMAMAXRAD,prfix,NULL,ptrgeom,-1);
5921  limit_gamma(0,OUTERDEATHGAMMAMAX,OUTERDEATHGAMMAMAXRAD,prfix,NULL,ptrgeom,-1);
5922 
5923  }
5924  }
5925 
5926 
5927  if(VARTOINTERP==PRIMTOINTERP_GDETFULLVERSION_WALD){
5928 
5929  }
5930  else{
5931 
5932  FTYPE rbr=500.0;
5933  FTYPE localrbr=rbr; //500.0; // rbr;
5934  FTYPE localrbrr0=MAX(100.0,localrbr/2.0);
5935 
5936  FTYPE switch0 = 0.5+1.0/M_PI*atan((V[1]-localrbr)/localrbrr0); // large r
5937  FTYPE switch2 = 1.0-switch0; // small r
5938  // FTYPE myhslope1=h0 + pow( (V[1]-rsjet3)/r0jet3 , njet);
5939  // FTYPE myhslope2=h0 + pow( (localrbr-rsjet3)/r0jet3 , njet);
5940  // myhslope = myhslope1*switch2 + myhslope2*switch0;
5941 
5942  if(V[1]<100.0) switch2=1.0;
5943  //if(V[1]>500.0) switch2=0.0;
5944 
5945  if(prfix[U1]<0.0){
5946  // switch2=0.0;
5947  prfix[U1]*=switch2;
5948  if(DOENOFLUX != NOENOFLUX){
5949  ufix[U1]*=switch2;
5950  }
5951  }
5952  if(DOENOFLUX != NOENOFLUX){
5953  // ufix[U2]=ufix[U3]=0.0;
5954  // SLOOPA(jjj) ufix[U1+jjj-1]=prfix[U1+jjj-1] = 0.0;
5955  }
5956 
5957  if(URAD0>=0){
5958  if(prfix[URAD1]<0.0){
5959  //switch2=0.0;
5960  prfix[URAD1]*=switch2;
5961  if(DOENOFLUX != NOENOFLUX){
5962  ufix[URAD1]*=switch2;
5963  }
5964  }
5965  }
5966  }
5967  }// end OUTERDEATHDAMP
5968  }// end FULLLOOP
5969  }
5970  }// end if(0)
5971  else{
5972 
5973  FTYPE prfixtry[NPR],ufixtry[NPR];
5974  FTYPE prfixtryb[NPR],ufixtryb[NPR];
5975  int pliter,pl;
5976 
5977 
5978  FULLLOOP{
5979  prfix=&MACP0A1(prim,i,j,k,0);//&GLOBALMACP0A1(pglobal,i,j,k,0);
5980  ufix=&GLOBALMACP0A1(unewglobal,i,j,k,0);
5981 
5982  PLOOP(pliter,pl) prfixtry[pl]=prfix[pl];
5983  if(DOENOFLUX != NOENOFLUX){
5984  PLOOP(pliter,pl) ufixtry[pl]=ufix[pl];
5985  }
5986 
5987  PLOOP(pliter,pl) prfixtryb[pl]=prfix[pl];
5988  if(DOENOFLUX != NOENOFLUX){
5989  PLOOP(pliter,pl) ufixtryb[pl]=ufix[pl];
5990  }
5991 
5992  // get geometry for center pre-interpolated values
5993  get_geometry(i, j, k, CENT, ptrgeom);
5994  bl_coord_ijk_2(i,j,k,CENT,X, V) ;
5995 
5996  FTYPE bsq=0.0;
5997  bsq_calc(prfix,ptrgeom,&bsq);
5998 
5999 #define BEGINDEATHTRANSITION (OUTERDEATHRADIUS*0.8)
6000 
6001  if(V[1]>OUTERDEATHRADIUS|| V[1]>BEGINDEATHTRANSITION){
6002 
6003  FTYPE rbr=OUTERDEATHRADIUS;
6004  FTYPE localrbr=rbr; //500.0; // rbr;
6005  FTYPE localrbrr0=BEGINDEATHTRANSITION;
6006  FTYPE localrbrdr=localrbr/3.0;
6007 
6008 #define cr(x) (exp(-1.0/(x)))
6009 #define tr(x) (cr(x)/(cr(x) + cr(1.0-(x))))
6010 #define trans(x,L,R) ((x)<=(L) ? 0.0 : ((x)>=(R) ? 1.0 : tr(((x)-(L))/((R)-(L)))) )
6011 
6012  FTYPE switch0,switch2;
6013  FTYPE switch0b,switch2b;
6014  FTYPE rinner;
6015 
6016  if(0){
6017  switch0 = 0.5+1.0/M_PI*atan((V[1]-localrbr)/localrbrdr); // large r
6018  switch2 = 1.0-switch0; // small r
6019 
6020  switch0b = 0.5+1.0/M_PI*atan((V[1]-localrbr)/localrbrdr); // large r
6021  switch2b = 1.0-switch0; // small r
6022  }
6023  else{
6024  rinner=MAX(rbr-localrbrdr,BEGINDEATHTRANSITION);
6025  switch0 = trans(V[1],rinner,rbr+localrbrdr);
6026  switch2 = 1.0-switch0;
6027 
6028 #define FACTORNEXTSWTICH (2.0)
6029  //rinner=MAX(FACTORNEXTSWTICH*rbr-localrbrdr,BEGINDEATHTRANSITION);
6030  //switch0b = trans(V[1],rinner,FACTORNEXTSWTICH*rbr+localrbrdr);
6031  switch0b = trans(V[1],rbr+localrbrdr,FACTORNEXTSWTICH*rbr+localrbrdr);
6032  switch2b = 1.0-switch0b;
6033  }
6034 
6035  if(V[1]<BEGINDEATHTRANSITION){
6036  switch0=0.0;
6037  switch2=1.0;
6038 
6039  switch0b=0.0;
6040  switch2b=1.0;
6041  }
6042  //if(V[1]>500.0) switch2=0.0;
6043 
6044 
6045 
6046 
6048 
6049  // FTYPE gamma,qsq;
6050  // gamma_calc(prfixtry,ptrgeom,&gamma,&qsq);
6051  // FTYPE gammanew=OUTERDEATHGAMMAMAX;
6052  // limit_gamma(0,gammanew,OUTERDEATHGAMMAMAXRAD,prfixtry,NULL,ptrgeom,-1);
6053  //prfixtry[U1]=0.0;
6054  //dualfprintf(fail_file,"problem: gamma=%g gammanew=%g\n",gamma,gammanew);
6055 
6056 
6057  //prfixtry[RHO] = 1E-10*pow(V[1]/500.0,-1.5);
6058 
6059  prfixtry[UU]= MIN(1.0*prfixtry[RHO],prfixtry[UU]); // no more than u/rho=1
6060  // prfixtry[UU] = MIN(prfixtry[UU],
6061  if(DOENOFLUX != NOENOFLUX){
6062  ufixtry[UU]=MAX(-prfixtry[UU],ufixtry[UU]);
6063  ufixtry[ENTROPY] = ufixtry[UU];
6064  ufixtry[ENTROPY] = MAX(0.0001,MIN(ufixtry[ENTROPY],1.0)); // like u/rho=1
6065  }
6066 
6067  if(URAD0>=0){
6068  //prfixtry[URAD0] = 2E-10*pow(V[1]/500.0,-1.5);
6069  //prfixtry[URAD0] = MIN(MIN(prfixtry[RHO],prfixtry[URAD0]),prfixtry[UU]); // no more than Erf/rho=1 and Erf/u=1
6070  prfixtry[URAD0] = MIN(10.0*prfixtry[RHO],prfixtry[URAD0]); // no more than Erf/rho=10. No limit on u/prad0 that can be very large in optically thin or thick parts
6071  if(DOENOFLUX != NOENOFLUX){
6072  ufixtry[URAD0]=MAX(-prfixtry[URAD0],ufixtry[URAD0]);
6073  }
6074  }
6075 
6076  // limit_gamma(0,1.5,GAMMAMAXRAD,prfixtry,NULL,ptrgeom,-1);
6077  limit_gamma(0,OUTERDEATHGAMMAMAX,OUTERDEATHGAMMAMAXRAD,prfixtry,NULL,ptrgeom,-1);
6078 
6079 
6080 
6081 #if(1)
6082  if(prfixtry[U1]<0.0){
6083  prfixtry[U1]=0.0;
6084  if(DOENOFLUX != NOENOFLUX){
6085  ufixtry[U1]=0.0;
6086  }
6087  }
6088  if(DOENOFLUX != NOENOFLUX){
6089  // ufixtry[U2]=ufixtry[U3]=0.0;
6090  // SLOOPA(jjj) ufixtry[U1+jjj-1]=prfixtry[U1+jjj-1] = 0.0;
6091  }
6092 
6093  if(URAD0>=0){
6094  if(prfixtry[URAD1]<0.0){
6095  prfixtry[URAD1]=0.0;
6096  if(DOENOFLUX != NOENOFLUX){
6097  ufixtry[URAD1]=0.0;
6098  }
6099  }
6100  }
6101 #endif
6102 
6104  // set b
6105  PLOOP(pliter,pl){
6106  prfixtryb[pl] = prfixtry[pl];
6107  if(DOENOFLUX != NOENOFLUX){
6108  ufixtryb[pl] = ufixtry[pl];
6109  }
6110  }
6111 
6112 #if(1)
6113  // stronger restrictions for b on prad0
6114  if(URAD0>=0){
6115  //FTYPE Tgaslocal=compute_temp_simple(i,j,k,CENT,prfixtryb[RHO],prfixtryb[UU]);
6116  // prfixtryb[URAD0] = MIN(MIN(prfixtryb[RHO],prfixtryb[URAD0]),10.0*prfixtryb[UU]); // no more than Erf/rho=1 and Erf/u=10
6117  prfixtryb[URAD0] = MIN(prfixtryb[RHO],prfixtryb[URAD0]);
6118  // prfixtryb[URAD0] = MAX(prfixtryb[URAD0],prfixtryb[UU]); // avoid drop-outs in Erf
6119  // try thermal equilibrium (arad Tgas^4 = urad )
6120  // prfixtryb[URAD0] = calc_LTE_EfromT(Tgaslocal);
6121 
6122  if(DOENOFLUX != NOENOFLUX){
6123  ufixtryb[URAD0]=MAX(-prfixtryb[URAD0],ufixtryb[URAD0]);
6124  }
6125  }
6126 #endif
6127 
6128  // final setting of fixed up quantities
6129  PLOOP(pliter,pl){
6130  prfix[pl] = prfixtryb[pl]*switch0b + prfixtry[pl]*switch0*switch2b + prfix[pl]*switch2*switch2b;
6131  if(DOENOFLUX != NOENOFLUX){
6132  ufix[pl] = ufixtryb[pl]*switch0b + ufixtry[pl]*switch0*switch2b + ufix[pl]*switch2*switch2b;
6133  }
6134  }
6135 
6136  // dualfprintf(fail_file,"V[1]=%g ijk=%d %d %d switch0=%g switch2=%g\n",V[1],i,j,k,switch0,switch2);
6137 
6138  } // end if beyond some radius
6139  }
6140 
6141  }
6142 
6143 }
6144 
6145 
6146 
6148 void debugfixup(void)
6149 {
6150  int i,j,k;
6151  // struct of_geom geomdontuse;
6152  // struct of_geom *ptrgeom=&geomdontuse;
6153  FTYPE X[NDIM],V[NDIM];
6154  FTYPE *prfix,*ufix;
6155  int jjj;
6156  FULLLOOP{
6157  prfix=&GLOBALMACP0A1(pglobal,i,j,k,0);
6158  ufix=&GLOBALMACP0A1(unewglobal,i,j,k,0);
6159 
6160  // get geometry for center pre-interpolated values
6161  //get_geometry(i, j, k, CENT, ptrgeom);
6162  bl_coord_ijk_2(i,j,k,CENT,X, V) ;
6163 
6164  if(V[1]>300.0 && V[2]<0.075){
6165 
6166  prfix[RHO] = RHOMIN*pow(V[1],-2.0);
6167 
6168  prfix[UU]= MIN(prfix[RHO],prfix[UU]); // no more than u/rho=1
6169  prfix[URAD0] = MIN(MIN(prfix[RHO],prfix[URAD0]),prfix[UU]); // no more than Erf/rho=1 and Erf/u=1
6170 
6171  if(DOENOFLUX != NOENOFLUX){
6172  ufix[UU]=MAX(-prfix[UU],ufix[UU]);
6173  ufix[ENTROPY] = ufix[UU];
6174  ufix[ENTROPY] = MAX(0.0001,MIN(ufix[ENTROPY],1.0)); // like u/rho=1
6175  ufix[URAD0]=MAX(-prfix[URAD0],ufix[URAD0]);
6176  }
6177 
6178  SLOOPA(jjj) prfix[U1+jjj-1] = 0.0;
6179  SLOOPA(jjj) prfix[URAD1+jjj-1] = 0.0;
6180 
6181  if(DOENOFLUX != NOENOFLUX){
6182  SLOOPA(jjj) ufix[U1+jjj-1]= 0.0;
6183  SLOOPA(jjj) ufix[URAD1+jjj-1]=0.0;
6184  }
6185  }
6186  }
6187 
6188 
6189 }
6190 
6192 void debugfixupnonmad(void)
6193 {
6194  // hack to get rid of bad region at large distances when restarted
6195  int i,j,k;
6196  // struct of_geom geomdontuse;
6197  // struct of_geom *ptrgeom=&geomdontuse;
6198  FTYPE X[NDIM],V[NDIM];
6199  FTYPE *prfix,*ufix;
6200  int jjj;
6201  FULLLOOP{
6202  prfix=&GLOBALMACP0A1(pglobal,i,j,k,0);
6203  ufix=&GLOBALMACP0A1(unewglobal,i,j,k,0);
6204 
6205  // get geometry for center pre-interpolated values
6206  //get_geometry(i, j, k, CENT, ptrgeom);
6207  bl_coord_ijk_2(i,j,k,CENT,X, V) ;
6208 
6209  // if(V[1]>300.0 && V[2]<0.075){
6210  // if(V[1]>500.0 && V[2]<0.075){
6211  //if(V[1]>500.0 && V[2]>M_PI-0.08){
6212  // if(V[1]>800.0 && V[2]>M_PI-0.08){
6213  // if(V[1]>800.0 && V[2]>M_PI*0.5*1.1){
6214  if(V[1]>1E3 && V[2]>M_PI*0.5*1.1){
6215 
6216  prfix[RHO] = RHOMIN*pow(V[1],-2.0);
6217 
6218  prfix[UU]= MIN(prfix[RHO],prfix[UU]); // no more than u/rho=1
6219  if(DOENOFLUX != NOENOFLUX){
6220  ufix[UU]=MAX(-prfix[UU],ufix[UU]);
6221  ufix[ENTROPY] = ufix[UU];
6222  ufix[ENTROPY] = MAX(0.0001,MIN(ufix[ENTROPY],1.0)); // like u/rho=1
6223  }
6224 
6225  prfix[URAD0] = MIN(MIN(prfix[RHO],prfix[URAD0]),prfix[UU]); // no more than Erf/rho=1 and Erf/u=1
6226  if(DOENOFLUX != NOENOFLUX){
6227  ufix[URAD0]=MAX(-prfix[URAD0],ufix[URAD0]);
6228  }
6229 
6230  SLOOPA(jjj) prfix[U1+jjj-1] = 0.0;
6231  SLOOPA(jjj) prfix[URAD1+jjj-1] = 0.0;
6232 
6233  if(DOENOFLUX != NOENOFLUX){
6234  SLOOPA(jjj) ufix[U1+jjj-1]= 0.0;
6235  SLOOPA(jjj) ufix[URAD1+jjj-1]=0.0;
6236  }
6237  }
6238  }
6239 
6240 
6241 }
6242 
6244 void debugfixupmad(void)
6245 {
6246  // hack to get rid of bad region at large distances when restarted
6247  int i,j,k;
6248  struct of_geom geomdontuse;
6249  struct of_geom *ptrgeom=&geomdontuse;
6250  FTYPE X[NDIM],V[NDIM];
6251  FTYPE *prfix,*ufix;
6252  int jjj;
6253  FULLLOOP{
6254  prfix=&GLOBALMACP0A1(pglobal,i,j,k,0);
6255  ufix=&GLOBALMACP0A1(unewglobal,i,j,k,0);
6256 
6257  // get geometry for center pre-interpolated values
6258  get_geometry(i, j, k, CENT, ptrgeom);
6259  bl_coord_ijk_2(i,j,k,CENT,X, V) ;
6260 
6261  FTYPE bsq=0.0;
6262  bsq_calc(prfix,ptrgeom,&bsq);
6263 
6264 
6265  if(V[1]>60.0){
6266  limit_gamma(0,1.5,1.5,prfix,NULL,ptrgeom,-1);
6267  }
6268 
6269  if(V[1]>60.0 && bsq/prfix[RHO]>1.0){
6270  prfix[UU]= MIN(prfix[RHO],prfix[UU]); // no more than u/rho=1
6271 
6272  if(DOENOFLUX != NOENOFLUX){
6273  ufix[UU]=MAX(-prfix[UU],ufix[UU]);
6274  ufix[ENTROPY] = ufix[UU];
6275  ufix[ENTROPY] = MAX(0.0001,MIN(ufix[ENTROPY],1.0)); // like u/rho=1
6276  }
6277 
6278  prfix[URAD0] = MIN(MIN(prfix[RHO],prfix[URAD0]),prfix[UU]); // no more than Erf/rho=1 and Erf/u=1
6279  if(DOENOFLUX != NOENOFLUX){
6280  ufix[URAD0]=MAX(-prfix[URAD0],ufix[URAD0]);
6281  }
6282 
6283  limit_gamma(0,1.5,GAMMAMAXRAD,prfix,NULL,ptrgeom,-1);
6284 
6285  }
6286 
6287  if(V[1]>900.0 && (V[2]>M_PI*0.5*1.1 || V[2]<M_PI*0.5*0.9) ){
6288 
6289  prfix[RHO] = 1.5E-9*pow(V[1]/500.0,-1.5);
6290 
6291  prfix[UU]= MIN(prfix[RHO],prfix[UU]); // no more than u/rho=1
6292  if(DOENOFLUX != NOENOFLUX){
6293  ufix[UU]=MAX(-prfix[UU],ufix[UU]);
6294  ufix[ENTROPY] = ufix[UU];
6295  ufix[ENTROPY] = MAX(0.0001,MIN(ufix[ENTROPY],1.0)); // like u/rho=1
6296  }
6297 
6298 
6299  prfix[URAD0] = 1E-9*pow(V[1]/900.0,-1.0);
6300  prfix[URAD0] = MIN(MIN(prfix[RHO],prfix[URAD0]),prfix[UU]); // no more than Erf/rho=1 and Erf/u=1
6301  if(DOENOFLUX != NOENOFLUX){
6302  ufix[URAD0]=MAX(-prfix[URAD0],ufix[URAD0]);
6303  }
6304 
6305  limit_gamma(0,1.5,GAMMAMAXRAD,prfix,NULL,ptrgeom,-1);
6306 
6307  }
6308  }
6309 
6310 
6311 }
6312 
6314 void debugfixupalt1(void)
6315 {
6316  // hack to get rid of bad region at large distances when restarted
6317  int i,j,k;
6318  struct of_geom geomdontuse;
6319  struct of_geom *ptrgeom=&geomdontuse;
6320  FTYPE X[NDIM],V[NDIM];
6321  FTYPE *prfix,*ufix;
6322  int jjj;
6323  FULLLOOP{
6324  prfix=&GLOBALMACP0A1(pglobal,i,j,k,0);
6325  ufix=&GLOBALMACP0A1(unewglobal,i,j,k,0);
6326 
6327  // get geometry for center pre-interpolated values
6328  get_geometry(i, j, k, CENT, ptrgeom);
6329  bl_coord_ijk_2(i,j,k,CENT,X, V) ;
6330 
6331  FTYPE bsq=0.0;
6332  bsq_calc(prfix,ptrgeom,&bsq);
6333 
6334  if(0){
6335  if(V[1]>60.0){
6336  limit_gamma(0,1.5,GAMMAMAXRAD,prfix,NULL,ptrgeom,-1);
6337 
6338  prfix[UU]= MIN(prfix[RHO],prfix[UU]); // no more than u/rho=1
6339  if(DOENOFLUX != NOENOFLUX){
6340  ufix[UU]=MAX(-prfix[UU],ufix[UU]);
6341  }
6342  }
6343  }
6344 
6345  if(0){
6346  if(V[1]>60.0 && bsq/prfix[RHO]>1.0){
6347  prfix[UU]= MIN(prfix[RHO],prfix[UU]); // no more than u/rho=1
6348 
6349  if(DOENOFLUX != NOENOFLUX){
6350  ufix[UU]=MAX(-prfix[UU],ufix[UU]);
6351  ufix[ENTROPY] = ufix[UU];
6352  ufix[ENTROPY] = MAX(0.0001,MIN(ufix[ENTROPY],1.0)); // like u/rho=1
6353  }
6354 
6355  prfix[URAD0] = MIN(MIN(prfix[RHO],prfix[URAD0]),prfix[UU]); // no more than Erf/rho=1 and Erf/u=1
6356  if(DOENOFLUX != NOENOFLUX){
6357  ufix[URAD0]=MAX(-prfix[URAD0],ufix[URAD0]);
6358  }
6359 
6360  limit_gamma(0,1.5,GAMMAMAXRAD,prfix,NULL,ptrgeom,-1);
6361 
6362  }
6363  }
6364  if(0){
6365  // if(V[1]>600.0 && (V[2]>M_PI*0.5*1.1 || V[2]<M_PI*0.5*0.9) ){
6366  if(V[1]>600.0 ){
6367 
6368  //prfix[RHO] = 1E-10*pow(V[1]/500.0,-1.5);
6369 
6370  prfix[UU]= MIN(prfix[RHO],prfix[UU]); // no more than u/rho=1
6371  if(DOENOFLUX != NOENOFLUX){
6372  ufix[UU]=MAX(-prfix[UU],ufix[UU]);
6373  ufix[ENTROPY] = ufix[UU];
6374  ufix[ENTROPY] = MAX(0.0001,MIN(ufix[ENTROPY],1.0)); // like u/rho=1
6375  }
6376 
6377  //prfix[URAD0] = 2E-10*pow(V[1]/500.0,-1.5);
6378  prfix[URAD0] = MIN(MIN(prfix[RHO],prfix[URAD0]),prfix[UU]); // no more than Erf/rho=1 and Erf/u=1
6379  if(DOENOFLUX != NOENOFLUX){
6380  ufix[URAD0]=MAX(-prfix[URAD0],ufix[URAD0]);
6381  }
6382 
6383  limit_gamma(0,1.5,GAMMAMAXRAD,prfix,NULL,ptrgeom,-1);
6384 
6385  if(DOENOFLUX != NOENOFLUX){
6386  if(ufix[U1]<0.0) ufix[U1]=0.0;
6387  ufix[U2]=ufix[U3]=0.0;
6388  if(ufix[URAD1]<0.0) ufix[URAD1]=0.0;
6389  ufix[URAD2]=ufix[URAD3]=0.0;
6390  }
6391 
6392  }
6393  }
6394  }
6395 
6396 
6397 }
6398 
6399 
6401 void debugfixupaltdeath(void)
6402 {
6403  // hack to get rid of bad region at large distances when restarted
6404  int i,j,k;
6405  struct of_geom geomdontuse;
6406  struct of_geom *ptrgeom=&geomdontuse;
6407  FTYPE X[NDIM],V[NDIM];
6408  FTYPE *prfix,*ufix;
6409  int jjj;
6410  FULLLOOP{
6411  prfix=&GLOBALMACP0A1(pglobal,i,j,k,0);
6412  ufix=&GLOBALMACP0A1(unewglobal,i,j,k,0);
6413 
6414  // get geometry for center pre-interpolated values
6415  get_geometry(i, j, k, CENT, ptrgeom);
6416  bl_coord_ijk_2(i,j,k,CENT,X, V) ;
6417 
6418  FTYPE bsq=0.0;
6419  bsq_calc(prfix,ptrgeom,&bsq);
6420 
6421  if(0){
6422  if(V[1]>60.0){
6423  limit_gamma(0,1.5,GAMMAMAXRAD,prfix,NULL,ptrgeom,-1);
6424 
6425  prfix[UU]= MIN(prfix[RHO],prfix[UU]); // no more than u/rho=1
6426  if(DOENOFLUX != NOENOFLUX){
6427  ufix[UU]=MAX(-prfix[UU],ufix[UU]);
6428  }
6429  }
6430  }
6431 
6432  if(0){
6433  if(V[1]>60.0 && bsq/prfix[RHO]>1.0){
6434  prfix[UU]= MIN(prfix[RHO],prfix[UU]); // no more than u/rho=1
6435 
6436  if(DOENOFLUX != NOENOFLUX){
6437  ufix[UU]=MAX(-prfix[UU],ufix[UU]);
6438  ufix[ENTROPY] = ufix[UU];
6439  ufix[ENTROPY] = MAX(0.0001,MIN(ufix[ENTROPY],1.0)); // like u/rho=1
6440  }
6441 
6442  prfix[URAD0] = MIN(MIN(prfix[RHO],prfix[URAD0]),prfix[UU]); // no more than Erf/rho=1 and Erf/u=1
6443  if(DOENOFLUX != NOENOFLUX){
6444  ufix[URAD0]=MAX(-prfix[URAD0],ufix[URAD0]);
6445  }
6446 
6447  limit_gamma(0,1.5,GAMMAMAXRAD,prfix,NULL,ptrgeom,-1);
6448 
6449  }
6450  }
6451  if(1){
6452  // if(V[1]>600.0 && (V[2]>M_PI*0.5*1.1 || V[2]<M_PI*0.5*0.9) ){
6453  if(V[1]>1E3 ){
6454 
6455  //prfix[RHO] = 1E-10*pow(V[1]/500.0,-1.5);
6456 
6457  prfix[UU]= MIN(prfix[RHO],prfix[UU]); // no more than u/rho=1
6458 
6459  if(DOENOFLUX != NOENOFLUX){
6460  ufix[UU]=MAX(-prfix[UU],ufix[UU]);
6461  ufix[ENTROPY] = ufix[UU];
6462  ufix[ENTROPY] = MAX(0.0001,MIN(ufix[ENTROPY],1.0)); // like u/rho=1
6463  }
6464 
6465  //prfix[URAD0] = 2E-10*pow(V[1]/500.0,-1.5);
6466  prfix[URAD0] = MIN(MIN(prfix[RHO],prfix[URAD0]),prfix[UU]); // no more than Erf/rho=1 and Erf/u=1
6467  if(DOENOFLUX != NOENOFLUX){
6468  ufix[URAD0]=MAX(-prfix[URAD0],ufix[URAD0]);
6469  }
6470 
6471  // limit_gamma(0,1.5,GAMMAMAXRAD,prfix,NULL,ptrgeom,-1);
6472  limit_gamma(0,6.0,GAMMAMAXRAD,prfix,NULL,ptrgeom,-1);
6473 
6474  if(DOENOFLUX != NOENOFLUX){
6475  if(ufix[U1]<0.0) ufix[U1]=0.0;
6476  // ufix[U2]=ufix[U3]=0.0;
6477 
6478  if(ufix[URAD1]<0.0) ufix[URAD1]=0.0;
6479  // ufix[URAD2]=ufix[URAD3]=0.0;
6480  }
6481 
6482  }
6483  }
6484  }
6485 
6486 
6487 }
6488