1 | /***************************************
2 | $Header$
3 |
4 | This file contains the functions
5 | +latex+{\tt IlluMultiSave()}, {\tt IlluMultiLoad()} and {\tt IlluMultiRead()}
6 | +html+ <tt>IlluMultiSave()</tt>, <tt>IlluMultiLoad()</tt> and
7 | +html+ <tt>IlluMultiRead()</tt>
8 | designed to handle distributed storage and retrieval of data on local drives
9 | of machines in a Beowulf cluster. This should allow rapid loading of
10 | timestep data for "playback" from what is essentially a giant RAID-0 array of
11 | distributed disks, at enormously higher speeds than via NFS from a hard drive
12 | or RAID array on the head node. The danger of course is that if one node's
13 | disk goes down, you don't have a valid data set any more, but that's the
14 | nature of RAID-0, right?
15 | +latex+\par
16 | +html+ <p>
17 | The filenames saved are:
18 | +latex+\begin{itemize} \item {\tt $<$basename$>$.cpu\#\#\#\#.meta}, where
19 | +latex+{\tt \#\#\#\#}
20 | +html+ <ul><li><tt><basename>.cpu####.meta</tt>, where <tt>####</tt>
21 | is replaced by the CPU number (more than four digits if there are more than
22 | 9999 CPUs :-), which has the metadata for the whole thing in XML format
23 | (written by
24 | +latex+{\tt GNOME libxml}),
25 | +html+ <tt>GNOME libxml</tt>),
26 | as described in the notes on the
27 | +latex+{\tt IlluMultiStoreXML()} function.
28 | +html+ <tt>IlluMultiStoreXML()</tt> function.
29 | +latex+\item {\tt $<$basename$>$.cpu\#\#\#\#.data}
30 | +html+ </li><li><tt><basename>.cpu####.data</tt>
31 | which is simply a stream of the raw data, optionally compressed by
32 | +latex+{\tt gzip}. \end{itemize}
33 | +html+ <tt>gzip</tt>.</li></ul>
34 | If one were saving timesteps, one might include a timestep number in the
35 | basename, and also timestep and simulation time in the metadata. The
36 | metadata can also hold simulation parameters, etc.
37 | +latex+\par
38 | +html+ <p>
39 | This supports 1-D, 2-D and 3-D distributed arrays. As an extra feature, you
40 | can load a multi-CPU distributed array scattered over lots of files into a
41 | single CPU, to facilitate certain modes of data visualization.
42 | ***************************************/
43 |
44 |
45 | #include "illuminator.h" /* Also includes petscda.h */
46 | #include <glib.h> /* for guint*, GUINT*_SWAP_LE_BE */
47 | #include <petscblaslapack.h> /* for BLAScopy_() */
48 | #include <libxml/tree.h> /* To build the XML tree to store */
49 | #include <libxml/parser.h> /* XML parsing */
50 | #include <stdio.h> /* FILE, etc. */
51 | #include <errno.h> /* errno */
52 | #include <string.h> /* strncpy(), etc. */
53 | #include <stdlib.h> /* malloc() and free() */
54 |
55 | /* Build with -DDEBUG for debugging output */
56 | #undef DPRINTF
57 | #ifdef DEBUG
58 | #define DPRINTF(fmt, args...) ierr = PetscPrintf (comm, "%s: " fmt, __FUNCT__, args); CHKERRQ (ierr)
59 | #else
60 | #define DPRINTF(fmt, args...)
61 | #endif
62 |
63 |
64 | #undef __FUNCT__
65 | #define __FUNCT__ "IlluMultiParseXML"
66 |
67 | /*++++++++++++++++++++++++++++++++++++++
68 | This function reads in the XML metadata document and returns the various
69 | parameter values in the addresses pointed to by the arguments. It is called
70 | by IlluMultiLoad() and IlluMultiRead().
71 |
72 | int IlluMultiParseXML It returns zero or an error code.
73 |
74 | MPI_Comm comm MPI communicator, if NULL it uses PETSC_COMM_WORLD.
75 |
76 | char *basename Base file name.
77 |
78 | int rank CPU number to read data for.
79 |
80 | int *compressed Data compression: if zero then no compression (fastest), 1-9
81 | then gzip compression level, 10-15 then gzip --best. If 16-31 then save
82 | guint32s representing relative values between min and max for each field,
83 | compressed according to this value minus 16. Likewise for 32-47 and
84 | guint16s, and 48-63 for guint8s. Yes, these alternative formats lose
85 | information and can't be used for accurate checkpointing, but they should
86 | retain enough data for visualization (except perhaps for the guint8s, which
87 | are possibly acceptable for vectors but likely not contours).
88 |
89 | int *wrongendian Tells whether the data are stored in the opposite endian
90 | format from this platform, and thus must be switched when the data are
91 | streamed in.
92 |
93 | int *dim Dimensionality of the space.
94 |
95 | int *px Number of processors in the
96 | +latex+$x$-direction.
97 | +html+ <i>x</i>-direction.
98 |
99 | int *py Number of processors in the
100 | +latex+$y$-direction.
101 | +html+ <i>y</i>-direction.
102 |
103 | int *pz Number of processors in the
104 | +latex+$z$-direction.
105 | +html+ <i>z</i>-direction.
106 |
107 | int *nx Number of grid points over the entire array in the
108 | +latex+$x$-direction.
109 | +html+ <i>x</i>-direction.
110 |
111 | int *ny Number of grid points over the entire array in the
112 | +latex+$y$-direction.
113 | +html+ <i>y</i>-direction.
114 |
115 | int *nz Number of grid points over the entire array in the
116 | +latex+$z$-direction.
117 | +html+ <i>z</i>-direction.
118 |
119 | PetscScalar *wx Physical overall width in the
120 | +latex+$x$-direction, {\tt PETSC\_NULL} if not needed.
121 | +html+ <i>x</i>-direction, <tt>PETSC_NULL</tt> if not needed.
122 |
123 | PetscScalar *wy Physical overall width in the
124 | +latex+$y$-direction, {\tt PETSC\_NULL} if not needed.
125 | +html+ <i>y</i>-direction, <tt>PETSC_NULL</tt> if not needed.
126 |
127 | PetscScalar *wz Physical overall width in the
128 | +latex+$z$-direction, {\tt PETSC\_NULL} if not needed.
129 | +html+ <i>z</i>-direction, <tt>PETSC_NULL</tt> if not needed.
130 |
131 | int *xm Number of grid points over the local part of the array in the
132 | +latex+$x$-direction.
133 | +html+ <i>x</i>-direction.
134 |
135 | int *ym Number of grid points over the local part of the array in the
136 | +latex+$y$-direction.
137 | +html+ <i>y</i>-direction.
138 |
139 | int *zm Number of grid points over the local part of the array in the
140 | +latex+$z$-direction.
141 | +html+ <i>z</i>-direction.
142 |
143 | int *dof Degrees of freedom at each node, a.k.a. number of field variables.
144 |
145 | int *sw Stencil width.
146 |
147 | DAStencilType *st Stencil type, given by the
148 | +latex+{\tt PETSc enum} values.
149 | +html+ <tt>PETSc enum</tt> values.
150 |
151 | DAPeriodicType *wrap Periodic type, given by the
152 | +latex+{\tt PETSc enum} values.
153 | +html+ <tt>PETSc enum</tt> values.
154 |
155 | char ***fieldnames Names of the field variables.
156 |
157 | field_plot_type **fieldtypes Data (plot) types for field variables,
158 | +latex+{\tt PETSC\_NULL} if not needed.
159 | +html+ <tt>PETSC_NULL</tt> if not needed.
160 |
161 | PetscScalar **fieldmin Minimum value of each field variable.
162 |
163 | PetscScalar **fieldmax Maximum value of each field variable.
164 |
165 | int *usermetacount Number of user metadata parameters.
166 |
167 | char ***usermetanames User metadata parameter names.
168 |
169 | char ***usermetadata User metadata parameter strings.
170 | ++++++++++++++++++++++++++++++++++++++*/
171 |
172 | static int IlluMultiParseXML
173 | (MPI_Comm comm, char *basename, int rank, int *compressed, int *wrongendian,
174 | int *dim, int *px,int *py,int *pz, int *nx,int *ny,int *nz,
175 | PetscScalar *wx,PetscScalar *wy,PetscScalar *wz, int *xm,int *ym,int *zm,
176 | int *dof, int *sw, DAStencilType *st, DAPeriodicType *wrap,
177 | char ***fieldnames, field_plot_type **fieldtypes, PetscScalar **fieldmin,
178 | PetscScalar **fieldmax,
179 | int *usermetacount, char ***usermetanames, char ***usermetadata)
180 | {
181 | xmlDocPtr thedoc;
182 | xmlNodePtr thenode;
183 | char filename [1000];
184 | xmlChar *buffer, *version;
185 | int field=0, ierr;
186 |
187 | if (snprintf (filename, 999, "%s.cpu%.4d.meta", basename, rank) > 999)
188 | SETERRQ1 (PETSC_ERR_ARG_SIZ, "Base file name too long (%d chars, 960 max)",
189 | strlen (basename));
190 |
191 | DPRINTF ("Parsing XML file %s\n", filename);
192 | thedoc = xmlParseFile (filename);
193 | if (thedoc == NULL)
194 | SETERRQ (PETSC_ERR_FILE_OPEN, "xmlParseFile returned NULL\n");
195 |
196 | version = xmlGetProp (thedoc -> children, BAD_CAST "version");
197 | if (strncmp ((char *) version, "0.1", 3) &&
198 | strncmp ((char *) version, "0.2", 3))
199 | {
200 | free (version);
201 | SETERRQ (PETSC_ERR_FILE_UNEXPECTED,
202 | "Version is neither 0.1 nor 0.2, can't parse\n");
203 | }
204 |
205 | *wrongendian = 0;
206 | buffer = xmlGetProp (thedoc -> children, BAD_CAST "endian");
207 | #ifdef WORDS_BIGENDIAN
208 | if (strncasecmp ((char *) buffer, "big", 3))
209 | *wrongendian = 1;
210 | #else
211 | if (!strncasecmp ((char *) buffer, "big", 3))
212 | *wrongendian = 1;
213 | #endif
214 | free (buffer);
215 |
216 | buffer = xmlGetProp (thedoc -> children, BAD_CAST "compression");
217 | /* Have to handle all of the various compressed string values */
218 | /* TODO: make this deal with "float" and "complex" */
219 | if (buffer[0] == 'l' || buffer[0] == 'L')
220 | {
221 | if (buffer[4] == 'n' || buffer[4] == 'N')
222 | *compressed = COMPRESS_INT_LONG | COMPRESS_GZIP_NONE;
223 | else if (buffer[4] >= '1' && buffer[4] <= '9')
224 | {
225 | sscanf ((char *) buffer+4, "%d", compressed);
226 | *compressed |= COMPRESS_INT_LONG;
227 | }
228 | else if (buffer[4] == 'b' || buffer[4] == 'B')
229 | *compressed = COMPRESS_INT_LONG | COMPRESS_GZIP_BEST;
230 | else
231 | SETERRQ1 (PETSC_ERR_ARG_OUTOFRANGE,
232 | "Unrecognized compression type %s in .meta file", buffer);
233 | }
234 | else if (buffer[0] == 's' || buffer[0] == 'S')
235 | {
236 | if (buffer[5] == 'n' || buffer[5] == 'N')
237 | *compressed = COMPRESS_INT_SHORT | COMPRESS_GZIP_NONE;
238 | else if (buffer[5] >= '1' && buffer[5] <= '9')
239 | {
240 | sscanf ((char *) buffer+5, "%d", compressed);
241 | *compressed |= COMPRESS_INT_SHORT;
242 | }
243 | else if (buffer[5] == 'b' || buffer[5] == 'B')
244 | *compressed = COMPRESS_INT_SHORT | COMPRESS_GZIP_BEST;
245 | else
246 | SETERRQ1 (PETSC_ERR_ARG_OUTOFRANGE,
247 | "Unrecognized compression type %s in .meta file", buffer);
248 | }
249 | else if (buffer[0] == 'c' || buffer[0] == 'C')
250 | {
251 | if (buffer[4] == 'n' || buffer[4] == 'N')
252 | *compressed = COMPRESS_INT_CHAR | COMPRESS_GZIP_NONE;
253 | else if (buffer[4] >= '1' && buffer[4] <= '9')
254 | {
255 | sscanf ((char *) buffer+4, "%d", compressed);
256 | *compressed |= COMPRESS_INT_CHAR;
257 | }
258 | else if (buffer[4] == 'b' || buffer[4] == 'B')
259 | *compressed = COMPRESS_INT_CHAR | COMPRESS_GZIP_BEST;
260 | else
261 | SETERRQ1 (PETSC_ERR_ARG_OUTOFRANGE,
262 | "Unrecognized compression type %s in .meta file", buffer);
263 | }
264 | else
265 | {
266 | if (buffer[0] == 'n' || buffer[0] == 'N')
267 | *compressed = COMPRESS_INT_NONE | COMPRESS_GZIP_NONE;
268 | else if (buffer[0] >= '1' && buffer[0] <= '9')
269 | {
270 | sscanf ((char *) buffer, "%d", compressed);
271 | *compressed |= COMPRESS_INT_NONE;
272 | }
273 | else if (buffer[0] == 'b' || buffer[0] == 'B')
274 | *compressed = COMPRESS_INT_NONE | COMPRESS_GZIP_BEST;
275 | else
276 | SETERRQ1 (PETSC_ERR_ARG_OUTOFRANGE,
277 | "Unrecognized compression type %s in .meta file", buffer);
278 | }
279 | free (buffer);
280 | DPRINTF ("Root node: version %s, %s endian, compression %d|%d\n",
281 | version, *wrongendian ? "switched" : "native",
282 | *compressed & COMPRESS_INT_MASK, *compressed & COMPRESS_GZIP_MASK);
283 |
284 | *fieldnames = PETSC_NULL;
285 | if (fieldtypes != PETSC_NULL)
286 | *fieldtypes = PETSC_NULL;
287 | *fieldmin = *fieldmax = PETSC_NULL;
288 | *usermetacount = 0;
289 | *usermetanames = *usermetadata = PETSC_NULL;
290 |
291 | /* The main parsing loop; because xmlStringGetNodeList() doesn't work. */
292 | for (thenode = thedoc->children->xmlChildrenNode;
293 | thenode;
294 | thenode = thenode->next)
295 | {
296 | DPRINTF ("Looping through nodes, this is %s, fieldnames 0x%lx\n",
297 | thenode->name, *fieldnames);
298 | if (!strncasecmp ((char *) thenode->name, "GlobalCPUs", 10))
299 | {
300 | buffer = xmlGetProp (thenode, BAD_CAST "dimensions");
301 | sscanf ((char *) buffer, "%d", dim);
302 | free (buffer);
303 |
304 | buffer = xmlGetProp (thenode, BAD_CAST "xwidth");
305 | sscanf ((char *) buffer, "%d", px);
306 | free (buffer);
307 |
308 | if (*dim > 1)
309 | {
310 | buffer = xmlGetProp (thenode, BAD_CAST "ywidth");
311 | sscanf ((char *) buffer, "%d", py);
312 | free (buffer);
313 |
314 | if (*dim > 2)
315 | {
316 | buffer = xmlGetProp (thenode, BAD_CAST "zwidth");
317 | sscanf ((char *) buffer, "%d", pz);
318 | free (buffer);
319 | }
320 | else
321 | *pz = 1;
322 | }
323 | else
324 | *py = 1;
325 | }
326 |
327 | else if (!strncasecmp ((char *) thenode->name, "GlobalSize", 10))
328 | {
329 | buffer = xmlGetProp (thenode, BAD_CAST "xwidth");
330 | sscanf ((char *) buffer, "%d", nx);
331 | free (buffer);
332 | /*+ For GlobalSize, since there's no *size attribute (for an 0.1
333 | version document), assume 1. +*/
334 | buffer = xmlGetProp (thenode, BAD_CAST "xsize");
335 | if (wx != PETSC_NULL)
336 | {
337 | if (buffer)
338 | #ifdef PETSC_USE_SINGLE
339 | sscanf ((char *) buffer, "%f", wx);
340 | #else
341 | sscanf ((char *) buffer, "%lf", wx);
342 | #endif
343 | else
344 | *wx = 1.;
345 | }
346 | if (buffer)
347 | free (buffer);
348 |
349 | if (*dim > 1)
350 | {
351 | buffer = xmlGetProp (thenode, BAD_CAST "ywidth");
352 | sscanf ((char *) buffer, "%d", ny);
353 | free (buffer);
354 | buffer = xmlGetProp (thenode, BAD_CAST "ysize");
355 | if (wy != PETSC_NULL)
356 | {
357 | if (buffer)
358 | #ifdef PETSC_USE_SINGLE
359 | sscanf ((char *) buffer, "%f", wy);
360 | #else
361 | sscanf ((char *) buffer, "%lf", wy);
362 | #endif
363 | else
364 | *wy = 1.;
365 | }
366 | if (buffer)
367 | free (buffer);
368 |
369 | if (*dim > 2)
370 | {
371 | buffer = xmlGetProp (thenode, BAD_CAST "zwidth");
372 | sscanf ((char *) buffer, "%d", nz);
373 | free (buffer);
374 | buffer = xmlGetProp (thenode, BAD_CAST "zsize");
375 | if (wz != PETSC_NULL)
376 | {
377 | if (buffer)
378 | #ifdef PETSC_USE_SINGLE
379 | sscanf ((char *) buffer, "%f", wz);
380 | #else
381 | sscanf ((char *) buffer, "%lf", wz);
382 | #endif
383 | else
384 | *wz = 1.;
385 | }
386 | if (buffer)
387 | free (buffer);
388 | }
389 | else
390 | *nz = 1;
391 | }
392 | else
393 | *ny = *nz = 1;
394 |
395 | buffer = xmlGetProp (thenode, BAD_CAST "fields");
396 | sscanf ((char *) buffer, "%d", dof);
397 | free (buffer);
398 |
399 | if (*dof > field)
400 | {
401 | if (!(*fieldnames = (char **) realloc
402 | (*fieldnames, *dof * sizeof (char *))))
403 | SETERRQ (PETSC_ERR_MEM, "Can't allocate field names");
404 | if (fieldtypes != PETSC_NULL)
405 | if (!(*fieldtypes = (field_plot_type *) realloc
406 | (*fieldtypes, *dof * sizeof(field_plot_type))))
407 | SETERRQ (PETSC_ERR_MEM, "Can't allocate field types");
408 | if (!(*fieldmin = (PetscScalar *) realloc
409 | (*fieldmin, *dof * sizeof(PetscScalar))))
410 | SETERRQ (PETSC_ERR_MEM, "Can't allocate field params");
411 | if (!(*fieldmax = (PetscScalar *) realloc
412 | (*fieldmax, *dof * sizeof(PetscScalar))))
413 | SETERRQ (PETSC_ERR_MEM, "Can't allocate field params");
414 | }
415 | }
416 |
417 | else if (!strncasecmp ((char *) thenode->name, "LocalSize", 9))
418 | {
419 | buffer = xmlGetProp (thenode, BAD_CAST "xwidth");
420 | sscanf ((char *) buffer, "%d", xm);
421 | free (buffer);
422 |
423 | if (*dim > 1)
424 | {
425 | buffer = xmlGetProp (thenode, BAD_CAST "ywidth");
426 | sscanf ((char *) buffer, "%d", ym);
427 | free (buffer);
428 |
429 | if (*dim > 2)
430 | {
431 | buffer = xmlGetProp (thenode, BAD_CAST "zwidth");
432 | sscanf ((char *) buffer, "%d", zm);
433 | free (buffer);
434 | }
435 | else
436 | *zm = 1;
437 | }
438 | else
439 | *ym = 1;
440 | }
441 |
442 | else if (!strncasecmp ((char *) thenode->name, "Stencil", 7))
443 | {
444 | buffer = xmlGetProp (thenode, BAD_CAST "width");
445 | sscanf ((char *) buffer, "%d", sw);
446 | free (buffer);
447 |
448 | buffer = xmlGetProp (thenode, BAD_CAST "type");
449 | sscanf ((char *) buffer, "%d", st);
450 | free (buffer);
451 |
452 | buffer = xmlGetProp (thenode, BAD_CAST "periodic");
453 | sscanf ((char *) buffer, "%d", wrap);
454 | free (buffer);
455 | }
456 |
457 | else if (!strncasecmp ((char *) thenode->name, "Field", 5))
458 | {
459 | if (!*fieldnames || !*fieldmax || !*fieldmin)
460 | {
461 | printf ("Warning: reading a Field, but the number of them is unknown!\n");
462 | }
463 |
464 | if (field >= *dof)
465 | {
466 | printf ("Warning: more <Field> tags than declared degrees of freedom.\nThis may be because tags are out of order.\n");
467 | *fieldnames = (char **) realloc
468 | (*fieldnames, (field+1) * sizeof (char *));
469 | if (fieldtypes != PETSC_NULL)
470 | *fieldtypes = (field_plot_type *) realloc
471 | (*fieldtypes, (field+1) * sizeof (field_plot_type));
472 | *fieldmin = (PetscScalar *) realloc
473 | (*fieldmin, (field+1) * sizeof (PetscScalar));
474 | *fieldmax = (PetscScalar *) realloc
475 | (*fieldmax, (field+1) * sizeof (PetscScalar));
476 | }
477 |
478 | (*fieldnames) [field] =
479 | (char *) xmlGetProp (thenode, BAD_CAST "name");
480 | /*+ If the type attribute is missing from the Field node (as it is in
481 | version 0.1 documents), assume
482 | +latex+{\tt FIELD\_SCALAR}.
483 | +html+ <tt>FIELD_SCALAR</tt>. +*/
484 | buffer = xmlGetProp (thenode, BAD_CAST "type");
485 | if (fieldtypes)
486 | {
487 | if (buffer)
488 | sscanf ((char *) buffer, "0x%x", *fieldtypes + field);
489 | else
490 | (*fieldtypes) [field] = FIELD_SCALAR;
491 | }
492 | if (buffer)
493 | free(buffer);
494 | buffer = xmlGetProp (thenode, BAD_CAST "min");
495 | #ifdef PETSC_USE_SINGLE
496 | sscanf ((char *) buffer, "%f", *fieldmin + field);
497 | #else
498 | sscanf ((char *) buffer, "%lf", *fieldmin + field);
499 | #endif
500 | free (buffer);
501 | buffer = xmlGetProp (thenode, BAD_CAST "max");
502 | #ifdef PETSC_USE_SINGLE
503 | sscanf ((char *) buffer, "%f", *fieldmax + field);
504 | #else
505 | sscanf ((char *) buffer, "%lf", *fieldmax + field);
506 | #endif
507 | free (buffer);
508 | field++;
509 | }
510 |
511 | else if (!strncasecmp ((char *) thenode->name, "User", 4))
512 | {
513 | (*usermetacount)++;
514 | (*usermetanames) = (char **) realloc
515 | (*usermetanames, (*usermetacount) * sizeof (char *));
516 | (*usermetadata) = (char **) realloc
517 | (*usermetadata, (*usermetacount) * sizeof (char *));
518 | (*usermetanames) [(*usermetacount)-1] =
519 | (char *) xmlGetProp (thenode, BAD_CAST "name");
520 | (*usermetadata) [(*usermetacount)-1] =
521 | (char *) xmlGetProp (thenode, BAD_CAST "value");
522 | }
523 | }
524 |
525 | /* This is another way to do things, which would be a whole lot easier, if it
526 | worked... (See Debian's libxml bug list.)
527 | thenode = xmlStringGetNodeList (thedoc, "GlobalCPUs");
528 | printf ("thenode = 0x%lx\n", thenode);
529 | if (!thenode)
530 | return 1;
531 | thenode = xmlStringGetNodeList (thedoc, "GlobalSize");
532 | if (!thenode)
533 | return 1;
534 | thenode = xmlStringGetNodeList (thedoc, "LocalSize");
535 | if (!thenode)
536 | return 1;
537 | thenode = xmlStringGetNodeList (thedoc, "Stencil");
538 | if (!thenode)
539 | return 1;
540 |
541 | i = 0;
542 | thenode = xmlStringGetNodeList (thedoc, "Field");
543 | if (!thenode)
544 | return 1;
545 | while (thenode)
546 | {
547 | i++;
548 | thenode = thenode -> next;
549 | }
550 |
551 | thenode = xmlStringGetNodeList (thedoc, "User");
552 | while (thenode)
553 | {
554 | thenode = thenode -> next;
555 | } */
556 |
557 | free (version);
558 | xmlFreeDoc (thedoc);
559 |
560 | return 0;
561 | }
562 |
563 |
564 | #undef __FUNCT__
565 | #define __FUNCT__ "IlluMultiParseData"
566 |
567 | /*++++++++++++++++++++++++++++++++++++++
568 | This function reads in the data stored by IlluMultiStoreData(), complete with
569 | int/gzip compression.
570 |
571 | int IlluMultiParseData It returns zero or an error code.
572 |
573 | MPI_Comm comm MPI communicator, if NULL it uses PETSC_COMM_WORLD.
574 |
575 | PetscScalar *globalarray Array into which to load the (local) data.
576 |
577 | char *basename Base file name.
578 |
579 | int rank CPU number to read data for.
580 |
581 | int compressed Data compression: if zero then no compression (fastest), 1-9
582 | then gzip compression level, 10-15 then gzip --best. If 16-31 then save
583 | guint32s representing relative values between min and max for each field,
584 | compressed according to this value minus 16. Likewise for 32-47 and
585 | guint16s, and 48-63 for guint8s. Yes, these alternative formats lose
586 | information and can't be used for accurate checkpointing, but they should
587 | retain enough data for visualization (except perhaps for the guint8s, which
588 | are possibly acceptable for vectors but likely not contours).
589 |
590 | int gridpoints Number of gridpoints to read data for.
591 |
592 | int dof Degrees of freedom at each node, a.k.a. number of field variables.
593 |
594 | int wrongendian Tells whether the data are stored in the opposite endian
595 | format from this platform, and thus must be switched when the data are
596 | streamed in.
597 |
598 | PetscScalar *fieldmin Minimum value of each field variable.
599 |
600 | PetscScalar *fieldmax Maximum value of each field variable.
601 | ++++++++++++++++++++++++++++++++++++++*/
602 |
603 | static int IlluMultiParseData
604 | (MPI_Comm comm, PetscScalar *globalarray, char *basename, int rank,
605 | int compressed, int gridpoints, int dof, int wrongendian,
606 | PetscScalar *fieldmin, PetscScalar *fieldmax)
607 | {
608 | int i;
609 | char filename [1000];
610 | FILE *infile;
611 | size_t readoutput;
612 |
613 | if (compressed & COMPRESS_GZIP_MASK)
614 | {
615 | if (snprintf (filename, 999, "gunzip -c < %s.cpu%.4d.data",
616 | basename, rank) > 999)
617 | SETERRQ1 (PETSC_ERR_ARG_SIZ,
618 | "Base file name too long (%d chars, 960 max)",
619 | strlen (basename));
620 | if (!(infile = popen (filename, "r")))
621 | SETERRQ4 (PETSC_ERR_FILE_OPEN,
622 | "CPU %d: error %d (%s) opening input pipe \"%s\"",
623 | rank, errno, strerror (errno), filename);
624 | }
625 | else
626 | {
627 | if (snprintf (filename, 999, "%s.cpu%.4d.data", basename, rank) > 999)
628 | SETERRQ1 (PETSC_ERR_ARG_SIZ,
629 | "Base file name too long (%d chars, 960 max)",
630 | strlen (basename));
631 | if (!(infile = fopen (filename, "r")))
632 | SETERRQ4 (PETSC_ERR_FILE_OPEN,
633 | "CPU %d: error %d (%s) opening input file \"%s\"",
634 | rank, errno, strerror (errno), filename);
635 | }
636 |
637 | /* Read and store the data */
638 | /* TODO: make this deal with float and complex */
639 | switch (compressed & COMPRESS_INT_MASK)
640 | {
641 | /* Yes, this way of doing it is a bit redundant, but it seems better than
642 | the multitude of subconditionals which would be required if I did it
643 | the other way. */
644 | case COMPRESS_INT_LONG:
645 | {
646 | guint32 *storage;
647 | if (!(storage = (guint32 *) malloc
648 | (gridpoints * dof * sizeof (guint32))))
649 | SETERRQ (PETSC_ERR_MEM,
650 | "Can't allocate data decompression buffer");
651 | readoutput = fread (storage, sizeof (guint32), gridpoints*dof, infile);
652 |
653 | if (wrongendian)
654 | {
655 | for (i=0; i<gridpoints*dof; i++)
656 | storage[i] = GUINT32_SWAP_LE_BE_CONSTANT (storage[i]);
657 | }
658 | for (i=0; i<gridpoints*dof; i++)
659 | globalarray[i] = fieldmin [i%dof] +
660 | ((fieldmax [i%dof] - fieldmin [i%dof]) * storage[i]) / 4294967295.;
661 | free (storage);
662 | break;
663 | }
664 | case COMPRESS_INT_SHORT:
665 | {
666 | guint16 *storage;
667 | if (!(storage = (guint16 *) malloc
668 | (gridpoints * dof * sizeof (guint16))))
669 | SETERRQ (PETSC_ERR_MEM,
670 | "Can't allocate data decompression buffer");
671 | readoutput = fread (storage, sizeof (guint16), gridpoints*dof, infile);
672 |
673 | if (wrongendian)
674 | {
675 | for (i=0; i<gridpoints*dof; i++)
676 | storage[i] = GUINT16_SWAP_LE_BE_CONSTANT (storage[i]);
677 | }
678 | for (i=0; i<gridpoints*dof; i++)
679 | globalarray[i] = fieldmin [i%dof] +
680 | ((fieldmax [i%dof] - fieldmin [i%dof]) * storage[i]) / 65535.;
681 | free (storage);
682 | break;
683 | }
684 | case COMPRESS_INT_CHAR:
685 | {
686 | guint8 *storage;
687 | if (!(storage = (guint8 *) malloc
688 | (gridpoints * dof * sizeof (guint8))))
689 | SETERRQ (PETSC_ERR_MEM,
690 | "Can't allocate data decompression buffer");
691 | readoutput = fread (storage, sizeof (guint8), gridpoints*dof, infile);
692 |
693 | for (i=0; i<gridpoints*dof; i++)
694 | globalarray[i] = fieldmin [i%dof] +
695 | ((fieldmax [i%dof] - fieldmin [i%dof]) * storage[i]) / 255.;
696 | free (storage);
697 | break;
698 | }
699 | default:
700 | {
701 | /* TODO: check for different size data, complex */
702 | readoutput = fread (globalarray, sizeof (PetscScalar), gridpoints*dof,
703 | infile);
704 |
705 | if (wrongendian)
706 | {
707 | if (sizeof (PetscScalar) == 8)
708 | {
709 | for (i=0; i<gridpoints*dof; i++)
710 | globalarray[i]= GUINT64_SWAP_LE_BE_CONSTANT (globalarray[i]);
711 | }
712 | else if (sizeof (PetscScalar) == 4)
713 | {
714 | for (i=0; i<gridpoints*dof; i++)
715 | globalarray[i]= GUINT32_SWAP_LE_BE_CONSTANT (globalarray[i]);
716 | }
717 | else
718 | SETERRQ1 (PETSC_ERR_ARG_UNKNOWN_TYPE,
719 | "Unknown PetscScalar size %d (4-8 supported)",
720 | sizeof (PetscScalar));
721 | }
722 |
723 | }
724 | }
725 |
726 | if (compressed & COMPRESS_GZIP_MASK)
727 | pclose (infile);
728 | else
729 | fclose (infile);
730 |
731 | if (readoutput < gridpoints*dof)
732 | {
733 | int readerr = ferror (infile);
734 | SETERRQ5 (PETSC_ERR_FILE_READ,
735 | "CPU %d: error %d (%s) during read, got %d of %d items",
736 | rank, readerr, strerror (readerr), readoutput, gridpoints*dof);
737 | }
738 |
739 | return 0;
740 | }
741 |
742 |
743 | #undef __FUNCT__
744 | #define __FUNCT__ "IlluMultiStoreXML"
745 |
746 | /*++++++++++++++++++++++++++++++++++++++
747 | This function opens, stores and closes the XML metadata file for IlluMulti
748 | format storage. It is called by IlluMultiSave().
749 |
750 | int IlluMultiStoreXML It returns zero or an error code.
751 |
752 | MPI_Comm comm MPI communicator, if NULL it uses PETSC_COMM_WORLD.
753 |
754 | char *basename Base file name.
755 |
756 | int rank CPU number to store data for.
757 |
758 | int compressed Data compression: if zero then no compression (fastest), 1-9
759 | then gzip compression level, 10-15 then gzip --best. If 16-31 then save
760 | guint32s representing relative values between min and max for each field,
761 | compressed according to this value minus 16. Likewise for 32-47 and
762 | guint16s, and 48-63 for guint8s. Yes, these alternative formats lose
763 | information and can't be used for accurate checkpointing, but they should
764 | retain enough data for visualization (except perhaps for the guint8s, which
765 | are possibly acceptable for vectors but certainly not contours).
766 |
767 | int dim Dimensionality of the space.
768 |
769 | int px Number of processors in the
770 | +latex+$x$-direction.
771 | +html+ <i>x</i>-direction.
772 |
773 | int py Number of processors in the
774 | +latex+$y$-direction.
775 | +html+ <i>y</i>-direction.
776 |
777 | int pz Number of processors in the
778 | +latex+$z$-direction.
779 | +html+ <i>z</i>-direction.
780 |
781 | int nx Number of grid points over the entire array in the
782 | +latex+$x$-direction.
783 | +html+ <i>x</i>-direction.
784 |
785 | int ny Number of grid points over the entire array in the
786 | +latex+$y$-direction.
787 | +html+ <i>y</i>-direction.
788 |
789 | int nz Number of grid points over the entire array in the
790 | +latex+$z$-direction.
791 | +html+ <i>z</i>-direction.
792 |
793 | PetscScalar wx Physical overall width in the
794 | +latex+$x$-direction.
795 | +html+ <i>x</i>-direction.
796 |
797 | PetscScalar wy Physical overall width in the
798 | +latex+$y$-direction.
799 | +html+ <i>y</i>-direction.
800 |
801 | PetscScalar wz Physical overall width in the
802 | +latex+$z$-direction.
803 | +html+ <i>z</i>-direction.
804 |
805 | int xm Number of grid points over the local part of the array in the
806 | +latex+$x$-direction.
807 | +html+ <i>x</i>-direction.
808 |
809 | int ym Number of grid points over the local part of the array in the
810 | +latex+$y$-direction.
811 | +html+ <i>y</i>-direction.
812 |
813 | int zm Number of grid points over the local part of the array in the
814 | +latex+$z$-direction.
815 | +html+ <i>z</i>-direction.
816 |
817 | int dof Degrees of freedom at each node, a.k.a. number of field variables.
818 |
819 | int sw Stencil width.
820 |
821 | int st Stencil type, given by the
822 | +latex+{\tt PETSc enum} values.
823 | +html+ <tt>PETSc enum</tt> values.
824 |
825 | int wrap Periodic type, given by the
826 | +latex+{\tt PETSc enum} values.
827 | +html+ <tt>PETSc enum</tt> values.
828 |
829 | char **fieldnames Names of the field variables.
830 |
831 | field_plot_type *fieldtypes Data (plot) types for field variables.
832 |
833 | PetscReal *fieldmin Minimum value of each field variable.
834 |
835 | PetscReal *fieldmax Maximum value of each field variable.
836 |
837 | int usermetacount Number of user metadata parameters.
838 |
839 | char **usermetanames User metadata parameter names.
840 |
841 | char **usermetadata User metadata parameter strings.
842 | ++++++++++++++++++++++++++++++++++++++*/
843 |
844 | static int IlluMultiStoreXML
845 | (MPI_Comm comm, char *basename, int rank, int compressed,
846 | int dim, int px,int py,int pz, int nx,int ny,int nz,
847 | PetscScalar wx,PetscScalar wy,PetscScalar wz, int xm,int ym,int zm,
848 | int dof, int sw, int st, int wrap, char **fieldnames,
849 | field_plot_type *fieldtypes, PetscReal *fieldmin, PetscReal *fieldmax,
850 | int usermetacount, char **usermetanames, char **usermetadata)
851 | {
852 | xmlDocPtr thedoc;
853 | xmlNodePtr thenode;
854 | FILE *outfile;
855 | char filename [1000];
856 | xmlChar buffer [1000];
857 | int i;
858 |
859 | /*+The XML tags in the .meta file consist of:
860 | +latex+\begin{center} \begin{tabular}{|l|l|} \hline
861 | +html+ <center><table BORDER>
862 | +latex+{\tt IlluMulti} & Primary tag, with attributes {\tt version},
863 | +latex+\\ & {\tt endian} ({\tt big} or {\tt little}) and
864 | +latex+{\tt compression} \\ & (none, 1-9, best, float*, long*, short* or
865 | +latex+char*\footnote{longnone, long1-long9 or longbest compresses each
866 | +latex+double to a 32-bit unsigned int, with 0 representing the minimum for
867 | +latex+that field and 2$^{32}-1$ representing the maximum; likewise
868 | +latex+shortnone, short1-short9, shortbest uses 16-bit unsigned ints, and
869 | +latex+char* uses 8-bit unsigned ints. float is there to indicate that
870 | +latex+PetscScalar is 4 bytes at save time, loading should adjust
871 | +latex+accordingly; that's not fully supported just yet. At some point
872 | +latex+I'll have to figure out how to represent complex.}). \\ \hline
873 | +html+ <tr><td><tt>IlluMulti</tt></td><td>Primary tag, with attributes
874 | +html+ <tt>version</tt>,<br><tt>endian</tt> (<tt>big</tt> or
875 | +html+ <tt>little</tt>) and <tt>compression</tt><br>(none, 1-9, best;
876 | +html+ long*, short* or char*)</td></tr>
877 | +*/
878 |
879 | thedoc = xmlNewDoc (BAD_CAST "1.0");
880 | thedoc->children = xmlNewDocNode (thedoc,NULL, BAD_CAST"IlluMulti", NULL);
881 | xmlSetProp (thedoc->children, BAD_CAST "version", BAD_CAST "0.2");
882 | #ifdef WORDS_BIGENDIAN
883 | xmlSetProp (thedoc->children, BAD_CAST "endian", BAD_CAST "big");
884 | #else
885 | xmlSetProp (thedoc->children, BAD_CAST "endian", BAD_CAST "little");
886 | #endif
887 | if ((compressed & COMPRESS_INT_MASK) == COMPRESS_INT_LONG)
888 | strcpy ((char *) buffer, "long");
889 | else if ((compressed & COMPRESS_INT_MASK) == COMPRESS_INT_SHORT)
890 | strcpy ((char *) buffer, "short");
891 | else if ((compressed & COMPRESS_INT_MASK) == COMPRESS_INT_CHAR)
892 | strcpy ((char *) buffer, "char");
893 | /* Note: this doesn't really support complex types yet... */
894 | else if (sizeof (PetscScalar) == 4)
895 | strcpy ((char *) buffer, "float");
896 | else
897 | *buffer = '\0';
898 | if ((compressed & COMPRESS_GZIP_MASK) == 0)
899 | strcat ((char *) buffer, "none");
900 | else if ((compressed & COMPRESS_GZIP_MASK) >= 1 &&
901 | (compressed & COMPRESS_GZIP_MASK) <= 9)
902 | sprintf ((char *) buffer + strlen ((char *) buffer), "%d",
903 | compressed & COMPRESS_GZIP_MASK);
904 | else
905 | strcat ((char *) buffer, "best");
906 | xmlSetProp (thedoc->children, BAD_CAST "compression", buffer);
907 |
908 | /*+
909 | +latex+{\tt GlobalCPUs} & Number of CPUs in each direction, with
910 | +latex+\\ & attributes {\tt dimensions}, {\tt xwidth}, {\tt ywidth} and
911 | +latex+{\tt zwidth} \\ \hline
912 | +html+ <tr><td><tt>GlobalCPUs</tt></td><td>Number of CPUs in each
913 | +html+ direction, with<br>attributes <tt>dimensions</tt>, <tt>xwidth</tt>,
914 | +html+ <tt>ywidth</tt> and <tt>zwidth</tt></td></tr>
915 | +*/
916 |
917 | thenode = xmlNewChild (thedoc->children, NULL, BAD_CAST "GlobalCPUs", NULL);
918 | snprintf ((char *) buffer, 999, "%d", dim);
919 | xmlSetProp (thenode, BAD_CAST "dimensions", buffer);
920 | snprintf ((char *) buffer, 999, "%d", px);
921 | xmlSetProp (thenode, BAD_CAST "xwidth", buffer);
922 | if (dim > 1)
923 | {
924 | snprintf ((char *) buffer, 999, "%d", py);
925 | xmlSetProp (thenode, BAD_CAST "ywidth", buffer);
926 | if (dim > 2)
927 | {
928 | snprintf ((char *) buffer, 999, "%d", pz);
929 | xmlSetProp (thenode, BAD_CAST "zwidth", buffer);
930 | }
931 | }
932 |
933 | /*+
934 | +latex+{\tt GlobalSize} & Size of the entire distributed array, with
935 | +latex+\\ & attributes {\tt xwidth}, {\tt ywidth},
936 | +latex+{\tt zwidth}, {\tt xsize}**, {\tt ysize}**, {\tt zsize}** and
937 | +latex+{\tt fields} \\ \hline
938 | +html+ <tr><td><tt>GlobalSize</tt></td><td>Size of the entire distributed
939 | +html+ array, with<br>attributes <tt>xwidth</tt>, <tt>ywidth</tt>,
940 | +html+ <tt>zwidth</tt>, <tt>xsize</tt>**, <tt>ysize</tt>**,
941 | +html+ <tt>zsize</tt>** and <tt>fields</tt></td></tr>
942 | +*/
943 |
944 | thenode = xmlNewChild (thedoc->children, NULL, BAD_CAST "GlobalSize", NULL);
945 | snprintf ((char *) buffer, 999, "%d", nx);
946 | xmlSetProp (thenode, BAD_CAST "xwidth", buffer);
947 | snprintf ((char *) buffer, 999, "%g", wx);
948 | xmlSetProp (thenode, BAD_CAST "xsize", buffer);
949 | if (dim > 1)
950 | {
951 | snprintf ((char *) buffer, 999, "%d", ny);
952 | xmlSetProp (thenode, BAD_CAST "ywidth", buffer);
953 | snprintf ((char *) buffer, 999, "%g", wy);
954 | xmlSetProp (thenode, BAD_CAST "ysize", buffer);
955 | if (dim > 2)
956 | {
957 | snprintf ((char *) buffer, 999, "%d", nz);
958 | xmlSetProp (thenode, BAD_CAST "zwidth", buffer);
959 | snprintf ((char *) buffer, 999, "%g", wz);
960 | xmlSetProp (thenode, BAD_CAST "zsize", buffer);
961 | }
962 | }
963 | snprintf ((char *) buffer, 999, "%d", dof);
964 | xmlSetProp (thenode, BAD_CAST "fields", buffer);
965 |
966 | /*+
967 | +latex+{\tt LocalSize} & Size of the entire local part of the array,
968 | +latex+\\ & with attributes {\tt xwidth}, {\tt ywidth} and {\tt zwidth}
969 | +latex+\\ \hline
970 | +html+ <tr><td><tt>LocalSize</tt></td><td>Size of the local part of the
971 | +html+ array, with<br>attributes <tt>xwidth</tt>, <tt>ywidth</tt> and
972 | +html+ <tt>zwidth</tt></td></tr>
973 | +*/
974 |
975 | thenode = xmlNewChild (thedoc->children, NULL, BAD_CAST "LocalSize", NULL);
976 | snprintf ((char *) buffer, 999, "%d", xm);
977 | xmlSetProp (thenode, BAD_CAST "xwidth", buffer);
978 | if (dim > 1)
979 | {
980 | snprintf ((char *) buffer, 999, "%d", ym);
981 | xmlSetProp (thenode, BAD_CAST "ywidth", buffer);
982 | if (dim > 2)
983 | {
984 | snprintf ((char *) buffer, 999, "%d", zm);
985 | xmlSetProp (thenode, BAD_CAST "zwidth", buffer);
986 | }
987 | }
988 |
989 | /*+
990 | +latex+{\tt Stencil} & Stencil and periodic data, with attributes
991 | +latex+\\ & {\tt width}, {\tt type} and {\tt periodic} (using {\tt PETSc
992 | +latex+enum} values) \\ \hline
993 | +html+ <tr><td><tt>Stencil</tt></td><td>Stencil and periodic data, with
994 | +html+ attributes<br><tt>width</tt>, <tt>type</tt> and <tt>periodic</tt>
995 | +html+ (using <tt>PETSc enum</tt> values)</td></tr>
996 | +*/
997 |
998 | thenode = xmlNewChild (thedoc->children, NULL, BAD_CAST "Stencil", NULL);
999 | snprintf ((char *) buffer, 999, "%d", sw);
1000 | xmlSetProp (thenode, BAD_CAST "width", buffer);
1001 | snprintf ((char *) buffer, 999, "%d", (int) st);
1002 | xmlSetProp (thenode, BAD_CAST "type", buffer);
1003 | snprintf ((char *) buffer, 999, "%d", (int) wrap);
1004 | xmlSetProp (thenode, BAD_CAST "periodic", buffer);
1005 |
1006 | /*+
1007 | +latex+{\tt Field} & Data on each field, attributes {\tt name},
1008 | +latex+{\tt type}**, {\tt min} and {\tt max} \\ \hline
1009 | +html+ <tr><td><tt>Field</tt></td><td>Data on each field, attributes
1010 | +html+ <tt>name</tt>, <tt>type</tt>*, <tt>min</tt> and
1011 | +html+ <tt>max</tt></td></tr>
1012 | +*/
1013 |
1014 | for (i=0; i<dof; i++)
1015 | {
1016 | thenode = xmlNewChild (thedoc->children, NULL, BAD_CAST "Field", NULL);
1017 | xmlSetProp (thenode, BAD_CAST "name", BAD_CAST fieldnames [i]);
1018 | snprintf ((char *) buffer, 999, "0x%.2x", fieldtypes [i]);
1019 | xmlSetProp (thenode, BAD_CAST "type", buffer);
1020 | snprintf ((char *) buffer, 999, "%g", fieldmin [i]);
1021 | xmlSetProp (thenode, BAD_CAST "min", buffer);
1022 | snprintf ((char *) buffer, 999, "%g", fieldmax [i]);
1023 | xmlSetProp (thenode, BAD_CAST "max", buffer);
1024 | }
1025 |
1026 | /*+
1027 | +latex+{\tt User} & User parameters, attributes {\tt name} and {\tt value}
1028 | +latex+\\ \hline \end{tabular} \end{center}
1029 | +html+ <tr><td><tt>User</tt></td><td>User parameters, attributes
1030 | +html+ <tt>name</tt> and <tt>value</tt></td></tr></table></center>
1031 | *Lossy compression to smaller data types.
1032 | +latex+\par
1033 | +html+ <br>
1034 | **Represents new attribute for IlluMulti 0.2 file format.
1035 | +*/
1036 |
1037 | for (i=0; i<usermetacount; i++)
1038 | {
1039 | thenode = xmlNewChild (thedoc->children, NULL, BAD_CAST "User", NULL);
1040 | xmlSetProp (thenode, BAD_CAST "name", BAD_CAST usermetanames [i]);
1041 | xmlSetProp (thenode, BAD_CAST "value", BAD_CAST usermetadata [i]);
1042 | }
1043 |
1044 | if (snprintf ((char *)filename, 999, "%s.cpu%.4d.meta", basename, rank) >999)
1045 | SETERRQ1 (PETSC_ERR_ARG_SIZ, "Base file name too long (%d chars, 960 max)",
1046 | strlen (basename));
1047 | if (!(outfile = fopen (filename, "w")))
1048 | SETERRQ4 (PETSC_ERR_FILE_OPEN,
1049 | "CPU %d: error %d (%s) opening output file \"%s\"",
1050 | rank, errno, strerror (errno), filename);
1051 | xmlDocDump (outfile, thedoc);
1052 | xmlFreeDoc (thedoc);
1053 | fclose(outfile);
1054 |
1055 | return 0;
1056 | }
1057 |
1058 |
1059 | #undef __FUNCT__
1060 | #define __FUNCT__ "IlluMultiStoreData"
1061 |
1062 | /*++++++++++++++++++++++++++++++++++++++
1063 | This function stores the data file.
1064 |
1065 | int IlluMultiStoreData It returns zero or an error code.
1066 |
1067 | MPI_Comm comm MPI communicator, if NULL it uses PETSC_COMM_WORLD.
1068 |
1069 | PetscScalar *globalarray Array from which to save the (local) data.
1070 |
1071 | char *basename Base file name.
1072 |
1073 | int rank CPU number to read data for.
1074 |
1075 | int compressed Data compression: if zero then no compression (fastest), 1-9
1076 | then gzip compression level, 10-15 then gzip --best. If 16-31 then save
1077 | guint32s representing relative values between min and max for each field,
1078 | compressed according to this value minus 16. Likewise for 32-47 and
1079 | guint16s, and 48-63 for guint8s. Yes, these alternative formats lose
1080 | information and can't be used for accurate checkpointing, but they should
1081 | retain enough data for visualization (except perhaps for the guint8s, which
1082 | are possibly acceptable for vectors but likely not contours).
1083 |
1084 | int gridpoints Number of gridpoints to store data for.
1085 |
1086 | int dof Degrees of freedom at each node, a.k.a. number of field variables.
1087 |
1088 | PetscScalar *fieldmin Minimum value of each field variable.
1089 |
1090 | PetscScalar *fieldmax Maximum value of each field variable.
1091 | ++++++++++++++++++++++++++++++++++++++*/
1092 |
1093 | static int IlluMultiStoreData
1094 | (MPI_Comm comm, PetscScalar *globalarray, char *basename, int rank,
1095 | int compressed, int gridpoints, int dof, PetscScalar *fieldmin,
1096 | PetscScalar *fieldmax)
1097 | {
1098 | int i;
1099 | char filename [1000];
1100 | FILE *outfile;
1101 | size_t writeout;
1102 |
1103 | if (compressed & COMPRESS_GZIP_MASK)
1104 | {
1105 | if ((compressed & COMPRESS_GZIP_MASK) >= 1 &&
1106 | (compressed & COMPRESS_GZIP_MASK) <= 9)
1107 | {
1108 | if (snprintf (filename, 999, "gzip -c -%d > %s.cpu%.4d.data",
1109 | compressed & COMPRESS_GZIP_MASK, basename, rank) > 999)
1110 | SETERRQ1 (PETSC_ERR_ARG_SIZ,
1111 | "Base file name too long (%d chars, 960 max)",
1112 | strlen (basename));
1113 | }
1114 | else
1115 | {
1116 | if (snprintf (filename,999, "gzip -c --best > %s.cpu%.4d.data",
1117 | basename,rank) > 999)
1118 | SETERRQ1 (PETSC_ERR_ARG_SIZ,
1119 | "Base file name too long (%d chars, 960 max)",
1120 | strlen (basename));
1121 | }
1122 | if (!(outfile = popen (filename, "w")))
1123 | SETERRQ4 (PETSC_ERR_FILE_OPEN,
1124 | "CPU %d: error %d (%s) opening output pipe \"%s\"",
1125 | rank, errno, strerror (errno), filename);
1126 | }
1127 | else
1128 | {
1129 | if (snprintf (filename, 999, "%s.cpu%.4d.data", basename, rank) > 999)
1130 | SETERRQ1 (PETSC_ERR_ARG_SIZ,
1131 | "Base file name too long (%d chars, 960 max)",
1132 | strlen (basename));
1133 | if (!(outfile = fopen (filename, "w")))
1134 | SETERRQ4 (PETSC_ERR_FILE_OPEN,
1135 | "CPU %d: error %d (%s) opening output file \"%s\"",
1136 | rank, errno, strerror (errno), filename);
1137 | }
1138 |
1139 | switch (compressed & COMPRESS_INT_MASK)
1140 | {
1141 | case COMPRESS_INT_LONG:
1142 | {
1143 | guint32 *storage;
1144 | if (!(storage = (guint32 *) malloc
1145 | (gridpoints * dof * sizeof (guint32))))
1146 | SETERRQ (PETSC_ERR_MEM, "Can't allocate data compression buffer");
1147 | for (i=0; i<gridpoints*dof; i++)
1148 | storage [i] = (guint32)
1149 | (4294967295. * (globalarray [i] - fieldmin [i%dof]) /
1150 | (fieldmax [i%dof] - fieldmin [i%dof]));
1151 | writeout = fwrite (storage, sizeof (guint32), gridpoints*dof, outfile);
1152 | free (storage);
1153 | break;
1154 | }
1155 | case COMPRESS_INT_SHORT:
1156 | {
1157 | guint16 *storage;
1158 | if (!(storage = (guint16 *) malloc
1159 | (gridpoints * dof * sizeof (guint16))))
1160 | SETERRQ (PETSC_ERR_MEM, "Can't allocate data compression buffer");
1161 | for (i=0; i<gridpoints*dof; i++)
1162 | storage [i] = (guint16)
1163 | (65535. * (globalarray [i] - fieldmin [i%dof]) /
1164 | (fieldmax [i%dof] - fieldmin [i%dof]));
1165 | writeout = fwrite (storage, sizeof (guint16), gridpoints*dof,outfile);
1166 | free (storage);
1167 | break;
1168 | }
1169 | case COMPRESS_INT_CHAR:
1170 | {
1171 | guint8 *storage;
1172 | if (!(storage = (guint8 *) malloc
1173 | (gridpoints * dof * sizeof (guint8))))
1174 | SETERRQ (PETSC_ERR_MEM, "Can't allocate data compression buffer");
1175 | for (i=0; i<gridpoints*dof; i++)
1176 | storage [i] = (guint8)
1177 | (255. * (globalarray [i] - fieldmin [i%dof]) /
1178 | (fieldmax [i%dof] - fieldmin [i%dof]));
1179 | writeout = fwrite (storage, sizeof (guint8), gridpoints*dof, outfile);
1180 | free (storage);
1181 | break;
1182 | }
1183 | default:
1184 | writeout = fwrite (globalarray, sizeof (PetscScalar), gridpoints*dof,
1185 | outfile);
1186 | }
1187 |
1188 | if (compressed & COMPRESS_GZIP_MASK)
1189 | pclose (outfile);
1190 | else
1191 | fclose (outfile);
1192 |
1193 | if (writeout < gridpoints*dof)
1194 | {
1195 | int writerr = ferror (outfile);
1196 | SETERRQ5 (PETSC_ERR_FILE_READ,
1197 | "CPU %d: error %d (%s) during write, sent %d of %d items",
1198 | rank, writerr, strerror (writerr), writeout, gridpoints*dof);
1199 | }
1200 |
1201 | return 0;
1202 | }
1203 |
1204 |
1205 | #undef __FUNCT__
1206 | #define __FUNCT__ "checkagree"
1207 |
1208 | /*++++++++++++++++++++++++++++++++++++++
1209 | Ancillary routine for
1210 | +latex+{\tt IlluMultiRead()}:
1211 | +html+ <tt>IlluMultiRead()</tt>:
1212 | checks agreement of parameters and reports disagreement if necessary.
1213 |
1214 | int checkagree Returns 0 if they agree, 1 otherwise.
1215 |
1216 | MPI_Comm comm MPI communicator, if NULL it uses PETSC_COMM_WORLD.
1217 |
1218 | int da Integer parameter from the existing DA.
1219 |
1220 | int file Integer parameter read from the file.
1221 |
1222 | char *parameter Parameter name for reporting.
1223 | ++++++++++++++++++++++++++++++++++++++*/
1224 |
1225 | static inline int checkagree (MPI_Comm comm, int da, int file, char *parameter)
1226 | {
1227 | int ierr;
1228 |
1229 | if (da == file)
1230 | return 0;
1231 | SETERRQ3 (PETSC_ERR_ARG_OUTOFRANGE,
1232 | "Disagreement in parameter %s: da has %d, file %d\n",
1233 | parameter, da, file);
1234 | }
1235 |
1236 |
1237 | #undef __FUNCT__
1238 | #define __FUNCT__ "IlluMultiRead"
1239 |
1240 | /*++++++++++++++++++++++++++++++++++++++
1241 | This reads the data into an existing distributed array and vector, checking
1242 | that the sizes are right etc.
1243 |
1244 | int IlluMultiRead It returns zero or an error code.
1245 |
1246 | MPI_Comm comm MPI communicator, if NULL it uses PETSC_COMM_WORLD.
1247 |
1248 | DA theda Distributed array object controlling the data to read.
1249 |
1250 | Vec X Vector into which to read the data.
1251 |
1252 | char *basename Base file name.
1253 |
1254 | int *usermetacount Pointer to an int where we put the number of user metadata
1255 | parameters loaded.
1256 |
1257 | char ***usermetanames Pointer to a char ** where the loaded parameter names
1258 | are stored. This is
1259 | +latex+{\tt malloc}ed by this function, so a call to {\tt free()}
1260 | +html+ <tt>malloc</tt>ed by this function, so a call to <tt>free()</tt>
1261 | is needed to free up its data.
1262 |
1263 | char ***usermetadata Pointer to a char ** where the loaded parameter strings
1264 | are stored. This is
1265 | +latex+{\tt malloc}ed by this function, so a call to {\tt free()}
1266 | +html+ <tt>malloc</tt>ed by this function, so a call to <tt>free()</tt>
1267 | is needed to free up its data.
1268 | ++++++++++++++++++++++++++++++++++++++*/
1269 |
1270 | int IlluMultiRead
1271 | (MPI_Comm comm, DA theda, Vec X, char *basename,
1272 | int *usermetacount, char ***usermetanames, char ***usermetadata)
1273 | {
1274 | int dim, px,py,pz, nx,ny,nz, xm,ym,zm, dof, sw, size, rank, ierr;
1275 | DAPeriodicType wrap, fwrap;
1276 | DAStencilType st, fst;
1277 | int fdim, fpx,fpy,fpz, fnx,fny,fnz, fxm,fym,fzm, fdof, fsw, fsize = 1;
1278 | int compressed, wrongendian;
1279 | char filename[1000], **fieldnames;
1280 | FILE *infile;
1281 | PetscScalar *globalarray, *fieldmin, *fieldmax;
1282 |
1283 | if (!comm)
1284 | comm = PETSC_COMM_WORLD;
1285 |
1286 | /*+ First it gets the properties of the distributed array for comparison with
1287 | the metadata. +*/
1288 | ierr = DAGetInfo (theda, &dim, &nx,&ny,&nz, &px,&py,&pz, &dof, &sw, &wrap,
1289 | &st); CHKERRQ (ierr);
1290 | ierr = DAGetCorners (theda, PETSC_NULL,PETSC_NULL,PETSC_NULL, &xm,&ym,&zm);
1291 | CHKERRQ (ierr);
1292 |
1293 | /*+ Next it parses the XML metadata file into the document tree, and reads
1294 | its content into the appropriate structures, comparing parameters with
1295 | those of the existing distributed array structure. +*/
1296 | ierr = MPI_Comm_size (comm, &size); CHKERRQ (ierr);
1297 | ierr = MPI_Comm_rank (comm, &rank); CHKERRQ (ierr);
1298 |
1299 | DPRINTF ("Parsing XML metadata file from basename %s\n", basename);
1300 | if (ierr = IlluMultiParseXML
1301 | (comm, basename, rank, &compressed, &wrongendian, &fdim, &fpx,&fpy,&fpz,
1302 | &fnx,&fny,&fnz, PETSC_NULL,PETSC_NULL,PETSC_NULL, &fxm,&fym,&fzm, &fdof,
1303 | &fsw, &fst, &fwrap, &fieldnames, PETSC_NULL, &fieldmin, &fieldmax,
1304 | usermetacount, usermetanames, usermetadata))
1305 | return ierr;
1306 |
1307 | /* The size>1 checks are because we support loading multiple CPUs' data into
1308 | a 1-CPU distributed array, as long as the global dimensions are right. */
1309 | DPRINTF ("Checking size agreement between DA and data to read\n",0);
1310 | fsize = fpx * ((dim > 1) ? fpy : 1) * ((dim > 2) ? fpz : 1);
1311 | if (ierr = checkagree (comm, dim, fdim, "dimensions")) return ierr;
1312 | if (size > 1 || fsize == 1) {
1313 | if (ierr = checkagree (comm, px, fpx, "cpu xwidth")) return ierr;
1314 | if (dim > 1) {
1315 | if (ierr = checkagree (comm, py, fpy, "cpu ywidth")) return ierr;
1316 | if (dim > 2) {
1317 | if (ierr = checkagree (comm, pz, fpz, "cpu zwidth")) return ierr; }}}
1318 | if (ierr = checkagree (comm, nx, fnx, "global xwidth")) return ierr;
1319 | if (dim > 1) {
1320 | if (ierr = checkagree (comm, ny, fny, "global ywidth")) return ierr;
1321 | if (dim > 2) {
1322 | if (ierr = checkagree (comm, nz, fnz, "global zwidth")) return ierr; }}
1323 | if (size > 1 || fsize == 1) {
1324 | if (ierr = checkagree (comm, xm, fxm, "local xwidth")) return ierr;
1325 | if (dim > 1) {
1326 | if (ierr = checkagree (comm, ym, fym, "local ywidth")) return ierr;
1327 | if (dim > 2) {
1328 | if (ierr = checkagree (comm, zm, fzm, "local zwidth")) return ierr; }}}
1329 | if (ierr = checkagree (comm, dof, fdof, "degrees of freedom")) return ierr;
1330 | if (ierr = checkagree (comm, sw, fsw, "stencil width")) return ierr;
1331 | if (ierr = checkagree (comm, (int)st, (int)fst, "stencil type")) return ierr;
1332 | if (ierr = checkagree (comm, (int)wrap, (int)fwrap, "periodic type"))
1333 | return ierr;
1334 |
1335 | /*+ Then it streams in the data in one big slurp. +*/
1336 | ierr = VecGetArray (X, &globalarray); CHKERRQ (ierr);
1337 | if (size > 1 || fsize == 1) /* Default condition */
1338 | {
1339 | DPRINTF ("Reading data\n",0);
1340 | ierr = IlluMultiParseData
1341 | (comm, globalarray, basename, rank, compressed,
1342 | xm * ((dim>1)?ym:1) * ((dim>2)?zm:1), dof, wrongendian, fieldmin,
1343 | fieldmax);
1344 | if (ierr)
1345 | {
1346 | xm = VecRestoreArray (X, &globalarray); CHKERRQ (xm);
1347 | return ierr;
1348 | }
1349 | }
1350 | else /* Getting distributed data into a single array */
1351 | {
1352 | int i,j,k, xs,ys,zs;
1353 | PetscScalar *storage;
1354 |
1355 | /* Incrementally read in data to storage, store it in globalarray */
1356 | /* First loop through the stored CPUs */
1357 | for (pz=0, zs=0; pz<((dim>2)?fpz:1); pz++, zs+=fzm)
1358 | for (py=0, ys=0; py<((dim>1)?fpy:1); py++, ys+=fym)
1359 | for (px=0, xs=0; px<fpx; px++, xs+=fxm, rank++)
1360 | {
1361 | int gridpoints;
1362 |
1363 | DPRINTF ("Parsing XML metadata file for CPU %d\n", rank);
1364 | if (ierr = IlluMultiParseXML
1365 | (comm, basename, rank, &compressed, &wrongendian, &fdim,
1366 | &fpx,&fpy,&fpz, &fnx,&fny,&fnz,
1367 | PETSC_NULL,PETSC_NULL,PETSC_NULL,
1368 | &fxm,&fym,&fzm, &fdof, &fsw, &fst, &fwrap,
1369 | &fieldnames, PETSC_NULL, &fieldmin, &fieldmax,
1370 | usermetacount, usermetanames, usermetadata))
1371 | return ierr;
1372 | gridpoints = fxm * ((dim>1)?fym:1) * ((dim>2)?fzm:1);
1373 |
1374 | /* This allocate/read/copy/free storage business is annoying.
1375 | Replacing it with "zero-copy" would involve changing the
1376 | IlluMultiParseData() API to allow loading into part of the
1377 | local array. But I'll leave that for future expansion, when
1378 | this will be able to do arbitrary m->n CPU loading as long as
1379 | the global sizes match. */
1380 | storage = (PetscScalar *) malloc
1381 | (gridpoints * dof * sizeof (PetscScalar));
1382 | DPRINTF ("Reading data for CPU %d\n", rank);
1383 | ierr = IlluMultiParseData
1384 | (comm, storage, basename, rank, compressed, gridpoints, dof,
1385 | wrongendian, fieldmin, fieldmax);
1386 |
1387 | fxm *= dof; /* so fxm is the number of PetscScalars to copy */
1388 | i=1;
1389 | /* Now distribute */
1390 | for (k=0; k<((dim>2)?fzm:1); k++)
1391 | for (j=0; j<((dim>1)?fym:1); j++)
1392 | BLAScopy_ (&fxm, storage + (k*fym + j)*fxm /* *dof */, &i,
1393 | globalarray + (((zs+k)*ny + ys+j)*nx + xs)*dof, &i);
1394 |
1395 | free (storage);
1396 | fxm /= dof;
1397 | }
1398 | }
1399 |
1400 | ierr = VecRestoreArray (X, &globalarray); CHKERRQ (ierr);
1401 |
1402 | return 0;
1403 | }
1404 |
1405 |
1406 | #undef __FUNCT__
1407 | #define __FUNCT__ "IlluMultiLoad"
1408 |
1409 | /*++++++++++++++++++++++++++++++++++++++
1410 | This creates a new distributed array of the appropriate size and loads the
1411 | data into the vector contained in it (as retrieved by
1412 | +latex+{\tt DAGetVector()}).
1413 | +html+ <tt>DAGetVector()</tt>).
1414 | It also reads the user metadata parameters into arrays stored at the supplied
1415 | pointers.
1416 |
1417 | int IlluMultiLoad It returns zero or an error code.
1418 |
1419 | MPI_Comm comm MPI communicator, if NULL it uses PETSC_COMM_WORLD.
1420 |
1421 | char *basename Base file name.
1422 |
1423 | DA *theda Pointer to a DA object (to be created by this function).
1424 |
1425 | PetscScalar *wx Physical overall width in the
1426 | +latex+$x$-direction.
1427 | +html+ <i>x</i>-direction.
1428 |
1429 | PetscScalar *wy Physical overall width in the
1430 | +latex+$y$-direction.
1431 | +html+ <i>y</i>-direction.
1432 |
1433 | PetscScalar *wz Physical overall width in the
1434 | +latex+$z$-direction.
1435 | +html+ <i>z</i>-direction.
1436 |
1437 | field_plot_type **fieldtypes Data (plot) types for field variables.
1438 |
1439 | int *usermetacount Pointer to an int where we put the number of user metadata
1440 | parameters loaded.
1441 |
1442 | char ***usermetanames Pointer to a char ** where the loaded parameter names
1443 | are stored. This is
1444 | +latex+{\tt malloc}ed by this function, so a call to {\tt free()}
1445 | +html+ <tt>malloc</tt>ed by this function, so a call to <tt>free()</tt>
1446 | is needed to free up its data.
1447 |
1448 | char ***usermetadata Pointer to a char ** where the loaded parameter strings
1449 | are stored. This is
1450 | +latex+{\tt malloc}ed by this function, so a call to {\tt free()}
1451 | +html+ <tt>malloc</tt>ed by this function, so a call to <tt>free()</tt>
1452 | is needed to free up its data.
1453 | ++++++++++++++++++++++++++++++++++++++*/
1454 |
1455 | int IlluMultiLoad
1456 | (MPI_Comm comm, char *basename, DA *theda, PetscScalar *wx,PetscScalar *wy,
1457 | PetscScalar *wz, field_plot_type **fieldtypes, int *usermetacount,
1458 | char ***usermetanames, char ***usermetadata)
1459 | {
1460 | int dim, px,py,pz, nx,ny,nz, dof, sw, fxm,fym,fzm, rank,size, ierr;
1461 | int compressed, wrongendian, fpx,fpy,fpz, fsize, xm,ym,zm;
1462 | DAPeriodicType wrap;
1463 | DAStencilType st;
1464 | char filename[1000], **fieldnames;
1465 | xmlChar *buffer;
1466 | FILE *infile;
1467 | xmlDocPtr thedoc;
1468 | xmlNodePtr thenode;
1469 | PetscScalar *globalarray, *fieldmin, *fieldmax;
1470 | Vec X;
1471 |
1472 | if (!comm)
1473 | comm = PETSC_COMM_WORLD;
1474 |
1475 | ierr = MPI_Comm_rank (comm, &rank); CHKERRQ (ierr);
1476 | ierr = MPI_Comm_size (comm, &size); CHKERRQ (ierr);
1477 |
1478 | /*+ First it gets the parameters from the XML file. +*/
1479 | DPRINTF ("Parsing XML metadata file from basename %s\n", basename);
1480 | if (ierr = IlluMultiParseXML
1481 | (comm, basename, rank, &compressed, &wrongendian, &dim, &fpx,&fpy,&fpz,
1482 | &nx,&ny,&nz, wx,wy,wz, &fxm,&fym,&fzm, &dof, &sw, &st, &wrap,
1483 | &fieldnames, fieldtypes, &fieldmin, &fieldmax,
1484 | usermetacount, usermetanames, usermetadata))
1485 | return ierr;
1486 |
1487 | fsize = fpx * ((dim > 1) ? fpy : 1) * ((dim > 2) ? fpz : 1);
1488 | if (size == fsize) { /* Default condition: n->n */
1489 | px = fpx; py = fpy; pz = fpz; }
1490 | else if (size == 1) /* Special condition: n->1 */
1491 | px = py = pz = 1;
1492 | else /* Error: incorrect number of CPUs */
1493 | SETERRQ5 (PETSC_ERR_ARG_SIZ,
1494 | "Incorrect CPU count: current %d, stored %dx%dx%d=%d\n",
1495 | size, fpx,fpy,fpz, fsize);
1496 |
1497 | /*+ Next it creates a distributed array based on those parameters, and sets
1498 | the names of its fields. +*/
1499 | switch (dim)
1500 | {
1501 | case 1:
1502 | {
1503 | DPRINTF ("Creating %d-%d 1-D DA\n", nx, dof);
1504 | ierr = DACreate1d (comm, wrap, nx, dof, sw, PETSC_NULL,
1505 | theda); CHKERRQ (ierr);
1506 | break;
1507 | }
1508 | case 2:
1509 | {
1510 | DPRINTF ("Creating %dx%d-%d 2-D DA\n", nx,ny, dof);
1511 | ierr = DACreate2d (comm, wrap, st, nx,ny, px,py, dof, sw,
1512 | PETSC_NULL, PETSC_NULL, theda); CHKERRQ (ierr);
1513 | break;
1514 | }
1515 | case 3:
1516 | {
1517 | DPRINTF ("Creating %dx%dx%d-%d 3-D DA\n", nx,ny,nz, dof);
1518 | ierr = DACreate3d (comm, wrap, st, nx,ny,nz, px,py,pz, dof,
1519 | sw, PETSC_NULL, PETSC_NULL, PETSC_NULL, theda);
1520 | CHKERRQ (ierr);
1521 | break;
1522 | }
1523 | default:
1524 | SETERRQ1 (PETSC_ERR_ARG_OUTOFRANGE,
1525 | "Unsupported number of dimensions %d\n", dim);
1526 | }
1527 |
1528 | ierr = DAGetCorners (*theda, PETSC_NULL, PETSC_NULL, PETSC_NULL,
1529 | &xm,&ym,&zm); CHKERRQ (ierr);
1530 | if(size==fsize && (xm != fxm || ym != fym || zm != fzm))
1531 | SETERRQ6 (PETSC_ERR_ARG_SIZ,
1532 | "Local array size mismatch: DA %dx%dx%d, stored %dx%dx%d\n",
1533 | xm,ym,zm, fxm,fym,fzm);
1534 |
1535 | for (px=0; px<dof; px++)
1536 | DASetFieldName (*theda, px, fieldnames [px]);
1537 |
1538 | /*+ Then it streams the data into the distributed array's vector in one big
1539 | slurp. +*/
1540 | DPRINTF ("Getting global vector and array\n",0);
1541 | ierr = DAGetGlobalVector (*theda, &X); CHKERRQ (ierr);
1542 | ierr = VecGetArray (X, &globalarray); CHKERRQ (ierr);
1543 | if (size > 1 || fsize == 1) /* Default condition */
1544 | {
1545 | DPRINTF ("Loading data in parallel\n",0);
1546 | ierr = IlluMultiParseData
1547 | (comm, globalarray, basename, rank, compressed,
1548 | fxm * ((dim>1)?fym:1) * ((dim>2)?fzm:1), dof, wrongendian, fieldmin,
1549 | fieldmax);
1550 | if (ierr)
1551 | {
1552 | fxm = VecRestoreArray (X, &globalarray); CHKERRQ (fxm);
1553 | return ierr;
1554 | }
1555 | }
1556 | else /* Getting distributed data into a single array */
1557 | {
1558 | int i,j,k, xs,ys,zs;
1559 | PetscScalar *storage;
1560 |
1561 | /* Incrementally read in data to storage, store it in globalarray */
1562 | /* First loop through the stored CPUs */
1563 | DPRINTF ("Loading data into a single array on 1 CPU\n",0);
1564 | for (pz=0, zs=0; pz<((dim>2)?fpz:1); pz++, zs+=fzm)
1565 | for (py=0, ys=0; py<((dim>1)?fpy:1); py++, ys+=fym)
1566 | for (px=0, xs=0; px<fpx; px++, xs+=fxm, rank++)
1567 | {
1568 | int gridpoints, fdim, fnx,fny,fnz, fdof, fsw;
1569 | DAPeriodicType fwrap;
1570 | DAStencilType fst;
1571 |
1572 | if (ierr = IlluMultiParseXML
1573 | (comm, basename, rank, &compressed, &wrongendian, &fdim,
1574 | &fpx,&fpy,&fpz, &fnx,&fny,&fnz,
1575 | PETSC_NULL,PETSC_NULL,PETSC_NULL,
1576 | &fxm,&fym,&fzm, &fdof, &fsw, &fst, &fwrap,
1577 | &fieldnames, PETSC_NULL, &fieldmin, &fieldmax,
1578 | usermetacount, usermetanames, usermetadata))
1579 | return ierr;
1580 | gridpoints = fxm * ((dim>1)?fym:1) * ((dim>2)?fzm:1);
1581 |
1582 | /* This allocate/read/copy/free storage business is annoying.
1583 | Replacing it with "zero-copy" would involve changing the
1584 | IlluMultiParseData() API to allow loading into part of the
1585 | local array. But I'll leave that for future expansion, when
1586 | this will be able to do arbitrary m->n CPU loading as long as
1587 | the global sizes match. */
1588 | storage = (PetscScalar *) malloc
1589 | (gridpoints * dof * sizeof (PetscScalar));
1590 | ierr = IlluMultiParseData
1591 | (comm, storage, basename, rank, compressed, gridpoints, dof,
1592 | wrongendian, fieldmin, fieldmax);
1593 |
1594 | fxm *= dof; /* so fxm is the number of PetscScalars to copy */
1595 | i=1;
1596 | /* Now distribute */
1597 | for (k=0; k<((dim>2)?fzm:1); k++)
1598 | for (j=0; j<((dim>1)?fym:1); j++)
1599 | BLAScopy_ (&fxm, storage + (k*fym + j)*fxm /* *dof */, &i,
1600 | globalarray + (((zs+k)*ny + ys+j)*nx + xs)*dof, &i);
1601 |
1602 | free (storage);
1603 | fxm /= dof;
1604 | }
1605 | }
1606 |
1607 | DPRINTF ("Done.\n",0);
1608 | ierr = VecRestoreArray (X, &globalarray); CHKERRQ (ierr);
1609 | ierr = DARestoreGlobalVector (*theda, &X); CHKERRQ (ierr);
1610 |
1611 | return 0;
1612 | }
1613 |
1614 |
1615 | #undef __FUNCT__
1616 | #define __FUNCT__ "IlluMultiSave"
1617 |
1618 | /*++++++++++++++++++++++++++++++++++++++
1619 | This saves the vector X in multiple files, two per process.
1620 |
1621 | int IlluMultiSave it returns zero or an error code.
1622 |
1623 | MPI_Comm comm MPI communicator, if NULL it uses PETSC_COMM_WORLD.
1624 |
1625 | DA theda Distributed array object controlling data saved.
1626 |
1627 | Vec X Vector whose data are actually being saved.
1628 |
1629 | char *basename Base file name.
1630 |
1631 | int usermetacount Number of user metadata parameters.
1632 |
1633 | char **usermetanames User metadata parameter names.
1634 |
1635 | char **usermetadata User metadata parameter strings.
1636 |
1637 | int compressed Data compression: if zero then no compression (fastest), 1-9
1638 | then gzip compression level, 10-15 then gzip --best. If 16-31 then save
1639 | guint32s representing relative values between min and max for each field,
1640 | compressed according to this value minus 16. Likewise for 32-47 and
1641 | guint16s, and 48-63 for guint8s. Yes, these alternative formats lose
1642 | information and can't be used for accurate checkpointing, but they should
1643 | retain enough data for visualization (except perhaps for the guint8s, which
1644 | are possibly acceptable for vectors but certainly not contours).
1645 | ++++++++++++++++++++++++++++++++++++++*/
1646 |
1647 | int IlluMultiSave
1648 | (MPI_Comm comm, DA theda, Vec X, char *basename, PetscScalar wx,PetscScalar wy,
1649 | PetscScalar wz, field_plot_type *fieldtypes, int usermetacount,
1650 | char **usermetanames, char **usermetadata, int compressed)
1651 | {
1652 | int dim, nx,ny,nz, px,py,pz, dof, sw, xs,ys,zs, xm,ym,zm, rank, i, ierr;
1653 | DAPeriodicType wrap;
1654 | DAStencilType st;
1655 | FILE *outfile;
1656 | char filename[1000], **fieldnames;
1657 | PetscReal *fieldmin, *fieldmax;
1658 | PetscScalar *globalarray;
1659 | guint64 *array64;
1660 | guint32 *array32;
1661 | guint16 *array16;
1662 | guint8 *array8;
1663 |
1664 | if (!comm)
1665 | comm = PETSC_COMM_WORLD;
1666 |
1667 | /*+First a check to verify a supported value of
1668 | +latex+{\tt compressed},
1669 | +html+ <tt>compressed</tt>,
1670 | but no fancy guint* compression for complex!
1671 | +*/
1672 | #ifdef PETSC_USE_COMPLEX
1673 | if (compressed < 0 || compressed > COMPRESS_GZIP_MASK)
1674 | SETERRQ1(PETSC_ERR_ARG_OUTOFRANGE,"Unsupported compression type code 0x%x",
1675 | compressed);
1676 | #else
1677 | if (compressed < 0 || compressed > (COMPRESS_GZIP_MASK | COMPRESS_INT_MASK))
1678 | SETERRQ1(PETSC_ERR_ARG_OUTOFRANGE,"Unsupported compression type code 0x%x",
1679 | compressed);
1680 | #endif
1681 |
1682 | /*+Then get the distributed array parameters and processor number, and store
1683 | all this data in the XML .meta file. +*/
1684 | ierr = DAGetInfo (theda, &dim, &nx,&ny,&nz, &px,&py,&pz, &dof, &sw, &wrap,
1685 | &st); CHKERRQ (ierr);
1686 | ierr = DAGetCorners (theda, &xs,&ys,&zs, &xm,&ym,&zm); CHKERRQ (ierr);
1687 | ierr = MPI_Comm_rank (comm, &rank); CHKERRQ (ierr);
1688 |
1689 | if (!(fieldnames = (char **) malloc (dof * sizeof (char *))))
1690 | SETERRQ (PETSC_ERR_MEM, "Can't allocate field names");
1691 | if (!(fieldmin = (PetscScalar *) malloc (2 * dof * sizeof (PetscScalar))))
1692 | SETERRQ (PETSC_ERR_MEM, "Can't allocate field params");
1693 | fieldmax = fieldmin + dof;
1694 | for (i=0; i<dof; i++)
1695 | {
1696 | ierr = DAGetFieldName (theda, i, fieldnames + i); CHKERRQ (ierr);
1697 | ierr = VecStrideMin (X, i, PETSC_NULL, fieldmin + i); CHKERRQ (ierr);
1698 | ierr = VecStrideMax (X, i, PETSC_NULL, fieldmax + i); CHKERRQ (ierr);
1699 | }
1700 |
1701 | if (ierr = IlluMultiStoreXML
1702 | (comm, basename, rank, compressed, dim, px,py,pz, nx,ny,nz, wx,wy,wz,
1703 | xm,ym,zm, dof, sw, (int) st, (int) wrap, fieldnames,fieldtypes,
1704 | fieldmin,fieldmax, usermetacount, usermetanames, usermetadata))
1705 | return ierr;
1706 |
1707 | /*+ Finally, the data just stream out to the data file or gzip pipe in one
1708 | big lump. +*/
1709 | ierr = VecGetArray (X, &globalarray); CHKERRQ (ierr);
1710 | IlluMultiStoreData (comm, globalarray, basename, rank, compressed,
1711 | xm * ((dim>1)?ym:1) * ((dim>2)?zm:1), dof,
1712 | fieldmin, fieldmax);
1713 | ierr = VecRestoreArray (X, &globalarray); CHKERRQ (ierr);
1714 |
1715 | free (fieldmin);
1716 | free (fieldnames);
1717 |
1718 | return 0;
1719 | }