HARM
harm and utilities
 All Data Structures Files Functions Variables Typedefs Macros Pages
metric.tools.c
Go to the documentation of this file.
1 
7 #include "decs.h"
8 
9 
10 
11 #if(0) // uniformly low values although not always lower than original version
12 #define MAXITER 5
13 #define NRANSI
14 #define CON 1.1
15 #define CON2 (CON*CON)
16 //#define NTAB 130 // number of function evaluations is 2XNTAB
17 #define NTAB 30 // number of function evaluations is 2XNTAB
18 
19 #define SAFE 2.0
20 #define TRYTOL (trytollocal) // attempted tolerance
21 #define OKTOL 1e-5 // error must be below this to avoid report
22 #define FAILTOL 1e-1
23 //#define HSTARTCHANGE 10.0
24 #endif
25 
26 #if(1) // original version (gets pretty damn low for many, but not all, derivatives -- some 1E-6 rather than 1E-13)
27 #define MAXITER 15 //Increased by Sasha up from 5 to avoid non-convergence errors (>1e-5)
28 #define NRANSI
29 #define CON 1.3
30 #define CON2 (CON*CON)
31 #define NTAB 10 // number of function evaluations is 2XNTAB
32 
33 #define SAFE 2.0
34 #define TRYTOL 1E-10 // attempted tolerance
35 #define OKTOL 1e-5 // error must be below this to avoid report
36 #define FAILTOL 1e-3
37 #endif
38 
39 // whether to turn on extensive recent debugging
40 #define DEBUGDF 0
41 
42 // whether to output matrix of derivative and extrapolations of various orders
43 #define DEBUGOUTPUTAMATRIX 0
44 
45 // whether to use Jon's additional starting "h" checks [about 2X-3X more expensive in typical cases]
46 #define USEJONEXTENSION 1
47 
48 
49 
50 // jon's version of NR's dfridr modified to accept more general, needed, function
51 int dfridr(FTYPE (*func)(struct of_geom *, FTYPE*,int,int), struct of_geom *ptrgeom, FTYPE *X,int ii, int jj, int kk, FTYPE *finalanswer)
52 {
53  int i,j,k;
54  FTYPE errt,fac,hh,**a,ans;
55  FTYPE dX[NDIM],Xh[NDIM],Xl[NDIM];
56  FTYPE h,err;
57  FTYPE firsthstart,hstart;
58  FTYPE newdx,temp;
59  int iter;
60  FTYPE errlist[MAXITER];
61  FTYPE hhlist[MAXITER];
62  FTYPE anslist[MAXITER];
63  FTYPE minerror,minerrorhstart;
64  int miniter;
65  int iterdebug;
66  FTYPE minans;
67  FTYPE trytollocal;
68  int lasti;
69  int Nvec[NDIM];
70  int goodi,goodj;
71  FTYPE goodhh;
72  FTYPE truecon,truecon2;
73  int nrlasti,nrgoodi,nrgoodj;
74  FTYPE nrerr,nrans,nrgoodhh;
75  int iterfailed;
76 
77 
78  trytollocal=NUMEPSILONPOW23;
79  // minhlocal=NUMEPSILONPOW23;
80 
81  // allocate memory
82  a=dmatrix(1,NTAB,1,NTAB);
83 
84  // starting delta shouldn't be so large to cross interesting boundaries (like pole or r=0), but should be large enough to allow convergence to answer at smaller intervals.
85  Nvec[0]=0;
86  Nvec[1]=N1;
87  Nvec[2]=N2;
88  Nvec[3]=N3;
89  if(kk==TT) hstart=CONNDELTA;
90  else if(Nvec[kk]==1) hstart=CONNDELTA;
91  else hstart=0.5*dx[kk];
92 
93 
94 
95  firsthstart=hstart;
96  miniter=0; // didn't find minimum is assumed, which means first is minimum!
97  iter=0;
98  minerror=BIG;
99  minans=BIG;
100  truecon=CON;
101  truecon2=CON2;
102  goodhh=-1.0; // -1.0 indicates if ever found goodhh
103 
104 
106  //
107  // START BIG LOOP over Jon's iterations
108  //
110  iterfailed=0;
111  while(1){
112 
113 
114 
115 #if(DEBUGDF)
116  dualfprintf(fail_file,"iter=%d\n",iter);
117 #endif
118 
119  // Set initial differential size
120  h=hstart;
121  if (h <= 0.0) nrerror("h must be positive definite in dfridr.");
122  hh=h;
123 
124 
125  // HARM STUFF
126  for(k=0;k<NDIM;k++) dX[k]=0.0; // other components will remains 0 for this function
127  dX[kk]=hh;
128  for(k=0;k<NDIM;k++){
129  // Xl[k]=X[k]-dX[k];
130  temp=X[k]-dX[k];
131  newdx=temp-X[k];
132  Xl[k]=X[k]+newdx;
133  // dualfprintf(fail_file,"newdx[k=%d] for low = %21.15g : %21.15g\n",k,newdx,Xl[k]);
134  }
135  for(k=0;k<NDIM;k++){
136  // Xh[k]=X[k]+dX[k];
137  temp=X[k]+dX[k];
138  newdx=temp-X[k];
139  Xh[k]=X[k]+newdx;
140  // dualfprintf(fail_file,"newdx[k=%d] for high = %21.15g : %21.15g\n",k,newdx,Xh[k]);
141  }
142  // dualfprintf(fail_file,"a[1][1]=%21.15g dX[kk=%d]=%21.15g\n",a[1][1],kk,dX[kk]);
143  // for(k=0;k<NDIM;k++) dualfprintf(fail_file,"X[%d]=%21.15g Xh[%d]=%21.15g Xl[%d]=%21.15g dXhl=%21.15g\n",k,X[k],k,Xh[k],k,Xl[k],Xh[k]-Xl[k]);
144  // for(k=0;k<NDIM;k++) dualfprintf(fail_file,"funch=%21.15g funcl=%21.15g\n",k,(*func)(ptrgeom,Xh,ii,jj),k,(*func)(ptrgeom,Xl,ii,jj));
145  // end HARM STUFF
146 
147 
148  // compute df/dx
149  a[1][1]=((*func)(ptrgeom,Xh,ii,jj)-(*func)(ptrgeom,Xl,ii,jj))/(2.0*hh);
150  err=BIG;
151  for (i=2;i<=NTAB;i++) {
152  hh /= truecon;
153 
154  // HARM STUFF
155  dX[kk]=hh;
156  for(k=0;k<NDIM;k++) Xl[k]=X[k]-dX[k];
157  for(k=0;k<NDIM;k++) Xh[k]=X[k]+dX[k];
158  // for(k=0;k<NDIM;k++) dualfprintf(fail_file,"i=%d k=%d %21.15g %21.15g\n",i,k,X[k],dX[k]);
159  // end HARM STUFF
160 
161  // compute df/dx with dx here as dx*dx in dx size
162  a[1][i]=((*func)(ptrgeom,Xh,ii,jj)-(*func)(ptrgeom,Xl,ii,jj))/(2.0*hh);
163  fac=truecon2;
164 
165  // compute Neville table for each smaller step size
166  // Contains extrapolation to dx->0 for each possible order of extrapolation
167  for (j=2;j<=i;j++) {
168  // extrapolate for dx->0
169  a[j][i]=(a[j-1][i]*fac-a[j-1][i-1])/(fac-1.0);
170  fac=truecon2*fac;
171  //unnormalized error
172  // errt=MAX(fabs(a[j][i]-a[j-1][i]),fabs(a[j][i]-a[j-1][i-1]));
173  //normalized error
174  // errt=MAX(fabs(a[j][i]-a[j-1][i]),fabs(a[j][i]-a[j-1][i-1]))/((*func)(ptrgeom,X,ii,jj)+SMALL);
175  // normalized error
176  errt=MAX(fabs(a[j][i]-a[j-1][i]),fabs(a[j][i]-a[j-1][i-1]));
177  errt/=MAX(MAX(MAX(fabs(a[j][i]),fabs(a[j-1][i])),MAX(fabs(a[j][i]),fabs(a[j-1][i-1]))),SMALL);
178 
179  if (errt <= err) {
180  err=errt;
181  ans=a[j][i];
182  goodj=j;
183  goodi=i;
184  goodhh=hh;
185  }
186  }// end over j derivatives
187 
188  // unnormalize error
189  // if (fabs((a[i][i]-a[i-1][i-1])) >= SAFE*(err)) break;
190  // normalized error
191  // if (fabs((a[i][i]-a[i-1][i-1])/( (*func)(ptrgeom,X,ii,jj)+SMALL)) >= SAFE*(err)) break;
192  // if (fabs((a[i][i]-a[i-1][i-1]))/(fabs(ans)+SMALL) >= SAFE*(err)) break;
193  // normalized error
194  if (fabs((a[i][i]-a[i-1][i-1]))/MAX(SMALL,MAX(fabs(a[i][i]),fabs(a[i-1][i-1]))) >= SAFE*(err)){
195  nrlasti=i;
196  nrerr=err;
197  nrans=ans;
198  nrgoodj=goodj;
199  nrgoodi=goodi;
200  nrgoodhh=goodhh;
201 #if(0) // NR breaks a bit early -- found solution can sometimes be better if wait -- so just do whole NTAB
202  break;
203 #endif
204  }// end NR early termination check
205 
206  lasti=i;
207 
208  }// end over NTAB entries for i
209 
210 
211 
212 
213 
214 
215 
216 
217 
218 #if(DEBUGOUTPUTAMATRIX) // debug report to see exactly what table contained
219  dualfprintf(fail_file,"BEGIN DFRIDR table for i=%d j=%d k=%d :: ii=%d jj=%d kk=%d hstart=%21.15g lasti=%d nrlasti=%d\n",ptrgeom->i,ptrgeom->j,ptrgeom->k,ii,jj,kk,hstart,lasti,nrlasti);
220  dualfprintf(fail_file," %21s ","ORDERS");
221  for(i=1;i<=lasti;i++){
222  dualfprintf(fail_file,"%21d",i);
223  }
224  dualfprintf(fail_file,"\n");
225  for(i=1;i<=lasti;i++){
226  dualfprintf(fail_file,"hh=%21.15g ",hstart/pow(truecon,i-1));
227  for(j=1;j<=lasti;j++){
228  if(j<=i){
229  if(i==goodi && j==goodj){
230  dualfprintf(fail_file,"*%21.15g*",a[j][i]);
231  }
232  else if(i==nrgoodi && j==nrgoodj){
233  dualfprintf(fail_file,"?%21.15g?",a[j][i]);
234  }
235  else{
236  dualfprintf(fail_file," %21.15g ",a[j][i]);
237  }
238  }
239  else{
240  dualfprintf(fail_file,"%21s"," ");
241  }
242  if(j==lasti) dualfprintf(fail_file,"\n");
243  else dualfprintf(fail_file," ");
244  }
245  if(i==lasti) dualfprintf(fail_file,"\n");
246  }
247  dualfprintf(fail_file,"final error=%21.15g goodhh=%21.15g ans=%21.15g\n",err,goodhh,ans);
248  dualfprintf(fail_file,"NR error=%21.15g goodhh=%21.15g ans=%21.15g\n",nrerr,nrgoodhh,nrans);
249  dualfprintf(fail_file,"END DFRIDR table for i=%d j=%d k=%d ii=%d jj=%d kk=%d\n",ptrgeom->i,ptrgeom->j,ptrgeom->k,ii,jj,kk);
250 #endif
251 
252 
253 
254 
255 
256 
257 
258 #if(USEJONEXTENSION==0)
259  break;
260 #else
261 
263  //
264  // Entire below section is new JCM additions
265  //
267 
269  //
270  // now check error is not crazy with the starting large h, decrease h if crazy until not crazy and get good error
271  //
273 
274  errlist[iter]=err; // store final error from NR method
275  // hhlist[iter]=hstart;
276  hhlist[iter]=goodhh; // current hh to focus around
277  anslist[iter]=ans; // store result
278 
279  if(err>TRYTOL){ // TRYTOL is error we are attempting to reach
280 
281  if(errlist[iter]<minerror){
282  // store min error event
283  minerror=errlist[iter];
284  minerrorhstart=hhlist[iter];
285  miniter=iter;
286  minans=ans;
287 #if(DEBUGDF)
288  dualfprintf(fail_file,"minerr=%21.15g minhstart=%21.15g miniter=%d minans=%21.15g\n",minerror,minerrorhstart,miniter,minans);
289 #endif
290  }
291  else{
292  // if did no better through bisecting hhgood, then probably done
293  ans=minans;
294 #if(DEBUGDF)
295  if(iter==1) dualfprintf(fail_file,"Done with no better error (iter=%d)\n",iter);
296  else dualfprintf(fail_file,"Did better on multiple iterations (iter=%d)\n",iter);
297 #endif
298  break;
299  }
300 
301  // if here, then not done yet with bisecting on hh
302  // try h around goodh, but narrow down the factor by which hh changes so don't skip over better hh than goodhh
303  // NTAB entries, and want to go back to prior hh before goodhh, but need to cross to goodhh/truecon
304  // so change on truecon is fixed to be:
305  FTYPE htop,hbottom;
306  htop=goodhh*truecon;
307  hbottom=goodhh/truecon;
308  truecon= pow(hbottom,-1.0/NTAB) * pow(htop,1.0/NTAB);
309  truecon2=truecon*truecon; // reset truecon2
310  // now set new hstart:
311  hstart=htop;
312  // so start at htop and will exponentially approach hbottom after NTAB entries.
313 #if(DEBUGDF)
314  dualfprintf(fail_file,"iter=%d err=%21.15g goodhh=%21.15g lasti=%d truecon=%21.15g htop=%21.15g hbottom=%21.15g\n",iter,err,goodhh,lasti,truecon,htop,hbottom);
315 #endif
316 
317  if(truecon==1.0){
318  ans=minans;
319  break; // can't do any better
320  }
321 
322  }// end if err > TRYTOL
323  else break; // then error<TRYTOL! GOOD!
324 
325 
326 
327 
328  iter++;
329  if(iter>=MAXITER){
330  if(err<OKTOL) break;
331  else{ // then maybe problems
332  ans=minans;
333 
334  if((minerror<OKTOL || err<OKTOL)){
335  break; // then accept as answer
336  }
337  else{
338  // then must fail
339  dualfprintf(fail_file,"iter=%d>=MAXITER=%d: never found error below %21.15g: err=%21.15g : ii=%d jj=%d kk=%d\n",iter,MAXITER,OKTOL,err,ii,jj,kk);
340  dualfprintf(fail_file,"gi=%d gj=%d gk=%d\n",ptrgeom->i,ptrgeom->j,ptrgeom->k);
341  dualfprintf(fail_file,"ti=%d tj=%d tk=%d\n",startpos[1]+ptrgeom->i,startpos[2]+ptrgeom->j,startpos[3]+ptrgeom->k);
342  dualfprintf(fail_file,"miniter=%d errlist[miniter]=%21.15g hhlist[miniter]=%21.15g\n",miniter,errlist[miniter],hhlist[miniter]);
343  dualfprintf(fail_file,"minerror=%21.15g minans=%21.15g\n",minerror,minans);
344  for(iterdebug=0;iterdebug<iter;iterdebug++){
345  dualfprintf(fail_file,"h[%d]=%21.15g err[%d]=%21.15g ans[%d]=%21.15g\n",iterdebug,hhlist[iterdebug],iterdebug,errlist[iterdebug],iterdebug,anslist[iterdebug]);
346  }
347 
348  iterfailed=1;
349  break;
350  }//end if not OKTOL
351  }// end if not OKTOL for only err
352  }// end if iter>=MAXITER
353 #endif // end JCM addition
354 
355 
356 
357  }// end while loop
358 
359 
360  // done
361  free_dmatrix(a,1,NTAB,1,NTAB);
362 
363 
364  if(err<OKTOL){
365  // then good, nothing to report
366  }
367  else{
368  if(debugfail>=2) dualfprintf(fail_file,"Bad NUMREC error at i=%d j=%d k=%d ii=%d jj=%d kk=%d error=%21.15g ans=%21.15g :: iter=%d lasti=%d minerrorhstart=%21.15g firsthstart=%21.15g\n",ptrgeom->i,ptrgeom->j,ptrgeom->k,ii,jj,kk,err,ans,iter,lasti,minerrorhstart,firsthstart);
369  }
370 
371 
372  if(err>=FAILTOL){
373  return(1); // indicate failure
374  }
375 
376 
377  *finalanswer=ans;
378  return(0); // indicate no failure
379 
380 
381 
382 }
383 #undef CON
384 #undef CON2
385 #undef BIG
386 #undef NTAB
387 #undef SAFE
388 #undef NRANSI
389 // (C) Copr. 1986-92 Numerical Recipes Software *1.@Q..
390 
391 
392 
393 
394 
395 
396 /*
397  FTYPE delta(int i, int j) { if(i == j) return(1.) ; else return(0.)
398  ; } */
399 
400 /* Minkowski metric; signature +2 */
401 /*
402  FTYPE mink(int i, int j) { if(i == j) { if(i == 0) return(-1.) ;
403  else return(1.) ; } else return(0.) ; } */
404 
405 
406 
407 
408 
409 // continuous mod function like Mathematica
410 int contmod(int x, int y)
411 {
412  int ans;
413 
414  ans = x%y;
415  if(ans<0) ans+=y;
416 
417  return(ans);
418 }
419 
420 // below 2 mysin/mycos preserve symmetry for arg=-M_PI/2 .. 2M_PI by restricting arg to 0..PI/2
421 FTYPE mysin(FTYPE th)
422 {
423  int contmod(int x, int y);
424  int nmod4;
425  FTYPE ans;
426 
427 #if(ACCURATESINCOS)
428  nmod4=contmod((int)(th/(M_PI*0.5)),4);
429 
430  switch(nmod4){
431  case 0:
432  ans=sin(th);
433  break;
434  case 1:
435  ans=sin(M_PI-th);
436  break;
437  case 2:
438  ans=-sin(th-M_PI);
439  break;
440  case 3:
441  ans=-sin(2.0*M_PI-th);
442  break;
443  default:
444  dualfprintf(fail_file,"No such case for mysin with nmod4=%d\n",nmod4);
445  ans=sqrt(-1.0);
446  myexit(206983462);
447  break;
448  }
449 #else
450  ans=sin(th);
451 #endif
452 
453  return(ans);
454 
455 }
456 
457 
458 FTYPE mycos(FTYPE th)
459 {
460 
461 
462 #if(ACCURATESINCOS)
463  return(mysin(th+M_PI*0.5));
464 #else
465  return(cos(th));
466 #endif
467 
468 }
469 
470 
471 
472 #if(SUPERLONGDOUBLE==0)
473 //#ifdef WIN32
474 // cot = 1/tan = cos/sin
475 FTYPE cot(FTYPE arg)
476 {
477 
478  if(fabs(fmod(arg,M_PI))<1E-14){
479  return(0.0); // avoid singularity - assume handled elsewhere
480  }
481  else{
482 #if(ACCURATESINCOS)
483  return(mycos(arg)/mysin(arg));
484 #else
485  return(1.0/tan(arg));
486 #endif
487  }
488 
489 }
490 //#endif
491 #endif
492 
493 FTYPE csc(FTYPE arg)
494 {
495  if(fabs(fmod(arg,M_PI))<1E-14){
496  return(0.0); // avoid singularity - assume handled elsewhere
497  }
498  else{
499 #if(ACCURATESINCOS)
500  return(1.0/mysin(arg));
501 #else
502  return(1.0/sin(arg));
503 #endif
504  }
505 }
506 
507 
508 FTYPE sec(FTYPE arg)
509 {
510  if(fabs(fmod(arg,0.5*M_PI))<1E-14){
511  return(0.0); // avoid singularity - assume handled elsewhere
512  }
513  else{
514 #if(ACCURATESINCOS)
515  return(1.0/mycos(arg));
516 #else
517  return(1.0/cos(arg));
518 #endif
519  }
520 }
521 
522 
523 
524 
525 // fromwhere==0 or 1, then assume horizoni on downside of actual location
526 // fromwhere==2, then assume horizoni on upside of actual location
527 int find_horizon(int fromwhere)
528 {
529  int i, j, k, ii;
530  FTYPE r1, r2;
531  FTYPE X[NDIM],V[NDIM];
532  int gotit;
533  FTYPE horizonvalue;
534  // called after grid is setup for all cpus
535 
536 
537  // first compute where horizon is
538 
539  // these 2 below are used prior, but not initialized otherwised on restart
540  // some calculations
541  // these 2 below are also used by find_horizon() below
542  Rhor=rhor_calc(0);
543  Risco=rmso_calc(PROGRADERISCO);
544 
545 
546  // need to find horizon and place horizoni on right-hand-side of location
547 
548 
549 
550  // was testing to make sure if held horizoni fixed that conserved mass if using conservation of baryon number as basis for masses in Mvsr and accretion of energy
551  // GODMARK DEBUG DEBUG DEBUG
552  // if(fromwhere!=2){
553  // return(0);
554  // }
555 
556 
557  // definition of horizoni must be consistent so fluxes are consistent and have conservation
558  fromwhere=2; // force to be on upside unless Rhor=0, which is caught first
559 
560 
561 
562 
563 
564  if(fromwhere==0) trifprintf("begin: find_horizon ... ");
565 
566  // find cpu column that brackets the horizon and determine the
567  // i-offset of horizon
568  // notice that only 1 CPU will get horizon since stop process once found
569  // notice that radius(horizoni) is below or equal to actual horizon radius
570 
571 
572 
573 
574  horizonvalue = Rhor;
575  horizoni = -100;
576  horizoncpupos1=-1;
577  gotit = 0;
578  for (ii = numprocs - 1; ii >= 0; ii--) { // should get done by first row
579  if (ii == myid) {
580  for (i = N1-1; i >= 0; i--) {
581 
582 
583  j = N2 / 2; // doesn't matter (spherical polar assumed)
584  k = N3 / 2; // doesn't matter (spherical polar assumed)
585  coord_ijk(i, j, k, FACE1, X);
586  bl_coord_ijk(i, j, k, FACE1, V);
587  r1=V[1];
588  coord_ijk(ip1mac(i), j, k, FACE1, X);
589  bl_coord_ijk(ip1mac(i), j, k, FACE1, V);
590  r2=V[1];
591  // looking between FACE1's r value and upper FACE1's r value, so loop is from i=N1-1..i=0
592 
593  if(ii==myid && myid==0 && i==0){
594  // special check in case horizon inside inner-most radial grid
595  if(horizonvalue<=r1 || horizonvalue<SMALL){ // GODMARK: this means horizon can't be chosen to <SMALL and mean there is a black hole there
596  // then horizon off grid or right on edge, but still ok
597  // treat as if horizon is off grid if right on edge
598  horizoni = 0;
599  horizoncpupos1=mycpupos[1];
600  break;
601  }
602  }
603 
604 
605  // if (fabs(r1 - horizonvalue) <= (r2 - r1)) { // find horizon
606  if (fromwhere!=2){
607  if(horizonvalue >= r1 && horizonvalue < r2){ // note that if strictly on r2, then next CPU should pick it up
608  horizoni = i;
609  horizoncpupos1 = mycpupos[1];
610  break;
611  }
612  }
613  else if (fromwhere==2){
614  if(horizonvalue >= r1 && horizonvalue < r2){
615  horizoni = ip1mac(i);
616  horizoncpupos1 = mycpupos[1];
617  if(horizoni>=N1){
618  horizoni=0;
619  horizoncpupos1++;
620  }
621  else{
622  // then on original CPU
623  horizoncpupos1 = mycpupos[1];
624  }
625  //dualfprintf(fail_file,"horizon: %d %d\n",horizoni,horizoncpupos1);
626  break;
627  }
628  // dualfprintf(fail_file,"horizonnot: %d %d :: %21.15g %21.15g %21.15g\n",horizoni,horizoncpupos1,r1,Rhor,r2);
629  }
630  }
631  }
632 
633  if (numprocs > 0) {
634 #if(USEMPI)
635  MPI_Bcast(&horizoni, 1, MPI_INT, MPIid[ii], MPI_COMM_GRMHD);
636  MPI_Bcast(&horizoncpupos1, 1, MPI_INT, MPIid[ii], MPI_COMM_GRMHD);
637 #endif
638  }
639  if (horizoni >= 0) gotit = 1; // can stop entire process
640 
641  // keep horizoni as relative to CPU with horizon so any CPU knows where horizon is
642  // if (mycpupos[1] != horizoncpupos1) {
643  // horizoni = -100;
644  //} // reset if not right cpu group
645  if (gotit) break;
646  }
647 
648 
649 
650 
651  if(gotit==0){
652  dualfprintf(fail_file,"Never found horizon: fromwhere=%d :: MBH=%21.15g a=%21.15g :: Rhor=%21.15g Risco=%21.15g\n",fromwhere,MBH,a,Rhor,Risco);
653  myexit(6246);
654  }
655 
656 
658  //
659  // report some information
660  if(fromwhere==0) {
661  trifprintf("horizoni: %d horizoncpupos1: %d\n", horizoni, horizoncpupos1);
662  // just a check
663  dualfprintf(log_file,"horizoni: %d mycpupos[1]: %d horizoncpupos1: %d\n", horizoni, mycpupos[1], horizoncpupos1);
664 
665  trifprintf("end: find_horizon\n");
666  }
667 
668 
669  return(0);
670 }
671 
672 
673 
674 
675 
676 
677 
678 
679 
680 
681 
682 int find_RinRout(FTYPE *localRin, FTYPE *localRout)
683 {
684  FTYPE X[NDIM],V[NDIM];
685  int whichcpu;
686 
687 
688  // assume outer radius is on outer CPU
689  // only 1 CPU needs to get this
690  whichcpu=ncpux1-1;
691 
692  if(myid==whichcpu){
693 
694  coord_ijk(N1, 0, 0, FACE1, X);
695  bl_coord_ijk(N1, 0, 0, FACE1, V);
696  *localRout=V[1];
697  }
698 
699  if (numprocs > 0) {
700 #if(USEMPI)
701  MPI_Bcast(localRout, 1, MPI_FTYPE, MPIid[whichcpu], MPI_COMM_GRMHD);
702 #endif
703  }
704 
705 
706  // assume inner radius is on inner CPU
707  // only 1 CPU needs to get this
708  whichcpu=0;
709 
710  if(myid==whichcpu){
711 
712  coord_ijk(0, 0, 0, FACE1, X);
713  bl_coord_ijk(0, 0, 0, FACE1, V);
714  *localRin=V[1];
715  }
716 
717  if (numprocs > 0) {
718 #if(USEMPI)
719  MPI_Bcast(localRin, 1, MPI_FTYPE, MPIid[whichcpu], MPI_COMM_GRMHD);
720 #endif
721  }
722 
723 
724  return(0);
725 }
726 
727 
728 
729 
730 
731 
732 
733 
734 // get dr(i=0)
735 // assumes black hole is at r=0
736 void set_drsing(void)
737 {
738  FTYPE dxdxp[NDIM][NDIM];
739  FTYPE V[NDIM],X[NDIM];
740  FTYPE dr;
741 
742  // assume BH r=0 is at inner radial boundary
743  if(mycpupos[1]==0){
744  coord_ijk(0,0,0,CENT,X);
745  bl_coord_ijk(0,0,0,CENT,V);
746  dxdxprim_ijk(0,0,0,CENT,dxdxp);
747 
748  dr = (dxdxp[1][1]*dx[1] + dxdxp[1][2]*dx[2])/10.0; // divide by 10 so doesn't dominate
749 
750  // dualfprintf(fail_file,"%21.15g %21.15g %21.15g %21.15g\n",dxdxp[1][1],dx[1],dxdxp[1][2],dx[2]);
751  }
752  else dr=0;
753 
754 #if(USEMPI)
755  MPI_Allreduce(&dr, &drsing,1, MPI_FTYPE, MPI_MAX,MPI_COMM_GRMHD);
756 #else
757  drsing=dr;
758 #endif
759 
760 }
761 
762 // get 1-D line over all CPUs that has the radius (only applicable if r(x_1) and not r(x_1,x_2)
763 void set_rvsr(void)
764 {
765  FTYPE X[NDIM],V[NDIM];
766  int i,j,k;
767  int ii;
768  FTYPE r;
769 
770 
772  //
773  // get radius for this CPU
774 
775  // initialize full rcent for all CPUs, choosing value so that can use MAX over all CPUs and get answer we want
776  GRAVLOOP(ii) rcent[ii]=-1E30;
777 
778  j=0; // assumes j==0 has same radial dependence as all other j (true if r(x_1))
779  k=0; // as with j
780  COMPLOOPFP11{
781  coord_ijk(i, j, k, CENT, X);
782  bl_coord_ijk(i, j, k, CENT, V);
783  r=V[1];
784  ii=startpos[1]+i;
785  rcent[ii] = r;
786  }
787 
788  // send information to myid=0 since only this processor needs rcent
789 #if(USEMPI)
790  MPI_Reduce(&(rcent[-N1BND]),&(rcent_tot[-N1BND]),NUMGRAVPOS,MPI_FTYPE,MPI_MAX,MPIid[0], MPI_COMM_GRMHD);
791  // MPI_Reduce(rcent,rcent_tot,ncpux1*N1,MPI_FTYPE,MPI_MAX,MPIid[0], MPI_COMM_GRMHD);
792 #else
793  GRAVLOOP(ii) rcent_tot[ii]=rcent[ii];
794 #endif
795 
796 
797 }
798 
799 
800 // determinant not simply transformed from analytic function -- so no analytic form possible yet
801 int gdet_func_metric(int whichcoord, FTYPE *V,FTYPE *gcov, FTYPE *gdet)
802 {
803  int jj,kk;
804  FTYPE generalmatrixlower[NDIM][NDIM];
805  int toreturn;
806 
807  // copy to full 2D space for NR functions
808  DLOOP(jj,kk) generalmatrixlower[jj][kk] = gcov[GIND(jj,kk)];
809 
810  toreturn=gdet_func_singcheck(whichcoord, V,generalmatrixlower,gdet);
811 
812  if(FORCEGDETPOSITIVE==1){
813  *gdet=fabs(*gdet);
814  }
815 
816 
817  return(toreturn);
818 
819 
820 
821 }
822 
823 // determinant not simply transformed from analytic function -- so no analytic form possible yet
824 int gdet_func_singcheck(int whichcoord, FTYPE *V,FTYPE (*generalmatrixlower)[NDIM], FTYPE *gdet)
825 {
826  int gdet_func_orig(int whichcoord,FTYPE (*generalmatrixlower)[NDIM], FTYPE *gdet);
827  int toreturn;
828 
829 
830  toreturn=gdet_func_orig(whichcoord, generalmatrixlower,gdet);
831 
832 
833 #if(FLIPGDETAXIS)
834  if(ISSPCMCOORD(whichcoord)){
835  if(V[2]<0.0) *gdet*=-1.0;
836  if(V[2]>M_PI) *gdet*=-1.0;
837  }
838 #endif
839 
840 
841  return(toreturn);
842 
843 }
844 
845 // find determinant in general of a metric
846 /* assumes gcov has been set first; returns determinant */
847 
848 // determinant not simply transformed from analytic function -- so no analytic form possible yet
849 int gdet_func_orig(int whichcoord, FTYPE (*generalmatrixlower)[NDIM], FTYPE *gdet)
850 {
851  FTYPE d;
852  int j, k, indx[NDIM];
853  int singfix;
854  int anglesing,centersing,truedim;
855  FTYPE finalvalue;
856 
857 
858 
859 #if(USEOPENMP)
860  // maintain thread safety
861  FTYPE **tmp;
862  tmp = dmatrix(1, NDIM, 1, NDIM);
863 #else
864  static int firstc = 1;
865  static FTYPE **tmp;
866  if (firstc) {
867  tmp = dmatrix(1, NDIM, 1, NDIM);
868  firstc = 0;
869  }
870 #endif
871 
872 
873  singfix=0;
874 
875 #if(1)
876  // check for coordinate singularity (using this avoids bad input to ludcmp())
877  metric_sing_check(whichcoord, generalmatrixlower, &anglesing, &centersing, &truedim);
878 #else
879  truedim=NDIM;
880 #endif
881 
882 
883 
884 
885  DLOOP(j,k) tmp[j + 1][k + 1] = generalmatrixlower[j][k];
886 
887 
888  if(ludcmp(tmp, truedim, indx - 1, &d)>=1){
889  if(debugfail>=2) dualfprintf(fail_file,"ludcmp failure in gdet_func_orig\n");
890 #if(1)
891  if(ISSPCMCOORD(whichcoord)){
892  // super hack
893  if(debugfail>=2) dualfprintf(fail_file,"Assuming on polar axis\n");
894  singfix=1;
895  }
896  else{
897  dualfprintf(fail_file,"ludcmp failure 2: whichcoord=%d\n",whichcoord);
898  myexit(3477);
899  }
900 #else
901  dualfprintf(fail_file,"ludcmp failure: Can't assume on polar axis for whichcoord=%d\n",whichcoord);
902  myexit(246);
903 #endif
904 
905  }
906 
907  if(singfix==0 || truedim<NDIM){
908  // below from 1..NDIM due to ludcmp requiring 1..N
909  for (j = 1; j <= NDIM; j++) d *= tmp[j][j];
910 
911  if(d>0.0 && d<SMALL){
912  finalvalue=0.0;
913  }
914  else if(d>0.0){
915  dualfprintf(fail_file,"Metric has bad signature: d=%21.15g\n",d);
916  DLOOP(j,k) dualfprintf(fail_file,"generalmatrixlower[%d][%d]=%21.15g\n",j,k,generalmatrixlower[j][k]); // generalmatrixlower[][] is 2d array since using NR routines.
917  myexit(3478);
918  }
919  else{
920  finalvalue=sqrt(-d);
921  }
922  }
923  else{
924  finalvalue=0.0;
925  }
926 
927 
928 #if(USEOPENMP)
929  // maintain thread safety
930  free_dmatrix(tmp, 1, NDIM, 1, NDIM);
931 #endif
932 
933 
934 
935  *gdet=finalvalue;
936 
937  if(singfix) return(-1); // indicates some problem, may want to report a bit
938  else return(0);
939 }
940 
941 
942 
943 
944 // t-r inverse only (for r=0 coordinate singularity)
945 // NDIM in size, but invert submatrix that's 2x2 in size
946 void matrix_inverse_2d(FTYPE (*genmatrixlower)[NDIM], FTYPE (*genmatrixupper)[NDIM])
947 {
948  int jj,kk;
949 
950  DLOOP(jj,kk) genmatrixupper[jj][kk] = 0.0;
951 
952 
953  genmatrixupper[TT][TT] = 1.0/(-genmatrixlower[TT][RR]*genmatrixlower[TT][RR]/genmatrixlower[RR][RR] + genmatrixlower[TT][TT]);
954 
955  genmatrixupper[TT][RR] = genmatrixupper[RR][TT] = -genmatrixlower[TT][RR]/(-genmatrixlower[TT][RR]*genmatrixlower[TT][RR] + genmatrixlower[RR][RR]*genmatrixlower[TT][TT]);
956 
957  genmatrixupper[RR][RR] = 1.0/(-genmatrixlower[TT][RR]*genmatrixlower[TT][RR]/genmatrixlower[TT][TT] + genmatrixlower[RR][RR]);
958 
959 }
960 
961 // t-r-\theta inverse only (for \theta={0,\pi} coordinate singularities)
962 // NDIM in size, but invert submatrix that's 3x3 in size
963 void matrix_inverse_3d(FTYPE (*genmatrixlower)[NDIM], FTYPE (*genmatrixupper)[NDIM])
964 {
965  int jj,kk;
966  FTYPE genmatrixlowerrhsq,genmatrixlowertrsq,genmatrixlowerthsq;
967 
968 
969  DLOOP(jj,kk) genmatrixupper[jj][kk] = 0.0;
970 
971 
972  genmatrixlowerrhsq = genmatrixlower[RR][TH]*genmatrixlower[RR][TH];
973  genmatrixlowertrsq = genmatrixlower[TT][RR]*genmatrixlower[TT][RR];
974  genmatrixlowerthsq = genmatrixlower[TT][TH]*genmatrixlower[TT][TH];
975 
976  genmatrixupper[TT][TT]=(genmatrixlowerrhsq - genmatrixlower[TH][TH]*genmatrixlower[RR][RR])/(genmatrixlower[RR][RR]*genmatrixlowerthsq - 2.0*genmatrixlower[RR][TH]*genmatrixlower[TT][TH]*genmatrixlower[TT][RR] + genmatrixlower[TH][TH]*genmatrixlowertrsq + genmatrixlowerrhsq*genmatrixlower[TT][TT] - genmatrixlower[TH][TH]*genmatrixlower[RR][RR]*genmatrixlower[TT][TT]);
977  genmatrixupper[TT][RR]=genmatrixupper[RR][TT]=(-(genmatrixlower[RR][TH]*genmatrixlower[TT][TH]) + genmatrixlower[TH][TH]*genmatrixlower[TT][RR])/(-2.0*genmatrixlower[RR][TH]*genmatrixlower[TT][TH]*genmatrixlower[TT][RR] + genmatrixlower[TH][TH]*genmatrixlowertrsq + genmatrixlowerrhsq*genmatrixlower[TT][TT] + genmatrixlower[RR][RR]*(genmatrixlowerthsq - genmatrixlower[TH][TH]*genmatrixlower[TT][TT]));
978  genmatrixupper[TT][TH]=genmatrixupper[TH][TT]=(genmatrixlower[RR][RR]*genmatrixlower[TT][TH] - genmatrixlower[RR][TH]*genmatrixlower[TT][RR])/(genmatrixlower[RR][RR]*genmatrixlowerthsq - 2.0*genmatrixlower[RR][TH]*genmatrixlower[TT][TH]*genmatrixlower[TT][RR] + genmatrixlower[TH][TH]*genmatrixlowertrsq + genmatrixlowerrhsq*genmatrixlower[TT][TT] - genmatrixlower[TH][TH]*genmatrixlower[RR][RR]*genmatrixlower[TT][TT]);
979  genmatrixupper[RR][RR]=(genmatrixlowerthsq - genmatrixlower[TH][TH]*genmatrixlower[TT][TT])/(genmatrixlower[RR][RR]*genmatrixlowerthsq - 2.0*genmatrixlower[RR][TH]*genmatrixlower[TT][TH]*genmatrixlower[TT][RR] + genmatrixlower[TH][TH]*genmatrixlowertrsq + genmatrixlowerrhsq*genmatrixlower[TT][TT] - genmatrixlower[TH][TH]*genmatrixlower[RR][RR]*genmatrixlower[TT][TT]);
980  genmatrixupper[RR][TH]=genmatrixupper[TH][RR]=(-(genmatrixlower[TT][TH]*genmatrixlower[TT][RR]) + genmatrixlower[RR][TH]*genmatrixlower[TT][TT])/(-2.0*genmatrixlower[RR][TH]*genmatrixlower[TT][TH]*genmatrixlower[TT][RR] + genmatrixlower[TH][TH]*genmatrixlowertrsq + genmatrixlowerrhsq*genmatrixlower[TT][TT] + genmatrixlower[RR][RR]*(genmatrixlowerthsq - genmatrixlower[TH][TH]*genmatrixlower[TT][TT]));
981  genmatrixupper[TH][TH]=(genmatrixlowertrsq - genmatrixlower[RR][RR]*genmatrixlower[TT][TT])/(genmatrixlower[RR][RR]*genmatrixlowerthsq - 2.0*genmatrixlower[RR][TH]*genmatrixlower[TT][TH]*genmatrixlower[TT][RR] + genmatrixlower[TH][TH]*genmatrixlowertrsq + genmatrixlowerrhsq*genmatrixlower[TT][TT] - genmatrixlower[TH][TH]*genmatrixlower[RR][RR]*genmatrixlower[TT][TT]);
982 
983 
984 }
985 
986 
987 // wrapper for symmetric matrix
988 void matrix_inverse_metric(int whichcoord, FTYPE *gcov, FTYPE *gcon)
989 {
990  FTYPE genmatrixlower[NDIM][NDIM];
991  FTYPE genmatrixupper[NDIM][NDIM];
992  int jj,kk;
993 
994  // copy both in case both needed for some reason
995  DLOOP(jj,kk) genmatrixlower[jj][kk]=gcov[GIND(jj,kk)];
996  DLOOP(jj,kk) genmatrixupper[jj][kk]=gcon[GIND(jj,kk)];
997 
998  matrix_inverse(whichcoord, genmatrixlower, genmatrixupper);
999 
1000  // translate back
1001  DLOOP(jj,kk) gcon[GIND(jj,kk)]=genmatrixupper[jj][kk];
1002 
1003 }
1004 
1005 /* invert genmatrixlower to get genmatrixupper */
1006 // can be used to invert any 2nd rank tensor (symmetric or not)
1007 // actually returns the inverse transpose, so if
1008 // genmatrixlower=T^j_k then out pops (iT)^k_j such that T^j_k (iT)^k_l = \delta^j_l
1009 void matrix_inverse(int whichcoord, FTYPE (*genmatrixlower)[NDIM], FTYPE (*genmatrixupper)[NDIM])
1010 {
1011  int pl,pliter;
1012  int j, k;
1013  int anglesing,centersing,truedim;
1014  void metric_sing_check(int whichcoord, FTYPE (*genmatrixlower)[NDIM], int *anglesing, int*centersing, int *truedim);
1015  void matrix_inverse_2d(FTYPE (*genmatrixlower)[NDIM], FTYPE (*genmatrixupper)[NDIM]);
1016  void matrix_inverse_3d(FTYPE (*genmatrixlower)[NDIM], FTYPE (*genmatrixupper)[NDIM]);
1017 
1018 
1019 #if(USEOPENMP)
1020  // maintain thread safety
1021  FTYPE **tmp;
1022  tmp = dmatrix(1, NDIM, 1, NDIM);
1023 #else
1024  static int firstc = 1;
1025  static FTYPE **tmp;
1026  if (firstc) {
1027  tmp = dmatrix(1, NDIM, 1, NDIM);
1028  firstc = 0;
1029  }
1030 #endif
1031 
1032 
1033 
1034 
1035  DLOOP(j,k) tmp[j + 1][k + 1] = genmatrixlower[j][k];
1036 
1037 #if(1)
1038  // check for singularities
1039  // only truedim used
1040  // allow avoiding of gaussj fail as detection of singularity
1041  metric_sing_check(whichcoord, genmatrixlower, &anglesing, &centersing, &truedim);
1042  // dualfprintf(fail_file,"anglesing=%d centersing=%d truedim=%d\n",anglesing,centersing,truedim);
1043 #else
1044  truedim=NDIM;
1045 #endif
1046 
1047  // 0-out all genmatrixupper
1048  DLOOP(j,k) genmatrixupper[j][k]=0.0;
1049 
1050 
1051  // DLOOP(j,k) dualfprintf(fail_file,"tmp[%d][%d]=%21.15g\n",j+1,k+1,tmp[j+1][k+1]);
1052 
1053  // GODMARK: Feeding in nan's results in gaussj segfaulting with erroneous access. Seems bad behavior!
1054  if(gaussj(tmp, truedim, NULL, 0)){
1055  // then singular
1056 #if(0) // new singularity check before if(gaussj) should work
1057  if(ISSPCMCOORD(whichcoord)){
1058  // super hack
1059  //if(centersing) matrix_inverse_2d(genmatrixlower,genmatrixupper);
1060  // else if(anglesing) matrix_inverse_3d(genmatrixlower,genmatrixupper);
1061  matrix_inverse_2d(genmatrixlower,genmatrixupper);
1062  }
1063  else{
1064  dualfprintf(fail_file,"whichcoord=%d\n",whichcoord);
1065  myexit(6243);
1066  }
1067 #else
1068  dualfprintf(fail_file,"Singularity check didn't work\n");
1069  dualfprintf(fail_file,"whichcoord=%d anglesing=%d centersing=%d truedim=%d\n",whichcoord,anglesing,centersing,truedim);
1070  DLOOP(j,k) dualfprintf(fail_file,"inputmatrix[%d][%d]=%21.15g\n",j,k,genmatrixlower[j][k]);
1071  myexit(2714);
1072 #endif
1073 
1074  }
1075  else{
1076  // assign but also transpose (shouldn't do in general, confusing)
1077  //DLOOP(j,k) genmatrixupper[j][k] = tmp[k + 1][j + 1];
1078  DLOOP(j,k) genmatrixupper[j][k] = tmp[j + 1][k + 1];
1079  }
1080 
1081 
1082 #if(1) // check for nan's
1083  DLOOP(j,k) if(!finite(genmatrixupper[j][k])){
1084  dualfprintf(fail_file,"Came out of matrix_inverse with inf/nan for genmatrixupper at j=%d k=%d\n",j,k);
1085  myexit(5);
1086  }
1087 #endif
1088 
1089 
1090 
1091 
1092 #if(USEOPENMP)
1093  // maintain thread safety
1094  free_dmatrix(tmp, 1, NDIM, 1, NDIM);
1095 #endif
1096 
1097 
1098 
1099 }
1100 
1101 /* invert genmatrixlower to get genmatrixupper */
1102 // can be used to invert any 2nd rank tensor (symmetric or not)
1103 // actually returns the inverse transpose, so if
1104 // genmatrixlower=T^j_k then out pops (iT)^k_j such that T^j_k (iT)^k_l = \delta^j_l
1105 void matrix_inverse_gen(int truedim, FTYPE (*genmatrixlower)[NDIM], FTYPE (*genmatrixupper)[NDIM])
1106 {
1107  int pl,pliter;
1108  int j, k;
1109 
1110 
1111 #if(USEOPENMP)
1112  // maintain thread safety
1113  FTYPE **tmp;
1114  tmp = dmatrix(1, NDIM, 1, NDIM);
1115 #else
1116  static int firstc = 1;
1117  static FTYPE **tmp;
1118  if (firstc) {
1119  tmp = dmatrix(1, NDIM, 1, NDIM);
1120  firstc = 0;
1121  }
1122 #endif
1123 
1124 
1125 
1126 
1127  DLOOP(j,k) tmp[j + 1][k + 1] = genmatrixlower[j][k];
1128 
1129 
1130  // 0-out all genmatrixupper
1131  DLOOP(j,k) genmatrixupper[j][k]=0.0;
1132 
1133 
1134  if(gaussj(tmp, truedim, NULL, 0)){
1135  // then singular
1136  dualfprintf(fail_file,"Singularity\n");
1137  DLOOP(j,k) dualfprintf(fail_file,"inputmatrix[%d][%d]=%21.15g\n",j,k,genmatrixlower[j][k]);
1138  myexit(2715);
1139  }
1140  else{
1141  // assign but also transpose (shouldn't do in general, confusing)
1142  //DLOOP(j,k) genmatrixupper[j][k] = tmp[k + 1][j + 1];
1143  DLOOP(j,k) genmatrixupper[j][k] = tmp[j + 1][k + 1];
1144  }
1145 
1146 
1147 #if(1) // check for nan's
1148  DLOOP(j,k) if(!finite(genmatrixupper[j][k])){
1149  dualfprintf(fail_file,"Came out of matrix_inverse_gen with inf/nan for genmatrixupper at j=%d k=%d\n",j,k);
1150  myexit(5);
1151  }
1152 #endif
1153 
1154 
1155 
1156 
1157 #if(USEOPENMP)
1158  // maintain thread safety
1159  free_dmatrix(tmp, 1, NDIM, 1, NDIM);
1160 #endif
1161 
1162 
1163 
1164 }
1165 
1166 
1167 /* invert genmatrixlower to get genmatrixupper */
1168 // can be used to invert any 2nd rank tensor (symmetric or not)
1169 // actually returns the inverse transpose, so if
1170 // genmatrixlower=T^j_k then out pops (iT)^k_j such that T^j_k (iT)^k_l = \delta^j_l
1171 void matrix_inverse_4d(FTYPE (*genmatrixlower)[NDIM], FTYPE (*genmatrixupper)[NDIM])
1172 {
1173  int pl,pliter;
1174  int j, k;
1175  void matrix_inverse_2d(FTYPE (*genmatrixlower)[NDIM], FTYPE (*genmatrixupper)[NDIM]);
1176  void matrix_inverse_3d(FTYPE (*genmatrixlower)[NDIM], FTYPE (*genmatrixupper)[NDIM]);
1177 
1178 
1179 #if(USEOPENMP)
1180  // maintain thread safety
1181  FTYPE **tmp;
1182  tmp = dmatrix(1, NDIM, 1, NDIM);
1183 #else
1184  static int firstc = 1;
1185  static FTYPE **tmp;
1186  if (firstc) {
1187  tmp = dmatrix(1, NDIM, 1, NDIM);
1188  firstc = 0;
1189  }
1190 #endif
1191 
1192 
1193 
1194 
1195  DLOOP(j,k) tmp[j + 1][k + 1] = genmatrixlower[j][k];
1196 
1197  int truedim=NDIM;
1198 
1199  // 0-out all genmatrixupper
1200  DLOOP(j,k) genmatrixupper[j][k]=0.0;
1201 
1202 
1203  // DLOOP(j,k) dualfprintf(fail_file,"tmp[%d][%d]=%21.15g\n",j+1,k+1,tmp[j+1][k+1]);
1204 
1205  // GODMARK: Feeding in nan's results in gaussj segfaulting with erroneous access. Seems bad behavior!
1206  if(gaussj(tmp, truedim, NULL, 0)){
1207  // then singular
1208  dualfprintf(fail_file,"Singular\n");
1209  DLOOP(j,k) dualfprintf(fail_file,"inputmatrix[%d][%d]=%21.15g\n",j,k,genmatrixlower[j][k]);
1210  myexit(2714);
1211  }
1212  else{
1213  // assign but also transpose (shouldn't do in general, confusing)
1214  //DLOOP(j,k) genmatrixupper[j][k] = tmp[k + 1][j + 1];
1215  DLOOP(j,k) genmatrixupper[j][k] = tmp[j + 1][k + 1];
1216  }
1217 
1218 
1219 #if(1) // check for nan's
1220  DLOOP(j,k) if(!finite(genmatrixupper[j][k])){
1221  dualfprintf(fail_file,"Came out of matrix_inverse_4d with inf/nan for genmatrixupper at j=%d k=%d\n",j,k);
1222  myexit(5);
1223  }
1224 #endif
1225 
1226 
1227 
1228 
1229 #if(USEOPENMP)
1230  // maintain thread safety
1231  free_dmatrix(tmp, 1, NDIM, 1, NDIM);
1232 #endif
1233 
1234 
1235 
1236 }
1237 
1238 
1239 
1240 
1241 
1242 void metric_sing_check(int whichcoord, FTYPE (*genmatrixlower)[NDIM], int *anglesing, int*centersing, int *truedim)
1243 {
1244 
1245 
1246  // check for singularity (always removes highest coordinates -- so can use same matrix inverse with lower dimension when on coordinate singularities)
1247  *anglesing=0;
1248  *centersing=0;
1249  if(ISSPCMCOORD(whichcoord)){
1250  if(fabs(genmatrixlower[PH][PH])<10.0*SMALL){
1251  *anglesing=1;
1252  }
1253  if(fabs(genmatrixlower[TH][TH])<10.0*SMALL){
1254  *centersing=1;
1255  }
1256  if(*centersing){
1257  *truedim=NDIM-2;
1258  }
1259  else if(*anglesing){
1260  *truedim=NDIM-1;
1261  }
1262  else *truedim=NDIM;
1263  }
1264  else *truedim=NDIM;
1265 
1266 
1267 }
1268 
1269 
1270 
1271 
1272 
1273 
1274 // compute the radius of the inner most stable circular orbit
1275 FTYPE rmso_calc(int which)
1276 {
1277  FTYPE rmso,Z1,Z2,sign ;
1278  FTYPE j;
1279 
1280  j=a/(MBH+SMALL); // so doesn't nan for MBH=0 since "a" should be "0" if MBH=0
1281 
1282  if(which==PROGRADERISCO) sign=1; else sign=-1;
1283 
1284  Z1 = 1. + pow(1. - j*j,1./3.)*(pow(1. + j,1./3.) +
1285  pow(1. - j, 1./3.)) ;
1286  Z2 = sqrt(fabs(3.*j*j + Z1*Z1)) ;
1287  rmso=3. + Z2-sign*sqrt(fabs(3. - Z1)*fabs(3. + Z1 + 2.*Z2)) ;
1288 
1289  return(MBH*rmso) ;
1290 }
1291 
1292 FTYPE uphi_isco_calc(int which,FTYPE rold)
1293 {
1294  FTYPE uphi;
1295  FTYPE sign;
1296  FTYPE Z1,Z2;
1297  FTYPE j,rnew;
1298 
1299  rnew=rold/MBH;
1300  j=a/MBH;
1301 
1302  if(which==PROGRADERISCO) sign=1; else sign=-1;
1303 
1304  Z1=rnew*rnew-sign*2.*j*sqrt(rnew)+j*j;
1305  Z2=rnew*(rnew*rnew-3.*rnew+sign*2.*j*sqrt(rnew));
1306 
1307  uphi=sign*Z1/sqrt(Z2);
1308 
1309  return(MBH*uphi);
1310 
1311 }
1312 
1313 FTYPE rhor_calc(int which)
1314 {
1315  FTYPE sign,rhor;
1316  FTYPE j,jsq,disc;
1317 
1318  j=a/(MBH+SMALL); // so doesn't nan for MBH=0 since "a" should be "0" if MBH=0
1319  jsq=j*j;
1320 
1321  if(which==0) sign=1; else sign=-1;
1322 
1323  disc=MAX(1.0 - jsq,0.0);
1324 
1325 
1326 
1327  rhor=1. +sign*sqrt(disc);
1328 
1329  return(rhor*MBH);
1330 }
1331 
1332 
1333 
1335 //
1336 // below are independent of user choice of metric/coords/grid
1337 //
1339 
1340 // presume get_geometry() only necessarily feeds back pointer where geometry is located
1341 // still required to have pointer point to physical allocated memory in general
1342 void get_and_copy_geometry(int i, int j, int k, int loc, struct of_geom *ptrgeom)
1343 {
1344  struct of_geom *ptrgeomorig;
1345 
1346  ptrgeomorig=ptrgeom;
1347 
1348  get_geometry(i,j,k,loc,ptrgeom); // potentially overwrites ptrgeom
1349 
1350  if(ptrgeom!=ptrgeomorig){
1351  // direct copy of geometry structure
1352  *ptrgeomorig=*ptrgeom;
1353  }
1354 
1355 }
1356 
1357 
1358 
1359 // set igdet part of geometry since more expensive and not always needed
1360 void set_igdet_old(struct of_geom *geom)
1361 {
1362  int pl,pliter;
1363 
1365  // avoids 0.0 for any sign of ptrgeom->e[pl]
1366 #if(GDETVOLDIFF==0)
1367 
1368  geom->igdetnosing = sign(geom->gdet)/(fabs(geom->gdet)+SMALL);
1369 
1370  // use PALLLOOP so compiler can optimize
1371 #if(WHICHEOM!=WITHGDET)
1372  PALLLOOP(pl) geom->IEOMFUNCNOSINGMAC(pl) = sign(geom->EOMFUNCMAC(pl))/(fabs(geom->EOMFUNCMAC(pl))+SMALL);
1373 #else
1374  // required to set to something since in general refer to this
1375  PALLLOOP(pl) geom->IEOMFUNCNOSINGMAC(pl)=geom->igdetnosing;
1376 #endif
1377 
1378 
1379 #else
1380 
1381  // volume regularization (correct to second order) // GODMARK: NOT FOR FINITE VOLUME WENO METHOD
1382  igdetnosing = sign(geom->gdetvol)/(fabs(geom->gdetvol)+SMALL);
1383  geom->igdetnosing = igdetnosing;
1384  // use PALLLOOP so compiler can optimize
1385  PALLLOOP(pl) geom->IEOMFUNCNOSINGMAC(pl) = igdetnosing;
1386 #endif
1387 
1388 
1389 }
1390 // set igdet part of geometry since more expensive and not always needed
1391 void set_igdetsimple_old(struct of_geom *geom)
1392 {
1393  int pl,pliter;
1394 
1395 
1396 #if(WHICHEOM!=WITHGDET)
1397  dualfprintf(fail_file,"Using set_igdetsimple() but WHICHEOM!=WITHGDET\n");
1398  myexit(342968347);
1399 #endif
1400 
1401 
1403  // avoids 0.0 for any sign of ptrgeom->e[pl]
1404 #if(GDETVOLDIFF==0)
1405  geom->igdetnosing = sign(geom->gdet)/(fabs(geom->gdet)+SMALL);
1406 #else
1407 
1408  // volume regularization (correct to second order) // GODMARK: NOT FOR FINITE VOLUME WENO METHOD
1409  igdetnosing = sign(geom->gdetvol)/(fabs(geom->gdetvol)+SMALL);
1410  geom->igdetnosing = igdetnosing;
1411 #endif
1412 
1413 }
1414 
1415 
1416 
1417 // obtain prime or non-prime alphalapse
1418 void alphalapse_func(struct of_geom *ptrgeom, int getprim, int whichcoord, FTYPE *X, FTYPE *gcov, FTYPE *gcon, FTYPE *alphalapse)
1419 {
1420 
1421  // set alpha -- fabs just for roundoff error
1422  // Note ptrgeom only has i,j,k,loc at this point
1423  *alphalapse = 1.0/sqrt(fabs(-gcon[GIND(TT,TT)]));
1424 
1425 }
1426 
1427 // obtain \beta^2
1428 void betasqoalphasq_func(struct of_geom *ptrgeom, int getprim, int whichcoord, FTYPE *X, FTYPE *gcov, FTYPE *gcon, FTYPE *betasqoalphasq)
1429 {
1430  int j;
1431 
1432  // \beta^i \beta_i / \alpha^2 = g^{ti} g_{ti}
1433  *betasqoalphasq = 0.0;
1434  SLOOPA(j) *betasqoalphasq += (gcov[GIND(TT,j)])*(gcon[GIND(TT,j)]);
1435 
1436 }
1437 
1438 // obtain \beta^i
1439 void beta_func(struct of_geom *ptrgeom, int getprim, int whichcoord, FTYPE *X, FTYPE *gcov, FTYPE *gcon, FTYPE alphalapse, FTYPE *beta)
1440 {
1441  int j;
1442  FTYPE alphasq;
1443 
1444  alphasq=alphalapse*alphalapse;
1445 
1446  // \beta^\mu = {0,\alpha^2 g^{ti}}
1447  beta[TT]=0.0;
1448  SLOOPA(j) beta[j] = gcon[GIND(TT,j)]*alphasq ;
1449 
1450 
1451 }
1452 
1453 
1454 
1455 
1456 
1457 // find the con/cov forms of the chosen metric
1458 void gset(int getprim, int whichcoord, int i, int j, int k, struct of_geom *ptrgeom)
1459 {
1460  // assumes loc=CENT
1461  int loc;
1462 
1463  loc=CENT;
1464  gset_genloc(getprim, whichcoord, i, j, k, loc, ptrgeom);
1465 
1466 }
1467 
1468 
1469 // find the con/cov forms of the chosen metric
1470 // fills in information like get_geometry but for arbitrary metric not just PRIMECOORDS
1471 // GODMARK: doesn't yet set igdet's
1472 // NOTE:
1473 // This function returns contents of ptrgeom when not necessarily internal coordinate system
1474 // So this function should be quite similar to set_grid.c for that part setting ptrgeom contents, but really setting contents rather than looking them up
1475 void gset_genloc(int getprim, int whichcoord, int i, int j, int k, int loc, struct of_geom *ptrgeom)
1476 {
1477  FTYPE X[NDIM];
1478 
1479  if(whichcoord>=0){
1480  coord_ijk(i, j, k, loc, X);
1481  }
1482  else if(whichcoord==PRIMECOORDS){ // special case
1483  // won't neeed X . Assumes user never requests PRIMECOORDS unless on grid. User should call getprim=1 with whichcoord>=0 for normal coordinates in prime coords
1484  }
1485  else{
1486  dualfprintf(fail_file,"gset(): no such whichcoord=%d\n",whichcoord);
1487  myexit(3466);
1488  }
1489 
1490  gset_X(getprim, whichcoord, i, j, k, loc, X, ptrgeom);
1491 
1492 }
1493 
1494 
1495 
1496 
1497 // find the con/cov forms of the chosen metric
1498 // fills in information like get_geometry but for arbitrary metric not just PRIMECOORDS
1499 // GODMARK: doesn't yet set igdet's
1500 // NOTE:
1501 // This function returns contents of ptrgeom when not necessarily internal coordinate system
1502 // So this function should be quite similar to set_grid.c for that part setting ptrgeom contents, but really setting contents rather than looking them up
1503 // i,j,k,loc can be filled or loc=NOWHERE then recomputes bl_coord()
1504 void gset_X(int getprim, int whichcoord, int i, int j, int k, int loc, FTYPE *X, struct of_geom *ptrgeom)
1505 {
1506  FTYPE V[NDIM];
1507  struct of_geom tempgeom;
1508  extern void assign_eomfunc(struct of_geom *geom, FTYPE *EOMFUNCNAME);
1509  void gcov_func(struct of_geom *ptrgeom, int getprim, int whichcoord, FTYPE *X, FTYPE *gcov, FTYPE *gcovpert);
1510  void gcon_func(struct of_geom *ptrgeom, int getprim, int whichcoord, FTYPE *X, FTYPE *gcov, FTYPE *gcon);
1511 
1512  FTYPE *gcovptr;
1513  FTYPE *gconptr;
1514  FTYPE *gcovpertptr;
1515 
1516 
1517  ptrgeom->i=i;
1518  ptrgeom->j=j;
1519  ptrgeom->k=k;
1520  ptrgeom->p=loc;
1521 
1522 
1523 #if(GETGEOMUSEPOINTER==0 || NEWMETRICSTORAGE==1)
1524  // then ptrgeom->gcov,gcon,gcovpert are real memory spaces
1525  gcovptr=ptrgeom->gcov;
1526  gconptr=ptrgeom->gcon;
1527  gcovpertptr=ptrgeom->gcovpert;
1528 #else
1529  // then need to use dummy pointer space that has real memory assigned
1530  gcovptr=ptrgeom->gengcov;
1531  gconptr=ptrgeom->gengcon;
1532  gcovpertptr=ptrgeom->gengcovpert;
1533 
1534  ptrgeom->gcov=ptrgeom->gengcov; // pointer
1535  ptrgeom->gcon=ptrgeom->gengcon; // pointer
1536  ptrgeom->gcovpert=ptrgeom->gengcovpert; // pointer
1537 #endif
1538 
1539 
1540  if(whichcoord>=0){
1541  if(X==NULL) coord_ijk(i, j, k, loc, X); // user passes X if X==NULL
1542  bl_coord_ijk_2(i, j, k, loc, X, V);
1543  gcov_func(ptrgeom,getprim,whichcoord,X,gcovptr,gcovpertptr);
1544  // must come after gcov_func() above
1545  if(gdet_func_metric(whichcoord,V,gcovptr,&(ptrgeom->gdet))!=0){
1546  if(debugfail>=2) dualfprintf(fail_file,"Caught gdet_func_metric() problem in gset_genloc()\n");
1547  }
1548  gcon_func(ptrgeom,getprim,whichcoord,X,gcovptr,gconptr); // must come after gcov_func() above
1549  alphalapse_func(ptrgeom,getprim,whichcoord,X,gcovptr,gconptr,&(ptrgeom->alphalapse));
1550  betasqoalphasq_func(ptrgeom,getprim,whichcoord,X,gcovptr,gconptr,&(ptrgeom->betasqoalphasq));
1551  beta_func(ptrgeom,getprim,whichcoord,X,gcovptr,gconptr,ptrgeom->alphalapse,ptrgeom->beta);
1552  eomfunc_func(ptrgeom, getprim, whichcoord,X,&(ptrgeom->EOMFUNCMAC(0)));
1553  assign_eomfunc(ptrgeom,&(ptrgeom->EOMFUNCMAC(0))); // must come after assigning ptrgeom->g above (won't use eomfuncgen if WHICHEOM==WITHGDET)
1554 #if(GDETVOLDIFF)
1555  // uses X,V (not det) from all locations
1556  gdetvol_func(ptrgeom,ptrgeom->gdet,&(ptrgeom->EOMFUNCMAC(0)),ptrgeom->gdetvol);
1557 #endif
1558  set_igdet_old(ptrgeom); // full 1/gdet and 1/EOMFUNCMAC(pl)
1559 
1560  }
1561  else if(whichcoord==PRIMECOORDS){ // special case
1562  get_and_copy_geometry(i,j,k,loc,ptrgeom);
1563  }
1564  else{
1565  dualfprintf(fail_file,"gset(): no such whichcoord=%d\n",whichcoord);
1566  myexit(3466);
1567  }
1568 
1569 }
1570 
1571 
1572 //void SHOULDNOTREACHHEREEVERBUGYOUHAVE(void)
1573 //{
1574 // exit(1);
1575 //}
1576 
1577 #define OPTMETRICLOOP 1 // whether to use highly optmized loop (assumes metric is symmetric)
1578 #define COMPUTEPERTURBEDMETRIC 0 // GODMARK: NOT CORRECT RIGHT NOW, so do NOT do it
1579 
1580 void gcov2gcovprim(struct of_geom *ptrgeom, FTYPE *X, FTYPE *V, FTYPE *gcov, FTYPE *gcovpert, FTYPE *gcovprim, FTYPE *gcovpertprim)
1581 {
1582  int j, k, l, m;
1583  FTYPE dxdxp[NDIM][NDIM];
1584  FTYPE tmpgcov[SYMMATRIXNDIM];
1585  FTYPE ftemp1,ftemp2;
1586  int q;
1587 
1588 
1589 
1590 
1591  // now take term by term:
1592  // g_{u v} = \vec{e_{\mu}}\cdot\vec{e_{\nu}}
1593  // * (dx/dx')_{mu} * (dx/dx')_{\nu} =
1594  // \vec{e'_{\mu}}\cdot\vec{e'_{\nu}}
1595 
1596  // dx/dx' where '=prim coords (i.e. nonuni coords)
1597  // dualfprintf(fail_file,"gcov2gcovprim: i=%d p=%d\n",ptrgeom->i,ptrgeom->p);
1598  dxdxprim_ijk_2(ptrgeom,X,V,dxdxp);
1599 
1600 #if(OPTMETRICLOOP==0)
1601  DLOOP(j,k){
1602  tmpgcov[GIND(j,k)] = 0.;
1603  for(l=0;l<NDIM;l++) for(m=0;m<NDIM;m++){
1604  // g_{mup nup} = g_{mu nu} T^mu_mup T^nu_nup
1605  // where T^mu_mup == dx^mu[BL]/dx^mup[KSP uni grid]
1606  tmpgcov[GIND(j,k)] += GINDASSIGNFACTOR(j,k)*gcov[GIND(l,m)] * dxdxp[l][j] * dxdxp[m][k];
1607  }
1608 
1609  }
1610  DLOOP(j,k){
1611  // also must be outside above DLOOP because tempgcov can be reset to zero if GINDASSIGNFACTOR is zero, and then gcovprim would be zero, but that normally points also to gcov!
1612  // use tmpgcov since gcon might be same address as gcovprim
1613  gcovprim[GIND(j,k)] = tmpgcov[GIND(j,k)];
1614  }
1615 #else
1616  transgcov(gcov,dxdxp,gcovprim);
1617 #endif
1618 
1619  // DLOOP(j,k) dualfprintf(fail_file,"prim gcov[%d][%d]=%21.15g\n",j,k,gcov[GIND(j,k)]);
1620  // DLOOP(j,k) dualfprintf(fail_file,"prim gcovprim[%d][%d]=%21.15g\n",j,k,gcovprim[GIND(j,k)]);
1621  // DLOOP(j,k) dualfprintf(fail_file,"prim dxdxp[%d][%d]=%21.15g\n",j,k,dxdxp[j][k]);
1622 
1623  get_gcovpert(gcovprim,gcovpert,gcovpertprim);
1624 
1625 
1626 }
1627 
1628 
1629 
1630 // get perturbed part of gcov
1631 void get_gcovpert(FTYPE *gcovprim, FTYPE *gcovpert, FTYPE *gcovpertprim)
1632 {
1633  int j, k, l, m;
1634  FTYPE ftemp1,ftemp2;
1635  int q;
1636 
1638  //
1639  // perturbed terms
1640  //
1642 
1643 #if(COMPUTEPERTURBEDMETRIC)
1644  // SUPERGODMARK GODMARK: This only works for non-rel gravity if dxdxp is diagonal. Not sure what to do in general
1645  DLOOPA(q){
1646  // dualfprintf(fail_file,"gcovpert[%d]=%21.15g\n",q,gcovpert[q]);
1647 
1648 
1649  // get q-q term
1650  // -1 + g_{q'q'} for q!=TT, else SOMECONSTANT + g_{q'q'} for q=TT
1651 
1652  // for spatial parts, only care about deviations from constant in getting connection coefficients
1653  // avoids catastrophic cancellation (GODMARK)
1654  if( (q!=TT)&&(defcoord==UNIFORMCOORDS)) ftemp1=0.0;
1655  // then unsure how to handle in general, so just leave alone (SUPERGODMARK GODMARK)
1656  // problem is with machine error in dxdxp leading to apparently large connection coefficients due to that machine error
1657  // also don't know ahead of time what constant to choose for prim coordinate quantities
1658  // and don't know how to handle variations in dxdxp, since then no constant to subtract off, but grid accelerations apparently wash out non-rel gravity accelerations?
1659  // except for time term, which is solid
1660  else ftemp1=(-1.0 + dxdxp[q][q] * dxdxp[q][q]);
1661  if(q==TT) ftemp1 *=-1.0;
1662  gcovpertprim[q] = (gcovpert[q] * dxdxp[q][q] * dxdxp[q][q]) + ftemp1;
1663 
1664  // now add 15 other terms
1665  ftemp2 = 0.;
1666  for(l=0;l<NDIM;l++) for(m=0;m<NDIM;m++){
1667  if((l!=q)&&(m!=q)) ftemp2+= gcov[GIND(l,m)] * dxdxp[l][q] * dxdxp[m][q];
1668  }
1669  // add other 15 terms to answer for total of 16 terms
1670  gcovpertprim[q]+=ftemp2;
1671 
1672 
1673  // dualfprintf(fail_file,"dxdxp[%d][%d]=%21.15g\n",q,q,dxdxp[q][q]);
1674  // dualfprintf(fail_file,"ftemp1[%d]=%21.15g ftemp2[%d]=%21.15g gcovpertprim[%d]=%21.15g\n",q,ftemp1,q,ftemp2,q,gcovpertprim[q]);
1675  }
1676 #elif(0)
1677  // override for now
1678  gcovpertprim[TT]=gcovprim[GIND(TT,TT)]-mink(TT,TT);
1679  SLOOPA(q){
1680  gcovpertprim[q]=gcovprim[GIND(q,q)]-mink(q,q);
1681  }
1682 #elif(1)
1683  // override for now
1684  gcovpertprim[TT]=gcovprim[GIND(TT,TT)]+1.0;
1685  gcovpertprim[RR]=gcovprim[GIND(RR,RR)]-1.0;
1686  gcovpertprim[TH]=gcovprim[GIND(TH,TH)]-1.0;
1687  gcovpertprim[PH]=gcovprim[GIND(PH,PH)]-1.0;
1688 #endif
1689 
1690 }
1691 
1692 
1693 void transgcov_old(FTYPE *gcov, FTYPE (*dxdxp)[NDIM], FTYPE *gcovprim)
1694 {
1695  int j, k, l, m;
1696  FTYPE tmpgcov[SYMMATRIXNDIM];
1697 
1698  DLOOP(j,k){ // OPTMARK: In places where deal with symmetric metric that's using GIND(), could introduce new DLOOPMET(j,k) that only goes over required elements. Only works for assignment to LHS, not RHS since need factors of two that arrive naturally in sum on RHS.
1699  tmpgcov[GIND(j,k)] = 0.;
1700  for(l=0;l<NDIM;l++) for(m=0;m<NDIM;m++){
1701  // g_{mup nup} = g_{mu nu} T^mu_mup T^nu_nup
1702  // where T^mu_mup == dx^mu[BL]/dx^mup[KSP uni grid]
1703  tmpgcov[GIND(j,k)] += gcov[GIND(l,m)] * dxdxp[l][j] * dxdxp[m][k]; // GINDASSIGNFACTOR(j,k) not needed because tmpgcov not += to itself. RHS is summed over as if entire metric there, as wanted.
1704  }
1705  }
1706  DLOOP(j,k){
1707  // use tmpgcov since gcov might be same address as gcovprim
1708  gcovprim[GIND(j,k)] = tmpgcov[GIND(j,k)];
1709  }
1710 
1711 }
1712 
1713 
1714 
1715 // used to transform gcov and put back into gcov
1716 void transgcovself(FTYPE *gcov, FTYPE (*trans)[NDIM])
1717 {
1718  FTYPE gcovprim[SYMMATRIXNDIM];
1719  int j,k;
1720 
1721  transgcov(gcov,trans,gcovprim);
1722  DLOOP(j,k) gcov[GIND(j,k)] = gcovprim[GIND(j,k)];
1723 
1724 }
1725 
1726 
1727 // used to transform gcov&gcovprim and put back into gcov&gcovprim
1728 void transgcovgcovpertself(FTYPE *gcov, FTYPE *gcovpert, FTYPE (*trans)[NDIM])
1729 {
1730  FTYPE gcovprim[SYMMATRIXNDIM];
1731  FTYPE gcovpertprim[NDIM];
1732  int j,k;
1733 
1734  transgcov(gcov,trans,gcovprim);
1735  DLOOP(j,k) gcov[GIND(j,k)] = gcovprim[GIND(j,k)];
1736 
1737  get_gcovpert(gcovprim, gcovpert, gcovpertprim);
1738  DLOOPA(j) gcovpert[j] = gcovpertprim[j];
1739 
1740 
1741 }
1742 
1743 
1744 
1745 
1746 // gcov might be same memory address as gcovprim, so use tmp
1747 // assumes metric is symmetric 2nd rank tensor
1748 void transgcov(FTYPE *gcov, FTYPE (*trans)[NDIM], FTYPE *gcovprim)
1749 {
1750  int j, k, l, m;
1751  FTYPE tmp[NDIM][NDIM];
1752 
1753  // g_{\alpha \beta} = g_{\mu \nu} \Lambda^\mu_\alpha \Lambda^\nu_\beta
1754 
1755 
1756  /*
1757  // 4 along diagonal and 6 off-diagonal with 6 other identical values
1758  #define GCOV_DOT_TRANS_DOT_TRANS(a,b)\
1759  gcov[GIND(0,0)] * trans[0][a]* trans[0][b]\
1760  + gcov[GIND(1,1)] * trans[1][a]* trans[1][b]\
1761  + gcov[GIND(2,2)] * trans[2][a]* trans[2][b]\
1762  + gcov[GIND(3,3)] * trans[3][a]* trans[3][b]\
1763  + 2.0*gcov[GIND(0,1)] * trans[0][a]* trans[1][b]\
1764  + 2.0*gcov[GIND(0,2)] * trans[0][a]* trans[2][b]\
1765  + 2.0*gcov[GIND(0,3)] * trans[0][a]* trans[3][b]\
1766  + 2.0*gcov[GIND(1,2)] * trans[1][a]* trans[2][b]\
1767  + 2.0*gcov[GIND(1,3)] * trans[1][a]* trans[3][b]\
1768  + 2.0*gcov[GIND(2,3)] * trans[2][a]* trans[3][b]
1769  */
1770 
1771 
1772  // 4 along diagonal and 6 off-diagonal with 6 other identical values
1773 #define GCOV_DOT_TRANS_DOT_TRANS(a,b) \
1774  gcov[GIND(0,0)] * trans[0][a]* trans[0][b] \
1775  + gcov[GIND(1,1)] * trans[1][a]* trans[1][b] \
1776  + gcov[GIND(2,2)] * trans[2][a]* trans[2][b] \
1777  + gcov[GIND(3,3)] * trans[3][a]* trans[3][b] \
1778  + gcov[GIND(0,1)] * (trans[0][a]* trans[1][b] + trans[1][a]* trans[0][b]) \
1779  + gcov[GIND(0,2)] * (trans[0][a]* trans[2][b] + trans[2][a]* trans[0][b]) \
1780  + gcov[GIND(0,3)] * (trans[0][a]* trans[3][b] + trans[3][a]* trans[0][b]) \
1781  + gcov[GIND(1,2)] * (trans[1][a]* trans[2][b] + trans[2][a]* trans[1][b]) \
1782  + gcov[GIND(1,3)] * (trans[1][a]* trans[3][b] + trans[3][a]* trans[1][b]) \
1783  + gcov[GIND(2,3)] * (trans[2][a]* trans[3][b] + trans[3][a]* trans[2][b])
1784 
1785 
1786  // first do 4 along diagonal
1787  DLOOPA(j) tmp[j][j] = GCOV_DOT_TRANS_DOT_TRANS(j,j);
1788 
1789  // now do 6 others off-diagonal
1790  j=0; k=1; tmp[j][k] = GCOV_DOT_TRANS_DOT_TRANS(j,k);
1791  j=0; k=2; tmp[j][k] = GCOV_DOT_TRANS_DOT_TRANS(j,k);
1792  j=0; k=3; tmp[j][k] = GCOV_DOT_TRANS_DOT_TRANS(j,k);
1793  j=1; k=2; tmp[j][k] = GCOV_DOT_TRANS_DOT_TRANS(j,k);
1794  j=1; k=3; tmp[j][k] = GCOV_DOT_TRANS_DOT_TRANS(j,k);
1795  j=2; k=3; tmp[j][k] = GCOV_DOT_TRANS_DOT_TRANS(j,k);
1796 
1797 
1798  // copy over result assuming tmp on upper-diagonal only, but filling gcovprim fully
1799  // 4 along diagonal (00,11,22,33) and 6 off-diagonal (01, 02, 03, 12, 13, 23), 6 more as copies of off-diagonal
1800  DLOOPA(j) gcovprim[GIND(j,j)] = tmp[j][j];
1801  gcovprim[GIND(0,1)]=gcovprim[GIND(1,0)]=tmp[0][1];
1802  gcovprim[GIND(0,2)]=gcovprim[GIND(2,0)]=tmp[0][2];
1803  gcovprim[GIND(0,3)]=gcovprim[GIND(3,0)]=tmp[0][3];
1804  gcovprim[GIND(1,2)]=gcovprim[GIND(2,1)]=tmp[1][2];
1805  gcovprim[GIND(1,3)]=gcovprim[GIND(3,1)]=tmp[1][3];
1806  gcovprim[GIND(2,3)]=gcovprim[GIND(3,2)]=tmp[2][3];
1807 
1808 }
1809 
1810 
1811 void gcon2gconprim(struct of_geom *ptrgeom, FTYPE *X, FTYPE *V, FTYPE *gcon,FTYPE *gconprim)
1812 {
1813  int j, k, l, m;
1814  //FTYPE dxdxp[NDIM][NDIM],idxdxp[NDIM][NDIM];
1815  FTYPE idxdxp[NDIM][NDIM];
1816  FTYPE tmpgcon[SYMMATRIXNDIM];
1817 
1818  // see transforms.c and mettometp() and see gcov2gcovprim()
1819  idxdxprim_ijk_2(ptrgeom, X, V, idxdxp);
1820 
1821  // dualfprintf(fail_file,"mi in gcon2gconprim\n");
1822  // PRIMECOORDS indicates no special considerations if fails to get inverse
1823  // matrix_inverse(PRIMECOORDS, dxdxp,idxdxp);
1824 
1825  DLOOP(j,k) tmpgcon[GIND(j,k)] = 0.;
1826  DLOOP(j,k){
1827  for(l=0;l<NDIM;l++){
1828  for(m=0;m<NDIM;m++){
1829  tmpgcon[GIND(j,k)] += GINDASSIGNFACTOR(j,k)*idxdxp[j][l] * idxdxp[k][m] * gcon[GIND(l,m)] ;
1830  }
1831  }
1832  }
1833  // use tmpgcon since gcon might be same address as gconprim
1834  DLOOP(j,k) gconprim[GIND(j,k)] = tmpgcon[GIND(j,k)];
1835 
1836 }
1837 
1838 
1839 
1840 
1841 
1842 
1843 
1844 
1845 void setup_delta(int whichfun,int whichdifftype, FTYPE defaultdelta, struct of_geom *geom, struct of_geom (*localptrgeoml)[NDIM], struct of_geom (*localptrgeomh)[NDIM], FTYPE *truedelta)
1846 {
1847  int jj;
1848  int localwhichdifftype;
1849  int N[NDIM];
1850 
1851  N[0]=2; // if evolving metric then 2, but doesn't matter since at CENT always
1852  N[1]=N1;
1853  N[2]=N2;
1854  N[3]=N3;
1855 
1856 
1857 
1858  if(whichfun==0 || whichfun==1){ // connection or dxdxp
1859  if(whichdifftype==DIFFGAMMIE || whichdifftype==DIFFNUMREC){
1860  localwhichdifftype=0; // infinitesimal
1861  }
1862  else if(whichdifftype==DIFFFINITE){
1863  localwhichdifftype=1; // finite
1864  }
1865  }
1866  else{
1867  dualfprintf(fail_file,"no such whichfun=%d for whichdifftype=%d",whichfun,whichdifftype);
1868  myexit(915);
1869  }
1870 
1871 
1872 
1873  if(localwhichdifftype==0){ // infinitesimal
1874  DLOOPA(jj){
1875  // specify NOWHERE so that won't use gridded position
1876  ((*localptrgeoml)[jj]).p=((*localptrgeomh)[jj]).p=NOWHERE;
1877 
1878  ((*localptrgeoml)[jj]).i=((*localptrgeomh)[jj]).i=(geom->i);
1879  ((*localptrgeoml)[jj]).j=((*localptrgeomh)[jj]).j=(geom->j);
1880  ((*localptrgeoml)[jj]).k=((*localptrgeomh)[jj]).k=(geom->k);
1881 
1882  // an infinitesimal fraction of total grid difference (so scales with artificial startx and endx)
1883  truedelta[jj]=defaultdelta*Diffx[jj];
1884  }
1885 #if(DOEVOLVEMETRIC)
1886  // if evolving metric, then use CENT and gcov_func() will use X[0] to choose if present time or not
1887  // GODMARK: assume didstoremetricdata==1 is set when presenttime==1 is set so doesn't reach bl_coord_ijk_2() and do coord(i,j,k...) that would be wrong
1888  jj=TT;
1889  ((*localptrgeoml)[jj]).p=((*localptrgeomh)[jj]).p=CENT; // when taking temporal differences, values are always spatially located at loc=CENT
1890 #endif
1891 
1892  }
1893  else if(localwhichdifftype==1){ // finite
1894  // use gridded data
1895  // assumes connection at CENT
1896  if(geom->p==CENT){
1897 
1898 
1899  DLOOPA(jj){
1900  ((*localptrgeoml)[jj]).p=((*localptrgeomh)[jj]).p=CENT + jj*(N[jj]!=1); // jj=0 -> CENT , jj=1 -> FACE1 , jj=2 -> FACE2 , jj=3 -> FACE3
1901 
1902  ((*localptrgeoml)[jj]).i=(geom->i);
1903  ((*localptrgeoml)[jj]).j=(geom->j);
1904  ((*localptrgeoml)[jj]).k=(geom->k);
1905 
1906  // GODMARK: Note that infinitesimal version obtains correct connection even in reduced dimensions, while the finite version just reduces to no connection if reduced dimension (so this assumes something about the position or may not even be right in general)
1907  ((*localptrgeomh)[jj]).i=(geom->i) + (jj==1 ? SHIFT1 : 0);
1908  ((*localptrgeomh)[jj]).j=(geom->j) + (jj==2 ? SHIFT2 : 0);
1909  ((*localptrgeomh)[jj]).k=(geom->k) + (jj==3 ? SHIFT3 : 0);
1910 
1911  truedelta[jj]=dx[jj];
1912  }
1913  }
1914  else{
1915  dualfprintf(fail_file,"No localptrgeom defined for p=%d\n",geom->p);
1916  myexit(916);
1917  }
1918  }
1919  else{
1920  dualfprintf(fail_file,"No such difftype=%d\n",localwhichdifftype);
1921  myexit(917);
1922  }
1923 
1924 #if(0)
1925  dualfprintf(fail_file,"whichfun=%d whichdifftype=%d localwhichdifftype=%d\n",whichfun,whichdifftype,localwhichdifftype);
1926  DLOOPA(jj) dualfprintf(fail_file,"%d :: truedelta=%21.15g\n",jj,truedelta[jj]);
1927  DLOOPA(jj) dualfprintf(fail_file,"%d :: low: i=%d j=%d k=%d p=%d :: high: i=%d j=%d k=%d p=%d\n",jj,((*localptrgeoml)[jj]).i,((*localptrgeoml)[jj]).j,((*localptrgeoml)[jj]).k,((*localptrgeoml)[jj]).p,((*localptrgeomh)[jj]).i,((*localptrgeomh)[jj]).j,((*localptrgeomh)[jj]).k,((*localptrgeomh)[jj]).p);
1928 #endif
1929 
1930 }
1931 
1932 
1933 
1934 
1935 
1936 
1937 
1938 
1939 
1940 
1941 
1942 
1943 /*
1944  this gives the connection coefficient \Gamma^{i}_{j,k} =
1945  GLOBALMETMACP1A0(conn,..,i,j,k) where i,j,k = {0,1,2,3} corresponds to {t,r,theta,phi}
1946 */
1947 
1948 /*
1949  we also compute the 2nd connection:
1950  -d/dj(ln(\detg))
1951 */
1952 
1953 
1954 
1955 
1956 /*
1957 
1958  Based upon EOM:
1959 
1960  $$
1961  (f_{(\nu)} T^t_\nu)_{,t}
1962  = -(f_{(\nu)}T^j_\nu)_{,j}
1963  + f_{(\nu)} [
1964  T^\lambda_\nu[\ln(f_{(\nu)}/\detg)]_{,\lambda}
1965  + T^\mu_\lambda \Gamma^\lambda_{\mu\nu}
1966  + \ln(f_{(\nu)})_{,t} T^t_\nu
1967  ]
1968  $$
1969 
1970  More compactly (JCM: 05/07/08):
1971 
1972  $$
1973  d_\mu (f T^\mu_\nu) = f[ T^\lambda_\kappa \Gamma^\kappa_{\nu\lambda} - d_\mu ln (\detg/f) T^\mu_\nu]
1974  $$
1975 
1976  SUPERGODMARK: for if(WHICHEOM!=WITHGDET) don't yet compute correct conn2 I believe
1977 
1978  // to avoid body forces in general, must compute (e.g. through correction):
1979 
1980  $$
1981  \Gamma^\kappa_{\nu\lambda}[new] = \Gamma^\kappa_{\nu\lambda}[old] - (1/4)[\Gamma^\alpha_{\nu\alpha} - ( (d_\nu f)/f + d_\nu \ln(\detg/f) )]\delta^\kappa_\lambda
1982  $$
1983 
1984  or:
1985 
1986  $$
1987  \Gamma^\kappa_{\nu\lambda} += - (1/4)[ \Gamma^\alpha_{\nu\alpha} - ( (d_\nu f)/f - conn2_\nu )]\delta^\kappa_\lambda
1988  $$
1989 
1990  For the WHICHEOM!=WITHGDET version, in general this would require a different \Gamma for each EOM. For simplicity just assume f=\detg and don't allow this non-body version unless WHICHEOM==WITHGDET, and so then one has:
1991 
1992  $$
1993  \Gamma^\kappa_{\nu\lambda} += - (1/4)[ \Gamma^\alpha_{\nu\alpha} - (d_\nu \detg)/\detg ]\delta^\kappa_\lambda
1994  $$
1995 
1996  The above is only true for second order scheme. For higher order FLUXRECON scheme one has:
1997 
1998  $$
1999  \Gamma^\kappa_{\nu\lambda} += - (1/4)[ \Gamma^\alpha_{\nu\alpha} - (d_\nu a2c_\nu \detg)/\detg ]\delta^\kappa_\lambda
2000  $$
2001 
2002  Problems:
2003  1) But a2c_\nu would normally be chosen adaptively and so one wouldn't have good cancellation.
2004  2) Fact that \detg is there means difficult to generally have cancellation otherwise unless use NOGDET
2005  3) NOGDET not yet setup for EVOLVEMETRIC, but otherwise for spherical polar coordinates should just use NOGDET for r and \theta.
2006  4) Alternatively, can perform correction every timestep and use same a2c as used for each flux -- expensive
2007  Problem is b^2/2 and p will be treated differently for SPLITMAEM unless constant all weights
2008  5) So seems only solution apart from NOGDET is to use constant all weights that is unstable
2009 
2010 */
2011 
2012 
2013 
2014 // if evolving in time, then metric changed from Xmetricold[TT] to now (t)
2015 // by shifting time in past we tell metric to use old metric
2016 // if t==Xmetricold[0], then feeding the below tells gcov to assume standard difference at present time
2017 #define DELTAFORTIMEh(DELTA) (Xmetricnew[TT]!=Xmetricold[TT] ? (0.0) : 2.0*DELTA)
2018 #define DELTAFORTIMEl(DELTA) (Xmetricnew[TT]!=Xmetricold[TT] ? (-(Xmetricnew[TT]-Xmetricold[TT])) : DELTA)
2019 
2020 #define MYDELTAh(DELTA,k) ( k==TT ? DELTAFORTIMEh(DELTA) : +DELTA*0.5 )
2021 #define MYDELTAl(DELTA,k) ( k==TT ? DELTAFORTIMEl(DELTA) : -DELTA*0.5 )
2022 
2023 
2024 /* NOTE: parameter hides global variable */
2025 // note that inputted geom is used not only for i,j,k but for actual CENT gcon, etc.
2026 void conn_func_numerical1(FTYPE DELTA, FTYPE *X, struct of_geom *geom,
2027  FTYPE (*conn)[NDIM][NDIM],FTYPE *conn2)
2028 {
2029  int i, j, k, l;
2030  int kk,jj;
2031  FTYPE conndiag[NDIM],conndiag2[NDIM];
2032  FTYPE gdethgen[NDIM],gdetlgen[NDIM];
2033  FTYPE lngdethgen[NDIM],lngdetlgen[NDIM];
2034  // FTYPE *gdeth,*gdetl;
2035  // FTYPE *lngdeth,*lngdetl;
2036  FTYPE gdetmid;
2037  FTYPE tmp[NDIM][NDIM][NDIM];
2038  FTYPE Xhgen[NDIM][NDIM];
2039  FTYPE Xlgen[NDIM][NDIM];
2040  FTYPE signdXgen[NDIM];
2041  // FTYPE *Xh, *Xl;
2042  // FTYPE signdX;
2043  FTYPE V[NDIM];
2044  FTYPE gmid[SYMMATRIXNDIM];
2045  FTYPE ghgen[NDIM][SYMMATRIXNDIM];
2046  FTYPE glgen[NDIM][SYMMATRIXNDIM];
2047  // FTYPE (*gh)[NDIM];
2048  // FTYPE (*gl)[NDIM];
2049  FTYPE gcovpertmid[NDIM];
2050  FTYPE gcovperthgen[NDIM][NDIM],gcovpertlgen[NDIM][NDIM];
2051  // FTYPE *gcovperth,*gcovpertl;
2052  void gcov_func(struct of_geom *ptrgeom, int getprim, int whichcoord, FTYPE *X, FTYPE *gcov, FTYPE *gcovpert);
2053  int dfridr(FTYPE (*func)(struct of_geom *,FTYPE*,int,int), struct of_geom *geom, FTYPE *X,int ii, int jj, int kk, FTYPE *ans);
2054  FTYPE gcov_func_mcoord(struct of_geom *ptrgeom, FTYPE* X, int i, int j);
2055  FTYPE lngdet_func_mcoord(struct of_geom *ptrgeom, FTYPE* X, int i, int j);
2056  //
2057  // below variables used for setup_delta()
2058  // and used for defining whether infinitesimal or finite difference
2059  FTYPE truedelta[NDIM];
2060  struct of_geom localgeoml[NDIM];
2061  struct of_geom localgeomh[NDIM];
2062  struct of_geom (*localptrgeoml)[NDIM];
2063  struct of_geom (*localptrgeomh)[NDIM];
2064  struct of_geom localgeom;
2065  struct of_geom *localptrgeom=&localgeom;
2066  void setup_delta(int whichfun,int whichdifftype, FTYPE defaultdelta, struct of_geom *geom, struct of_geom (*localptrgeoml)[NDIM], struct of_geom (*localptrgeomh)[NDIM], FTYPE *truedelta);
2067  FTYPE gdet_func_mcoord_usegcov(FTYPE *gcovmcoord, struct of_geom *ptrgeom, FTYPE* X, int i, int j);
2068  int doingmachinebody;
2069  int setup_XlXh(FTYPE *X,FTYPE *truedelta, FTYPE (*Xlgen)[NDIM],FTYPE (*Xhgen)[NDIM],FTYPE *signdXgen);
2070  int compute_metricquantities_midlh(int donormal, int domachinebody
2071  ,struct of_geom *geom, FTYPE *X, FTYPE *gmid, FTYPE *gcovpertmid,FTYPE *gdetmid,FTYPE *gdetlgen,FTYPE *gdethgen
2072  ,struct of_geom (*localptrgeoml)[NDIM],struct of_geom (*localptrgeomh)[NDIM],FTYPE (*Xlgen)[NDIM],FTYPE (*Xhgen)[NDIM]
2073  ,FTYPE *lngdetlgen, FTYPE *lngdethgen, FTYPE (*glgen)[SYMMATRIXNDIM], FTYPE (*ghgen)[SYMMATRIXNDIM], FTYPE (*gcovpertlgen)[NDIM], FTYPE (*gcovperthgen)[NDIM]
2074  );
2075  int failreturn;
2076  int conndertypelocal;
2077  FTYPE value;
2078 
2079 
2080 
2081 
2082  // setup conditional
2083  doingmachinebody = CONNMACHINEBODY && WHICHEOM==WITHGDET;
2084 
2085 
2086 
2087 
2088 
2089  localptrgeoml=&localgeoml;
2090  localptrgeomh=&localgeomh;
2091 
2092 
2093 
2094 
2095 #if(DOEVOLVEMETRIC==1 && CONNDERTYPE==DIFFNUMREC)
2096  // no choice in sense that NUMREC is too slow
2097  dualfprintf(fail_file,"Not good idea to use DOEVOLVEMETRIC==1 with CONNDERTYPE==DIFFNUMREC\n");
2098  myexit(7225);
2099 #endif
2100 
2101 
2102 
2104  //
2105  // gabc_{ijk}=dg_{ij}/dx^k
2106  // gammie derivative (attempt to get analytical value) or finite difference (true finite difference with given cell size)
2107  //
2109  if(CONNDERTYPE==DIFFGAMMIE || CONNDERTYPE==DIFFFINITE || CONNDERTYPE==DIFFNUMREC){ // get default value even for DIFFNUMREC
2110 
2111  if(CONNDERTYPE==DIFFNUMREC) conndertypelocal=DIFFGAMMIE; // assume ok to use gammie diff as default
2112  else conndertypelocal=CONNDERTYPE; // normal
2113 
2114 
2116  //
2117  // Setup delta
2118  //
2120  // 0 indicates connection type
2121  setup_delta(0,conndertypelocal,DELTA,geom,localptrgeoml,localptrgeomh,truedelta);
2122 
2123 
2124  // setup Xl and Xh and signdX
2125  setup_XlXh(X,truedelta,Xlgen,Xhgen,signdXgen);
2126 
2127  // get metric quantities
2128  compute_metricquantities_midlh(1,doingmachinebody
2129  ,geom, X, gmid, gcovpertmid,&gdetmid,gdetlgen,gdethgen
2130  ,localptrgeoml,localptrgeomh,Xlgen,Xhgen
2131  ,lngdetlgen, lngdethgen, glgen, ghgen, gcovpertlgen, gcovperthgen
2132  );
2133 
2134  // resetup Xl and Xh and signdX since above "compute" function will have overwritten X positions for non-existent dimensions
2135  // if do this, must set truedelta correctly consistent in setup_delta()
2136  setup_XlXh(X,truedelta,Xlgen,Xhgen,signdXgen);
2137 
2138 
2140  //
2141  // Now compute derivatives
2142  //
2144 
2145  for (k = 0; k < NDIM; k++) {
2146 
2147 
2148  if(WHICHEOM!=WITHGDET){
2149  //$$
2150  //d_\mu (f T^\mu_\nu) = f[ T^\lambda_\kappa \Gamma^\kappa_{\nu\lambda} - d_\mu ln (\detg/f) T^\mu_\nu]
2151  //$$
2152  conn2[k]= signdXgen[k]*(lngdethgen[k] - lngdetlgen[k]) / (Xhgen[k][k] - Xlgen[k][k]);
2153  }
2154  else{
2155  conn2[k]=0.0; // no 2nd connection then
2156  }
2157 
2158 
2159  // answer is symmetric on i,j since uses g_{ij}, so only do part of work
2160  for (i = 0; i < NDIM; i++){
2161  for (j = 0; j <=i; j++){
2162  // d(1+g_{tt}) -> dg_{tt}, so can use 1+g_{tt} for accurate non-relativistic gravity
2163  //if(i==j) conn[i][j][k] = (gcovperth[i] - gcovpertl[i]) / (Xh[k] - Xl[k]);
2164  // else
2165  conn[i][j][k] = signdXgen[k]*(ghgen[k][GIND(i,j)] - glgen[k][GIND(i,j)]) / (Xhgen[k][k] - Xlgen[k][k]);
2166 
2167  // dualfprintf(fail_file,"ii=%d jj=%d kk=%d :: i=%d j=%d k=%d c=%21.15g gh=%21.15g gl=%21.15g Xh=%21.15g Xl=%21.15g\n",localptrgeom->i,localptrgeom->j,localptrgeom->k,i,j,k,conn[i][j][k],ghgen[k][GIND(i,j)],glgen[k][GIND(i,j)],Xhgen[k][k],Xlgen[k][k]);
2168 
2169  }
2170  }
2171 
2172  }// end over k
2173 
2174 
2175  }
2176 
2177 
2178 
2179 
2180 
2181 
2182  // DIFFNUMREC version has DIFFGAMMIE as default value
2183  if(CONNDERTYPE==DIFFNUMREC){
2184 
2185 
2187  //
2188  // Setup delta
2189  //
2191  // 0 indicates connection type
2192  setup_delta(0,CONNDERTYPE,DELTA,geom,localptrgeoml,localptrgeomh,truedelta);
2193 
2194 
2195  // use local copy so don't overwrite original, which could be pointing to global storage
2196  localptrgeom->i=geom->i;
2197  localptrgeom->j=geom->j;
2198  localptrgeom->k=geom->k;
2199  localptrgeom->p=NOWHERE; // informs rest of calls that X will generally be arbitrary
2200 
2201  // dualfprintf(fail_file,"DIFFNUMREC: doing i=%d j=%d\n",geom->i,geom->j);
2202  for (k = 0; k < NDIM; k++) {
2203 
2204 
2205  if(WHICHEOM!=WITHGDET){
2206  failreturn = dfridr(lngdet_func_mcoord,localptrgeom,X,0,0,k,&value); // 0,0 not used
2207  if(failreturn==0) conn2[k]=value; // else leave as default
2208  }
2209  else conn2[k]=0.0; // then no 2nd connection
2210 
2211  // answer is symmetric on i,j since uses g_{ij}, so only do part of work
2212  for (i = 0; i < NDIM; i++){
2213  for (j = 0; j <=i; j++){
2214  failreturn = dfridr(gcov_func_mcoord,localptrgeom,X,i,j,k,&value);
2215  if(failreturn==0) conn[i][j][k]=value; // else leave as default
2216  }
2217  }
2218 
2219  }// end over k
2220 
2221  }// end if CONNDERTYPE==DIFFNUMREC
2222 
2223 
2224 
2225 
2226 
2227 
2228 
2229 
2231  //
2232  // fill in rest of conn[i][j][k] (or enforce symmetry of connection)
2233  //
2235  for (k = 0; k < NDIM; k++) {
2236  for (i = 0; i < NDIM; i++){
2237  for (j = i+1; j <NDIM; j++){
2238  conn[i][j][k] = conn[j][i][k];
2239  }
2240  }
2241  }
2242 
2243 
2245  //
2246  // now rearrange to find \Gamma_{ijk}=1/2*(gabc_{jik}+gabc_{kij}-gabc_{kji})
2247  //
2249  for (i = 0; i < NDIM; i++)
2250  for (j = 0; j < NDIM; j++)
2251  for (k = 0; k < NDIM; k++)
2252  tmp[i][j][k] =
2253  0.5 * (conn[j][i][k] + conn[k][i][j] - conn[k][j][i]);
2254 
2255 
2257  //
2258  // finally, raise first index
2259  //
2261  for (i = 0; i < NDIM; i++)
2262  for (j = 0; j < NDIM; j++)
2263  for (k = 0; k < NDIM; k++) {
2264  conn[i][j][k] = 0.;
2265  for (l = 0; l < NDIM; l++){
2266  conn[i][j][k] += geom->gcon[GIND(i,l)] * tmp[l][j][k];
2267  }
2268  }
2269 
2270 
2272  //
2273  // now correct for accurate body forces
2274  // only makes sense for no a2c on flux right now
2275  //
2276  // Idea is that pressure term in stress-energy tensor does not cancel between flux and source term, leading to lack of cancellation leading to secular errors -- especially near pole where flux must vanish so errors in flux are not removed by flux differencing.
2277  //
2278  // So by looking at $T^\mu_\nu += p_{\rm tot} \delta^\mu_\nu$ term, one can see how to correct the connection so this pressure term (for constant total pressure) cancels between flux and source
2279  //
2280  // --> $d_t (\detg \delta^t_\nu) = -d_j(\detg \delta^j_\nu) + \detg \delta^k_\lambda \Gamma^\lambda_{\nu\kappa}$
2281  //
2282  // assume $d_t(\detg)\sim 0$ and $d_j(ptot)\sim 0$, then:
2283  //
2284  // $0 = -d_\nu(\detg) + \detg \Gamma^\kappa_{\nu\kappa}$
2285  // or:
2286  // $\Gamma^\kappa_{\nu\kappa} = d_\nu(\detg) / \detg$
2287  //
2288  // Since source is at center while flux is at faces, we need to subtract off face-related values and add center-related values. The $d_\nu(\detg)$ term is really the only face part related to the flux calculation. That needs to be removed and then we should add back on the center version
2289  //
2290  // 1) $[each \kappa]\Gamma^\kappa_{\nu\kappa} -= (1/4) (d_\nu(\detg) @ cent / (\detg @ cent)) = [sum over \kappa] (1/4) \Gamma^\kappa_{\nu\kappa}$
2291  // 2) $[each \kappa]\Gamma^\kappa_{\nu\kappa} += (1/4) (\delta_\nu(\detg)/(\Delta_\nu) / (\detg @ cent))$
2292  //
2293  // No, cannot just completely change each \kappa term like this, since could change dramatically the value...
2294  //
2295  // Instead, Form Q_\nu = ([wanted version, with sum over \kappa] \Gamma^\kappa_{\nu\kappa}) / ([original, with sum over \kappa] \Gamma^\kappa_{\nu\kappa})
2296  //
2297  // Then to get final \Gamma^\kappa_{\nu\kappa}, multiply *each* \kappa term by Q_\nu
2298  // Then one has minimal multiplicative factor that multiplies each term so that sum will be desired value
2299  //
2301  if(doingmachinebody){
2302 
2304  //
2305  // First repeat setup for connection calculation but use DIFFFINITE so that metric quantities are evaluated at FACES rather than CENT
2306  //
2308  if(CONNDERTYPE!=DIFFFINITE){
2309  // Setup finite difference for correction
2310  // 0 indicates connection type
2311  // correction always uses DIFFFINITE
2312  setup_delta(0,DIFFFINITE,DELTA,geom,localptrgeoml,localptrgeomh,truedelta);
2313 
2314  // setup Xl and Xh and signdX
2315  setup_XlXh(X,truedelta,Xlgen,Xhgen,signdXgen);
2316 
2317  // 0 means not normal calculations
2318  compute_metricquantities_midlh(0,doingmachinebody
2319  ,geom, X, gmid, gcovpertmid,&gdetmid,gdetlgen,gdethgen
2320  ,localptrgeoml,localptrgeomh,Xlgen,Xhgen
2321  ,lngdetlgen, lngdethgen, glgen, ghgen, gcovpertlgen, gcovperthgen
2322  );
2323 
2324  // resetup Xl and Xh and signdX since above "compute" function will have overwritten X positions for non-existent dimensions
2325  // if do this, must set truedelta correctly consistent in setup_delta()
2326  setup_XlXh(X,truedelta,Xlgen,Xhgen,signdXgen);
2327 
2328  }
2329 
2331  // form original contracted connection
2333  DLOOPA(kk){
2334  conndiag[kk]=0.0;
2335  DLOOPA(jj) conndiag[kk] += conn[jj][kk][jj];
2336  }
2337 
2339  // form new finite differential \detg divided by centered \detg
2341  for (k = 0; k < NDIM; k++) {
2342  conndiag2[k] = signdXgen[k]*(gdethgen[k] - gdetlgen[k]) / (Xhgen[k][k] - Xlgen[k][k]);
2343  conndiag2[k] /= gdetmid;
2344  }
2345 
2347  // now obtain correction
2348  // for Xh-Xl->0 this vanishes as required
2349  // Plugging this new conn into EOM for T^\mu_\nu = p \delta^\mu_\nu gives exactly cancellation between source and flux differencing of pressure
2350 #if(0)
2351  // Note that correction applies to \Gamma^\kappa_{\nu\kappa}, which is a contracted quantity. So we spread correction across each \kappa=0,1,2,3. Hence the 0.25
2352  // Once later contractions operate and pressure and flux source terms appear, the contracted term involving the pressure will cancel correctly for constant pressure
2354  // for (i = 0; i < NDIM; i++) for (j = 0; j < NDIM; j++) for (k = 0; k < NDIM; k++) conn[i][j][k] += -0.25*(conndiag[j] - conndiag2[j])*delta(i,k);
2355  // apply delta(i,k) directly by setting i=k
2356  // for (i = 0; i < NDIM; i++) for (j = 0; j < NDIM; j++) conn[i][j][i] += -0.25*(conndiag[j] - conndiag2[j]);
2357  // for (i = 0; i < NDIM; i++) for (j = 0; j < NDIM; j++) conn[i][j][i] *= pow(fabs(conndiag2[j])/(fabs(conndiag[j])+SMALL),0.25);
2358  FTYPE ftemp;
2359  for (i = 0; i < NDIM; i++){
2360  for (j = 0; j < NDIM; j++){
2361  if(fabs(conndiag2[j])==0.0) ftemp=1.0;
2362  else ftemp=(fabs(conndiag2[j])+SMALL)/(fabs(conndiag[j])+SMALL);
2363  // Note, there is difficulty when sum passes through zero. Won't matter for pressure term, but each individual term in connection may be far from zero and simply different terms cancel.
2364  if(fabs(ftemp-1.0)>0.5){
2365  dualfprintf(fail_file,"WARNING: Large correction for machinebody: i=%d j=%d :: i=%d j=%d ftemp=%21.15g :: %21.15g %21.15g :: %21.15g\n",geom->i,geom->j,i,j,ftemp,conndiag[j],conndiag2[j],gdetmid);
2366  }
2367  else{
2368  // otherwise do correction
2369  conn[i][j][i] *= ftemp;
2370  // if(i==0) dualfprintf(fail_file,"machinebody: i=%d j=%d :: j=%d ftemp=%21.15g :: %21.15g %21.15g :: %21.15g\n",geom->i,geom->j,j,ftemp,conndiag[j],conndiag2[j],gdetmid);
2371  }
2372  }
2373  }
2374 #else
2375  // Sasha's weighted method for correction
2376  FTYPE sumabsconn[NDIM];
2377  DLOOPA(kk){
2378  sumabsconn[kk]=SMALL; // to avoid 0/0 in weight
2379  DLOOPA(jj) sumabsconn[kk] += fabs(conn[jj][kk][jj]);
2380  }
2381 
2382  FTYPE dS,weight;
2383  for (i = 0; i < NDIM; i++){ // over traced terms
2384  for (j = 0; j < NDIM; j++){ // over each j
2385  // correction to sum of trace:
2386  dS=conndiag2[j]-conndiag[j];
2387  // weight for this connection term
2388  weight=fabs(conn[i][j][i])/sumabsconn[j];
2389  // weighted correction
2390  conn[i][j][i] += dS*weight;
2391  }
2392  }
2393 #endif
2394 
2395 
2396  } // end if correcting body forces
2397 
2398 
2399 
2400 
2401 #if(0) // DEBUG
2402  if(geom->i==0 || geom->i==-1){
2403  for (i = 0; i < NDIM; i++) for (l = 0; l < NDIM; l++){
2404  dualfprintf(fail_file,"i=%d gcon[%d][%d]=%21.15g\n",geom->i,i,l,geom->gcon[GIND(i,l)]);
2405  }
2406  dualfprintf(fail_file,"i=%d conn[0][0][0]=%21.15g\n",geom->i,conn[0][0][0]);
2407  }
2408 #endif
2409 
2410  /* done! */
2411 }
2412 
2413 
2414 
2415 int setup_XlXh(FTYPE *X,FTYPE *truedelta, FTYPE (*Xlgen)[NDIM],FTYPE (*Xhgen)[NDIM],FTYPE *signdXgen)
2416 {
2417  int k,l;
2418 
2420  // first form high and low positions for function locations
2422  for (k = 0; k < NDIM; k++) {
2423 
2424  for (l = 0; l < NDIM; l++){
2425  Xhgen[k][l] = X[l];
2426  Xlgen[k][l] = X[l];
2427  }
2428 
2429  // normal case
2430  Xhgen[k][k] += MYDELTAh(truedelta[k],k);
2431  Xlgen[k][k] += MYDELTAl(truedelta[k],k);
2432  signdXgen[k]=1.0;
2433 
2434 #if(0)// debug
2435  if(k==TT){
2436  if(Xlgen[k][k]>X[k]){
2437  dualfprintf(fail_file,"Xl in future! Xl=%21.15g Xh=%21.15g\n",Xlgen[k][k],Xhgen[k][k]);
2438  }
2439  }
2440 #endif
2441 
2442 
2443 #if(0)
2444  // not really wanted since want "force" to be in same "radial" direction so sign SHOULD flip
2445  if(k==1 && ISSPCMCOORD(MCOORD)){
2446  // then check if r<0 and invert Xh and Xl if so
2447  bl_coord(X, V);
2448  if(V[k]<0.0){
2449  signdXgen[k]=-1.0;
2450  }
2451  }
2452 #endif
2453  // if(k==TT) dualfprintf(fail_file,"k=%d DELTAl=%21.15g DELTAh=%21.15g : true=%21.15g\n",k,MYDELTAl(truedelta[k],k),MYDELTAh(truedelta[k],k),truedelta[k]);
2454  // DEBUG
2455  //Xhgen[k][k] += DELTA;
2456  // Xlgen[k][k] -= DELTA;
2457 
2458  }
2459 
2460  return(0);
2461 
2462 }
2463 
2464 
2465 
2466 // compute low/high metric quantities
2467 int compute_metricquantities_midlh(int donormal, int domachinebody
2468  ,struct of_geom *geom, FTYPE *X, FTYPE *gmid, FTYPE *gcovpertmid,FTYPE *gdetmid,FTYPE *gdetlgen,FTYPE *gdethgen
2469  ,struct of_geom (*localptrgeoml)[NDIM],struct of_geom (*localptrgeomh)[NDIM],FTYPE (*Xlgen)[NDIM],FTYPE (*Xhgen)[NDIM]
2470  ,FTYPE *lngdetlgen, FTYPE *lngdethgen, FTYPE (*glgen)[SYMMATRIXNDIM], FTYPE (*ghgen)[SYMMATRIXNDIM], FTYPE (*gcovpertlgen)[NDIM], FTYPE (*gcovperthgen)[NDIM]
2471  )
2472 {
2473  int k;
2474  FTYPE gcov_func_mcoord(struct of_geom *ptrgeom, FTYPE* X, int i, int j);
2475  FTYPE lngdet_func_mcoord(struct of_geom *ptrgeom, FTYPE* X, int i, int j);
2476  FTYPE gdet_func_mcoord_usegcov(FTYPE *gcovmcoord, struct of_geom *ptrgeom, FTYPE* X, int i, int j);
2477 
2478  // get gdet in cell center
2479  if(domachinebody){
2480  gcov_func(geom,1,MCOORD,X, gmid,gcovpertmid);
2481  *gdetmid=gdet_func_mcoord_usegcov(gmid, geom, X, 0,0);
2482  }
2483 
2484  // get k-dependent things
2485  for (k = 0; k < NDIM; k++) {
2486 
2487  if(donormal){
2488  if(WHICHEOM!=WITHGDET){
2489  lngdethgen[k]=lngdet_func_mcoord(&((*localptrgeomh)[k]),Xhgen[k],0,0); // doesn't use 0,0
2490  lngdetlgen[k]=lngdet_func_mcoord(&((*localptrgeoml)[k]),Xlgen[k],0,0); // doesn't use 0,0
2491  }
2492  }
2493  if(donormal || domachinebody){
2494  // must come before below gdet_func_mcoord_usegcov() call
2495  gcov_func(&((*localptrgeomh)[k]),1,MCOORD,Xhgen[k], ghgen[k],gcovperthgen[k]);
2496  gcov_func(&((*localptrgeoml)[k]),1,MCOORD,Xlgen[k], glgen[k],gcovpertlgen[k]);
2497  }
2498 
2499  if(domachinebody){
2500  gdethgen[k]=gdet_func_mcoord_usegcov(ghgen[k], &((*localptrgeomh)[k]), Xhgen[k], 0,0);
2501  gdetlgen[k]=gdet_func_mcoord_usegcov(glgen[k], &((*localptrgeoml)[k]), Xlgen[k], 0,0);
2502  }
2503 
2504  }
2505 
2506  return(0);
2507 
2508 }
2509 
2510 
2511 
2512 
2513 
2514 // returns MCOORD gcov value for i,j element
2515 // excessive to compute other elements, but ok for now
2516 // input ptrgeom is only used for i,j,k,loc positions
2517 FTYPE gcov_func_mcoord(struct of_geom *ptrgeom, FTYPE* X, int i, int j)
2518 {
2519  FTYPE gcovmcoord[SYMMATRIXNDIM];
2520  FTYPE gcovpert[NDIM];
2521  void gcov_func(struct of_geom *ptrgeom, int getprim, int whichcoord, FTYPE *X, FTYPE *gcov, FTYPE *gcovpert);
2522 
2523  gcov_func(ptrgeom,1,MCOORD,X,gcovmcoord,gcovpert);
2524  // if(i==j) return(gcovpert[i]); // for accurate non-rel gravity (works since d(1+g_{tt}) = dg_{tt})
2525  // else
2526  return(gcovmcoord[GIND(i,j)]);
2527 }
2528 
2529 
2530 // returns MCOORD value for log(f/gdet). Doesn't use i,j (these are not grid locations)
2531 // input ptrgeom is only used for i,j,k,loc positions
2532 FTYPE lngdet_func_mcoord(struct of_geom *ptrgeom, FTYPE* X, int i, int j)
2533 {
2534  FTYPE gcovmcoord[SYMMATRIXNDIM];
2535  FTYPE gcovpertcoord[NDIM];
2536  FTYPE toreturn;
2537  FTYPE eomfunc[NPR]; // don't change this as related to WHICHEOM==WITHGDET
2538  FTYPE *eomfuncptr=eomfunc; // don't change this as related to WHICHEOM==WITHGDET
2539  FTYPE V[NDIM];
2540  void gcov_func(struct of_geom *ptrgeom, int getprim, int whichcoord, FTYPE *X, FTYPE *gcov, FTYPE *gcovpert);
2541  FTYPE gdet;
2542 
2543  gcov_func(ptrgeom,1,MCOORD,X,gcovmcoord,gcovpertcoord);
2544  eomfunc_func(ptrgeom,1,MCOORD,X,EOMFUNCPTR);
2545 
2546  bl_coord_ijk_2(ptrgeom->i,ptrgeom->j,ptrgeom->k,ptrgeom->p, X, V);
2547  // GODMARK: assumes all RHO, etc. use same eomfunc if using 2nd connection
2548  if(gdet_func_metric(MCOORD,V,gcovmcoord,&gdet)!=0){
2549  if(debugfail>=2) dualfprintf(fail_file,"Caught gdet_func_metric() issue in lngdet_func_mcoord()\n");
2550  }
2551  toreturn=log(fabs(EOMFUNCMAC(RHO)/gdet)); // can't take log of negative (FLIPGDETAXIS note about how this method wouldn't work)
2552 
2553  return(toreturn);
2554 }
2555 
2556 // returns MCOORD value for gdet. Doesn't use i,j (these are not grid locations)
2557 FTYPE gdet_func_mcoord(struct of_geom *ptrgeom, FTYPE* X, int i, int j)
2558 {
2559  FTYPE gcovmcoord[SYMMATRIXNDIM];
2560  FTYPE gcovpertcoord[NDIM];
2561  FTYPE toreturn;
2562  FTYPE V[NDIM];
2563  void gcov_func(struct of_geom *ptrgeom, int getprim, int whichcoord, FTYPE *X, FTYPE *gcov, FTYPE *gcovpert);
2564  FTYPE gdet;
2565 
2566  gcov_func(ptrgeom,1,MCOORD,X,gcovmcoord,gcovpertcoord);
2567 
2568  bl_coord_ijk_2(ptrgeom->i,ptrgeom->j,ptrgeom->k,ptrgeom->p, X, V);
2569  if(gdet_func_metric(MCOORD,V,gcovmcoord,&gdet)!=0){
2570  if(debugfail>=2) dualfprintf(fail_file,"Caught gdet_func_metric() issue in gdet_func_mcoord()\n");
2571  }
2572 
2573  return(gdet);
2574 }
2575 
2576 // returns MCOORD value for gdet using gcovmcoord as input to avoid repeated computations of gcovmcoord. Doesn't use i,j (these are not grid locations)
2577 FTYPE gdet_func_mcoord_usegcov(FTYPE *gcovmcoord, struct of_geom *ptrgeom, FTYPE* X, int i, int j)
2578 {
2579  // FTYPE gcovmcoord[SYMMATRIXNDIM];
2580  // FTYPE gcovpertcoord[NDIM];
2581  FTYPE toreturn;
2582  FTYPE V[NDIM];
2583  //void gcov_func(struct of_geom *ptrgeom, int getprim, int whichcoord, FTYPE *X, FTYPE *gcov, FTYPE *gcovpert);
2584  FTYPE gdet;
2585 
2586  // gcov_func(ptrgeom,1,MCOORD,X,gcovmcoord,gcovpertcoord);
2587 
2588  bl_coord_ijk_2(ptrgeom->i,ptrgeom->j,ptrgeom->k,ptrgeom->p, X, V);
2589  if(gdet_func_metric(MCOORD,V,gcovmcoord,&gdet)!=0){
2590  if(debugfail>=2) dualfprintf(fail_file,"Caught gdet_func_metric() issue in gdet_func_mcoord_usegcov()\n");
2591  }
2592 
2593  return(gdet);
2594 }
2595 
2596 
2597 
2598 
2599 
2600 
2601 
2602 
2603 
2604 
2605 
2606 
2607 
2608 
2609 
2610 
2611 
2612 
2613 
2614 
2615 
2616 int metric_checks(struct of_geom *ptrgeom)
2617 {
2618  FTYPE delta[NDIM][NDIM];
2619  FTYPE V[NDIM],X[NDIM];
2620  int i,j,k,loc;
2621  int jj,kk,pp;
2622 
2623 
2624  i=ptrgeom->i;
2625  j=ptrgeom->j;
2626  k=ptrgeom->k;
2627  loc=ptrgeom->p;
2628 
2630  //
2631  // delta^jj_kk = gcov_{kk pp} gcon^{pp jj}
2632  //
2634  DLOOP(jj,kk){
2635  delta[jj][kk]=0.0;
2636  DLOOPA(pp){
2637  delta[jj][kk]+= ptrgeom->gcov[GIND(jj,pp)]* ptrgeom->gcon[GIND(pp,kk)];
2638  }
2639 
2640  if(PRODUCTION==0){
2641  if(fabs(delta[jj][kk]-delta(jj,kk))>NUMEPSILON*100.0){
2642  if(ISSPCMCOORD(MCOORD)==0 || (ISSPCMCOORD(MCOORD)==1 && j!=totalsize[2] && j!=0 && loc!=FACE2 && loc!=CORN1 && loc!=CORN3 && loc!=CORNT) ){
2643  dualfprintf(fail_file,"Problem with metric at i=%d j=%d k=%d loc=%d delta[%d][%d]=%21.15g should be %21.15g\n",ptrgeom->i,ptrgeom->j,ptrgeom->k,ptrgeom->p,jj,kk,delta[jj][kk],delta(jj,kk));
2644  }
2645  }
2646  }
2647 
2648  if(fabs(delta[jj][kk]-delta(jj,kk))>NUMEPSILON*1000.0){
2649  if(ISSPCMCOORD(MCOORD)==0 || (ISSPCMCOORD(MCOORD)==1 && j!=totalsize[2] && j!=0 && loc!=FACE2 && loc!=CORN1 && loc!=CORN3 && loc!=CORNT) ){
2650  dualfprintf(fail_file,"MAJOR Problem with metric at i=%d j=%d k=%d loc=%d delta[%d][%d]=%21.15g should be %21.15g\n",ptrgeom->i,ptrgeom->j,ptrgeom->k,ptrgeom->p,jj,kk,delta[jj][kk],delta(jj,kk));
2651  }
2652  }
2653  }
2654 
2655 
2656 
2657  // check if near static limit since can't divide by the below in ucon_calc
2658  // GODMARK
2659  if (fabs(ptrgeom->gcon[GIND(TT,TT)]) < SMALL) {
2660  bl_coord_ijk_2(i,j,k,loc,X,V);
2661  dualfprintf(fail_file, "grid location too near g_{tt}==0: %d %d %d : r=%21.15g th=%21.15g phi=%21.15g : Rin=%21.15g %21.15g\n", i,j,k,V[1],V[2],V[3],Rin,ptrgeom->gcon[GIND(TT,TT)]);
2662  myexit(1);
2663  }
2664  if (0 && fabs(ptrgeom->gcon[GIND(RR,RR)]) < SMALL) {
2665  bl_coord_ijk_2(i,j,k,loc,X,V);
2666  dualfprintf(fail_file, "grid location too near g^{rr}==0: %d %d %d : r=%21.15g th=%21.15g phi=%21.15g : Rin=%21.15g %21.15g\n", i,j,k,V[1],V[2],V[3],Rin,ptrgeom->gcon[GIND(RR,RR)]);
2667  myexit(1);
2668  }
2669  if (0 && fabs(ptrgeom->gcon[GIND(TH,TH)]) < SMALL) {
2670  bl_coord_ijk_2(i,j,k,loc,X,V);
2671  dualfprintf(fail_file,"grid location too near g^{\\theta\\theta}==0: %d %d %d : r=%21.15g th=%21.15g phi=%21.15g : Rin=%21.15g %21.15g\n", i,j,k,V[1],V[2],V[3],Rin,ptrgeom->gcon[GIND(TH,TH)]);
2672  myexit(1);
2673  }
2674  if (0 && fabs(ptrgeom->gcon[GIND(PH,PH)]) < SMALL) {
2675  bl_coord_ijk_2(i,j,k,loc,X,V);
2676  dualfprintf(fail_file,"grid location too near g^{\\phi\\phi}==0: %d %d %d : r=%21.15g th=%21.15g phi=%21.15g : Rin=%21.15g %21.15g\n", i,j,k,V[1],V[2],V[3],Rin,ptrgeom->gcon[GIND(PH,PH)]);
2677  myexit(1);
2678  }
2679 
2680  // what about g_{tt}==0? Do I ever divide by g_{tt}?
2681  // Yes, for ucon[TT] for 4 velocity, which is done if WHICHVEL==VEL4 or init.c
2682  // what about g^{rr}==0? Do I ever divide by g^{rr}?
2683  // Doesn't appear so
2684  // g^{pp} inf taken care of in metric.c by avoiding theta=0,Pi
2685 
2686 
2687  return(0);
2688 
2689 }
2690 
2691 
2692 
2693 void check_rmin(void)
2694 {
2695  int i,j,k;
2696  FTYPE X[NDIM],V[NDIM],r;
2697 
2698 
2699  // diagnostic
2700  // determine nature of inner radial edge (assumes myid==0 is always there)
2701  if(myid==0){
2702  i=INFULL1;
2703  j=k=0;
2704  coord(i,j,k, FACE1, X);
2705  bl_coord(X, V);
2706  r=V[1];
2707  trifprintf("rmin(i=%d,X=%21.15g) = %21.15g\n", i,X[1],r);
2708  trifprintf("r=%21.15g Rhor=%21.15g :: rmin/rh: %21.15g\n",r,Rhor, r / (fabs(Rhor)+SMALL) );
2709  // trifprintf("rmin/rsing: %21.15g\n", r / (a+SMALL));
2710  if(r/(fabs(Rhor)+SMALL)<=1.0){
2711  trifprintf("inner grid is inside horizon\n");
2712  }
2713  else{
2714  trifprintf("inner grid is outside horizon\n");
2715  }
2716  }
2717 
2718  // show cells that are inside horizon
2719  if(mycpupos[2]==ncpux2/2 && mycpupos[3]==0){
2720  j=0;
2721  k=0;
2722  LOOPF1{
2723  coord(i,j,k, FACE1, X);
2724  bl_coord(X, V);
2725  r=V[1];
2726  if(r<Rhor){
2727  logfprintf("INSIDE Horizon (r=%g): i=%d r=%21.15g\n",Rhor, i, r);
2728  }
2729 
2730  }
2731  }
2732 
2733 
2734 
2735 
2736 }
Generated on Fri May 20 2016 15:52:34 for HARM by   doxygen 1.8.3.1