HARM
harm and utilities
 All Data Structures Files Functions Variables Typedefs Macros Pages
set_grid.c
Go to the documentation of this file.
1 
7 #include "decs.h"
8 
9 
10 static void set_position_stores(void);
11 static void set_grid_metrics(void);
12 static void set_grid_metrics_withsymmetrization(void);
13 static void set_connection(void);
14 static void set_idxvol(void);
15 static void symmetrize_connection(void);
16 static void symmetrize_gcov(void);
17 static void set_grid_metrics_gcov(void);
18 static void set_grid_metrics_others(void);
19 static void symmetrize_X_V_dxdxp_idxdxp(void);
20 static void set_tlab2ortho(void);
21 
22 
24 #define ATTEMPTSYMMETRIZATION 0
25 
27 int assignmetricstorage_new(struct of_compgeom *mygeom, FTYPE **localgcov, FTYPE **localgcon, FTYPE **localgcovpert, FTYPE **localgdet, FTYPE **localgdetvol, FTYPE **localalphalapse, FTYPE **localbetasqoalphasq, FTYPE **localbeta, FTYPE **localeomfunc)
28 {
29 
30 #if(NEWMETRICSTORAGE==1)
31 
32 
33  *localgcov=mygeom->gcov;
34  *localgcon=mygeom->gcon;
35  *localgcovpert=mygeom->gcovpert;
36  *localgdet=&(mygeom->gdet);
37 #if(GDETVOLDIFF)
38  *localgdetvol=&(mygeom->gdetvol);
39 #endif
40  *localalphalapse=&(mygeom->alphalapse);
41  *localbetasqoalphasq=&(mygeom->betasqoalphasq);
42  *localbeta=mygeom->beta;
43  *localeomfunc=&(mygeom->EOMFUNCMAC(0));
44 
45 
46 #endif
47  return(0);
48 }
49 
50 
53 int assignmetricstorage_old(int loc, int i, int j, int k, FTYPE **localgcov, FTYPE **localgcon, FTYPE **localgcovpert, FTYPE **localgdet, FTYPE **localgdetvol, FTYPE **localalphalapse, FTYPE **localbetasqoalphasq, FTYPE **localbeta, FTYPE **localeomfunc)
54 {
55 #if(NEWMETRICSTORAGE==0)
56 
57 
58  *localgcov=GLOBALMETMACP1A0(gcov,loc,i,j,k);
59  *localgcon=GLOBALMETMACP1A0(gcon,loc,i,j,k);
60  *localgcovpert=GLOBALMETMACP1A0(gcovpert,loc,i,j,k);
61  *localgdet=&GLOBALMETMACP1A0(gdet,loc,i,j,k);
62 #if(GDETVOLDIFF)
63  *localgdetvol=&GLOBALMETMACP1A0(gdetvol,loc,i,j,k);
64 #endif
65  *localalphalapse=&GLOBALMETMACP1A0(alphalapse,loc,i,j,k);
66  *localbetasqoalphasq=&GLOBALMETMACP1A0(betasqoalphasq,loc,i,j,k);
67  *localbeta=GLOBALMETMACP1A0(beta,loc,i,j,k);
68 
69 #if(WHICHEOM!=WITHGDET)
70  *localeomfunc=GLOBALMETMACP1A0(eomfunc,loc,i,j,k);
71 #else
72  *localeomfunc=&GLOBALMETMACP1A0(gdet,loc,i,j,k)
73 #endif
74 
75 
76 #endif
77 
78  return(0);
79 }
80 
83 int assignmetricstorage_oldlast(int loc, int i, int j, int k, FTYPE **localgcov, FTYPE **localgcon, FTYPE **localgcovpert, FTYPE **localgdet, FTYPE **localgdetvol, FTYPE **localalphalapse, FTYPE **localbetasqoalphasq, FTYPE **localbeta, FTYPE **localeomfunc)
84 {
85 #if(NEWMETRICSTORAGE==0)
86 
87  *localgcovlast=GLOBALMETMACP1A0(gcovlast,loc,i,j,k);
88  // *localgconlast=GLOBALMETMACP1A0(gconlast,loc,i,j,k);
89  *localgcovpertlast=GLOBALMETMACP1A0(gcovpertlast,loc,i,j,k);
90  // *localgdetlast=&GLOBALMETMACP1A0(gdetlast,loc,i,j,k);
91 #if(GDETVOLDIFF)
92  // *localgdetvollast=&GLOBALMETMACP1A0(gdetvollast,loc,i,j,k);
93 #endif
94  *localalphalapselast=&GLOBALMETMACP1A0(alphalapselast,loc,i,j,k);
95  // *localbetasqoalphasqlast=&GLOBALMETMACP1A0(betasqoalphasqlast,loc,i,j,k);
96  // *localbetalast=GLOBALMETMACP1A0(betalast,loc,i,j,k);
97  //#if(WHICHEOM!=WITHGDET)
98  // *localeomfunclast=GLOBALMETMACP1A0(eomfunclast,loc,i,j,k);
99  //#else
100  // *localeomfunclast=&GLOBALMETMACP1A0(gdetlast,loc,i,j,k)
101  //#endif
102 
103 
104 #endif
105 
106  return(0);
107 }
108 
109 
110 
111 
120 void set_grid(int whichtime, FTYPE *CUf, FTYPE *Cunew)
121 {
122  extern void set_rvsr(void);
123  extern void control_time_store_metric(int whichtime, FTYPE *Cunew);
124  extern void set_drsing(void);
125  extern void set_rvsr(void);
126  extern int bound_spacetime_inside_horizon(void);
127  extern int store_old_metric(void);
128 
129  int pliter,pl;
130 
131  // NOTE:
132  // don't assume we enter here with whichtime==0 only once since may have to "iteratively" obtain metric from stress-energy tensor and metric
133  // this happens at t=0 in init() when DOSELFGRAVVSR==1
134 
135 
136 
137  // choose time of metric and whether to store metric before overwritten
138  control_time_store_metric(whichtime,Cunew);
139 
140 
141  /* set up boundaries, steps in coordinate grid */
142  if(whichtime==0 && dt==0.0){
143  // && DOEVOLVEMETRIC){ // no, can't add this conditional. need dt to be certainly non-zero in general for connection.
144  dt=1.0; // dummy value that should lead to 0 connection coefficient for derivatives in time
145  // just avoids 0/0 that should really be 0
146  }
147  // set_points() is purely dealing with coordinates, not metric quantities
148  // but sets dx[0]=dt
149  set_points();
150 
151 
152 
154  //
155  // first compute coordinate labels, which never change in time right now
156  //
158  if(whichtime==0){
159  // set minimum dr used to smooth metric (requires set_points())
160  // assume drsing doesn't change (i.e. coordinates don't change)
161  trifprintf("set_drsing() START\n");
162  set_drsing();
163  trifprintf("set_drsing() END\n");
164  trifprintf("drsing=%21.15g\n",drsing);
165 
166  if(DOSELFGRAVVSR){
167  // set rvsr
168  // assume r doesn't change with time
169  trifprintf("Setting r(i)\n");
170  set_rvsr();
171  trifprintf("Done setting r(i)\n");
172  }
173 
174 
175  if(ATTEMPTSYMMETRIZATION){
176  trifprintf("Will attempt to symmetrize coordinate/metric/connection type quantities\n");
177  }
178 
179 
181  //
182  // set positional stores
183  //
185 #if(DOSTOREPOSITIONDATA)
186  trifprintf("set_position_stores() BEGIN\n");
187  set_position_stores();
188  didstorepositiondata=1;
189  trifprintf("set_position_stores() END\n");
190 #endif // end if DOSTOREPOSITIONDATA==1
191 
192  }// end if whichtime==0
193 
194 
195 
196 
197 
199  //
200  // set metric
201  //
203 
204  // if here, then doing over again
205  if(whichtime==0) trifprintf("set_grid_metrics() BEGIN\n");
206  didstoremetricdata=0;
207  if(ATTEMPTSYMMETRIZATION==0){
208  set_grid_metrics();
209  }
210  else{
211  set_grid_metrics_withsymmetrization();
212  }
213  // if here, then did store new metric data
214  didstoremetricdata=1;
215  if(whichtime==0) trifprintf("set_grid_metrics() END\n");
216 
217 
218 
220  //
221  // set connection
222  //
224  if(whichtime==0){
225  trifprintf("set_connection() BEGIN\n");
226  if(CONNDERTYPE==DIFFNUMREC){
227  trifprintf("set_connection() will take a while with CONNDERTYPE==DIFFNUMREC. Be patient ...\n");
228  }
229  }
230  set_connection();
231  if(whichtime==0) trifprintf("set_connection() END\n");
232 
233 
235  //
236  // set idxvol
237  //
239 
240  if(VOLUMEDIFF){
241  if(whichtime==0) trifprintf("set_idxvol() BEGIN\n");
242  set_idxvol();
243  if(whichtime==0) trifprintf("set_idxvol() END\n");
244  }
245 
246 
247  // set boundary conditions on metric
248  // this only modifies outside unused computational regions, so can come last. One might imagine that bounding metric would change connection calculation...true but connection bounded too. +-1 issues? GODMARK
249  if(DOEVOLVEMETRIC){
250  if(whichtime==0) trifprintf("bound_spacetime_inside_horizon() BEGIN\n");
251  bound_spacetime_inside_horizon();
252  if(whichtime==0) trifprintf("bound_spacetime_inside_horizon() END\n");
253  }
254 
255 
256  if(DOEVOLVEMETRIC && whichtime==0){
257  // if first time to call set_grid() (i.e. set_grid(0)), then store present metric as old metric
258  // if evolving metric, then store old metric before computing new one
259  // store metric, needed to have dg/dt terms in connection coefficients
260  // initially dg/dt = 0
261  // This should come after metric calculation and after connection calculation so that conneciton has old and new metric rather than new and new metric
262  if(whichtime==0) trifprintf("store_old_metric() START\n");
263  trifprintf("Storing old metric initially\n");
264  store_old_metric();
265  if(whichtime==0) trifprintf("store_old_metric() END\n");
266  }
267 
268 
269 
270  if(EOMRADTYPE!=EOMRADNONE && STORETLAB2ORTHO==1){
271  if(whichtime==0) trifprintf("set_tlab2ortho() BEGIN\n");
272  set_tlab2ortho();
273  if(whichtime==0) trifprintf("set_tlab2ortho() END\n");
274  }
275 
276 
277 
278  /* done! */
279 }
280 
281 
282 
283 
284 
285 
287 static void set_position_stores(void)
288 {
289  extern void dxdxprim(FTYPE *X, FTYPE *V, FTYPE (*dxdxp)[NDIM]);
290 
291 
292 
293  // separately store X since over different domain in i,j,k and no MCOORD dependency
294 #pragma omp parallel
295  {
296  int i, j, k;
297  int loc;
298 
299 
300  // over grid locations needing these quantities
301  for (loc = NPG - 1; loc >= 0; loc--) {
302 
303 
304  OPENMP3DLOOPVARSDEFINE; OPENMP3DLOOPSETUPFULLP2;
305 #pragma omp for schedule(OPENMPSCHEDULE(),OPENMPCHUNKSIZE(blocksize)) nowait
306  OPENMP3DLOOPBLOCK{
307  OPENMP3DLOOPBLOCK2IJK(i,j,k);
308 
309  // Want X vs. i,j,k
310 
311  // store X since can be expensive to keep recomputing these things, esp. if bl_coord() involves alot of complicated functions
312  // X must be always allowed to vary in space
313  coord(i, j, k, loc, GLOBALMACP1A0(Xstore,loc,i,j,k));
314 
315  }// end internal loop block
316  }// end over locations
317  }// end parallel region (and implied barrier)
318 
319 
320 
321 
322  // separately store V since over different domain in i,j,k
323 #pragma omp parallel // no ucon or EOS stuff needed
324  {
325  int i, j, k;
326  int loc;
327 
328 
329  // over grid locations needing these quantities
330  for (loc = NPG - 1; loc >= 0; loc--) {
331 
332  OPENMP3DLOOPVARSDEFINE; OPENMP3DLOOPSETUPFULLP2;
333 #pragma omp for schedule(OPENMPSCHEDULE(),OPENMPCHUNKSIZE(blocksize)) nowait
334  OPENMP3DLOOPBLOCK{
335  OPENMP3DLOOPBLOCK2IJK(i,j,k);
336 
337 
338  // Want V vs. i,j,k
339 #if(MCOORD==CARTMINKMETRIC)
340  // doesn't depend on position, only store/use 1 value
341  // if(i!=0 || j!=0 || k!=0) continue; // simple way to avoid other i,j,k when doing OpenMP
342 #endif
343 
344 
345 
346  // store V since can be expensive to keep recomputing these things, esp. if bl_coord() involves alot of complicated functions
347  bl_coord(GLOBALMACP1A0(Xstore,loc,i,j,k),GLOBALMACP1A0(Vstore,loc,i,j,k));
348  // SLOOPA(jj) dualfprintf(fail_file,"loc=%d i=%d j=%d k=%d :: V[%d]=%21.15g\n",loc,i,j,k,jj,GLOBALMACP1A0(Vstore,loc,i,j,k)[jj]);
349  }// end internal loop block
350  }// end over locations
351  }// end parallel region (and implied barrier)
352 
353 
354 
355 
356 
357  // separately do other non-X stores since X has no MCOORD==CARTMINKMETRIC conditional
358 #pragma omp parallel // no ucon or EOS stuff needed
359  {
360  int i, j, k;
361  int loc;
362 
363 
364  // over grid locations needing these quantities
365  for (loc = NPG - 1; loc >= 0; loc--) {
366 
367  OPENMP3DLOOPVARSDEFINE; OPENMP3DLOOPSETUPFULLP1;
368 #pragma omp for schedule(OPENMPSCHEDULE(),OPENMPCHUNKSIZE(blocksize)) nowait
369  OPENMP3DLOOPBLOCK{
370  OPENMP3DLOOPBLOCK2IJK(i,j,k);
371 
372 
373  // dxdxp and idxdxp are only different if non-Minkowski non-Cartesian
374 #if(MCOORD==CARTMINKMETRIC)
375  // doesn't depend on position, only store/use 1 value
376  if(i!=0 || j!=0 || k!=0) continue; // simple way to avoid other i,j,k when doing OpenMP
377 #endif
378 
379 
380 
381  // store dxdxp since can be expensive to keep recomputing these things, esp. if bl_coord() involves alot of complicated functions
382  dxdxprim(GLOBALMACP1A0(Xstore,loc,i,j,k),GLOBALMACP1A0(Vstore,loc,i,j,k),GLOBALMETMACP1A0(dxdxpstore,loc,i,j,k));
383 
384  // SLOOPA(jj) dualfprintf(fail_file,"loc=%d i=%d j=%d k=%d :: V[%d]=%21.15g\n",loc,i,j,k,jj,GLOBALMACP1A0(Vstore,loc,i,j,k)[jj]);
385 
386  matrix_inverse(PRIMECOORDS, GLOBALMETMACP1A0(dxdxpstore,loc,i,j,k),GLOBALMETMACP1A0(idxdxpstore,loc,i,j,k));
387  }// end internal loop block
388  }// end over locations
389  }// end parallel region (and implied barrier)
390 
391 
392 
393 
394  if(ATTEMPTSYMMETRIZATION) symmetrize_X_V_dxdxp_idxdxp();
395 
396 }
397 
398 
399 
400 
402 static void symmetrize_X_V_dxdxp_idxdxp(void)
403 {
404 
405  // symmetrize gcov
406  if(numprocs==1 && ISSPCMCOORD(MCOORD) && BCtype[X2DN]==POLARAXIS && BCtype[X2UP]==POLARAXIS && N2>1){
407 
408 
409 
410  // deal with X,V
411 #pragma omp parallel
412  {
413  int i, j, k;
414  int jj,kk;
415  int loc;
416 
417 
418 
419  // over grid locations needing these quantities
420  for (loc = NPG - 1; loc >= 0; loc--) {
421 
422 
423  if(loc==FACE2 || loc==CORN1 || loc==CORN3 || loc==CORNT){
424 
425  OPENMP3DLOOPVARSDEFINE; OPENMP3DLOOPSETUPFULLP2;
426 #pragma omp for schedule(OPENMPSCHEDULE(),OPENMPCHUNKSIZE(blocksize))
427  OPENMP3DLOOPBLOCK{
428  OPENMP3DLOOPBLOCK2IJK(i,j,k);
429 
430 
431 #if(MCOORD==CARTMINKMETRIC)
432  // doesn't depend on position, only store/use 1 value
433  if(i!=0 || j!=0 || k!=0) continue; // simple way to avoid other i,j,k when doing OpenMP
434 #endif
435 
436 
437  if(j>=N2/2){
438 
439  GLOBALMACP1A1(Vstore,loc,i,j,k,TH)=M_PI-GLOBALMACP1A1(Vstore,loc,i,N2-j,k,TH);
440  GLOBALMACP1A1(Xstore,loc,i,j,k,TH)=endx[TH]-(GLOBALMACP1A1(Xstore,loc,i,N2-j,k,TH)-startx[TH]); //Sasha corrected to work correctly in case startx[TH]!=0
441 
442  }// end over upper hemisphere
443 
444  }// end 3D LOOP
445  }// end if symmetrizing FACE2-like positions
446  else{
447 
448  OPENMP3DLOOPVARSDEFINE; OPENMP3DLOOPSETUPFULLP1EXCEPTX2;
449 #pragma omp for schedule(OPENMPSCHEDULE(),OPENMPCHUNKSIZE(blocksize))
450  OPENMP3DLOOPBLOCK{
451  OPENMP3DLOOPBLOCK2IJK(i,j,k);
452 
453 
454 #if(MCOORD==CARTMINKMETRIC)
455  // doesn't depend on position, only store/use 1 value
456  if(i!=0 || j!=0 || k!=0) continue; // simple way to avoid other i,j,k when doing OpenMP
457 #endif
458 
459  if(j>=N2/2){
460 
461  GLOBALMACP1A1(Vstore,loc,i,j,k,TH)=M_PI-GLOBALMACP1A1(Vstore,loc,i,N2-1-j,k,TH);
462  GLOBALMACP1A1(Xstore,loc,i,j,k,TH)=endx[TH]-(GLOBALMACP1A1(Xstore,loc,i,N2-1-j,k,TH)-startx[TH]); //Sasha corrected to work correctly in case startx[TH]!=0
463 
464 
465  }// end over upper hemisphere
466  }// end 3D LOOP
467  }// end if over CENT-like w.r.t. FACE2
468 
469  }// end over locations
470  }// end parallel region
471 
472 
473 
474 
475  // deal with dxdxp and idxdxp
476 #pragma omp parallel
477  {
478  int i, j, k;
479  int jj,kk;
480  int loc;
481 
482 
483 
484  // over grid locations needing these quantities
485  for (loc = NPG - 1; loc >= 0; loc--) {
486 
487 
488  if(loc==FACE2 || loc==CORN1 || loc==CORN3 || loc==CORNT){
489 
491  OPENMP3DLOOPVARSDEFINE; OPENMP3DLOOPSETUPFULLP1;
492 #pragma omp for schedule(OPENMPSCHEDULE(),OPENMPCHUNKSIZE(blocksize))
493  OPENMP3DLOOPBLOCK{
494  OPENMP3DLOOPBLOCK2IJK(i,j,k);
495 
496 
497 #if(MCOORD==CARTMINKMETRIC)
498  // doesn't depend on position, only store/use 1 value
499  if(i!=0 || j!=0 || k!=0) continue; // simple way to avoid other i,j,k when doing OpenMP
500 #endif
501 
502 
503  if(j>=N2/2){
504  DLOOP(jj,kk){
505  GLOBALMETMACP1A2(dxdxpstore,loc,i,j,k,jj,kk) = GLOBALMETMACP1A2(dxdxpstore,loc,i,N2-j,k,jj,kk);
506  GLOBALMETMACP1A2(idxdxpstore,loc,i,j,k,jj,kk) = GLOBALMETMACP1A2(idxdxpstore,loc,i,N2-j,k,jj,kk);
507  }
508 
509  DLOOP(jj,kk){
510  if(jj==2 && kk!=2 || jj!=2 && kk==2){
511  GLOBALMETMACP1A2(dxdxpstore,loc,i,j,k,jj,kk) *= -1.0;
512  GLOBALMETMACP1A2(idxdxpstore,loc,i,j,k,jj,kk) *= -1.0;
513  }
514  }
515 
516  }// end over upper hemisphere
517 
518  }// end 3D LOOP
519  }// end if symmetrizing FACE2-like positions
520  else{
521 
522  OPENMP3DLOOPVARSDEFINE; OPENMP3DLOOPSETUPFULLP1EXCEPTX2;
523 #pragma omp for schedule(OPENMPSCHEDULE(),OPENMPCHUNKSIZE(blocksize))
524  OPENMP3DLOOPBLOCK{
525  OPENMP3DLOOPBLOCK2IJK(i,j,k);
526 
527 
528 #if(MCOORD==CARTMINKMETRIC)
529  // doesn't depend on position, only store/use 1 value
530  if(i!=0 || j!=0 || k!=0) continue; // simple way to avoid other i,j,k when doing OpenMP
531 #endif
532 
533  if(j>=N2/2){
534  DLOOP(jj,kk){
535  GLOBALMETMACP1A2(dxdxpstore,loc,i,j,k,jj,kk) = GLOBALMETMACP1A2(dxdxpstore,loc,i,N2-1-j,k,jj,kk);
536  GLOBALMETMACP1A2(idxdxpstore,loc,i,j,k,jj,kk) = GLOBALMETMACP1A2(idxdxpstore,loc,i,N2-1-j,k,jj,kk);
537  }
538 
539  DLOOP(jj,kk){
540  if(jj==2 && kk!=2 || jj!=2 && kk==2){
541  GLOBALMETMACP1A2(dxdxpstore,loc,i,j,k,jj,kk) *= -1.0;
542  GLOBALMETMACP1A2(idxdxpstore,loc,i,j,k,jj,kk) *= -1.0;
543  }
544  }
545 
546 
547  }// end over upper hemisphere
548  }// end 3D LOOP
549  }// end if over CENT-like w.r.t. FACE2
550 
551  }// end over locations
552  }// end parallel region
553 
554 
555  }// end if can symmetrize
556 
557 
558 
559 }
560 
561 
562 
563 
564 
565 
566 
567 
580 static void set_grid_metrics(void)
581 {
582 
583 
584 
585  // dualfprintf(fail_file,"Computing metric stuff\n");
586 
587 
588 #pragma omp parallel
589  {
590  int i, j, k, l, m;
591  int ii, jj, kk, ll;
592  int pl,pliter;
593  FTYPE X[NDIM];
594  FTYPE V[NDIM];
595  struct of_geom geomdontuse;
596  struct of_geom *ptrgeom=&geomdontuse;
597  struct of_geom geomdontusetest;
598  struct of_geom *ptrgeomtest=&geomdontusetest;
599  int loc;
600  extern void gcov_func(struct of_geom *ptrgeom, int getprim, int whichcoord, FTYPE *X, FTYPE *gcov, FTYPE *gcovpert);
601  extern void gcon_func(struct of_geom *ptrgeom, int getprim, int whichcoord, FTYPE *X, FTYPE *gcov, FTYPE *gcon);
602  extern void eomfunc_func(struct of_geom *ptrgeom, int getprim, int whichcoord, FTYPE *X, FTYPE *EOMFUNCNAME);
603  void assign_eomfunc(struct of_geom *geom, FTYPE *EOMFUNCNAME);
604  LOCALMETRICTEMPVARS;
605 
606 
607 
608 
610  OPENMP3DLOOPVARSDEFINE; OPENMP3DLOOPSETUPFULLP1;
611 #pragma omp for schedule(OPENMPSCHEDULE(),OPENMPCHUNKSIZE(blocksize))
612  OPENMP3DLOOPBLOCK{
613  OPENMP3DLOOPBLOCK2IJK(i,j,k);
614 
615 
616 #if(MCOORD==CARTMINKMETRIC)
617  // doesn't depend on position, only store/use 1 value
618  if(i!=0 || j!=0 || k!=0) continue; // simple way to avoid other i,j,k when doing OpenMP
619 #endif
620 
621  // over grid locations needing these quantities
622  for (loc = NPG - 1; loc >= 0; loc--) {
623 
624  bl_coord_ijk_2(i,j,k,loc,X, V);
626  //
627  // (1,MCOORD) here actually means PRIMCOORDS since the "1" means convert MCOORD to PRIMCOORDS.
628 
629  // fake geom just for indices
630  ptrgeom->i=i;
631  ptrgeom->j=j;
632  ptrgeom->k=k;
633  ptrgeom->p=loc;
634 
635  // get local metric quantities for this loc,i,j,k
636  GETLOCALMETRIC(loc,i,j,k);
637 
638 
639  // get metric terms and anything in struct of_geom
640  gcov_func(ptrgeom,1,MCOORD,X, localgcov,localgcovpert);
641  if(gdet_func_metric(MCOORD,V,localgcov,&(localgdet[0]))!=0){
642  if(debugfail>=2) dualfprintf(fail_file,"Caught gdet_func_metric() problem in set_grid.c:set_grid_metrics() i=%d j=%d k=%d\n",i,j,k);
643  }
644  gcon_func(ptrgeom,1,MCOORD,X,localgcov,localgcon);
645 #if(WHICHEOM!=WITHGDET)
646  // don't need (want) to recompute gdet that is just where localeomfunc points to if WHICHEOM==WITHGDET
647  eomfunc_func(ptrgeom,1,MCOORD,X,localeomfunc);
648 #endif
649  alphalapse_func(ptrgeom,1,MCOORD,X,localgcov,localgcon,localalphalapse);
650  betasqoalphasq_func(ptrgeom,1,MCOORD,X,localgcov,localgcon,localbetasqoalphasq);
651  beta_func(ptrgeom,1,MCOORD,X,localgcov,localgcon,*localalphalapse,localbeta);
652  if(GDETVOLDIFF){
653  // uses X,V (not det) from all locations
654  gdetvol_func(ptrgeom,localgdet,localeomfunc,localgdetvol);
655  }
656 
657 
658 #if(NEWMETRICSTORAGE)
659  GLOBALMETMACP1A0(compgeom,loc,i,j,k).i=ptrgeom->i;
660  GLOBALMETMACP1A0(compgeom,loc,i,j,k).j=ptrgeom->j;
661  GLOBALMETMACP1A0(compgeom,loc,i,j,k).k=ptrgeom->k;
662  GLOBALMETMACP1A0(compgeom,loc,i,j,k).p=ptrgeom->p;
663 
664 
665  // go ahead and assign eomfunc as per variable as either original eomfunc or gdet
666  // Must come *before* set_igdet_old() that sets 1/eomfunc
667  assign_eomfunc(&GLOBALMETMACP1A0(compgeom,loc,i,j,k), &(GLOBALMETMACP1A0(compgeom,loc,i,j,k).EOMFUNCMAC(0)));
668 
669  // then store extra things in global geom structure
670  // refer to "_old" version so really computs it
671  set_igdet_old(&GLOBALMETMACP1A0(compgeom,loc,i,j,k)); // full 1/gdet and 1/EOMFUNCMAC(pl)
672 
673 
674 
675  GLOBALMETMACP0A1(gdetgeom,i,j,k,loc).gdet=GLOBALMETMACP1A0(gdetgeomnormal,loc,i,j,k).gdet=GLOBALMETMACP1A0(compgeom,loc,i,j,k).gdet;
676 
677 
678  GLOBALMETMACP0A1(gdetgeom,i,j,k,loc).igdetnosing=GLOBALMETMACP1A0(gdetgeomnormal,loc,i,j,k).igdetnosing=GLOBALMETMACP1A0(compgeom,loc,i,j,k).igdetnosing;
679 
680 
681  PLOOP(pliter,pl){
682  GLOBALMETMACP0A1(gdetgeom,i,j,k,loc).EOMFUNCMAC(pl)=GLOBALMETMACP1A0(gdetgeomnormal,loc,i,j,k).EOMFUNCMAC(pl)=GLOBALMETMACP1A0(compgeom,loc,i,j,k).EOMFUNCMAC(pl);
683  GLOBALMETMACP0A1(gdetgeom,i,j,k,loc).IEOMFUNCNOSINGMAC(pl)=GLOBALMETMACP1A0(gdetgeomnormal,loc,i,j,k).IEOMFUNCNOSINGMAC(pl)=GLOBALMETMACP1A0(compgeom,loc,i,j,k).IEOMFUNCNOSINGMAC(pl);
684  }
685  // exit(0);
686 
687 #else
688  // GODMARK: then unless want to create independent storage area, computed each time needed
689 #endif
690 
691 
692 
694  //
695  // metric checks
696  //
698  get_geometry(i,j,k,loc,ptrgeomtest);
699  metric_checks(ptrgeomtest);
700 
701  }// end over location
702  }// end 3D LOOP
703  }// end parallel region
704 }
705 
706 
707 
708 
709 
711 static void set_grid_metrics_withsymmetrization(void)
712 {
713 
714 
715 
716  // dualfprintf(fail_file,"Computing metric stuff\n");
717 
718 
719  set_grid_metrics_gcov();
720 
721  if(ATTEMPTSYMMETRIZATION) symmetrize_gcov();
722 
723  set_grid_metrics_others();
724 
725 }
726 
727 
728 
729 
730 
731 
732 
734 static void set_grid_metrics_gcov(void)
735 {
736 
737 
738  // dualfprintf(fail_file,"Computing metric stuff\n");
739 
740 
741 #pragma omp parallel
742  {
743  int i, j, k, l, m;
744  int ii, jj, kk, ll;
745  int pl,pliter;
746  FTYPE X[NDIM];
747  FTYPE V[NDIM];
748  struct of_geom geomdontuse;
749  struct of_geom *ptrgeom=&geomdontuse;
750  int loc;
751  extern void gcov_func(struct of_geom *ptrgeom, int getprim, int whichcoord, FTYPE *X, FTYPE *gcov, FTYPE *gcovpert);
752  LOCALMETRICTEMPVARS;
753 
754 
755 
756 
758  OPENMP3DLOOPVARSDEFINE; OPENMP3DLOOPSETUPFULLP1;
759 #pragma omp for schedule(OPENMPSCHEDULE(),OPENMPCHUNKSIZE(blocksize))
760  OPENMP3DLOOPBLOCK{
761  OPENMP3DLOOPBLOCK2IJK(i,j,k);
762 
763 
764 #if(MCOORD==CARTMINKMETRIC)
765  // doesn't depend on position, only store/use 1 value
766  if(i!=0 || j!=0 || k!=0) continue; // simple way to avoid other i,j,k when doing OpenMP
767 #endif
768 
769  // over grid locations needing these quantities
770  for (loc = NPG - 1; loc >= 0; loc--) {
771 
772  bl_coord_ijk_2(i,j,k,loc,X, V);
774  //
775  // (1,MCOORD) here actually means PRIMCOORDS since the "1" means convert MCOORD to PRIMCOORDS.
776 
777  // fake geom just for indices
778  ptrgeom->i=i;
779  ptrgeom->j=j;
780  ptrgeom->k=k;
781  ptrgeom->p=loc;
782 
783  // get local metric quantities for this loc,i,j,k
784  GETLOCALMETRIC(loc,i,j,k);
785 
786  // get metric terms and anything in struct of_geom
787  gcov_func(ptrgeom,1,MCOORD,X, localgcov,localgcovpert);
788  if(gdet_func_metric(MCOORD,V,localgcov,&(localgdet[0]))!=0){
789  if(debugfail>=2) dualfprintf(fail_file,"Caught gdet_func_metric() problem in set_grid.c:set_grid_metrics_gcov() i=%d j=%d k=%d gcovtt=%21.15g gcovperttt=%21.15g\n",i,j,k,localgcov[GIND(TT,TT)],localgcovpert[TT]);
790  }
791  }// end over location
792  }// end 3D LOOP
793  }// end parallel region
794 
795 
796 }
797 
798 
799 
800 
801 
802 
803 
804 
805 static void set_grid_metrics_others(void)
806 {
807 
808 
809 
810 #pragma omp parallel
811  {
812  int i, j, k, l, m;
813  int ii, jj, kk, ll;
814  int pl,pliter;
815  FTYPE X[NDIM];
816  FTYPE V[NDIM];
817  struct of_geom geomdontuse;
818  struct of_geom *ptrgeom=&geomdontuse;
819  struct of_geom geomdontusetest;
820  struct of_geom *ptrgeomtest=&geomdontusetest;
821  int loc;
822  extern void gcov_func(struct of_geom *ptrgeom, int getprim, int whichcoord, FTYPE *X, FTYPE *gcov, FTYPE *gcovpert);
823  extern void gcon_func(struct of_geom *ptrgeom, int getprim, int whichcoord, FTYPE *X, FTYPE *gcov, FTYPE *gcon);
824  extern void eomfunc_func(struct of_geom *ptrgeom, int getprim, int whichcoord, FTYPE *X, FTYPE *EOMFUNCNAME);
825  void assign_eomfunc(struct of_geom *geom, FTYPE *EOMFUNCNAME);
826  LOCALMETRICTEMPVARS;
827 
828 
829 
830 
832  OPENMP3DLOOPVARSDEFINE; OPENMP3DLOOPSETUPFULLP1;
833 #pragma omp for schedule(OPENMPSCHEDULE(),OPENMPCHUNKSIZE(blocksize))
834  OPENMP3DLOOPBLOCK{
835  OPENMP3DLOOPBLOCK2IJK(i,j,k);
836 
837 
838 #if(MCOORD==CARTMINKMETRIC)
839  // doesn't depend on position, only store/use 1 value
840  if(i!=0 || j!=0 || k!=0) continue; // simple way to avoid other i,j,k when doing OpenMP
841 #endif
842 
843  // over grid locations needing these quantities
844  for (loc = NPG - 1; loc >= 0; loc--) {
845 
846  bl_coord_ijk_2(i,j,k,loc,X, V);
848  //
849  // (1,MCOORD) here actually means PRIMCOORDS since the "1" means convert MCOORD to PRIMCOORDS.
850 
851  // fake geom just for indices
852  ptrgeom->i=i;
853  ptrgeom->j=j;
854  ptrgeom->k=k;
855  ptrgeom->p=loc;
856 
857  // get local metric quantities for this loc,i,j,k
858  GETLOCALMETRIC(loc,i,j,k);
859 
860 
861  // get metric terms and anything in struct of_geom
862  gcon_func(ptrgeom,1,MCOORD,X,localgcov,localgcon);
863 #if(WHICHEOM!=WITHGDET)
864  eomfunc_func(ptrgeom,1,MCOORD,X,localeomfunc);
865 #endif
866  alphalapse_func(ptrgeom,1,MCOORD,X,localgcov,localgcon,localalphalapse);
867  betasqoalphasq_func(ptrgeom,1,MCOORD,X,localgcov,localgcon,localbetasqoalphasq);
868  beta_func(ptrgeom,1,MCOORD,X,localgcov,localgcon,*localalphalapse,localbeta);
869  if(GDETVOLDIFF){
870  // uses X,V (not det) from all locations
871  gdetvol_func(ptrgeom,localgdet,localeomfunc,localgdetvol);
872  }
873 
874 
875 #if(NEWMETRICSTORAGE)
876  GLOBALMETMACP1A0(compgeom,loc,i,j,k).i=ptrgeom->i;
877  GLOBALMETMACP1A0(compgeom,loc,i,j,k).j=ptrgeom->j;
878  GLOBALMETMACP1A0(compgeom,loc,i,j,k).k=ptrgeom->k;
879  GLOBALMETMACP1A0(compgeom,loc,i,j,k).p=ptrgeom->p;
880 
881 
882  // then store extra things in global geom structure
883  // refer to "_old" version so really computs it
884  set_igdet_old(&GLOBALMETMACP1A0(compgeom,loc,i,j,k)); // full 1/gdet and 1/EOMFUNCMAC(pl)
885 
886  // go ahead and assign eomfunc as per variable as either original eomfunc or gdet
887  assign_eomfunc(&GLOBALMETMACP1A0(compgeom,loc,i,j,k), &(GLOBALMETMACP1A0(compgeom,loc,i,j,k).EOMFUNCMAC(0)));
888 
889 
890  GLOBALMETMACP0A1(gdetgeom,i,j,k,loc).gdet=GLOBALMETMACP1A0(gdetgeomnormal,loc,i,j,k).gdet=GLOBALMETMACP1A0(compgeom,loc,i,j,k).gdet;
891 
892 
893  GLOBALMETMACP0A1(gdetgeom,i,j,k,loc).igdetnosing=GLOBALMETMACP1A0(gdetgeomnormal,loc,i,j,k).igdetnosing=GLOBALMETMACP1A0(compgeom,loc,i,j,k).igdetnosing;
894 
895 
896  PLOOP(pliter,pl){
897  GLOBALMETMACP0A1(gdetgeom,i,j,k,loc).EOMFUNCMAC(pl)=GLOBALMETMACP1A0(gdetgeomnormal,loc,i,j,k).EOMFUNCMAC(pl)=GLOBALMETMACP1A0(compgeom,loc,i,j,k).EOMFUNCMAC(pl);
898  GLOBALMETMACP0A1(gdetgeom,i,j,k,loc).IEOMFUNCNOSINGMAC(pl)=GLOBALMETMACP1A0(gdetgeomnormal,loc,i,j,k).IEOMFUNCNOSINGMAC(pl)=GLOBALMETMACP1A0(compgeom,loc,i,j,k).IEOMFUNCNOSINGMAC(pl);
899  }
900  // exit(0);
901 
902 #else
903  // GODMARK: then unless want to create independent storage area, computed each time needed
904 #endif
905 
906 
907 
909  //
910  // metric checks
911  //
913  get_geometry(i,j,k,loc,ptrgeomtest);
914  metric_checks(ptrgeomtest);
915 #if(0)
916  if(startpos[1]+i>totalsize[1]/2){
917 
918  //ptrgeomtest->gcon[GIND(PH,TT)]=0.0; // doesn't matter
919  ptrgeomtest->gcon[GIND(PH,RR)]=0.0; // matters
920  //ptrgeomtest->gcon[GIND(TT,RR)]=0.0; // doesn't matter
921  ptrgeomtest->gcon[GIND(PH,TH)]=0.0; // matters
922 
923  ptrgeomtest->gcov[GIND(PH,TT)]=0.0; // matters
924  ptrgeomtest->gcov[GIND(PH,RR)]=0.0; // matters
926 
927  // ptrgeomtest->beta[TT]=0.0; // doesn't matter
928  // ptrgeomtest->beta[RR]=0.0; // doesn't matter
929  // ptrgeomtest->beta[TH]=0.0; // doesn't matter
930  // ptrgeomtest->beta[PH]=0.0; // doesn't matter as long as gcon analytic!
931 
932  dualfprintf(fail_file,"i=%d j=%d beta[PH]=%21.15g\n",i,j,ptrgeomtest->beta[PH]);
933  }
934 #endif
935 
936  }// end over location
937  }// end 3D LOOP
938  }// end parallel region
939 
940 
941 
942 
943 }
944 
945 
946 
947 
948 
949 
952 static void symmetrize_gcov(void)
953 {
954 
955 
956 #if(NEWMETRICSTORAGE==1)
957 
958  // symmetrize gcov
959  if(numprocs==1 && ISSPCMCOORD(MCOORD) && BCtype[X2DN]==POLARAXIS && BCtype[X2UP]==POLARAXIS && N2>1){
960 
961 #pragma omp parallel
962  {
963  int i, j, k;
964  int jj,kk;
965  int loc;
966 
967 
968 
969  // over grid locations needing these quantities
970  for (loc = NPG - 1; loc >= 0; loc--) {
971 
972 
973  if(loc==FACE2 || loc==CORN1 || loc==CORN3 || loc==CORNT){
974 
976  OPENMP3DLOOPVARSDEFINE; OPENMP3DLOOPSETUPFULLP1;
977 #pragma omp for schedule(OPENMPSCHEDULE(),OPENMPCHUNKSIZE(blocksize))
978  OPENMP3DLOOPBLOCK{
979  OPENMP3DLOOPBLOCK2IJK(i,j,k);
980 
981 
982 #if(MCOORD==CARTMINKMETRIC)
983  // doesn't depend on position, only store/use 1 value
984  if(i!=0 || j!=0 || k!=0) continue; // simple way to avoid other i,j,k when doing OpenMP
985 #endif
986 
987 
988  if(j>=N2/2){
989  DLOOP(jj,kk) GLOBALMETMACP1A0(compgeom,loc,i,j,k).gcov[GIND(jj,kk)] = GLOBALMETMACP1A0(compgeom,loc,i,N2-j,k).gcov[GIND(jj,kk)];
990 
991  for(jj=0;jj<NDIM;jj++){
992  for(kk=0;kk<=jj;kk++){ // must avoid duplication
993  if(jj==2 && kk!=2 || jj!=2 && kk==2){
994  GLOBALMETMACP1A0(compgeom,loc,i,j,k).gcov[GIND(jj,kk)] *= -1.0;
995  }
996  }
997  }
998 
999  GLOBALMETMACP1A0(compgeom,loc,i,j,k).gdet = GLOBALMETMACP1A0(compgeom,loc,i,N2-j,k).gdet;
1000 
1001  }// end over upper hemisphere
1002 
1003  }// end 3D LOOP
1004  }// end if symmetrizing FACE2-like positions
1005  else{
1006 
1007  OPENMP3DLOOPVARSDEFINE; OPENMP3DLOOPSETUPFULLP1EXCEPTX2;
1008 #pragma omp for schedule(OPENMPSCHEDULE(),OPENMPCHUNKSIZE(blocksize))
1009  OPENMP3DLOOPBLOCK{
1010  OPENMP3DLOOPBLOCK2IJK(i,j,k);
1011 
1012 
1013 #if(MCOORD==CARTMINKMETRIC)
1014  // doesn't depend on position, only store/use 1 value
1015  if(i!=0 || j!=0 || k!=0) continue; // simple way to avoid other i,j,k when doing OpenMP
1016 #endif
1017 
1018 
1019  if(j>=N2/2){
1020  DLOOP(jj,kk) GLOBALMETMACP1A0(compgeom,loc,i,j,k).gcov[GIND(jj,kk)] = GLOBALMETMACP1A0(compgeom,loc,i,N2-1-j,k).gcov[GIND(jj,kk)];
1021 
1022  for(jj=0;jj<NDIM;jj++){
1023  for(kk=0;kk<=jj;kk++){ // must avoid duplication
1024  if(jj==2 && kk!=2 || jj!=2 && kk==2){
1025  GLOBALMETMACP1A0(compgeom,loc,i,j,k).gcov[GIND(jj,kk)] *= -1.0;
1026  }
1027  }
1028  }
1029 
1030  GLOBALMETMACP1A0(compgeom,loc,i,j,k).gdet = GLOBALMETMACP1A0(compgeom,loc,i,N2-1-j,k).gdet;
1031 
1032  }// end over upper hemisphere
1033  }// end 3D LOOP
1034  }// end if over CENT-like w.r.t. FACE2
1035 
1036  }// end over locations
1037  }// end parallel region
1038  }// end if can symmetrize
1039 
1040 #endif
1041 
1042 }
1043 
1044 
1051 static void set_connection(void)
1052 {
1053 
1054 
1055  // dualfprintf(fail_file,"Computing connection\n");
1056 
1057 
1058 #pragma omp parallel
1059  {
1060  int i, j, k;
1061  FTYPE X[NDIM];
1062  struct of_geom geomdontuse;
1063  struct of_geom *ptrgeom=&geomdontuse;
1064  int loc;
1065  int dim1, dim2, dim3;
1066 
1067 
1069  OPENMP3DLOOPVARSDEFINE; OPENMP3DLOOPSETUPFULL;
1070 #pragma omp for schedule(OPENMPSCHEDULE(),OPENMPCHUNKSIZE(blocksize))
1071  OPENMP3DLOOPBLOCK{
1072  OPENMP3DLOOPBLOCK2IJK(i,j,k);
1073 
1074 
1075 #if(MCOORD==CARTMINKMETRIC)
1076  // doesn't depend on position, only store/use 1 value
1077  if(i!=0 || j!=0 || k!=0) continue; // simple way to avoid other i,j,k when doing OpenMP
1078 #endif
1079 
1080 
1081 #if( CONNAXISYMM == 1 )
1082  //compute connection only for k = 0 (assumes it is axisymmetric)
1083  if( k != 0 ) continue;
1084 #endif
1085 
1086  loc=CENT;
1087  coord_ijk(i, j, k, loc, X);
1088  // in geom, only metric thing used is gcon to raise the lower connection.
1089  get_geometry(i, j, k, loc, ptrgeom);
1090 
1091  // dualfprintf(fail_file,"conn: i=%d j=%d k=%d\n",i,j,k);
1092  conn_func(MCOORD,X, ptrgeom, GLOBALMETMAC(conn,i,j,k),GLOBALMETMAC(conn2,i,j,k));
1093 
1094 
1095 
1096 
1097 
1098 
1099  // checks
1100  if(CONNMACHINEBODY){
1101  struct of_geom geomf1dontuse;
1102  struct of_geom *ptrgeomf1=&geomf1dontuse;
1103  struct of_geom geomf2dontuse;
1104  struct of_geom *ptrgeomf2=&geomf2dontuse;
1105 
1106  if(i!=N1+N1BND){
1107  // in geom, only metric thing used is gcon to raise the lower connection.
1108  get_geometry(i, j, k, FACE1, ptrgeomf1);
1109  get_geometry(ip1mac(i), j, k, FACE1, ptrgeomf2);
1110 
1111  FTYPE shouldbe = (ptrgeomf2->gdet - ptrgeomf1->gdet)/dx[1]/(ptrgeom->gdet);
1112  FTYPE isbe;
1113  int jj;
1114  isbe=0.0;
1115  DLOOPA(jj) isbe+=GLOBALMETMACP0A3(conn,i,j,k,jj,1,jj);
1116  FTYPE diffbe = fabs(shouldbe-isbe)/(fabs(shouldbe)+fabs(isbe)+SMALL);
1117 
1118  if(diffbe>NUMEPSILON*100){
1119  dualfprintf(fail_file,"Connection will lead to errors in constant pressure regions: dir=%d i=%d j=%d k=%d : %21.15g %21.15g : %21.15g\n",1,i,j,k,shouldbe,isbe,diffbe);
1120  }
1121  }
1122  if(j!=N2+N2BND){
1123  // in geom, only metric thing used is gcon to raise the lower connection.
1124  get_geometry(i, j, k, FACE2, ptrgeomf1);
1125  get_geometry(i, jp1mac(j), k, FACE2, ptrgeomf2);
1126 
1127  FTYPE shouldbe = (ptrgeomf2->gdet - ptrgeomf1->gdet)/dx[2]/(ptrgeom->gdet);
1128  FTYPE isbe;
1129  int jj;
1130  isbe=0.0;
1131  DLOOPA(jj) isbe+=GLOBALMETMACP0A3(conn,i,j,k,jj,2,jj);
1132  FTYPE diffbe = fabs(shouldbe-isbe)/(fabs(shouldbe)+fabs(isbe)+SMALL);
1133 
1134  if(diffbe>NUMEPSILON*100){
1135  dualfprintf(fail_file,"Connection will lead to errors in constant pressure regions: dir=%d i=%d j=%d k=%d : %21.15g %21.15g : %21.15g\n",2,i,j,k,shouldbe,isbe,diffbe);
1136  }
1137  }
1138  }
1139 
1140 
1141 
1142 
1143  }// end 3D LOOP
1144  }// end parallel region
1145 
1146 
1147 #if( CONNAXISYMM == 1 )
1148  //assuming connection axisymmetric, copy it from k = 0 cells to the rest of the cells (with k!=0)
1149 #pragma omp parallel
1150  {
1151  int i, j, k;
1152  FTYPE X[NDIM];
1153  struct of_geom geomdontuse;
1154  struct of_geom *ptrgeom=&geomdontuse;
1155  int loc;
1156  int dim1, dim2, dim3;
1157 
1158 
1160  OPENMP3DLOOPVARSDEFINE; OPENMP3DLOOPSETUPFULL;
1161 #pragma omp for schedule(OPENMPSCHEDULE(),OPENMPCHUNKSIZE(blocksize))
1162  OPENMP3DLOOPBLOCK{
1163  OPENMP3DLOOPBLOCK2IJK(i,j,k);
1164  if( k == 0 ) continue; //already computed connection for k = 0
1165  //now, assuming axisymmetry, reuse those values
1166  for( dim1 = 0; dim1 < NDIM; dim1++ ) {
1167  for( dim2 = 0; dim2 < NDIM; dim2++ ) {
1168  for( dim3 = 0; dim3 < NDIM; dim3++ ) {
1169  GLOBALMETMAC(conn,i,j,k)[dim1][dim2][dim3] = GLOBALMETMAC(conn,i,j,0)[dim1][dim2][dim3];
1170  }
1171  }
1172  }
1173  for( dim1 = 0; dim1 < NDIM; dim1++ ) {
1174  GLOBALMETMAC(conn2,i,j,k)[dim1] = GLOBALMETMAC(conn2,i,j,0)[dim1];
1175  }
1176  }// end 3D LOOP
1177  }// end parallel region
1178 #endif
1179 
1180  if(ATTEMPTSYMMETRIZATION) symmetrize_connection();
1181 
1182 }
1183 
1184 
1185 
1186 
1187 static void symmetrize_connection(void)
1188 {
1189 
1190 
1191 
1192  if(numprocs==1 && ISSPCMCOORD(MCOORD) && BCtype[X2DN]==POLARAXIS && BCtype[X2UP]==POLARAXIS && N2>1){
1193  // symmetrize
1194 #pragma omp parallel
1195  {
1196  int i, j, k;
1197  int jj,kk,pp;
1198 
1200  OPENMP3DLOOPVARSDEFINE; OPENMP3DLOOPSETUPFULL;
1201 #pragma omp for schedule(OPENMPSCHEDULE(),OPENMPCHUNKSIZE(blocksize))
1202  OPENMP3DLOOPBLOCK{
1203  OPENMP3DLOOPBLOCK2IJK(i,j,k);
1204 
1205 
1206 #if(MCOORD==CARTMINKMETRIC)
1207  // doesn't depend on position, only store/use 1 value
1208  if(i!=0 || j!=0 || k!=0) continue; // simple way to avoid other i,j,k when doing OpenMP
1209 #endif
1210 
1211 
1212  if(j>=N2/2){
1213  DLOOP(jj,kk) DLOOPA(pp) GLOBALMETMACP0A3(conn,i,j,k,jj,kk,pp) = GLOBALMETMACP0A3(conn,i,N2-1-j,k,jj,kk,pp); // N2-1-j because at loc=CENT
1214  DLOOP(jj,kk) DLOOPA(pp){
1215  if(jj==2 && kk!=2 && pp!=2 || jj!=2 && kk==2 && pp!=2 || jj!=2 && kk!=2 && pp==2 || jj==2 && kk==2 && pp==2){
1216  GLOBALMETMACP0A3(conn,i,j,k,jj,kk,pp) *= -1.0;
1217  }
1218  }
1219  }
1220 
1221  }// end 3D LOOP
1222  }// end parallel region
1223  } // end symmetrization
1224 
1225 
1226 }
1227 
1228 
1229 
1230 
1231 
1237 static void set_idxvol(void)
1238 {
1239 
1240 
1241 #pragma omp parallel
1242  {
1243  int i, j, k;
1244 
1246  OPENMP3DLOOPVARSDEFINE; OPENMP3DLOOPSETUPZLOOP;
1247 #pragma omp for schedule(OPENMPSCHEDULE(),OPENMPCHUNKSIZE(blocksize))
1248  OPENMP3DLOOPBLOCK{
1249  OPENMP3DLOOPBLOCK2IJK(i,j,k);
1250 
1251 
1252 #if(MCOORD==CARTMINKMETRIC)
1253  // doesn't depend on position, only store/use 1 value
1254  if(i!=0 || j!=0 || k!=0) continue; // simple way to avoid other i,j,k when doing OpenMP
1255 #endif
1256 
1257  // only at 1 location, centered, using surrounding edge values
1258  if((defcoord==LOGRUNITH)&&(MCOORD==KSCOORDS)){
1259  mks_unitheta_idxvol_func(i,j,k,GLOBALMETMAC(idxvol,i,j,k));
1260  }
1261  else{
1262  GLOBALMETMACP0A1(idxvol,i,j,k,TT)=1.0; // really 1/dt, but changes in time
1263  GLOBALMETMACP0A1(idxvol,i,j,k,RR)=1.0/dx[1];
1264  GLOBALMETMACP0A1(idxvol,i,j,k,TH)=1.0/dx[2];
1265  GLOBALMETMACP0A1(idxvol,i,j,k,PH)=1.0/dx[3];
1266  }
1267  }// end 3D LOOP
1268  }// end parallel region
1269 
1270 }
1271 
1272 
1273 
1274 
1280 static void set_tlab2ortho(void)
1281 {
1282  struct of_geom geomdontuse;
1283  struct of_geom *ptrgeom=&geomdontuse;
1284  int primcoord=1; // tells calc_ORTHOes() can use dxdxp to simplify metric before getting tetrad
1285 
1286 #pragma omp parallel
1287  {
1288  int i, j, k;
1289 
1291  OPENMP3DLOOPVARSDEFINE; OPENMP3DLOOPSETUPFULLP1;
1292 #pragma omp for schedule(OPENMPSCHEDULE(),OPENMPCHUNKSIZE(blocksize))
1293  OPENMP3DLOOPBLOCK{
1294  OPENMP3DLOOPBLOCK2IJK(i,j,k);
1295 
1296 
1297 #if(MCOORD==CARTMINKMETRIC)
1298  // doesn't depend on position, only store/use 1 value
1299  if(i!=0 || j!=0 || k!=0) continue; // simple way to avoid other i,j,k when doing OpenMP
1300 #endif
1301  int ll;
1302  for(ll=CENT;ll<CENT+BOOSTGRIDPOS;ll++){
1303  get_geometry(i, j, k, ll, ptrgeom);
1304  calc_ORTHOes(primcoord, ptrgeom, GLOBALMETMACP2A0(tlab2ortho,ll,LAB2ORTHO,i,j,k),GLOBALMETMACP2A0(tlab2ortho,ll,ORTHO2LAB,i,j,k));// pass [4][4] array
1305 
1306  // DEBUG:
1307  // int jj,kk;
1308  // dualfprintf(fail_file,"TEST ijkll=%d %d %d : %d\n",i,j,k,ll);
1309  // DLOOP(jj,kk) dualfprintf(fail_file,"jj=%d kk=%d : boostup=%g boostdown=%g\n",jj,kk,GLOBALMETMACP2A2(tlab2ortho,ll,LAB2ORTHO,i,j,k,jj,kk),GLOBALMETMACP2A2(tlab2ortho,ll,ORTHO2LAB,i,j,k,jj,kk));
1310 
1311 
1312  }
1313  }// end 3D LOOP
1314  }// end parallel region
1315 
1316 
1317 }
1318 
1319 
Generated on Fri May 20 2016 15:52:35 for HARM by   doxygen 1.8.3.1