OSHUN  beta
Arbitrary Order Spherical-Harmonic 1D-3P Vlasov-Fokker-Planck-Maxwell code
export.cpp
Go to the documentation of this file.
1 
8 // Standard libraries
9 #include <mpi.h>
10 #include <iostream>
11 #include <vector>
12 #include <valarray>
13 #include <complex>
14 #include <algorithm>
15 #include <fstream>
16 #include <iomanip>
17 #include <cstdlib>
18 #include <sstream>
19 #include <string>
20 #include <cstring>
21 
22 #include <math.h>
23 #include <map>
24 
25 #include <sys/stat.h>
26 #include <sys/types.h>
27 
28 // My libraries
29 #include "lib-array.h"
30 #include "lib-algorithms.h"
31 #include "H5Cpp.h"
32 #include "exprtk.hpp"
33 
34 // Declarations
35 #include "input.h"
36 #include "state.h"
37 #include "formulary.h"
38 #include "setup.h"
39 #include "nmethods.h"
40 #include "parallel.h"
41 #include "nmethods.h"
42 #include "export.h"
43 
44 
45 
46 //------------------------------------------------------------------------------
47 // Create a folder
48 //
49 // @param[in] _name The name of the folder
50 //
51 int Export_Files::Makefolder(string _name){
52 
53  mode_t _permissions(S_IRWXU | S_IRGRP | S_IXGRP | S_IROTH | S_IXOTH);
54  // char* foldername = new char [_name.length() + 1];
55 
56 
57  const char* foldername = _name.data();
58 
59 
60  // strcpy(foldername, _name.c_str());
61 
62  // vector<char> foldername(_name.c_str(), _name.c_str() + _name.size() + 1);
63 
64  int status(mkdir(foldername,_permissions));
65 
66  // if (status != 0) cout << "Warning: Folder "<< _name<< " exists\n";
67 
68  // delete[] foldername;
69 
70  return status;
71 }
72 //--------------------------------------------------------------
73 //---------------------------------------------------------------
74 void Export_Files::Folders(){
75 //--------------------------------------------------------------
76 // create the directory tree
77 //--------------------------------------------------------------
78 
79  if (Makefolder("RESTART") != 0) cout << "Warning: Folder 'RESTART' exists" << endl;
80 
81  if (Makefolder("OUTPUT") != 0) cout << "Warning: Folder 'OUTPUT' exists" << endl;
82 
83  if ( Input::List().o_EHist ) {
84  if ( Makefolder("OUTPUT/NUM") != 0) cout << "Warning: Folder 'OUTPUT/NUM' exists" << endl;
85  }
86 
87  if ( Input::List().o_Ex || Input::List().o_Ey || Input::List().o_Ez ||
88  Input::List().o_Bx || Input::List().o_By || Input::List().o_Bz ) {
89  if (Makefolder("OUTPUT/FLD") != 0)
90  cout<<"Warning: Folder 'OUTPUT/FLD' exists" << endl;
91 
92  if (Input::List().o_Ex) {
93  if (Makefolder("OUTPUT/FLD/Ex") != 0)
94  cout<<"Warning: Folder 'OUTPUT/FLD/EX' exists" << endl;
95  }
96  if (Input::List().o_Ey) {
97  if (Makefolder("OUTPUT/FLD/Ey") != 0)
98  cout<<"Warning: Folder 'OUTPUT/FLD/EY' exists" << endl;
99  }
100  if (Input::List().o_Ez) {
101  if (Makefolder("OUTPUT/FLD/Ez") != 0)
102  cout<<"Warning: Folder 'OUTPUT/FLD/EZ' exists" << endl;
103  }
104  if (Input::List().o_Bx) {
105  if (Makefolder("OUTPUT/FLD/Bx") != 0)
106  cout<<"Warning: Folder 'OUTPUT/FLD/BX' exists" << endl;
107  }
108  if (Input::List().o_By) {
109  if (Makefolder("OUTPUT/FLD/By") != 0)
110  cout<<"Warning: Folder 'OUTPUT/FLD/BY' exists" << endl;
111  }
112  if (Input::List().o_Bz) {
113  if (Makefolder("OUTPUT/FLD/Bz") != 0)
114  cout<<"Warning: Folder 'OUTPUT/FLD/BZ' exists" << endl;
115  }
116  }
117 
118  if ( Input::List().o_x1x2 || Input::List().o_pth ||
119  Input::List().o_G ||
120  Input::List().o_Px || Input::List().o_PxPx ||
121  Input::List().o_Py || Input::List().o_PxPy || Input::List().o_PyPy ||
122  Input::List().o_Pz || Input::List().o_PxPz || Input::List().o_PyPz || Input::List().o_PzPz ||
123  Input::List().o_Vx || Input::List().o_VxVx ||
124  Input::List().o_Vy || Input::List().o_VxVy || Input::List().o_VyVy ||
125  Input::List().o_VxVz || Input::List().o_VyVz || Input::List().o_VzVz ||
126  Input::List().o_Vsq || Input::List().o_Temperature || Input::List().o_Pressure ||
127  Input::List().o_Qx || Input::List().o_Qy || Input::List().o_Qz ||
128  Input::List().o_p1x1 ) {
129 
130  if (Makefolder("OUTPUT/MOM") != 0)
131  cout<<"Warning: Folder 'OUTPUT/MOM' exists" << endl;
132 
133 // Relativistic Energy - Momentum Tensor
134 // - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
135  if (Input::List().o_x1x2) {
136  if (Makefolder("OUTPUT/MOM/n") != 0)
137  cout<<"Warning: Folder 'OUTPUT/MOM/N' exists" << endl;
138  }
139  if (Input::List().o_pth) {
140  if (Makefolder("OUTPUT/MOM/Pth") != 0)
141  cout<<"Warning: Folder 'OUTPUT/MOM/Pth' exists" << endl;
142  }
143  if (Input::List().o_G) {
144  if (Makefolder("OUTPUT/MOM/Gam") != 0)
145  cout<<"Warning: Folder 'OUTPUT/MOM/Gam' exists" << endl;
146  }
147  if (Input::List().o_Px) {
148  if (Makefolder("OUTPUT/MOM/Px") != 0)
149  cout<<"Warning: Folder 'OUTPUT/MOM/Px' exists" << endl;
150  }
151  if (Input::List().o_PxPx) {
152  if (Makefolder("OUTPUT/MOM/PxPx") != 0)
153  cout<<"Warning: Folder 'OUTPUT/MOM/PxPx' exists" << endl;
154  }
155  if (Input::List().o_Py) {
156  if (Makefolder("OUTPUT/MOM/Py") != 0)
157  cout<<"Warning: Folder 'OUTPUT/MOM/Py' exists" << endl;
158  }
159  if (Input::List().o_PxPy) {
160  if (Makefolder("OUTPUT/MOM/PxPy") != 0)
161  cout<<"Warning: Folder 'OUTPUT/MOM/PxPy' exists" << endl;
162  }
163  if (Input::List().o_PyPy) {
164  if (Makefolder("OUTPUT/MOM/PyPy") != 0)
165  cout<<"Warning: Folder 'OUTPUT/MOM/PyPy' exists" << endl;
166  }
167  if (Input::List().o_Pz) {
168  if (Makefolder("OUTPUT/MOM/Pz") != 0)
169  cout<<"Warning: Folder 'OUTPUT/MOM/Pz' exists" << endl;
170  }
171  if (Input::List().o_PxPz) {
172  if (Makefolder("OUTPUT/MOM/PxPz") != 0)
173  cout<<"Warning: Folder 'OUTPUT/MOM/PxPz' exists" << endl;
174  }
175  if (Input::List().o_PyPz) {
176  if (Makefolder("OUTPUT/MOM/PyPz") != 0)
177  cout<<"Warning: Folder 'OUTPUT/MOM/PyPz' exists" << endl;
178  }
179  if (Input::List().o_PzPz) {
180  if (Makefolder("OUTPUT/MOM/PzPz") != 0)
181  cout<<"Warning: Folder 'OUTPUT/MOM/PzPz' exists" << endl;
182  }
183 
184 // Current
185 // - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
186  if (Input::List().o_Jx) {
187  if (Makefolder("OUTPUT/MOM/Jx") != 0)
188  cout<<"Warning: Folder 'OUTPUT/MOM/Jx' exists" << endl;
189  }
190  if (Input::List().o_Jy) {
191  if (Makefolder("OUTPUT/MOM/Jy") != 0)
192  cout<<"Warning: Folder 'OUTPUT/MOM/Jy' exists" << endl;
193  }
194  if (Input::List().o_Jz) {
195  if (Makefolder("OUTPUT/MOM/Jz") != 0)
196  cout<<"Warning: Folder 'OUTPUT/MOM/Jz' exists" << endl;
197  }
198 
199 
200 // Nonrelativistic output
201 // - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
202  if (Input::List().o_Vx) {
203  if (Makefolder("OUTPUT/MOM/Vx") != 0)
204  cout<<"Warning: Folder 'OUTPUT/MOM/Vx' exists" << endl;
205  }
206  if (Input::List().o_VxVx) {
207  if (Makefolder("OUTPUT/MOM/VxVx") != 0)
208  cout<<"Warning: Folder 'OUTPUT/MOM/VxVx' exists" << endl;
209  }
210  if (Input::List().o_Vy) {
211  if (Makefolder("OUTPUT/MOM/Vy") != 0)
212  cout<<"Warning: Folder 'OUTPUT/MOM/Vy' exists" << endl;
213  }
214  if (Input::List().o_VxVy) {
215  if (Makefolder("OUTPUT/MOM/VxVy") != 0)
216  cout<<"Warning: Folder 'OUTPUT/MOM/VxVy' exists" << endl;
217  }
218  if (Input::List().o_VyVy) {
219  if (Makefolder("OUTPUT/MOM/VyVy") != 0)
220  cout<<"Warning: Folder 'OUTPUT/MOM/VyVy' exists" << endl;
221  }
222  if (Input::List().o_VxVz) {
223  if (Makefolder("OUTPUT/MOM/VxVz") != 0)
224  cout<<"Warning: Folder 'OUTPUT/MOM/VxVz' exists" << endl;
225  }
226  if (Input::List().o_VyVz) {
227  if (Makefolder("OUTPUT/MOM/VyVz") != 0)
228  cout<<"Warning: Folder 'OUTPUT/MOM/VyVz' exists" << endl;
229  }
230  if (Input::List().o_VzVz) {
231  if (Makefolder("OUTPUT/MOM/VzVz") != 0)
232  cout<<"Warning: Folder 'OUTPUT/MOM/VzVz' exists" << endl;
233  }
234  if (Input::List().o_Vsq) {
235  if (Makefolder("OUTPUT/MOM/Vsq") != 0)
236  cout<<"Warning: Folder 'OUTPUT/MOM/Vsq' exists" << endl;
237  }
238  if (Input::List().o_Temperature) {
239 // if (Makefolder("OUTPUT/MOM/T_eV") != 0)
240 // cout<<"Warning: Folder 'OUTPUT/MOM/T_eV' exists" << endl;
241  if (Makefolder("OUTPUT/MOM/T") != 0)
242  cout<<"Warning: Folder 'OUTPUT/MOM/T' exists" << endl;
243  }
244  if (Input::List().o_Pressure) {
245  if (Makefolder("OUTPUT/MOM/P_Mbar") != 0)
246  cout<<"Warning: Folder 'OUTPUT/MOM/P_Mbar' exists" << endl;
247  }
248  if (Input::List().o_ND) {
249  if (Makefolder("OUTPUT/MOM/ND") != 0)
250  cout<<"Warning: Folder 'OUTPUT/MOM/ND' exists" << endl;
251  }
252  if (Input::List().o_Nu) {
253  if (Makefolder("OUTPUT/MOM/Nu") != 0)
254  cout<<"Warning: Folder 'OUTPUT/MOM/Nu' exists" << endl;
255  }
256  if (Input::List().o_Qx) {
257  if (Makefolder("OUTPUT/MOM/Qx") != 0)
258  cout<<"Warning: Folder 'OUTPUT/MOM/Qx' exists" << endl;
259  }
260  if (Input::List().o_Qy) {
261  if (Makefolder("OUTPUT/MOM/Qy") != 0)
262  cout<<"Warning: Folder 'OUTPUT/MOM/Qy' exists" << endl;
263  }
264  if (Input::List().o_Qz) {
265  if (Makefolder("OUTPUT/MOM/Qz") != 0)
266  cout<<"Warning: Folder 'OUTPUT/MOM/Qz' exists" << endl;
267  }
268  if (Input::List().o_vNx) {
269  if (Makefolder("OUTPUT/MOM/vNx") != 0)
270  cout<<"Warning: Folder 'OUTPUT/MOM/vNx' exists" << endl;
271  }
272  if (Input::List().o_vNy) {
273  if (Makefolder("OUTPUT/MOM/vNy") != 0)
274  cout<<"Warning: Folder 'OUTPUT/MOM/vNy' exists" << endl;
275  }
276  if (Input::List().o_vNz) {
277  if (Makefolder("OUTPUT/MOM/vNz") != 0)
278  cout<<"Warning: Folder 'OUTPUT/MOM/vNz' exists" << endl;
279  }
280  if (Input::List().o_Ux) {
281  if (Makefolder("OUTPUT/MOM/Ux") != 0)
282  cout<<"Warning: Folder 'OUTPUT/MOM/Ux' exists" << endl;
283  }
284  if (Input::List().o_Uy) {
285  if (Makefolder("OUTPUT/MOM/Uy") != 0)
286  cout<<"Warning: Folder 'OUTPUT/MOM/Uy' exists" << endl;
287  }
288  if (Input::List().o_Uz) {
289  if (Makefolder("OUTPUT/MOM/Uz") != 0)
290  cout<<"Warning: Folder 'OUTPUT/MOM/Ux' exists" << endl;
291  }
292  if (Input::List().o_Z) {
293  if (Makefolder("OUTPUT/MOM/Z") != 0)
294  cout<<"Warning: Folder 'OUTPUT/MOM/Z' exists" << endl;
295  }
296  if (Input::List().o_ni) {
297  if (Makefolder("OUTPUT/MOM/ni") != 0)
298  cout<<"Warning: Folder 'OUTPUT/MOM/ni' exists" << endl;
299  }
300  if (Input::List().o_Ti) {
301  if (Makefolder("OUTPUT/MOM/Ti") != 0)
302  cout<<"Warning: Folder 'OUTPUT/MOM/Ti' exists" << endl;
303  }
304 //><-><-><-><-><-><-><-><-><-><-><-><-><-><-><-><-><-><-><-><-><-><-><-><-><-><-><
305 // TWO STREAM STUFF
306 //><-><-><-><-><-><-><-><-><-><-><-><-><-><-><-><-><-><-><-><-><-><-><-><-><-><-><
307 // if (Input::List().o_p1x1) {
308 // if (Makefolder("OUTPUT/MOM/p1x1") != 0)
309 // cout<<"Warning: Folder 'OUTPUT/MOM/p1x1' exists\n";
310 // }
311 //><-><-><-><-><-><-><-><-><-><-><-><-><-><-><-><-><-><-><-><-><-><-><-><-><-><-><
312 
313  }
314 
315  if ( Input::List().o_f0x1 || Input::List().o_p1x1
316  || Input::List().o_f10x1 || Input::List().o_f11x1 || Input::List().o_f20x1
317  || Input::List().o_fl0x1) {
318  if (Makefolder("OUTPUT/DISTR") != 0)
319  cout<<"Warning: Folder 'OUTPUT/DISTR' exists" << endl;
320 
321  }
322 
323 
324 }
325 //--------------------------------------------------------------
326 //**************************************************************
327 
328 //--------------------------------------------------------------
329 // List of the acceptable Tags
330 // Constructor
331 //------------------------------------------------------------------------------
332 // List of the acceptable Tags Constructor
333 //
334 // @param[in] species The species
335 //
336 Export_Files::DefaultTags::DefaultTags(size_t species){
337 
338 // Time
339  time.push_back( "Time_cgs");
340  time.push_back( "Time_si" );
341  time.push_back( "Time_fs" );
342  time.push_back( "Time_ps" );
343  time.push_back( "Time" ); // Default tag
344 
345 // Space
346  space.push_back( "Space_cgs");
347  space.push_back( "Space_si" );
348  space.push_back( "Space");
349 
350 // Fields
351  fld.push_back( "Ex" );
352  fld.push_back( "Ex_cgs" );
353  fld.push_back( "Ex_si" );
354  fld.push_back( "Ey" );
355  fld.push_back( "Ey_cgs" );
356  fld.push_back( "Ey_si" );
357  fld.push_back( "Ez" );
358  fld.push_back( "Ez_cgs" );
359  fld.push_back( "Ez_si" );
360  fld.push_back( "Bx" );
361  fld.push_back( "Bx_cgs" );
362  fld.push_back( "Bx_si" );
363  fld.push_back( "By" );
364  fld.push_back( "By_cgs" );
365  fld.push_back( "By_si" );
366  fld.push_back( "Bz" );
367  fld.push_back( "Bz_cgs" );
368  fld.push_back( "Bz_si" );
369 
370 // Moments
371  mom.push_back( "P" );
372  mom.push_back( "P_cgs" );
373  mom.push_back( "P_si" );
374  mom.push_back( "P_Mbar" );
375  mom.push_back( "T" );
376  mom.push_back( "T_cgs" );
377  mom.push_back( "T_si" );
378  mom.push_back( "T_eV" );
379  mom.push_back( "n" );
380  mom.push_back( "n_cgs" );
381  mom.push_back( "n_si" );
382  mom.push_back( "Qx" );
383  mom.push_back( "Qy" );
384  mom.push_back( "Qz" );
385  mom.push_back( "vNx" );
386  mom.push_back( "vNy" );
387  mom.push_back( "vNz" );
388  mom.push_back( "Jx" );
389  mom.push_back( "Jx_cgs" );
390  mom.push_back( "Jx_si" );
391  mom.push_back( "Jy" );
392  mom.push_back( "Jy_cgs" );
393  mom.push_back( "Jy_si" );
394  mom.push_back( "Jz" );
395  mom.push_back( "Jz_cgs" );
396  mom.push_back( "Jz_si" );
397 
398  mom.push_back( "Ux" );
399  mom.push_back( "Uy" );
400  mom.push_back( "Uz" );
401  mom.push_back( "Z" );
402  mom.push_back( "ni" );
403  mom.push_back( "Ti" );
404 
405 
406 
407 
408 // p-x
409  for (size_t s(0); s < species; ++s) {
410  pvsx.push_back( "px-x");
411  pvsx.push_back( "py-x");
412  pvsx.push_back( "pz-x");
413  }
414 
415 // f-x
416  for (size_t s(0); s < species; ++s) {
417  fvsx.push_back( "f0-x");//+stringify(s) );
418  fvsx.push_back( "f10-x");
419  fvsx.push_back( "f11-x");
420  fvsx.push_back( "f20-x");
421  fvsx.push_back( "fl0-x");
422  }
423 
424 // p-p-x
425  for (size_t s(0); s < species; ++s) {
426  pvspvsx.push_back( "pxpy-x");
427  pvspvsx.push_back( "pypz-x");
428  pvspvsx.push_back( "pxpz-x");
429  }
430 
431 }
432 
433 
434 
435 //--------------------------------------------------------------
436 
437 
438 // Definition of the output axis
439 // Constructor
440 Export_Files::oAxis::oAxis() : label(""), units(""), min(0.0), max(1.0), sz(3) {}
441 Export_Files::oAxis::oAxis( const float _m, const float _M,
442  const size_t _sz)
443  : label(""), units(""), min(_m), max(_M), sz(_sz) {}
444 Export_Files::oAxis::oAxis(const string _l, const string _u, const float _m, const float _M,
445  const size_t _sz)
446  : label(_l), units(_u), min(_m), max(_M), sz(_sz) {}
447 // Copy constructor
448 Export_Files::oAxis::oAxis(const oAxis& other) {
449  label = other.label;
450  units = other.units;
451  min = other.min;
452  max = other.max;
453  sz = other.sz;
454 }
455 //--------------------------------------------------------------
456 
457 //--------------------------------------------------------------
458 // 1D header constructor
459 Export_Files::Header::Header(oAxis _x,
460  string _Ql, float _Qc,
461  string _tl, string _tu, float _tc,
462  string _oD)
463  : title(_Ql), titleC(_Qc),
464  time(_tl), timeU(_tu), timeC(_tc),
465  oDir(_oD) {
466  xyz_axis.push_back(_x);
467 }
468 
469 // 2D header constructor
470 Export_Files::Header::Header(oAxis _x, oAxis _y,
471  string _Ql, float _Qc,
472  string _tl, string _tu, float _tc,
473  string _oD)
474  : title(_Ql), time(_tl), timeU(_tu),
475  titleC(_Qc), timeC(_tc),
476  oDir(_oD) {
477  xyz_axis.push_back(_x);
478  xyz_axis.push_back(_y);
479 }
480 
481 // 3D header constructor
482 Export_Files::Header::Header(oAxis _x, oAxis _y, oAxis _z,
483  string _Ql, float _Qc,
484  string _tl, string _tu, float _tc,
485  string _oD)
486  : title(_Ql), time(_tl), timeU(_tu),
487  titleC(_Qc), timeC(_tc),
488  oDir(_oD) {
489  xyz_axis.push_back(_x);
490  xyz_axis.push_back(_y);
491  xyz_axis.push_back(_z);
492 }
493 
494 // xD header constructor
495 Export_Files::Header::Header(vector< oAxis > _xyz,
496  string _Ql, float _Qc,
497  string _tl, string _tu, float _tc,
498  string _oD)
499  : xyz_axis(_xyz), title(_Ql), time(_tl), timeU(_tu),
500  titleC(_Qc), timeC(_tc),
501  oDir(_oD) {}
502 
503 // number of header dimensions
504 size_t Export_Files::Header::dim() {
505  return xyz_axis.size();
506 }
507 //- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
508 
509 valarray<float> Export_Files::Header::axis(const size_t i) {
510 // return Algorithms::MakeAxis(xyz_axis[i].min, xyz_axis[i].max, xyz_axis[i].sz);
511  return Algorithms::MakeCAxis(xyz_axis[i].min, xyz_axis[i].max, xyz_axis[i].sz);
512 }
513 string Export_Files::Header::label(const size_t i) {
514  return xyz_axis[i].label;
515 }
516 string Export_Files::Header::units(const size_t i) {
517  return xyz_axis[i].units;
518 }
519 
520 string Export_Files::Header::Title_label() { return title; }
521 float Export_Files::Header::Title_conv() { return titleC; }
522 string Export_Files::Header::Time_label() { return time; }
523 float Export_Files::Header::Time_conv() { return timeC; }
524 string Export_Files::Header::Directory() { return oDir; }
525 //--------------------------------------------------------------
526 
527 
528 //--------------------------------------------------------------
529 // Constructor of the export facility for data structures
530 Export_Files::Xport::Xport(const Algorithms::AxisBundle<double>& _axis,
531  const vector< string > oTags,
532  string homedir){
533 
534  size_t species(_axis.pdim());
535  DefaultTags dTags(species);
536 
537  vector< oAxis > xyz, pxyz, imre, pr;
538  xyz.push_back( Export_Files::oAxis(_axis.xgmin(0), _axis.xgmax(0), _axis.Nxg(0)) );
539 
540  for (size_t s(0); s < species; ++s) {
541  pr.push_back( Export_Files::oAxis(_axis.pmin(s), _axis.pmax(s), _axis.Np(s)));
542  pxyz.push_back( Export_Files::oAxis( (-1.0)*float(_axis.pmax(s)), float(_axis.pmax(s)), _axis.Npx(s)) );
543  pxyz.push_back( Export_Files::oAxis( (-1.0)*float(_axis.pmax(s)), float(_axis.pmax(s)), _axis.Npy(s)) );
544  pxyz.push_back( Export_Files::oAxis( (-1.0)*float(_axis.pmax(s)), float(_axis.pmax(s)), _axis.Npz(s)) );
545  }
546  imre.push_back( Export_Files::oAxis(0,1,2) );
547 // Time
548  size_t tloc(0); // Find the location of the right tag
549  while ( ( tloc < dTags.time.size()-1 ) &&
550  ( find(oTags.begin(),oTags.end(), dTags.time[tloc]) == oTags.end() ) ) {
551  ++tloc;
552  }
553  string tlabel = "t[" +formulary().Label(dTags.time[tloc])+"]";
554  string tunits = formulary().Label(dTags.time[tloc]);
555  float tconv = formulary().Uconv(dTags.time[tloc]);
556 
557 // xyz Axis
558  size_t xloc(0); // Find the location of the right tag
559  while ( ( xloc < dTags.space.size()-1 ) &&
560  ( find(oTags.begin(),oTags.end(), dTags.space[xloc]) == oTags.end() ) ) {
561  ++xloc;
562  }
563  xyz[0].label = "x["+ formulary().Label(dTags.space[xloc]) +"]";
564  if ( xyz.size() > 1 ) xyz[1].label = "y["+ formulary().Label(dTags.space[xloc]) +"]";
565  if ( xyz.size() > 2 ) xyz[2].label = "z["+ formulary().Label(dTags.space[xloc]) +"]";
566  xyz[0].units = formulary().Label(dTags.space[xloc]);
567  if ( xyz.size() > 1 ) xyz[1].units = formulary().Label(dTags.space[xloc]);
568  if ( xyz.size() > 2 ) xyz[2].units = formulary().Label(dTags.space[xloc]);
569 
570 
571 // pxyz Axis
572  for (size_t i(0); i < species; ++i) {
573  pr[i].label = "p";
574  pr[i].units = "m_e c";
575 
576  pxyz[i].label = "px[mc]";
577  pxyz[i].units = "m_e c";
578  // if ( pxyz.size()/species > 1 )
579  pxyz[1+i].label = "py[mc]";
580  pxyz[1+i].units = "m_e c";
581  // if ( pxyz.size()/species > 2 )
582  pxyz[2+i].label = "pz[mc]";
583  pxyz[2+i].units = "m_e c";
584 
585  // if ( pxyz.size()/species > 1 )
586  // if ( pxyz.size()/species > 2 ) pxyz[2*species+i].units = "m_e c";
587  }
588 
589 // Tags for Fields -->
590  for (size_t i(0); i < dTags.fld.size(); ++i) {
591 
592  // If this tag is an output tag
593  if ( find(oTags.begin(),oTags.end(), dTags.fld[i]) != oTags.end() ) {
594 
595  string nounits = dTags.fld[i].substr(0, dTags.fld[i].find("_"));
596  string folder = homedir + "OUTPUT/FLD/" + nounits + "/";
597  // Generate a header file for this tag
598  Hdr[dTags.fld[i]] = Header(xyz,
599  nounits+"["+formulary().Label(dTags.fld[i])+"]",
600  formulary().Uconv(dTags.fld[i]),
601  tlabel, tunits, tconv, folder);
602  }
603  } // <--
604 
605 // Tags for Moments -->
606  for (size_t i(0); i < dTags.mom.size(); ++i) {
607 
608  // If this tag is an output tag
609  if ( find(oTags.begin(),oTags.end(), dTags.mom[i]) != oTags.end() ) {
610 
611  string nounits = dTags.mom[i].substr(0, dTags.mom[i].find("_"));
612  string folder = homedir + "OUTPUT/MOM/" + nounits + "/";
613  // Generate a header file for this tag
614  Hdr[dTags.mom[i]] = Header(xyz,
615  nounits+"["+formulary().Label(dTags.mom[i])+"]",
616  formulary().Uconv(dTags.mom[i]),
617  tlabel, tunits, tconv, folder);
618  }
619  } // <--
620 
621 // Tags for p-x -->
622  for (size_t i(0); i < dTags.pvsx.size(); ++i) {
623 
624  // for (size_t k(0); k < oTags.size(); ++k)
625  // std::cout << " \n \n k = " << oTags[k];
626 
627  // If this tag is an output tag
628  if ( find(oTags.begin(),oTags.end(), dTags.pvsx[i]) != oTags.end() ) {
629 
630  string folder = homedir + "OUTPUT/DISTR/" + dTags.pvsx[i] + "/";
631  Makefolder(folder);
632 
633 
634 // Generate a header file for this tag
635 // For each 9 you have a different species
636  Hdr[dTags.pvsx[i]] = Header( pr[i/3], xyz[0],
637  "f"+stringify(i/3), 1.0, tlabel, tunits, tconv, folder);
638  // std::cout << " \n k = " << dTags.pvsx[i] << "\n";
639  // std::cout << " \n \n 11 ";
640  }
641  } //<--
642 
643 // Tags for f-x -->
644  for (size_t i(0); i < dTags.fvsx.size(); ++i) {
645 
646  // If this tag is an output tag
647  if ( find(oTags.begin(),oTags.end(), dTags.fvsx[i]) != oTags.end() ) {
648 
649  string folder = homedir + "OUTPUT/DISTR/" + dTags.fvsx[i] + "/";
650  Makefolder(folder);
651 
652 // Generate a header file for this tag
653  Hdr[dTags.fvsx[i]] = Header( pr[i/5], xyz[0], imre[0],
654  "f", 1.0, tlabel, tunits, tconv, folder);
655 
656  }
657  } //<--
658 
659  // Tags for p-p-x -->
660  for (size_t i(0); i < dTags.pvspvsx.size(); ++i) {
661 
662  // If this tag is an output tag
663  if ( find(oTags.begin(),oTags.end(), dTags.pvspvsx[i]) != oTags.end() ) {
664 
665  string folder = homedir + "OUTPUT/DISTR/" + dTags.pvspvsx[i] + "/";
666  Makefolder(folder);
667 
668 // Generate a header file for this tag
669  Hdr[dTags.pvspvsx[i]] = Header( pxyz[i/1+species],pxyz[(i+1)/1+species], xyz[0],
670  "f"+stringify(i/1), 1.0, tlabel, tunits, tconv, folder);
671 
672  // std::cout << " \n \n 11 ";
673  }
674  } //<--
675 }
676 //--------------------------------------------------------------
677 
678 
679 //--------------------------------------------------------------
680 //--------------------------------------------------------------
681 // Adjust H5 filenames with zeros to reach some prescribed length.
682 // Add the H5 filename extension.
683 string Export_Files::Xport::oH5Fextension(size_t step, int species){
684 
685  stringstream sFilename;
686 
687  if(species >= 0) sFilename << "_s" << species;
688 
689  sFilename << "_";
690 
691  // Number of zeros to add to the filename
692  int Nzeros(ofconventions::ofile_digits - stringify(step).length());
693  while (Nzeros-- > 0) {
694  sFilename << "0";
695  }
696 
697  sFilename << step << ofconventions::h5file_extension;
698 
699  return sFilename.str();
700 }
701 //--------------------------------------------------------------
702 //**************************************************************
703 
704 
705 
706 //**************************************************************
707 //--------------------------------------------------------------
708 Export_Files::Restart_Facility::Restart_Facility(const int rank, string homedir) {
709  hdir = homedir;
710 
711  if (!rank) Makefolder(hdir+"RESTART/");
712 
713 // if (Makefolder(hdir+"RESTART/") != 0) cout<<"Warning: Folder "<< hdir+"RESTART/"<<" exists\n";
714 }
715 
716 //--------------------------------------------------------------
717 // Read restart file
718 void Export_Files::Restart_Facility::Read(const int rank, const size_t re_step, State1D& Y) {
719 
720 // Generate filename
721  string filename(hdir+"RESTART/re_1D_");
722  filename.append(rFextension(rank,re_step));
723 
724 // Open file
725  ifstream fin(filename.c_str(), ios::binary);
726 
727  if (fin)
728  {
729  // Read distribution functions
730  for(size_t s(0); s < Y.Species(); ++s) {
731  for(size_t nh(0); nh < Y.DF(s).dim(); ++nh) {
732  for(size_t i(0); i < (Y.DF(s))(nh).dim(); ++i) {
733  fin.read((char *)&(Y.DF(s))(nh)(i), sizeof((Y.DF(s))(nh)(i)));
734  }
735  }
736  }
737  }
738  else {
739 
740  if (!rank) std::cout << "\n\n ERROR :: No files to read! \n\n";
741  exit(1);
742 
743  }
744 
745 // Read Ex
746  for(size_t i(0); i < Y.EMF().Ex().numx(); ++i){
747  fin.read((char *)(&Y.EMF().Ex()(i)), sizeof(Y.EMF().Ex()(i)));
748  }
749 
750  fin.close();
751 }
752 //--------------------------------------------------------------
753 
754 //--------------------------------------------------------------
755 // Write restart file
756 void Export_Files::Restart_Facility::Write(const int rank, const size_t re_step, State1D& Y) {
757 
758 // Generate filename
759  string filename(hdir+"RESTART/re_1D_");
760  filename.append(rFextension(rank,re_step));
761 
762 // Open file
763  ofstream fout(filename.c_str(), ios::binary);
764 
765 // Write distribution functions
766  for(size_t s(0); s < Y.Species(); ++s) {
767  for (size_t s(0); s < Y.Species(); ++s) {
768  for (size_t nh(0); nh < Y.DF(s).dim(); ++nh) {
769  for (size_t i(0); i < (Y.DF(s))(nh).dim(); ++i) {
770  fout.write((char *) &(Y.DF(s))(nh)(i), sizeof((Y.DF(s))(nh)(i)));
771  }
772  }
773  }
774  }
775 
776 // Write Ex
777  for(size_t i(0); i < Y.EMF().Ex().numx(); ++i) {
778  fout.write((char *)(&Y.EMF().Ex()(i)), sizeof(Y.EMF().Ex()(i)));
779  }
780 
781  fout.flush();
782  fout.close();
783 }
784 //--------------------------------------------------------------
785 
786 //--------------------------------------------------------------
787 // Adjust filenames with zeros to reach some prescribed length.
788 // Add the filename extension.
789 string Export_Files::Restart_Facility::rFextension(const int rank, const size_t rstep){
790  stringstream sFilename;
791 
792  // Number of zeros to add to the filename
793  int Nzeros(ofconventions::rank_digits - stringify(rank).length());
794  while (Nzeros-- > 0) {
795  sFilename << "0";
796  }
797 
798  sFilename << rank << "_";
799 
800  // Number of zeros to add to the filename
801  Nzeros = ofconventions::rfile_digits - stringify(rstep).length();
802  while (Nzeros-- > 0) {
803  sFilename << "0";
804  }
805 
806  sFilename << rstep << ofconventions::rfile_extension;
807 
808  return sFilename.str();
809 }
810 //--------------------------------------------------------------
811 //**************************************************************
812 
813 
814 //**************************************************************
815 //--------------------------------------------------------------
816 Output_Data::PLegendre1D::PLegendre1D( size_t Nl, size_t Nm, size_t Np,
817  float pmin, float pmax, size_t Npx ) {
818 
819  size_t sz(0);
820 
821  valarray<float> p( Algorithms::MakeAxis(pmin, pmax, Np)),
822  px(Algorithms::MakeAxis(float(-1.0)*pmax, pmax, Npx));
823 
824 // Generate the structure to save the polynomials
825  sz = ((Nm+1)*(2*Nl-Nm+2))/2;
826 //
827  plegendre = new vector< Array2D<float>>(sz,Array2D<float>(Npx,Np)) ;
828  // for (size_t m(0); m < Nm+1; ++m) {
829  // for (size_t l(m); l < Nl+1; ++l) {
830  // (*plegendre).push_back( Array2D<float>(Npx,Np) );
831  // ++sz;
832  // }
833  // }
834  size_t k(0);
835 // Generate the polynomial values for each cos(theta) = px/p
836  for (size_t j(0); j < p.size(); ++j) {
837  float invp(1.0/p[j]);
838  for (size_t i(0); i < px.size(); ++i) {
839 
840 // For given px/p generate all the polynomial values ...
841 // valarray<float> vL( Algorithms::Legendre( px[i]*invp, Nl, Nm ) );
842  Array2D<float> vL( Algorithms::Legendre( px[i]*invp, Nl, Nm ) );
843  // std::cout << "\n x = " << px[i]*invp << ", sqrt(1-x^2) = " << sqrt(1.0-px[i]*invp*px[i]*invp) ;
844 // ... and save them
845  // std::cout << "dim = " << vL.dim() << "\n";
846  // std::cout << "dim1:" << vL.dim1() << "\n";
847  // std::cout << "dim2:" << vL.dim1() << "\n\n";
848  for (size_t l(0); l < Nl+1; ++l){
849  for (size_t m=0; m<((Nm<l)?Nm:l)+1; ++m){
850  // for (size_t k(0); k < sz; ++k) {
851  k = ((l < Nm+1)?((l*(l+1))/2+m):(l*(Nm+1)-(Nm*(Nm+1))/2 + m));
852 
853  (*plegendre)[k](i,j) = vL(l,m);
854 
855  std::cout << "\n LP(" << l << "," << m << "," << k << ") = " << vL(l,m);
856 
857  }
858  // (*plegendre)[k](i,j) = vL(k);
859  }
860 
861  // exit(1);
862  }
863  }
864 
865 }
866 //- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
867 
868 Output_Data::PLegendre1D::PLegendre1D( const PLegendre1D& other ) {
869 
870 // Generate the structure to save the polynomials
871  plegendre = new vector< Array2D<float> > ;
872  for (size_t i(0); i < other.dim(); ++i) {
873  (*plegendre).push_back( other(i) );
874  }
875 }
876 //- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
877 
878 Output_Data::PLegendre1D::~PLegendre1D(){
879  delete plegendre;
880 }
881 //--------------------------------------------------------------
882 //**************************************************************
883 
884 //**************************************************************
885 //--------------------------------------------------------------
886 Output_Data::PLegendre2D::PLegendre2D( size_t Nl, size_t Nm, size_t Np,
887  float pmin, float pmax, size_t Npx, size_t Npy ) {
888 
889  size_t sz(0);
890 
891  valarray<float> p( Algorithms::MakeAxis(pmin, pmax, Np)),
892  px(Algorithms::MakeAxis(float(-1.0)*pmax, pmax, Npx)),
893  py(Algorithms::MakeAxis(float(-1.0)*pmax, pmax, Npy));
894 
895 // Generate the structure to save the polynomials
896  sz = ((Nm+1)*(2*Nl-Nm+2))/2;
897 //
898  plegendre = new vector< Array2D<float>>(sz,Array2D<float>(Npx,Npy)) ;
899  // for (size_t m(0); m < Nm+1; ++m) {
900  // for (size_t l(m); l < Nl+1; ++l) {
901  // (*plegendre).push_back( Array2D<float>(Npx,Np) );
902  // ++sz;
903  // }
904  // }
905  size_t k(0);
906  // float invpmax(1.0/pmax);
907 // Generate the polynomial values for each cos(theta) = px/p
908  for (size_t j(0); j < py.size(); ++j) {
909 
910  for (size_t i(0); i < px.size(); ++i) {
911  float invp(sqrt(py[j]*py[j]+px[i]+px[i]));
912  if (invp > 0.0 ) invp=1.0/invp;
913  else if (invp > pmax) invp *= 2;
914  else invp = 0.0;
915 // For given px/p generate all the polynomial values ...
916 // valarray<float> vL( Algorithms::Legendre( px[i]*invp, Nl, Nm ) );
917 
918  Array2D<float> vL( Algorithms::Legendre( px[i]*invp, Nl, Nm ) );
919  // std::cout << "\n x = " << px[i]*invp << ", sqrt(1-x^2) = " << sqrt(1.0-px[i]*invp*px[i]*invp) ;
920  // ... and save them
921  // std::cout << "dim = " << vL.dim() << "\n";
922  // std::cout << "dim1:" << vL.dim1() << "\n";
923  // std::cout << "dim2:" << vL.dim1() << "\n\n";
924  for (size_t l(0); l < Nl+1; ++l){
925  for (size_t m(0); m<((Nm<l)?Nm:l)+1; ++m){
926  // for (size_t k(0); k < sz; ++k) {
927  k = ((l < Nm+1)?((l*(l+1))/2+m):(l*(Nm+1)-(Nm*(Nm+1))/2 + m));
928 
929  (*plegendre)[k](i,j) = vL(l,m);
930  // if (invp == 0.0) (*plegendre)[k](i,j) = 0.0;
931 
932  // std::cout << "\n LP(" << l << "," << m << "," << k << ") = " << vL(l,m);
933 
934  }
935  // (*plegendre)[k](i,j) = vL(k);
936  }
937  }
938  // exit(1);
939  }
940 
941 }
942 //- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
943 
944 Output_Data::PLegendre2D::PLegendre2D( const PLegendre2D& other ) {
945 
946 // Generate the structure to save the polynomials
947  plegendre = new vector< Array2D<float> > ;
948  for (size_t i(0); i < other.dim(); ++i) {
949  (*plegendre).push_back( other(i) );
950  }
951 }
952 //- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
953 
954 Output_Data::PLegendre2D::~PLegendre2D(){
955  delete plegendre;
956 }
957 //--------------------------------------------------------------
958 //**************************************************************
959 
960 //**************************************************************
961 //--------------------------------------------------------------
962 Output_Data::p1x1_1D::p1x1_1D( const Grid_Info& _G) {
963 
964 // Generate the required structures
965  for (size_t s(0); s < _G.axis.pdim(); ++s) {
966  pmax.push_back( static_cast<float>(_G.axis.pmax(s)) );
967  pmin.push_back( static_cast<float>(-1.0*_G.axis.pmax(s)) );
968 // Pl.push_back( PLegendre1D( _G.l0[s], _G.m0[s], _G.axis.Np(s), static_cast<float>(_G.axis.pmin(s)),
969 // static_cast<float>(_G.axis.pmax(s)), _G.axis.Npx(s) ) );
970 //
971 // polarR.push_back( Array2D<float>( _G.axis.Npx(s), _G.axis.Np(s) ) );
972 
973  p1x1.push_back( valarray<float>( _G.axis.Npx(s)) );
974 
975  }
976 
977 /*// Calculate a 2D array for the polar radius for each px and pr
978  for (size_t s(0); s < _G.axis.pdim(); ++s) {
979  valarray<float> p( Algorithms::MakeAxis( pmin[s], pmax[s], _G.axis.Np(s) )),
980  px(Algorithms::MakeAxis( float(-1.0)* pmax[s], pmax[s], _G.axis.Npx(s) ));
981 // --->
982  for (size_t i(0); i < _G.axis.Npx(s); ++i){
983  for (size_t j(0); j < _G.axis.Np(s); ++j){
984  float polarR_sq = p[j] * p[j] - px[i] * px[i];
985  if (polarR_sq < 0.0) {
986  polarR[s](i,j) = -1.0;
987  }
988  else {
989  polarR[s](i,j) = sqrt(abs(polarR_sq));
990  }
991  // std::cout << "\n polarR[s](" << px[i] << ',' << p[j] << ") = " << polarR[s](i,j) << "\n";
992  }
993  }
994 // <---
995  }*/
996 // p1x1 = new valarray<float>(Npx);
997 }
998 //- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
999 
1000 Output_Data::p1x1_1D::p1x1_1D( const p1x1_1D& other) {
1001 
1002  for (size_t s(0); s < other.Species(); ++s) {
1003  pmin.push_back( other.Pmin(s) );
1004  pmax.push_back( other.Pmax(s) );
1005 // Pl.push_back( other.PL(s) );
1006 // polarR.push_back( other.PolarR(s) );
1007  p1x1.push_back( other.p1_x1(s) );
1008  }
1009 }
1010 //- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
1011 
1012 Output_Data::p1x1_1D::~p1x1_1D(){
1013  // delete p1x1;
1014 }
1015 //--------------------------------------------------------------
1016 
1017 //--------------------------------------------------------------
1018 // Turn the Distribution function at some spatial location x0
1019 // into a cartesian grid.
1020 valarray<float> Output_Data::p1x1_1D::operator()(DistFunc1D& df, size_t x0, size_t s) {
1021 
1022  p1x1[s] = 0.0; // this is a valarray<float>
1023  vector<complex<double> > dummyvec;
1024 
1025 /* float p_im1(0.0), p_ip1(0.0);
1026  float integrant_low(0.0),
1027  integrant_high(0.0);
1028 
1029 // Calculate the integral for each harmonic separately
1030 // for (size_t l(0); l < Nl(s); ++l) {
1031  for (size_t i(0); i < (Pl[s]).dim(); ++i) {
1032 
1033  valarray<float> InPx( Npx(s) );
1034  for (size_t ipx(0); ipx < Npx(s)-1; ++ipx) { // at each location in px
1035 
1036 // // At each location in |p|
1037  size_t ip(0);
1038  while (polarR[s](ipx,ip) < 0 ) { ++ip; }
1039 
1040  p_ip1 = polarR[s](ipx,ip);
1041  integrant_high =(float( (df(i))(ip,x0).real() ))*Pl[s](i)(ipx, ip) ;
1042 
1043 
1044 
1045  InPx[ipx] += p_ip1 * (0.5*p_ip1) * integrant_high;
1046  ++ip;
1047  while ( (ip < Np(s) ) && (polarR[s](ipx,ip) > 0) ) {
1048  p_im1 = p_ip1;
1049  integrant_low = integrant_high;
1050 
1051  p_ip1 = polarR[s](ipx,ip);
1052  integrant_high =(float( (df(i))(ip,x0).real() )) * Pl[s](i)(ipx, ip);
1053 
1054  InPx[ipx] += p_im1 * (0.5*(p_ip1-p_im1)) * integrant_low;
1055  InPx[ipx] += p_ip1 * (0.5*(p_ip1-p_im1)) * integrant_high;
1056 
1057  ++ip;
1058  }
1059  // if (i==3) std::cout << "\n p1x1[s](" << ipx << ',' << ip << ") = " << InPx[ipx] << "\n";
1060  (p1x1[s])[ipx] = integrant_high;
1061  }
1062  }*/
1063 
1064  int sgn(1);
1065  size_t midpoint,ipxp,ipxm;
1066 
1067  for (size_t i(0); i < df.dim(); ++i)
1068  {
1069  dummyvec = df(i).xVec(x0);
1070 
1071  midpoint = df(0,0).nump();
1072  p1x1[s][midpoint] = (float) ((dummyvec[0]).real());
1073  for (size_t ip(0); ip < df(0,0).nump(); ++ip)
1074  {
1075  ipxm = midpoint - 1 - ip;
1076  ipxp = midpoint + 1 + ip;
1077 
1078  (p1x1[s])[ipxm] += (float) (sgn*(dummyvec[ip]).real());
1079  (p1x1[s])[ipxp] += (float) ((dummyvec[ip]).real());
1080  }
1081  sgn *= -1;
1082  }
1083 
1084 // p1x1[s] *= 2.0 * M_PI;
1085  // exit(1);
1086  return p1x1[s];
1087 }
1088 //--------------------------------------------------------------
1089 //**************************************************************
1090 
1091 //**************************************************************
1092 //--------------------------------------------------------------
1093 Output_Data::fx1_1D::fx1_1D(const Grid_Info& _G) {
1094 
1095 // Generate the required structures
1096  for (size_t s(0); s < _G.axis.pdim(); ++s) {
1097  pmax.push_back( static_cast<float>(_G.axis.pmax(s)) );
1098  pmin.push_back( static_cast<float>(_G.axis.pmin(s)) );
1099  nump.push_back( static_cast<float>(_G.axis.Np(s)));
1100  // f0x1.push_back( valarray<float>( _G.axis.Npx(s)) );
1101  }
1102 // f0x1 = new valarray<float>(Npx);
1103 }
1104 //- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
1105 
1106 Output_Data::fx1_1D::fx1_1D( const fx1_1D& other) {
1107 
1108  for (size_t s(0); s < other.Species(); ++s) {
1109  pmin.push_back( other.Pmin(s) );
1110  pmax.push_back( other.Pmax(s) );
1111  nump.push_back( other.Np(s));
1112  // fx1.push_back( other.f_x1(s) );
1113  }
1114 }
1115 //- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
1116 
1117 Output_Data::fx1_1D::~fx1_1D(){
1118 
1119 }
1120 //-------------------------------------------------------------
1121 //--------------------------------------------------------------
1122 // Turn the Distribution function at some spatial location x0
1123 // into a cartesian grid.
1124 Array2D<float> Output_Data::fx1_1D::operator()(DistFunc1D& df, size_t l, size_t m, size_t x0, size_t s) {
1125 
1126  // valarray<float> fx1(0.0,Npx(s)); // this is a valarray<float>
1127  Array2D<float> fx1(Np(s),2); // this is a valarray<float>//
1128 
1129  for (size_t ip(0); ip < Np(s); ++ip) {
1130  // At each location in |p|
1131  // fx1[0.5*(Npx(s)-1)+ip] = (float( (df(whichdist))(ip,x0).real() ));
1132  // fx1[0.5*(Npx(s)-1)-ip] = (float( (df(whichdist))(ip,x0).real() ));
1133 
1134  // fx1(0.5*(Npx(s)-1)+1+ip,1) = (float( (df(whichdist))(ip,x0).real() ));
1135  // fx1(0.5*(Npx(s)-1)-1-ip,1) = (float( (df(whichdist))(ip,x0).real() ));
1136 
1137  // fx1(0.5*(Npx(s)-1)+1+ip,2) = (float( (df(whichdist))(ip,x0).imag() ));
1138  // fx1(0.5*(Npx(s)-1)-1-ip,2) = (float( (df(whichdist))(ip,x0).imag() ));
1139  //
1140  fx1(ip,0) = (float( (df(l,m))(ip,x0).real() ));
1141  fx1(ip,1) = (float( (df(l,m))(ip,x0).imag() ));
1142 
1143  // std::cout << "fx1(" << ip << ") = " << (df(whichdist))(ip,x0).real() << "\n";
1144  // << fx1(0.5*(Npx(s)-1)+1+ip,1) << "," << fx1(0.5*(Npx(s)-1)+1+ip,2) << "\n";
1145  }
1146  // fx1[Npx(s)] = (float( (df(whichdist))(0,x0).real() ));
1147  // fx1(,1) = (float( (df(whichdist))(0,x0).real() ));
1148  // fx1(0.5*(Npx(s)-1),2) = (float( (df(whichdist))(0,x0).imag() ));
1149  return fx1;
1150 }
1151 //--------------------------------------------------------------
1152 //**************************************************************
1153 
1154 //**************************************************************
1155 //--------------------------------------------------------------
1156 Output_Data::p2p1x1_1D::p2p1x1_1D( const Grid_Info& _G) {
1157 
1158 // Generate the required structures
1159  for (size_t s(0); s < _G.axis.pdim(); ++s) {
1160  pmax.push_back( static_cast<float>(_G.axis.pmax(s)) );
1161  pmin.push_back( static_cast<float>(_G.axis.pmin(s)) );
1162  numl.push_back( _G.l0[s]);
1163  numm.push_back( _G.m0[s]);
1164 
1165  Pl.push_back( PLegendre2D( _G.l0[s], _G.m0[s], _G.axis.Np(s), static_cast<float>(_G.axis.pmin(s)),
1166  static_cast<float>(_G.axis.pmax(s)), _G.axis.Npx(s), _G.axis.Npy(s) ) );
1167  pr.push_back( Array2D<float>( _G.axis.Npx(s),_G.axis.Npy(s)));
1168 
1169 
1170  nextpcell.push_back( Array2D<size_t>( _G.axis.Npx(s),_G.axis.Npy(s)));
1171  distancetothatcell.push_back( Array2D<float>( _G.axis.Npx(s),_G.axis.Npy(s)));
1172 
1173 
1174  p2p1x1.push_back( Array2D<float>( _G.axis.Npx(s),_G.axis.Npy(s)) );
1175  }
1176 
1177  size_t k(0);
1178  // Calculate a 2D array for the pr for each px and py
1179  for (size_t s(0); s < _G.axis.pdim(); ++s) {
1180 
1181  valarray<float> p( Algorithms::MakeAxis( pmin[s], pmax[s], _G.axis.Np(s) )),
1182  px(Algorithms::MakeAxis( float(-1.0)* pmax[s], pmax[s], _G.axis.Npx(s) )),
1183  py(Algorithms::MakeAxis( float(-1.0)* pmax[s], pmax[s], _G.axis.Npy(s) ));
1184 // --->
1185  for (size_t i(0); i < _G.axis.Npx(s); ++i){
1186  for (size_t j(0); j < _G.axis.Npy(s); ++j){
1187  pr[s](i,j) = sqrt(px[i]*px[i]+py[j]*py[j]);
1188 
1189  while ((pr[s](i,j) > p[k]) && (k < _G.axis.Np(s))) ++k;
1190  if (k == 0) ++k;
1191  nextpcell[s](i,j) = k-1;
1192  distancetothatcell[s](i,j) = (p[k] - pr[s](i,j))/(p[1]-p[0]);
1193  // std::cout << "\n\n (px,py,pr,p[k],nextpcell(i,j)) = " << px[i]
1194  // << ", " << py[j]
1195  // << ", " << pr[s](i,j)
1196  // << ", " << p[k]
1197  // << ", " << nextpcell[s](i,j);
1198 
1199  k = 0;
1200 
1201  }
1202  }
1203 // <---
1204  }
1205 
1206 
1207 // f0x1 = new valarray<float>(Npx);
1208 }
1209 //- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
1210 
1211 Output_Data::p2p1x1_1D::p2p1x1_1D( const p2p1x1_1D& other) {
1212 
1213  for (size_t s(0); s < other.Species(); ++s) {
1214  pmin.push_back( other.Pmin(s) );
1215  pmax.push_back( other.Pmax(s) );
1216  numl.push_back( other.Nl(s));
1217  numm.push_back( other.Nm(s));
1218  Pl.push_back(other.PL(s));
1219  pr.push_back(other.prad(s));
1220  nextpcell.push_back(other.npcell(s));
1221  distancetothatcell.push_back(other.dcell(s));
1222  p2p1x1.push_back( other.p2p1_x1(s) );
1223  }
1224 }
1225 //- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
1226 
1227 Output_Data::p2p1x1_1D::~p2p1x1_1D(){
1228 
1229 }
1230 //-------------------------------------------------------------
1231 //--------------------------------------------------------------
1232 // Turn the Distribution function at some spatial location x0
1233 // into a cartesian grid.
1234 Array2D<float> Output_Data::p2p1x1_1D::operator()(DistFunc1D& df, size_t x0, size_t s) {
1235 
1236  p2p1x1[s] = 0.0; // this is a Array2D<float>
1237  // std::cout << "\n\n x0 = " << x0 << " \n";
1238  size_t i(0), temploc(0);
1239  float tempdist(0.0);
1240 
1241  // valarray<float> p( Algorithms::MakeAxis(pmin, pmax, Np));
1242 
1243  for (size_t l(0); l < Nl(s)+1; ++l){
1244  for (size_t m(0); m<((Nm(s)<l)?Nm(s):l)+1; ++m){
1245 
1246  i = ((l < Nm(s)+1)?((l*(l+1))/2+m):(l*(Nm(s)+1)-(Nm(s)*(Nm(s)+1))/2 + m));
1247 
1248  for (size_t ipy(0); ipy < Npy(s); ++ipy) {
1249  for (size_t ipx(0); ipx < Npx(s); ++ipx) {
1250  // At each location in px,py
1251  if (nextpcell[s](ipx,ipy) == 0)
1252  {
1253  temploc = 1;
1254  tempdist = 1.0;
1255  }
1256  else
1257  {
1258  temploc = nextpcell[s](ipx,ipy);
1259  tempdist = distancetothatcell[s](ipx,ipy);
1260  }
1261 
1262  p2p1x1[s](ipx,ipy) += Pl[s](i)(ipx, ipy) *
1263  ((float((df(i))(temploc ,x0).real())*(1.0-tempdist))+
1264  (float((df(i))(temploc-1,x0).real())*( tempdist)));
1265  }
1266  }
1267  }
1268  }
1269  return p2p1x1[s];
1270 }
1271 //--------------------------------------------------------------
1272 //**************************************************************
1273 
1274 //**************************************************************
1275 //--------------------------------------------------------------
1276 
1277 void Output_Data::Output_Preprocessor_1D::operator()(const State1D& Y, const Grid_Info& grid, const size_t tout,
1278  const Parallel_Environment_1D& PE) {
1279 
1280  if (Input::List().o_Ex) {
1281  Ex( Y, grid, tout, PE );
1282  }
1283 
1284  if (Input::List().o_Ey) {
1285  Ey( Y, grid, tout, PE );
1286  }
1287 
1288  if (Input::List().o_Ez) {
1289  Ez( Y, grid, tout, PE );
1290  }
1291 
1292  if (Input::List().o_Bx) {
1293  Bx( Y, grid, tout, PE );
1294  }
1295 
1296  if (Input::List().o_By) {
1297  By( Y, grid, tout, PE );
1298  }
1299 
1300  if (Input::List().o_Bz) {
1301  Bz( Y, grid, tout, PE );
1302  }
1303 
1304  if (Input::List().o_x1x2) {
1305  n( Y, grid, tout, PE );
1306  }
1307 
1308  if (Input::List().o_Temperature) {
1309  T( Y, grid, tout, PE );
1310  }
1311 
1312  if (Input::List().o_Jx) {
1313  Jx( Y, grid, tout, PE );
1314  }
1315 
1316  if (Input::List().o_Jy) {
1317  Jy( Y, grid, tout, PE );
1318  }
1319 
1320  if (Input::List().o_Jz) {
1321  Jz( Y, grid, tout, PE );
1322  }
1323 
1324  if (Input::List().o_Qx) {
1325  Qx( Y, grid, tout, PE );
1326  }
1327 
1328  if (Input::List().o_Qy) {
1329  Qy( Y, grid, tout, PE );
1330  }
1331 
1332  if (Input::List().o_Qz) {
1333  Qz( Y, grid, tout, PE );
1334  }
1335  if (Input::List().o_vNx) {
1336  vNx( Y, grid, tout, PE );
1337  }
1338  if (Input::List().o_vNy) {
1339  vNy( Y, grid, tout, PE );
1340  }
1341  if (Input::List().o_vNz) {
1342  vNz( Y, grid, tout, PE );
1343  }
1344 
1345 // if (Input::List().o_p1x1){
1346 // px( Y, grid, tout, PE );
1347 // }
1348 // if (Input::List().o_f0x1){
1349 // f0( Y, grid, tout, PE );
1350 // }
1351 // if (Input::List().o_f10x1){
1352 // f10( Y, grid, tout, PE );
1353 // }
1354 // if (Input::List().o_f11x1){
1355 // f11( Y, grid, tout, PE );
1356 // }
1357 // if (Input::List().o_p2p1x1){
1358 // pxpy( Y, grid, tout, PE );
1359 // }
1360 
1361  if (Input::List().o_Ux) {
1362  Ux( Y, grid, tout, PE );
1363  }
1364 
1365  if (Input::List().o_Uy) {
1366  Uy( Y, grid, tout, PE );
1367  }
1368 
1369  if (Input::List().o_Uz) {
1370  Uz( Y, grid, tout, PE );
1371  }
1372 
1373  if (Input::List().o_Z) {
1374  Z( Y, grid, tout, PE );
1375  }
1376 
1377  if (Input::List().o_ni) {
1378  ni( Y, grid, tout, PE );
1379  }
1380 
1381  if (Input::List().o_Ux) {
1382  Ti( Y, grid, tout, PE );
1383  }
1384 
1385 }
1386 //--------------------------------------------------------------
1387 //--------------------------------------------------------------
1388 
1389 void Output_Data::Output_Preprocessor_1D::distdump(const State1D& Y, const Grid_Info& grid, const size_t tout,
1390  const Parallel_Environment_1D& PE) {
1391 
1392  if (Input::List().o_p1x1){
1393  px( Y, grid, tout, PE );
1394  }
1395  if (Input::List().o_f0x1){
1396  f0( Y, grid, tout, PE );
1397  }
1398  if (Input::List().o_f10x1){
1399  f10( Y, grid, tout, PE );
1400  }
1401  if (Input::List().o_f11x1){
1402  f11( Y, grid, tout, PE );
1403  }
1404  if (Input::List().o_f20x1){
1405  f20( Y, grid, tout, PE );
1406  }
1407  if (Input::List().o_fl0x1){
1408  fl0( Y, grid, tout, PE );
1409  }
1410 // if (Input::List().o_p2p1x1){
1411 // pxpy( Y, grid, tout, PE );
1412 // }
1413 
1414 }
1415 //--------------------------------------------------------------
1416 //--------------------------------------------------------------
1417 // void Parallel_Output::parallel_Ex(size_t step, State1D& Y) {
1418 //--------------------------------------------------------------
1419 // Parallel output for Ex
1420 //--------------------------------------------------------------
1421 void Output_Data::Output_Preprocessor_1D::Ex(const State1D& Y, const Grid_Info& grid, const size_t tout,
1422  const Parallel_Environment_1D& PE) {
1423 
1424  size_t Nbc = Input::List().BoundaryCells;
1425  MPI_Status status;
1426  size_t st(0), bi(0);
1427  size_t outNxLocal(grid.axis.Nx(0) - 2*Nbc);
1428  size_t outNxGlobal(grid.axis.Nxg(0));
1429 
1430  int msg_sz(outNxLocal); //*szy);
1431  float* Exbuf = new float[msg_sz];
1432  valarray<float> ExGlobal(outNxGlobal); //, yglob_axis.dim());
1433 
1434  for(size_t i(0); i < msg_sz; ++i) {
1435  Exbuf[i] = static_cast<float>( Y.EMF().Ex()(Nbc+i).real() );
1436  }
1437 
1438  if (PE.NODES() > 1) {
1439  if (PE.RANK()!=0) {
1440  MPI_Send(Exbuf, msg_sz, MPI_FLOAT, 0, PE.RANK(), MPI_COMM_WORLD);
1441  }
1442  else {
1443  // Fill data for rank = 0
1444  for(size_t i(0); i < outNxLocal; i++) {
1445  ExGlobal[i] = Exbuf[i];
1446  }
1447  // Fill data for rank > 0
1448  for (int rr = 1; rr < PE.NODES(); ++rr){
1449  MPI_Recv(Exbuf, msg_sz, MPI_FLOAT, rr, rr, MPI_COMM_WORLD, &status);
1450  for(size_t i(0); i < outNxLocal; i++) {
1451  ExGlobal[i + outNxLocal*rr] = Exbuf[i];
1452  }
1453  }
1454  }
1455  }
1456  // Fill data for Nodes = 0
1457  else {
1458  for(size_t i(0); i < outNxGlobal; i++) {
1459  ExGlobal[i] = Exbuf[i];
1460  }
1461  }
1462 
1463  if (PE.RANK() == 0) expo.Export_h5("Ex", ExGlobal, tout);
1464 
1465  delete[] Exbuf;
1466 
1467 }
1468 //--------------------------------------------------------------
1469 //--------------------------------------------------------------
1471 void Output_Data::Output_Preprocessor_1D::Ey(const State1D& Y, const Grid_Info& grid, const size_t tout,
1472  const Parallel_Environment_1D& PE) {
1473 
1474  size_t Nbc = Input::List().BoundaryCells;
1475  MPI_Status status;
1476  size_t st(0), bi(0);
1477  size_t outNxLocal(grid.axis.Nx(0) - 2*Nbc);
1478  size_t outNxGlobal(grid.axis.Nxg(0));
1479 
1480  int msg_sz(outNxLocal); //*szy);
1481  float* Eybuf = new float[msg_sz];
1482  valarray<float> EyGlobal(outNxGlobal); //, yglob_axis.dim());
1483 
1484  for(size_t i(0); i < msg_sz; ++i) {
1485  Eybuf[i] = static_cast<float>( Y.EMF().Ey()(Nbc+i).real() );
1486  }
1487 
1488  if (PE.NODES() > 1) {
1489  if (PE.RANK()!=0) {
1490  MPI_Send(Eybuf, msg_sz, MPI_FLOAT, 0, PE.RANK(), MPI_COMM_WORLD);
1491  }
1492  else {
1493  // Fill data for rank = 0
1494  for(size_t i(0); i < outNxLocal; i++) {
1495  EyGlobal[i] = Eybuf[i];
1496  }
1497  // Fill data for rank > 0
1498  for (int rr = 1; rr < PE.NODES(); ++rr){
1499  MPI_Recv(Eybuf, msg_sz, MPI_FLOAT, rr, rr, MPI_COMM_WORLD, &status);
1500  for(size_t i(0); i < outNxLocal; i++) {
1501  EyGlobal[i + outNxLocal*rr] = Eybuf[i];
1502  }
1503  }
1504  }
1505  }
1506  // Fill data for Nodes = 0
1507  else {
1508  for(size_t i(0); i < outNxGlobal; i++) {
1509  EyGlobal[i] = Eybuf[i];
1510  }
1511  }
1512 
1513  if (PE.RANK() == 0) expo.Export_h5("Ey", EyGlobal, tout);
1514 
1515  delete[] Eybuf;
1516 
1517 }
1518 //--------------------------------------------------------------
1519 //--------------------------------------------------------------
1521 void Output_Data::Output_Preprocessor_1D::Ez(const State1D& Y, const Grid_Info& grid, const size_t tout,
1522  const Parallel_Environment_1D& PE) {
1523 
1524  size_t Nbc = Input::List().BoundaryCells;
1525  MPI_Status status;
1526  size_t st(0), bi(0);
1527  size_t outNxLocal(grid.axis.Nx(0) - 2*Nbc);
1528  size_t outNxGlobal(grid.axis.Nxg(0));
1529 
1530  int msg_sz(outNxLocal); //*szy);
1531  float* Ezbuf = new float[msg_sz];
1532  valarray<float> EzGlobal(outNxGlobal); //, yglob_axis.dim());
1533 
1534  for(size_t i(0); i < msg_sz; ++i) {
1535  Ezbuf[i] = static_cast<float>( Y.EMF().Ez()(Nbc+i).real() );
1536  }
1537 
1538  if (PE.NODES() > 1) {
1539  if (PE.RANK()!=0) {
1540  MPI_Send(Ezbuf, msg_sz, MPI_FLOAT, 0, PE.RANK(), MPI_COMM_WORLD);
1541  }
1542  else {
1543  // Fill data for rank = 0
1544  for(size_t i(0); i < outNxLocal; i++) {
1545  EzGlobal[i] = Ezbuf[i];
1546  }
1547  // Fill data for rank > 0
1548  for (int rr = 1; rr < PE.NODES(); ++rr){
1549  MPI_Recv(Ezbuf, msg_sz, MPI_FLOAT, rr, rr, MPI_COMM_WORLD, &status);
1550  for(size_t i(0); i < outNxLocal; i++) {
1551  EzGlobal[i + outNxLocal*rr] = Ezbuf[i];
1552  }
1553  }
1554  }
1555  }
1556  // Fill data for Nodes = 0
1557  else {
1558  for(size_t i(0); i < outNxGlobal; i++) {
1559  EzGlobal[i] = Ezbuf[i];
1560  }
1561  }
1562 
1563  if (PE.RANK() == 0) expo.Export_h5("Ez", EzGlobal, tout);
1564 
1565  delete[] Ezbuf;
1566 
1567 }
1568 //--------------------------------------------------------------
1569 //--------------------------------------------------------------
1571 void Output_Data::Output_Preprocessor_1D::Bx(const State1D& Y, const Grid_Info& grid, const size_t tout,
1572  const Parallel_Environment_1D& PE) {
1573 
1574  size_t Nbc = Input::List().BoundaryCells;
1575  MPI_Status status;
1576  size_t st(0), bi(0);
1577  size_t outNxLocal(grid.axis.Nx(0) - 2*Nbc);
1578  size_t outNxGlobal(grid.axis.Nxg(0));
1579 
1580  int msg_sz(outNxLocal); //*szy);
1581  float* Bxbuf = new float[msg_sz];
1582  valarray<float> BxGlobal(outNxGlobal); //, yglob_axis.dim());
1583 
1584  for(size_t i(0); i < msg_sz; ++i) {
1585  Bxbuf[i] = static_cast<float>( Y.EMF().Bx()(Nbc+i).real() );
1586  }
1587 
1588  if (PE.NODES() > 1) {
1589  if (PE.RANK()!=0) {
1590  MPI_Send(Bxbuf, msg_sz, MPI_FLOAT, 0, PE.RANK(), MPI_COMM_WORLD);
1591  }
1592  else {
1593  // Fill data for rank = 0
1594  for(size_t i(0); i < outNxLocal; i++) {
1595  BxGlobal[i] = Bxbuf[i];
1596  }
1597  // Fill data for rank > 0
1598  for (int rr = 1; rr < PE.NODES(); ++rr){
1599  MPI_Recv(Bxbuf, msg_sz, MPI_FLOAT, rr, rr, MPI_COMM_WORLD, &status);
1600  for(size_t i(0); i < outNxLocal; i++) {
1601  BxGlobal[i + outNxLocal*rr] = Bxbuf[i];
1602  }
1603  }
1604  }
1605  }
1606  // Fill data for Nodes = 0
1607  else {
1608  for(size_t i(0); i < outNxGlobal; i++) {
1609  BxGlobal[i] = Bxbuf[i];
1610  }
1611  }
1612 
1613  if (PE.RANK() == 0) expo.Export_h5("Bx", BxGlobal, tout);
1614 
1615  delete[] Bxbuf;
1616 
1617 }
1618 //--------------------------------------------------------------
1619 //--------------------------------------------------------------
1621 void Output_Data::Output_Preprocessor_1D::By(const State1D& Y, const Grid_Info& grid, const size_t tout,
1622  const Parallel_Environment_1D& PE) {
1623 
1624  size_t Nbc = Input::List().BoundaryCells;
1625  MPI_Status status;
1626  size_t st(0), bi(0);
1627  size_t outNxLocal(grid.axis.Nx(0) - 2*Nbc);
1628  size_t outNxGlobal(grid.axis.Nxg(0));
1629 
1630  int msg_sz(outNxLocal); //*szy);
1631  float* Bybuf = new float[msg_sz];
1632  valarray<float> ByGlobal(outNxGlobal); //, yglob_axis.dim());
1633 
1634  for(size_t i(0); i < msg_sz; ++i) {
1635  Bybuf[i] = static_cast<float>( Y.EMF().By()(Nbc+i).real() );
1636  }
1637 
1638  if (PE.NODES() > 1) {
1639  if (PE.RANK()!=0) {
1640  MPI_Send(Bybuf, msg_sz, MPI_FLOAT, 0, PE.RANK(), MPI_COMM_WORLD);
1641  }
1642  else {
1643  // Fill data for rank = 0
1644  for(size_t i(0); i < outNxLocal; i++) {
1645  ByGlobal[i] = Bybuf[i];
1646  }
1647  // Fill data for rank > 0
1648  for (int rr = 1; rr < PE.NODES(); ++rr){
1649  MPI_Recv(Bybuf, msg_sz, MPI_FLOAT, rr, rr, MPI_COMM_WORLD, &status);
1650  for(size_t i(0); i < outNxLocal; i++) {
1651  ByGlobal[i + outNxLocal*rr] = Bybuf[i];
1652  }
1653  }
1654  }
1655  }
1656  // Fill data for Nodes = 0
1657  else {
1658  for(size_t i(0); i < outNxGlobal; i++) {
1659  ByGlobal[i] = Bybuf[i];
1660  }
1661  }
1662 
1663  if (PE.RANK() == 0) expo.Export_h5("By", ByGlobal, tout);
1664 
1665  delete[] Bybuf;
1666 
1667 }
1668 //--------------------------------------------------------------
1669 //--------------------------------------------------------------
1671 void Output_Data::Output_Preprocessor_1D::Bz(const State1D& Y, const Grid_Info& grid, const size_t tout,
1672  const Parallel_Environment_1D& PE) {
1673 
1674  size_t Nbc = Input::List().BoundaryCells;
1675  MPI_Status status;
1676  size_t st(0), bi(0);
1677  size_t outNxLocal(grid.axis.Nx(0) - 2*Nbc);
1678  size_t outNxGlobal(grid.axis.Nxg(0));
1679 
1680  int msg_sz(outNxLocal); //*szy);
1681  float* Bzbuf = new float[msg_sz];
1682  valarray<float> BzGlobal(outNxGlobal); //, yglob_axis.dim());
1683 
1684  for(size_t i(0); i < msg_sz; ++i) {
1685  Bzbuf[i] = static_cast<float>( Y.EMF().Bz()(Nbc+i).real() );
1686  }
1687 
1688  if (PE.NODES() > 1) {
1689  if (PE.RANK()!=0) {
1690  MPI_Send(Bzbuf, msg_sz, MPI_FLOAT, 0, PE.RANK(), MPI_COMM_WORLD);
1691  }
1692  else {
1693  // Fill data for rank = 0
1694  for(size_t i(0); i < outNxLocal; i++) {
1695  BzGlobal[i] = Bzbuf[i];
1696  }
1697  // Fill data for rank > 0
1698  for (int rr = 1; rr < PE.NODES(); ++rr){
1699  MPI_Recv(Bzbuf, msg_sz, MPI_FLOAT, rr, rr, MPI_COMM_WORLD, &status);
1700  for(size_t i(0); i < outNxLocal; i++) {
1701  BzGlobal[i + outNxLocal*rr] = Bzbuf[i];
1702  }
1703  }
1704  }
1705  }
1706  // Fill data for Nodes = 0
1707  else {
1708  for(size_t i(0); i < outNxGlobal; i++) {
1709  BzGlobal[i] = Bzbuf[i];
1710  }
1711  }
1712 
1713  if (PE.RANK() == 0) expo.Export_h5("Bz", BzGlobal, tout);
1714 
1715  delete[] Bzbuf;
1716 
1717 }
1718 //--------------------------------------------------------------
1719 //--------------------------------------------------------------
1720 void Output_Data::Output_Preprocessor_1D::px(const State1D& Y, const Grid_Info& grid, const size_t tout,
1721  const Parallel_Environment_1D& PE) {
1722  // std::cout << "0 \n";
1723  size_t Nbc = Input::List().BoundaryCells;
1724  MPI_Status status;
1725  size_t st(0), bi(0);
1726  size_t outNxLocal(grid.axis.Nx(0) - 2*Nbc);
1727  size_t outNxGlobal(grid.axis.Nxg(0));
1728 
1729 
1730 
1731  for(int s(0); s < Y.Species(); ++s) {
1732  size_t Npx(2*f_x.Np(s)+1);
1733  int msg_sz(outNxLocal*Npx);
1734  Array2D<float> p1x1Global(Npx,outNxGlobal); //, yglob_axis.dim());
1735  float* pxbuf = new float[Npx*outNxLocal];
1736 
1737  for (size_t i(0); i < outNxLocal; ++i) {
1738 
1739  valarray<float> data1D = px_x( Y.DF(s), i+Nbc, s);
1740 
1741  for (size_t j(0); j < Npx; ++j) {
1742  pxbuf[j+i*Npx]=data1D[j];
1743  }
1744  }
1745 
1746  if (PE.NODES() > 1) {
1747  if (PE.RANK()!=0) {
1748  MPI_Send(pxbuf, msg_sz, MPI_FLOAT, 0, PE.RANK(), MPI_COMM_WORLD);
1749  }
1750  else {
1751  // Fill data for rank = 0
1752  for(size_t i(0); i < outNxLocal; i++) {
1753  for (size_t j(0); j < Npx; ++j) {
1754  p1x1Global(j,i) = pxbuf[j+i*Npx];
1755  }
1756  }
1757  // Fill data for rank > 0
1758  for (int rr = 1; rr < PE.NODES(); ++rr){
1759  MPI_Recv(pxbuf, msg_sz, MPI_FLOAT, rr, rr, MPI_COMM_WORLD, &status);
1760  for(size_t i(0); i < outNxLocal; i++) {
1761  for (size_t j(0); j < Npx; ++j) {
1762  p1x1Global(j,i + outNxLocal*rr) = pxbuf[j+i*Npx];
1763  }
1764  }
1765  }
1766  }
1767  }
1768  else {
1769  for(size_t i(0); i < outNxGlobal; i++) {
1770  for (size_t j(0); j < Npx; ++j) {
1771  p1x1Global(j,i) = pxbuf[j+i*Npx];
1772  }
1773  }
1774  }
1775 
1776  if (PE.RANK() == 0) expo.Export_h5("px-x", p1x1Global, tout, s);
1777 
1778  delete[] pxbuf;
1779  }
1780 
1781 }
1782 //--------------------------------------------------------------
1784 //void Output_Data::Output_Preprocessor_1D::pxpy(const State1D& Y, const Grid_Info& grid, const size_t tout,
1785 // const Parallel_Environment_1D& PE) {
1786 //
1787 // size_t Nbc = Input::List().BoundaryCells;
1788 // MPI_Status status;
1789 // size_t st(0), bi(0);
1790 // size_t outNxLocal(grid.axis.Nx(0) - 2*Nbc);
1791 // size_t outNxGlobal(grid.axis.Nxg(0));
1792 // size_t ind(0);
1793 // for(int s(0); s < Y.Species(); ++s) {
1794 //
1795 // int msg_sz(outNxLocal*pxpy_x.Npx(s)*pxpy_x.Npy(s));
1796 // Array3D<float> pxpyGlobal(pxpy_x.Npx(s),pxpy_x.Npy(s),outNxGlobal); //, yglob_axis.dim());
1797 // float* pbuf = new float[pxpy_x.Npx(s)*pxpy_x.Npy(s)*outNxGlobal];
1798 //
1799 // for (size_t i(0); i < outNxLocal; ++i) {
1800 //
1801 // Array2D<float> data2D = pxpy_x( Y.DF(s), i+Nbc, s);
1802 //
1803 // for (size_t j(0); j < pxpy_x.Npx(s); ++j) {
1804 // for (size_t k(0); k < pxpy_x.Npy(s); ++k) {
1805 // pbuf[ind+i*(pxpy_x.Npx(s)*pxpy_x.Npy(s))]=data2D(j,k);
1806 // ++ind;
1807 // }
1808 // }
1809 // ind = 0;
1810 // }
1811 //
1812 // if (PE.NODES() > 1) {
1813 // if (PE.RANK()!=0) {
1814 // MPI_Send(pbuf, msg_sz, MPI_FLOAT, 0, PE.RANK(), MPI_COMM_WORLD);
1815 // }
1816 // else {
1817 // // Fill data for rank = 0
1818 // for(size_t i(0); i < outNxLocal; i++) {
1819 // ind = 0;
1820 // for (size_t j(0); j < pxpy_x.Npx(s); ++j) {
1821 // for (size_t k(0); k < pxpy_x.Npy(s); ++k) {
1822 // pxpyGlobal(k,j,i) = pbuf[ind+i*(pxpy_x.Npx(s)*pxpy_x.Npy(s))];
1823 // ++ind;
1824 // }
1825 // }
1826 // }
1827 // // Fill data for rank > 0
1828 // for (int rr = 1; rr < PE.NODES(); ++rr){
1829 // MPI_Recv(pbuf, msg_sz, MPI_FLOAT, rr, rr, MPI_COMM_WORLD, &status);
1830 // for(size_t i(0); i < outNxLocal; i++) {
1831 // ind = 0;
1832 // for (size_t j(0); j < pxpy_x.Npx(s); ++j) {
1833 // for (size_t k(0); k < pxpy_x.Npy(s); ++k) {
1834 // pxpyGlobal(k,j,i + outNxLocal*rr) = pbuf[ind+i*(pxpy_x.Npx(s)*pxpy_x.Npy(s))];
1835 // // p1x1Global(j,i + outNxLocal*rr) = pxbuf[j+i*f_x.Np(s)];
1836 // }
1837 // }
1838 // }
1839 // }
1840 // }
1841 // }
1842 // else {
1843 // for(size_t i(0); i < outNxLocal; i++) {
1844 // ind = 0;
1845 // for (size_t j(0); j < pxpy_x.Npx(s); ++j) {
1846 // for (size_t k(0); k < pxpy_x.Npy(s); ++k) {
1847 // pxpyGlobal(k,j,i) = pbuf[ind+i*(pxpy_x.Npx(s)*pxpy_x.Npy(s))];
1848 // ++ind;
1849 // }
1850 // }
1851 // }
1852 // }
1853 //
1854 // if (PE.RANK() == 0) expo.Export_h5("pxpy-x", pxpyGlobal, tout, s);
1855 //
1856 // delete[] pbuf;
1857 // }
1858 //
1859 //}
1860 //--------------------------------------------------------------
1861 //--------------------------------------------------------------
1862 void Output_Data::Output_Preprocessor_1D::f0(const State1D& Y, const Grid_Info& grid, const size_t tout,
1863  const Parallel_Environment_1D& PE) {
1864 
1865  size_t Nbc = Input::List().BoundaryCells;
1866  MPI_Status status;
1867 
1868  size_t outNxLocal(grid.axis.Nx(0) - 2*Nbc);
1869  size_t outNxGlobal(grid.axis.Nxg(0));
1870 
1871  for(int s(0); s < Y.Species(); ++s) {
1872 
1873  int msg_sz(2*outNxLocal*f_x.Np(s));
1874  Array3D<float> f0x1Global(f_x.Np(s),outNxGlobal,2); //, yglob_axis.dim());
1875  float* f0xbuf = new float[msg_sz];
1876 
1877  for (size_t i(0); i < outNxLocal; ++i) {
1878 
1879  Array2D<float> data2D = f_x( Y.DF(s),0,0, i+Nbc, s);
1880 
1881  for (size_t j(0); j < f_x.Np(s); ++j) {
1882  // std::cout << "\n f0(" << i << "," << j << ") = (" << data2D(j,1) << "," << data2D(j,2) <<")";
1883  f0xbuf[2*j+ 2*i*f_x.Np(s)]=data2D(j,0);
1884 
1885  f0xbuf[2*j+1+ 2*i*f_x.Np(s)]=data2D(j,1);
1886  }
1887  }
1888 
1889  if (PE.NODES() > 1) {
1890  if (PE.RANK()!=0) {
1891  MPI_Send(f0xbuf, msg_sz, MPI_FLOAT, 0, PE.RANK(), MPI_COMM_WORLD);
1892  }
1893  else {
1894  // Fill data for rank = 0
1895  for(size_t i(0); i < outNxLocal; i++) {
1896 
1897  for (size_t j(0); j < f_x.Np(s); ++j) {
1898  f0x1Global(j,i,0) = f0xbuf[2*j+ 2*i*f_x.Np(s)];
1899  f0x1Global(j,i,1) = f0xbuf[2*j+1+ 2*i*f_x.Np(s)];
1900  }
1901  }
1902  // Fill data for rank > 0
1903  for (int rr = 1; rr < PE.NODES(); ++rr){
1904  MPI_Recv(f0xbuf, msg_sz, MPI_FLOAT, rr, rr, MPI_COMM_WORLD, &status);
1905 
1906  for(size_t i(0); i < outNxLocal; i++) {
1907  for (size_t j(0); j < f_x.Np(s); ++j) {
1908  f0x1Global(j,i + outNxLocal*rr,0) = f0xbuf[2*j+ 2*i*f_x.Np(s)];
1909  f0x1Global(j,i + outNxLocal*rr,1) = f0xbuf[2*j+1+ 2*i*f_x.Np(s)];
1910  }
1911  }
1912  }
1913  }
1914  }
1915  else {
1916  for(size_t i(0); i < outNxGlobal; i++) {
1917 
1918  for (size_t j(0); j < f_x.Np(s); ++j) {
1919  f0x1Global(j,i,0) = f0xbuf[2*j+ 2*i*f_x.Np(s)];
1920  f0x1Global(j,i,1) = f0xbuf[2*j+1+ 2*i*f_x.Np(s)];
1921  }
1922  }
1923  }
1924 
1925  if (PE.RANK() == 0) expo.Export_h5("f0-x", f0x1Global, tout, s);
1926 
1927  delete[] f0xbuf;
1928  }
1929 
1930 }
1931 //--------------------------------------------------------------
1932 //--------------------------------------------------------------
1933 //--------------------------------------------------------------
1934 // void Output_Data::Output_Preprocessor_1D::f0(const State1D& Y, const Grid_Info& grid, const size_t tout,
1935 // const Parallel_Environment_1D& PE) {
1936 // // std::cout << "0 \n";
1937 // size_t Nbc = Input::List().BoundaryCells;
1938 // MPI_Status status;
1939 
1940 // size_t outNxLocal(grid.axis.Nx(0) - 2*Nbc);
1941 // size_t outNxGlobal(grid.axis.Nxg(0));
1942 
1943 // for(int s(0); s < Y.Species(); ++s) {
1944 
1945 // int msg_sz(outNxLocal*f_x.Npx(s));
1946 // Array2D<float> f0x1Global(f_x.Npx(s),outNxGlobal); //, yglob_axis.dim());
1947 // float* f0xbuf = new float[f_x.Npx(s)*outNxLocal];
1948 
1949 // for (size_t i(0); i < outNxLocal; ++i) {
1950 
1951 // valarray<float> data1D = f_x( Y.DF(s),0, i+Nbc, s);
1952 
1953 // for (size_t j(0); j < f_x.Npx(s); ++j) {
1954 // f0xbuf[j+i*f_x.Npx(s)]=data1D[j];
1955 // }
1956 // }
1957 
1958 // if (PE.NODES() > 1) {
1959 // if (PE.RANK()!=0) {
1960 // MPI_Send(f0xbuf, msg_sz, MPI_FLOAT, 0, PE.RANK(), MPI_COMM_WORLD);
1961 // }
1962 // else {
1963 // // Fill data for rank = 0
1964 // for(size_t i(0); i < outNxLocal; i++) {
1965 // for (size_t j(0); j < f_x.Npx(s); ++j) {
1966 // f0x1Global(j,i) = f0xbuf[j+i*f_x.Npx(s)];
1967 // }
1968 // }
1969 // // Fill data for rank > 0
1970 // for (int rr = 1; rr < PE.NODES(); ++rr){
1971 // MPI_Recv(f0xbuf, msg_sz, MPI_FLOAT, rr, rr, MPI_COMM_WORLD, &status);
1972 // for(size_t i(0); i < outNxLocal; i++) {
1973 // for (size_t j(0); j < f_x.Npx(s); ++j) {
1974 // f0x1Global(j,i + outNxLocal*rr) = f0xbuf[j+i*f_x.Npx(s)];
1975 // }
1976 // }
1977 // }
1978 // }
1979 // }
1980 // else {
1981 // for(size_t i(0); i < outNxGlobal; i++) {
1982 // for (size_t j(0); j < f_x.Npx(s); ++j) {
1983 // f0x1Global(j,i) = f0xbuf[j+i*f_x.Npx(s)];
1984 // }
1985 // }
1986 // }
1987 
1988 // if (PE.RANK() == 0) expo.Export_h5("f0-x", f0x1Global, tout, s);
1989 
1990 // delete[] f0xbuf;
1991 // }
1992 
1993 // }
1994 // //--------------------------------------------------------------
1995 // //--------------------------------------------------------------
1996 //--------------------------------------------------------------
1997 void Output_Data::Output_Preprocessor_1D::f10(const State1D& Y, const Grid_Info& grid, const size_t tout,
1998  const Parallel_Environment_1D& PE) {
1999 
2000  size_t Nbc = Input::List().BoundaryCells;
2001  MPI_Status status;
2002 
2003  size_t outNxLocal(grid.axis.Nx(0) - 2*Nbc);
2004  size_t outNxGlobal(grid.axis.Nxg(0));
2005 
2006  for(int s(0); s < Y.Species(); ++s) {
2007 
2008  int msg_sz(2*outNxLocal*f_x.Np(s));
2009  Array3D<float> f0x1Global(f_x.Np(s),outNxGlobal,2); //, yglob_axis.dim());
2010  float* f0xbuf = new float[msg_sz];
2011 
2012  for (size_t i(0); i < outNxLocal; ++i) {
2013 
2014  Array2D<float> data2D = f_x( Y.DF(s), 1, 0, i+Nbc, s);
2015 
2016  for (size_t j(0); j < f_x.Np(s); ++j) {
2017  // std::cout << "\n f0(" << i << "," << j << ") = (" << data2D(j,1) << "," << data2D(j,2) <<")";
2018  f0xbuf[2*j+ 2*i*f_x.Np(s)]=data2D(j,0);
2019 
2020  f0xbuf[2*j+1+ 2*i*f_x.Np(s)]=data2D(j,1);
2021  }
2022  }
2023 
2024  if (PE.NODES() > 1) {
2025  if (PE.RANK()!=0) {
2026  MPI_Send(f0xbuf, msg_sz, MPI_FLOAT, 0, PE.RANK(), MPI_COMM_WORLD);
2027  }
2028  else {
2029  // Fill data for rank = 0
2030  for(size_t i(0); i < outNxLocal; i++) {
2031 
2032  for (size_t j(0); j < f_x.Np(s); ++j) {
2033  f0x1Global(j,i,0) = f0xbuf[2*j+ 2*i*f_x.Np(s)];
2034  f0x1Global(j,i,1) = f0xbuf[2*j+1+ 2*i*f_x.Np(s)];
2035  }
2036  }
2037  // Fill data for rank > 0
2038  for (int rr = 1; rr < PE.NODES(); ++rr){
2039  MPI_Recv(f0xbuf, msg_sz, MPI_FLOAT, rr, rr, MPI_COMM_WORLD, &status);
2040 
2041  for(size_t i(0); i < outNxLocal; i++) {
2042  for (size_t j(0); j < f_x.Np(s); ++j) {
2043  f0x1Global(j,i + outNxLocal*rr,0) = f0xbuf[2*j+ 2*i*f_x.Np(s)];
2044  f0x1Global(j,i + outNxLocal*rr,1) = f0xbuf[2*j+1+ 2*i*f_x.Np(s)];
2045  }
2046  }
2047  }
2048  }
2049  }
2050  else {
2051  for(size_t i(0); i < outNxGlobal; i++) {
2052 
2053  for (size_t j(0); j < f_x.Np(s); ++j) {
2054  f0x1Global(j,i,0) = f0xbuf[2*j+ 2*i*f_x.Np(s)];
2055  f0x1Global(j,i,1) = f0xbuf[2*j+1+ 2*i*f_x.Np(s)];
2056  }
2057  }
2058  }
2059 
2060  if (PE.RANK() == 0) expo.Export_h5("f10-x", f0x1Global, tout, s);
2061 
2062  delete[] f0xbuf;
2063  }
2064 
2065 }
2066 //--------------------------------------------------------------
2067 //--------------------------------------------------------------
2068 void Output_Data::Output_Preprocessor_1D::f11(const State1D& Y, const Grid_Info& grid, const size_t tout,
2069  const Parallel_Environment_1D& PE) {
2070 
2071  size_t Nbc = Input::List().BoundaryCells;
2072  MPI_Status status;
2073 
2074  size_t outNxLocal(grid.axis.Nx(0) - 2*Nbc);
2075  size_t outNxGlobal(grid.axis.Nxg(0));
2076 
2077  for(int s(0); s < Y.Species(); ++s) {
2078 
2079  int msg_sz(2*outNxLocal*f_x.Np(s));
2080  Array3D<float> f0x1Global(f_x.Np(s),outNxGlobal,2); //, yglob_axis.dim());
2081  float* f0xbuf = new float[msg_sz];
2082 
2083  for (size_t i(0); i < outNxLocal; ++i) {
2084 
2085  Array2D<float> data2D = f_x( Y.DF(s),1,1, i+Nbc, s);
2086 
2087  for (size_t j(0); j < f_x.Np(s); ++j) {
2088  // std::cout << "\n f0(" << i << "," << j << ") = (" << data2D(j,1) << "," << data2D(j,2) <<")";
2089  f0xbuf[2*j+ 2*i*f_x.Np(s)]=data2D(j,0);
2090 
2091  f0xbuf[2*j+1+ 2*i*f_x.Np(s)]=data2D(j,1);
2092  }
2093  }
2094 
2095  if (PE.NODES() > 1) {
2096  if (PE.RANK()!=0) {
2097  MPI_Send(f0xbuf, msg_sz, MPI_FLOAT, 0, PE.RANK(), MPI_COMM_WORLD);
2098  }
2099  else {
2100  // Fill data for rank = 0
2101  for(size_t i(0); i < outNxLocal; i++) {
2102 
2103  for (size_t j(0); j < f_x.Np(s); ++j) {
2104  f0x1Global(j,i,0) = f0xbuf[2*j+ 2*i*f_x.Np(s)];
2105  f0x1Global(j,i,1) = f0xbuf[2*j+1+ 2*i*f_x.Np(s)];
2106  }
2107  }
2108  // Fill data for rank > 0
2109  for (int rr = 1; rr < PE.NODES(); ++rr){
2110  MPI_Recv(f0xbuf, msg_sz, MPI_FLOAT, rr, rr, MPI_COMM_WORLD, &status);
2111 
2112  for(size_t i(0); i < outNxLocal; i++) {
2113  for (size_t j(0); j < f_x.Np(s); ++j) {
2114  f0x1Global(j,i + outNxLocal*rr,0) = f0xbuf[2*j+ 2*i*f_x.Np(s)];
2115  f0x1Global(j,i + outNxLocal*rr,1) = f0xbuf[2*j+1+ 2*i*f_x.Np(s)];
2116  }
2117  }
2118  }
2119  }
2120  }
2121  else {
2122  for(size_t i(0); i < outNxGlobal; i++) {
2123 
2124  for (size_t j(0); j < f_x.Np(s); ++j) {
2125  f0x1Global(j,i,0) = f0xbuf[2*j+ 2*i*f_x.Np(s)];
2126  f0x1Global(j,i,1) = f0xbuf[2*j+1+ 2*i*f_x.Np(s)];
2127  }
2128  }
2129  }
2130 
2131  if (PE.RANK() == 0) expo.Export_h5("f11-x", f0x1Global, tout, s);
2132 
2133  delete[] f0xbuf;
2134  }
2135 
2136 
2137 }
2138 //--------------------------------------------------------------
2139 //--------------------------------------------------------------
2140 void Output_Data::Output_Preprocessor_1D::f20(const State1D& Y, const Grid_Info& grid, const size_t tout,
2141  const Parallel_Environment_1D& PE) {
2142 
2143  size_t Nbc = Input::List().BoundaryCells;
2144  MPI_Status status;
2145 
2146  size_t outNxLocal(grid.axis.Nx(0) - 2*Nbc);
2147  size_t outNxGlobal(grid.axis.Nxg(0));
2148 
2149  for(int s(0); s < Y.Species(); ++s) {
2150 
2151  int msg_sz(2*outNxLocal*f_x.Np(s));
2152  Array3D<float> f0x1Global(f_x.Np(s),outNxGlobal,2); //, yglob_axis.dim());
2153  float* f0xbuf = new float[msg_sz];
2154 
2155  for (size_t i(0); i < outNxLocal; ++i) {
2156 
2157  Array2D<float> data2D = f_x( Y.DF(s), 2, 0, i+Nbc, s);
2158 
2159  for (size_t j(0); j < f_x.Np(s); ++j) {
2160  // std::cout << "\n f0(" << i << "," << j << ") = (" << data2D(j,1) << "," << data2D(j,2) <<")";
2161  f0xbuf[2*j+ 2*i*f_x.Np(s)]=data2D(j,0);
2162 
2163  f0xbuf[2*j+1+ 2*i*f_x.Np(s)]=data2D(j,1);
2164  }
2165  }
2166 
2167  if (PE.NODES() > 1) {
2168  if (PE.RANK()!=0) {
2169  MPI_Send(f0xbuf, msg_sz, MPI_FLOAT, 0, PE.RANK(), MPI_COMM_WORLD);
2170  }
2171  else {
2172  // Fill data for rank = 0
2173  for(size_t i(0); i < outNxLocal; i++) {
2174 
2175  for (size_t j(0); j < f_x.Np(s); ++j) {
2176  f0x1Global(j,i,0) = f0xbuf[2*j+ 2*i*f_x.Np(s)];
2177  f0x1Global(j,i,1) = f0xbuf[2*j+1+ 2*i*f_x.Np(s)];
2178  }
2179  }
2180  // Fill data for rank > 0
2181  for (int rr = 1; rr < PE.NODES(); ++rr){
2182  MPI_Recv(f0xbuf, msg_sz, MPI_FLOAT, rr, rr, MPI_COMM_WORLD, &status);
2183 
2184  for(size_t i(0); i < outNxLocal; i++) {
2185  for (size_t j(0); j < f_x.Np(s); ++j) {
2186  f0x1Global(j,i + outNxLocal*rr,0) = f0xbuf[2*j+ 2*i*f_x.Np(s)];
2187  f0x1Global(j,i + outNxLocal*rr,1) = f0xbuf[2*j+1+ 2*i*f_x.Np(s)];
2188  }
2189  }
2190  }
2191  }
2192  }
2193  else {
2194  for(size_t i(0); i < outNxGlobal; i++) {
2195 
2196  for (size_t j(0); j < f_x.Np(s); ++j) {
2197  f0x1Global(j,i,0) = f0xbuf[2*j+ 2*i*f_x.Np(s)];
2198  f0x1Global(j,i,1) = f0xbuf[2*j+1+ 2*i*f_x.Np(s)];
2199  }
2200  }
2201  }
2202 
2203  if (PE.RANK() == 0) expo.Export_h5("f20-x", f0x1Global, tout, s);
2204 
2205  delete[] f0xbuf;
2206  }
2207 
2208 }
2209 //--------------------------------------------------------------
2210 //--------------------------------------------------------------
2211 //--------------------------------------------------------------
2212 void Output_Data::Output_Preprocessor_1D::fl0(const State1D& Y, const Grid_Info& grid, const size_t tout,
2213  const Parallel_Environment_1D& PE) {
2214 
2215  size_t Nbc = Input::List().BoundaryCells;
2216  MPI_Status status;
2217 
2218  size_t outNxLocal(grid.axis.Nx(0) - 2*Nbc);
2219  size_t outNxGlobal(grid.axis.Nxg(0));
2220 
2221  for(int s(0); s < Y.Species(); ++s) {
2222 
2223  int msg_sz(2*outNxLocal*f_x.Np(s));
2224  Array3D<float> f0x1Global(f_x.Np(s),outNxGlobal,2); //, yglob_axis.dim());
2225  float* f0xbuf = new float[msg_sz];
2226 
2227  for (size_t i(0); i < outNxLocal; ++i) {
2228 
2229  Array2D<float> data2D = f_x( Y.DF(s), Y.DF(s).l0()-1,0, i+Nbc, s);
2230 
2231  for (size_t j(0); j < f_x.Np(s); ++j) {
2232  // std::cout << "\n f0(" << i << "," << j << ") = (" << data2D(j,1) << "," << data2D(j,2) <<")";
2233  f0xbuf[2*j+ 2*i*f_x.Np(s)]=data2D(j,0);
2234 
2235  f0xbuf[2*j+1+ 2*i*f_x.Np(s)]=data2D(j,1);
2236  }
2237  }
2238 
2239  if (PE.NODES() > 1) {
2240  if (PE.RANK()!=0) {
2241  MPI_Send(f0xbuf, msg_sz, MPI_FLOAT, 0, PE.RANK(), MPI_COMM_WORLD);
2242  }
2243  else {
2244  // Fill data for rank = 0
2245  for(size_t i(0); i < outNxLocal; i++) {
2246 
2247  for (size_t j(0); j < f_x.Np(s); ++j) {
2248  f0x1Global(j,i,0) = f0xbuf[2*j+ 2*i*f_x.Np(s)];
2249  f0x1Global(j,i,1) = f0xbuf[2*j+1+ 2*i*f_x.Np(s)];
2250  }
2251  }
2252  // Fill data for rank > 0
2253  for (int rr = 1; rr < PE.NODES(); ++rr){
2254  MPI_Recv(f0xbuf, msg_sz, MPI_FLOAT, rr, rr, MPI_COMM_WORLD, &status);
2255 
2256  for(size_t i(0); i < outNxLocal; i++) {
2257  for (size_t j(0); j < f_x.Np(s); ++j) {
2258  f0x1Global(j,i + outNxLocal*rr,0) = f0xbuf[2*j+ 2*i*f_x.Np(s)];
2259  f0x1Global(j,i + outNxLocal*rr,1) = f0xbuf[2*j+1+ 2*i*f_x.Np(s)];
2260  }
2261  }
2262  }
2263  }
2264  }
2265  else {
2266  for(size_t i(0); i < outNxGlobal; i++) {
2267 
2268  for (size_t j(0); j < f_x.Np(s); ++j) {
2269  f0x1Global(j,i,0) = f0xbuf[2*j+ 2*i*f_x.Np(s)];
2270  f0x1Global(j,i,1) = f0xbuf[2*j+1+ 2*i*f_x.Np(s)];
2271  }
2272  }
2273  }
2274 
2275  if (PE.RANK() == 0) expo.Export_h5("fl0-x", f0x1Global, tout, s);
2276 
2277  delete[] f0xbuf;
2278  }
2279 
2280 }
2281 //--------------------------------------------------------------
2282 //--------------------------------------------------------------
2283 //--------------------------------------------------------------
2284 //--------------------------------------------------------------
2285 void Output_Data::Output_Preprocessor_1D::n(const State1D& Y, const Grid_Info& grid, const size_t tout,
2286  const Parallel_Environment_1D& PE) {
2287 
2288 
2289  size_t Nbc = Input::List().BoundaryCells;
2290  MPI_Status status;
2291  size_t st(0), bi(0);
2292  size_t outNxLocal(grid.axis.Nx(0) - 2*Nbc);
2293  size_t outNxGlobal(grid.axis.Nxg(0));
2294 
2295  int msg_sz(outNxLocal); //*szy);
2296  float* nbuf = new float[msg_sz];
2297  valarray<float> nGlobal(outNxGlobal); //, yglob_axis.dim());
2298 
2299  for(int s(0); s < Y.Species(); ++s) {
2300 
2301 // valarray<float> pra( Algorithms::MakeAxis( f_x.Pmin(s), f_x.Pmax(s), f_x.Np(s) ) );
2302 
2303  valarray<float> pra( Algorithms::MakeCAxis(float(0.0), f_x.Pmax(s), f_x.Np(s) ) );
2304  for(size_t i(0); i < msg_sz; ++i) {
2305 
2306  // std::cout << "f00n[" << i << "]=" << (Y.SH(s,0,0)).xVec(i+Nbc)[0] << "\n";
2307 
2308  nbuf[i] = 4.0*M_PI*Algorithms::moment( vfloat( (Y.SH(s,0,0)).xVec(i+Nbc) ), pra, 2);
2309  }
2310 
2311  if (PE.NODES() > 1) {
2312  if (PE.RANK()!=0) {
2313  MPI_Send(nbuf, msg_sz, MPI_FLOAT, 0, PE.RANK(), MPI_COMM_WORLD);
2314  }
2315  else {
2316  // Fill data for rank = 0
2317  for(size_t i(0); i < outNxLocal; i++) {
2318  nGlobal[i] = nbuf[i];
2319  }
2320  // Fill data for rank > 0
2321  for (int rr = 1; rr < PE.NODES(); ++rr){
2322  MPI_Recv(nbuf, msg_sz, MPI_FLOAT, rr, rr, MPI_COMM_WORLD, &status);
2323  for(size_t i(0); i < outNxLocal; i++) {
2324  nGlobal[i + outNxLocal*rr] = nbuf[i];
2325  }
2326  }
2327  }
2328  }
2329  // Fill data for Nodes = 0
2330  else {
2331  for(size_t i(0); i < outNxGlobal; i++) {
2332  nGlobal[i] = nbuf[i];
2333  }
2334  }
2335 
2336  if (PE.RANK() == 0) expo.Export_h5("n", nGlobal, tout, s);
2337 
2338  }
2339 
2340  delete[] nbuf;
2341 
2342 }
2343 //--------------------------------------------------------------
2344 
2345 void Output_Data::Output_Preprocessor_1D::T(const State1D& Y, const Grid_Info& grid, const size_t tout,
2346  const Parallel_Environment_1D& PE) {
2347 
2348 
2349  size_t Nbc = Input::List().BoundaryCells;
2350  MPI_Status status;
2351  // size_t st(0), bi(0);
2352  size_t outNxLocal(grid.axis.Nx(0) - 2*Nbc);
2353  size_t outNxGlobal(grid.axis.Nxg(0));
2354 
2355  int msg_sz(outNxLocal); //*szy);
2356  // float* nbuf = new float[msg_sz];
2357  float* tbuf = new float[msg_sz];
2358  valarray<float> tGlobal(outNxGlobal); //, yglob_axis.dim());
2359 
2360  double convert_factor = (2.99792458e8)*(2.99792458e8)*(9.1093829e-31)/(1.602176565e-19);
2361 
2362  for(int s(0); s < Y.Species(); ++s) {
2363 
2364 // valarray<float> pra( Algorithms::MakeAxis( f_x.Pmin(s), f_x.Pmax(s), f_x.Np(s) ) );
2365  valarray<float> pra( Algorithms::MakeCAxis(float(0.0), f_x.Pmax(s), f_x.Np(s) ) );
2366  // float deltapra = pra[1] - pra[0];
2367 
2368  // need nbuf for density normalization
2369  for(size_t i(0); i < msg_sz; ++i) {
2370  tbuf[i] = 4.0*M_PI*Algorithms::moment( vfloat( (Y.SH(s,0,0)).xVec(i+Nbc) ), pra, 4);
2371  tbuf[i] /= 3.0*4.0*M_PI*Algorithms::moment( vfloat( (Y.SH(s,0,0)).xVec(i+Nbc) ), pra, 2);
2372 
2373  // tbuf[i] *= convert_factor/Y.DF(s).mass();
2374 
2375  tbuf[i] *= 1.0/Y.DF(s).mass();
2376 
2377  }
2378 
2379 // for (int ix(0); ix < outNxLocal; ++ix){
2380 // double accumulator = 0.0;
2381 // for(int ip(0); ip < pra.size(); ++ip) {
2382 // double p_squared = pra[ip]*pra[ip];
2383 // double val = Y.SH(s,0,0)(ip,ix+Nbc).real();
2384 // // pr/dx() can be moved outside loop if anyone cares.
2385 // val *= (deltapra) *p_squared*p_squared;
2386 // accumulator += val;
2387 // }
2388 // tbuf[ix] = accumulator;
2389 // // do the normalization
2390 // tbuf[ix] *= ( 4.0*M_PI*2.0/3.0 );
2391 // // do the density part of the normalization...
2392 // tbuf[ix] /= nbuf[ix];
2393 
2394 // // double pth2 = ( TemperatureFromPthSq[ix] );
2395 // tbuf[ix] *= convert_factor;//pth2 / (1.0 + pth2) * convert_factor;
2396 // }
2397 
2398  if (PE.NODES() > 1) {
2399  if (PE.RANK()!=0) {
2400  MPI_Send(tbuf, msg_sz, MPI_FLOAT, 0, PE.RANK(), MPI_COMM_WORLD);
2401  }
2402  else {
2403  // Fill data for rank = 0
2404  for(size_t i(0); i < outNxLocal; i++) {
2405  tGlobal[i] = tbuf[i];
2406  }
2407  // Fill data for rank > 0
2408  for (int rr = 1; rr < PE.NODES(); ++rr){
2409  MPI_Recv(tbuf, msg_sz, MPI_FLOAT, rr, rr, MPI_COMM_WORLD, &status);
2410  for(size_t i(0); i < outNxLocal; i++) {
2411  tGlobal[i + outNxLocal*rr] = tbuf[i];
2412  }
2413  }
2414  }
2415  }
2416  // Fill data for Nodes = 0
2417  else {
2418  for(size_t i(0); i < outNxGlobal; i++) {
2419  tGlobal[i] = tbuf[i];
2420  }
2421  }
2422 
2423  if (PE.RANK() == 0) expo.Export_h5("T", tGlobal, tout, s);
2424 
2425  }
2426 
2427  // delete[] nbuf;
2428  delete[] tbuf;
2429 
2430 }
2431 //--------------------------------------------------------------
2432 //--------------------------------------------------------------
2433 //--------------------------------------------------------------
2434 void Output_Data::Output_Preprocessor_1D::Jx(const State1D& Y, const Grid_Info& grid, const size_t tout,
2435  const Parallel_Environment_1D& PE) {
2436 
2437 
2438  size_t Nbc = Input::List().BoundaryCells;
2439  MPI_Status status;
2440  // size_t st(0), bi(0);
2441  size_t outNxLocal(grid.axis.Nx(0) - 2*Nbc);
2442  size_t outNxGlobal(grid.axis.Nxg(0));
2443 
2444  int msg_sz(outNxLocal); //*szy);
2445  float* Jxbuf = new float[msg_sz];
2446  valarray<float> JxGlobal(outNxGlobal); //, yglob_axis.dim());
2447 
2448  for(int s(0); s < Y.Species(); ++s) {
2449 
2450 // valarray<float> pra( Algorithms::MakeAxis( f_x.Pmin(s), f_x.Pmax(s), f_x.Np(s) ) );
2451  valarray<float> pra( Algorithms::MakeCAxis(float(0.0), f_x.Pmax(s), f_x.Np(s) ) );
2452 
2453  for(size_t i(0); i < msg_sz; ++i) {
2454  Jxbuf[i] = Y.DF(s).q()/Y.DF(s).mass()*4.0/3.0*M_PI*Algorithms::moment( vfloat( (Y.SH(s,1,0)).xVec(i+Nbc) ), pra, 3);
2455  }
2456 
2457  if (PE.NODES() > 1) {
2458  if (PE.RANK()!=0) {
2459  MPI_Send(Jxbuf, msg_sz, MPI_FLOAT, 0, PE.RANK(), MPI_COMM_WORLD);
2460  }
2461  else {
2462  // Fill data for rank = 0
2463  for(size_t i(0); i < outNxLocal; i++) {
2464  JxGlobal[i] = Jxbuf[i];
2465  }
2466  // Fill data for rank > 0
2467  for (int rr = 1; rr < PE.NODES(); ++rr){
2468  MPI_Recv(Jxbuf, msg_sz, MPI_FLOAT, rr, rr, MPI_COMM_WORLD, &status);
2469  for(size_t i(0); i < outNxLocal; i++) {
2470  JxGlobal[i + outNxLocal*rr] = Jxbuf[i];
2471  }
2472  }
2473  }
2474  }
2475  // Fill data for Nodes = 0
2476  else {
2477  for(size_t i(0); i < outNxGlobal; i++) {
2478  JxGlobal[i] = Jxbuf[i];
2479  }
2480  }
2481 
2482  if (PE.RANK() == 0) expo.Export_h5("Jx", JxGlobal, tout, s);
2483 
2484  }
2485 
2486  delete[] Jxbuf;
2487 
2488 }
2489 //--------------------------------------------------------------
2490 //--------------------------------------------------------------
2491 //--------------------------------------------------------------
2492 //--------------------------------------------------------------
2493 void Output_Data::Output_Preprocessor_1D::Jy(const State1D& Y, const Grid_Info& grid, const size_t tout,
2494  const Parallel_Environment_1D& PE) {
2495 
2496 
2497  size_t Nbc = Input::List().BoundaryCells;
2498  MPI_Status status;
2499  // size_t st(0), bi(0);
2500  size_t outNxLocal(grid.axis.Nx(0) - 2*Nbc);
2501  size_t outNxGlobal(grid.axis.Nxg(0));
2502 
2503  int msg_sz(outNxLocal); //*szy);
2504  float* Jybuf = new float[msg_sz];
2505  valarray<float> JyGlobal(outNxGlobal); //, yglob_axis.dim());
2506 
2507  for(int s(0); s < Y.Species(); ++s) {
2508 
2509 // valarray<float> pra( Algorithms::MakeAxis( f_x.Pmin(s), f_x.Pmax(s), f_x.Np(s) ) );
2510  valarray<float> pra( Algorithms::MakeCAxis(float(0.0), f_x.Pmax(s), f_x.Np(s) ) );
2511 
2512  for(size_t i(0); i < msg_sz; ++i) {
2513  Jybuf[i] = Y.DF(s).q()/Y.DF(s).mass()*8.0/3.0*M_PI*Algorithms::moment( vfloat( (Y.SH(s,1,1)).xVec(i+Nbc) ), pra, 3);
2514  }
2515 
2516  if (PE.NODES() > 1) {
2517  if (PE.RANK()!=0) {
2518  MPI_Send(Jybuf, msg_sz, MPI_FLOAT, 0, PE.RANK(), MPI_COMM_WORLD);
2519  }
2520  else {
2521  // Fill data for rank = 0
2522  for(size_t i(0); i < outNxLocal; i++) {
2523  JyGlobal[i] = Jybuf[i];
2524  }
2525  // Fill data for rank > 0
2526  for (int rr = 1; rr < PE.NODES(); ++rr){
2527  MPI_Recv(Jybuf, msg_sz, MPI_FLOAT, rr, rr, MPI_COMM_WORLD, &status);
2528  for(size_t i(0); i < outNxLocal; i++) {
2529  JyGlobal[i + outNxLocal*rr] = Jybuf[i];
2530  }
2531  }
2532  }
2533  }
2534  // Fill data for Nodes = 0
2535  else {
2536  for(size_t i(0); i < outNxGlobal; i++) {
2537  JyGlobal[i] = Jybuf[i];
2538  }
2539  }
2540 
2541  if (PE.RANK() == 0) expo.Export_h5("Jy", JyGlobal, tout, s);
2542 
2543  }
2544 
2545  delete[] Jybuf;
2546 
2547 }
2548 //--------------------------------------------------------------
2549 //--------------------------------------------------------------
2550 //--------------------------------------------------------------
2551 void Output_Data::Output_Preprocessor_1D::Jz(const State1D& Y, const Grid_Info& grid, const size_t tout,
2552  const Parallel_Environment_1D& PE) {
2553 
2554 
2555  size_t Nbc = Input::List().BoundaryCells;
2556  MPI_Status status;
2557  // size_t st(0), bi(0);
2558  size_t outNxLocal(grid.axis.Nx(0) - 2*Nbc);
2559  size_t outNxGlobal(grid.axis.Nxg(0));
2560 
2561  int msg_sz(outNxLocal); //*szy);
2562  float* Jzbuf = new float[msg_sz];
2563  valarray<float> JzGlobal(outNxGlobal); //, yglob_axis.dim());
2564 
2565  for(int s(0); s < Y.Species(); ++s) {
2566 
2567 // valarray<float> pra( Algorithms::MakeAxis( f_x.Pmin(s), f_x.Pmax(s), f_x.Np(s) ) );
2568  valarray<float> pra( Algorithms::MakeCAxis(float(0.0), f_x.Pmax(s), f_x.Np(s) ) );
2569  for(size_t i(0); i < msg_sz; ++i) {
2570  Jzbuf[i] = Y.DF(s).q()/Y.DF(s).mass()*-8.0/3.0*M_PI*Algorithms::moment( vfloat_complex( (Y.SH(s,1,1)).xVec(i+Nbc) ), pra, 3);
2571  }
2572 
2573  if (PE.NODES() > 1) {
2574  if (PE.RANK()!=0) {
2575  MPI_Send(Jzbuf, msg_sz, MPI_FLOAT, 0, PE.RANK(), MPI_COMM_WORLD);
2576  }
2577  else {
2578  // Fill data for rank = 0
2579  for(size_t i(0); i < outNxLocal; i++) {
2580  JzGlobal[i] = Jzbuf[i];
2581  }
2582  // Fill data for rank > 0
2583  for (int rr = 1; rr < PE.NODES(); ++rr){
2584  MPI_Recv(Jzbuf, msg_sz, MPI_FLOAT, rr, rr, MPI_COMM_WORLD, &status);
2585  for(size_t i(0); i < outNxLocal; i++) {
2586  JzGlobal[i + outNxLocal*rr] = Jzbuf[i];
2587  }
2588  }
2589  }
2590  }
2591  // Fill data for Nodes = 0
2592  else {
2593  for(size_t i(0); i < outNxGlobal; i++) {
2594  JzGlobal[i] = Jzbuf[i];
2595  }
2596  }
2597 
2598  if (PE.RANK() == 0) expo.Export_h5("Jz", JzGlobal, tout, s);
2599 
2600  }
2601 
2602  delete[] Jzbuf;
2603 
2604 }
2605 //--------------------------------------------------------------
2606 //--------------------------------------------------------------
2607 void Output_Data::Output_Preprocessor_1D::Qx(const State1D& Y, const Grid_Info& grid, const size_t tout,
2608  const Parallel_Environment_1D& PE) {
2609 
2610 
2611  size_t Nbc = Input::List().BoundaryCells;
2612  MPI_Status status;
2613 
2614  size_t outNxLocal(grid.axis.Nx(0) - 2*Nbc);
2615  size_t outNxGlobal(grid.axis.Nxg(0));
2616 
2617  int msg_sz(outNxLocal); //*szy);
2618 
2619  double Vxbuf;
2620  double Vybuf;
2621  double Vzbuf;
2622  double VxVxbuf;
2623  double VyVybuf;
2624  double VzVzbuf;
2625  double VxVybuf;
2626  double VxVzbuf;
2627  double VyVzbuf;
2628  double nbuf;
2629  double Ubuf;
2630 
2631  float* Qxbuf = new float[msg_sz];
2632  valarray<float> QxGlobal(outNxGlobal); //, yglob_axis.dim());
2633 
2634  for(int s(0); s < Y.Species(); ++s) {
2635 
2636 // valarray<float> pra( Algorithms::MakeAxis( f_x.Pmin(s), f_x.Pmax(s), f_x.Np(s) ) );
2637  valarray<float> pra( Algorithms::MakeCAxis(float(0.0), f_x.Pmax(s), f_x.Np(s) ) );
2638  for(size_t i(0); i < msg_sz; ++i) {
2639 // nbuf = 4.0*M_PI*Algorithms::moment( vfloat( (Y.SH(s,0,0)).xVec(i+Nbc) ), pra, 2);
2640 // Ubuf = 4.0*M_PI*Algorithms::moment( vfloat( (Y.SH(s,0,0)).xVec(i+Nbc) ), pra, 4);
2641 //
2642 // Vxbuf = 4.0/3.0*M_PI*Algorithms::moment( vfloat( (Y.SH(s,1,0)).xVec(i+Nbc) ), pra, 3);
2643 // Vybuf = 4.0/3.0*M_PI*Algorithms::moment( vfloat( (Y.SH(s,1,1)).xVec(i+Nbc) ), pra, 3);
2644 // Vzbuf = 4.0/3.0*M_PI*Algorithms::moment( vfloat_complex( (Y.SH(s,1,1)).xVec(i+Nbc) ), pra, 3);
2645 // Vxbuf /= nbuf;
2646 // Vybuf /= nbuf;
2647 // Vzbuf /= nbuf;
2648 //
2649 // VxVxbuf = 4.0/3.0*2.0/5.0*M_PI*Algorithms::moment( vfloat( (Y.SH(s,2,0)).xVec(i+Nbc) ), pra, 4);
2650 // VyVybuf = 4.0*4.0/5.0*M_PI*Algorithms::moment( vfloat( (Y.SH(s,2,2)).xVec(i+Nbc) ), pra, 4);
2651 // VzVzbuf = -1.0*VyVybuf;
2652 //
2653 // VyVybuf -= VxVxbuf/2.0;
2654 // VzVzbuf -= VxVxbuf/2.0;
2655 //
2656 // VxVxbuf += Ubuf/2.0;
2657 // VyVybuf += Ubuf/2.0;
2658 // VzVzbuf += Ubuf/2.0;
2659 //
2660 // VxVybuf = 4.0*2.0/5.0*M_PI*Algorithms::moment( vfloat( (Y.SH(s,2,1)).xVec(i+Nbc) ), pra, 4);
2661 // VxVzbuf = -4.0*2.0/5.0*M_PI*Algorithms::moment( vfloat_complex( (Y.SH(s,2,1)).xVec(i+Nbc) ), pra, 4);
2662 // VyVzbuf = -4.0*4.0/5.0*M_PI*Algorithms::moment( vfloat_complex( (Y.SH(s,2,2)).xVec(i+Nbc) ), pra, 4);
2663 //
2664 // VxVxbuf /= nbuf;
2665 // VyVybuf /= nbuf;
2666 // VzVzbuf /= nbuf;
2667 // VxVybuf /= nbuf;
2668 // VxVzbuf /= nbuf;
2669 // VyVzbuf /= nbuf;
2670 
2671  Qxbuf[i] = 4.0*M_PI/3.0*Algorithms::moment( vfloat( (Y.SH(s,1,0)).xVec(i+Nbc) ), pra, 5);
2672 // Qxbuf[i] -= Vxbuf *(VxVxbuf +VyVybuf +VzVzbuf )*nbuf ;
2673 // Qxbuf[i] -= 2.0 *( VxVxbuf *Vxbuf + VxVybuf *Vybuf + VxVzbuf*Vzbuf)*nbuf ;
2674 // Qxbuf[i] += 2.0 *( Vxbuf *Vxbuf + Vybuf *Vybuf +Vzbuf*Vzbuf ) * Vxbuf * nbuf ;
2675  Qxbuf[i] *= 0.5;
2676  }
2677 
2678  if (PE.NODES() > 1) {
2679  if (PE.RANK()!=0) {
2680  MPI_Send(Qxbuf, msg_sz, MPI_FLOAT, 0, PE.RANK(), MPI_COMM_WORLD);
2681  }
2682  else {
2683  // Fill data for rank = 0
2684  for(size_t i(0); i < outNxLocal; i++) {
2685  QxGlobal[i] = Qxbuf[i];
2686  }
2687  // Fill data for rank > 0
2688  for (int rr = 1; rr < PE.NODES(); ++rr){
2689  MPI_Recv(Qxbuf, msg_sz, MPI_FLOAT, rr, rr, MPI_COMM_WORLD, &status);
2690  for(size_t i(0); i < outNxLocal; i++) {
2691  QxGlobal[i + outNxLocal*rr] = Qxbuf[i];
2692  }
2693  }
2694  }
2695  }
2696  // Fill data for Nodes = 0
2697  else {
2698  for(size_t i(0); i < outNxGlobal; i++) {
2699  QxGlobal[i] = Qxbuf[i];
2700  }
2701  }
2702 
2703  if (PE.RANK() == 0) expo.Export_h5("Qx", QxGlobal, tout, s);
2704 
2705  }
2706 
2707  delete[] Qxbuf;
2708 
2709 }
2710 //--------------------------------------------------------------
2711 //--------------------------------------------------------------
2712 //--------------------------------------------------------------
2713 void Output_Data::Output_Preprocessor_1D::Qy(const State1D& Y, const Grid_Info& grid, const size_t tout,
2714  const Parallel_Environment_1D& PE) {
2715 
2716 
2717  size_t Nbc = Input::List().BoundaryCells;
2718  MPI_Status status;
2719 
2720  size_t outNxLocal(grid.axis.Nx(0) - 2*Nbc);
2721  size_t outNxGlobal(grid.axis.Nxg(0));
2722 
2723  int msg_sz(outNxLocal); //*szy);
2724 
2725  double Vxbuf;
2726  double Vybuf;
2727  double Vzbuf;
2728  double VxVxbuf;
2729  double VyVybuf;
2730  double VzVzbuf;
2731  double VxVybuf;
2732  double VxVzbuf;
2733  double VyVzbuf;
2734  double nbuf;
2735  double Ubuf;
2736 
2737  float* Qxbuf = new float[msg_sz];
2738  valarray<float> QxGlobal(outNxGlobal); //, yglob_axis.dim());
2739 
2740  for(int s(0); s < Y.Species(); ++s) {
2741 
2742 // valarray<float> pra( Algorithms::MakeAxis( f_x.Pmin(s), f_x.Pmax(s), f_x.Np(s) ) );
2743  valarray<float> pra( Algorithms::MakeCAxis(float(0.0), f_x.Pmax(s), f_x.Np(s) ) );
2744  for(size_t i(0); i < msg_sz; ++i) {
2745 // nbuf = 4.0*M_PI*Algorithms::moment( vfloat( (Y.SH(s,0,0)).xVec(i+Nbc) ), pra, 2);
2746 // Ubuf = 4.0*M_PI*Algorithms::moment( vfloat( (Y.SH(s,0,0)).xVec(i+Nbc) ), pra, 4);
2747 //
2748 // Vxbuf = 4.0/3.0*M_PI*Algorithms::moment( vfloat( (Y.SH(s,1,0)).xVec(i+Nbc) ), pra, 3);
2749 // Vybuf = 4.0/3.0*M_PI*Algorithms::moment( vfloat( (Y.SH(s,1,1)).xVec(i+Nbc) ), pra, 3);
2750 // Vzbuf = 4.0/3.0*M_PI*Algorithms::moment( vfloat_complex( (Y.SH(s,1,1)).xVec(i+Nbc) ), pra, 3);
2751 // Vxbuf /= nbuf;
2752 // Vybuf /= nbuf;
2753 // Vzbuf /= nbuf;
2754 //
2755 // VxVxbuf = 4.0/3.0*2.0/5.0*M_PI*Algorithms::moment( vfloat( (Y.SH(s,2,0)).xVec(i+Nbc) ), pra, 4);
2756 // VyVybuf = 4.0*4.0/5.0*M_PI*Algorithms::moment( vfloat( (Y.SH(s,2,2)).xVec(i+Nbc) ), pra, 4);
2757 // VzVzbuf = -1.0*VyVybuf;
2758 //
2759 // VyVybuf -= VxVxbuf/2.0;
2760 // VzVzbuf -= VxVxbuf/2.0;
2761 //
2762 // VxVxbuf += Ubuf/2.0;
2763 // VyVybuf += Ubuf/2.0;
2764 // VzVzbuf += Ubuf/2.0;
2765 //
2766 // VxVybuf = 4.0*2.0/5.0*M_PI*Algorithms::moment( vfloat( (Y.SH(s,2,1)).xVec(i+Nbc) ), pra, 4);
2767 // VxVzbuf = -4.0*2.0/5.0*M_PI*Algorithms::moment( vfloat_complex( (Y.SH(s,2,1)).xVec(i+Nbc) ), pra, 4);
2768 // VyVzbuf = -4.0*4.0/5.0*M_PI*Algorithms::moment( vfloat_complex( (Y.SH(s,2,2)).xVec(i+Nbc) ), pra, 4);
2769 //
2770 // VxVxbuf /= nbuf;
2771 // VyVybuf /= nbuf;
2772 // VzVzbuf /= nbuf;
2773 // VxVybuf /= nbuf;
2774 // VxVzbuf /= nbuf;
2775 // VyVzbuf /= nbuf;
2776 
2777  Qxbuf[i] = 8.0*M_PI/3.0*Algorithms::moment( vfloat( (Y.SH(s,1,1)).xVec(i+Nbc) ), pra, 5);
2778 // Qxbuf[i] -= Vxbuf *(VxVxbuf +VyVybuf +VzVzbuf )*nbuf ;
2779 // Qxbuf[i] -= 2.0 *( VxVxbuf *Vxbuf + VxVybuf *Vybuf + VxVzbuf*Vzbuf)*nbuf ;
2780 // Qxbuf[i] += 2.0 *( Vxbuf *Vxbuf + Vybuf *Vybuf +Vzbuf*Vzbuf ) * Vxbuf * nbuf ;
2781  Qxbuf[i] *= 0.5;
2782  }
2783 
2784  if (PE.NODES() > 1) {
2785  if (PE.RANK()!=0) {
2786  MPI_Send(Qxbuf, msg_sz, MPI_FLOAT, 0, PE.RANK(), MPI_COMM_WORLD);
2787  }
2788  else {
2789  // Fill data for rank = 0
2790  for(size_t i(0); i < outNxLocal; i++) {
2791  QxGlobal[i] = Qxbuf[i];
2792  }
2793  // Fill data for rank > 0
2794  for (int rr = 1; rr < PE.NODES(); ++rr){
2795  MPI_Recv(Qxbuf, msg_sz, MPI_FLOAT, rr, rr, MPI_COMM_WORLD, &status);
2796  for(size_t i(0); i < outNxLocal; i++) {
2797  QxGlobal[i + outNxLocal*rr] = Qxbuf[i];
2798  }
2799  }
2800  }
2801  }
2802  // Fill data for Nodes = 0
2803  else {
2804  for(size_t i(0); i < outNxGlobal; i++) {
2805  QxGlobal[i] = Qxbuf[i];
2806  }
2807  }
2808 
2809  if (PE.RANK() == 0) expo.Export_h5("Qy", QxGlobal, tout, s);
2810 
2811  }
2812 
2813  delete[] Qxbuf;
2814 
2815 }
2816 //--------------------------------------------------------------
2817 //--------------------------------------------------------------
2818 //--------------------------------------------------------------
2819 void Output_Data::Output_Preprocessor_1D::Qz(const State1D& Y, const Grid_Info& grid, const size_t tout,
2820  const Parallel_Environment_1D& PE) {
2821 
2822 
2823  size_t Nbc = Input::List().BoundaryCells;
2824  MPI_Status status;
2825 
2826  size_t outNxLocal(grid.axis.Nx(0) - 2*Nbc);
2827  size_t outNxGlobal(grid.axis.Nxg(0));
2828 
2829  int msg_sz(outNxLocal); //*szy);
2830 
2831  double Vxbuf;
2832  double Vybuf;
2833  double Vzbuf;
2834  double VxVxbuf;
2835  double VyVybuf;
2836  double VzVzbuf;
2837  double VxVybuf;
2838  double VxVzbuf;
2839  double VyVzbuf;
2840  double nbuf;
2841  double Ubuf;
2842 
2843  float* Qxbuf = new float[msg_sz];
2844  valarray<float> QxGlobal(outNxGlobal); //, yglob_axis.dim());
2845 
2846  for(int s(0); s < Y.Species(); ++s) {
2847 
2848 // valarray<float> pra( Algorithms::MakeAxis( f_x.Pmin(s), f_x.Pmax(s), f_x.Np(s) ) );
2849  valarray<float> pra( Algorithms::MakeCAxis( float(0.0),f_x.Pmax(s), f_x.Np(s) ) );
2850  for(size_t i(0); i < msg_sz; ++i) {
2851 // nbuf = 4.0*M_PI*Algorithms::moment( vfloat( (Y.SH(s,0,0)).xVec(i+Nbc) ), pra, 2);
2852 // Ubuf = 4.0*M_PI*Algorithms::moment( vfloat( (Y.SH(s,0,0)).xVec(i+Nbc) ), pra, 4);
2853 //
2854 // Vxbuf = 4.0/3.0*M_PI*Algorithms::moment( vfloat( (Y.SH(s,1,0)).xVec(i+Nbc) ), pra, 3);
2855 // Vybuf = 4.0/3.0*M_PI*Algorithms::moment( vfloat( (Y.SH(s,1,1)).xVec(i+Nbc) ), pra, 3);
2856 // Vzbuf = 4.0/3.0*M_PI*Algorithms::moment( vfloat_complex( (Y.SH(s,1,1)).xVec(i+Nbc) ), pra, 3);
2857 // Vxbuf /= nbuf;
2858 // Vybuf /= nbuf;
2859 // Vzbuf /= nbuf;
2860 //
2861 // VxVxbuf = 4.0/3.0*2.0/5.0*M_PI*Algorithms::moment( vfloat( (Y.SH(s,2,0)).xVec(i+Nbc) ), pra, 4);
2862 // VyVybuf = 4.0*4.0/5.0*M_PI*Algorithms::moment( vfloat( (Y.SH(s,2,2)).xVec(i+Nbc) ), pra, 4);
2863 // VzVzbuf = -1.0*VyVybuf;
2864 //
2865 // VyVybuf -= VxVxbuf/2.0;
2866 // VzVzbuf -= VxVxbuf/2.0;
2867 //
2868 // VxVxbuf += Ubuf/2.0;
2869 // VyVybuf += Ubuf/2.0;
2870 // VzVzbuf += Ubuf/2.0;
2871 //
2872 // VxVybuf = 4.0*2.0/5.0*M_PI*Algorithms::moment( vfloat( (Y.SH(s,2,1)).xVec(i+Nbc) ), pra, 4);
2873 // VxVzbuf = -4.0*2.0/5.0*M_PI*Algorithms::moment( vfloat_complex( (Y.SH(s,2,1)).xVec(i+Nbc) ), pra, 4);
2874 // VyVzbuf = -4.0*4.0/5.0*M_PI*Algorithms::moment( vfloat_complex( (Y.SH(s,2,2)).xVec(i+Nbc) ), pra, 4);
2875 //
2876 // VxVxbuf /= nbuf;
2877 // VyVybuf /= nbuf;
2878 // VzVzbuf /= nbuf;
2879 // VxVybuf /= nbuf;
2880 // VxVzbuf /= nbuf;
2881 // VyVzbuf /= nbuf;
2882 
2883  Qxbuf[i] = -8.0*M_PI/3.0*Algorithms::moment( vfloat_complex( (Y.SH(s,1,1)).xVec(i+Nbc) ), pra, 5);
2884 // Qxbuf[i] -= Vxbuf *(VxVxbuf +VyVybuf +VzVzbuf )*nbuf ;
2885 // Qxbuf[i] -= 2.0 *( VxVxbuf *Vxbuf + VxVybuf *Vybuf + VxVzbuf*Vzbuf)*nbuf ;
2886 // Qxbuf[i] += 2.0 *( Vxbuf *Vxbuf + Vybuf *Vybuf +Vzbuf*Vzbuf ) * Vxbuf * nbuf ;
2887  Qxbuf[i] *= 0.5;
2888  }
2889 
2890  if (PE.NODES() > 1) {
2891  if (PE.RANK()!=0) {
2892  MPI_Send(Qxbuf, msg_sz, MPI_FLOAT, 0, PE.RANK(), MPI_COMM_WORLD);
2893  }
2894  else {
2895  // Fill data for rank = 0
2896  for(size_t i(0); i < outNxLocal; i++) {
2897  QxGlobal[i] = Qxbuf[i];
2898  }
2899  // Fill data for rank > 0
2900  for (int rr = 1; rr < PE.NODES(); ++rr){
2901  MPI_Recv(Qxbuf, msg_sz, MPI_FLOAT, rr, rr, MPI_COMM_WORLD, &status);
2902  for(size_t i(0); i < outNxLocal; i++) {
2903  QxGlobal[i + outNxLocal*rr] = Qxbuf[i];
2904  }
2905  }
2906  }
2907  }
2908  // Fill data for Nodes = 0
2909  else {
2910  for(size_t i(0); i < outNxGlobal; i++) {
2911  QxGlobal[i] = Qxbuf[i];
2912  }
2913  }
2914 
2915  if (PE.RANK() == 0) expo.Export_h5("Qz", QxGlobal, tout, s);
2916 
2917  }
2918 
2919  delete[] Qxbuf;
2920 
2921 }
2922 //--------------------------------------------------------------
2923 //--------------------------------------------------------------
2924 //--------------------------------------------------------------
2925 //--------------------------------------------------------------
2926 void Output_Data::Output_Preprocessor_1D::vNx(const State1D& Y, const Grid_Info& grid, const size_t tout,
2927  const Parallel_Environment_1D& PE) {
2928 
2929 
2930  size_t Nbc = Input::List().BoundaryCells;
2931  MPI_Status status;
2932  // size_t st(0), bi(0);
2933  size_t outNxLocal(grid.axis.Nx(0) - 2*Nbc);
2934  size_t outNxGlobal(grid.axis.Nxg(0));
2935 
2936  int msg_sz(outNxLocal); //*szy);
2937  float* vNxbuf = new float[msg_sz];
2938  valarray<float> vNxGlobal(outNxGlobal); //, yglob_axis.dim());
2939 
2940  for(int s(0); s < Y.Species(); ++s) {
2941 
2942 // valarray<float> pra( Algorithms::MakeAxis( f_x.Pmin(s), f_x.Pmax(s), f_x.Np(s) ) );
2943 
2944  valarray<float> pra( Algorithms::MakeCAxis( float(0.0),f_x.Pmax(s), f_x.Np(s) ) );
2945 
2946  for(size_t i(0); i < msg_sz; ++i) {
2947  vNxbuf[i] = (float) (1.0 / 6.0 * (Algorithms::moment(vfloat((Y.SH(s, 1, 0)).xVec(i + Nbc) ), pra, 6)
2948  / Algorithms::moment( vfloat( (Y.SH(s,0,0)).xVec(i+Nbc) ), pra, 5)));
2949  }
2950 
2951  if (PE.NODES() > 1) {
2952  if (PE.RANK()!=0) {
2953  MPI_Send(vNxbuf, msg_sz, MPI_FLOAT, 0, PE.RANK(), MPI_COMM_WORLD);
2954  }
2955  else {
2956  // Fill data for rank = 0
2957  for(size_t i(0); i < outNxLocal; i++) {
2958  vNxGlobal[i] = vNxbuf[i];
2959  }
2960  // Fill data for rank > 0
2961  for (int rr = 1; rr < PE.NODES(); ++rr){
2962  MPI_Recv(vNxbuf, msg_sz, MPI_FLOAT, rr, rr, MPI_COMM_WORLD, &status);
2963  for(size_t i(0); i < outNxLocal; i++) {
2964  vNxGlobal[i + outNxLocal*rr] = vNxbuf[i];
2965  }
2966  }
2967  }
2968  }
2969  // Fill data for Nodes = 0
2970  else {
2971  for(size_t i(0); i < outNxGlobal; i++) {
2972  vNxGlobal[i] = vNxbuf[i];
2973  }
2974  }
2975 
2976  if (PE.RANK() == 0) expo.Export_h5("vNx", vNxGlobal, tout, 0);
2977 
2978  }
2979 
2980  delete[] vNxbuf;
2981 
2982 }
2983 //--------------------------------------------------------------
2984 //--------------------------------------------------------------
2985 void Output_Data::Output_Preprocessor_1D::vNy(const State1D& Y, const Grid_Info& grid, const size_t tout,
2986  const Parallel_Environment_1D& PE) {
2987 
2988 
2989  size_t Nbc = Input::List().BoundaryCells;
2990  MPI_Status status;
2991  size_t outNxLocal(grid.axis.Nx(0) - 2*Nbc);
2992  size_t outNxGlobal(grid.axis.Nxg(0));
2993 
2994  int msg_sz(outNxLocal); //*szy);
2995  float* vNxbuf = new float[msg_sz];
2996  valarray<float> vNxGlobal(outNxGlobal);
2997 
2998  for(int s(0); s < Y.Species(); ++s) {
2999 // valarray<float> pra( Algorithms::MakeAxis( f_x.Pmin(s), f_x.Pmax(s), f_x.Np(s) ) );
3000  valarray<float> pra( Algorithms::MakeCAxis( float(0.0), f_x.Pmax(s), f_x.Np(s) ) );
3001 
3002  for(size_t i(0); i < msg_sz; ++i) {
3003  vNxbuf[i] = (float) (2.0 / 6.0 * (Algorithms::moment(vfloat((Y.SH(s, 1, 1)).xVec(i + Nbc) ), pra, 6)
3004  / Algorithms::moment( vfloat( (Y.SH(s,0,0)).xVec(i+Nbc) ), pra, 5)));
3005  }
3006  if (PE.NODES() > 1) {
3007  if (PE.RANK()!=0) {
3008  MPI_Send(vNxbuf, msg_sz, MPI_FLOAT, 0, PE.RANK(), MPI_COMM_WORLD);
3009  }
3010  else {
3011  // Fill data for rank = 0
3012  for(size_t i(0); i < outNxLocal; i++) {
3013  vNxGlobal[i] = vNxbuf[i];
3014  }
3015  // Fill data for rank > 0
3016  for (int rr = 1; rr < PE.NODES(); ++rr){
3017  MPI_Recv(vNxbuf, msg_sz, MPI_FLOAT, rr, rr, MPI_COMM_WORLD, &status);
3018  for(size_t i(0); i < outNxLocal; i++) {
3019  vNxGlobal[i + outNxLocal*rr] = vNxbuf[i];
3020  }
3021  }
3022  }
3023  }
3024  // Fill data for Nodes = 0
3025  else {
3026  for(size_t i(0); i < outNxGlobal; i++) {
3027  vNxGlobal[i] = vNxbuf[i];
3028  }
3029  }
3030  if (PE.RANK() == 0) expo.Export_h5("vNy", vNxGlobal, tout, 0);
3031  }
3032  delete[] vNxbuf;
3033 }
3034 // --------------------------------------------------------------
3035 // --------------------------------------------------------------
3036 //--------------------------------------------------------------
3037 void Output_Data::Output_Preprocessor_1D::vNz(const State1D& Y, const Grid_Info& grid, const size_t tout,
3038  const Parallel_Environment_1D& PE) {
3039 
3040 
3041  size_t Nbc = Input::List().BoundaryCells;
3042  MPI_Status status;
3043  size_t outNxLocal(grid.axis.Nx(0) - 2*Nbc);
3044  size_t outNxGlobal(grid.axis.Nxg(0));
3045 
3046  int msg_sz(outNxLocal);
3047  float* vNxbuf = new float[msg_sz];
3048  valarray<float> vNxGlobal(outNxGlobal);
3049 
3050  for(int s(0); s < Y.Species(); ++s) {
3051 // valarray<float> pra( Algorithms::MakeAxis( f_x.Pmin(s), f_x.Pmax(s), f_x.Np(s) ) );
3052 
3053  valarray<float> pra( Algorithms::MakeCAxis( float(0.0), f_x.Pmax(s), f_x.Np(s) ) );
3054 
3055  for(size_t i(0); i < msg_sz; ++i) {
3056  vNxbuf[i] = (float) (-2.0 / 6.0 * (Algorithms::moment(vfloat_complex((Y.SH(s, 1, 1)).xVec(i + Nbc) ), pra, 6)
3057  / Algorithms::moment( vfloat( (Y.SH(s,0,0)).xVec(i+Nbc) ), pra, 5)));
3058  }
3059  if (PE.NODES() > 1) {
3060  if (PE.RANK()!=0) {
3061  MPI_Send(vNxbuf, msg_sz, MPI_FLOAT, 0, PE.RANK(), MPI_COMM_WORLD);
3062  }
3063  else {
3064  // Fill data for rank = 0
3065  for(size_t i(0); i < outNxLocal; i++) {
3066  vNxGlobal[i] = vNxbuf[i];
3067  }
3068  // Fill data for rank > 0
3069  for (int rr = 1; rr < PE.NODES(); ++rr){
3070  MPI_Recv(vNxbuf, msg_sz, MPI_FLOAT, rr, rr, MPI_COMM_WORLD, &status);
3071  for(size_t i(0); i < outNxLocal; i++) {
3072  vNxGlobal[i + outNxLocal*rr] = vNxbuf[i];
3073  }
3074  }
3075  }
3076  }
3077  // Fill data for Nodes = 0
3078  else {
3079  for(size_t i(0); i < outNxGlobal; i++) {
3080  vNxGlobal[i] = vNxbuf[i];
3081  }
3082  }
3083  if (PE.RANK() == 0) expo.Export_h5("vNz", vNxGlobal, tout, 0);
3084  }
3085  delete[] vNxbuf;
3086 }
3087 // --------------------------------------------------------------
3088 // --------------------------------------------------------------
3089 
3090 
3091 //--------------------------------------------------------------
3092 void Output_Data::Output_Preprocessor_1D::Ux(const State1D& Y, const Grid_Info& grid, const size_t tout,
3093  const Parallel_Environment_1D& PE) {
3094 
3095 
3096  size_t Nbc = Input::List().BoundaryCells;
3097  MPI_Status status;
3098  // size_t st(0), bi(0);
3099  size_t outNxLocal(grid.axis.Nx(0) - 2*Nbc);
3100  size_t outNxGlobal(grid.axis.Nxg(0));
3101 
3102  int msg_sz(outNxLocal); //*szy);
3103  float* Uxbuf = new float[msg_sz];
3104  valarray<float> UxGlobal(outNxGlobal); //, yglob_axis.dim());
3105 
3106  for(size_t i(0); i < msg_sz; ++i) {
3107  Uxbuf[i] = static_cast<float>(Y.HYDRO().vx(i+Nbc));
3108  }
3109 
3110  if (PE.NODES() > 1) {
3111  if (PE.RANK()!=0) {
3112  MPI_Send(Uxbuf, msg_sz, MPI_FLOAT, 0, PE.RANK(), MPI_COMM_WORLD);
3113  }
3114  else {
3115  // Fill data for rank = 0
3116  for(size_t i(0); i < outNxLocal; i++) {
3117  UxGlobal[i] = Uxbuf[i];
3118  }
3119  // Fill data for rank > 0
3120  for (int rr = 1; rr < PE.NODES(); ++rr){
3121  MPI_Recv(Uxbuf, msg_sz, MPI_FLOAT, rr, rr, MPI_COMM_WORLD, &status);
3122  for(size_t i(0); i < outNxLocal; i++) {
3123  UxGlobal[i + outNxLocal*rr] = Uxbuf[i];
3124  }
3125  }
3126  }
3127  }
3128  // Fill data for Nodes = 0
3129  else {
3130  for(size_t i(0); i < outNxGlobal; i++) {
3131  UxGlobal[i] = Uxbuf[i];
3132  }
3133  }
3134 
3135  if (PE.RANK() == 0) expo.Export_h5("Ux", UxGlobal, tout, 0);
3136 
3137 
3138 
3139  delete[] Uxbuf;
3140 
3141 }
3142 //--------------------------------------------------------------
3143 //--------------------------------------------------------------
3144 void Output_Data::Output_Preprocessor_1D::Uy(const State1D& Y, const Grid_Info& grid, const size_t tout,
3145  const Parallel_Environment_1D& PE) {
3146 
3147 
3148  size_t Nbc = Input::List().BoundaryCells;
3149  MPI_Status status;
3150  // size_t st(0), bi(0);
3151  size_t outNxLocal(grid.axis.Nx(0) - 2*Nbc);
3152  size_t outNxGlobal(grid.axis.Nxg(0));
3153 
3154  int msg_sz(outNxLocal); //*szy);
3155  float* Uxbuf = new float[msg_sz];
3156  valarray<float> UxGlobal(outNxGlobal); //, yglob_axis.dim());
3157 
3158  for(size_t i(0); i < msg_sz; ++i) {
3159  Uxbuf[i] = static_cast<float>(Y.HYDRO().vy(i+Nbc));
3160  }
3161 
3162  if (PE.NODES() > 1) {
3163  if (PE.RANK()!=0) {
3164  MPI_Send(Uxbuf, msg_sz, MPI_FLOAT, 0, PE.RANK(), MPI_COMM_WORLD);
3165  }
3166  else {
3167  // Fill data for rank = 0
3168  for(size_t i(0); i < outNxLocal; i++) {
3169  UxGlobal[i] = Uxbuf[i];
3170  }
3171  // Fill data for rank > 0
3172  for (int rr = 1; rr < PE.NODES(); ++rr){
3173  MPI_Recv(Uxbuf, msg_sz, MPI_FLOAT, rr, rr, MPI_COMM_WORLD, &status);
3174  for(size_t i(0); i < outNxLocal; i++) {
3175  UxGlobal[i + outNxLocal*rr] = Uxbuf[i];
3176  }
3177  }
3178  }
3179  }
3180  // Fill data for Nodes = 0
3181  else {
3182  for(size_t i(0); i < outNxGlobal; i++) {
3183  UxGlobal[i] = Uxbuf[i];
3184  }
3185  }
3186 
3187  if (PE.RANK() == 0) expo.Export_h5("Uy", UxGlobal, tout, 0);
3188 
3189 
3190 
3191  delete[] Uxbuf;
3192 
3193 }
3194 //--------------------------------------------------------------
3195 //--------------------------------------------------------------
3196 void Output_Data::Output_Preprocessor_1D::Uz(const State1D& Y, const Grid_Info& grid, const size_t tout,
3197  const Parallel_Environment_1D& PE) {
3198 
3199 
3200  size_t Nbc = Input::List().BoundaryCells;
3201  MPI_Status status;
3202  // size_t st(0), bi(0);
3203  size_t outNxLocal(grid.axis.Nx(0) - 2*Nbc);
3204  size_t outNxGlobal(grid.axis.Nxg(0));
3205 
3206  int msg_sz(outNxLocal); //*szy);
3207  float* Uxbuf = new float[msg_sz];
3208  valarray<float> UxGlobal(outNxGlobal); //, yglob_axis.dim());
3209 
3210 
3211  for(size_t i(0); i < msg_sz; ++i) {
3212  Uxbuf[i] = static_cast<float>(Y.HYDRO().vz(i+Nbc));
3213  }
3214 
3215  if (PE.NODES() > 1) {
3216  if (PE.RANK()!=0) {
3217  MPI_Send(Uxbuf, msg_sz, MPI_FLOAT, 0, PE.RANK(), MPI_COMM_WORLD);
3218  }
3219  else {
3220  // Fill data for rank = 0
3221  for(size_t i(0); i < outNxLocal; i++) {
3222  UxGlobal[i] = Uxbuf[i];
3223  }
3224  // Fill data for rank > 0
3225  for (int rr = 1; rr < PE.NODES(); ++rr){
3226  MPI_Recv(Uxbuf, msg_sz, MPI_FLOAT, rr, rr, MPI_COMM_WORLD, &status);
3227  for(size_t i(0); i < outNxLocal; i++) {
3228  UxGlobal[i + outNxLocal*rr] = Uxbuf[i];
3229  }
3230  }
3231  }
3232  }
3233  // Fill data for Nodes = 0
3234  else {
3235  for(size_t i(0); i < outNxGlobal; i++) {
3236  UxGlobal[i] = Uxbuf[i];
3237  }
3238  }
3239 
3240  if (PE.RANK() == 0) expo.Export_h5("Uz", UxGlobal, tout, 0);
3241 
3242 
3243 
3244  delete[] Uxbuf;
3245 
3246 }
3247 //--------------------------------------------------------------
3248 //--------------------------------------------------------------
3249 void Output_Data::Output_Preprocessor_1D::Z(const State1D& Y, const Grid_Info& grid, const size_t tout,
3250  const Parallel_Environment_1D& PE) {
3251 
3252 
3253  size_t Nbc = Input::List().BoundaryCells;
3254  MPI_Status status;
3255  // size_t st(0), bi(0);
3256  size_t outNxLocal(grid.axis.Nx(0) - 2*Nbc);
3257  size_t outNxGlobal(grid.axis.Nxg(0));
3258 
3259  int msg_sz(outNxLocal); //*szy);
3260  float* Uxbuf = new float[msg_sz];
3261  valarray<float> UxGlobal(outNxGlobal); //, yglob_axis.dim());
3262 
3263 
3264  for(size_t i(0); i < msg_sz; ++i) {
3265  Uxbuf[i] = static_cast<float>(Y.HYDRO().Z(i+Nbc));
3266  }
3267 
3268  if (PE.NODES() > 1) {
3269  if (PE.RANK()!=0) {
3270  MPI_Send(Uxbuf, msg_sz, MPI_FLOAT, 0, PE.RANK(), MPI_COMM_WORLD);
3271  }
3272  else {
3273  // Fill data for rank = 0
3274  for(size_t i(0); i < outNxLocal; i++) {
3275  UxGlobal[i] = Uxbuf[i];
3276  }
3277  // Fill data for rank > 0
3278  for (int rr = 1; rr < PE.NODES(); ++rr){
3279  MPI_Recv(Uxbuf, msg_sz, MPI_FLOAT, rr, rr, MPI_COMM_WORLD, &status);
3280  for(size_t i(0); i < outNxLocal; i++) {
3281  UxGlobal[i + outNxLocal*rr] = Uxbuf[i];
3282  }
3283  }
3284  }
3285  }
3286  // Fill data for Nodes = 0
3287  else {
3288  for(size_t i(0); i < outNxGlobal; i++) {
3289  UxGlobal[i] = Uxbuf[i];
3290  }
3291  }
3292 
3293  if (PE.RANK() == 0) expo.Export_h5("Z", UxGlobal, tout, 0);
3294 
3295 
3296 
3297  delete[] Uxbuf;
3298 
3299 }
3300 //--------------------------------------------------------------
3301 //--------------------------------------------------------------
3302 //
3303 //--------------------------------------------------------------
3304 //--------------------------------------------------------------
3305 void Output_Data::Output_Preprocessor_1D::ni(const State1D& Y, const Grid_Info& grid, const size_t tout,
3306  const Parallel_Environment_1D& PE) {
3307 
3308 
3309  size_t Nbc = Input::List().BoundaryCells;
3310  MPI_Status status;
3311  // size_t st(0), bi(0);
3312  size_t outNxLocal(grid.axis.Nx(0) - 2*Nbc);
3313  size_t outNxGlobal(grid.axis.Nxg(0));
3314 
3315  int msg_sz(outNxLocal); //*szy);
3316  float* nibuf = new float[msg_sz];
3317  valarray<float> niGlobal(outNxGlobal); //, yglob_axis.dim());
3318 
3319  for(size_t i(0); i < msg_sz; ++i) {
3320  nibuf[i] = static_cast<float>(Y.HYDRO().density(i+Nbc));
3321  }
3322 
3323  if (PE.NODES() > 1) {
3324  if (PE.RANK()!=0) {
3325  MPI_Send(nibuf, msg_sz, MPI_FLOAT, 0, PE.RANK(), MPI_COMM_WORLD);
3326  }
3327  else {
3328  // Fill data for rank = 0
3329  for(size_t i(0); i < outNxLocal; i++) {
3330  niGlobal[i] = nibuf[i];
3331  }
3332  // Fill data for rank > 0
3333  for (int rr = 1; rr < PE.NODES(); ++rr){
3334  MPI_Recv(nibuf, msg_sz, MPI_FLOAT, rr, rr, MPI_COMM_WORLD, &status);
3335  for(size_t i(0); i < outNxLocal; i++) {
3336  niGlobal[i + outNxLocal*rr] = nibuf[i];
3337  }
3338  }
3339  }
3340  }
3341  // Fill data for Nodes = 0
3342  else {
3343  for(size_t i(0); i < outNxGlobal; i++) {
3344  niGlobal[i] = nibuf[i];
3345  }
3346  }
3347 
3348  if (PE.RANK() == 0) expo.Export_h5("ni", niGlobal, tout, 0);
3349 
3350 
3351 
3352  delete[] nibuf;
3353 
3354 }
3355 //--------------------------------------------------------------
3356 //--------------------------------------------------------------
3357 //
3358 //--------------------------------------------------------------
3359 //--------------------------------------------------------------
3360 void Output_Data::Output_Preprocessor_1D::Ti(const State1D& Y, const Grid_Info& grid, const size_t tout,
3361  const Parallel_Environment_1D& PE) {
3362 
3363 
3364  size_t Nbc = Input::List().BoundaryCells;
3365  MPI_Status status;
3366  // size_t st(0), bi(0);
3367  size_t outNxLocal(grid.axis.Nx(0) - 2*Nbc);
3368  size_t outNxGlobal(grid.axis.Nxg(0));
3369 
3370  int msg_sz(outNxLocal); //*szy);
3371  float* Thydrobuf = new float[msg_sz];
3372  valarray<float> ThydroGlobal(outNxGlobal); //, yglob_axis.dim());
3373 
3374  for(size_t i(0); i < msg_sz; ++i) {
3375  Thydrobuf[i] = static_cast<float>(511000.0/3.0*Y.HYDRO().temperature(i+Nbc));
3376  }
3377 
3378  if (PE.NODES() > 1) {
3379  if (PE.RANK()!=0) {
3380  MPI_Send(Thydrobuf, msg_sz, MPI_FLOAT, 0, PE.RANK(), MPI_COMM_WORLD);
3381  }
3382  else {
3383  // Fill data for rank = 0
3384  for(size_t i(0); i < outNxLocal; i++) {
3385  ThydroGlobal[i] = Thydrobuf[i];
3386  }
3387  // Fill data for rank > 0
3388  for (int rr = 1; rr < PE.NODES(); ++rr){
3389  MPI_Recv(Thydrobuf, msg_sz, MPI_FLOAT, rr, rr, MPI_COMM_WORLD, &status);
3390  for(size_t i(0); i < outNxLocal; i++) {
3391  ThydroGlobal[i + outNxLocal*rr] = Thydrobuf[i];
3392  }
3393  }
3394  }
3395  }
3396  // Fill data for Nodes = 0
3397  else {
3398  for(size_t i(0); i < outNxGlobal; i++) {
3399  ThydroGlobal[i] = Thydrobuf[i];
3400  }
3401  }
3402 
3403  if (PE.RANK() == 0) expo.Export_h5("Ti", ThydroGlobal, tout, 0);
3404 
3405 
3406 
3407  delete[] Thydrobuf;
3408 
3409 }
3410 //--------------------------------------------------------------
3411 //--------------------------------------------------------------
3412 
3413 //**************************************************************
3414 //--------------------------------------------------------------
3415 //--------------------------------------------------------------
3416 void Export_Files::Xport:: Export_h5(const std::string tag,
3417  std::valarray<float> data,
3418  const size_t& step,
3419  const int spec){
3420 //--------------------------------------------------------------
3421 // Export data to H5 file
3422 //--------------------------------------------------------------
3423 
3424  string filename(Hdr[tag].Directory());
3425 
3426 // Check Header file correctness
3427  if (Hdr[tag].dim() != 1) {
3428  cout << "ERROR "<< tag <<" : " << Hdr[tag].dim() << " dimensions != 1D structure\n";
3429  exit(1);
3430  }
3431 
3432 // Open File
3433  filename.append(tag).append(oH5Fextension(step,spec));
3434  H5::H5File file = hmake_file(filename);
3435 
3436  hsize_t dimsf[1] = {data.size()}; // dataset dimensions
3437  H5::DataSpace dataspace( 1, dimsf );
3438  H5::DataSet dataset = file.createDataSet( tag, H5::PredType::NATIVE_FLOAT, dataspace );
3439 
3440  float xmin = { static_cast<float>( Hdr[tag].min(0)) };
3441  float xmax = { static_cast<float>( Hdr[tag].max(0)) };
3442 
3443  hinit_attr(file, tag, step, xmax, xmin);
3444 
3445  H5::Group group = file.createGroup("/AXIS");
3446 
3447  std::size_t pos = Hdr[tag].label(0).find("[");
3448 
3449  string axismainname = "AXIS1";
3450  float axisrange[2];
3451  axisrange[0] = xmin;
3452  axisrange[1] = xmax;
3453  string axislongname = Hdr[tag].label(0).substr(0,pos);
3454  string axisname = Hdr[tag].label(0).substr(0,pos);
3455  string axistype = "linear";
3456  string axisunits = Hdr[tag].units(0);
3457 
3458  haxis(group, axismainname, axisrange, axislongname, axisname, axistype, axisunits);
3459 
3460  group.close();
3461 
3462  dataset.write( &(data[0]) , H5::PredType::NATIVE_FLOAT );
3463  hfile_add_attr_todataset(dataset, "LONG_NAME", tag);
3464  hfile_add_attr_todataset(dataset, "UNITS", formulary().Label(tag));
3465 
3466  dataset.close();
3467 
3468  hclose_file(file);
3469 
3470 }
3471 //--------------------------------------------------------------
3472 
3473 
3474 //--------------------------------------------------------------
3475 void Export_Files::Xport:: Export_h5(const std::string tag,
3476  Array2D<float> data,
3477  const size_t& step,
3478  const int spec){
3479 //--------------------------------------------------------------
3480 // Export data to H5 file
3481 //--------------------------------------------------------------
3482 
3483  string filename(Hdr[tag].Directory());
3484 
3485 // Check Header file correctness
3486  if (Hdr[tag].dim() != 2) {
3487  cout << "ERROR "<< tag <<" : " << Hdr[tag].dim() << " dimensions != 2D structure\n";
3488  exit(1);
3489  }
3490 
3491 // Open File
3492  filename.append(tag).append(oH5Fextension(step,spec));
3493  H5::H5File file = hmake_file(filename);
3494 
3495  hsize_t dimsf[2] = { data.dim2(), data.dim1() }; // dataset dimensions
3496  H5::DataSpace dataspace( 2, dimsf );
3497  H5::DataSet dataset = file.createDataSet( tag, H5::PredType::NATIVE_FLOAT, dataspace );
3498 
3499  float xmin[2] = { static_cast<float>( Hdr[tag].min(0)),
3500  static_cast<float>( Hdr[tag].min(1)) };
3501 
3502  float xmax[2] = { static_cast<float>( Hdr[tag].max(0)),
3503  static_cast<float>( Hdr[tag].max(1)) };
3504 
3505  hinit_attr2(file, tag, step, xmax, xmin);
3506 
3507  H5::Group group = file.createGroup("/AXIS");
3508 
3509  std::size_t pos = Hdr[tag].label(1).find("[");
3510 
3511  string axismainname = "AXIS1";
3512  float axisrange[2];
3513  axisrange[0] = xmin[1];
3514  axisrange[1] = xmax[1];
3515  string axislongname = Hdr[tag].label(1).substr(0,pos);
3516  string axisname = Hdr[tag].label(1).substr(0,pos);
3517  string axistype = "linear";
3518  string axisunits = Hdr[tag].units(1);
3519 
3520  haxis(group, axismainname, axisrange, axislongname, axisname, axistype, axisunits);
3521 
3522  pos = Hdr[tag].label(0).find("[");
3523 
3524  axismainname = "AXIS2";
3525  axisrange[0] = xmin[0];
3526  axisrange[1] = xmax[0];
3527  axislongname = Hdr[tag].label(0).substr(0,pos);
3528  axisname = Hdr[tag].label(0).substr(0,pos);
3529  axisunits = Hdr[tag].units(0);
3530 
3531  haxis(group, axismainname, axisrange, axislongname, axisname, axistype, axisunits);
3532 
3533  group.close();
3534 
3535  dataset.write( &(data.array())[0] , H5::PredType::NATIVE_FLOAT );
3536  hfile_add_attr_todataset(dataset, "LONG_NAME", tag);
3537  hfile_add_attr_todataset(dataset, "UNITS", formulary().Label(tag));
3538 
3539  dataset.close();
3540 
3541  hclose_file(file);
3542 
3543 }
3544 //--------------------------------------------------------------
3545 
3546 
3547 //--------------------------------------------------------------
3548 void Export_Files::Xport:: Export_h5(const std::string tag,
3549  Array3D<float> data,
3550  const size_t& step,
3551  const int spec){
3552 //--------------------------------------------------------------
3553 // Export data to H5 file
3554 //--------------------------------------------------------------
3555 
3556  string filename(Hdr[tag].Directory());
3557 
3558 // Check Header file correctness
3559  if (Hdr[tag].dim() != 3) {
3560  cout << "ERROR "<< tag <<" : " << Hdr[tag].dim() << " dimensions != 3D structure\n";
3561  exit(1);
3562  }
3563 
3564 // Open File
3565  filename.append(tag).append(oH5Fextension(step,spec));
3566  H5::H5File file = hmake_file(filename);
3567 
3568  hsize_t dimsf[3] = { data.dim3(), data.dim2(), data.dim1() }; // dataset dimensions
3569  H5::DataSpace dataspace( 3, dimsf );
3570  H5::DataSet dataset = file.createDataSet( tag, H5::PredType::NATIVE_FLOAT, dataspace );
3571 
3572  float xmin[3] = { static_cast<float>( Hdr[tag].min(0)),
3573  static_cast<float>( Hdr[tag].min(1)),
3574  static_cast<float>( Hdr[tag].min(2))};
3575  float xmax[3] = { static_cast<float>( Hdr[tag].max(0)),
3576  static_cast<float>( Hdr[tag].max(1)),
3577  static_cast<float>( Hdr[tag].max(2))};
3578 
3579  hinit_attr2(file, tag, step, xmax, xmin);
3580 
3581  H5::Group group = file.createGroup("/AXIS");
3582 
3583  std::size_t pos = Hdr[tag].label(1).find("[");
3584 
3585  string axismainname = "AXIS1";
3586  float axisrange[2];
3587  axisrange[0] = xmin[1];
3588  axisrange[1] = xmax[1];
3589  string axislongname = Hdr[tag].label(1).substr(0,pos);
3590  string axisname = Hdr[tag].label(1).substr(0,pos);
3591  string axistype = "linear";
3592  string axisunits = Hdr[tag].units(1);
3593 
3594  haxis(group, axismainname, axisrange, axislongname, axisname, axistype, axisunits);
3595 
3596  pos = Hdr[tag].label(0).find("[");
3597 
3598  axismainname = "AXIS2";
3599  axisrange[0] = xmin[0];
3600  axisrange[1] = xmax[0];
3601  axislongname = Hdr[tag].label(0).substr(0,pos);
3602  axisname = Hdr[tag].label(0).substr(0,pos);
3603  axisunits = Hdr[tag].units(0);
3604 
3605  haxis(group, axismainname, axisrange, axislongname, axisname, axistype, axisunits);
3606 
3607  pos = Hdr[tag].label(2).find("[");
3608 
3609  axismainname = "AXIS3";
3610  axisrange[0] = xmin[2];
3611  axisrange[1] = xmax[2];
3612  axislongname = Hdr[tag].label(2).substr(0,pos);
3613  axisname = Hdr[tag].label(2).substr(0,pos);
3614  axisunits = Hdr[tag].units(2);
3615 
3616  haxis(group, axismainname, axisrange, axislongname, axisname, axistype, axisunits);
3617 
3618  group.close();
3619 
3620  dataset.write( &(data.array())[0] , H5::PredType::NATIVE_FLOAT );
3621  hfile_add_attr_todataset(dataset, "LONG_NAME", tag);
3622  hfile_add_attr_todataset(dataset, "UNITS", formulary().Label(tag));
3623 
3624  dataset.close();
3625 
3626  hclose_file(file);
3627 
3628 }
3629 //--------------------------------------------------------------
3630 //
3631 //--------------------------------------------------------------
3632 H5::H5File Export_Files::Xport:: hmake_file(string ofilename) {
3633 //--------------------------------------------------------------
3634 // Create HDF5 file
3635 //--------------------------------------------------------------
3636 
3637  const H5std_string FILE_NAME( ofilename );
3638  H5::H5File file( FILE_NAME, H5F_ACC_TRUNC );
3639  return file;
3640 
3641 }
3642 //--------------------------------------------------------------
3643 
3644 
3645 //--------------------------------------------------------------
3646 void Export_Files::Xport:: hclose_file(H5::H5File &file) {
3647 //--------------------------------------------------------------
3648 // Close HDF5 file
3649 //--------------------------------------------------------------
3650 
3651  file.close();
3652 
3653 }
3654 //--------------------------------------------------------------
3655 
3656 
3657 //--------------------------------------------------------------
3658 void Export_Files::Xport:: hinit_attr(H5::H5File &hfilehandle, const std::string tag,
3659  size_t step,
3660  float xmax, float xmin) {
3661 //--------------------------------------------------------------
3662 // Add initial attributes:
3663 // time step, iteration number, name of diagnostic,
3664 // physical time, time units, type, and max and min of axes ranges
3665 //--------------------------------------------------------------
3666 
3667  double dt_out(Input::List().t_stop / Input::List().n_outsteps);
3668 
3669  float dt = { static_cast<float>(
3670  dt_out/static_cast<double>(
3671  size_t(static_cast<int>(dt_out
3672  /Input::List().clf_dp))+1
3673  ))};
3674 
3675  hfile_add_attr(hfilehandle, "DT", dt);
3676 
3677  float time = { static_cast<float>(dt_out*step) };
3678  hfile_add_attr(hfilehandle, "TIME", time);
3679 
3680  hfile_add_attr(hfilehandle, "NAME", tag);
3681 
3682  // float time = { static_cast<float>(dt_out*step*dt) };
3683  // hfile_add_attr(hfilehandle, "TIME", time);
3684 
3685  string timeunits = "1/\\omega_p";
3686  hfile_add_attr(hfilehandle, "TIME UNITS", timeunits);
3687 
3688  string typegrid = "grid";
3689  hfile_add_attr(hfilehandle, "TYPE", typegrid);
3690 
3691  hfile_add_attr(hfilehandle, "XMAX", xmax);
3692 
3693  hfile_add_attr(hfilehandle, "XMIN", xmin);
3694 
3695 }
3696 //--------------------------------------------------------------
3697 //--------------------------------------------------------------
3698 void Export_Files::Xport:: hinit_attr2(H5::H5File &hfilehandle, const std::string tag,
3699  size_t step,
3700  float xmax[], float xmin[]) {
3701 //--------------------------------------------------------------
3702 // Add initial attributes:
3703 // time step, iteration number, name of diagnostic,
3704 // physical time, time units, type, and max and min of axes ranges
3705 //--------------------------------------------------------------
3706 
3707  double dt_out(Input::List().t_stop / Input::List().n_outsteps);
3708 
3709  float dt = { static_cast<float>(
3710  dt_out/static_cast<double>(
3711  size_t(static_cast<int>(dt_out
3712  /Input::List().clf_dp))+1
3713  ))};
3714  hfile_add_attr(hfilehandle, "DT", dt);
3715 
3716  float time = { static_cast<float>(dt_out*step) };
3717  hfile_add_attr(hfilehandle, "TIME", time);
3718 
3719  hfile_add_attr(hfilehandle, "NAME", tag);
3720 
3721  string timeunits = "1/\\omega_p";
3722  hfile_add_attr(hfilehandle, "TIME UNITS", timeunits);
3723 
3724  string typegrid = "grid";
3725  hfile_add_attr(hfilehandle, "TYPE", typegrid);
3726 
3727  hfile_add_attr2(hfilehandle, "XMAX", xmax);
3728 
3729  hfile_add_attr2(hfilehandle, "XMIN", xmin);
3730 
3731 }
3732 //--------------------------------------------------------------
3733 
3734 
3735 //--------------------------------------------------------------
3736 void Export_Files::Xport:: haxis(H5::Group &hgrouphandle, string axismainname, float axisrange[2],
3737  string axislongname, string axisname,
3738  string axistype, string axisunits) {
3739 //--------------------------------------------------------------
3740 // Add axis attributes
3741 //--------------------------------------------------------------
3742 
3743  hsize_t dims2[1] = { 2 };
3744  H5::DataSpace dataspaceg1( 1, dims2 );
3745  H5::DataSet datasetg1 = hgrouphandle.createDataSet(axismainname, H5::PredType::NATIVE_FLOAT,
3746  dataspaceg1);
3747  datasetg1.write( axisrange, H5::PredType::NATIVE_FLOAT );
3748  hfile_add_attr_todataset(datasetg1, "LONG_NAME", axislongname);
3749  hfile_add_attr_todataset(datasetg1, "NAME", axisname);
3750  hfile_add_attr_todataset(datasetg1, "TYPE", axistype);
3751  hfile_add_attr_todataset(datasetg1, "UNITS", axisunits);
3752 
3753  dataspaceg1.close();
3754  datasetg1.close();
3755 
3756 }
3757 //--------------------------------------------------------------
3758 
3759 
3760 //--------------------------------------------------------------
3761 void Export_Files::Xport:: hfile_add_attr(H5::H5File &hfilehandle, string attrname,
3762  int attrdata) {
3763 //--------------------------------------------------------------
3764 // Add 1D integer attribute
3765 //--------------------------------------------------------------
3766 
3767  hsize_t dims1[1] = { 1 };
3768  H5::DataSpace attr_dataspace_1 = H5::DataSpace (1, dims1 );
3769  H5::Attribute attribute_int = hfilehandle.createAttribute( attrname, H5::PredType::NATIVE_INT,
3770  attr_dataspace_1);
3771  attribute_int.write( H5::PredType::NATIVE_INT, &attrdata);
3772  attr_dataspace_1.close();
3773  attribute_int.close();
3774 
3775 }
3776 
3777 //--------------------------------------------------------------
3778 
3779 
3780 //--------------------------------------------------------------
3781 void Export_Files::Xport:: hfile_add_attr2(H5::H5File &hfilehandle, string attrname,
3782  int attrdata[2]) {
3783 //--------------------------------------------------------------
3784 // Add 2D integer attribute
3785 //--------------------------------------------------------------
3786 
3787  hsize_t dims1[1] = { 2 };
3788  H5::DataSpace attr_dataspace_2 = H5::DataSpace (1, dims1 );
3789  H5::Attribute attribute_int = hfilehandle.createAttribute( attrname, H5::PredType::NATIVE_INT,
3790  attr_dataspace_2);
3791  attribute_int.write( H5::PredType::NATIVE_INT, attrdata);
3792  attr_dataspace_2.close();
3793  attribute_int.close();
3794 
3795 }
3796 
3797 //--------------------------------------------------------------
3798 
3799 
3800 //--------------------------------------------------------------
3801 void Export_Files::Xport:: hfile_add_attr(H5::H5File &hfilehandle, string attrname,
3802  float attrdata) {
3803 //--------------------------------------------------------------
3804 // Add 1D float attribute
3805 //--------------------------------------------------------------
3806 
3807  hsize_t dims1[1] = { 1 };
3808  H5::DataSpace attr_dataspace_1 = H5::DataSpace (1, dims1 );
3809  H5::Attribute attribute_float = hfilehandle.createAttribute( attrname, H5::PredType::NATIVE_FLOAT,
3810  attr_dataspace_1);
3811  attribute_float.write( H5::PredType::NATIVE_FLOAT, &attrdata);
3812  attr_dataspace_1.close();
3813  attribute_float.close();
3814 
3815 }
3816 
3817 //--------------------------------------------------------------
3818 
3819 
3820 //--------------------------------------------------------------
3821 void Export_Files::Xport:: hfile_add_attr2(H5::H5File &hfilehandle, string attrname,
3822  float attrdata[2]) {
3823 //--------------------------------------------------------------
3824 // Add 2D float attribute
3825 //--------------------------------------------------------------
3826 
3827  hsize_t dims1[1] = { 2 };
3828  H5::DataSpace attr_dataspace_2 = H5::DataSpace (1, dims1 );
3829  H5::Attribute attribute_float = hfilehandle.createAttribute( attrname, H5::PredType::NATIVE_FLOAT,
3830  attr_dataspace_2);
3831  attribute_float.write( H5::PredType::NATIVE_FLOAT, attrdata);
3832  attr_dataspace_2.close();
3833  attribute_float.close();
3834 
3835 }
3836 
3837 //--------------------------------------------------------------
3838 
3839 
3840 //--------------------------------------------------------------
3841 void Export_Files::Xport:: hfile_add_attr(H5::H5File &hfilehandle, string attrname,
3842  string attrdata) {
3843 //--------------------------------------------------------------
3844 // Add string attribute
3845 //--------------------------------------------------------------
3846 
3847  hsize_t dims1[1] = { 1 };
3848  H5::DataSpace attr_dataspace_1 = H5::DataSpace (1, dims1 );
3849  H5::StrType strdatatype(H5::PredType::C_S1, 256); // of length 256 characters
3850  H5std_string strwritebuf ( attrdata );
3851  H5::Attribute attribute_string = hfilehandle.createAttribute( attrname, strdatatype, attr_dataspace_1);
3852  attribute_string.write(strdatatype, strwritebuf);
3853  attr_dataspace_1.close();
3854  attribute_string.close();
3855 
3856 }
3857 
3858 //--------------------------------------------------------------
3859 
3860 
3861 //--------------------------------------------------------------
3862 void Export_Files::Xport:: hfile_add_attr_todataset(H5::DataSet &hdatasethandle, string attrname,
3863  string attrdata) {
3864 //--------------------------------------------------------------
3865 // Add string attribute to a specific dataset
3866 //--------------------------------------------------------------
3867 
3868  hsize_t dims1[1] = { 1 };
3869  H5::DataSpace attr_dataspace_1 = H5::DataSpace (1, dims1 );
3870  H5::StrType strdatatype(H5::PredType::C_S1, 256); // of length 256 characters
3871  H5std_string strwritebuf ( attrdata );
3872  H5::Attribute attribute_string = hdatasethandle.createAttribute( attrname, strdatatype, attr_dataspace_1);
3873  attribute_string.write(strdatatype, strwritebuf);
3874  attr_dataspace_1.close();
3875  attribute_string.close();
3876 
3877 }
3878 
3879 //-------------------------------------------------------------
3880 
3881 
3882 
3883 //**************************************************************
3884 //**************************************************************
Export_Files::oAxis::min
float min
Definition: export.h:117
Export_Files::Header::timeU
string timeU
Definition: export.h:165
Output_Data::p2p1x1_1D::dcell
Array2D< float > dcell(size_t s) const
Definition: export.h:354
Export_Files::Header::Time_label
string Time_label()
Definition: export.cpp:522
lib-algorithms.h
Algorithms
Export_Files::oAxis::units
string units
Definition: export.h:116
formulary.h
Plasma Formulary and Units - Declarations.
Output_Data::p1x1_1D::Pmin
float Pmin(size_t s) const
Definition: export.h:307
Output_Data::Output_Preprocessor_1D::f10
void f10(const State1D &Y, const Grid_Info &grid, const size_t tout, const Parallel_Environment_1D &PE)
Definition: export.cpp:1997
Output_Data::p1x1_1D::Pmax
float Pmax(size_t s) const
Definition: export.h:308
Array3D::array
valarray< T > & array() const
Definition: lib-array.h:1056
Field1D::numx
size_t numx() const
Definition: state.h:197
Output_Data::p1x1_1D::~p1x1_1D
~p1x1_1D()
Definition: export.cpp:1012
Output_Data::p2p1x1_1D::Species
size_t Species() const
Definition: export.h:340
Export_Files::Header::min
float min(const size_t i)
Definition: export.h:154
Export_Files::Makefolder
int Makefolder(string _name)
Definition: export.cpp:51
Output_Data::p2p1x1_1D::Pmax
float Pmax(size_t s) const
Definition: export.h:347
Export_Files::Header::label
string label(const size_t i)
Definition: export.cpp:513
EMF1D::By
Field1D & By()
Definition: state.h:294
Output_Data::Output_Preprocessor_1D::Qx
void Qx(const State1D &Y, const Grid_Info &grid, const size_t tout, const Parallel_Environment_1D &PE)
Definition: export.cpp:2607
vfloat_complex
vector< float > vfloat_complex(const vector< complex< double > > vDouble)
Definition: nmethods.cpp:412
input.h
Input reader - Declarations.
Export_Files::stringify
std::string stringify(T const &x)
Definition: export.h:30
Export_Files::Xport::hfile_add_attr
void hfile_add_attr(H5::H5File &hfilehandle, string attrname, int attrdata)
Definition: export.cpp:3761
Export_Files::DefaultTags
Definition: export.h:44
Export_Files::Xport::hinit_attr2
void hinit_attr2(H5::H5File &hfilehandle, const std::string tag, size_t step, float xmax[], float xmin[])
Definition: export.cpp:3698
Output_Data::Output_Preprocessor_1D::Bx
void Bx(const State1D &Y, const Grid_Info &grid, const size_t tout, const Parallel_Environment_1D &PE)
Definition: export.cpp:1571
Export_Files::Header::oDir
string oDir
Definition: export.h:167
lib-array.h
Underlying data structures.
vfloat
valarray< float > vfloat(const valarray< double > &vDouble)
Definition: nmethods.cpp:381
Algorithms::AxisBundle::pmin
T pmin(size_t i) const
Definition: lib-algorithms.h:114
Grid_Info::l0
const vector< size_t > l0
Definition: setup.h:39
Export_Files::Header::time
string time
Definition: export.h:165
Output_Data::p2p1x1_1D::Pmin
float Pmin(size_t s) const
Definition: export.h:346
Output_Data::Output_Preprocessor_1D::Jy
void Jy(const State1D &Y, const Grid_Info &grid, const size_t tout, const Parallel_Environment_1D &PE)
Definition: export.cpp:2493
Algorithms::Legendre
valarray< T > Legendre(const T x, const size_t Nl)
Definition: lib-algorithms.h:227
export.h
Export Files - Declarations.
Export_Files::Header::units
string units(const size_t i)
Definition: export.cpp:516
Export_Files::Restart_Facility::Read
void Read(const int rank, const size_t re_step, State1D &Y)
Definition: export.cpp:718
Export_Files::oAxis::label
string label
Definition: export.h:115
Output_Data::Output_Preprocessor_1D::Uy
void Uy(const State1D &Y, const Grid_Info &grid, const size_t tout, const Parallel_Environment_1D &PE)
Definition: export.cpp:3144
Algorithms::AxisBundle::Npy
size_t Npy(size_t i) const
Definition: lib-algorithms.h:128
Export_Files::DefaultTags::fvsx
vector< string > fvsx
Definition: export.h:52
Output_Data::p2p1x1_1D::~p2p1x1_1D
~p2p1x1_1D()
Definition: export.cpp:1227
Array3D::dim3
size_t dim3() const
Definition: lib-array.h:1055
Formulary::Uconv
double Uconv(string key)
Definition: formulary.h:78
Output_Data::Output_Preprocessor_1D::Qy
void Qy(const State1D &Y, const Grid_Info &grid, const size_t tout, const Parallel_Environment_1D &PE)
Definition: export.cpp:2713
DistFunc1D::mass
double mass() const
Definition: state.h:400
Export_Files::Restart_Facility::Restart_Facility
Restart_Facility(const int rank, string homedir="")
Definition: export.cpp:708
Array2D::dim1
size_t dim1() const
Definition: lib-array.h:289
Output_Data::Output_Preprocessor_1D::f0
void f0(const State1D &Y, const Grid_Info &grid, const size_t tout, const Parallel_Environment_1D &PE)
Definition: export.cpp:1862
Export_Files::oAxis
Definition: export.h:112
Export_Files::Xport::hinit_attr
void hinit_attr(H5::H5File &hfilehandle, const std::string tag, size_t step, float xmax, float xmin)
Definition: export.cpp:3658
Export_Files::Header::max
float max(const size_t i)
Definition: export.h:155
Output_Data::p1x1_1D::operator()
valarray< float > operator()(DistFunc1D &df, size_t x0, size_t s)
Definition: export.cpp:1020
Export_Files::Folders
void Folders()
Definition: export.cpp:74
Output_Data::Output_Preprocessor_1D::Jz
void Jz(const State1D &Y, const Grid_Info &grid, const size_t tout, const Parallel_Environment_1D &PE)
Definition: export.cpp:2551
Output_Data::Output_Preprocessor_1D::vNx
void vNx(const State1D &Y, const Grid_Info &grid, const size_t tout, const Parallel_Environment_1D &PE)
Definition: export.cpp:2926
Output_Data::fx1_1D::Pmax
float Pmax(size_t s) const
Definition: export.h:390
Export_Files::Xport::hfile_add_attr2
void hfile_add_attr2(H5::H5File &hfilehandle, string attrname, int attrdata[2])
Definition: export.cpp:3781
Export_Files::DefaultTags::mom
vector< string > mom
Definition: export.h:50
Export_Files::DefaultTags::pvsx
vector< string > pvsx
Definition: export.h:51
Array2D::dim2
size_t dim2() const
Definition: lib-array.h:290
Export_Files::DefaultTags::fld
vector< string > fld
Definition: export.h:49
Output_Data::Output_Preprocessor_1D::Qz
void Qz(const State1D &Y, const Grid_Info &grid, const size_t tout, const Parallel_Environment_1D &PE)
Definition: export.cpp:2819
Output_Data::p1x1_1D
Definition: export.h:290
Array3D::dim2
size_t dim2() const
Definition: lib-array.h:1054
Output_Data::fx1_1D::Np
size_t Np(size_t s) const
Definition: export.h:388
Export_Files::Header::titleC
float titleC
Definition: export.h:166
Grid_Info::m0
const vector< size_t > m0
Definition: setup.h:42
Algorithms::AxisBundle::xgmin
T xgmin(size_t i) const
Definition: lib-algorithms.h:107
Output_Data::p2p1x1_1D
Definition: export.h:327
Formulary::Label
string Label(string key)
Definition: formulary.h:76
Array3D::dim1
size_t dim1() const
Definition: lib-array.h:1053
Output_Data::PLegendre2D::PLegendre2D
PLegendre2D(size_t Nl, size_t Nm, size_t Np, float pmin, float pmax, size_t Npx, size_t Npy)
Definition: export.cpp:886
Export_Files::oAxis::oAxis
oAxis()
Definition: export.cpp:440
Output_Data::fx1_1D::fx1_1D
fx1_1D(const Grid_Info &_G)
Definition: export.cpp:1093
Hydro1D::vz
double & vz(size_t i)
Definition: state.h:536
Export_Files::Xport::hmake_file
H5::H5File hmake_file(string ofilename)
Definition: export.cpp:3632
Output_Data::Output_Preprocessor_1D::Uz
void Uz(const State1D &Y, const Grid_Info &grid, const size_t tout, const Parallel_Environment_1D &PE)
Definition: export.cpp:3196
DistFunc1D::dim
size_t dim() const
Definition: state.h:395
EMF1D::Ez
Field1D & Ez()
Definition: state.h:292
Algorithms::AxisBundle::xgmax
T xgmax(size_t i) const
Definition: lib-algorithms.h:108
Output_Data::Output_Preprocessor_1D::Ex
void Ex(const State1D &Y, const Grid_Info &grid, const size_t tout, const Parallel_Environment_1D &PE)
Definition: export.cpp:1421
Output_Data::Output_Preprocessor_1D::f20
void f20(const State1D &Y, const Grid_Info &grid, const size_t tout, const Parallel_Environment_1D &PE)
Definition: export.cpp:2140
Output_Data::Output_Preprocessor_1D::vNy
void vNy(const State1D &Y, const Grid_Info &grid, const size_t tout, const Parallel_Environment_1D &PE)
Definition: export.cpp:2985
Output_Data::Output_Preprocessor_1D::fl0
void fl0(const State1D &Y, const Grid_Info &grid, const size_t tout, const Parallel_Environment_1D &PE)
Definition: export.cpp:2212
Output_Data::Output_Preprocessor_1D::distdump
void distdump(const State1D &Y, const Grid_Info &grid, const size_t tout, const Parallel_Environment_1D &PE)
Definition: export.cpp:1389
State1D
Definition: state.h:577
Output_Data::Output_Preprocessor_1D::T
void T(const State1D &Y, const Grid_Info &grid, const size_t tout, const Parallel_Environment_1D &PE)
Definition: export.cpp:2345
Export_Files::Header::Header
Header()
Definition: export.h:138
EMF1D::Bx
Field1D & Bx()
Definition: state.h:293
Export_Files::ofconventions::rank_digits
const int rank_digits
Definition: export.h:22
State1D::HYDRO
Hydro1D & HYDRO()
Definition: state.h:613
Output_Data::fx1_1D::Species
size_t Species() const
Definition: export.h:386
Array2D::array
valarray< T > & array() const
Definition: lib-array.h:291
Hydro1D::temperature
double & temperature(size_t i)
Definition: state.h:539
Hydro1D::Z
double & Z(size_t i)
Definition: state.h:542
Algorithms::AxisBundle::Np
size_t Np(size_t i) const
Definition: lib-algorithms.h:113
Output_Data::Output_Preprocessor_1D::f11
void f11(const State1D &Y, const Grid_Info &grid, const size_t tout, const Parallel_Environment_1D &PE)
Definition: export.cpp:2068
Output_Data::Output_Preprocessor_1D::px
void px(const State1D &Y, const Grid_Info &grid, const size_t tout, const Parallel_Environment_1D &PE)
Definition: export.cpp:1720
Algorithms::MakeAxis
valarray< T > MakeAxis(const T min, const T max, const size_t N)
Definition: lib-algorithms.h:25
Parallel_Environment_1D::NODES
int NODES() const
Definition: parallel.cpp:1092
Output_Data::PLegendre2D
Definition: export.h:266
Export_Files::Header::timeC
float timeC
Definition: export.h:166
Output_Data::Output_Preprocessor_1D::operator()
void operator()(const State1D &Y, const Grid_Info &grid, const size_t tout, const Parallel_Environment_1D &PE)
Definition: export.cpp:1277
Export_Files::Xport::hclose_file
void hclose_file(H5::H5File &file)
Definition: export.cpp:3646
Output_Data::PLegendre1D
Definition: export.h:241
Export_Files::Xport::Xport
Xport(const Algorithms::AxisBundle< double > &_axis, const vector< string > oTags, string homedir="")
Definition: export.cpp:530
Parallel_Environment_1D::RANK
int RANK() const
Definition: parallel.cpp:1090
DistFunc1D::q
double q() const
Definition: state.h:399
Output_Data::fx1_1D
Definition: export.h:373
EMF1D::Bz
Field1D & Bz()
Definition: state.h:295
State1D::DF
DistFunc1D & DF(size_t s)
Definition: state.h:602
State1D::SH
SHarmonic1D & SH(size_t s, size_t lh, size_t mh)
Definition: state.h:605
setup.h
Grid Setup - Declaration.
Export_Files::Header::axis
valarray< float > axis(const size_t i)
Definition: export.cpp:509
Export_Files::ofconventions::rfile_digits
const int rfile_digits
Definition: export.h:21
Export_Files::Restart_Facility::Write
void Write(const int rank, const size_t re_step, State1D &Y)
Definition: export.cpp:756
Output_Data::PLegendre1D::~PLegendre1D
~PLegendre1D()
Definition: export.cpp:878
state.h
Fields, Distributions, Harmonics, States - Declarations.
EMF1D::Ex
Field1D & Ex()
Definition: state.h:290
Export_Files::Xport::haxis
void haxis(H5::Group &hgrouphandle, string axismainname, float axisrange[2], string axislongname, string axisname, string axistype, string axisunits)
Definition: export.cpp:3736
formulary
Formulary & formulary()
Definition: formulary.cpp:328
Output_Data::p2p1x1_1D::p2p1x1_1D
p2p1x1_1D(const Grid_Info &_G)
Definition: export.cpp:1156
Output_Data::Output_Preprocessor_1D::Ux
void Ux(const State1D &Y, const Grid_Info &grid, const size_t tout, const Parallel_Environment_1D &PE)
Definition: export.cpp:3092
Export_Files::Restart_Facility::rFextension
string rFextension(const int rank, const size_t rstep)
Definition: export.cpp:789
Output_Data::Output_Preprocessor_1D::Ez
void Ez(const State1D &Y, const Grid_Info &grid, const size_t tout, const Parallel_Environment_1D &PE)
Definition: export.cpp:1521
Output_Data::PLegendre2D::dim
size_t dim() const
Definition: export.h:275
Hydro1D::vx
double & vx(size_t i)
Definition: state.h:530
Grid_Info::axis
const Algorithms::AxisBundle< double > axis
Definition: setup.h:47
Grid_Info
Definition: setup.h:15
Export_Files::DefaultTags::pvspvsx
vector< string > pvspvsx
Definition: export.h:53
Output_Data::p2p1x1_1D::PL
PLegendre2D PL(size_t s) const
Definition: export.h:350
Export_Files::Header::Time_conv
float Time_conv()
Definition: export.cpp:523
Output_Data::fx1_1D::operator()
Array2D< float > operator()(DistFunc1D &df, size_t l, size_t m, size_t x0, size_t s)
Definition: export.cpp:1124
Output_Data::Output_Preprocessor_1D::n
void n(const State1D &Y, const Grid_Info &grid, const size_t tout, const Parallel_Environment_1D &PE)
Definition: export.cpp:2285
Export_Files::DefaultTags::time
vector< string > time
Definition: export.h:47
Output_Data::p2p1x1_1D::npcell
Array2D< size_t > npcell(size_t s) const
Definition: export.h:353
Output_Data::p2p1x1_1D::prad
Array2D< float > prad(size_t s) const
Definition: export.h:352
Algorithms::AxisBundle::Npx
size_t Npx(size_t i) const
Definition: lib-algorithms.h:121
Algorithms::MakeCAxis
valarray< T > MakeCAxis(const T min, const T max, const size_t N)
Definition: lib-algorithms.h:36
Input::List
Input_List & List()
Definition: input.cpp:1585
Export_Files::DefaultTags::space
vector< string > space
Definition: export.h:48
Export_Files::DefaultTags::DefaultTags
DefaultTags(size_t species)
Definition: export.cpp:336
Algorithms::AxisBundle::pmax
T pmax(size_t i) const
Definition: lib-algorithms.h:115
Output_Data::p1x1_1D::Species
size_t Species() const
Definition: export.h:303
Export_Files::Xport::hfile_add_attr_todataset
void hfile_add_attr_todataset(H5::DataSet &hdatasethandle, string attrname, string attrdata)
Definition: export.cpp:3862
Hydro1D::density
double & density(size_t i)
Definition: state.h:527
Parallel_Environment_1D
Definition: parallel.h:130
Export_Files::Header::Title_conv
float Title_conv()
Definition: export.cpp:521
Output_Data::PLegendre2D::~PLegendre2D
~PLegendre2D()
Definition: export.cpp:954
EMF1D::Ey
Field1D & Ey()
Definition: state.h:291
Algorithms::AxisBundle::Npz
size_t Npz(size_t i) const
Definition: lib-algorithms.h:135
Export_Files::Header::Title_label
string Title_label()
Definition: export.cpp:520
Output_Data::fx1_1D::Pmin
float Pmin(size_t s) const
Definition: export.h:389
Export_Files::oAxis::max
float max
Definition: export.h:117
Export_Files::Xport::oH5Fextension
string oH5Fextension(size_t step, int species=-1)
Definition: export.cpp:683
Output_Data::Output_Preprocessor_1D::vNz
void vNz(const State1D &Y, const Grid_Info &grid, const size_t tout, const Parallel_Environment_1D &PE)
Definition: export.cpp:3037
DistFunc1D::l0
size_t l0() const
Definition: state.h:396
Output_Data::Output_Preprocessor_1D::Ti
void Ti(const State1D &Y, const Grid_Info &grid, const size_t tout, const Parallel_Environment_1D &PE)
Definition: export.cpp:3360
Output_Data::p1x1_1D::p1x1_1D
p1x1_1D(const Grid_Info &_G)
Definition: export.cpp:962
Output_Data::p2p1x1_1D::p2p1_x1
Array2D< float > p2p1_x1(size_t s) const
Definition: export.h:351
Hydro1D::vy
double & vy(size_t i)
Definition: state.h:533
Output_Data::fx1_1D::~fx1_1D
~fx1_1D()
Definition: export.cpp:1117
Algorithms::AxisBundle::Nx
size_t Nx(size_t i) const
Definition: lib-algorithms.h:98
Output_Data::PLegendre1D::PLegendre1D
PLegendre1D(size_t Nl, size_t Nm, size_t Np, float pmin, float pmax, size_t Npx)
Definition: export.cpp:816
Output_Data::Output_Preprocessor_1D::Z
void Z(const State1D &Y, const Grid_Info &grid, const size_t tout, const Parallel_Environment_1D &PE)
Definition: export.cpp:3249
Algorithms::AxisBundle::Nxg
size_t Nxg(size_t i) const
Definition: lib-algorithms.h:106
Export_Files::oAxis::sz
size_t sz
Definition: export.h:118
Output_Data::PLegendre1D::dim
size_t dim() const
Definition: export.h:250
Export_Files::ofconventions::ofile_digits
const int ofile_digits
Definition: export.h:17
Export_Files::Xport::Export_h5
void Export_h5(const std::string tag, std::valarray< float > ex, const size_t &step, const int spec=-1)
Definition: export.cpp:3416
Output_Data::p2p1x1_1D::operator()
Array2D< float > operator()(DistFunc1D &df, size_t x0, size_t s)
Definition: export.cpp:1234
Output_Data::Output_Preprocessor_1D::By
void By(const State1D &Y, const Grid_Info &grid, const size_t tout, const Parallel_Environment_1D &PE)
Definition: export.cpp:1621
Export_Files::ofconventions::h5file_extension
const string h5file_extension
Definition: export.h:25
Input::Input_List::BoundaryCells
int BoundaryCells
Definition: input.h:50
Export_Files::Header::title
string title
Definition: export.h:165
Export_Files::Header::xyz_axis
vector< oAxis > xyz_axis
Definition: export.h:164
State1D::Species
size_t Species() const
Definition: state.h:596
Output_Data::Output_Preprocessor_1D::ni
void ni(const State1D &Y, const Grid_Info &grid, const size_t tout, const Parallel_Environment_1D &PE)
Definition: export.cpp:3305
Output_Data::Output_Preprocessor_1D::Jx
void Jx(const State1D &Y, const Grid_Info &grid, const size_t tout, const Parallel_Environment_1D &PE)
Definition: export.cpp:2434
Output_Data::p2p1x1_1D::Nm
size_t Nm(size_t s) const
Definition: export.h:342
Output_Data::Output_Preprocessor_1D::Ey
void Ey(const State1D &Y, const Grid_Info &grid, const size_t tout, const Parallel_Environment_1D &PE)
Definition: export.cpp:1471
Output_Data::p2p1x1_1D::Nl
size_t Nl(size_t s) const
Definition: export.h:341
DistFunc1D
The 1D distribution function is the container for all SHarmonic1D per species.
Definition: state.h:376
Algorithms::moment
T moment(const vector< T > q, const vector< T > x, const int p)
Definition: lib-algorithms.h:305
Export_Files::Header::Directory
string Directory()
Definition: export.cpp:524
Export_Files::Header::dim
size_t dim()
Definition: export.cpp:504
Output_Data::p1x1_1D::p1_x1
valarray< float > p1_x1(size_t s) const
Definition: export.h:311
Output_Data::Output_Preprocessor_1D::Bz
void Bz(const State1D &Y, const Grid_Info &grid, const size_t tout, const Parallel_Environment_1D &PE)
Definition: export.cpp:1671
State1D::EMF
EMF1D & EMF() const
Definition: state.h:610
Export_Files::ofconventions::rfile_extension
const string rfile_extension
Definition: export.h:23
nmethods.h
Numerical Methods - Declarations.