HARM
harm and utilities
 All Data Structures Files Functions Variables Typedefs Macros Pages
dump.c
Go to the documentation of this file.
1 #include "decs.h"
2 
55 
56 int init_dumps(void)
57 {
58 
59  trifprintf("begin: init_dumps\n");
60 
62  //
63  // Output nprlist information to SM-readable file
64  //
66  output_nprlist_info();
67 
69  //
70  // setup number of columns per dump file
71  // (see dumpgen.c or dump.c for how used)
72  //
74  init_dnumcolumns_dnumversion();
75 
76 
77  if(mpicombine==1 && mpicombinetype==MPICOMBINEMINMEM){
79  //
80  // setup link list (only used for MINMEM method)
81  //
83  init_linklists();
84  }
85 
86  trifprintf("end: init_dumps\n");
87 
88 
89  return(0);
90 }
91 
93 void output_nprlist_info(void)
94 {
95  int pliter,pl;
96  int numversion;
97  int numlines;
98  FILE *out;
99 
100  // only CPU=0
101  if(myid==0){
102 
103  out=fopen("nprlistinfo.dat","wt");
104  if(out==NULL){
105  dualfprintf(fail_file,"Couldn't open nprlistinfo.dat\n");
106  myexit(12358235);
107  }
108 
109  numlines=7;
110  numversion=0; // version number of this file
111 
112  myfprintf(out,"%d %d\n",numlines,numversion);
113 
114  // NPR: (conserved quantities: U0-U? in SM)
115  PLOOP(pliter,pl){
116  myfprintf(out,"%d ",pl);
117  }
118  if(pliter==0) myfprintf(out,"-1"); // nothing in this list
119  myfprintf(out,"\n");
120 
121  // NPR2INTERP:
122  PINTERPLOOP(pliter,pl){
123  myfprintf(out,"%d ",pl);
124  }
125  if(pliter==0) myfprintf(out,"-1"); // nothing in this list
126  myfprintf(out,"\n");
127 
128  // NPR2NOTINTERP:
129  PNOTINTERPLOOP(pliter,pl){
130  myfprintf(out,"%d ",pl);
131  }
132  if(pliter==0) myfprintf(out,"-1"); // nothing in this list
133  myfprintf(out,"\n");
134 
135  // NPRBOUND:
136  PBOUNDLOOP(pliter,pl){
137  myfprintf(out,"%d ",pl);
138  }
139  if(pliter==0) myfprintf(out,"-1"); // nothing in this list
140  myfprintf(out,"\n");
141 
142  // NPRFLUXBOUND:
143  PFLUXBOUNDLOOP(pliter,pl){
144  myfprintf(out,"%d ",pl);
145  }
146  if(pliter==0) myfprintf(out,"-1"); // nothing in this list
147  myfprintf(out,"\n");
148 
149  // NPRDUMP:
150  PDUMPLOOP(pliter,pl){
151  myfprintf(out,"%d ",pl);
152  }
153  if(pliter==0) myfprintf(out,"-1"); // nothing in this list
154  myfprintf(out,"\n");
155 
156  // NPRINVERT:
157  PINVERTLOOP(pliter,pl){
158  myfprintf(out,"%d ",pl);
159  }
160  if(pliter==0) myfprintf(out,"-1"); // nothing in this list
161  myfprintf(out,"\n");
162 
163 
164  fclose(out);
165 
166  }// end if CPU==0
167 
168 
169 }
170 
171 
172 
174 void init_dnumcolumns_dnumversion(void)
175 {
176  char tempdumpnamelist[NUMDUMPTYPES][MAXFILENAME]=MYDUMPNAMELIST;
177  int i;
178 
179  extern void set_image_content_dnumcolumns_dnumversion(int *numcolumns, int *numversion);
180  extern void set_dump_content_dnumcolumns_dnumversion(int *numcolumns, int *numversion);
181  extern void set_gdump_content_dnumcolumns_dnumversion(int *numcolumns, int *numversion);
182  extern void set_avg_content_dnumcolumns_dnumversion(int *numcolumns, int *numversion);
183  extern void set_avg2_content_dnumcolumns_dnumversion(int *numcolumns, int *numversion);
184  extern void set_debug_content_dnumcolumns_dnumversion(int *numcolumns, int *numversion);
185  extern void set_enodebug_content_dnumcolumns_dnumversion(int *numcolumns, int *numversion);
186  extern void set_fieldline_content_dnumcolumns_dnumversion(int *numcolumns, int *numversion);
187  extern void set_dissdump_content_dnumcolumns_dnumversion(int *numcolumns, int *numversion);
188  extern void set_dumpother_content_dnumcolumns_dnumversion(int *numcolumns, int *numversion);
189  extern void set_fluxdump_content_dnumcolumns_dnumversion(int *numcolumns, int *numversion);
190  extern void set_eosdump_content_dnumcolumns_dnumversion(int *numcolumns, int *numversion);
191  extern void set_raddump_content_dnumcolumns_dnumversion(int *numcolumns, int *numversion);
192  extern void set_vpotdump_content_dnumcolumns_dnumversion(int *numcolumns, int *numversion);
193  extern void set_failfloordudump_content_dnumcolumns_dnumversion(int *numcolumns, int *numversion);
194  extern void set_dissmeasuredump_content_dnumcolumns_dnumversion(int *numcolumns, int *numversion);
195  extern void set_fluxsimpledump_content_dnumcolumns_dnumversion(int *numcolumnsvar, int *numversion);
196 
197  extern void set_rupperpoledump_read_content_dnumcolumns_dnumversion(int *numcolumns, int *numversion);
198  extern void set_rupperpoledump_content_dnumcolumns_dnumversion(int *numcolumns, int *numversion);
199 
200  extern void set_rdump_content_dnumcolumns_dnumversion(int *numcolumns, int *numversion);
201 
202  extern void set_rmetricdump_read_content_dnumcolumns_dnumversion(int *numcolumns, int *numversion);
203  extern void set_rmetricdump_content_dnumcolumns_dnumversion(int *numcolumns, int *numversion);
204 
205 
206 
207  // assign local to global -- do this since global not easily assigned to value with defs and decs approach
208  int dumpiter;
209  for(dumpiter=0;dumpiter<NUMDUMPTYPES;dumpiter++){
210  strcpy(dumpnamelist[dumpiter],tempdumpnamelist[dumpiter]);
211  }
212 
213 
214 
215  // always numcolumns=0 for fake dump
216  // version=0 shouldn't matter for fake dump
217  dnumcolumns[FAKEDUMPTYPE]=0; dnumversion[FAKEDUMPTYPE]=0;
218 
219 
220  // image
221  set_image_content_dnumcolumns_dnumversion(&dnumcolumns[IMAGEDUMPTYPE],&dnumversion[IMAGEDUMPTYPE]);
222  // dump
223  set_dump_content_dnumcolumns_dnumversion(&dnumcolumns[MAINDUMPTYPE],&dnumversion[MAINDUMPTYPE]);
224  // gdump
225  set_gdump_content_dnumcolumns_dnumversion(&dnumcolumns[GRIDDUMPTYPE],&dnumversion[GRIDDUMPTYPE]);
226  // avg
227  set_avg_content_dnumcolumns_dnumversion(&dnumcolumns[AVG1DUMPTYPE],&dnumversion[AVG1DUMPTYPE]);
228  // avg2
229  set_avg2_content_dnumcolumns_dnumversion(&dnumcolumns[AVG2DUMPTYPE],&dnumversion[AVG2DUMPTYPE]);
230  // debug
231  set_debug_content_dnumcolumns_dnumversion(&dnumcolumns[DEBUGDUMPTYPE],&dnumversion[DEBUGDUMPTYPE]);
232  // enodebug
233  set_enodebug_content_dnumcolumns_dnumversion(&dnumcolumns[ENODEBUGDUMPTYPE],&dnumversion[ENODEBUGDUMPTYPE]);
234  // fieldline
235  set_fieldline_content_dnumcolumns_dnumversion(&dnumcolumns[FIELDLINEDUMPTYPE],&dnumversion[FIELDLINEDUMPTYPE]);
236  // dissdump
237  set_dissdump_content_dnumcolumns_dnumversion(&dnumcolumns[DISSDUMPTYPE],&dnumversion[DISSDUMPTYPE]);
238  // dumpother
239  set_dumpother_content_dnumcolumns_dnumversion(&dnumcolumns[OTHERDUMPTYPE],&dnumversion[OTHERDUMPTYPE]);
240  // fluxdump
241  set_fluxdump_content_dnumcolumns_dnumversion(&dnumcolumns[FLUXDUMPTYPE],&dnumversion[FLUXDUMPTYPE]);
242  // eosdump
243  set_eosdump_content_dnumcolumns_dnumversion(&dnumcolumns[EOSDUMPTYPE],&dnumversion[EOSDUMPTYPE]);
244  // raddump
245  set_raddump_content_dnumcolumns_dnumversion(&dnumcolumns[RADDUMPTYPE],&dnumversion[RADDUMPTYPE]);
246  // vpotdump
247  set_vpotdump_content_dnumcolumns_dnumversion(&dnumcolumns[VPOTDUMPTYPE],&dnumversion[VPOTDUMPTYPE]);
248  // failfloordudump
249  set_failfloordudump_content_dnumcolumns_dnumversion(&dnumcolumns[FAILFLOORDUDUMPTYPE],&dnumversion[FAILFLOORDUDUMPTYPE]);
250  // dissmeasuredump
251  set_dissmeasuredump_content_dnumcolumns_dnumversion(&dnumcolumns[DISSMEASUREDUMPTYPE],&dnumversion[DISSMEASUREDUMPTYPE]);
252  // fluxsimpledump
253  set_fluxsimpledump_content_dnumcolumns_dnumversion(&dnumcolumns[FLUXSIMPLEDUMPTYPE],&dnumversion[FLUXSIMPLEDUMPTYPE]);
254 
255  // rdump (must come after all normal dumps since dnumcolumns used by restart to store other things needed up restart that are dealt with also above)
256  // rupperpoledump
257  set_rupperpoledump_content_dnumcolumns_dnumversion(&dnumcolumns[RESTARTUPPERPOLEDUMPTYPE],&dnumversion[RESTARTUPPERPOLEDUMPTYPE]);
258  set_rdump_content_dnumcolumns_dnumversion(&dnumcolumns[RESTARTDUMPTYPE],&dnumversion[RESTARTDUMPTYPE]);
259  // rmetricdump
260  set_rmetricdump_content_dnumcolumns_dnumversion(&dnumcolumns[RESTARTMETRICDUMPTYPE],&dnumversion[RESTARTMETRICDUMPTYPE]);
261 
262 
263 
264  trifprintf("dump number of columns(see global.nondepmnemonics.h)\n");
265  for(i=0;i<NUMDUMPTYPES;i++){
266  trifprintf("%s dnumcolumns[%d]=%d\n",dumpnamelist[i],i,dnumcolumns[i]);
267  }
268  trifprintf("\n");
269 
270 
271 }
272 
273 
274 
275 
276 
277 
279 int dump(long dump_cnt)
280 {
281  MPI_Datatype datatype;
282  int whichdump;
283  char fileprefix[MAXFILENAME]={'\0'};
284  char filesuffix[MAXFILENAME]={'\0'};
285  char fileformat[MAXFILENAME]={'\0'};
286 
287 
288  trifprintf("begin dumping dump# %ld ... ",dump_cnt);
289 
290  whichdump=MAINDUMPTYPE;
291  datatype=MPI_FTYPE;
292  strcpy(fileprefix,"dumps/dump");
293  strcpy(fileformat,"%04ld"); //atch adjust dump every substep
294  strcpy(filesuffix,"");
295 
296  if(dump_gen(WRITEFILE,dump_cnt,binaryoutput,whichdump,datatype,fileprefix,fileformat,filesuffix,dump_header,dump_content)>=1) return(1);
297 
299  //writesyminfo();
301 
302  trifprintf("end dumping dump# %ld ... ",dump_cnt);
303 
304 
305  return(0);
306 
307 }
308 
309 
310 
312 void writesyminfo( void )
313 {
314  int i, j;
315 
316  for( i = 0; i < N1; i++ ) {
317  }
318 
319 }
320 
322 int dump_header(int whichdump, int whichdumpversion, int numcolumnsvar, int bintxt, FILE *headerptr)
323 {
324  int dump_header_general(int whichdump, int whichdumpversion, int numcolumnsvar, long localrealnstep, SFTYPE localdt, int bintxt, FILE *headerptr);
325  int retval;
326 
327  retval=dump_header_general(whichdump, whichdumpversion, numcolumnsvar, realnstep,dt, bintxt, headerptr);
328 
329  return(retval);
330 
331 }
332 
334 int fluxdump_header(int whichdump, int whichdumpversion, int numcolumnsvar, int bintxt, FILE *headerptr)
335 {
336  int dump_header_general(int whichdump, int whichdumpversion, int numcolumnsvar, long localrealnstep, SFTYPE localdt, int bintxt, FILE *headerptr);
337  int retval;
338 
339  retval=dump_header_general(whichdump, whichdumpversion, numcolumnsvar, fluxdumprealnstep,fluxdumpdt, bintxt, headerptr);
340 
341  return(retval);
342 
343 }
344 
346 int dump_header_general(int whichdump, int whichdumpversion, int numcolumnsvar, long localrealnstep, SFTYPE localdt, int bintxt, FILE *headerptr)
347 {
348  int realtotalsize[NDIM];
349  FTYPE realstartx[NDIM];
350  FTYPE X[NDIM];
351  int is,ie,js,je,ks,ke;
352 
353  // set computational domain for output header
354  is=Uconsevolveloop[FIS]+SHIFTX1DN;
355  ie=Uconsevolveloop[FIE]+SHIFTX1UP;
356  js=Uconsevolveloop[FJS]+SHIFTX2DN;
357  je=Uconsevolveloop[FJE]+SHIFTX2UP;
358  ks=Uconsevolveloop[FKS]+SHIFTX3DN;
359  ke=Uconsevolveloop[FKE]+SHIFTX3UP;
360 
361 
362  // get real startx's (assumes rectangular grid)
363  realtotalsize[1]=totalsize[1]+2*EXTRADUMP1;
364  if(EXTRADUMP1!=0){
365  coord(0-EXTRADUMP1,0,0,FACE1,X);
366  realstartx[1]=X[1];
367  }
368  else realstartx[1]=startx[1];
369 
370  realtotalsize[2]=totalsize[2]+2*EXTRADUMP2;
371  if(EXTRADUMP2!=0){
372  coord(0,0-EXTRADUMP2,0,FACE2,X);
373  realstartx[2]=X[2];
374  }
375  else realstartx[2]=startx[2];
376 
377  realtotalsize[3]=totalsize[3]+2*EXTRADUMP3;
378  if(EXTRADUMP3!=0){
379  coord(0,0,0-EXTRADUMP3,FACE3,X);
380  realstartx[3]=X[3];
381  }
382  else realstartx[3]=startx[3];
383 
384 
385  // dx is the same (constant)
386 
387  // 15+3=18 elements total
388  if(bintxt==BINARYOUTPUT){
389  fwrite(&tsteppartf,sizeof(FTYPE),1,headerptr);
390  fwrite(&realtotalsize[1],sizeof(int),1,headerptr);
391  fwrite(&realtotalsize[2],sizeof(int),1,headerptr);
392  fwrite(&realtotalsize[3],sizeof(int),1,headerptr);
393  fwrite(&realstartx[1],sizeof(FTYPE),1,headerptr);
394  fwrite(&realstartx[2],sizeof(FTYPE),1,headerptr);
395  fwrite(&realstartx[3],sizeof(FTYPE),1,headerptr);
396  fwrite(&dx[1],sizeof(FTYPE),1,headerptr);
397  fwrite(&dx[2],sizeof(FTYPE),1,headerptr);
398  fwrite(&dx[3],sizeof(FTYPE),1,headerptr);
399  fwrite(&localrealnstep,sizeof(long),1,headerptr);
400  fwrite(&gam,sizeof(FTYPE),1,headerptr);
401  fwrite(&a,sizeof(FTYPE),1,headerptr);
402  fwrite(&R0,sizeof(FTYPE),1,headerptr);
403  fwrite(&Rin,sizeof(FTYPE),1,headerptr);
404  fwrite(&Rout,sizeof(FTYPE),1,headerptr);
405  fwrite(&hslope,sizeof(FTYPE),1,headerptr);
406  fwrite(&localdt,sizeof(FTYPE),1,headerptr);
407  fwrite(&defcoord,sizeof(int),1,headerptr);
408  // SUPERGODMARK: For jon_interp to not need init.h specific stuff, need to at least add MCOORD so knows if spherical polar coordinates or not. defcoord is insufficient in general.
409  // Maybe also: COMPDIM, WHICHVEL, WHICHEOM, REMOVERESTMASSFROMUU, EOMTYPE, NPR, DOENTROPY,
410  fwrite(&MBH,sizeof(FTYPE),1,headerptr);
411  fwrite(&QBH,sizeof(FTYPE),1,headerptr);
412  fwrite(&EP3,sizeof(FTYPE),1,headerptr);
413  fwrite(&THETAROT,sizeof(FTYPE),1,headerptr);
414  fwrite(&is,sizeof(int),1,headerptr);
415  fwrite(&ie,sizeof(int),1,headerptr);
416  fwrite(&js,sizeof(int),1,headerptr);
417  fwrite(&je,sizeof(int),1,headerptr);
418  fwrite(&ks,sizeof(int),1,headerptr);
419  fwrite(&ke,sizeof(int),1,headerptr);
420  fwrite(&whichdump,sizeof(int),1,headerptr);
421  fwrite(&whichdumpversion,sizeof(int),1,headerptr);
422  fwrite(&numcolumnsvar,sizeof(int),1,headerptr);
423 
424  // init.h type stuff related to physical settings or what data is outputted
425  int var;
426  var=TRACKVPOT; fwrite(&var,sizeof(int),1,headerptr);
427  var=MCOORD; fwrite(&var,sizeof(int),1,headerptr);
428  var=DODISS; fwrite(&var,sizeof(int),1,headerptr);
429  var=DOEVOLVEMETRIC; fwrite(&var,sizeof(int),1,headerptr);
430  var=WHICHVEL; fwrite(&var,sizeof(int),1,headerptr);
431  var=WHICHEOM; fwrite(&var,sizeof(int),1,headerptr);
432  var=REMOVERESTMASSFROMUU; fwrite(&var,sizeof(int),1,headerptr);
433  var=RELTYPE; fwrite(&var,sizeof(int),1,headerptr);
434  var=EOMTYPE; fwrite(&var,sizeof(int),1,headerptr);
435  var=WHICHEOS; fwrite(&var,sizeof(int),1,headerptr);
436  var=DOENTROPY; fwrite(&var,sizeof(int),1,headerptr);
437  var=WHICHENTROPYEVOLVE; fwrite(&var,sizeof(int),1,headerptr);
438  var=CALCFARADAYANDCURRENTS; fwrite(&var,sizeof(int),1,headerptr);
439  var=DOPOLEDEATH; fwrite(&var,sizeof(int),1,headerptr);
440  var=DOPOLESMOOTH; fwrite(&var,sizeof(int),1,headerptr);
441  var=DOPOLEGAMMADEATH; fwrite(&var,sizeof(int),1,headerptr);
442  var=IF3DSPCTHENMPITRANSFERATPOLE; fwrite(&var,sizeof(int),1,headerptr);
443  var=EOMRADTYPE; fwrite(&var,sizeof(int),1,headerptr);
444  var=WHICHRADSOURCEMETHOD; fwrite(&var,sizeof(int),1,headerptr);
445  var=OUTERDEATH; fwrite(&var,sizeof(int),1,headerptr);
446  var=OUTERDEATHRADIUS; fwrite(&var,sizeof(int),1,headerptr);
447 
448 
449 
450  }
451  else{
452 
453 #define DUMPHEADERLIST tsteppartf, realtotalsize[1], realtotalsize[2], realtotalsize[3], realstartx[1], realstartx[2], realstartx[3], dx[1], dx[2], dx[3], localrealnstep,gam,a,R0,Rin,Rout,hslope,localdt,defcoord,MBH,QBH,EP3,THETAROT,is,ie,js,je,ks,ke,whichdump,whichdumpversion,numcolumnsvar, TRACKVPOT ,MCOORD ,DODISS ,DOEVOLVEMETRIC ,WHICHVEL ,WHICHEOM ,REMOVERESTMASSFROMUU ,RELTYPE ,EOMTYPE ,WHICHEOS ,DOENTROPY ,WHICHENTROPYEVOLVE ,CALCFARADAYANDCURRENTS ,DOPOLEDEATH ,DOPOLESMOOTH ,DOPOLEGAMMADEATH ,IF3DSPCTHENMPITRANSFERATPOLE ,EOMRADTYPE ,WHICHRADSOURCEMETHOD ,OUTERDEATH ,OUTERDEATHRADIUS
454 
455 #if(REALTYPE==DOUBLETYPE)
456  fprintf(headerptr, "%21.15g %d %d %d %21.15g %21.15g %21.15g %21.15g %21.15g %21.15g %ld %21.15g %21.15g %21.15g %21.15g %21.15g %21.15g %21.15g %d %21.15g %21.15g %21.15g %21.15g %d %d %d %d %d %d %d %d %d %d %d %d %d %d %d %d %d %d %d %d %d %d %d %d %d %d %d %d %d %21.15g\n", DUMPHEADERLIST);
457 #elif(REALTYPE==LONGDOUBLETYPE)
458  fprintf(headerptr, "%31.25Lg %d %d %d %31.25Lg %31.25Lg %31.25Lg %31.25Lg %31.25Lg %31.25Lg %ld %31.25Lg %31.25Lg %31.25Lg %31.25Lg %31.25Lg %31.25Lg %31.25Lg %d %31.25Lg %31.25Lg %31.25Lg %31.25Lg %d %d %d %d %d %d %d %d %d %d %d %d %d %d %d %d %d %d %d %d %d %d %d %d %d %d %d %d %d %31.25Lg\n", DUMPHEADERLIST);
459 #endif
460  }
461  fflush(headerptr);
462  return(0);
463 }
464 
465 
466 
468 void set_dump_content_dnumcolumns_dnumversion(int *numcolumnsvar, int *numversion)
469 {
470 
471 
472 
473  if(DOMAINDUMPDIAG){
474  // always NPRDUMP
475  if(GAMMIEDUMP) *numcolumnsvar=2*3 + NPRDUMP+NPR + 3 + 1 + 4 * NDIM + 6 + 1
476 #if(CALCFARADAYANDCURRENTS)
477  + NDIM*2
478  + 2*6
479 #endif
480  ;
481  else{
482  *numcolumnsvar=3*3 + NPRDUMP + 3 + (nprend+1) + 1 + 4 * NDIM + 6 + 1 //replace NPR -> (nprend+1) since nprend, not NPR, controls dumping. Fixes: DOEXTRAINTERP=1 case
483 #if(CALCFARADAYANDCURRENTS)
484  + NDIM*2
485  + 2*6
486 #endif
487  ; // 61 total if also doing currents and faraday, 41 otherwise
488 
489  if(FLUXB==FLUXCTSTAG && 0){ // DEBUG (change corresponding code in dump.c)
490  *numcolumnsvar+= NPR2INTERP*COMPDIM*2 + NPR + COMPDIM*3*2 + COMPDIM*3*2*2;
491  }
492  }
493  }
494  else{
495  *numcolumnsvar=0;
496  }
497 
498  // Version number:
499  *numversion=0;
500 
501 }
502 
503 
504 
506 int dump_content(int i, int j, int k, MPI_Datatype datatype,void *writebuf)
507 {
508  int pl,pliter;
509  FTYPE r, th, vmin[NDIM], vmax[NDIM];
510  int ignorecourant;
511  struct of_state qdontuse;
512  struct of_state *qptr=&qdontuse;
513  FTYPE X[NDIM],V[NDIM];
514  FTYPE divb;
515  FTYPE b[NDIM],ucon[NDIM];
516  FTYPE U[NPR];
517  FTYPE ftemp;
518  FTYPE jcov[NDIM];
519  FTYPE fcov[NUMFARADAY];
520  FTYPE rho,u,pressure,cs2,Sden;
521  int dir,l,m,n,o;
522  struct of_geom geomdontuse;
523  struct of_geom *ptrgeom=&geomdontuse;
524  int loc=CENT;
525 
527  //
528  // some calculations
529  //
530 
531  bl_coord_ijk_2(i, j, k, loc, X,V);
532  // if failed, then data output for below invalid, but columns still must exist
533 
534  get_geometry(i, j, k, loc, ptrgeom);
535 
536  if (!failed) {
537  if (get_state(GLOBALMAC(pdump,i,j,k), ptrgeom, qptr) >= 1)
538  FAILSTATEMENT("dump.c:dump()", "get_state() dir=0", 1);
539  if(N1NOT1){
540  if (vchar_all(GLOBALMAC(pdump,i,j,k), qptr, 1, ptrgeom, &vmax[1], &vmin[1],&ignorecourant) >= 1) FAILSTATEMENT("dump.c:dump()", "vchar_all() dir=1or2", 1);
541  }
542  else{
543  vmax[1]=vmin[1]=0;
544  }
545  if(N2NOT1){
546  if (vchar_all(GLOBALMAC(pdump,i,j,k), qptr, 2, ptrgeom, &vmax[2], &vmin[2],&ignorecourant) >= 1) FAILSTATEMENT("dump.c:dump()", "vchar_all() dir=1or2", 2);
547  }
548  else{
549  vmax[2]=vmin[2]=0;
550  }
551  if(N3NOT1){
552  if (vchar_all(GLOBALMAC(pdump,i,j,k), qptr, 3, ptrgeom, &vmax[3], &vmin[3],&ignorecourant) >= 1) FAILSTATEMENT("dump.c:dump()", "vchar_all() dir=1or2", 3);
553  }
554  else{
555  vmax[3]=vmin[3]=0;
556  }
557 
558  }
559  else {// do a per zone check, otherwise set to 0
560  whocalleducon=1; // force no failure mode, just return like failure, and don't return if failure, just set to 0 and continue
561  if (get_state(GLOBALMAC(pdump,i,j,k), ptrgeom, qptr) >= 1){
562  for (pl = 0; pl < NDIM; pl++)
563  qptr->ucon[pl]=0;
564  for (pl = 0; pl < NDIM; pl++)
565  qptr->ucov[pl]=0;
566  for (pl = 0; pl < NDIM; pl++)
567  qptr->bcon[pl]=0;
568  for (pl = 0; pl < NDIM; pl++)
569  qptr->bcov[pl]=0;
570  }
571  if (vchar_all(GLOBALMAC(pdump,i,j,k), qptr, 1, ptrgeom, &vmax[1], &vmin[1],&ignorecourant) >= 1){
572  vmax[1]=vmin[1]=0;
573  }
574 
575  if (vchar_all(GLOBALMAC(pdump,i,j,k), qptr, 2, ptrgeom, &vmax[2], &vmin[2],&ignorecourant) >= 1){
576  vmax[2]=vmin[2]=0;
577  }
578 
579  if (vchar_all(GLOBALMAC(pdump,i,j,k), qptr, 3, ptrgeom, &vmax[3], &vmin[3],&ignorecourant) >= 1){
580  vmax[3]=vmin[3]=0;
581  }
582 
583  whocalleducon=0; // return to normal state
584 
585  }
586 
587 
588  // GODMARK: use of globals: Ok for dumps since probably always globally dumped, not dumped per section computed or something similar
589  setfdivb(&divb, GLOBALPOINT(pdump), GLOBALPOINT(pstagdump), GLOBALPOINT(udump), GLOBALPOINT(Bhatdump), i, j, k); // udump also set externally GODMARK
590 
592  //
593  // do the assignments
594  //
595  // if you change # of outputted vars, remember to change numcolumnsvar
596 
597 
598  //static
599  if(!GAMMIEDUMP){
600  ftemp=(FTYPE)(i+startpos[1]);
601  myset(datatype,&ftemp,0,1,writebuf); //ti
602  ftemp=(FTYPE)(j+startpos[2]);
603  myset(datatype,&ftemp,0,1,writebuf); //tj
604  ftemp=(FTYPE)(k+startpos[3]);
605  myset(datatype,&ftemp,0,1,writebuf); //tk
606  }
607  myset(datatype,X,1,3,writebuf); //x1, x2, x3
608  myset(datatype,V,1,3,writebuf); //r, h, ph
609  // 9
610 
612  //
613  // rest dynamic
614 
615  // primitives
616  // must use PDUMPLOOP() since may be any order unlike NPR loop
617  PDUMPLOOP(pliter,pl) myset(datatype,&(GLOBALMACP0A1(pdump,i,j,k,pl)),0,1,writebuf); // NPRDUMP //rho u v1 v2 v3 B1 B2 B3 ??
618 
620  //
621  // output some EOS stuff since in general not simple function of rho0,u
622  rho = GLOBALMACP0A1(pdump,i,j,k,RHO);
623  u = GLOBALMACP0A1(pdump,i,j,k,UU);
624 
625 
626  pressure = pressure_rho0_u_simple(i,j,k,loc,rho,u);
627  cs2 = cs2_compute_simple(i,j,k,loc,rho,u);
628  Sden = compute_entropy_simple(i,j,k,loc,rho,u);
629 
630  myset(datatype,&pressure,0,1,writebuf); // 1 //p
631  myset(datatype,&cs2,0,1,writebuf); // 1 //cs2
632  myset(datatype,&Sden,0,1,writebuf); // 1 //Sden
633 
635  //
636  // output the conserved quantities since not easily inverted and at higher order aren't invertable from point primitives
637  // PLOOP() used since conserved quantities always fill full PLOOP, while PDUMPLOOP is for primitives that may be duplicate among conserved quantities
638  PLOOP(pliter,pl) myset(datatype,&(GLOBALMACP0A1(udump,i,j,k,pl)),0,1,writebuf); // NPR //U0 U1 U2 U3 U4 U5 U6 U7 ??
639  myset(datatype,&divb,0,1,writebuf); // 1 //divb
640 
641  for (pl = 0; pl < NDIM; pl++)
642  myset(datatype,&(qptr->ucon[pl]),0,1,writebuf); //uu0 uu1 uu2 uu3
643  for (pl = 0; pl < NDIM; pl++)
644  myset(datatype,&(qptr->ucov[pl]),0,1,writebuf); //ud0 ud1 ud2 ud3
645  for (pl = 0; pl < NDIM; pl++)
646  myset(datatype,&(qptr->bcon[pl]),0,1,writebuf); //bu0 bu1 bu2 bu3
647  for (pl = 0; pl < NDIM; pl++)
648  myset(datatype,&(qptr->bcov[pl]),0,1,writebuf); //bd0 bd1 bd2 bd3
649  // 4*4
650 
651  myset(datatype,&vmin[1],0,1,writebuf); //v1m v1p v2m v2p v3m v3p
652  myset(datatype,&vmax[1],0,1,writebuf);
653  myset(datatype,&vmin[2],0,1,writebuf);
654  myset(datatype,&vmax[2],0,1,writebuf);
655  myset(datatype,&vmin[3],0,1,writebuf);
656  myset(datatype,&vmax[3],0,1,writebuf);
657  // 6
658 
659  // one static term
660  myset(datatype,&(ptrgeom->gdet),0,1,writebuf); // 1 //gdet //end of default read
661 
662 
663  if(CALCFARADAYANDCURRENTS){ // NIM*2+6*2 = 8+12=20
664  // updated 11/16/2003
665  // new 10/23/2003
666  // current density
667  lower_vec(GLOBALMAC(jcon,i,j,k),ptrgeom,jcov);
668  myset(datatype,GLOBALMAC(jcon,i,j,k),0,NDIM,writebuf); // (NDIM)
669  myset(datatype,jcov,0,NDIM,writebuf);// (NDIM)
670  // faraday (2*6)
671  lowerf(GLOBALMAC(fcon,i,j,k),ptrgeom,fcov);
672  myset(datatype,GLOBALMAC(fcon,i,j,k),0,NUMFARADAY,writebuf); // (6)
673  myset(datatype,fcov,0,NUMFARADAY,writebuf); // (6)
674  }
675 
676 
677  // DEBUG: Also add +3 to numcolumnsvar for this to work
678  if(0){
679  if(FLUXB==FLUXCTSTAG) myset(datatype,GLOBALMAC(pstagdump,i,j,k),B1,3,writebuf);
680  else{
681  FTYPE plblob[NPR]={0};
682  myset(datatype,plblob,B1,3,writebuf);
683  }
684  }
685 
686 
687  if(FLUXB==FLUXCTSTAG && 0){ // DEBUG (change corresponding code in dump.c)
688  // uses jrdp3dudebug in gtwod.m that assumes CALCFARADAYANDCURRENTS==0
689  for(l=1;l<=COMPDIM;l++) myset(datatype,GLOBALMACP1A0(pl_ctdump,l,i,j,k),0,NPR2INTERP,writebuf); // 3*8 = 24
690  for(l=1;l<=COMPDIM;l++) myset(datatype,GLOBALMACP1A0(pr_ctdump,l,i,j,k),0,NPR2INTERP,writebuf); // 3*8 = 24
691  myset(datatype,GLOBALMAC(pstagdump,i,j,k),0,NPR,writebuf); // 8 // GODMARK: use of globals
692  // for(dir=1;dir<=COMPDIM;dir++) for(pl=B1;pl<=B3;pl++) for(n=0;n<=1;n++) myset(datatype,&GLOBALMACP3A0(pbcorninterp,dir,pl,n,i,j,k),0,1,writebuf); // 3*3*2 = 18
693  for(dir=1;dir<=COMPDIM;dir++) for(pl=1;pl<=COMPDIM;pl++) for(n=0;n<NUMCS+1;n++) for(o=0;o<NUMCS;o++) myset(datatype,&GLOBALMACP1A3(pvbcorninterp,dir,i,j,k,pl,n,o),0,1,writebuf); // 3*3*(2+1)*2 = 54 (was 36)
694  }
695 
696  return (0);
697 }
698 
699 
700 
701 
702 
704 int debugdump(long dump_cnt)
705 {
706  MPI_Datatype datatype;
707  int whichdump;
708  char fileprefix[MAXFILENAME]={'\0'};
709  char filesuffix[MAXFILENAME]={'\0'};
710  char fileformat[MAXFILENAME]={'\0'};
711 
712 
713  trifprintf("begin dumping debug dump# %ld ... ",dump_cnt);
714 
715  whichdump=DEBUGDUMPTYPE;
716  datatype=MPI_CTYPE;
717  strcpy(fileprefix,"dumps/debugdump");
718  strcpy(fileformat,"%04ld");
719  strcpy(filesuffix,"");
720 
721  // same header as dump
722  if(dump_gen(WRITEFILE,dump_cnt,binaryoutput,whichdump,datatype,fileprefix,fileformat,filesuffix,dump_header,debug_content)>=1) return(1);
723 
724  trifprintf("end dumping debug# %ld ... ",dump_cnt);
725 
726  return(0);
727 
728 }
729 
731 extern void set_debug_content_dnumcolumns_dnumversion(int *numcolumnsvar, int *numversion)
732 {
733 
734  if(DODEBUG && DODEBUGDUMP){ //different than doing just DODEBUG
735  *numcolumnsvar=2*NUMFAILFLOORFLAGS*NUMTSCALES;
736  }
737  else *numcolumnsvar=0;
738 
739  // Version number:
740  *numversion=0;
741 
742 }
743 
745 int debug_content(int i, int j, int k, MPI_Datatype datatype,void *writebuf)
746 {
747  // could also make everything FTYPE and convert like for normal i,j dump file
748 
749 
750  // NOTEMARK: see also restart.c since this is added to restart
751  myset(datatype,GLOBALMAC(failfloorcount,i,j,k),0,2*NUMTSCALES*NUMFAILFLOORFLAGS,writebuf);
752 
753 
754  return(0);
755 }
756 
757 
758 
759 
761 int enodebugdump(long dump_cnt)
762 {
763  MPI_Datatype datatype;
764  int whichdump;
765  char fileprefix[MAXFILENAME]={'\0'};
766  char filesuffix[MAXFILENAME]={'\0'};
767  char fileformat[MAXFILENAME]={'\0'};
768 
769 
770  trifprintf("begin dumping enodebug dump# %ld ... ",dump_cnt);
771 
772  whichdump=ENODEBUGDUMPTYPE;
773  // datatype=MPI_FTYPE;
774  datatype=MPI_CTYPE;
775  strcpy(fileprefix,"dumps/enodebugdump");
776  strcpy(fileformat,"%04ld");
777  strcpy(filesuffix,"");
778 
779  // same header as dump
780  if(dump_gen(WRITEFILE,dump_cnt,binaryoutput,whichdump,datatype,fileprefix,fileformat,filesuffix,eno_dump_header,enodebug_content)>=1) return(1);
781 
782  trifprintf("end dumping enodebug# %ld ... ",dump_cnt);
783 
784  return(0);
785 
786 }
787 
788 
790 void set_enodebug_content_dnumcolumns_dnumversion(int *numcolumnsvar, int *numversion)
791 {
792 
793  if(DOENODEBUG){
794  //*numcolumnsvar=NUMENODEBUGS;
795  *numcolumnsvar=(3-1)* NUMENOINTERPTYPES * (NPR-4) * NUMENODEBUGS; //SASMARK2
796  }
797  else *numcolumnsvar=0;
798 
799  // Version number:
800  *numversion=0;
801 
802 
803 }
804 
806 int enodebug_content(int i, int j, int k, MPI_Datatype datatype,void *writebuf)
807 {
808  // could also make everything FTYPE and convert like for normal i,j dump file
809  //myset(datatype,GLOBALMAC(enodebugarray,i,j,k),0,3*NUMENOINTERPTYPES*NPR*NUMENODEBUGS,writebuf);
810  myset(datatype,GLOBALMAC(enodebugarray,i,j,k),0,(3-1)*NUMENOINTERPTYPES*(NPR-4)*NUMENODEBUGS,writebuf); //atch corrected
811 
812  return(0);
813 }
814 
816 int eno_dump_header(int whichdump, int whichdumpversion, int numcolumnsvar, int bintxt, FILE *headerptr)
817 {
818  int dump_header_general(int whichdump, int whichdumpversion, int numcolumnsvar, long localrealnstep, SFTYPE localdt, int bintxt, FILE *headerptr);
819  int retval;
820 
821  retval=dump_header_general(whichdump, whichdumpversion, numcolumnsvar, realnstep,dt, bintxt, headerptr);
822 
823  return(retval);
824 }
825 
826 
827 
828 
829 
830 
831 
833 int avgdump(long dump_cnt)
834 {
835  MPI_Datatype datatype;
836  int whichdump;
837  char fileprefix[MAXFILENAME]={'\0'};
838  char filesuffix[MAXFILENAME]={'\0'};
839  char fileformat[MAXFILENAME]={'\0'};
840 
841 
842 
843  trifprintf("begin dumping avgdump# %ld ... ",dump_cnt);
844 
845  whichdump=AVG1DUMPTYPE;
846  datatype=MPI_FTYPE;
847  strcpy(fileprefix,"dumps/avg");
848  strcpy(fileformat,"%04ld");
849  strcpy(filesuffix,"");
850 
851  if(dump_gen(WRITEFILE,dump_cnt,binaryoutput,whichdump,datatype,fileprefix,fileformat,filesuffix,dump_header,avg_content)>=1) return(1);
852 
853  trifprintf("end dumping avgdump# %ld ... ",dump_cnt);
854 
855 
856  return(0);
857 
858 }
859 
861 void set_avg_content_dnumcolumns_dnumversion(int *numcolumnsvar, int *numversion)
862 {
863 
864  if(DOAVG){
865  // 36+29+8*2+4*2+2+12*2+96*2=339
866  *numcolumnsvar=3*3 + 1 + NUMNORMDUMP // (6+1+29=36)
867  + NUMNORMDUMP // |normal terms| (29)
868 #if(CALCFARADAYANDCURRENTS)
869  + NDIM*2 // jcon/jcov (8)
870  + NDIM*2 // |jcon|/|jcov| (8)
871 #endif
872  + NDIM*2 // massflux/|massflux|
873  + NUMOTHER*2 // other stuff and fabs of each
874 #if(CALCFARADAYANDCURRENTS)
875  +6*2 // fcon/fcov (12)
876  +6*2 // |fcon|,|fcov| (12)
877 #endif
878  +7*16 // Tud all 7 parts, all 4x4 terms (112)
879  +7*16 // |Tud| all 7 parts, all 4x4 terms (112)
880  ;
881 
882  if(DOAVG2){
883  *numcolumnsvar-=224;
884  }
885  }
886  else{
887  *numcolumnsvar=0;
888  }
889 
890  // Version number:
891  *numversion=0;
892 
893 
894 
895 }
896 
897 
899 int avg_content(int i, int j, int k, MPI_Datatype datatype,void *writebuf)
900 {
901  int pl = 0, l = 0, col = 0;
902  FTYPE X[NDIM],V[NDIM];
903  FTYPE ftemp;
904  struct of_geom geomdontuse;
905  struct of_geom *ptrgeom=&geomdontuse;
906  int loc=CENT;
907 
908 
909  bl_coord_ijk_2(i, j, k, loc, X,V);
910  get_geometry(i, j, k, loc, ptrgeom);
911 
912  if(!GAMMIEDUMP){
913  ftemp=(FTYPE)(i+startpos[1]);
914  myset(datatype,&ftemp,0,1,writebuf);
915  ftemp=(FTYPE)(j+startpos[2]);
916  myset(datatype,&ftemp,0,1,writebuf);
917  ftemp=(FTYPE)(k+startpos[3]);
918  myset(datatype,&ftemp,0,1,writebuf);
919  }
920  myset(datatype,X,1,3,writebuf);
921  myset(datatype,V,1,3,writebuf);
922 
923 
924  myset(datatype,&(ptrgeom->gdet),0,1,writebuf);
925 
926  // now do time average stuff
927  myset(datatype,GLOBALMAC(normalvarstavg,i,j,k),0,NUMNORMDUMP,writebuf);
928  myset(datatype,GLOBALMAC(anormalvarstavg,i,j,k),0,NUMNORMDUMP,writebuf);
929 
930 #if(CALCFARADAYANDCURRENTS)
931  myset(datatype,GLOBALMAC(jcontavg,i,j,k),0,NDIM,writebuf);
932  myset(datatype,GLOBALMAC(jcovtavg,i,j,k),0,NDIM,writebuf);
933  myset(datatype,GLOBALMAC(ajcontavg,i,j,k),0,NDIM,writebuf);
934  myset(datatype,GLOBALMAC(ajcovtavg,i,j,k),0,NDIM,writebuf);
935 #endif
936  myset(datatype,GLOBALMAC(massfluxtavg,i,j,k),0,NDIM,writebuf);
937  myset(datatype,GLOBALMAC(amassfluxtavg,i,j,k),0,NDIM,writebuf);
938 
939  myset(datatype,GLOBALMAC(othertavg,i,j,k),0,NUMOTHER,writebuf);
940  myset(datatype,GLOBALMAC(aothertavg,i,j,k),0,NUMOTHER,writebuf);
941 
942 #if(CALCFARADAYANDCURRENTS)
943  myset(datatype,GLOBALMAC(fcontavg,i,j,k),0,NUMFARADAY,writebuf);
944  myset(datatype,GLOBALMAC(fcovtavg,i,j,k),0,NUMFARADAY,writebuf);
945  myset(datatype,GLOBALMAC(afcontavg,i,j,k),0,NUMFARADAY,writebuf);
946  myset(datatype,GLOBALMAC(afcovtavg,i,j,k),0,NUMFARADAY,writebuf);
947 #endif
948 
949 #if(DOAVG2==0)
950  myset(datatype,GLOBALMAC(tudtavg,i,j,k),0,NUMSTRESSTERMS,writebuf);
951  myset(datatype,GLOBALMAC(atudtavg,i,j,k),0,NUMSTRESSTERMS,writebuf);
952 #endif
953 
954  return(0);
955 
956 }
957 
959 int avg2dump(long dump_cnt)
960 {
961  MPI_Datatype datatype;
962  int whichdump;
963  char fileprefix[MAXFILENAME]={'\0'};
964  char filesuffix[MAXFILENAME]={'\0'};
965  char fileformat[MAXFILENAME]={'\0'};
966 
967 
968 
969  trifprintf("begin dumping avg2dump# %ld ... ",dump_cnt);
970 
971  whichdump=AVG2DUMPTYPE;
972  datatype=MPI_FTYPE;
973  strcpy(fileprefix,"dumps/avg2");
974  strcpy(fileformat,"%04ld");
975  strcpy(filesuffix,"");
976 
977  if(dump_gen(WRITEFILE,dump_cnt,binaryoutput,whichdump,datatype,fileprefix,fileformat,filesuffix,dump_header,avg2_content)>=1) return(1);
978 
979  trifprintf("end dumping avg2dump# %ld ... ",dump_cnt);
980 
981 
982  return(0);
983 
984 }
985 
987 void set_avg2_content_dnumcolumns_dnumversion(int *numcolumnsvar, int *numversion)
988 {
989 
990 
991  if(DOAVG2){
992  *numcolumnsvar=10 + 224; // otherwise doesn't exist so don't need to set
993  }
994  else *numcolumnsvar=0;
995 
996  // Version number:
997  *numversion=0;
998 
999 
1000 
1001 }
1002 
1004 int avg2_content(int i, int j, int k, MPI_Datatype datatype,void *writebuf)
1005 {
1006  int pl = 0, l = 0, col = 0;
1007  FTYPE X[NDIM],V[NDIM];
1008  FTYPE ftemp;
1009  struct of_geom geomdontuse;
1010  struct of_geom *ptrgeom=&geomdontuse;
1011  int loc=CENT;
1012 
1013  bl_coord_ijk_2(i, j, k, loc, X,V);
1014  get_geometry(i, j, k, loc, ptrgeom);
1015  // if you change # of outputted vars, remember to change numcolumnsvar above
1016 
1017  if(!GAMMIEDUMP){
1018  ftemp=(FTYPE)(i+startpos[1]);
1019  myset(datatype,&ftemp,0,1,writebuf);
1020  ftemp=(FTYPE)(j+startpos[2]);
1021  myset(datatype,&ftemp,0,1,writebuf);
1022  ftemp=(FTYPE)(k+startpos[3]);
1023  myset(datatype,&ftemp,0,1,writebuf);
1024  }
1025  myset(datatype,X,1,3,writebuf);
1026  myset(datatype,V,1,3,writebuf);
1027 
1028  myset(datatype,&(ptrgeom->gdet),0,1,writebuf);
1029  // 10
1030 
1031  myset(datatype,GLOBALMAC(tudtavg,i,j,k),0,NUMSTRESSTERMS,writebuf);
1032  myset(datatype,GLOBALMAC(atudtavg,i,j,k),0,NUMSTRESSTERMS,writebuf);
1033  // 112*2
1034 
1035  // total=10+112*2=234
1036 
1037  return(0);
1038 }
1039 
1040 
1041 
1043 int gdump(long dump_cnt)
1044 {
1045  MPI_Datatype datatype;
1046  int whichdump;
1047  char fileprefix[MAXFILENAME]={'\0'};
1048  char filesuffix[MAXFILENAME]={'\0'};
1049  char fileformat[MAXFILENAME]={'\0'};
1050 
1051 
1052  trifprintf("begin dumping gdump# %ld ... ",dump_cnt);
1053 
1054  whichdump=GRIDDUMPTYPE;
1055  datatype=MPI_FTYPE;
1056  strcpy(fileprefix,"dumps/gdump");
1057  strcpy(fileformat,"%04ld");
1058  strcpy(filesuffix,"");
1059 
1060  // dump_cnt==-1 means no file number
1061  if(dump_gen(WRITEFILE,dump_cnt,binaryoutput,whichdump,datatype,fileprefix,fileformat,filesuffix,dump_header,gdump_content)>=1) return(1);
1062 
1063  trifprintf("end dumping gdump# %ld ... ",dump_cnt);
1064 
1065  return(0);
1066 }
1067 
1068 
1070 extern void set_gdump_content_dnumcolumns_dnumversion(int *numcolumnsvar, int *numversion)
1071 {
1072 
1073  if(DOGDUMP){
1074  // 205+4+4*4 currently
1075  //*numcolumnsvar=3*3+NDIM*NDIM*NDIM+NPG*NDIM*NDIM*2+NPG+4+4*4;
1076  //NPG was replaced with unity in order to avoid excessive dumping of info (only center info now available)
1077  *numcolumnsvar=3*3 + NDIM*NDIM*NDIM + 1*NDIM*NDIM*2 + 1 + NDIM + NDIM*NDIM;
1078  //t^i x^i V^i, \Gamma^\mu_{\nu\tau}, g^{\mu\nu} g_{\mu\nu}, \sqrt{-g}, \gamma_\mu, dx^\mu/dx^\nu
1079  }
1080  else{
1081  *numcolumnsvar=0;
1082  }
1083  // Version number:
1084  *numversion=0;
1085 
1086 }
1087 
1088 
1090 int gdump_content(int i, int j, int k, MPI_Datatype datatype, void *writebuf)
1091 {
1092  int pl = 0, l = 0, m = 0, n = 0, col = 0;
1093  FTYPE X[NDIM],V[NDIM];
1094  FTYPE ftemp;
1095  FTYPE *ptrftemp;
1096  FTYPE dxdxp[NDIM][NDIM];
1097  int myii,myjj,mykk;
1098  LOCALMETRICTEMPVARS;
1099  int loc=CENT;
1100  FTYPE generalmatrixlower[NDIM][NDIM];
1101  FTYPE generalmatrixupper[NDIM][NDIM];
1102  int jj,kk;
1103 
1104 
1105 
1106  bl_coord_ijk_2(i, j, k, loc, X, V);
1107  dxdxprim_ijk(i, j, k, loc, dxdxp);
1108 
1109 
1110 
1111  ftemp=(FTYPE)(i+startpos[1]);
1112  myset(datatype,&ftemp,0,1,writebuf);
1113  ftemp=(FTYPE)(j+startpos[2]);
1114  myset(datatype,&ftemp,0,1,writebuf);
1115  ftemp=(FTYPE)(k+startpos[3]);
1116  myset(datatype,&ftemp,0,1,writebuf);
1117  // 3
1118  myset(datatype,X,1,3,writebuf);
1119  myset(datatype,V,1,3,writebuf);
1120  // 6
1121 
1122 
1123 
1124 
1125 #if(MCOORD!=CARTMINKMETRIC)
1126  myii=i;
1127  myjj=j;
1128  mykk=k;
1129 #else
1130  myii=0;
1131  myjj=0;
1132  mykk=0;
1133 #endif
1134 
1135 
1136  ptrftemp=(FTYPE*)(&GLOBALMETMACP0A3(conn,myii,myjj,mykk,0,0,0));
1137  myset(datatype,ptrftemp,0,NDIM*NDIM*NDIM,writebuf);
1138 
1139  // get local metric quantities for this loc,i,j,k
1140  GETLOCALMETRIC(loc,myii,myjj,mykk);
1141 
1142  DLOOP(jj,kk){
1143  generalmatrixlower[jj][kk]=localgcov[GIND(jj,kk)];
1144  generalmatrixupper[jj][kk]=localgcon[GIND(jj,kk)];
1145  }
1146 
1147 
1148  ptrftemp=&generalmatrixupper[0][0];
1149  myset(datatype,ptrftemp,0,NDIM*NDIM,writebuf);
1150  ptrftemp=&generalmatrixlower[0][0];
1151  myset(datatype,ptrftemp,0,NDIM*NDIM,writebuf);
1152  ptrftemp=&localgdet[0];
1153  myset(datatype,ptrftemp,0,1,writebuf);
1154  //ptrftemp=(FTYPE*)(&localgdetvol); // can take a peek if GDETVOLDIFF==1
1155  // myset(datatype,ptrftemp,0,1,writebuf);
1156 
1157 
1158  ptrftemp=&GLOBALMETMACP0A1(conn2,myii,myjj,mykk,0);
1159  myset(datatype,ptrftemp,0,NDIM,writebuf);
1160 
1161  // 4*4
1162  ptrftemp=&dxdxp[0][0];
1163  myset(datatype,ptrftemp,0,NDIM*NDIM,writebuf);
1164 
1165 
1166  return(0);
1167 
1168 }
1169 
1170 
1172 int fieldlinedump(long dump_cnt)
1173 {
1174  MPI_Datatype datatype;
1175  int whichdump;
1176  char fileprefix[MAXFILENAME]={'\0'};
1177  char filesuffix[MAXFILENAME]={'\0'};
1178  char fileformat[MAXFILENAME]={'\0'};
1179 
1180 
1181 
1182  trifprintf("begin dumping fieldlinedump# %ld ... ",dump_cnt);
1183 
1184  whichdump=FIELDLINEDUMPTYPE;
1185  datatype=MPI_FLOAT; // don't need good precision
1186  strcpy(fileprefix,"dumps/fieldline");
1187  strcpy(fileformat,"%04ld");
1188  strcpy(filesuffix,"");
1189 
1190  // MIXEDOUTPUT means text header and forced binary data
1191  if(dump_gen(WRITEFILE,dump_cnt,MIXEDOUTPUT,whichdump,datatype,fileprefix,fileformat,filesuffix,dump_header,fieldline_content)>=1) return(1);
1192 
1193  trifprintf("end dumping fieldlinedump# %ld ... ",dump_cnt);
1194 
1195 
1196  return(0);
1197 
1198 }
1199 
1201 extern void set_fieldline_content_dnumcolumns_dnumversion(int *numcolumnsvar, int *numversion)
1202 {
1203 
1204 
1208 #if( FIELDLINEGDETB == 1)
1209 #define NUMFIELDLINEQUANTITIES (14-2 + NUMYFL*(DOYFL!=0) + (DOYL!=0) + (DOYNU!=0) + (1+NDIM+11)*(EOMRADTYPE!=EOMRADNONE))
1210 
1211 
1212 #else
1213 #define NUMFIELDLINEQUANTITIES (11-2 + NUMYFL*(DOYFL!=0) + (DOYL!=0) + (DOYNU!=0) + (1+NDIM+11)*(EOMRADTYPE!=EOMRADNONE))
1214 
1215 
1216 #endif
1217 
1218 
1219 
1220  if(DOFIELDLINE){
1221  *numcolumnsvar=NUMFIELDLINEQUANTITIES;
1222  }
1223  else *numcolumnsvar=0;
1224 
1225 
1226 
1227 
1228  // Version number:
1229  *numversion=0;
1230 
1231 }
1232 
1234 int fieldline_content(int i, int j, int k, MPI_Datatype datatype,void *writebuf)
1235 {
1236  int pl = 0, l = 0, col = 0;
1237  struct of_state q;
1238  //FTYPE U[NPR];
1239  FTYPE FL[NPR];
1240  // must be same precision as written content
1241  float ftemp;
1242  struct of_geom geomdontuse;
1243  struct of_geom *ptrgeom=&geomdontuse;
1244  int loc=CENT;
1245 
1246 
1247  FTYPE *pr=GLOBALMAC(pdump,i,j,k);
1248 
1250  //
1251  // some calculations
1252  //
1253 
1254  // if failed, then data output for below invalid, but columns still must exist
1255  get_geometry(i, j, k, loc, ptrgeom);
1256  if (!failed) {
1257  if (get_state(GLOBALMAC(pdump,i,j,k), ptrgeom, &q) >= 1)
1258  FAILSTATEMENT("dump.c:dump()", "get_state() dir=0", 1);
1259  }
1260  else {// do a per zone check, otherwise set to 0
1261  whocalleducon=1; // force no failure mode, just return like failure, and don't return if failure, just set to 0 and continue
1262  if (get_state(GLOBALMAC(pdump,i,j,k), ptrgeom, &q) >= 1){
1263  for (pl = 0; pl < NDIM; pl++) q.ucon[pl]=0;
1264  for (pl = 0; pl < NDIM; pl++) q.ucov[pl]=0;
1265  for (pl = 0; pl < NDIM; pl++) q.bcon[pl]=0;
1266  for (pl = 0; pl < NDIM; pl++) q.bcov[pl]=0;
1267 
1268  if(EOMRADTYPE!=EOMRADNONE){
1269  for (pl = 0; pl < NDIM; pl++) q.uradcon[pl]=0;
1270  for (pl = 0; pl < NDIM; pl++) q.uradcov[pl]=0;
1271  }
1272  }
1273  whocalleducon=0; // return to normal state
1274 
1275  }
1276 
1277  MYFUN(primtoflux(UDIAG,pr, &q, RR, ptrgeom, FL, NULL),"step_ch.c:fluxcalc()", "primtoflux_calc() dir=1/2 l", RR);
1278 
1279 
1281  //
1282  // do the assignments
1283  //
1284  // if you change # of outputted vars, remember to change numcolumnsvar
1285 
1286 
1287 
1289  //
1290  // 2 various things
1291 
1292  // rho (for various things)
1293  ftemp=(float)GLOBALMACP0A1(pdump,i,j,k,RHO);
1294  myset(datatype,&ftemp,0,1,writebuf); // 1
1295 
1296  // u (for various things)
1297  ftemp=(float)GLOBALMACP0A1(pdump,i,j,k,UU);
1298  myset(datatype,&ftemp,0,1,writebuf); // 1
1299 
1300  if(DOYFL){
1301  // Y_fl
1302  if(YFL1>=0){
1303  ftemp=(float)(GLOBALMACP0A1(pdump,i,j,k,YFL1));
1304  myset(datatype,&ftemp,0,1,writebuf); // 1
1305  }
1306  if(YFL2>=0){
1307  ftemp=(float)(GLOBALMACP0A1(pdump,i,j,k,YFL2));
1308  myset(datatype,&ftemp,0,1,writebuf); // 1
1309  }
1310  if(YFL3>=0){
1311  ftemp=(float)(GLOBALMACP0A1(pdump,i,j,k,YFL3));
1312  myset(datatype,&ftemp,0,1,writebuf); // 1
1313  }
1314  if(YFL4>=0){
1315  ftemp=(float)(GLOBALMACP0A1(pdump,i,j,k,YFL4));
1316  myset(datatype,&ftemp,0,1,writebuf); // 1
1317  }
1318  if(YFL5>=0){
1319  ftemp=(float)(GLOBALMACP0A1(pdump,i,j,k,YFL5));
1320  myset(datatype,&ftemp,0,1,writebuf); // 1
1321  }
1322  }
1323 
1324  if(DOYL){
1325  // Y_l
1326  ftemp=(float)(GLOBALMACP0A1(pdump,i,j,k,YL));
1327  myset(datatype,&ftemp,0,1,writebuf); // 1
1328  }
1329 
1330  if(DOYNU){
1331  // Y_nu
1332  ftemp=(float)(GLOBALMACP0A1(pdump,i,j,k,YNU));
1333  myset(datatype,&ftemp,0,1,writebuf); // 1
1334  }
1335 
1337  //
1338  // 2 things for jet/energy per baryon at infinity
1339 
1340  // -u_t (-hu_t can be found from this and rho/u/p above)
1341  // ftemp=(float)(-q.ucov[0]);
1342  // myset(datatype,&ftemp,0,1,writebuf); // 1
1343 
1344  // -T^t_t/(rho u^t)
1345  // ftemp=(float)(-U[UU]/(ptrgeom->gdet * GLOBALMACP0A1(pdump,i,j,k,RHO)*q.ucon[TT]));
1346  //myset(datatype,&ftemp,0,1,writebuf);
1347 
1348  // -T^r_t/(rho u^r)
1349  // if(q.ucon[RR]!=0.0){
1350  // ftemp=(float)(-FL[UU]/(ptrgeom->gdet * GLOBALMACP0A1(pdump,i,j,k,RHO)*q.ucon[RR]));
1351  // }
1352  // else ftemp=0.0;
1353  // myset(datatype,&ftemp,0,1,writebuf); // 1
1354 
1355 
1356  // 1 extra thing
1357 
1358  // u^t
1359  ftemp=(float)(q.ucon[0]);
1360  myset(datatype,&ftemp,0,1,writebuf); // 1
1361 
1362 
1364  //
1365  // 6 things for the field line stuff
1366 
1367  // v^r [ in grid frame]
1368  ftemp=(float)(q.ucon[1]/q.ucon[0]);
1369  myset(datatype,&ftemp,0,1,writebuf); // 1
1370 
1371  // v^\theta
1372  ftemp=(float)(q.ucon[2]/q.ucon[0]);
1373  myset(datatype,&ftemp,0,1,writebuf); // 1
1374 
1375  // v^\phi
1376  ftemp=(float)(q.ucon[3]/q.ucon[0]);
1377  myset(datatype,&ftemp,0,1,writebuf); // 1
1378 
1379  // B^r
1380  ftemp=(float)(GLOBALMACP0A1(pdump,i,j,k,B1));
1381  myset(datatype,&ftemp,0,1,writebuf); // 1
1382 
1383  // B^\theta
1384  ftemp=(float)(GLOBALMACP0A1(pdump,i,j,k,B2));
1385  myset(datatype,&ftemp,0,1,writebuf); // 1
1386 
1387  // B^\phi
1388  ftemp=(float)(GLOBALMACP0A1(pdump,i,j,k,B3));
1389  myset(datatype,&ftemp,0,1,writebuf); // 1
1390 
1391 #if( FIELDLINEGDETB == 1)
1392  //it is useful to have access to gdet*B^i at cell faces directly for plotting field lines
1393  ftemp=(float)(GLOBALMACP0A1(udump,i,j,k,B1));
1394  myset(datatype,&ftemp,0,1,writebuf); // 1
1395 
1396  ftemp=(float)(GLOBALMACP0A1(udump,i,j,k,B2));
1397  myset(datatype,&ftemp,0,1,writebuf); // 1
1398 
1399  ftemp=(float)(GLOBALMACP0A1(udump,i,j,k,B3));
1400  myset(datatype,&ftemp,0,1,writebuf); // 1
1401 #endif
1402 
1403  // see grmhd-dualfcon2omegaf.nb
1404  // below can be obtained from above set of v and B
1405  // \Omega_F_1
1406  // ftemp=(float)(v3-B3*v2/(B2+SMALL));
1407  //myset(datatype,&ftemp,0,1,writebuf); // 1
1408 
1409  // \Omega_F_2
1410  //ftemp=(float)(v3-B3*v1/(B1+SMALL));
1411  // myset(datatype,&ftemp,0,1,writebuf); // 1
1412 
1413 
1414  if(EOMRADTYPE!=EOMRADNONE){
1415 
1416  // Erf (for various things)
1417  ftemp=(float)GLOBALMACP0A1(pdump,i,j,k,PRAD0);
1418  myset(datatype,&ftemp,0,1,writebuf); // 1
1419 
1420  // urad^t [ in grid frame]
1421  ftemp=(float)(q.uradcon[0]);
1422  myset(datatype,&ftemp,0,1,writebuf); // 1
1423 
1424  // vrad^r [ in grid frame]
1425  ftemp=(float)(q.uradcon[1]/q.uradcon[0]);
1426  myset(datatype,&ftemp,0,1,writebuf); // 1
1427 
1428  // vrad^\theta
1429  ftemp=(float)(q.uradcon[2]/q.uradcon[0]);
1430  myset(datatype,&ftemp,0,1,writebuf); // 1
1431 
1432  // vrad^\phi
1433  ftemp=(float)(q.uradcon[3]/q.uradcon[0]);
1434  myset(datatype,&ftemp,0,1,writebuf); // 1
1435 
1436 
1437 
1438 
1439 
1440 
1441  // add kappaes, kappa, kappan, kappaemit, kappanemit so don't have to keep python scripts up to date with form or how opacity is dealth with.
1442 
1443  //radiative stress tensor in the lab frame
1444  FTYPE Rij[NDIM][NDIM];
1445 
1446  //this call returns R^i_j, i.e., the first index is contra-variant and the last index is co-variant
1447  mhdfull_calc_rad(pr, ptrgeom, &q, Rij);
1448 
1449  //the four-velocity of fluid in lab frame
1450  FTYPE *ucon,*ucov;
1451  ucon = q.ucon;
1452  ucov = q.ucov;
1453 
1454  //Eradff = R^a_b u_a u^b
1455  FTYPE Ruu=0.; DLOOP(i,j) Ruu+=Rij[i][j]*ucov[i]*ucon[j];
1456  // get relative Lorentz factor between gas and radiation
1457  FTYPE gammaradgas = 0.0;
1458  int jj;
1459  DLOOPA(jj) gammaradgas += - (q.ucov[jj] * q.uradcon[jj]);
1460 
1461  // get B
1462  FTYPE bsq = dot(q.bcon, q.bcov);
1463  FTYPE B=sqrt(bsq);
1464 
1465  FTYPE nradff=0;
1466  FTYPE expfactorradff=0;
1467  FTYPE kappa,kappan=0;
1468  FTYPE kappaemit,kappanemit=0;
1469  FTYPE kappaes;
1470  FTYPE lambda,nlambda=0;
1471  FTYPE Tgas=0,Tradff=0;
1472  calc_Tandopacityandemission(pr,ptrgeom,&q,Ruu,gammaradgas,B,&Tgas,&Tradff,&nradff,&expfactorradff,&kappa,&kappan,&kappaemit,&kappanemit,&kappaes, &lambda, &nlambda);
1473 
1474  myset(datatype,&Tgas,0,1,writebuf); // 1
1475  myset(datatype,&Tradff,0,1,writebuf); // 1
1476  myset(datatype,&nradff,0,1,writebuf); // 1
1477  myset(datatype,&expfactorradff,0,1,writebuf); // 1
1478  myset(datatype,&kappa,0,1,writebuf); // 1
1479  myset(datatype,&kappan,0,1,writebuf); // 1
1480  myset(datatype,&kappaemit,0,1,writebuf); // 1
1481  myset(datatype,&kappanemit,0,1,writebuf); // 1
1482  myset(datatype,&kappaes,0,1,writebuf); // 1
1483  myset(datatype,&lambda,0,1,writebuf); // 1
1484  myset(datatype,&nlambda,0,1,writebuf); // 1
1485 
1486  }
1487 
1488 
1489 
1490  return(0);
1491 
1492 }
1493 
1494 
1495 
1497 int dissdump(long dump_cnt)
1498 {
1499  MPI_Datatype datatype;
1500  int whichdump;
1501  char fileprefix[MAXFILENAME]={'\0'};
1502  char filesuffix[MAXFILENAME]={'\0'};
1503  char fileformat[MAXFILENAME]={'\0'};
1504 
1505 
1506  trifprintf("begin dumping dissdump# %ld ... ",dump_cnt);
1507 
1508  whichdump=DISSDUMPTYPE;
1509  datatype=MPI_FTYPE;
1510  strcpy(fileprefix,"dumps/dissdump");
1511  strcpy(fileformat,"%04ld");
1512  strcpy(filesuffix,"");
1513 
1514  if(dump_gen(WRITEFILE,dump_cnt,binaryoutput,whichdump,datatype,fileprefix,fileformat,filesuffix,dump_header,dissdump_content)>=1) return(1);
1515 
1516  trifprintf("end dumping dissdump# %ld ... ",dump_cnt);
1517 
1518 
1519  return(0);
1520 
1521 }
1522 
1524 void set_dissdump_content_dnumcolumns_dnumversion(int *numcolumnsvar, int *numversion)
1525 {
1526 
1527  if(DODISS){
1528  *numcolumnsvar=NUMDISSFUNPOS;
1529  }
1530  else *numcolumnsvar=0;
1531 
1532  // Version number:
1533  *numversion=0;
1534 
1535 }
1536 
1537 
1539 int dissdump_content(int i, int j, int k, MPI_Datatype datatype,void *writebuf)
1540 {
1541 
1542  // NOTEMARK: see also restart.c since this is added to restart
1543  myset(datatype,&GLOBALMAC(dissfunpos,i,j,k),0,NUMDISSFUNPOS,writebuf);
1544 
1545  return (0);
1546 }
1547 
1548 
1550 int dumpother(long dump_cnt)
1551 {
1552  MPI_Datatype datatype;
1553  int whichdump;
1554  char fileprefix[MAXFILENAME]={'\0'};
1555  char filesuffix[MAXFILENAME]={'\0'};
1556  char fileformat[MAXFILENAME]={'\0'};
1557 
1558 
1559  trifprintf("begin dumping dumpother# %ld ... ",dump_cnt);
1560 
1561  whichdump=OTHERDUMPTYPE;
1562  datatype=MPI_FTYPE;
1563  strcpy(fileprefix,"dumps/dumpother");
1564  strcpy(fileformat,"%04ld"); //atch adjust dump every substep
1565  strcpy(filesuffix,"");
1566 
1567  if(dump_gen(WRITEFILE,dump_cnt,binaryoutput,whichdump,datatype,fileprefix,fileformat,filesuffix,dump_header,dumpother_content)>=1) return(1);
1568 
1570  //writesyminfo();
1572 
1573  trifprintf("end dumping dumpother# %ld ... ",dump_cnt);
1574 
1575 
1576  return(0);
1577 
1578 }
1579 
1581 void set_dumpother_content_dnumcolumns_dnumversion(int *numcolumnsvar, int *numversion)
1582 {
1583 
1584  if(DODUMPOTHER){ // panalytic + numpother quantities
1585  *numcolumnsvar=NPR+NUMPOTHER;
1586  }
1587  else *numcolumnsvar=0;
1588 
1589  // Version number:
1590  *numversion=0;
1591 
1592 
1593 }
1594 
1596 int dumpother_content(int i, int j, int k, MPI_Datatype datatype,void *writebuf)
1597 {
1598  int pl,pliter;
1599 
1600 #if(ANALYTICMEMORY==0)
1601  dualfprintf(fail_file,"Need to set ANALYTICMEMORY==1 for using analytical memory in dumpother_content()\n");
1602  myexit(9857652);
1603 #endif
1604 
1605 
1606  myset(datatype,&GLOBALMAC(panalytic,i,j,k),0,NPR,writebuf);
1607 
1608  for(pl=0;pl<NUMPOTHER;pl++){
1609  myset(datatype,&GLOBALMACP1A0(pother,pl,i,j,k),0,1,writebuf);
1610  }
1611 
1612  return (0);
1613 }
1614 
1615 
1616 
1618 int fluxdumpdump(long dump_cnt)
1619 {
1620  MPI_Datatype datatype;
1621  int whichdump;
1622  char fileprefix[MAXFILENAME]={'\0'};
1623  char filesuffix[MAXFILENAME]={'\0'};
1624  char fileformat[MAXFILENAME]={'\0'};
1625 
1626 
1627  trifprintf("begin dumping fluxdump# %ld ... ",dump_cnt);
1628 
1629  whichdump=FLUXDUMPTYPE;
1630  datatype=MPI_FTYPE;
1631  strcpy(fileprefix,"dumps/fluxdump");
1632  strcpy(fileformat,"%04ld"); //atch adjust dump every substep
1633  strcpy(filesuffix,"");
1634 
1635  if(dump_gen(WRITEFILE,dump_cnt,binaryoutput,whichdump,datatype,fileprefix,fileformat,filesuffix,fluxdump_header,fluxdump_content)>=1) return(1);
1636 
1637  trifprintf("end dumping fluxdump# %ld ... ",dump_cnt);
1638 
1639 
1640  return(0);
1641 
1642 }
1643 
1645 void set_fluxdump_content_dnumcolumns_dnumversion(int *numcolumnsvar, int *numversion)
1646 {
1647 
1648  if(FLUXDUMP==1){ // dU, flux, and ppprimitives for flux
1649  *numcolumnsvar=NUMFLUXDUMP;
1650  }
1651  else *numcolumnsvar=0;
1652 
1653 
1654 
1655  // Version number:
1656  *numversion=0;
1657 
1658 
1659 }
1660 
1661 
1663 int fluxdump_content(int i, int j, int k, MPI_Datatype datatype,void *writebuf)
1664 {
1665  int pl,pliter;
1666 
1667  myset(datatype,&GLOBALMAC(fluxdump,i,j,k),0,NUMFLUXDUMP,writebuf);
1668 
1669  return (0);
1670 }
1671 
1672 
1673 
1675 int eosdump(long dump_cnt)
1676 {
1677  MPI_Datatype datatype;
1678  int whichdump;
1679  char fileprefix[MAXFILENAME]={'\0'};
1680  char filesuffix[MAXFILENAME]={'\0'};
1681  char fileformat[MAXFILENAME]={'\0'};
1682 
1683 
1684  trifprintf("begin dumping eosdump# %ld ... ",dump_cnt);
1685 
1686  whichdump=EOSDUMPTYPE;
1687  datatype=MPI_FTYPE;
1688  strcpy(fileprefix,"dumps/eosdump");
1689  strcpy(fileformat,"%04ld"); //atch adjust dump every substep
1690  strcpy(filesuffix,"");
1691 
1692  if(dump_gen(WRITEFILE,dump_cnt,binaryoutput,whichdump,datatype,fileprefix,fileformat,filesuffix,dump_header,eosdump_content)>=1) return(1);
1693 
1695  //writesyminfo();
1697 
1698  trifprintf("end dumping eosdump# %ld ... ",dump_cnt);
1699 
1700 
1701  return(0);
1702 
1703 }
1704 
1705 
1707 void set_eosdump_content_dnumcolumns_dnumversion(int *numcolumnsvar, int *numversion)
1708 {
1709 
1710  if(DOEOSDIAG==1){
1711  // all EOSs output same size data so uniform format
1712  // otherwise have to also put this condition in dump.c when outputting so don't overwrite memory!
1713  *numcolumnsvar=MAXPARLIST+1+MAXNUMEXTRAS+MAXPROCESSEDEXTRAS; // 1 is temperature
1714  }
1715  else{
1716  *numcolumnsvar=0;
1717  }
1718 
1719  // Version number:
1720  *numversion=0;
1721 
1722 
1723 }
1724 
1726 int eosdump_content(int i, int j, int k, MPI_Datatype datatype,void *writebuf)
1727 {
1728  FTYPE rho,u;
1729  FTYPE height,ye,ynu,temp;
1730  // NOTEMARK: have to fill with zero or else might not be number and if ldouble is used then seg faults
1731  FTYPE extras[MAXNUMEXTRAS]={0};
1732  FTYPE processed[MAXPROCESSEDEXTRAS]={0};
1733  int numextras;
1734  int extraiter;
1735  int numparms;
1736  FTYPE parlist[MAXPARLIST]={0};
1737  int loc=CENT;
1738 
1739 
1740  // dualfprintf(fail_file,"NUMEOSGLOBALS=%g LASTEOSGLOBAL=%g FIRSTEOSGLOBAL=%g\n",NUMEOSGLOBALS,LASTEOSGLOBAL,FIRSTEOSGLOBAL);
1741 
1743  //
1744  // do the assignments
1745  //
1746  // if you change # of outputted vars, remember to change numcolumnsvar
1747 
1748 
1750  //
1751  // output some EOS stuff since in general not simple function of rho0,u
1752  rho = GLOBALMACP0A1(pdump,i,j,k,RHO);
1753  u = GLOBALMACP0A1(pdump,i,j,k,UU);
1754 
1755  // compute EOS stuff
1756  get_EOS_parms_simple(&numparms, i,j,k,loc, parlist);
1757  temp = compute_temp_simple(i,j,k,loc,rho,u);
1758  // get extra EOS stuff
1759  get_extrasprocessed_simple(1, i, j, k, loc, GLOBALMAC(pdump,i,j,k), extras, processed);
1760  // compute_allextras(0,rho,u,&numextras,extras);
1761 
1762 
1763  // write EOS stuff
1764  // write out all stuff for all EOSs so uniform format to read for any run (otherwise need version information in file)
1765  myset(datatype,&parlist,0,MAXPARLIST,writebuf);
1766  myset(datatype,&temp,0,1,writebuf); // 1
1767  // write extras EOS stuff
1768  myset(datatype,extras,0,MAXNUMEXTRAS,writebuf); // numextras
1769  myset(datatype,processed,0,MAXPROCESSEDEXTRAS,writebuf); // numprocessedextras
1770 
1771 
1772  return (0);
1773 }
1774 
1775 
1776 
1777 
1778 
1779 
1780 
1781 
1782 
1783 
1784 
1785 
1786 
1787 
1789 int raddump(long dump_cnt)
1790 {
1791  MPI_Datatype datatype;
1792  int whichdump;
1793  char fileprefix[MAXFILENAME]={'\0'};
1794  char filesuffix[MAXFILENAME]={'\0'};
1795  char fileformat[MAXFILENAME]={'\0'};
1796 
1797 
1798  trifprintf("begin dumping raddump# %ld ... ",dump_cnt);
1799 
1800  whichdump=RADDUMPTYPE;
1801  datatype=MPI_FTYPE;
1802  strcpy(fileprefix,"dumps/raddump");
1803  strcpy(fileformat,"%04ld"); //atch adjust dump every substep
1804  strcpy(filesuffix,"");
1805 
1806  if(dump_gen(WRITEFILE,dump_cnt,binaryoutput,whichdump,datatype,fileprefix,fileformat,filesuffix,dump_header,raddump_content)>=1) return(1);
1807 
1809  //writesyminfo();
1811 
1812  trifprintf("end dumping raddump# %ld ... ",dump_cnt);
1813 
1814 
1815  return(0);
1816 
1817 }
1818 
1820 void set_raddump_content_dnumcolumns_dnumversion(int *numcolumnsvar, int *numversion)
1821 {
1822 
1823  if(EOMRADTYPE!=EOMRADNONE && DORADDIAG){
1824  *numcolumnsvar=NDIM*5 + 6 + 9 + NDIM+1 + 4*3;
1825  }
1826  else{
1827  *numcolumnsvar=0;
1828  }
1829 
1830  // Version number:
1831  *numversion=0;
1832 
1833 
1834 }
1835 
1837 int raddump_content(int i, int j, int k, MPI_Datatype datatype,void *writebuf)
1838 {
1839  int loc=CENT;
1840  FTYPE *pr=&GLOBALMACP0A1(pdump,i,j,k,0);
1841 
1842 
1843  // leave if shouldn't write anything
1844  if(EOMRADTYPE==EOMRADNONE) return(0);
1845 
1847  //
1848  // do the assignments
1849  //
1850  // if you change # of outputted vars, remember to change numcolumnsvar
1851 
1852 
1853  // get X,V
1854  FTYPE X[NDIM],V[NDIM];
1855  bl_coord_ijk_2(i, j, k, loc, X,V);
1856  // get geometry
1857  struct of_geom geomdontuse;
1858  struct of_geom *ptrgeom=&geomdontuse;
1859  get_geometry(i, j, k, loc, ptrgeom);
1860 
1861  // get full state
1862  struct of_state q;
1863  get_state(pr,ptrgeom,&q);
1864 
1865  myset(datatype,q.uradcon,0,NDIM,writebuf); // NDIM
1866  myset(datatype,q.uradcov,0,NDIM,writebuf); // NDIM
1867 
1868  // Transform these lab frame coordinate basis primitives to fluid frame E,F^i
1869  int whichvel=VEL4; // desired final ff version
1870  int whichcoord=MCOORD; // desired final ff version
1871  FTYPE pradffortho[NPR]={-1};
1872  FTYPE prff[NPR]; // not new compared to pr
1873  prad_fforlab(&whichvel, &whichcoord, LAB2FF, i,j,k,CENT,ptrgeom,pradffortho, pr, prff);
1874  myset(datatype,pradffortho,PRAD0,NDIM,writebuf); // NDIM
1875 
1876  // get 4-force in lab and fluid frame
1877  FTYPE Gdpl[NPR]={0},Gdabspl[NPR]={0},chi,Tgas,Trad;
1878  int computestate=0;// already computed above
1879  int computeentropy=1;
1880  koral_source_rad_calc(computestate,computeentropy,pr, ptrgeom, Gdpl, Gdabspl, &chi, &Tgas, &Trad, &q);
1881  myset(datatype,Gdpl,PRAD0,NDIM,writebuf); // NDIM
1882  myset(datatype,Gdabspl,PRAD0,NDIM,writebuf); // NDIM
1883 
1884  // get gas T
1885  myset(datatype,&Tgas,0,1,writebuf); // 1
1886 
1887  // get radiation T from actual Erf
1888  FTYPE Tradlte=calc_LTE_TfromE(pr[PRAD0]);
1889  myset(datatype,&Tradlte,0,1,writebuf); // 1
1890 
1891  // get radiation's fluid frame T from actual Erf in fluid frame
1892  FTYPE Tradff=calc_LTE_TfromE(pradffortho[PRAD0]);
1893  myset(datatype,&Tradff,0,1,writebuf); // 1
1894 
1895  // get Trad -- full fluid frame T
1896  myset(datatype,&Trad,0,1,writebuf); // 1
1897 
1898 
1899  //radiative stress tensor in the lab frame
1900  FTYPE Rij[NDIM][NDIM];
1901 
1902  //this call returns R^i_j, i.e., the first index is contra-variant and the last index is co-variant
1903  mhdfull_calc_rad(pr, ptrgeom, &q, Rij);
1904 
1905  //the four-velocity of fluid in lab frame
1906  FTYPE *ucon,*ucov;
1907  ucon = q.ucon;
1908  ucov = q.ucov;
1909 
1910  //Eradff = R^a_b u_a u^b
1911  FTYPE Ruu=0.; DLOOP(i,j) Ruu+=Rij[i][j]*ucov[i]*ucon[j];
1912  // get relative Lorentz factor between gas and radiation
1913  FTYPE gammaradgas = 0.0;
1914  int jj;
1915  DLOOPA(jj) gammaradgas += - (q.ucov[jj] * q.uradcon[jj]);
1916 
1917  // get Erf in fluid frame
1918  myset(datatype,&Ruu,0,1,writebuf); // 1
1919 
1920  // get Erf [assuming LTE]
1921  FTYPE Erf;
1922  Erf=calc_LTE_Efromurho(pr[UU],pr[RHO]);
1923  myset(datatype,&Erf,0,1,writebuf); // 1
1924 
1925  // get B
1926  FTYPE bsq = dot(q.bcon, q.bcov);
1927  FTYPE B=sqrt(bsq);
1928 
1929  // get lambda
1930  // FTYPE Tgas=calc_PEQ_Tfromurho(pr[UU],pr[RHO]);
1931  // calc_rad_lambda(pr, ptrgeom, Tgas, &lambda, &nlambda, &kappaemit, &kappanemit);
1932  // get absorption opacities
1933  // FTYPE Tgas;
1934  FTYPE nradff=0;
1935  FTYPE expfactorradff=0;
1936  FTYPE kappa,kappan=0;
1937  FTYPE kappaemit,kappanemit=0;
1938  FTYPE kappaes;
1939  FTYPE lambda,nlambda=0;
1940  calc_Tandopacityandemission(pr,ptrgeom,&q,Ruu,gammaradgas,B,&Tgas,&Tradff,&nradff,&expfactorradff,&kappa,&kappan,&kappaemit,&kappanemit,&kappaes, &lambda, &nlambda);
1941 
1942  myset(datatype,&nradff,0,1,writebuf); // 1
1943  myset(datatype,&expfactorradff,0,1,writebuf); // 1
1944  myset(datatype,&kappa,0,1,writebuf); // 1
1945  myset(datatype,&kappan,0,1,writebuf); // 1
1946  myset(datatype,&kappaemit,0,1,writebuf); // 1
1947  myset(datatype,&kappanemit,0,1,writebuf); // 1
1948  myset(datatype,&kappaes,0,1,writebuf); // 1
1949  myset(datatype,&lambda,0,1,writebuf); // 1
1950  myset(datatype,&nlambda,0,1,writebuf); // 1
1951 
1952  // get tau
1953  FTYPE tautot[NDIM]={0},tautotmax;
1954  calc_tautot(pr, ptrgeom, &q, tautot, &tautotmax);
1955  myset(datatype,tautot,0,NDIM,writebuf); // NDIM
1956  myset(datatype,&tautotmax,0,1,writebuf); // 1
1957 
1958 
1959 
1960  // get radiative vchar
1961  int ignorecourant;
1962  FTYPE vmin1=0,vmax1=0,vmax21=0,vmin21=0;
1963  vchar_rad(pr, &q, 1, ptrgeom, &vmax1, &vmin1, &vmax21, &vmin21, &ignorecourant);
1964  myset(datatype,&vmin1,0,1,writebuf); // 1
1965  myset(datatype,&vmax1,0,1,writebuf); // 1
1966  myset(datatype,&vmin21,0,1,writebuf); // 1
1967  myset(datatype,&vmax21,0,1,writebuf); // 1
1968  FTYPE vmin2=0,vmax2=0,vmax22=0,vmin22=0;
1969  vchar_rad(pr, &q, 2, ptrgeom, &vmax2, &vmin2, &vmax22, &vmin22, &ignorecourant);
1970  myset(datatype,&vmin2,0,1,writebuf); // 1
1971  myset(datatype,&vmax2,0,1,writebuf); // 1
1972  myset(datatype,&vmin22,0,1,writebuf); // 1
1973  myset(datatype,&vmax22,0,1,writebuf); // 1
1974  FTYPE vmin3=0,vmax3=0,vmax23=0,vmin23=0;
1975  vchar_rad(pr, &q, 3, ptrgeom, &vmax3, &vmin3, &vmax23, &vmin23, &ignorecourant);
1976  myset(datatype,&vmin3,0,1,writebuf); // 1
1977  myset(datatype,&vmax3,0,1,writebuf); // 1
1978  myset(datatype,&vmin23,0,1,writebuf); // 1
1979  myset(datatype,&vmax23,0,1,writebuf); // 1
1980 
1981  return (0);
1982 }
1983 
1984 
1985 
1986 
1987 
1988 
1989 
1990 
1991 
1992 
1993 
1995 int vpotdump(long dump_cnt)
1996 {
1997  MPI_Datatype datatype;
1998  int whichdump;
1999  char fileprefix[MAXFILENAME]={'\0'};
2000  char filesuffix[MAXFILENAME]={'\0'};
2001  char fileformat[MAXFILENAME]={'\0'};
2002 
2003 
2004  trifprintf("begin dumping vpotdump# %ld ... ",dump_cnt);
2005 
2006  whichdump=VPOTDUMPTYPE;
2007  datatype=MPI_FTYPE;
2008  strcpy(fileprefix,"dumps/vpotdump");
2009  strcpy(fileformat,"%04ld"); //atch adjust dump every substep
2010  strcpy(filesuffix,"");
2011 
2012  if(dump_gen(WRITEFILE,dump_cnt,binaryoutput,whichdump,datatype,fileprefix,fileformat,filesuffix,dump_header,vpotdump_content)>=1) return(1);
2013 
2014  trifprintf("end dumping vpotdump# %ld ... ",dump_cnt);
2015 
2016 
2017  return(0);
2018 
2019 }
2020 
2021 
2023 void set_vpotdump_content_dnumcolumns_dnumversion(int *numcolumnsvar, int *numversion)
2024 {
2025 
2026  if(DOVPOTDUMP){
2027  *numcolumnsvar=NUMVPOTDUMP;
2028  }
2029  else *numcolumnsvar=0;
2030 
2031  // Version number:
2032  *numversion=0;
2033 
2034 
2035 
2036 }
2037 
2039 int vpotdump_content(int i, int j, int k, MPI_Datatype datatype,void *writebuf)
2040 {
2041  int jj;
2042 
2043  // NOTEMARK: see also restart.c since this is added to restart
2044  for(jj=0;jj<NUMVPOTDUMP;jj++){
2045  myset(datatype,&GLOBALMACP1A0(vpotarraydump,jj,i,j,k),0,1,writebuf); // 1 each
2046  }
2047 
2048  return (0);
2049 }
2050 
2051 
2052 
2053 
2054 
2055 
2056 
2058 int failfloordudump(long dump_cnt)
2059 {
2060  MPI_Datatype datatype;
2061  int whichdump;
2062  char fileprefix[MAXFILENAME]={'\0'};
2063  char filesuffix[MAXFILENAME]={'\0'};
2064  char fileformat[MAXFILENAME]={'\0'};
2065 
2066 
2067  trifprintf("begin dumping failfloordudump# %ld ... ",dump_cnt);
2068 
2069  whichdump=FAILFLOORDUDUMPTYPE;
2070  datatype=MPI_FTYPE;
2071  strcpy(fileprefix,"dumps/failfloordudump");
2072  strcpy(fileformat,"%04ld"); //atch adjust dump every substep
2073  strcpy(filesuffix,"");
2074 
2075  if(dump_gen(WRITEFILE,dump_cnt,binaryoutput,whichdump,datatype,fileprefix,fileformat,filesuffix,dump_header,failfloordudump_content)>=1) return(1);
2076 
2077  trifprintf("end dumping failfloordudump# %ld ... ",dump_cnt);
2078 
2079 
2080  return(0);
2081 
2082 }
2083 
2084 
2086 void set_failfloordudump_content_dnumcolumns_dnumversion(int *numcolumnsvar, int *numversion)
2087 {
2088 
2089  if(DOFLOORDIAG){
2090  *numcolumnsvar=NPR; // normal dU from floor
2091  }
2092  else *numcolumnsvar=0;
2093 
2094  // Version number:
2095  *numversion=0;
2096 
2097 
2098 
2099 }
2100 
2102 int failfloordudump_content(int i, int j, int k, MPI_Datatype datatype,void *writebuf)
2103 {
2104  int pl;
2105 
2106  // NOTEMARK: see also restart.c since this is added to restart
2107  myset(datatype,&GLOBALMAC(failfloordu,i,j,k),0,NPR,writebuf); // NPR
2108 
2109  return (0);
2110 }
2111 
2112 
2113 
2114 
2116 int fluxsimpledump(long dump_cnt)
2117 {
2118  MPI_Datatype datatype;
2119  int whichdump;
2120  char fileprefix[MAXFILENAME]={'\0'};
2121  char filesuffix[MAXFILENAME]={'\0'};
2122  char fileformat[MAXFILENAME]={'\0'};
2123 
2124 
2125  trifprintf("begin dumping fluxsimpledump# %ld ... ",dump_cnt);
2126 
2127  whichdump=FLUXSIMPLEDUMPTYPE;
2128  datatype=MPI_FTYPE;
2129  strcpy(fileprefix,"dumps/fluxsimpledump");
2130  strcpy(fileformat,"%04ld"); //atch adjust dump every substep
2131  strcpy(filesuffix,"");
2132 
2133  if(dump_gen(WRITEFILE,dump_cnt,binaryoutput,whichdump,datatype,fileprefix,fileformat,filesuffix,dump_header,fluxsimpledump_content)>=1) return(1);
2134 
2135  trifprintf("end dumping fluxsimpledump# %ld ... ",dump_cnt);
2136 
2137 
2138  return(0);
2139 
2140 }
2141 
2142 
2144 void set_fluxsimpledump_content_dnumcolumns_dnumversion(int *numcolumnsvar, int *numversion)
2145 {
2146 
2147  if(DOFLOORDIAG && N1NOT1==1){
2148  *numcolumnsvar=NPR; // radial fluxes only
2149  }
2150  else *numcolumnsvar=0;
2151 
2152 
2153  if(FLUXDUMP==2) *numcolumnsvar+= (NUMFLUXESTOSAVE*NUMPHYSICALFLUXTERMS);
2154 
2155 
2156  // Version number:
2157  *numversion=0;
2158 
2159 
2160 
2161 }
2162 
2164 int fluxsimpledump_content(int i, int j, int k, MPI_Datatype datatype,void *writebuf)
2165 {
2166  int pl;
2167 
2168  myset(datatype,&GLOBALMACP0A1(F1,i,j,k,0),0,NPR,writebuf); // NPR
2169 
2170 
2171  FTYPE ftemp;
2172  if(FLUXDUMP==2){ // can reduce this alot (rho not in en or yflx's, etc.)
2173  int fluxterm;
2174  int pliter;
2175  for(fluxterm=0;fluxterm<NUMPHYSICALFLUXTERMS;fluxterm++){
2176  PLOOP(pliter,pl){
2177  if(FLUXESTOSAVEPL(pl)){
2178  ftemp=(GLOBALMACP0A1(fluxdump,i,j,k,fluxterm*NPR + pl)); // (float)
2179  myset(datatype,&ftemp,0,1,writebuf);
2180  }
2181  }
2182  }
2183  }
2184 
2185  return (0);
2186 }
2187 
2188 
2189 
2190 
2192 int dissmeasuredump(long dump_cnt)
2193 {
2194  MPI_Datatype datatype;
2195  int whichdump;
2196  char fileprefix[MAXFILENAME]={'\0'};
2197  char filesuffix[MAXFILENAME]={'\0'};
2198  char fileformat[MAXFILENAME]={'\0'};
2199 
2200 
2201  trifprintf("begin dumping dissmeasuredump# %ld ... ",dump_cnt);
2202 
2203  whichdump=DISSMEASUREDUMPTYPE;
2204  datatype=MPI_FTYPE;
2205  strcpy(fileprefix,"dumps/dissmeasuredump");
2206  strcpy(fileformat,"%04ld");
2207  strcpy(filesuffix,"");
2208 
2209  if(dump_gen(WRITEFILE,dump_cnt,binaryoutput,whichdump,datatype,fileprefix,fileformat,filesuffix,dump_header,dissmeasuredump_content)>=1) return(1);
2210 
2211  trifprintf("end dumping dissmeasuredump# %ld ... ",dump_cnt);
2212 
2213 
2214  return(0);
2215 
2216 }
2217 
2219 void set_dissmeasuredump_content_dnumcolumns_dnumversion(int *numcolumnsvar, int *numversion)
2220 {
2221 
2222  if(DODISSMEASURE && DODISSMEASUREDIAG){
2223  *numcolumnsvar=NSPECIAL+1; // dissmeasurepl and dissmeasure
2224  *numcolumnsvar+=3*2; // Fi for each direction and gas/rad
2225  }
2226  else *numcolumnsvar=0;
2227 
2228  // Version number:
2229  *numversion=0;
2230 
2231 }
2232 
2233 
2235 int dissmeasuredump_content(int i, int j, int k, MPI_Datatype datatype,void *writebuf)
2236 {
2237 
2238  // NOTEMARK: see also restart.c since this is added to restart
2239  myset(datatype,&GLOBALMAC(dissmeasurearray,i,j,k),0,NSPECIAL+1+3*2,writebuf);
2240 
2241  return (0);
2242 }
2243 
2244 
2245 
2248 int fakedump(long dump_cnt)// arg not used
2249 {
2250  MPI_Datatype datatype;
2251  int whichdump;
2252  char fileprefix[MAXFILENAME]={'\0'};
2253  char filesuffix[MAXFILENAME]={'\0'};
2254  char fileformat[MAXFILENAME]={'\0'};
2255 
2256 
2257  trifprintf("begin dumping fakedump");
2258 
2259  // return(0); // Kraken had problems with this when computing and dumping currents.
2260 
2261  whichdump=FAKEDUMPTYPE;
2262  datatype=MPI_FTYPE;
2263  strcpy(fileprefix,"dumps/fakedump");
2264  strcpy(fileformat,"%04ld"); //atch adjust dump every substep
2265  strcpy(filesuffix,"");
2266 
2267  if(dump_gen(WRITEFILE,-1,binaryoutput,whichdump,datatype,fileprefix,fileformat,filesuffix,fakedump_header,fakedump_content)>=1) return(1);
2268 
2269  trifprintf("end dumping fakedump");
2270 
2271 
2272  return(0);
2273 
2274 }
2275 
2276 
2277 // fake dump content number
2278 int fakedump_content(int i, int j, int k, MPI_Datatype datatype,void *writebuf)
2279 {
2280  int jj;
2281 
2282  // blank
2283 
2284  return (0);
2285 }
2286 
2287 
2288 // fake dump header
2289 int fakedump_header(int whichdump, int whichdumpversion, int numcolumnsvar,int bintxt, FILE *headerptr)
2290 {
2291  int fake;
2292 
2293  fake=0;
2294 
2295 
2296  if(bintxt==BINARYOUTPUT){
2297  fwrite(&fake,sizeof(int),1,headerptr);
2298  }
2299  else{
2300 #if(REALTYPE==DOUBLETYPE)
2301  fprintf(headerptr, "DONE: %d\n",fake);
2302 #elif(REALTYPE==LONGDOUBLETYPE)
2303  fprintf(headerptr, "DONE: %d\n",fake);
2304 #endif
2305  }
2306  fflush(headerptr);
2307  return(0);
2308 }
2309 
Generated on Fri May 20 2016 15:52:33 for HARM by   doxygen 1.8.3.1