HARM
harm and utilities
 All Data Structures Files Functions Variables Typedefs Macros Pages
phys.c
Go to the documentation of this file.
1 
8 #include "decs.h"
9 
10 
11 
15 //#define VARSTATIC static
16 #define VARSTATIC
17 
18 
20 #define MYII ptrgeom->i
21 #define MYJJ ptrgeom->j
22 #define MYKK ptrgeom->k
23 
24 
27 int get_stateforfluxcalc(int dimen, int isleftright, FTYPE *pr, struct of_geom *ptrgeom, struct of_state **qptr)
28 {
29  int pureget_stateforfluxcalc(FTYPE *pr, struct of_geom *ptrgeom, struct of_state *q);
30  struct of_state *qptrlocal;
31 
32 
33 #if(STOREFLUXSTATE)
34  // retrieve state from global array
35 
36 
37  // JCM: I see no point in copying the data, just should change pointer address
38  // pointer assign (overrides the existing default assignment of the pointer address in the calling function, so assumes calling function only uses pointer form subsequently)
39  if(isleftright==ISMIDDLE){
40  *qptr=&(GLOBALMAC(fluxstatecent,MYII,MYJJ,MYKK));
41  }
42  else{
43  *qptr=&(GLOBALMACP1A1(fluxstate,dimen,MYII,MYJJ,MYKK,isleftright));
44  }
45 
46 #else // else if not storing fluxstate
47  // compute state right now
48  pureget_stateforfluxcalc(pr, ptrgeom, *qptr);
49 #endif
50 
51  return (0);
52 }
53 
54 
57 int get_stateforsource(FTYPE *pr, struct of_geom *ptrgeom, struct of_state **qptr)
58 {
59  int pureget_stateforsource(FTYPE *pr, struct of_geom *ptrgeom, struct of_state *q);
60  // int i,j,k;
61 
62 
63 
64 #if(STOREFLUXSTATE)
65  // retrieve state from global array
66  // i=ptrgeom->i;
67  // j=ptrgeom->j;
68  // k=ptrgeom->k;
69 
70  *qptr=&(GLOBALMAC(fluxstatecent,MYII,MYJJ,MYKK)); // note that fluxstatecent for source has more than needed for source since was computed for flux+source and flux>source for all terms
71 #else
72  // compute state right now
73  pureget_stateforsource(pr, ptrgeom, *qptr);
74 #endif
75 
76 
77  return (0);
78 }
79 
82 int get_stateforinterpline(FTYPE *pr, struct of_geom *ptrgeom, struct of_state **qptr)
83 {
84  int pureget_stateforinterpline(FTYPE *pr, struct of_geom *ptrgeom, struct of_state *q);
85  // int i,j,k,jj;
86 
87 
88 
89 #if(STOREFLUXSTATE)
90  //if(ptrgeom->p==CENT) (GODMARK: assume interpline yrealin always at CENT, as it is now)
91 
92  // retrieve state from global array
93  // i=ptrgeom->i;
94  // j=ptrgeom->j;
95  // k=ptrgeom->k;
96 
97  *qptr=&(GLOBALMAC(fluxstatecent,MYII,MYJJ,MYKK)); // note that fluxstatecent for interpline has more than needed for interpline since was computed for flux+interpline and flux>interpline for all terms
98 #else
99  // compute state right now
100  pureget_stateforinterpline(pr, ptrgeom, *qptr);
101 #endif
102 
103 
104  return (0);
105 }
106 
107 
110 int get_stateforglobalwavespeeds(FTYPE *pr, struct of_geom *ptrgeom, struct of_state **qptr)
111 {
112  int pureget_stateforglobalwavespeeds(FTYPE *pr, struct of_geom *ptrgeom, struct of_state *q);
113  // int i,j,k,jj;
114 
115 
116 
117 #if(STOREFLUXSTATE)
118  //if(ptrgeom->p==CENT) (GODMARK: assume interpline yrealin always at CENT, as it is now)
119 
120  // retrieve state from global array
121  // i=ptrgeom->i;
122  // j=ptrgeom->j;
123  // k=ptrgeom->k;
124 
125  *qptr=&(GLOBALMAC(fluxstatecent,MYII,MYJJ,MYKK)); // note that fluxstatecent for globalwavespeds has more than needed for interpline since was computed for flux+(globalwavespeeds) and flux>interpline for all terms
126 #else
127  // compute state right now
128  pureget_stateforglobalwavespeeds(pr, ptrgeom, *qptr);
129 #endif
130 
131 
132  return (0);
133 }
134 
135 
136 
137 /* add in source terms related to geometry to equations of motion */
138 int source_conn(FTYPE *pr, struct of_geom *ptrgeom,
139  struct of_state *q,FTYPE *dU)
140 {
141  VARSTATIC int i, j, k, l;
142  VARSTATIC int pl,pliter;
143  VARSTATIC FTYPE dUconn[NPR];
144  VARSTATIC FTYPE todo[NPR];
145  VARSTATIC FTYPE ftemp;
146  VARSTATIC FTYPE mhd[NDIM][NDIM];
147  VARSTATIC FTYPE mhdrad[NDIM][NDIM];
148  VARSTATIC FTYPE flux[NDIM][NPR];
149  int primtofullflux(int returntype, FTYPE *pr, struct of_state *q, struct of_geom *ptrgeom, FTYPE (*flux)[NPR], FTYPE (*fluxabs)[NPR]);
150  void mhd_calc_0(FTYPE *pr, int dir, struct of_geom *geom, struct of_state *q, FTYPE *mhd, FTYPE *mhdabs);
151  void mhd_calc(FTYPE *pr, int dir, struct of_geom *geom, struct of_state *q, FTYPE *mhd, FTYPE *mhdabs);
152  VARSTATIC int myii,myjj,mykk;
153 
154 
155 
156  if((ANALYTICSOURCE)&&(defcoord==LOGRSINTH)&&(MCOORD==KSCOORDS)&&(EOMRADTYPE==EOMRADNONE)){
157  // have both WHICHEOM==WITHGDET and WHICHEOM==WITHNOGDET
158  // avoided if doing radiation since can't use mks_source_conn() for radiation without explicitly creating new function (can't just substitute q->uradcon/cov for q->ucon/cov)
159  mks_source_conn(pr,ptrgeom, q, dU); // returns without geometry prefactor
160  PLOOP(pliter,pl) dU[pl]*=ptrgeom->EOMFUNCMAC(pl);
161  }
162  else { // general non-analytic form, which uses an analytic or numerical origin for conn/conn2 (see set_grid.c)
163 
164 
166  //
167  // define first connection
168  //
169  // notice mhd_calc(), rather than primtoflux(), is used because this is a special connection only operated on the (at most) 4 energy-momentum EOMs.
170  //
171  // notice that mhd_calc_0 is used, which includes rest-mass since connection coefficient source term has rest-mass. The rest-mass was only subtracted out of conserved quantity and flux term.
172  // SUPERGODMARK: problem in non-rel gravity for below term
173  //DLOOPA(j) mhd_calc_0(pr, j, ptrgeom, q, mhd[j], NULL);
174  //
175 
176 
177  // Get MHD stress-energy tensor
178  DLOOPA(j) mhd_calc(pr, j, ptrgeom, q, mhd[j], NULL);
179 
180  // Get radiation stress-energy tensor (separate fluid from MHD fluid)
181  if(EOMRADTYPE!=EOMRADNONE){
182  DLOOPA(j) mhd_calc_rad(pr, j, ptrgeom, q, mhdrad[j], NULL);
183  }
184 
185 
186  /* contract mhd stress tensor with connection */
187  // mhd^{dir}_{comp} = mhd^j_k
188  // dU^{l} = T^j_k C^{k}_{l j}
189  PLOOP(pliter,pl) dUconn[pl]=0;
190 
191 
192 #if(MCOORD!=CARTMINKMETRIC)
193  myii=ptrgeom->i;
194  myjj=ptrgeom->j;
195  mykk=ptrgeom->k;
196 #else
197  myii=0;
198  myjj=0;
199  mykk=0;
200 #endif
201 
202 
203 
204 
205  // For MHD and radiation stress-energy tensor
206  // Group these since conn[] repeats so should be easier on memory accesses
207  if(EOMRADTYPE!=EOMRADNONE){
208  FTYPE connterm;
209  for(l=0;l<NDIM;l++) DLOOP(j,k){
210  connterm=GLOBALMETMACP0A3(conn,myii,myjj,mykk,k,l,j);
211  dUconn[UU+l] += mhd[j][k] * connterm;
212  dUconn[URAD0+l] += mhdrad[j][k] * connterm;
213  }
214  }
215  else{
216  // For MHD stress-energy tensor
217  for(l=0;l<NDIM;l++) DLOOP(j,k) dUconn[UU+l] += mhd[j][k] * GLOBALMETMACP0A3(conn,myii,myjj,mykk,k,l,j);
218  }
219 
220 
221 #if(REMOVERESTMASSFROMUU==2)
222  // then need to add-in density term to source (used to avoid non-rel problems)
223  // GODMARK: for no non-rel problems one must make sure conn^t_{lj} is computed accurately.
224  for(l=0;l<NDIM;l++) DLOOPA(j){
225  dUconn[UU+l] += - pr[RHO] * q->ucon[j] * GLOBALMETMACP0A3(conn,myii,myjj,mykk,TT,l,j);
226  }
227 #endif
228 
229 
231  //
232  // use first connection
233  //
234  // true as long as UU->U3 are not added/subtracted from other EOMs
235  // Note that UtoU is ignorant of this connection term and additions/subtractions of various EOMs to/from another must account for this connection here.
236  PLOOP(pliter,pl) dU[pl]+=ptrgeom->EOMFUNCMAC(pl)*dUconn[pl];
237 
238 
239 
240 
242  //
243  // second connection
244  //
245 
246 #if(WHICHEOM!=WITHGDET)
247 
248  // d_t(f T^t_\nu) = -d_x^j (f F^j_\nu) + (f F^j_\nu)(d_x^j(\ln(f/\detg))) + f T^\lambda_\kappa \Gamma^\kappa_{\nu\lambda}
249 
250  // conn2_j = (d_x^j(\ln(f/\detg)))
251 
252  // deal with second connection. Now can use general flux as defined by primtofullflux since all EOMs operate similarly with respect to second connection
253  primtofullflux(UEVOLVE,pr,q,ptrgeom,flux,NULL); // returns with geometry prefactor (f F^j_\nu)
254 
255 
256  // todo = whether that EOM has the NOGDET form. If so, then need 2nd connection. Do this instead of different connection2 for each EOM since each spatial component is actually the same.
257  PLOOP(pliter,pl) todo[pl] = (nogdetlist[pl]>0 ? 1 : 0);
258 
259  // conn2 is assumed to take care of sign
260  // conn2 has geom->e and normal d(ln(gdet)) factors combined
261 
262  if(REMOVERESTMASSFROMUU){
263  if(todo[RHO]!=todo[UU]){
264  dualfprintf(fail_file,"Mixed form of REMOVERESTMASSFROMUU and NOGDET's not allowed.\n");
265  myexit(1);
266  }
267  }
268 
270  //
271  // notice that we assume equations are differenced first, then one manipulates the connections.
272  // Thus, the manipulation of connections applies to final form of EOMs afer differencing or adding.
273  // This agrees with how code evolves conserved quantities, so simplest choice to make.
274  //
275  // Alternatively stated, UtoU() controls this ordering issue.
276 
277  PLOOP(pliter,pl) DLOOPA(j) dU[pl] += todo[pl]*(flux[j][pl] * GLOBALMETMACP0A1(conn2,myii,myjj,mykk,j)); // (f F^j_\nu) conn2_j
278 
279 
280 
281 #endif // end if (WHICHEOM!=WITHGDET)
282 
283  }
284 
285 
286 
287  /* done! */
288  return (0);
289 }
290 
291 
292 
295 #if(NEWMETRICSTORAGE==0)
296 void get_geometry_old(int ii, int jj, int kk, int pp, struct of_geom *geom)
297 {
298  VARSTATIC int j, k;
299  void assign_eomfunc(struct of_geom *geom, FTYPE *EOMFUNCNAME);
300  VARSTATIC int myii,myjj,mykk;
301  VARSTATIC int pl,pliter;
302 
303 
304 #if(MCOORD!=CARTMINKMETRIC)
305  myii=ii;
306  myjj=jj;
307  mykk=kk;
308 #else
309  myii=0;
310  myjj=0;
311  mykk=0;
312 #endif
313 
314 
315 
316 #if(GETGEOMUSEPOINTER==0)
317  // DLOOP(j,k) geom->gcov[GIND(j,k)] = GLOBALMETMACP1A2(gcov,pp,myii,myjj,mykk,j,k);
318  //DLOOP(j,k) geom->gcon[GIND(j,k)] = GLOBALMETMACP1A2(gcon,pp,myii,myjj,mykk,j,k);
319  // let's vectorize it
320  for(j=0;j<=NDIM*NDIM-1;j++){
321  geom->gcon[GIND(0,j)] = GLOBALMETMACP1A2(gcon,pp,myii,myjj,mykk,0,j);
322  geom->gcov[GIND(0,j)] = GLOBALMETMACP1A2(gcov,pp,myii,myjj,mykk,0,j);
323  }
324  for(j=0;j<NDIM;j++){
325  geom->gcovpert[j]=GLOBALMETMACP1A1(gcovpert,pp,myii,myjj,mykk,j);
326  }
327 #else
328  // let's pointer it, faster by a bit
329  geom->gcov=(FTYPE (*)[NDIM])(&(GLOBALMETMACP1A2(gcov,pp,myii,myjj,mykk,0,0))); // pointer
330 
331  geom->gcovpert=(FTYPE (*))(&(GLOBALMETMACP1A0(gcovpert,pp,myii,myjj,mykk))); // pointer
332 
333  geom->gcon=(FTYPE (*)[NDIM])(&(GLOBALMETMACP1A2(gcon,pp,myii,myjj,mykk,0,0))); // pointer
334 #endif
335 
336  geom->gdet = GLOBALMETMACP1A0(gdet,pp,myii,myjj,mykk);
337 
338  // get eomfunc (see eomfunc_func() and lngdet_func_mcoord() in metric.c)
339  assign_eomfunc(geom,GLOBALMETMACP1A0(eomfunc,pp,myii,myjj,mykk));
340 
341  geom->alphalapse = GLOBALMETMACP1A0(alphalapse,pp,myii,myjj,mykk);
342 
344  //
345  // these should be real grid i,j,k positions
346  geom->i = ii;
347  geom->j = jj;
348  geom->k = kk;
349  geom->p = pp;
350 }
351 
352 // load local geometry into structure geom
354 void get_geometry_gdetonly_old(int ii, int jj, int kk, int pp, struct of_geom *geom)
355 {
356  VARSTATIC int j, k;
357  void assign_eomfunc(struct of_geom *geom, FTYPE *EOMFUNCNAME);
358  VARSTATIC int myii,myjj,mykk;
359  VARSTATIC int pl,pliter;
360 
361 
362 #if(MCOORD!=CARTMINKMETRIC)
363  myii=ii;
364  myjj=jj;
365  mykk=kk;
366 #else
367  myii=0;
368  myjj=0;
369  mykk=0;
370 #endif
371 
372 
373 
374  geom->gdet = GLOBALMETMACP1A0(gdet,pp,myii,myjj,mykk);
375 
376 
378  //
379  // these should be real grid i,j,k positions
380  geom->i = ii;
381  geom->j = jj;
382  geom->k = kk;
383  geom->p = pp;
384 
385 
386 }
387 
388 
391 void get_geometry_geomeonly_old(int ii, int jj, int kk, int pp, struct of_geom *geom)
392 {
393  VARSTATIC int pl,pliter;
394  // void get_geometry(int ii, int jj, int kk, int pp, struct of_geom *geom);
395  // void get_geometry_gdetonly(int ii, int jj, int kk, int pp, struct of_geom *geom);
396 
397 
398 
399 #if(WHICHEOM!=WITHGDET)
400  // more complicated ptrgeom->e
401  get_geometry(ii, jj, kk, pp, geom);
402 
403 #else
404  // then ptrgeom->e is just ptrgeom->g
405  get_geometry_gdetonly(ii, jj, kk, pp, geom);
406  PALLLOOP(pl) geom->EOMFUNCMAC(pl)=geom->gdet;
407 #endif
408 
409 
410 }
411 
412 
413 #endif // end if old metric method
414 
415 
416 
417 
418 
419 
420 
421 /* load local geometry into structure geom */
422 void get_allgeometry(int ii, int jj, int kk, int pp, struct of_allgeom *allgeom, struct of_geom *ptrgeom)
423 {
424  VARSTATIC int j, k;
425  VARSTATIC int pl,pliter;
426 
427 
428  get_geometry(ii, jj, kk, pp, ptrgeom);
429 
430  allgeom->i=ptrgeom->i;
431  allgeom->j=ptrgeom->j;
432  allgeom->k=ptrgeom->k;
433  allgeom->p=ptrgeom->p;
434 
435 #if(GETGEOMUSEPOINTER==0 || NEWMETRICSTORAGE==1)
436  for(j=0;j<=NDIM*NDIM-1;j++){
437  allgeom->gcov[GIND(0,j)]=ptrgeom->gcov[GIND(0,j)];
438  allgeom->gcon[GIND(0,j)]=ptrgeom->gcon[GIND(0,j)];
439  }
440  for(j=0;j<NDIM;j++){
441  allgeom->gcovpert[j]=ptrgeom->gcovpert[j];
442  }
443 #else
444  allgeom->gcov=ptrgeom->gcov;
445  allgeom->gcon=ptrgeom->gcon;
446  allgeom->gcovpert=ptrgeom->gcovpert;
447 #endif
448 
449  allgeom->gdet=ptrgeom->gdet;
450  PLOOP(pliter,pl) allgeom->EOMFUNCMAC(pl)=ptrgeom->EOMFUNCMAC(pl);
451  allgeom->alphalapse = ptrgeom->alphalapse;
452 
453  PALLLOOP(pl) allgeom->IEOMFUNCNOSINGMAC(pl) = ptrgeom->IEOMFUNCNOSINGMAC(pl);
454  allgeom->igdetnosing = ptrgeom->igdetnosing;
455 
456 
457  coord_ijk(ii, jj, kk, pp, allgeom->X);
458  bl_coord_ijk( ii, jj, kk, pp , allgeom->V );
459  dxdxprim_ijk( ii, jj, kk, pp , allgeom->dxdxp);
460 
461 
462 }
463 
464 
467 void assign_eomfunc(struct of_geom *geom, FTYPE *EOMFUNCNAME)
468 {
469  int pliter,pl;
470 
471 #if(WHICHEOM==WITHGDET)
472  return; // now nothing to do since using new EOMFUNCMAC stuff
473 #endif
474 
475  // now set each EOM
476  PLOOP(pliter,pl){
477  if(nogdetlist[pl]==0) geom->EOMFUNCMAC(pl)=geom->gdet;
478  else geom->EOMFUNCMAC(pl)=EOMFUNCASSIGN(pl);
479  }
480 
481  if(NEWMETRICSTORAGE){
482  // assign inverses (should already be set if PRIMECOORDS)
483  // this only changes (overrides) to use igdetnosing since otherwise already should be set
484 
485  // now set each EOM
486  PLOOP(pliter,pl){
487  if(nogdetlist[pl]==0) geom->IEOMFUNCNOSINGMAC(pl)=geom->igdetnosing;
488  // else will use igdet
489  }
490  } // end if NEWMETRICSTORAGE==1
491 
492 }
493 
494 
495 
496 /*
497 Kraken stopped when generated core dump:
498 
499 Core was generated by `/lustre/scratch/jmckinne/rad1b//grmhd 32 8 4 1 0'.
500 Program terminated with signal 11, Segmentation fault.
501 #0 0x00000000005359fb in current_doprecalc (which=-131040, p=0xfffffffffffffffc) at phys.c:549
502 549 current_precalc(which,ptrgeom,&q,Dt,GLOBALMAC(cfaraday,i,j,k));
503 (gdb) bt
504 #0 0x00000000005359fb in current_doprecalc (which=-131040, p=0xfffffffffffffffc) at phys.c:549
505 #1 0x000000000054d57d in step_ch_full (truestep=-131040, prim=0xfffffffffffffffc, pstag=0xfffffffffffffffc, ucons=0xfffffffffffffffc, vpot=0xffffffffffffff40, Bhat=0xffffffffffff3400, pl_ct=0x302ca00, pr_ct=0x1415d20, F1=0x1f2f968, F2=0x2e4cd68, F3=0xb84cc8, Atemp=0x196dca0,
506  uconstemp=0xa60c60) at step_ch.c:51
507 #2 0x00000000005059a6 in main (argc=32700776, argv=0x302ca00) at main.c:115
508 
509 */
510 
511 
513 int current_doprecalc(int which, FTYPE (*p)[NSTORE2][NSTORE3][NPR])
514 {
515 
516 
517 #if(CALCFARADAYANDCURRENTS==0)
518  dualfprintf(fail_file,"Shouldn't be here in current_doprecalc()\n");
519  myexit(12);
520 #endif
521 
522 
523 
524 
525 #pragma omp parallel OPENMPGLOBALPRIVATEFORSTATEANDGEOM // requires full state (OPTMARK: But shouldn't have to since just need field part of state)
526  {
527  int i,j,k;
528  int idel, jdel,kdel;
529  struct of_geom geomdontuse;
530  struct of_geom *ptrgeom=&geomdontuse;
531  struct of_state q;
532  FTYPE Dt;
533  int face;
534  OPENMP3DLOOPVARSDEFINE;
535 
536 
537  if(WHICHCURRENTCALC==CURRENTCALC0){ // then need to save 2 times
538  // which==1,2 should be using face primitives, but probably not
539  if(which==CURTYPET){ face=CENT; idel=0; jdel=0; kdel=0; Dt=dt*0.5; }
540  else if(which==CURTYPEX){ face=FACE1; idel=1; jdel=0; kdel=0; Dt=dt; }
541  else if(which==CURTYPEY){ face=FACE2; idel=0; jdel=1; kdel=0; Dt=dt; }
542  else if(which==CURTYPEZ){ face=FACE3; idel=0; jdel=0; kdel=1; Dt=dt; }
543  else if(which==CURTYPEFARADAY){ face=CENT; idel=0; jdel=0; kdel=0; Dt=dt; }
544  }
545  else if(WHICHCURRENTCALC==CURRENTCALC1){ // time and space centered on present time (best method)
546  face=CENT;
547  idel=0;
548  jdel=0;
549  kdel=0;
550  Dt=dt;
551  }
552  else if(WHICHCURRENTCALC==CURRENTCALC2){ // then need to save 2 times
553  if(which==CURTYPET){ face=CENT; idel=0; jdel=0; kdel=0; Dt=dt*0.5; }
554  else if(which==CURTYPEX){ face=CENT; idel=1; jdel=0; kdel=0; Dt=dt; }
555  else if(which==CURTYPEY){ face=CENT; idel=0; jdel=1; kdel=0; Dt=dt; }
556  else if(which==CURTYPEZ){ face=CENT; idel=0; jdel=0; kdel=1; Dt=dt; }
557  else if(which==CURTYPEFARADAY){ face=CENT; idel=0; jdel=0; kdel=0; Dt=dt; }
558  }
559  else{
560  dualfprintf(fail_file,"No such WHICHCURRENTCALC=%d\n",WHICHCURRENTCALC);
561  }
562 
563  // assume no other conditions GODMARK
564 
565  // COMPFULLLOOP{
566  OPENMP3DLOOPSETUPFULL;
567 #pragma omp for schedule(OPENMPSCHEDULE(),OPENMPCHUNKSIZE(blocksize))
568  OPENMP3DLOOPBLOCK{
569  OPENMP3DLOOPBLOCK2IJK(i,j,k);
570 
571  get_geometry(i, j, k, face, ptrgeom);
572  MYFUN(get_state(MAC(p,i,j,k), ptrgeom, &q),"phys.c:current_doprecalc()", "get_state()", 1);
573  current_precalc(which,ptrgeom,&q,Dt,GLOBALMAC(cfaraday,i,j,k));
574  }
575  }// end parallel region
576 
577 
578 
579  return(0);
580 }
581 
582 
587 void current_precalc(int which, struct of_geom *geom, struct of_state *q, FTYPE Dt,FTYPE (*faraday)[3])
588 {
589  // assume outside loop is like flux, from 0..N in r for r, 0..N in h for h. And must wait till 2nd timestep before computing the current since need time differences
590 
591 #if(CALCFARADAYANDCURRENTS==0)
592  dualfprintf(fail_file,"Shouldn't be here in current_doprecalc()\n");
593  myexit(12);
594 #endif
595 
596 
597  if(which==CURTYPET){ // for d/dt
598 
599  // assume got here when DT==dt/2 and geom and state set at zone center
600  // first save old calculation
601  if((WHICHCURRENTCALC==CURRENTCALC0)||(WHICHCURRENTCALC==CURRENTCALC2)){ // then need to save 2 times
602  faraday[4][0]=faraday[0][0];
603  faraday[4][1]=faraday[0][1];
604  faraday[4][2]=faraday[0][2];
605  }
606  else if(WHICHCURRENTCALC==CURRENTCALC1){ // then need to save 3 times and have [4] as 2 times ago
607  // 2 times ago
608  faraday[4][0]=faraday[5][0];
609  faraday[4][1]=faraday[5][1];
610  faraday[4][2]=faraday[5][2];
611 
612  // 1 time ago
613  faraday[5][0]=faraday[0][0];
614  faraday[5][1]=faraday[0][1];
615  faraday[5][2]=faraday[0][2];
616  }
617  // now calculate new version
618 
619  // Need F^{xt} F^{yt} F^{zt}
620  faraday[0][0]=-1.0/geom->gdet * (-q->bcov[3]*q->ucov[2]+q->bcov[2]*q->ucov[3]); // f^{rt}
621  faraday[0][1]=-1.0/geom->gdet * (q->bcov[3]*q->ucov[1]-q->bcov[1]*q->ucov[3]); // f^{ht}
622  faraday[0][2]=-1.0/geom->gdet * (-q->bcov[2]*q->ucov[1]+q->bcov[1]*q->ucov[2]); // f^{pt}
623  }
624  else if(which==CURTYPEX){// for d/dx
625 
626  // Need F^{tx} F^{yx} F^{zx}
627 
628  // assume got here with DT=dt and geom and state at radial zone edge
629  faraday[1][0]=-1.0/geom->gdet * (q->bcov[3]*q->ucov[2]-q->bcov[2]*q->ucov[3]); // f^{tr}=-f^{rt}
630  faraday[1][1]=-1.0/geom->gdet * (-q->bcov[3]*q->ucov[0]+q->bcov[0]*q->ucov[3]); // f^{hr}
631  faraday[1][2]=-1.0/geom->gdet * (q->bcov[2]*q->ucov[0]-q->bcov[0]*q->ucov[2]); // f^{pr}
632  }
633  else if(which==CURTYPEY){ // for d/dy
634 
635  // Need F^{ty} F^{xy} F^{zy}
636 
637  // assume got here with DT=dt and geom and state at theta zone edge
638  faraday[2][0]=-1.0/geom->gdet * (-q->bcov[3]*q->ucov[1]+q->bcov[1]*q->ucov[3]); // f^{th}=-f^{ht}
639  faraday[2][1]=-1.0/geom->gdet * (q->bcov[3]*q->ucov[0]-q->bcov[0]*q->ucov[3]); // f^{rh}=-f^{hr}
640  faraday[2][2]=-1.0/geom->gdet * (-q->bcov[1]*q->ucov[0]+q->bcov[0]*q->ucov[1]); // f^{ph}
641  }
642  else if(which==CURTYPEZ){ // for d/dz
643 
644  // Need F^{tz} F^{xz} F^{yz}
645 
646  faraday[3][0]=1.0/geom->gdet * (-q->bcov[3]*q->ucov[2]+q->bcov[2]*q->ucov[3]); // f^{rt} // f^{tr}=-f^{rt}
647  faraday[3][1]=1.0/geom->gdet * (q->bcov[2]*q->ucov[0]-q->bcov[0]*q->ucov[2]); // f^{rp}=-f^{pr}
648  faraday[3][2]=1.0/geom->gdet * (-q->bcov[1]*q->ucov[0]+q->bcov[0]*q->ucov[1]); // f^{hp}= -f^{ph}
649 
650  }
651  else if(which==CURTYPEFARADAY){
652  // DT==dt, but zone center
653  GLOBALMACP0A1(fcon,geom->i,geom->j,geom->k,0)=-1.0/geom->gdet * (q->bcov[3]*q->ucov[2]-q->bcov[2]*q->ucov[3]); // f^{tr}
654  GLOBALMACP0A1(fcon,geom->i,geom->j,geom->k,1)=-1.0/geom->gdet * (-q->bcov[3]*q->ucov[1]+q->bcov[1]*q->ucov[3]); // f^{th}
655  GLOBALMACP0A1(fcon,geom->i,geom->j,geom->k,2)=-1.0/geom->gdet * (q->bcov[2]*q->ucov[1]-q->bcov[1]*q->ucov[2]); // f^{tp}
656  GLOBALMACP0A1(fcon,geom->i,geom->j,geom->k,3)=-1.0/geom->gdet * (q->bcov[3]*q->ucov[0]-q->bcov[0]*q->ucov[3]); // f^{rh}
657  GLOBALMACP0A1(fcon,geom->i,geom->j,geom->k,4)=-1.0/geom->gdet * (-q->bcov[2]*q->ucov[0]+q->bcov[0]*q->ucov[2]); // f^{rp}
658  GLOBALMACP0A1(fcon,geom->i,geom->j,geom->k,5)=-1.0/geom->gdet * (q->bcov[1]*q->ucov[0]-q->bcov[0]*q->ucov[1]); // f^{hp}
659  }
660  else{
661  dualfprintf(fail_file,"No such which in current_doprecalc()\n");
662  }
663 }
664 
665 
667 void current_calc(FTYPE (*cfaraday)[NSTORE2][NSTORE3][NUMCURRENTSLOTS][3])
668 {
669  void current_calc_0(FTYPE (*cfaraday)[NSTORE2][NSTORE3][NUMCURRENTSLOTS][3]);
670  void current_calc_1(int which, FTYPE (*cfaraday)[NSTORE2][NSTORE3][NUMCURRENTSLOTS][3]);
671 
672 
673 #if(CALCFARADAYANDCURRENTS==0)
674  dualfprintf(fail_file,"Shouldn't be here in current_doprecalc()\n");
675  myexit(12);
676 #endif
677 
678  if(WHICHCURRENTCALC==CURRENTCALC0){
679  current_calc_0(cfaraday);
680  }
681  else if((WHICHCURRENTCALC==CURRENTCALC1)||(WHICHCURRENTCALC==CURRENTCALC2)){
682  current_calc_1(WHICHCURRENTCALC,cfaraday);
683  }
684 
685 }
686 
687 
688 
692 void current_calc_0(FTYPE (*cfaraday)[NSTORE2][NSTORE3][NUMCURRENTSLOTS][3])
693 {
694  static FTYPE lastdt;
695  static int calls=0;
696 
697 
698 
699 #if(CALCFARADAYANDCURRENTS==0)
700  dualfprintf(fail_file,"Shouldn't be here in current_doprecalc()\n");
701  myexit(12);
702 #endif
703 
704 
705 #pragma omp parallel
706  {
707  int i,j,k;
708  struct of_geom geomtdontuse;
709  struct of_geom *ptrgeomt=&geomtdontuse;
710 
711  struct of_geom geomrdontuse;
712  struct of_geom *ptrgeomr=&geomrdontuse;
713 
714  struct of_geom geomhdontuse;
715  struct of_geom *ptrgeomh=&geomhdontuse;
716 
717  struct of_geom geompdontuse;
718  struct of_geom *ptrgeomp=&geompdontuse;
719 
720  struct of_geom geomrp1dontuse;
721  struct of_geom *ptrgeomrp1=&geomrp1dontuse;
722 
723  struct of_geom geomhp1dontuse;
724  struct of_geom *ptrgeomhp1=&geomhp1dontuse;
725 
726  struct of_geom geompp1dontuse;
727  struct of_geom *ptrgeompp1=&geompp1dontuse;
728 
729  FTYPE idtc,idx1,idx2,idx3;
730  FTYPE timeterm[NDIM];
731  OPENMP3DLOOPVARSDEFINE;
732 
733 
734  idtc=2.0/(lastdt+dt);
735  idx1=1.0/dx[1];
736  idx2=1.0/dx[2];
737  idx3=1.0/dx[3];
738 
740  OPENMP3DLOOPSETUP(INP11,OUTP11,INP12,OUTP12,INP13,OUTP13);
741 #pragma omp for schedule(OPENMPSCHEDULE(),OPENMPCHUNKSIZE(blocksize))
742  OPENMP3DLOOPBLOCK{
743  OPENMP3DLOOPBLOCK2IJK(i,j,k);
744 
745 
746  get_geometry(i,j,k,CENT,ptrgeomt);
747  get_geometry(i,j,k,FACE1,ptrgeomr);
748  get_geometry(i,j,k,FACE2,ptrgeomh);
749  get_geometry(i,j,k,FACE3,ptrgeomp);
750  // geomtp1 is same as geomt since d/dt( geometry) -> 0
751  get_geometry(ip1mac(i),j,k,FACE1,ptrgeomrp1);
752  get_geometry(i,jp1mac(j),k,FACE2,ptrgeomhp1);
753  get_geometry(i,j,kp1mac(k),FACE3,ptrgeompp1);
754 
755  if(calls>0){ // since need 2 times
756  timeterm[0]=0;
757  timeterm[1]=+(GLOBALMACP0A2(cfaraday,i,j,k,0,0)-GLOBALMACP0A2(cfaraday,i,j,k,4,0))*idtc; // F^{rt},t
758  timeterm[2]=+(GLOBALMACP0A2(cfaraday,i,j,k,0,1)-GLOBALMACP0A2(cfaraday,i,j,k,4,1))*idtc; // F^{ht},t
759  timeterm[3]=+(GLOBALMACP0A2(cfaraday,i,j,k,0,2)-GLOBALMACP0A2(cfaraday,i,j,k,4,2))*idtc; // F^{pt},t
760  }
761  else{
762  timeterm[0]=0;
763  timeterm[1]=0;
764  timeterm[2]=0;
765  timeterm[3]=0;
766  }
767 
768  // J^t = F^{tr},r + F^{th},h + F^{tp},p
769  GLOBALMACP0A1(jcon,i,j,k,0)=
770  +1./ptrgeomt->gdet*(ptrgeomrp1->gdet*GLOBALMACP0A2(cfaraday,ip1mac(i),j,k,1,0)-ptrgeomr->gdet*GLOBALMACP0A2(cfaraday,i,j,k,1,0))*idx1 // F^{tr},r
771  +1./ptrgeomt->gdet*(ptrgeomhp1->gdet*GLOBALMACP0A2(cfaraday,i,jp1mac(j),k,2,0)-ptrgeomh->gdet*GLOBALMACP0A2(cfaraday,i,j,k,2,0))*idx2 // F^{th},h
772  +1./ptrgeomt->gdet*(ptrgeompp1->gdet*GLOBALMACP0A2(cfaraday,i,j,kp1mac(k),3,0)-ptrgeomp->gdet*GLOBALMACP0A2(cfaraday,i,j,k,3,0))*idx3 // F^{tp},p
773  ;
774 
775  // J^r = F^{rt},t + F^{rh},h + F^{rp},p
776  GLOBALMACP0A1(jcon,i,j,k,1)=
777  +timeterm[1]
778  +1./ptrgeomt->gdet*(ptrgeomhp1->gdet*GLOBALMACP0A2(cfaraday,i,jp1mac(j),k,2,1)-ptrgeomh->gdet*GLOBALMACP0A2(cfaraday,i,j,k,2,1))*idx2 // F^{rh},h
779  +1./ptrgeomt->gdet*(ptrgeompp1->gdet*GLOBALMACP0A2(cfaraday,i,j,kp1mac(k),3,1)-ptrgeomp->gdet*GLOBALMACP0A2(cfaraday,i,j,k,3,1))*idx3 // F^{rp},p
780  ;
781 
782  // J^h = F^{ht},t + F^{hr},r + F^{hp},p
783  GLOBALMACP0A1(jcon,i,j,k,2)=
784  +timeterm[2]
785  +1./ptrgeomt->gdet*(ptrgeomrp1->gdet*GLOBALMACP0A2(cfaraday,ip1mac(i),j,k,1,1)-ptrgeomr->gdet*GLOBALMACP0A2(cfaraday,i,j,k,1,1))*idx1 // F^{hr},r
786  +1./ptrgeomt->gdet*(ptrgeompp1->gdet*GLOBALMACP0A2(cfaraday,i,j,kp1mac(k),3,2)-ptrgeomp->gdet*GLOBALMACP0A2(cfaraday,i,j,k,3,2))*idx3 // F^{hp},p
787  ;
788 
789  // J^p = F^{pt},t + F^{pr},r + F^{ph},h
790  GLOBALMACP0A1(jcon,i,j,k,3)=
791  +timeterm[3]
792  +1./ptrgeomt->gdet*(ptrgeomrp1->gdet*GLOBALMACP0A2(cfaraday,ip1mac(i),j,k,1,2)-ptrgeomr->gdet*GLOBALMACP0A2(cfaraday,i,j,k,1,2))*idx1 // F^{pr},r
793  +1./ptrgeomt->gdet*(ptrgeomhp1->gdet*GLOBALMACP0A2(cfaraday,i,jp1mac(j),k,2,2)-ptrgeomh->gdet*GLOBALMACP0A2(cfaraday,i,j,k,2,2))*idx2 // F^{ph},h
794  ;
795 
796  }// end loops
797 
798 
799  }// end over parallel region
800 
801 
802  // static (should be outside parallel region)
803  calls++;
804  lastdt=dt;
805 
806 }
807 
808 
809 
814 void current_calc_1(int which, FTYPE (*cfaraday)[NSTORE2][NSTORE3][NUMCURRENTSLOTS][3])
815 {
816  static FTYPE lastdt;
817  static long long calls=0;
818 
819 
820 #if(CALCFARADAYANDCURRENTS==0)
821  dualfprintf(fail_file,"Shouldn't be here in current_doprecalc()\n");
822  myexit(12);
823 #endif
824 
825 
826 
827 #pragma omp parallel
828  {
829  int i,j,k;
830  struct of_geom geomtdontuse;
831  struct of_geom *ptrgeomt=&geomtdontuse;
832 
833  struct of_geom geomcdontuse;
834  struct of_geom *ptrgeomc=&geomcdontuse;
835 
836  struct of_geom geomrp1dontuse;
837  struct of_geom *ptrgeomrp1=&geomrp1dontuse;
838 
839  struct of_geom geomhp1dontuse;
840  struct of_geom *ptrgeomhp1=&geomhp1dontuse;
841 
842  struct of_geom geompp1dontuse;
843  struct of_geom *ptrgeompp1=&geompp1dontuse;
844 
845  struct of_geom geomrm1dontuse;
846  struct of_geom *ptrgeomrm1=&geomrm1dontuse;
847 
848  struct of_geom geomhm1dontuse;
849  struct of_geom *ptrgeomhm1=&geomhm1dontuse;
850 
851  struct of_geom geompm1dontuse;
852  struct of_geom *ptrgeompm1=&geompm1dontuse;
853 
854  FTYPE idtc,idx1,idx2,idx3;
855  FTYPE AA,BB;
856  FTYPE dtl,dtr;
857  FTYPE fl,fr,f0,derf;
858  FTYPE timeterm[NDIM];
859  OPENMP3DLOOPVARSDEFINE;
860 
861 
862  if(which==CURRENTCALC1){
863  dtl=lastdt;
864  dtr=dt;
865  }
866  else if(which==CURRENTCALC2){
867  idtc=2.0/(lastdt+dt);
868  }
869 
870  idx1=1.0/(2.0*dx[1]);
871  idx2=1.0/(2.0*dx[2]);
872  idx3=1.0/(2.0*dx[3]);
873 
874 
875 
877  OPENMP3DLOOPSETUP(INP11,OUTP11,INP12,OUTP12,INP13,OUTP13);
878 #pragma omp for schedule(OPENMPSCHEDULE(),OPENMPCHUNKSIZE(blocksize))
879  OPENMP3DLOOPBLOCK{
880  OPENMP3DLOOPBLOCK2IJK(i,j,k);
881 
882 
883  get_geometry(i,j,k,CENT,ptrgeomt);
884 
885  get_geometry(i,j,k,CENT,ptrgeomc);
886 
887  // geomtp1 is same as geomt since d/dt( geometry) -> 0
888  get_geometry(ip1mac(i),j,k,CENT,ptrgeomrp1);
889  get_geometry(i,jp1mac(j),k,CENT,ptrgeomhp1);
890  get_geometry(i,j,kp1mac(k),CENT,ptrgeompp1);
891 
892  get_geometry(im1mac(i),j,k,CENT,ptrgeomrm1);
893  get_geometry(i,jm1mac(j),k,CENT,ptrgeomhm1);
894  get_geometry(i,j,km1mac(k),CENT,ptrgeompm1);
895 
896  if(which==CURRENTCALC1){
897 
898  if(calls>=2){ // since need 3 times to properly time center without having to worry about what RK is doing
899 
900  timeterm[0]=0;
901 
902  fl=GLOBALMACP0A2(cfaraday,i,j,k,4,0);
903  f0=GLOBALMACP0A2(cfaraday,i,j,k,5,0);
904  fr=GLOBALMACP0A2(cfaraday,i,j,k,0,0);
905  AA=(dtr*(fl-f0)+dtl*(fr-f0))/(dtl*dtr*(dtl+dtr));
906  BB=(-dtr*dtr*(fl-f0)+dtl*dtl*(fr-f0))/(dtl*dtr*(dtl+dtr));
907  derf=2.0*AA*dtr+BB;
908 
909  timeterm[1]=derf;
910 
911  fl=GLOBALMACP0A2(cfaraday,i,j,k,4,1);
912  f0=GLOBALMACP0A2(cfaraday,i,j,k,5,1);
913  fr=GLOBALMACP0A2(cfaraday,i,j,k,0,1);
914  AA=(dtr*(fl-f0)+dtl*(fr-f0))/(dtl*dtr*(dtl+dtr));
915  BB=(-dtr*dtr*(fl-f0)+dtl*dtl*(fr-f0))/(dtl*dtr*(dtl+dtr));
916  derf=2.0*AA*dtr+BB;
917 
918  timeterm[2]=derf;
919 
920  fl=GLOBALMACP0A2(cfaraday,i,j,k,4,2);
921  f0=GLOBALMACP0A2(cfaraday,i,j,k,5,2);
922  fr=GLOBALMACP0A2(cfaraday,i,j,k,0,2);
923  AA=(dtr*(fl-f0)+dtl*(fr-f0))/(dtl*dtr*(dtl+dtr));
924  BB=(-dtr*dtr*(fl-f0)+dtl*dtl*(fr-f0))/(dtl*dtr*(dtl+dtr));
925  derf=2.0*AA*dtr+BB;
926 
927  timeterm[3]=derf;
928  }
929  else{
930  timeterm[0]=0;
931  timeterm[1]=0;
932  timeterm[2]=0;
933  timeterm[3]=0;
934  }
935  }
936  else if(which==CURRENTCALC2){
937  // the current is calculated to end up at the zone and time edge
938  // point is to have J^t at same time as rest of J's, and same spatial points. time for J is one timestep back for all components.
939  // like current_calc_0 in time and current_calc_1 in space
940 
941  if(calls>0){ // since need 2 times
942  timeterm[0]=0;
943  timeterm[1]=(GLOBALMACP0A2(cfaraday,i,j,k,0,0)-GLOBALMACP0A2(cfaraday,i,j,k,4,0))*idtc; // F^{rt},t
944  timeterm[2]=(GLOBALMACP0A2(cfaraday,i,j,k,0,1)-GLOBALMACP0A2(cfaraday,i,j,k,4,1))*idtc; // F^{ht},t
945  timeterm[3]=(GLOBALMACP0A2(cfaraday,i,j,k,0,2)-GLOBALMACP0A2(cfaraday,i,j,k,4,2))*idtc; // F^{pt},t
946  }
947  else{
948  timeterm[0]=0;
949  timeterm[1]=0;
950  timeterm[2]=0;
951  timeterm[3]=0;
952  }
953 
954  }
955 
956  // similar to other current_calc's except position of current and spacing of derivatives
957 
958  // J^t = F^{tr},r + F^{th},h + F^{tp},p
959  GLOBALMACP0A1(jcon,i,j,k,0)=
960  +1./ptrgeomt->gdet*(ptrgeomrp1->gdet*GLOBALMACP0A2(cfaraday,ip1mac(i),j,k,1,0)-ptrgeomrm1->gdet*GLOBALMACP0A2(cfaraday,im1mac(i),j,k,1,0))*idx1 // F^{tr},r
961  +1./ptrgeomt->gdet*(ptrgeomhp1->gdet*GLOBALMACP0A2(cfaraday,i,jp1mac(j),k,2,0)-ptrgeomhm1->gdet*GLOBALMACP0A2(cfaraday,i,jm1mac(j),k,2,0))*idx2 // F^{th},h
962  +1./ptrgeomt->gdet*(ptrgeompp1->gdet*GLOBALMACP0A2(cfaraday,i,j,kp1mac(k),3,0)-ptrgeompm1->gdet*GLOBALMACP0A2(cfaraday,i,j,km1mac(k),3,0))*idx3 // F^{tp},p
963  ;
964 
965  // J^r = F^{rt},t + F^{rh},h + F^{rp},p
966  GLOBALMACP0A1(jcon,i,j,k,1)=
967  +timeterm[1]
968  +1./ptrgeomt->gdet*(ptrgeomhp1->gdet*GLOBALMACP0A2(cfaraday,i,jp1mac(j),k,2,1)-ptrgeomhm1->gdet*GLOBALMACP0A2(cfaraday,i,jm1mac(j),k,2,1))*idx2 // F^{rh},h
969  +1./ptrgeomt->gdet*(ptrgeompp1->gdet*GLOBALMACP0A2(cfaraday,i,j,kp1mac(k),3,1)-ptrgeompm1->gdet*GLOBALMACP0A2(cfaraday,i,j,km1mac(k),3,1))*idx3 // F^{rp},p
970  ;
971 
972  // J^h = F^{ht},t + F^{hr},r + F^{hp},p
973  GLOBALMACP0A1(jcon,i,j,k,2)=
974  +timeterm[2]
975  +1./ptrgeomt->gdet*(ptrgeomrp1->gdet*GLOBALMACP0A2(cfaraday,ip1mac(i),j,k,1,1)-ptrgeomrm1->gdet*GLOBALMACP0A2(cfaraday,im1mac(i),j,k,1,1))*idx1 // F^{hr},r
976  +1./ptrgeomt->gdet*(ptrgeompp1->gdet*GLOBALMACP0A2(cfaraday,i,j,kp1mac(k),3,2)-ptrgeompm1->gdet*GLOBALMACP0A2(cfaraday,i,j,km1mac(k),3,2))*idx3 // F^{hp},p
977  ;
978 
979  // J^p = F^{pt},t + F^{pr},r + F^{ph},h
980  GLOBALMACP0A1(jcon,i,j,k,3)=
981  +timeterm[3]
982  +1./ptrgeomt->gdet*(ptrgeomrp1->gdet*GLOBALMACP0A2(cfaraday,ip1mac(i),j,k,1,2)-ptrgeomrm1->gdet*GLOBALMACP0A2(cfaraday,im1mac(i),j,k,1,2))*idx1 // F^{pr},r
983  +1./ptrgeomt->gdet*(ptrgeomhp1->gdet*GLOBALMACP0A2(cfaraday,i,jp1mac(j),k,2,2)-ptrgeomhm1->gdet*GLOBALMACP0A2(cfaraday,i,jm1mac(j),k,2,2))*idx2 // F^{ph},h
984  ;
985 
986 
987  }// end loops
988 
989  }// end parallel region
990 
991 
992 
993  // static (should be outside parallel region)
994  calls++;
995  lastdt=dt;
996 
997 }
998 
999 
1000 
1001 
1005 int compute_vorticity(FTYPE (*p)[NSTORE2][NSTORE3][NPR],FTYPE (*pvort)[NSTORE2][NSTORE3][NPR],int whichpl)
1006 {
1007 
1008 
1009 #pragma omp parallel
1010  {
1011  int i,j,k;
1012  FTYPE X[NDIM],V[NDIM];
1013  struct of_geom geomdontuse;
1014  struct of_geom *ptrgeom=&geomdontuse;
1015  FTYPE dxdxp[NDIM][NDIM];
1016  FTYPE vxm,vxp,vym,vyp;
1017  OPENMP3DLOOPVARSDEFINE;
1018 
1019 
1020  // COMPLOOPVORT{
1021  OPENMP3DLOOPSETUP(INP11,OUTP11,INP12,OUTP12,INP13,OUTP13);
1022 #pragma omp for schedule(OPENMPSCHEDULE(),OPENMPCHUNKSIZE(blocksize))
1023  OPENMP3DLOOPBLOCK{
1024  OPENMP3DLOOPBLOCK2IJK(i,j,k);
1025 
1026  bl_coord_ijk_2(i, jm1mac(j), k, CENT, X, V);
1027  get_geometry(i,jm1mac(j),k,CENT,ptrgeom) ;
1028  // dx/dx' where '=prim coords (i.e. nonuni coords)
1029  vxm = MACP0A1(p,i,jm1mac(j),k,U1)*sqrt(ptrgeom->gcov[GIND(1,1)]);
1030 
1031  bl_coord_ijk_2(i, jp1mac(j), k, CENT, X, V);
1032  get_geometry(i,jp1mac(j),k,CENT,ptrgeom) ;
1033  // dx/dx' where '=prim coords (i.e. nonuni coords)
1034  vxp = MACP0A1(p,i,jp1mac(j),k,U1)*sqrt(ptrgeom->gcov[GIND(1,1)]);
1035 
1036 
1037 
1038  bl_coord_ijk_2(im1mac(i), j, k, CENT, X, V);
1039  get_geometry(im1mac(i),j,k,CENT,ptrgeom) ;
1040  // dx/dx' where '=prim coords (i.e. nonuni coords)
1041  vym = MACP0A1(p,im1mac(i),j,k,U2)*sqrt(ptrgeom->gcov[GIND(2,2)]);
1042 
1043  bl_coord_ijk_2(ip1mac(i), j, k, CENT, X, V);
1044  get_geometry(ip1mac(i),j,k,CENT,ptrgeom) ;
1045  // dx/dx' where '=prim coords (i.e. nonuni coords)
1046  vyp = MACP0A1(p,ip1mac(i),j,k,U2)*sqrt(ptrgeom->gcov[GIND(2,2)]);
1047 
1048  // center
1049  bl_coord_ijk_2(i, j, k, CENT, X, V);
1050  get_geometry(i,j,k,CENT,ptrgeom) ;
1051  dxdxprim_ijk(i,j,k, CENT, dxdxp);
1052 
1053  MACP0A1(pvort,i,j,k,whichpl) = (vyp-vym)/(2.0*dx[1])/dxdxp[1][1];
1054 
1055  MACP0A1(pvort,i,j,k,whichpl) += -(vxp-vxm)/(2.0*dx[2])/dxdxp[2][2];
1056 
1057  }
1058  }// end parallel region
1059 
1060 
1061  return(0);
1062 
1063 }
Generated on Fri May 20 2016 15:52:34 for HARM by   doxygen 1.8.3.1