tfile_io.cpp - sphere - GPU-based 3D discrete element method algorithm with optional fluid coupling
HTML git clone git://src.adamsgaard.dk/sphere
DIR Log
DIR Files
DIR Refs
DIR LICENSE
---
tfile_io.cpp (32869B)
---
1 #include <iostream>
2 #include <fstream>
3 #include <cstdio>
4 #include <cstdlib>
5
6 #include "typedefs.h"
7 #include "datatypes.h"
8 #include "constants.h"
9 #include "sphere.h"
10
11
12 // Get the address of the first byte of an object's representation
13 // See Stroustrup (2008) p. 388
14 template<class T>
15 char* as_bytes(T& i) // treat a T as a sequence of bytes
16 {
17 // get the address of the first byte of memory used
18 // to store the object
19 void* addr = &i;
20
21 // treat the object as bytes
22 return static_cast<char*>(addr);
23 }
24
25 // Read DEM data from binary file
26 // Note: Static-size arrays can be bulk-read with e.g.
27 // ifs.read(as_bytes(grid.L), sizeof(grid.L))
28 // while dynamic, and vector arrays (e.g. Float4) must
29 // be read one value at a time.
30 void DEM::readbin(const char *target)
31 {
32 using std::cout; // stdout
33 using std::cerr; // stderr
34 using std::endl; // endline. Implicitly flushes buffer
35 unsigned int i;
36
37 // Open input file
38 // if target is string:
39 // std::ifstream ifs(target.c_str(), std::ios_base::binary);
40 std::ifstream ifs(target, std::ios_base::binary);
41 if (!ifs) {
42 cerr << "Could not read input binary file '"
43 << target << "'" << endl;
44 exit(1);
45 }
46
47 Float version;
48 ifs.read(as_bytes(version), sizeof(Float));
49 if (version != VERSION) {
50 std::cerr << "Error: The input file '" << target << "' is written by "
51 "sphere version " << version << ", which is incompatible with this "
52 "version (" << VERSION << ")." << std::endl;
53 exit(1);
54 }
55
56 ifs.read(as_bytes(nd), sizeof(nd));
57 ifs.read(as_bytes(np), sizeof(np));
58
59 if (nd != ND) {
60 cerr << "Dimensionality mismatch between dataset and this SPHERE "
61 "program.\nThe dataset is " << nd
62 << "D, this SPHERE binary is " << ND << "D.\n"
63 << "This execution is terminating." << endl;
64 exit(-1); // Return unsuccessful exit status
65 }
66
67 // Check precision choice
68 if (sizeof(Float) != sizeof(double) && sizeof(Float) != sizeof(float)) {
69 cerr << "Error! Chosen precision not available. Check datatypes.h\n";
70 exit(1);
71 }
72
73 // Read time parameters
74 ifs.read(as_bytes(time.dt), sizeof(time.dt));
75 ifs.read(as_bytes(time.current), sizeof(time.current));
76 ifs.read(as_bytes(time.total), sizeof(time.total));
77 ifs.read(as_bytes(time.file_dt), sizeof(time.file_dt));
78 ifs.read(as_bytes(time.step_count), sizeof(time.step_count));
79
80 // For spatial vectors an array of Float4 vectors is chosen for best fit
81 // with GPU memory handling. Vector variable structure: ( x, y, z, <empty>).
82 // Indexing starts from 0.
83
84 // Allocate host arrays
85 if (verbose == 1)
86 cout << " Allocating host memory: ";
87 // Allocate more host arrays
88 k.x = new Float4[np];
89 k.xyzsum = new Float4[np];
90 k.vel = new Float4[np];
91 k.force = new Float4[np];
92 k.angpos = new Float4[np];
93 k.angvel = new Float4[np];
94 k.torque = new Float4[np];
95 k.color = new int[np];
96
97 e.es_dot = new Float[np];
98 e.es = new Float[np];
99 e.ev_dot = new Float[np];
100 e.ev = new Float[np];
101 e.p = new Float[np];
102
103 if (verbose == 1)
104 cout << "Done\n";
105
106 if (verbose == 1)
107 cout << " Reading remaining data from input binary: ";
108
109 // Read grid parameters
110 ifs.read(as_bytes(grid.origo), sizeof(grid.origo));
111 ifs.read(as_bytes(grid.L), sizeof(grid.L));
112 ifs.read(as_bytes(grid.num), sizeof(grid.num));
113 ifs.read(as_bytes(grid.periodic), sizeof(grid.periodic));
114 ifs.read(as_bytes(grid.adaptive), sizeof(grid.adaptive));
115
116 // Read kinematic values
117 for (i = 0; i<np; ++i) {
118 ifs.read(as_bytes(k.x[i].x), sizeof(Float));
119 ifs.read(as_bytes(k.x[i].y), sizeof(Float));
120 ifs.read(as_bytes(k.x[i].z), sizeof(Float));
121 ifs.read(as_bytes(k.x[i].w), sizeof(Float));
122 }
123 for (i = 0; i<np; ++i) {
124 ifs.read(as_bytes(k.xyzsum[i].x), sizeof(Float));
125 ifs.read(as_bytes(k.xyzsum[i].y), sizeof(Float));
126 ifs.read(as_bytes(k.xyzsum[i].z), sizeof(Float));
127 }
128 for (i = 0; i<np; ++i) {
129 ifs.read(as_bytes(k.vel[i].x), sizeof(Float));
130 ifs.read(as_bytes(k.vel[i].y), sizeof(Float));
131 ifs.read(as_bytes(k.vel[i].z), sizeof(Float));
132 ifs.read(as_bytes(k.vel[i].w), sizeof(Float));
133 }
134 for (i = 0; i<np; ++i) {
135 ifs.read(as_bytes(k.force[i].x), sizeof(Float));
136 ifs.read(as_bytes(k.force[i].y), sizeof(Float));
137 ifs.read(as_bytes(k.force[i].z), sizeof(Float));
138 //ifs.read(as_bytes(k.force[i].w), sizeof(Float));
139 }
140 for (i = 0; i<np; ++i) {
141 ifs.read(as_bytes(k.angpos[i].x), sizeof(Float));
142 ifs.read(as_bytes(k.angpos[i].y), sizeof(Float));
143 ifs.read(as_bytes(k.angpos[i].z), sizeof(Float));
144 //ifs.read(as_bytes(k.angpos[i].w), sizeof(Float));
145 }
146 for (i = 0; i<np; ++i) {
147 ifs.read(as_bytes(k.angvel[i].x), sizeof(Float));
148 ifs.read(as_bytes(k.angvel[i].y), sizeof(Float));
149 ifs.read(as_bytes(k.angvel[i].z), sizeof(Float));
150 //ifs.read(as_bytes(k.angvel[i].w), sizeof(Float));
151 }
152 for (i = 0; i<np; ++i) {
153 ifs.read(as_bytes(k.torque[i].x), sizeof(Float));
154 ifs.read(as_bytes(k.torque[i].y), sizeof(Float));
155 ifs.read(as_bytes(k.torque[i].z), sizeof(Float));
156 //ifs.read(as_bytes(k.torque[i].w), sizeof(Float));
157 }
158
159 // Read energies
160 for (i = 0; i<np; ++i)
161 ifs.read(as_bytes(e.es_dot[i]), sizeof(Float));
162 for (i = 0; i<np; ++i)
163 ifs.read(as_bytes(e.es[i]), sizeof(Float));
164 for (i = 0; i<np; ++i)
165 ifs.read(as_bytes(e.ev_dot[i]), sizeof(Float));
166 for (i = 0; i<np; ++i)
167 ifs.read(as_bytes(e.ev[i]), sizeof(Float));
168 for (i = 0; i<np; ++i)
169 ifs.read(as_bytes(e.p[i]), sizeof(Float));
170
171 // Read constant parameters
172 ifs.read(as_bytes(params.g), sizeof(params.g));
173 ifs.read(as_bytes(params.k_n), sizeof(params.k_n));
174 ifs.read(as_bytes(params.k_t), sizeof(params.k_t));
175 ifs.read(as_bytes(params.k_r), sizeof(params.k_r));
176 ifs.read(as_bytes(params.E), sizeof(params.E));
177 ifs.read(as_bytes(params.gamma_n), sizeof(params.gamma_n));
178 ifs.read(as_bytes(params.gamma_t), sizeof(params.gamma_t));
179 ifs.read(as_bytes(params.gamma_r), sizeof(params.gamma_r));
180 ifs.read(as_bytes(params.mu_s), sizeof(params.mu_s));
181 ifs.read(as_bytes(params.mu_d), sizeof(params.mu_d));
182 ifs.read(as_bytes(params.mu_r), sizeof(params.mu_r));
183 ifs.read(as_bytes(params.gamma_wn), sizeof(params.gamma_wn));
184 ifs.read(as_bytes(params.gamma_wt), sizeof(params.gamma_wt));
185 ifs.read(as_bytes(params.mu_ws), sizeof(params.mu_s));
186 ifs.read(as_bytes(params.mu_wd), sizeof(params.mu_d));
187 ifs.read(as_bytes(params.rho), sizeof(params.rho));
188 ifs.read(as_bytes(params.contactmodel), sizeof(params.contactmodel));
189 ifs.read(as_bytes(params.kappa), sizeof(params.kappa));
190 ifs.read(as_bytes(params.db), sizeof(params.db));
191 ifs.read(as_bytes(params.V_b), sizeof(params.V_b));
192
193 // Read wall parameters
194 ifs.read(as_bytes(walls.nw), sizeof(walls.nw));
195 if (walls.nw > MAXWALLS) {
196 cerr << "Error; MAXWALLS (" << MAXWALLS << ") in datatypes.h "
197 << "is smaller than the number of walls specified in the "
198 << "input file (" << walls.nw << ").\n";
199 exit(1);
200 }
201
202 // Allocate host memory for walls
203 // Wall normal (x,y,z), w: wall position on axis parallel to wall normal
204 // Wall mass (x), velocity (y), force (z), and deviatoric stress (w)
205 walls.nx = new Float4[walls.nw];
206 walls.mvfd = new Float4[walls.nw];
207 walls.tau_x = new Float[1];
208
209 for (i = 0; i<walls.nw; ++i)
210 ifs.read(as_bytes(walls.wmode[i]), sizeof(walls.wmode[i]));
211 for (i = 0; i<walls.nw; ++i) {
212 ifs.read(as_bytes(walls.nx[i].x), sizeof(Float));
213 ifs.read(as_bytes(walls.nx[i].y), sizeof(Float));
214 ifs.read(as_bytes(walls.nx[i].z), sizeof(Float));
215 ifs.read(as_bytes(walls.nx[i].w), sizeof(Float));
216 }
217 for (i = 0; i<walls.nw; ++i) {
218 ifs.read(as_bytes(walls.mvfd[i].x), sizeof(Float));
219 ifs.read(as_bytes(walls.mvfd[i].y), sizeof(Float));
220 ifs.read(as_bytes(walls.mvfd[i].z), sizeof(Float));
221 ifs.read(as_bytes(walls.mvfd[i].w), sizeof(Float));
222 }
223 ifs.read(as_bytes(params.sigma0_A), sizeof(params.sigma0_A));
224 ifs.read(as_bytes(params.sigma0_f), sizeof(params.sigma0_f));
225 ifs.read(as_bytes(walls.tau_x[0]), sizeof(walls.tau_x[0]));
226
227 // Read bond parameters
228 ifs.read(as_bytes(params.lambda_bar), sizeof(params.lambda_bar));
229 ifs.read(as_bytes(params.nb0), sizeof(params.nb0));
230 ifs.read(as_bytes(params.sigma_b), sizeof(params.sigma_b));
231 ifs.read(as_bytes(params.tau_b), sizeof(params.tau_b));
232 if (params.nb0 > 0)
233 k.bonds = new uint2[params.nb0];
234 k.bonds_delta = new Float4[np];
235 k.bonds_omega = new Float4[np];
236 for (i = 0; i<params.nb0; ++i) {
237 ifs.read(as_bytes(k.bonds[i].x), sizeof(unsigned int));
238 ifs.read(as_bytes(k.bonds[i].y), sizeof(unsigned int));
239 }
240 for (i = 0; i<params.nb0; ++i) // Normal component
241 ifs.read(as_bytes(k.bonds_delta[i].w), sizeof(Float));
242 for (i = 0; i<params.nb0; ++i) { // Tangential component
243 ifs.read(as_bytes(k.bonds_delta[i].x), sizeof(Float));
244 ifs.read(as_bytes(k.bonds_delta[i].y), sizeof(Float));
245 ifs.read(as_bytes(k.bonds_delta[i].z), sizeof(Float));
246 }
247 for (i = 0; i<params.nb0; ++i) // Normal component
248 ifs.read(as_bytes(k.bonds_omega[i].w), sizeof(Float));
249 for (i = 0; i<params.nb0; ++i) { // Tangential component
250 ifs.read(as_bytes(k.bonds_omega[i].x), sizeof(Float));
251 ifs.read(as_bytes(k.bonds_omega[i].y), sizeof(Float));
252 ifs.read(as_bytes(k.bonds_omega[i].z), sizeof(Float));
253 }
254
255 unsigned int x, y, z;
256
257 if (verbose == 1)
258 cout << "Done\n";
259
260 // Simulate fluid
261 if (fluid == 1) {
262
263 ifs.read(as_bytes(cfd_solver), sizeof(int));
264
265 if (cfd_solver < 0 || cfd_solver > 1) {
266 std::cerr << "Value of cfd_solver not understood ("
267 << cfd_solver << ")" << std::endl;
268 exit(1);
269 }
270
271 if (cfd_solver == 0) { // Navier Stokes flow
272
273 initNSmem();
274
275 ifs.read(as_bytes(ns.mu), sizeof(Float));
276
277 if (verbose == 1)
278 cout << " - Reading fluid values:\t\t\t ";
279
280 for (z = 0; z<grid.num[2]; ++z) {
281 for (y = 0; y<grid.num[1]; ++y) {
282 for (x = 0; x<grid.num[0]; ++x) {
283 i = idx(x,y,z);
284 ifs.read(as_bytes(ns.v[i].x), sizeof(Float));
285 ifs.read(as_bytes(ns.v[i].y), sizeof(Float));
286 ifs.read(as_bytes(ns.v[i].z), sizeof(Float));
287 ifs.read(as_bytes(ns.p[i]), sizeof(Float));
288 ifs.read(as_bytes(ns.phi[i]), sizeof(Float));
289 ifs.read(as_bytes(ns.dphi[i]), sizeof(Float));
290 }
291 }
292 }
293
294 ifs.read(as_bytes(ns.rho_f), sizeof(Float));
295 ifs.read(as_bytes(ns.p_mod_A), sizeof(Float));
296 ifs.read(as_bytes(ns.p_mod_f), sizeof(Float));
297 ifs.read(as_bytes(ns.p_mod_phi), sizeof(Float));
298
299 ifs.read(as_bytes(ns.bc_top), sizeof(int));
300 ifs.read(as_bytes(ns.bc_bot), sizeof(int));
301 ifs.read(as_bytes(ns.free_slip_bot), sizeof(int));
302 ifs.read(as_bytes(ns.free_slip_top), sizeof(int));
303 ifs.read(as_bytes(ns.bc_bot_flux), sizeof(Float));
304 ifs.read(as_bytes(ns.bc_top_flux), sizeof(Float));
305
306 for (z = 0; z<grid.num[2]; ++z)
307 for (y = 0; y<grid.num[1]; ++y)
308 for (x = 0; x<grid.num[0]; ++x)
309 ifs.read(as_bytes(ns.p_constant[idx(x,y,z)]),
310 sizeof(int));
311
312 ifs.read(as_bytes(ns.gamma), sizeof(Float));
313 ifs.read(as_bytes(ns.theta), sizeof(Float));
314 ifs.read(as_bytes(ns.beta), sizeof(Float));
315 ifs.read(as_bytes(ns.tolerance), sizeof(Float));
316 ifs.read(as_bytes(ns.maxiter), sizeof(unsigned int));
317 ifs.read(as_bytes(ns.ndem), sizeof(unsigned int));
318
319 ifs.read(as_bytes(ns.c_phi), sizeof(Float));
320 ifs.read(as_bytes(ns.c_v), sizeof(Float));
321 ifs.read(as_bytes(ns.dt_dem_fac), sizeof(Float));
322
323 for (i = 0; i<np; ++i) {
324 ifs.read(as_bytes(ns.f_d[i].x), sizeof(Float));
325 ifs.read(as_bytes(ns.f_d[i].y), sizeof(Float));
326 ifs.read(as_bytes(ns.f_d[i].z), sizeof(Float));
327 }
328 for (i = 0; i<np; ++i) {
329 ifs.read(as_bytes(ns.f_p[i].x), sizeof(Float));
330 ifs.read(as_bytes(ns.f_p[i].y), sizeof(Float));
331 ifs.read(as_bytes(ns.f_p[i].z), sizeof(Float));
332 }
333 for (i = 0; i<np; ++i) {
334 ifs.read(as_bytes(ns.f_v[i].x), sizeof(Float));
335 ifs.read(as_bytes(ns.f_v[i].y), sizeof(Float));
336 ifs.read(as_bytes(ns.f_v[i].z), sizeof(Float));
337 }
338 for (i = 0; i<np; ++i) {
339 ifs.read(as_bytes(ns.f_sum[i].x), sizeof(Float));
340 ifs.read(as_bytes(ns.f_sum[i].y), sizeof(Float));
341 ifs.read(as_bytes(ns.f_sum[i].z), sizeof(Float));
342 }
343
344 if (verbose == 1)
345 cout << "Done" << std::endl;
346
347 } else if (cfd_solver == 1) { // Darcy flow
348
349 initDarcyMem();
350
351 ifs.read(as_bytes(darcy.mu), sizeof(Float));
352
353 if (verbose == 1)
354 cout << " - Reading fluid values:\t\t\t ";
355
356 for (z = 0; z<darcy.nz; ++z) {
357 for (y = 0; y<darcy.ny; ++y) {
358 for (x = 0; x<darcy.nx; ++x) {
359 i = d_idx(x,y,z);
360 ifs.read(as_bytes(darcy.v[i].x), sizeof(Float));
361 ifs.read(as_bytes(darcy.v[i].y), sizeof(Float));
362 ifs.read(as_bytes(darcy.v[i].z), sizeof(Float));
363 ifs.read(as_bytes(darcy.p[i]), sizeof(Float));
364 ifs.read(as_bytes(darcy.phi[i]), sizeof(Float));
365 ifs.read(as_bytes(darcy.dphi[i]), sizeof(Float));
366 }
367 }
368 }
369
370 ifs.read(as_bytes(darcy.rho_f), sizeof(Float));
371 ifs.read(as_bytes(darcy.p_mod_A), sizeof(Float));
372 ifs.read(as_bytes(darcy.p_mod_f), sizeof(Float));
373 ifs.read(as_bytes(darcy.p_mod_phi), sizeof(Float));
374
375 ifs.read(as_bytes(darcy.bc_xn), sizeof(int));
376 ifs.read(as_bytes(darcy.bc_xp), sizeof(int));
377 ifs.read(as_bytes(darcy.bc_yn), sizeof(int));
378 ifs.read(as_bytes(darcy.bc_yp), sizeof(int));
379 ifs.read(as_bytes(darcy.bc_bot), sizeof(int));
380 ifs.read(as_bytes(darcy.bc_top), sizeof(int));
381 ifs.read(as_bytes(darcy.free_slip_bot), sizeof(int));
382 ifs.read(as_bytes(darcy.free_slip_top), sizeof(int));
383 ifs.read(as_bytes(darcy.bc_bot_flux), sizeof(Float));
384 ifs.read(as_bytes(darcy.bc_top_flux), sizeof(Float));
385
386 for (z = 0; z<darcy.nz; ++z)
387 for (y = 0; y<darcy.ny; ++y)
388 for (x = 0; x<darcy.nx; ++x)
389 ifs.read(as_bytes(darcy.p_constant[d_idx(x,y,z)]),
390 sizeof(int));
391
392 ifs.read(as_bytes(darcy.tolerance), sizeof(Float));
393 ifs.read(as_bytes(darcy.maxiter), sizeof(unsigned int));
394 ifs.read(as_bytes(darcy.ndem), sizeof(unsigned int));
395 ifs.read(as_bytes(darcy.c_phi), sizeof(Float));
396
397 for (i = 0; i<np; ++i) {
398 ifs.read(as_bytes(darcy.f_p[i].x), sizeof(Float));
399 ifs.read(as_bytes(darcy.f_p[i].y), sizeof(Float));
400 ifs.read(as_bytes(darcy.f_p[i].z), sizeof(Float));
401 }
402
403 ifs.read(as_bytes(darcy.beta_f), sizeof(Float));
404 ifs.read(as_bytes(darcy.k_c), sizeof(Float));
405
406 if (verbose == 1)
407 cout << "Done" << std::endl;
408 }
409 }
410
411 for (i = 0; i<np; ++i)
412 ifs.read(as_bytes(k.color[i]), sizeof(int));
413
414 // Close file if it is still open
415 if (ifs.is_open())
416 ifs.close();
417 }
418
419 // Write DEM data to binary file
420 void DEM::writebin(const char *target)
421 {
422 unsigned int i;
423
424 // Open output file
425 std::ofstream ofs(target, std::ios_base::binary);
426 if (!ofs) {
427 std::cerr << "could not create output binary file '"
428 << target << "'" << std::endl;
429 exit(1);
430 }
431
432 // If double precision: Values can be written directly
433 if (sizeof(Float) == sizeof(double)) {
434
435 double version = VERSION;
436 ofs.write(as_bytes(version), sizeof(Float));
437
438 ofs.write(as_bytes(nd), sizeof(nd));
439 ofs.write(as_bytes(np), sizeof(np));
440
441 // Write time parameters
442 ofs.write(as_bytes(time.dt), sizeof(time.dt));
443 ofs.write(as_bytes(time.current), sizeof(time.current));
444 ofs.write(as_bytes(time.total), sizeof(time.total));
445 ofs.write(as_bytes(time.file_dt), sizeof(time.file_dt));
446 ofs.write(as_bytes(time.step_count), sizeof(time.step_count));
447
448 // Write grid parameters
449 ofs.write(as_bytes(grid.origo), sizeof(grid.origo));
450 ofs.write(as_bytes(grid.L), sizeof(grid.L));
451 ofs.write(as_bytes(grid.num), sizeof(grid.num));
452 ofs.write(as_bytes(grid.periodic), sizeof(grid.periodic));
453 ofs.write(as_bytes(grid.adaptive), sizeof(grid.adaptive));
454
455 // Write kinematic values
456 for (i = 0; i<np; ++i) {
457 ofs.write(as_bytes(k.x[i].x), sizeof(Float));
458 ofs.write(as_bytes(k.x[i].y), sizeof(Float));
459 ofs.write(as_bytes(k.x[i].z), sizeof(Float));
460 ofs.write(as_bytes(k.x[i].w), sizeof(Float));
461 }
462 for (i = 0; i<np; ++i) {
463 ofs.write(as_bytes(k.xyzsum[i].x), sizeof(Float));
464 ofs.write(as_bytes(k.xyzsum[i].y), sizeof(Float));
465 ofs.write(as_bytes(k.xyzsum[i].z), sizeof(Float));
466 }
467 for (i = 0; i<np; ++i) {
468 ofs.write(as_bytes(k.vel[i].x), sizeof(Float));
469 ofs.write(as_bytes(k.vel[i].y), sizeof(Float));
470 ofs.write(as_bytes(k.vel[i].z), sizeof(Float));
471 ofs.write(as_bytes(k.vel[i].w), sizeof(Float));
472 }
473 for (i = 0; i<np; ++i) {
474 ofs.write(as_bytes(k.force[i].x), sizeof(Float));
475 ofs.write(as_bytes(k.force[i].y), sizeof(Float));
476 ofs.write(as_bytes(k.force[i].z), sizeof(Float));
477 //ofs.write(as_bytes(k.force[i].w), sizeof(Float));
478 }
479 for (i = 0; i<np; ++i) {
480 ofs.write(as_bytes(k.angpos[i].x), sizeof(Float));
481 ofs.write(as_bytes(k.angpos[i].y), sizeof(Float));
482 ofs.write(as_bytes(k.angpos[i].z), sizeof(Float));
483 //ofs.write(as_bytes(k.angpos[i].w), sizeof(Float));
484 }
485 for (i = 0; i<np; ++i) {
486 ofs.write(as_bytes(k.angvel[i].x), sizeof(Float));
487 ofs.write(as_bytes(k.angvel[i].y), sizeof(Float));
488 ofs.write(as_bytes(k.angvel[i].z), sizeof(Float));
489 //ofs.write(as_bytes(k.angvel[i].w), sizeof(Float));
490 }
491 for (i = 0; i<np; ++i) {
492 ofs.write(as_bytes(k.torque[i].x), sizeof(Float));
493 ofs.write(as_bytes(k.torque[i].y), sizeof(Float));
494 ofs.write(as_bytes(k.torque[i].z), sizeof(Float));
495 //ofs.write(as_bytes(k.torque[i].w), sizeof(Float));
496 }
497
498 // Write energies
499 for (i = 0; i<np; ++i)
500 ofs.write(as_bytes(e.es_dot[i]), sizeof(Float));
501 for (i = 0; i<np; ++i)
502 ofs.write(as_bytes(e.es[i]), sizeof(Float));
503 for (i = 0; i<np; ++i)
504 ofs.write(as_bytes(e.ev_dot[i]), sizeof(Float));
505 for (i = 0; i<np; ++i)
506 ofs.write(as_bytes(e.ev[i]), sizeof(Float));
507 for (i = 0; i<np; ++i)
508 ofs.write(as_bytes(e.p[i]), sizeof(Float));
509
510 // Write constant parameters
511 ofs.write(as_bytes(params.g), sizeof(params.g));
512 ofs.write(as_bytes(params.k_n), sizeof(params.k_n));
513 ofs.write(as_bytes(params.k_t), sizeof(params.k_t));
514 ofs.write(as_bytes(params.k_r), sizeof(params.k_r));
515 ofs.write(as_bytes(params.E), sizeof(params.E));
516 ofs.write(as_bytes(params.gamma_n), sizeof(params.gamma_n));
517 ofs.write(as_bytes(params.gamma_t), sizeof(params.gamma_t));
518 ofs.write(as_bytes(params.gamma_r), sizeof(params.gamma_r));
519 ofs.write(as_bytes(params.mu_s), sizeof(params.mu_s));
520 ofs.write(as_bytes(params.mu_d), sizeof(params.mu_d));
521 ofs.write(as_bytes(params.mu_r), sizeof(params.mu_r));
522 ofs.write(as_bytes(params.gamma_wn), sizeof(params.gamma_wn));
523 ofs.write(as_bytes(params.gamma_wt), sizeof(params.gamma_wt));
524 ofs.write(as_bytes(params.mu_ws), sizeof(params.mu_ws));
525 ofs.write(as_bytes(params.mu_wd), sizeof(params.mu_wd));
526 ofs.write(as_bytes(params.rho), sizeof(params.rho));
527 ofs.write(as_bytes(params.contactmodel), sizeof(params.contactmodel));
528 ofs.write(as_bytes(params.kappa), sizeof(params.kappa));
529 ofs.write(as_bytes(params.db), sizeof(params.db));
530 ofs.write(as_bytes(params.V_b), sizeof(params.V_b));
531
532 // Write wall parameters
533 ofs.write(as_bytes(walls.nw), sizeof(walls.nw));
534 ofs.write(as_bytes(walls.wmode), sizeof(walls.wmode[0])*walls.nw);
535 for (i = 0; i<walls.nw; ++i) {
536 ofs.write(as_bytes(walls.nx[i].x), sizeof(Float));
537 ofs.write(as_bytes(walls.nx[i].y), sizeof(Float));
538 ofs.write(as_bytes(walls.nx[i].z), sizeof(Float));
539 ofs.write(as_bytes(walls.nx[i].w), sizeof(Float));
540 }
541 for (i = 0; i<walls.nw; ++i) {
542 ofs.write(as_bytes(walls.mvfd[i].x), sizeof(Float));
543 ofs.write(as_bytes(walls.mvfd[i].y), sizeof(Float));
544 ofs.write(as_bytes(walls.mvfd[i].z), sizeof(Float));
545 ofs.write(as_bytes(walls.mvfd[i].w), sizeof(Float));
546 }
547 ofs.write(as_bytes(params.sigma0_A), sizeof(params.sigma0_A));
548 ofs.write(as_bytes(params.sigma0_f), sizeof(params.sigma0_f));
549 ofs.write(as_bytes(walls.tau_x[0]), sizeof(walls.tau_x[0]));
550
551 // Write bond parameters
552 ofs.write(as_bytes(params.lambda_bar), sizeof(params.lambda_bar));
553 ofs.write(as_bytes(params.nb0), sizeof(params.nb0));
554 ofs.write(as_bytes(params.sigma_b), sizeof(params.sigma_b));
555 ofs.write(as_bytes(params.tau_b), sizeof(params.tau_b));
556 for (i = 0; i<params.nb0; ++i) {
557 ofs.write(as_bytes(k.bonds[i].x), sizeof(unsigned int));
558 ofs.write(as_bytes(k.bonds[i].y), sizeof(unsigned int));
559 }
560 for (i = 0; i<params.nb0; ++i) // Normal component
561 ofs.write(as_bytes(k.bonds_delta[i].w), sizeof(Float));
562 for (i = 0; i<params.nb0; ++i) { // Tangential component
563 ofs.write(as_bytes(k.bonds_delta[i].x), sizeof(Float));
564 ofs.write(as_bytes(k.bonds_delta[i].y), sizeof(Float));
565 ofs.write(as_bytes(k.bonds_delta[i].z), sizeof(Float));
566 }
567 for (i = 0; i<params.nb0; ++i) // Normal component
568 ofs.write(as_bytes(k.bonds_omega[i].w), sizeof(Float));
569 for (i = 0; i<params.nb0; ++i) { // Tangential component
570 ofs.write(as_bytes(k.bonds_omega[i].x), sizeof(Float));
571 ofs.write(as_bytes(k.bonds_omega[i].y), sizeof(Float));
572 ofs.write(as_bytes(k.bonds_omega[i].z), sizeof(Float));
573 }
574
575 if (fluid == 1) {
576
577 ofs.write(as_bytes(cfd_solver), sizeof(int));
578
579 if (cfd_solver == 0) { // Navier Stokes flow
580
581 ofs.write(as_bytes(ns.mu), sizeof(Float));
582
583 int x, y, z;
584 for (z=0; z<ns.nz; z++) {
585 for (y=0; y<ns.ny; y++) {
586 for (x=0; x<ns.nx; x++) {
587 i = idx(x,y,z);
588
589 // Interpolated cell-center velocities
590 ofs.write(as_bytes(ns.v[i].x), sizeof(Float));
591 ofs.write(as_bytes(ns.v[i].y), sizeof(Float));
592 ofs.write(as_bytes(ns.v[i].z), sizeof(Float));
593
594 // Cell-face velocities
595 //ofs.write(as_bytes(ns.v_x[vidx(x,y,z)]), sizeof(Float));
596 //ofs.write(as_bytes(ns.v_y[vidx(x,y,z)]), sizeof(Float));
597 //ofs.write(as_bytes(ns.v_z[vidx(x,y,z)]), sizeof(Float));
598
599 ofs.write(as_bytes(ns.p[i]), sizeof(Float));
600 ofs.write(as_bytes(ns.phi[i]), sizeof(Float));
601 ofs.write(as_bytes(ns.dphi[i]), sizeof(Float));
602 }
603 }
604 }
605
606 ofs.write(as_bytes(ns.rho_f), sizeof(Float));
607 ofs.write(as_bytes(ns.p_mod_A), sizeof(Float));
608 ofs.write(as_bytes(ns.p_mod_f), sizeof(Float));
609 ofs.write(as_bytes(ns.p_mod_phi), sizeof(Float));
610
611 ofs.write(as_bytes(ns.bc_bot), sizeof(int));
612 ofs.write(as_bytes(ns.bc_top), sizeof(int));
613 ofs.write(as_bytes(ns.free_slip_bot), sizeof(int));
614 ofs.write(as_bytes(ns.free_slip_top), sizeof(int));
615 ofs.write(as_bytes(ns.bc_bot_flux), sizeof(Float));
616 ofs.write(as_bytes(ns.bc_top_flux), sizeof(Float));
617
618 for (z = 0; z<ns.nz; ++z)
619 for (y = 0; y<ns.ny; ++y)
620 for (x = 0; x<ns.nx; ++x)
621 ofs.write(as_bytes(ns.p_constant[idx(x,y,z)]),
622 sizeof(int));
623
624 ofs.write(as_bytes(ns.gamma), sizeof(Float));
625 ofs.write(as_bytes(ns.theta), sizeof(Float));
626 ofs.write(as_bytes(ns.beta), sizeof(Float));
627 ofs.write(as_bytes(ns.tolerance), sizeof(Float));
628 ofs.write(as_bytes(ns.maxiter), sizeof(unsigned int));
629 ofs.write(as_bytes(ns.ndem), sizeof(unsigned int));
630
631 ofs.write(as_bytes(ns.c_phi), sizeof(Float));
632 ofs.write(as_bytes(ns.c_v), sizeof(Float));
633 ofs.write(as_bytes(ns.dt_dem_fac), sizeof(Float));
634
635 for (i = 0; i<np; ++i) {
636 ofs.write(as_bytes(ns.f_d[i].x), sizeof(Float));
637 ofs.write(as_bytes(ns.f_d[i].y), sizeof(Float));
638 ofs.write(as_bytes(ns.f_d[i].z), sizeof(Float));
639 }
640 for (i = 0; i<np; ++i) {
641 ofs.write(as_bytes(ns.f_p[i].x), sizeof(Float));
642 ofs.write(as_bytes(ns.f_p[i].y), sizeof(Float));
643 ofs.write(as_bytes(ns.f_p[i].z), sizeof(Float));
644 }
645 for (i = 0; i<np; ++i) {
646 ofs.write(as_bytes(ns.f_v[i].x), sizeof(Float));
647 ofs.write(as_bytes(ns.f_v[i].y), sizeof(Float));
648 ofs.write(as_bytes(ns.f_v[i].z), sizeof(Float));
649 }
650 for (i = 0; i<np; ++i) {
651 ofs.write(as_bytes(ns.f_sum[i].x), sizeof(Float));
652 ofs.write(as_bytes(ns.f_sum[i].y), sizeof(Float));
653 ofs.write(as_bytes(ns.f_sum[i].z), sizeof(Float));
654 }
655
656 } else if (cfd_solver == 1) { // Darcy flow
657
658 ofs.write(as_bytes(darcy.mu), sizeof(Float));
659
660 int x, y, z;
661 for (z=0; z<darcy.nz; z++) {
662 for (y=0; y<darcy.ny; y++) {
663 for (x=0; x<darcy.nx; x++) {
664 i = d_idx(x,y,z);
665
666 // Interpolated cell-center velocities
667 ofs.write(as_bytes(darcy.v[i].x), sizeof(Float));
668 ofs.write(as_bytes(darcy.v[i].y), sizeof(Float));
669 ofs.write(as_bytes(darcy.v[i].z), sizeof(Float));
670
671 // Cell-face velocities
672 //ofs.write(as_bytes(darcy.v_x[vidx(x,y,z)]), sizeof(Float));
673 //ofs.write(as_bytes(darcy.v_y[vidx(x,y,z)]), sizeof(Float));
674 //ofs.write(as_bytes(darcy.v_z[vidx(x,y,z)]), sizeof(Float));
675
676 ofs.write(as_bytes(darcy.p[i]), sizeof(Float));
677 ofs.write(as_bytes(darcy.phi[i]), sizeof(Float));
678 ofs.write(as_bytes(darcy.dphi[i]), sizeof(Float));
679 }
680 }
681 }
682
683 ofs.write(as_bytes(darcy.rho_f), sizeof(Float));
684 ofs.write(as_bytes(darcy.p_mod_A), sizeof(Float));
685 ofs.write(as_bytes(darcy.p_mod_f), sizeof(Float));
686 ofs.write(as_bytes(darcy.p_mod_phi), sizeof(Float));
687
688 ofs.write(as_bytes(darcy.bc_xn), sizeof(int));
689 ofs.write(as_bytes(darcy.bc_xp), sizeof(int));
690 ofs.write(as_bytes(darcy.bc_yn), sizeof(int));
691 ofs.write(as_bytes(darcy.bc_yp), sizeof(int));
692 ofs.write(as_bytes(darcy.bc_bot), sizeof(int));
693 ofs.write(as_bytes(darcy.bc_top), sizeof(int));
694 ofs.write(as_bytes(darcy.free_slip_bot), sizeof(int));
695 ofs.write(as_bytes(darcy.free_slip_top), sizeof(int));
696 ofs.write(as_bytes(darcy.bc_bot_flux), sizeof(Float));
697 ofs.write(as_bytes(darcy.bc_top_flux), sizeof(Float));
698
699 for (z = 0; z<darcy.nz; ++z)
700 for (y = 0; y<darcy.ny; ++y)
701 for (x = 0; x<darcy.nx; ++x)
702 ofs.write(as_bytes(darcy.p_constant[d_idx(x,y,z)]),
703 sizeof(int));
704
705 ofs.write(as_bytes(darcy.tolerance), sizeof(Float));
706 ofs.write(as_bytes(darcy.maxiter), sizeof(unsigned int));
707 ofs.write(as_bytes(darcy.ndem), sizeof(unsigned int));
708 ofs.write(as_bytes(darcy.c_phi), sizeof(Float));
709
710 for (i = 0; i<np; ++i) {
711 ofs.write(as_bytes(darcy.f_p[i].x), sizeof(Float));
712 ofs.write(as_bytes(darcy.f_p[i].y), sizeof(Float));
713 ofs.write(as_bytes(darcy.f_p[i].z), sizeof(Float));
714 }
715 ofs.write(as_bytes(darcy.beta_f), sizeof(Float));
716 ofs.write(as_bytes(darcy.k_c), sizeof(Float));
717 }
718 }
719
720 for (i = 0; i<np; ++i)
721 ofs.write(as_bytes(k.color[i]), sizeof(int));
722
723 // Close file if it is still open
724 if (ofs.is_open())
725 ofs.close();
726
727 } else {
728 std::cerr << "Can't write output when in single precision mode.\n";
729 exit(1);
730 }
731 }
732
733 // Write image structure to PPM file
734 void DEM::writePPM(const char *target)
735 {
736 // Open output file
737 std::ofstream ofs(target);
738 if (!ofs) {
739 std::cerr << "Could not create output PPM file '"
740 << target << std::endl;
741 exit(1); // Return unsuccessful exit status
742 }
743
744 if (verbose == 1)
745 std::cout << " Saving image: " << target << std::endl;
746
747 // Write PPM header
748 ofs << "P6 " << width << " " << height << " 255\n";
749
750 // Write pixel array to ppm image file
751 for (unsigned int i=0; i<height*width; ++i)
752 ofs << img[i].r << img[i].g << img[i].b;
753
754 // Close file if it is still open
755 if (ofs.is_open())
756 ofs.close();
757 }
758
759 // Write write depth vs. porosity values to file
760 void DEM::writePorosities(
761 const char *target,
762 const int z_slices,
763 const Float *z_pos,
764 const Float *porosity)
765 {
766 // Open output file
767 std::ofstream ofs(target);
768 if (!ofs) {
769 std::cerr << "Could not create output porosity file '"
770 << target << std::endl;
771 exit(1); // Return unsuccessful exit status
772 }
773
774 if (verbose == 1)
775 std::cout << " Saving porosities: " << target << std::endl;
776
777 // Write pixel array to ppm image file
778 for (int i=0; i<z_slices; ++i)
779 ofs << z_pos[i] << '\t' << porosity[i] << '\n';
780
781 // Close file if it is still open
782 if (ofs.is_open())
783 ofs.close();
784 }
785
786
787
788 // vim: tabstop=8 expandtab shiftwidth=4 softtabstop=4