source: CIVL/examples/omp/HydroC/vtkfile.c@ 1aaefd4

main test-branch
Last change on this file since 1aaefd4 was ea777aa, checked in by Alex Wilton <awilton@…>, 3 years ago

Moved examples, include, build_default.properties, common.xml, and README out from dev.civl.com into the root of the repo.

git-svn-id: ​svn://vsl.cis.udel.edu/civl/trunk@5704 fb995dde-84ed-4084-dfe6-e5aef3e2452c
  • Property mode set to 100644
File size: 11.6 KB
Line 
1/*
2 A simple 2D hydro code
3 (C) Romain Teyssier : CEA/IRFU -- original F90 code
4 (C) Pierre-Francois Lavallee : IDRIS -- original F90 code
5 (C) Guillaume Colin de Verdiere : CEA/DAM -- for the C version
6 (C) Adèle Villiermet : CINES -- for FTI integration
7*/
8/*
9
10/*
11
12This software is governed by the CeCILL license under French law and
13abiding by the rules of distribution of free software. You can use,
14modify and/ or redistribute the software under the terms of the CeCILL
15license as circulated by CEA, CNRS and INRIA at the following URL
16"http://www.cecill.info".
17
18As a counterpart to the access to the source code and rights to copy,
19modify and redistribute granted by the license, users are provided only
20with a limited warranty and the software's author, the holder of the
21economic rights, and the successive licensors have only limited
22liability.
23
24In this respect, the user's attention is drawn to the risks associated
25with loading, using, modifying and/or developing or reproducing the
26software by the user in light of its specific status of free software,
27that may mean that it is complicated to manipulate, and that also
28therefore means that it is reserved for developers and experienced
29professionals having in-depth computer knowledge. Users are therefore
30encouraged to load and test the software's suitability as regards their
31requirements in conditions enabling the security of their systems and/or
32data to be ensured and, more generally, to use and operate it in the
33same conditions as regards security.
34
35The fact that you are presently reading this means that you have had
36knowledge of the CeCILL license and that you accept its terms.
37
38*/
39
40#ifdef MPI
41#include <mpi.h>
42#if FTI>0
43#include <fti.h>
44#endif
45#endif
46#include <stdio.h>
47#include <stdlib.h>
48#include <unistd.h>
49#include <string.h>
50#include <strings.h>
51#include <assert.h>
52#include <math.h>
53#include <sys/stat.h>
54#include <sys/types.h>
55
56
57#include "parametres.h"
58#include "utils.h"
59#include "vtkfile.h"
60#include "SplitSurface.h"
61
62typedef unsigned char byte;
63typedef int bool;
64
65const int false = 0;
66const int true = 1;
67
68const char s_CharPlusSign = '+';
69const char s_CharSlash = '/';
70
71static char SixBitToChar(byte b);
72static char *ToBase64(unsigned char *data, int length);
73
74char
75SixBitToChar(byte b) {
76 char c;
77 if (b < 26) {
78 c = (char) ((int) b + (int) 'A');
79 } else if (b < 52) {
80 c = (char) ((int) b - 26 + (int) 'a');
81 } else if (b < 62) {
82 c = (char) ((int) b - 52 + (int) '0');
83 } else if (b == 62) {
84 c = s_CharPlusSign;
85 } else {
86 c = s_CharSlash;
87 }
88 return c;
89}
90
91char *
92ToBase64(unsigned char *data, int length) {
93 int padding = length % 3;
94 int blocks = (length - 1) / 3 + 1;
95 size_t lalloc;
96 char *s;
97 int i;
98
99 if (length == 0)
100 return NULL;
101
102 if (padding > 0)
103 padding = 3 - padding;
104
105 // lalloc = (blocks * 4 + 1 + 16);
106 lalloc = blocks;
107 lalloc *= 4;
108 lalloc += 17;
109
110 s = malloc(lalloc);
111 if (s == NULL) {
112 fprintf(stderr, "Length=%d, blocks=%d lalloc=%ld\n", length, blocks, lalloc);
113 }
114 assert(s != NULL);
115
116 for (i = 0; i < blocks; i++) {
117 bool finalBlock = i == blocks - 1;
118 bool pad2 = false;
119 bool pad1 = false;
120 if (finalBlock) {
121 pad2 = padding == 2;
122 pad1 = padding > 0;
123 }
124
125 int index = i * 3;
126 byte b1 = data[index];
127 byte b2 = pad2 ? (byte) 0 : data[index + 1];
128 byte b3 = pad1 ? (byte) 0 : data[index + 2];
129
130 byte temp1 = (byte) ((b1 & 0xFC) >> 2);
131
132 byte temp = (byte) ((b1 & 0x03) << 4);
133 byte temp2 = (byte) ((b2 & 0xF0) >> 4);
134 temp2 += temp;
135
136 temp = (byte) ((b2 & 0x0F) << 2);
137 byte temp3 = (byte) ((b3 & 0xC0) >> 6);
138 temp3 += temp;
139
140 byte temp4 = (byte) (b3 & 0x3F);
141
142 index = i * 4;
143 s[index] = SixBitToChar(temp1);
144 s[index + 1] = SixBitToChar(temp2);
145 s[index + 2] = pad2 ? '=' : SixBitToChar(temp3);
146 s[index + 3] = pad1 ? '=' : SixBitToChar(temp4);
147 }
148 s[blocks * 4] = (byte) 0;
149 return s;
150}
151
152#define BINARY 1
153
154void
155vtkwpvd(int nout, char *r) {
156 char n[1024];
157 char vfname[1024];
158 int i;
159 FILE *vf = NULL;
160 char tmp[10];
161
162 vf = fopen("Hydro.pvd", "w");
163 fprintf(vf, "<?xml version=\"1.0\"?>\n");
164 fprintf(vf, " <VTKFile type=\"Collection\" version=\"0.1\" byte_order=\"LittleEndian\">\n");
165 fprintf(vf, " <Collection>\n");
166
167 for (i = 1; i <= nout; i++) {
168 sprintf(tmp, "%06d", i);
169 sprintf(n, "Dep/%c%c%c%c", tmp[0], tmp[1], tmp[2], tmp[3]);
170 sprintf(n, "%s/%c%c", n, tmp[4], tmp[5]);
171 sprintf(vfname, "%s/Hydro_%04d.pvtr", n, i);
172 fprintf(vf, " <DataSet timestep=\"%d\" part=\"0\" file=\"%s\" name=\"Asmb:FRAME\"/>\n", i, vfname);
173 }
174
175 fprintf(vf, " </Collection>\n");
176 fprintf(vf, "</VTKFile>\n");
177 fclose(vf);
178}
179
180void
181vtknm(char *n, int me, int nout) {
182 char tmp[10];
183
184 sprintf(tmp, "%06d", nout);
185 sprintf(n, "Dep");
186 if (me == 0) {
187 mkdir(n, 0777);
188 }
189 sprintf(n, "%s/%c%c%c%c", n, tmp[0], tmp[1], tmp[2], tmp[3]);
190 if (me == 0) {
191 mkdir(n, 0777);
192 }
193 sprintf(n, "%s/%c%c", n, tmp[4], tmp[5]);
194
195 if (me == 0) {
196 mkdir(n, 0777);
197 }
198}
199void
200vtkfile(int step, const hydroparam_t H, hydrovar_t * Hv) {
201 char name[1024];
202 char vfrname[1024];
203 FILE *fic, *vf;
204 int i, j, nv;
205
206 WHERE("vtkfile");
207
208 // First step : create the directory structure ONLY using PE0
209#ifdef MPI
210#if FTI==0
211 if (H.nproc > 1) MPI_Barrier(MPI_COMM_WORLD);
212#endif
213#if FTI>0
214 if (H.nproc > 1) MPI_Barrier(FTI_COMM_WORLD);
215#endif
216#endif
217 vtknm(vfrname, H.mype, step); // create the directory structure
218 // if (H.mype == 0) fprintf(stderr, "%s\n", vfrname);
219#ifdef MPI
220#if FTI==0
221 if (H.nproc > 1) MPI_Barrier(MPI_COMM_WORLD);
222#endif
223#if FTI>0
224 if (H.nproc > 1) MPI_Barrier(FTI_COMM_WORLD);
225#endif
226#endif
227
228 // Write a domain per PE
229 sprintf(name, "%s/Hydro_%05d_%04d.vtr", vfrname, H.mype, step);
230 fic = fopen(name, "w");
231 if (fic == NULL) {
232 fprintf(stderr, "Ouverture du fichier %s impossible\n", name);
233 exit(1);
234 }
235 fprintf(fic, "<?xml version=\"1.0\"?>\n");
236 fprintf(fic, "<VTKFile type=\"RectilinearGrid\" byte_order=\"LittleEndian\">\n");
237 fprintf(fic, " <RectilinearGrid WholeExtent=\" %d %d %d %d %d %d\">\n",
238 H.box[XMIN_BOX], H.box[XMAX_BOX], H.box[YMIN_BOX], H.box[YMAX_BOX], 0, 1);
239 fprintf(fic, " <Piece Extent=\" %d %d %d %d %d %d\" GhostLevel=\"0\">\n",
240 H.box[XMIN_BOX], H.box[XMAX_BOX], H.box[YMIN_BOX], H.box[YMAX_BOX], 0, 1);
241 fprintf(fic, " <Coordinates>\n");
242
243 fprintf(fic, " <DataArray type=\"Float32\" format=\"ascii\" NumberOfComponents=\"1\">\n");
244 for (i = H.box[XMIN_BOX]; i <= H.box[XMAX_BOX]; i++) {
245 fprintf(fic, "%f ", i * H.dx);
246 }
247 fprintf(fic, "\n");
248 fprintf(fic, " </DataArray>\n");
249 fprintf(fic, " <DataArray type=\"Float32\" format=\"ascii\" NumberOfComponents=\"1\">\n");
250 for (j = H.box[YMIN_BOX]; j <= H.box[YMAX_BOX]; j++) {
251 fprintf(fic, "%f ", j * H.dx);
252 }
253 fprintf(fic, "\n");
254 fprintf(fic, " </DataArray>\n");
255 fprintf(fic, " <DataArray type=\"Float32\" format=\"ascii\" NumberOfComponents=\"1\">\n");
256 fprintf(fic, "%f %f\n", 0., 1. * H.dx);
257 fprintf(fic, " </DataArray>\n");
258 fprintf(fic, " </Coordinates>\n");
259 name[0] = 0;
260 for (nv = 0; nv <= IP; nv++) {
261 if (nv == ID)
262 sprintf(name, "%s varID", name);
263 if (nv == IU)
264 sprintf(name, "%s varIU", name);
265 if (nv == IV)
266 sprintf(name, "%s varIV", name);
267 if (nv == IP)
268 sprintf(name, "%s varIP", name);
269 }
270
271 // declaration of the variable list
272 fprintf(fic, " <CellData Scalars=\"%s\">\n", name);
273 name[0] = 0;
274 for (nv = 0; nv <= IP; nv++) {
275 if (nv == ID)
276 sprintf(name, "varID");
277 if (nv == IU)
278 sprintf(name, "varIU");
279 if (nv == IV)
280 sprintf(name, "varIV");
281 if (nv == IP)
282 sprintf(name, "varIP");
283
284 //Definition of the cell values
285#if BINARY == 1
286 fprintf(fic,
287 " <DataArray Name=\"%s\" type=\"Float32\" format=\"binary\" encoding=\"base64\" NumberOfComponents=\"1\">\n",
288 name);
289 {
290 // float tuold[H.nxt * H.nyt];
291 float *tuold = NULL;
292 char *r64;
293 size_t p = 0, lst;
294
295 assert((H.nxt * H.nyt) > 0);
296 tuold = (float *) calloc(H.nxt * H.nyt + 16, sizeof(float));
297 assert(tuold != NULL);
298
299 for (j = H.jmin + ExtraLayer; j < H.jmax - ExtraLayer; j++) {
300 for (i = H.imin + ExtraLayer; i < H.imax - ExtraLayer; i++) {
301 tuold[p++] = (float) Hv->uold[IHv(i, j, nv)];
302 }
303 }
304 // Header = size of the following items
305 assert(p < H.nxt * H.nyt);
306
307 p *= sizeof(float);
308 r64 = ToBase64((byte *) & p, sizeof(int));
309 lst = strlen(r64);
310 fwrite(r64, 1, lst, fic);
311 free(r64);
312 r64 = ToBase64((byte *) tuold, p);
313 lst = strlen(r64);
314 fwrite(r64, 1, lst, fic);
315 free(r64);
316 free(tuold);
317 }
318#else
319 fprintf(fic, " <DataArray type=\"Float32\" Name=\"%s\" format=\"ascii\" NumberOfComponents=\"1\">\n", name);
320
321 // the image is the interior of the computed domain
322 for (j = H.jmin + ExtraLayer; j < H.jmax - ExtraLayer; j++) {
323 for (i = H.imin + ExtraLayer; i < H.imax - ExtraLayer; i++) {
324 fprintf(fic, "%lf ", Hv->uold[IHv(i, j, nv)]);
325 }
326 fprintf(fic, "\n");
327 }
328#endif
329 fprintf(fic, " </DataArray>\n");
330 }
331 fprintf(fic, " </CellData>\n");
332 fprintf(fic, " </Piece>\n");
333 fprintf(fic, " </RectilinearGrid>\n");
334 fprintf(fic, "</VTKFile>\n");
335 fclose(fic);
336
337 // At this stage we can write VTK containers. Since only one file is
338 // necessary even for multiple domains, it has to be written by one
339 // PE only.
340
341#ifdef MPI
342#if FTI==0
343 if (H.nproc > 1) MPI_Barrier(MPI_COMM_WORLD);
344#endif
345#if FTI>0
346 if (H.nproc > 1) MPI_Barrier(FTI_COMM_WORLD);
347#endif
348#endif
349 if (H.mype == 0) {
350 sprintf(name, "outputvtk_%05d.pvtr", step);
351 sprintf(name, "%s/Hydro_%04d.pvtr", vfrname, step);
352 vf = fopen(name, "w");
353 if (vf == NULL) {
354 fprintf(stderr, "Ouverture du fichier %s impossible\n", name);
355 exit(1);
356 }
357 fprintf(vf, "<?xml version=\"1.0\"?>\n");
358 fprintf(vf, "<VTKFile type=\"PRectilinearGrid\" byte_order=\"LittleEndian\">\n");
359 fprintf(vf, "<PRectilinearGrid WholeExtent=\"0 %d 0 %d 0 %d\" GhostLevel=\"0\" >\n", H.globnx, H.globny, 1);
360 fprintf(vf, " <PCellData>\n");
361 for (nv = 0; nv <= IP; nv++) {
362 name[0] = '\0';
363 if (nv == ID)
364 sprintf(name, "varID");
365 if (nv == IU)
366 sprintf(name, "varIU");
367 if (nv == IV)
368 sprintf(name, "varIV");
369 if (nv == IP)
370 sprintf(name, "varIP");
371
372#if BINARY == 1
373 fprintf(vf,
374 " <PDataArray Name=\"%s\" type=\"Float32\" format=\"binary\" encoding=\"base64\" NumberOfComponents=\"1\"/>\n",
375 name);
376#else
377 fprintf(vf, " <PDataArray Name=\"%s\" type=\"Float32\" format=\"ascii\" NumberOfComponents=\"1\"/>\n", name);
378#endif
379 }
380 fprintf(vf, " </PCellData>\n");
381 fprintf(vf, " <PCoordinates>\n");
382 fprintf(vf, " <PDataArray type=\"Float32\" format=\"ascii\" NumberOfComponents=\"1\"/>\n");
383 fprintf(vf, " <PDataArray type=\"Float32\" format=\"ascii\" NumberOfComponents=\"1\"/>\n");
384 fprintf(vf, " <PDataArray type=\"Float32\" format=\"ascii\" NumberOfComponents=\"1\"/>\n");
385 fprintf(vf, " </PCoordinates>\n");
386 for (i = 0; i < H.nproc; i++) {
387 int box[8];
388 memset(box, 0, 8 * sizeof(int));
389 CalcSubSurface(0, H.globnx, 0, H.globny, 0, H.nproc - 1, 0, box, i, 0);
390 sprintf(name, "Hydro_%05d_%04d.vtr", i, step);
391 fprintf(vf, " <Piece Extent=\"%d %d %d %d %d %d\" Source=\"%s\"/>\n", box[XMIN_BOX],
392 box[XMAX_BOX], box[YMIN_BOX], box[YMAX_BOX], 0, 1, name);
393 }
394 fprintf(vf, "</PRectilinearGrid>\n");
395 fprintf(vf, "</VTKFile>\n");
396 fclose(vf);
397
398 // We make the time step available only now to ensure consistency
399 vtkwpvd(step, "Dep");
400 }
401}
Note: See TracBrowser for help on using the repository browser.