MFC
Exascale flow solver
Loading...
Searching...
No Matches
m_data_output.fpp.f90
Go to the documentation of this file.
1# 1 "/home/runner/work/MFC/MFC/src/post_process/m_data_output.fpp"
2!>
3!! @file
4!! @brief Contains module m_data_output
5
6!> @brief Writes post-processed grid and flow-variable data to Silo-HDF5 or binary database files
7module m_data_output
8
9 use m_derived_types
10 use m_global_parameters
11 use m_derived_variables
12 use m_mpi_proxy
13 use m_compile_specific
14 use m_helper
15 use m_variables_conversion
16 use m_eos
17 use m_constants, only: model_eqns_gamma_law, model_eqns_5eq, model_eqns_6eq, format_silo, format_binary, precision_single
18
19 implicit none
20
21 private; public :: s_initialize_data_output_module, s_define_output_region, s_open_formatted_database_file, &
22 & s_open_intf_data_file, s_open_energy_data_file, s_write_grid_to_formatted_database_file, &
23 & s_write_variable_to_formatted_database_file, s_write_lag_bubbles_results_to_text, &
24 & s_write_lag_bubbles_to_formatted_database_file, s_write_ib_state_files, s_write_intf_data_file, &
25 & s_write_energy_data_file, s_write_ib_bodies_to_formatted_database_file, s_close_formatted_database_file, &
26 & s_close_intf_data_file, s_close_energy_data_file, s_finalize_data_output_module, out
27
28 ! Include Silo-HDF5 interface library
29 include 'silo_f9x.inc'
30
31 !> Output workspace: flow variable buffers, VisIt extents/offsets, directory paths, file handles, and variable count.
32 type(output_context) :: out
33
34 ! Generic error flag for Silo-HDF5 and Binary I/O operations
35 integer, private :: err
36
37contains
38
39 !> Allocate storage arrays, configure output directories, and count flow variables for formatted database output.
40 impure subroutine s_initialize_data_output_module()
41
42 character(LEN=len_trim(case_dir) + 2*name_len) :: file_loc
43 logical :: dir_check
44 integer :: i
45
46 allocate (out%q_sf(-offset_x%beg:m + offset_x%end,-offset_y%beg:n + offset_y%end,-offset_z%beg:p + offset_z%end))
47 if (grid_geometry == 3) then
48 allocate (out%cyl_q_sf(-offset_y%beg:n + offset_y%end,-offset_z%beg:p + offset_z%end,-offset_x%beg:m + offset_x%end))
49 end if
50
51 if (precision == precision_single) then
52 allocate (out%q_sf_s(-offset_x%beg:m + offset_x%end,-offset_y%beg:n + offset_y%end,-offset_z%beg:p + offset_z%end))
53 if (grid_geometry == 3) then
54 allocate (out%cyl_q_sf_s(-offset_y%beg:n + offset_y%end,-offset_z%beg:p + offset_z%end, &
55 & -offset_x%beg:m + offset_x%end))
56 end if
57 end if
58
59 if (n == 0) then
60 allocate (out%q_root_sf(0:m_root,0:0,0:0))
61 if (precision == precision_single) then
62 allocate (out%q_root_sf_s(0:m_root,0:0,0:0))
63 end if
64 end if
65
66 ! Allocating the spatial and data extents and also the variables for the offsets and the one bookkeeping the number of
67 ! cell-boundaries in each active coordinate direction. Note that all these variables are only needed by the Silo-HDF5 format
68 ! for multidimensional data.
69 if (format == format_silo) then
70 allocate (out%data_extents(1:2,0:num_procs - 1))
71
72 if (p > 0) then
73 allocate (out%spatial_extents(1:6,0:num_procs - 1))
74 allocate (out%lo_offset(1:3))
75 allocate (out%hi_offset(1:3))
76 allocate (out%dims(1:3))
77 else if (n > 0) then
78 allocate (out%spatial_extents(1:4,0:num_procs - 1))
79 allocate (out%lo_offset(1:2))
80 allocate (out%hi_offset(1:2))
81 allocate (out%dims(1:2))
82 else
83 allocate (out%spatial_extents(1:2,0:num_procs - 1))
84 allocate (out%lo_offset(1:1))
85 allocate (out%hi_offset(1:1))
86 allocate (out%dims(1:1))
87 end if
88 end if
89
90 ! The size of the ghost zone layer in each of the active coordinate directions was set in the module m_mpi_proxy.f90. The
91 ! results are now transferred to the local variables of this module when they are required by the Silo-HDF5 format, for
92 ! multidimensional data sets. With the same, latter, requirements, the variables bookkeeping the number of cell-boundaries
93 ! in each active coordinate direction are also set here.
94 if (format == format_silo) then
95 if (p > 0) then
96 if (grid_geometry == 3) then
97 out%lo_offset(:) = (/offset_y%beg, offset_z%beg, offset_x%beg/)
98 out%hi_offset(:) = (/offset_y%end, offset_z%end, offset_x%end/)
99 else
100 out%lo_offset(:) = (/offset_x%beg, offset_y%beg, offset_z%beg/)
101 out%hi_offset(:) = (/offset_x%end, offset_y%end, offset_z%end/)
102 end if
103
104 if (grid_geometry == 3) then
105 out%dims(:) = (/n + offset_y%beg + offset_y%end + 2, p + offset_z%beg + offset_z%end + 2, &
106 & m + offset_x%beg + offset_x%end + 2/)
107 else
108 out%dims(:) = (/m + offset_x%beg + offset_x%end + 2, n + offset_y%beg + offset_y%end + 2, &
109 & p + offset_z%beg + offset_z%end + 2/)
110 end if
111 else if (n > 0) then
112 out%lo_offset(:) = (/offset_x%beg, offset_y%beg/)
113 out%hi_offset(:) = (/offset_x%end, offset_y%end/)
114
115 out%dims(:) = (/m + offset_x%beg + offset_x%end + 2, n + offset_y%beg + offset_y%end + 2/)
116 else
117 out%lo_offset(:) = (/offset_x%beg/)
118 out%hi_offset(:) = (/offset_x%end/)
119 out%dims(:) = (/m + offset_x%beg + offset_x%end + 2/)
120 end if
121 end if
122
123 if (format == format_silo) then
124 out%dbdir = trim(case_dir) // '/silo_hdf5'
125
126 write (out%proc_rank_dir, '(A,I0)') '/p', proc_rank
127
128 out%proc_rank_dir = trim(out%dbdir) // trim(out%proc_rank_dir)
129
130 file_loc = trim(out%proc_rank_dir) // '/.'
131
132 call my_inquire(file_loc, dir_check)
133 if (dir_check .neqv. .true.) then
134 call s_create_directory(trim(out%proc_rank_dir))
135 end if
136
137 if (proc_rank == 0) then
138 out%rootdir = trim(out%dbdir) // '/root'
139
140 file_loc = trim(out%rootdir) // '/.'
141
142 call my_inquire(file_loc, dir_check)
143 if (dir_check .neqv. .true.) then
144 call s_create_directory(trim(out%rootdir))
145 end if
146 end if
147 else
148 out%dbdir = trim(case_dir) // '/binary'
149
150 write (out%proc_rank_dir, '(A,I0)') '/p', proc_rank
151
152 out%proc_rank_dir = trim(out%dbdir) // trim(out%proc_rank_dir)
153
154 file_loc = trim(out%proc_rank_dir) // '/.'
155
156 call my_inquire(file_loc, dir_check)
157
158 if (dir_check .neqv. .true.) then
159 call s_create_directory(trim(out%proc_rank_dir))
160 end if
161
162 if (n == 0 .and. proc_rank == 0) then
163 out%rootdir = trim(out%dbdir) // '/root'
164
165 file_loc = trim(out%rootdir) // '/.'
166
167 call my_inquire(file_loc, dir_check)
168
169 if (dir_check .neqv. .true.) then
170 call s_create_directory(trim(out%rootdir))
171 end if
172 end if
173 end if
174
175 if (bubbles_lagrange) then ! Lagrangian solver
176 if (lag_txt_wrt) then
177 out%dbdir = trim(case_dir) // '/lag_bubbles_post_process'
178 file_loc = trim(out%dbdir) // '/.'
179 call my_inquire(file_loc, dir_check)
180
181 if (dir_check .neqv. .true.) then
182 call s_create_directory(trim(out%dbdir))
183 end if
184 end if
185 end if
186
187 ! Contrary to the Silo-HDF5 database format, handles of the Binary database master/root and slave/local process files are
188 ! perfectly static throughout post-process. Hence, they are set here so that they do not have to be repetitively computed in
189 ! later procedures.
190 if (format == format_binary) then
191 if (n == 0 .and. proc_rank == 0) out%dbroot = 2
192 out%dbfile = 1
193 end if
194
195 if (format == format_binary) then
196 out%dbvars = 0
197
198 if ((model_eqns == model_eqns_5eq) .or. (model_eqns == model_eqns_6eq)) then
199 do i = 1, num_fluids
200 if (alpha_rho_wrt(i) .or. (cons_vars_wrt .or. prim_vars_wrt)) then
201 out%dbvars = out%dbvars + 1
202 end if
203 end do
204 end if
205
206 if ((rho_wrt .or. (model_eqns == model_eqns_gamma_law .and. (cons_vars_wrt .or. prim_vars_wrt))) &
207 & .and. (.not. relativity)) then
208 out%dbvars = out%dbvars + 1
209 end if
210
211 if (relativity .and. (rho_wrt .or. prim_vars_wrt)) out%dbvars = out%dbvars + 1
212 if (relativity .and. (rho_wrt .or. cons_vars_wrt)) out%dbvars = out%dbvars + 1
213
214 do i = 1, eqn_idx%E - eqn_idx%mom%beg
215 if (mom_wrt(i) .or. cons_vars_wrt) out%dbvars = out%dbvars + 1
216 end do
217
218 do i = 1, eqn_idx%E - eqn_idx%mom%beg
219 if (vel_wrt(i) .or. prim_vars_wrt) out%dbvars = out%dbvars + 1
220 end do
221
222 do i = 1, eqn_idx%E - eqn_idx%mom%beg
223 if (flux_wrt(i)) out%dbvars = out%dbvars + 1
224 end do
225
226 if (e_wrt .or. cons_vars_wrt) out%dbvars = out%dbvars + 1
227 if (pres_wrt .or. prim_vars_wrt) out%dbvars = out%dbvars + 1
228 if (hypoelasticity) out%dbvars = out%dbvars + (num_dims*(num_dims + 1))/2
229 if (cont_damage) out%dbvars = out%dbvars + 1
230 if (hyper_cleaning) out%dbvars = out%dbvars + 1
231
232 if (mhd) then
233 if (n == 0) then
234 out%dbvars = out%dbvars + 2
235 else
236 out%dbvars = out%dbvars + 3
237 end if
238 end if
239
240 if ((model_eqns == model_eqns_5eq) .or. (model_eqns == model_eqns_6eq)) then
241 do i = 1, num_fluids - 1
242 if (alpha_wrt(i) .or. (cons_vars_wrt .or. prim_vars_wrt)) then
243 out%dbvars = out%dbvars + 1
244 end if
245 end do
246
247 if (alpha_wrt(num_fluids) .or. (cons_vars_wrt .or. prim_vars_wrt)) then
248 out%dbvars = out%dbvars + 1
249 end if
250 end if
251
252 if (gamma_wrt .or. (model_eqns == model_eqns_gamma_law .and. (cons_vars_wrt .or. prim_vars_wrt))) then
253 out%dbvars = out%dbvars + 1
254 end if
255
256 if (heat_ratio_wrt) out%dbvars = out%dbvars + 1
257
258 if (pi_inf_wrt .or. (model_eqns == model_eqns_gamma_law .and. (cons_vars_wrt .or. prim_vars_wrt))) then
259 out%dbvars = out%dbvars + 1
260 end if
261
262 if (pres_inf_wrt) out%dbvars = out%dbvars + 1
263 if (c_wrt) out%dbvars = out%dbvars + 1
264
265 if (p > 0) then
266 do i = 1, num_vels
267 if (omega_wrt(i)) out%dbvars = out%dbvars + 1
268 end do
269 else if (n > 0) then
270 do i = 1, num_vels
271 if (omega_wrt(i)) out%dbvars = out%dbvars + 1
272 end do
273 end if
274
275 if (schlieren_wrt) out%dbvars = out%dbvars + 1
276 end if
277
278 end subroutine s_initialize_data_output_module
279
280 !> Compute the cell-index bounds for the user-specified partial output domain in each coordinate direction.
281 impure subroutine s_define_output_region
282
283 integer :: i
284 integer :: lower_bound, upper_bound
285
286# 286 "/home/runner/work/MFC/MFC/src/post_process/m_data_output.fpp"
287 if (m == 0) return ! Early return for y or z if simulation is 1D or 2D
288
289 lower_bound = -offset_x%beg
290 upper_bound = m + offset_x%end
291
292 do i = lower_bound, upper_bound
293 if (x_cc(i) > x_output%beg) then
294 x_output_idx%beg = i + offset_x%beg
295 exit
296 end if
297 end do
298
299 do i = upper_bound, lower_bound, -1
300 if (x_cc(i) < x_output%end) then
301 x_output_idx%end = i + offset_x%beg
302 exit
303 end if
304 end do
305
306 ! If no grid points are within the output region
307 if ((x_cc(lower_bound) > x_output%end) .or. (x_cc(upper_bound) < x_output%beg)) then
308 x_output_idx%beg = 0
309 x_output_idx%end = 0
310 end if
311# 286 "/home/runner/work/MFC/MFC/src/post_process/m_data_output.fpp"
312 if (n == 0) return ! Early return for y or z if simulation is 1D or 2D
313
314 lower_bound = -offset_y%beg
315 upper_bound = n + offset_y%end
316
317 do i = lower_bound, upper_bound
318 if (y_cc(i) > y_output%beg) then
319 y_output_idx%beg = i + offset_y%beg
320 exit
321 end if
322 end do
323
324 do i = upper_bound, lower_bound, -1
325 if (y_cc(i) < y_output%end) then
326 y_output_idx%end = i + offset_y%beg
327 exit
328 end if
329 end do
330
331 ! If no grid points are within the output region
332 if ((y_cc(lower_bound) > y_output%end) .or. (y_cc(upper_bound) < y_output%beg)) then
333 y_output_idx%beg = 0
334 y_output_idx%end = 0
335 end if
336# 286 "/home/runner/work/MFC/MFC/src/post_process/m_data_output.fpp"
337 if (p == 0) return ! Early return for y or z if simulation is 1D or 2D
338
339 lower_bound = -offset_z%beg
340 upper_bound = p + offset_z%end
341
342 do i = lower_bound, upper_bound
343 if (z_cc(i) > z_output%beg) then
344 z_output_idx%beg = i + offset_z%beg
345 exit
346 end if
347 end do
348
349 do i = upper_bound, lower_bound, -1
350 if (z_cc(i) < z_output%end) then
351 z_output_idx%end = i + offset_z%beg
352 exit
353 end if
354 end do
355
356 ! If no grid points are within the output region
357 if ((z_cc(lower_bound) > z_output%end) .or. (z_cc(upper_bound) < z_output%beg)) then
358 z_output_idx%beg = 0
359 z_output_idx%end = 0
360 end if
361# 311 "/home/runner/work/MFC/MFC/src/post_process/m_data_output.fpp"
362
363 end subroutine s_define_output_region
364
365 !> Open (or create) the Silo-HDF5 or Binary formatted database slave and master files for a given time step.
366 impure subroutine s_open_formatted_database_file(t_step)
367
368 integer, intent(in) :: t_step
369 character(LEN=len_trim(case_dir) + 3*name_len) :: file_loc
370 integer :: ierr
371
372 if (format == format_silo) then
373 write (file_loc, '(A,I0,A)') '/', t_step, '.silo'
374 file_loc = trim(out%proc_rank_dir) // trim(file_loc)
375
376 ierr = dbcreate(trim(file_loc), len_trim(file_loc), db_clobber, db_local, 'MFC v3.0', 8, db_hdf5, out%dbfile)
377
378 if (out%dbfile == -1) then
379 call s_mpi_abort('Unable to create Silo-HDF5 database ' // 'slave file ' // trim(file_loc) // '. ' // 'Exiting.')
380 end if
381
382 if (proc_rank == 0) then
383 write (file_loc, '(A,I0,A)') '/collection_', t_step, '.silo'
384 file_loc = trim(out%rootdir) // trim(file_loc)
385
386 ierr = dbcreate(trim(file_loc), len_trim(file_loc), db_clobber, db_local, 'MFC v3.0', 8, db_hdf5, out%dbroot)
387
388 if (out%dbroot == -1) then
389 call s_mpi_abort('Unable to create Silo-HDF5 database ' // 'master file ' // trim(file_loc) // '. ' &
390 & // 'Exiting.')
391 end if
392 end if
393 else
394 write (file_loc, '(A,I0,A)') '/', t_step, '.dat'
395 file_loc = trim(out%proc_rank_dir) // trim(file_loc)
396
397 open (out%dbfile, iostat=err, file=trim(file_loc), form='unformatted', status='replace')
398
399 if (err /= 0) then
400 call s_mpi_abort('Unable to create Binary database slave ' // 'file ' // trim(file_loc) // '. Exiting.')
401 end if
402
403 if (output_partial_domain) then
404 write (out%dbfile) x_output_idx%end - x_output_idx%beg, y_output_idx%end - y_output_idx%beg, &
405 & z_output_idx%end - z_output_idx%beg, out%dbvars
406 else
407 write (out%dbfile) m, n, p, out%dbvars
408 end if
409
410 if (n == 0 .and. proc_rank == 0) then
411 write (file_loc, '(A,I0,A)') '/', t_step, '.dat'
412 file_loc = trim(out%rootdir) // trim(file_loc)
413
414 open (out%dbroot, iostat=err, file=trim(file_loc), form='unformatted', status='replace')
415
416 if (err /= 0) then
417 call s_mpi_abort('Unable to create Binary database ' // 'master file ' // trim(file_loc) // '. Exiting.')
418 end if
419
420 if (output_partial_domain) then
421 write (out%dbroot) x_output_idx%end - x_output_idx%beg, 0, 0, out%dbvars
422 else
423 write (out%dbroot) m_root, 0, 0, out%dbvars
424 end if
425 end if
426 end if
427
428 end subroutine s_open_formatted_database_file
429
430 !> Open the interface data file for appending extracted interface coordinates.
431 impure subroutine s_open_intf_data_file()
432
433 character(LEN=path_len + 3*name_len) :: file_path
434
435 write (file_path, '(A)') '/intf_data.dat'
436 file_path = trim(case_dir) // trim(file_path)
437
438 open (211, file=trim(file_path), form='formatted', position='append', status='unknown')
439
440 end subroutine s_open_intf_data_file
441
442 !> Open the energy data file for appending volume-integrated energy budget quantities.
443 impure subroutine s_open_energy_data_file()
444
445 character(LEN=path_len + 3*name_len) :: file_path
446
447 write (file_path, '(A)') '/eng_data.dat'
448 file_path = trim(case_dir) // trim(file_path)
449
450 open (251, file=trim(file_path), form='formatted', position='append', status='unknown')
451
452 end subroutine s_open_energy_data_file
453
454 !> Write the computational grid (cell-boundary coordinates) to the formatted database slave and master files.
455 impure subroutine s_write_grid_to_formatted_database_file(t_step)
456
457 integer, intent(in) :: t_step
458
459 ! NAG compiler requires these to be statically sized
460 character(LEN=4*name_len), dimension(num_procs) :: meshnames
461 integer, dimension(num_procs) :: meshtypes
462 integer :: i
463 integer :: ierr
464 integer :: extents_size
465
466 if (format == format_silo) then
467 ! For multidimensional data sets, the spatial extents of all of the grid(s) handled by the local processor(s) are
468 ! recorded so that they may be written, by root processor, to the formatted database master file.
469 if (num_procs > 1) then
470 call s_mpi_gather_spatial_extents(out%spatial_extents)
471 else if (p > 0) then
472 if (grid_geometry == 3) then
473 out%spatial_extents(:,0) = (/minval(y_cb), minval(z_cb), minval(x_cb), maxval(y_cb), maxval(z_cb), &
474 & maxval(x_cb)/)
475 else
476 out%spatial_extents(:,0) = (/minval(x_cb), minval(y_cb), minval(z_cb), maxval(x_cb), maxval(y_cb), &
477 & maxval(z_cb)/)
478 end if
479 else if (n > 0) then
480 out%spatial_extents(:,0) = (/minval(x_cb), minval(y_cb), maxval(x_cb), maxval(y_cb)/)
481 else
482 out%spatial_extents(:,0) = (/minval(x_cb), maxval(x_cb)/)
483 end if
484
485 ! Next, the root processor proceeds to record all of the spatial extents in the formatted database master file. In
486 ! addition, it also records a sub-domain connectivity map so that the entire grid may be reassembled by looking at the
487 ! master file.
488 if (proc_rank == 0) then
489 do i = 1, num_procs
490 write (meshnames(i), '(A,I0,A,I0,A)') '../p', i - 1, '/', t_step, '.silo:rectilinear_grid'
491 end do
492
493 meshtypes = db_quad_rect
494
495 err = dbset2dstrlen(len(meshnames(1)))
496 err = dbmkoptlist(2, out%optlist)
497 extents_size = size(out%spatial_extents, 1)
498 err = dbaddiopt(out%optlist, dbopt_extents_size, extents_size)
499 err = dbadddopt(out%optlist, dbopt_extents, out%spatial_extents)
500 err = dbputmmesh(out%dbroot, 'rectilinear_grid', 16, num_procs, meshnames, len_trim(meshnames), meshtypes, &
501 & out%optlist, ierr)
502 err = dbfreeoptlist(out%optlist)
503 end if
504
505 ! Finally, the local quadrilateral mesh, either 2D or 3D, along with its offsets that indicate the presence and size of
506 ! ghost zone layer(s), are put in the formatted database slave file.
507
508 if (p > 0) then
509 err = dbmkoptlist(2, out%optlist)
510 err = dbaddiaopt(out%optlist, dbopt_lo_offset, size(out%lo_offset), out%lo_offset)
511 err = dbaddiaopt(out%optlist, dbopt_hi_offset, size(out%hi_offset), out%hi_offset)
512 if (grid_geometry == 3) then
513 err = dbputqm(out%dbfile, 'rectilinear_grid', 16, 'x', 1, 'y', 1, 'z', 1, y_cb, z_cb, x_cb, out%dims, 3, &
514 & db_double, db_collinear, out%optlist, ierr)
515 else
516 err = dbputqm(out%dbfile, 'rectilinear_grid', 16, 'x', 1, 'y', 1, 'z', 1, x_cb, y_cb, z_cb, out%dims, 3, &
517 & db_double, db_collinear, out%optlist, ierr)
518 end if
519 err = dbfreeoptlist(out%optlist)
520 else if (n > 0) then
521 err = dbmkoptlist(2, out%optlist)
522 err = dbaddiaopt(out%optlist, dbopt_lo_offset, size(out%lo_offset), out%lo_offset)
523 err = dbaddiaopt(out%optlist, dbopt_hi_offset, size(out%hi_offset), out%hi_offset)
524 err = dbputqm(out%dbfile, 'rectilinear_grid', 16, 'x', 1, 'y', 1, 'z', 1, x_cb, y_cb, db_f77null, out%dims, 2, &
525 & db_double, db_collinear, out%optlist, ierr)
526 err = dbfreeoptlist(out%optlist)
527 else
528 err = dbmkoptlist(2, out%optlist)
529 err = dbaddiaopt(out%optlist, dbopt_lo_offset, size(out%lo_offset), out%lo_offset)
530 err = dbaddiaopt(out%optlist, dbopt_hi_offset, size(out%hi_offset), out%hi_offset)
531 err = dbputqm(out%dbfile, 'rectilinear_grid', 16, 'x', 1, 'y', 1, 'z', 1, x_cb, db_f77null, db_f77null, out%dims, &
532 & 1, db_double, db_collinear, out%optlist, ierr)
533 err = dbfreeoptlist(out%optlist)
534 end if
535 else if (format == format_binary) then
536 ! Multidimensional local grid data is written to the formatted database slave file. Recall that no master file to
537 ! maintained in multidimensions.
538 if (p > 0) then
539 if (precision == precision_single) then
540 write (out%dbfile) real(x_cb, sp), real(y_cb, sp), real(z_cb, sp)
541 else
542 if (output_partial_domain) then
543 write (out%dbfile) x_cb(x_output_idx%beg - 1:x_output_idx%end), &
544 & y_cb(y_output_idx%beg - 1:y_output_idx%end), z_cb(z_output_idx%beg - 1:z_output_idx%end)
545 else
546 write (out%dbfile) x_cb, y_cb, z_cb
547 end if
548 end if
549 else if (n > 0) then
550 if (precision == precision_single) then
551 write (out%dbfile) real(x_cb, sp), real(y_cb, sp)
552 else
553 if (output_partial_domain) then
554 write (out%dbfile) x_cb(x_output_idx%beg - 1:x_output_idx%end), y_cb(y_output_idx%beg - 1:y_output_idx%end)
555 else
556 write (out%dbfile) x_cb, y_cb
557 end if
558 end if
559
560 ! One-dimensional local grid data is written to the formatted database slave file. In addition, the local grid data
561 ! is put together by the root process and written to the master file.
562 else
563 if (precision == precision_single) then
564 write (out%dbfile) real(x_cb, sp)
565 else
566 if (output_partial_domain) then
567 write (out%dbfile) x_cb(x_output_idx%beg - 1:x_output_idx%end)
568 else
569 write (out%dbfile) x_cb
570 end if
571 end if
572
573 if (num_procs > 1) then
574 call s_mpi_defragment_1d_grid_variable()
575 else
576 x_root_cb(:) = x_cb(:)
577 end if
578
579 if (proc_rank == 0) then
580 if (precision == precision_single) then
581 write (out%dbroot) real(x_root_cb, wp)
582 else
583 if (output_partial_domain) then
584 write (out%dbroot) x_root_cb(x_output_idx%beg - 1:x_output_idx%end)
585 else
586 write (out%dbroot) x_root_cb
587 end if
588 end if
589 end if
590 end if
591 end if
592
593 end subroutine s_write_grid_to_formatted_database_file
594
595 !> Write a single flow variable field to the formatted database slave and master files for a given time step.
596 impure subroutine s_write_variable_to_formatted_database_file(varname, t_step)
597
598 character(LEN=*), intent(in) :: varname
599 integer, intent(in) :: t_step
600
601 ! NAG compiler requires these to be statically sized
602 character(LEN=4*name_len), dimension(num_procs) :: varnames
603 integer, dimension(num_procs) :: vartypes
604 integer :: i, j, k
605 integer :: ierr
606 integer :: extents_size
607
608 if (format == format_silo) then
609 ! Determining the extents of the flow variable on each local process and gathering all this information on root process
610 if (num_procs > 1) then
611 call s_mpi_gather_data_extents(out%q_sf, out%data_extents)
612 else
613 out%data_extents(:,0) = (/minval(out%q_sf), maxval(out%q_sf)/)
614 end if
615
616 if (proc_rank == 0) then
617 do i = 1, num_procs
618 write (varnames(i), '(A,I0,A,I0,A)') '../p', i - 1, '/', t_step, '.silo:' // trim(varname)
619 end do
620
621 vartypes = db_quadvar
622
623 err = dbset2dstrlen(len(varnames(1)))
624 err = dbmkoptlist(2, out%optlist)
625 extents_size = size(out%data_extents, 1)
626 err = dbaddiopt(out%optlist, dbopt_extents_size, extents_size)
627 err = dbadddopt(out%optlist, dbopt_extents, out%data_extents)
628 err = dbputmvar(out%dbroot, trim(varname), len_trim(varname), num_procs, varnames, len_trim(varnames), vartypes, &
629 & out%optlist, ierr)
630 err = dbfreeoptlist(out%optlist)
631 end if
632
633 if (wp == dp) then
634 if (precision == precision_single) then
635 do i = -offset_x%beg, m + offset_x%end
636 do j = -offset_y%beg, n + offset_y%end
637 do k = -offset_z%beg, p + offset_z%end
638 out%q_sf_s(i, j, k) = real(out%q_sf(i, j, k), sp)
639 end do
640 end do
641 end do
642 if (grid_geometry == 3) then
643 do i = -offset_x%beg, m + offset_x%end
644 do j = -offset_y%beg, n + offset_y%end
645 do k = -offset_z%beg, p + offset_z%end
646 out%cyl_q_sf_s(j, k, i) = out%q_sf_s(i, j, k)
647 end do
648 end do
649 end do
650 end if
651 else
652 if (grid_geometry == 3) then
653 do i = -offset_x%beg, m + offset_x%end
654 do j = -offset_y%beg, n + offset_y%end
655 do k = -offset_z%beg, p + offset_z%end
656 out%cyl_q_sf(j, k, i) = out%q_sf(i, j, k)
657 end do
658 end do
659 end do
660 end if
661 end if
662 else if (wp == sp) then
663 do i = -offset_x%beg, m + offset_x%end
664 do j = -offset_y%beg, n + offset_y%end
665 do k = -offset_z%beg, p + offset_z%end
666 out%q_sf_s(i, j, k) = out%q_sf(i, j, k)
667 end do
668 end do
669 end do
670 if (grid_geometry == 3) then
671 do i = -offset_x%beg, m + offset_x%end
672 do j = -offset_y%beg, n + offset_y%end
673 do k = -offset_z%beg, p + offset_z%end
674 out%cyl_q_sf_s(j, k, i) = out%q_sf_s(i, j, k)
675 end do
676 end do
677 end do
678 end if
679 end if
680
681# 631 "/home/runner/work/MFC/MFC/src/post_process/m_data_output.fpp"
682 if (precision == 1) then
683 if (p > 0) then
684 if (grid_geometry == 3) then
685 err = dbputqv1(out%dbfile, trim(varname), len_trim(varname), 'rectilinear_grid', 16, &
686 & out%cyl_q_sf_s, out%dims - 1, 3, db_f77null, 0, db_float, db_zonecent, &
687 & db_f77null, ierr)
688 else
689 err = dbputqv1(out%dbfile, trim(varname), len_trim(varname), 'rectilinear_grid', 16, out%q_sf_s, &
690 & out%dims - 1, 3, db_f77null, 0, db_float, db_zonecent, db_f77null, ierr)
691 end if
692 else if (n > 0) then
693 err = dbputqv1(out%dbfile, trim(varname), len_trim(varname), 'rectilinear_grid', 16, out%q_sf_s, &
694 & out%dims - 1, 2, db_f77null, 0, db_float, db_zonecent, db_f77null, ierr)
695 else
696 err = dbputqv1(out%dbfile, trim(varname), len_trim(varname), 'rectilinear_grid', 16, out%q_sf_s, &
697 & out%dims - 1, 1, db_f77null, 0, db_float, db_zonecent, db_f77null, ierr)
698 end if
699 end if
700# 631 "/home/runner/work/MFC/MFC/src/post_process/m_data_output.fpp"
701 if (precision == 2) then
702 if (p > 0) then
703 if (grid_geometry == 3) then
704 err = dbputqv1(out%dbfile, trim(varname), len_trim(varname), 'rectilinear_grid', 16, &
705 & out%cyl_q_sf, out%dims - 1, 3, db_f77null, 0, db_double, db_zonecent, &
706 & db_f77null, ierr)
707 else
708 err = dbputqv1(out%dbfile, trim(varname), len_trim(varname), 'rectilinear_grid', 16, out%q_sf, &
709 & out%dims - 1, 3, db_f77null, 0, db_double, db_zonecent, db_f77null, ierr)
710 end if
711 else if (n > 0) then
712 err = dbputqv1(out%dbfile, trim(varname), len_trim(varname), 'rectilinear_grid', 16, out%q_sf, &
713 & out%dims - 1, 2, db_f77null, 0, db_double, db_zonecent, db_f77null, ierr)
714 else
715 err = dbputqv1(out%dbfile, trim(varname), len_trim(varname), 'rectilinear_grid', 16, out%q_sf, &
716 & out%dims - 1, 1, db_f77null, 0, db_double, db_zonecent, db_f77null, ierr)
717 end if
718 end if
719# 650 "/home/runner/work/MFC/MFC/src/post_process/m_data_output.fpp"
720 else
721 ! Writing the name of the flow variable and its data, associated with the local processor, to the formatted database
722 ! slave file
723 if (precision == precision_single) then
724 write (out%dbfile) varname, real(out%q_sf, wp)
725 else
726 write (out%dbfile) varname, out%q_sf
727 end if
728
729 ! In 1D, the root process also takes care of gathering the flow variable data from all of the local processor(s) and
730 ! writes it to the formatted database master file.
731 if (n == 0) then
732 if (num_procs > 1) then
733 call s_mpi_defragment_1d_flow_variable(out%q_sf, out%q_root_sf)
734 else
735 out%q_root_sf(:,:,:) = out%q_sf(:,:,:)
736 end if
737
738 if (proc_rank == 0) then
739 if (precision == precision_single) then
740 write (out%dbroot) varname, real(out%q_root_sf, wp)
741 else
742 write (out%dbroot) varname, out%q_root_sf
743 end if
744 end if
745 end if
746 end if
747
748 end subroutine s_write_variable_to_formatted_database_file
749
750 !> Write the post-processed results in the folder 'lag_bubbles_data'
751 impure subroutine s_write_lag_bubbles_results_to_text(t_step)
752
753 integer, intent(in) :: t_step
754 character(len=len_trim(case_dir) + 3*name_len) :: file_loc
755 integer :: id
756
757#ifdef MFC_MPI
758 real(wp), dimension(20) :: inputvals
759 real(wp) :: time_real
760 integer, dimension(MPI_STATUS_SIZE) :: status
761 integer(KIND=MPI_OFFSET_KIND) :: disp
762 integer :: view
763 logical :: file_exist
764 integer, dimension(2) :: gsizes, lsizes, start_idx_part
765 integer :: ifile
766 integer :: ierr
767 real(wp) :: file_time, file_dt
768 integer :: file_num_procs, file_tot_part
769 integer :: i
770 integer, dimension(:), allocatable :: proc_bubble_counts
771 real(wp), dimension(1:1,1:lag_io_vars) :: lag_io_null
772
773 lag_io_null = 0._wp
774
775 ! Construct file path
776 write (file_loc, '(A,I0,A)') 'lag_bubbles_', t_step, '.dat'
777 file_loc = trim(case_dir) // '/restart_data' // trim(mpiiofs) // trim(file_loc)
778
779 ! Check if file exists
780 inquire (file=trim(file_loc), exist=file_exist)
781 if (.not. file_exist) then
782 call s_mpi_abort('Restart file ' // trim(file_loc) // ' does not exist!')
783 end if
784
785 if (.not. parallel_io) return
786
787 if (proc_rank == 0) then
788 call mpi_file_open(mpi_comm_self, file_loc, mpi_mode_rdonly, mpi_info_int, ifile, ierr)
789
790 call mpi_file_read(ifile, file_tot_part, 1, mpi_integer, status, ierr)
791 call mpi_file_read(ifile, file_time, 1, mpi_p, status, ierr)
792 call mpi_file_read(ifile, file_dt, 1, mpi_p, status, ierr)
793 call mpi_file_read(ifile, file_num_procs, 1, mpi_integer, status, ierr)
794
795 call mpi_file_close(ifile, ierr)
796 end if
797
798 call mpi_bcast(file_tot_part, 1, mpi_integer, 0, mpi_comm_world, ierr)
799 call mpi_bcast(file_time, 1, mpi_p, 0, mpi_comm_world, ierr)
800 call mpi_bcast(file_dt, 1, mpi_p, 0, mpi_comm_world, ierr)
801 call mpi_bcast(file_num_procs, 1, mpi_integer, 0, mpi_comm_world, ierr)
802 time_real = file_time
803
804 allocate (proc_bubble_counts(file_num_procs))
805
806 if (proc_rank == 0) then
807 call mpi_file_open(mpi_comm_self, file_loc, mpi_mode_rdonly, mpi_info_int, ifile, ierr)
808
809 ! Skip to processor counts position
810 disp = int(sizeof(file_tot_part) + 2*sizeof(file_time) + sizeof(file_num_procs), mpi_offset_kind)
811 call mpi_file_seek(ifile, disp, mpi_seek_set, ierr)
812 call mpi_file_read(ifile, proc_bubble_counts, file_num_procs, mpi_integer, status, ierr)
813
814 call mpi_file_close(ifile, ierr)
815 end if
816
817 call mpi_bcast(proc_bubble_counts, file_num_procs, mpi_integer, 0, mpi_comm_world, ierr)
818
819 if (file_tot_part > 0) then
820 gsizes(1) = file_tot_part
821 gsizes(2) = lag_io_vars
822 lsizes(1) = file_tot_part
823 lsizes(2) = lag_io_vars
824 start_idx_part(1) = 0
825 start_idx_part(2) = 0
826
827 call mpi_type_create_subarray(2, gsizes, lsizes, start_idx_part, mpi_order_fortran, mpi_p, view, ierr)
828 call mpi_type_commit(view, ierr)
829
830 call mpi_file_open(mpi_comm_world, file_loc, mpi_mode_rdonly, mpi_info_int, ifile, ierr)
831
832 disp = int(sizeof(file_tot_part) + 2*sizeof(file_time) + sizeof(file_num_procs) &
833 & + file_num_procs*sizeof(proc_bubble_counts(1)), mpi_offset_kind)
834 call mpi_file_set_view(ifile, disp, mpi_p, view, 'native', mpi_info_null, ierr)
835
836 allocate (mpi_io_data_lg_bubbles(file_tot_part,1:lag_io_vars))
837
838 call mpi_file_read_all(ifile, mpi_io_data_lg_bubbles, lag_io_vars*file_tot_part, mpi_p, status, ierr)
839
840 write (file_loc, '(A,I0,A)') 'lag_bubbles_post_process_', t_step, '.dat'
841 file_loc = trim(case_dir) // '/lag_bubbles_post_process/' // trim(file_loc)
842
843 if (proc_rank == 0) then
844 open (unit=29, file=file_loc, form='formatted', position='rewind')
845
846 if (lag_header) then
847 write (29, '(A)', advance='no')
848 if (lag_id_wrt) write (29, '(A8)', advance='no') 'id, '
849 if (lag_pos_wrt) write (29, '(3(A17))', advance='no') 'px, ', 'py, ', 'pz, '
850 if (lag_pos_prev_wrt) write (29, '(3(A17))', advance='no') 'pvx, ', 'pvy, ', 'pvz, '
851 if (lag_vel_wrt) write (29, '(3(A17))', advance='no') 'vx, ', 'vy, ', 'vz, '
852 if (lag_rad_wrt) write (29, '(A17)', advance='no') 'radius, '
853 if (lag_rvel_wrt) write (29, '(A17)', advance='no') 'rvel, '
854 if (lag_r0_wrt) write (29, '(A17)', advance='no') 'r0, '
855 if (lag_rmax_wrt) write (29, '(A17)', advance='no') 'rmax, '
856 if (lag_rmin_wrt) write (29, '(A17)', advance='no') 'rmin, '
857 if (lag_dphidt_wrt) write (29, '(A17)', advance='no') 'dphidt, '
858 if (lag_pres_wrt) write (29, '(A17)', advance='no') 'pressure, '
859 if (lag_mv_wrt) write (29, '(A17)', advance='no') 'mv, '
860 if (lag_mg_wrt) write (29, '(A17)', advance='no') 'mg, '
861 if (lag_betat_wrt) write (29, '(A17)', advance='no') 'betaT, '
862 if (lag_betac_wrt) write (29, '(A17)', advance='no') 'betaC, '
863 write (29, '(A15)') 'time'
864 end if
865
866 do i = 1, file_tot_part
867 id = int(mpi_io_data_lg_bubbles(i, 1))
868 inputvals(1:20) = mpi_io_data_lg_bubbles(i,2:21)
869 if (id > 0) then
870 write (29, '(100(A))', advance='no') ''
871 if (lag_id_wrt) write (29, '(I6, A)', advance='no') id, ', '
872 if (lag_pos_wrt) write (29, '(3(E15.7, A))', advance='no') inputvals(1), ', ', inputvals(2), ', ', &
873 & inputvals(3), ', '
874 if (lag_pos_prev_wrt) write (29, '(3(E15.7, A))', advance='no') inputvals(4), ', ', inputvals(5), ', ', &
875 & inputvals(6), ', '
876 if (lag_vel_wrt) write (29, '(3(E15.7, A))', advance='no') inputvals(7), ', ', inputvals(8), ', ', &
877 & inputvals(9), ', '
878 if (lag_rad_wrt) write (29, '(E15.7, A)', advance='no') inputvals(10), ', '
879 if (lag_rvel_wrt) write (29, '(E15.7, A)', advance='no') inputvals(11), ', '
880 if (lag_r0_wrt) write (29, '(E15.7, A)', advance='no') inputvals(12), ', '
881 if (lag_rmax_wrt) write (29, '(E15.7, A)', advance='no') inputvals(13), ', '
882 if (lag_rmin_wrt) write (29, '(E15.7, A)', advance='no') inputvals(14), ', '
883 if (lag_dphidt_wrt) write (29, '(E15.7, A)', advance='no') inputvals(15), ', '
884 if (lag_pres_wrt) write (29, '(E15.7, A)', advance='no') inputvals(16), ', '
885 if (lag_mv_wrt) write (29, '(E15.7, A)', advance='no') inputvals(17), ', '
886 if (lag_mg_wrt) write (29, '(E15.7, A)', advance='no') inputvals(18), ', '
887 if (lag_betat_wrt) write (29, '(E15.7, A)', advance='no') inputvals(19), ', '
888 if (lag_betac_wrt) write (29, '(E15.7, A)', advance='no') inputvals(20), ', '
889 write (29, '(E15.7)') time_real
890 end if
891 end do
892 close (29)
893 end if
894
895 deallocate (mpi_io_data_lg_bubbles)
896
897 call s_mpi_barrier()
898
899 call mpi_file_close(ifile, ierr)
900 else
901 call mpi_type_contiguous(0, mpi_p, view, ierr)
902 call mpi_type_commit(view, ierr)
903
904 call mpi_file_open(mpi_comm_world, file_loc, mpi_mode_rdonly, mpi_info_int, ifile, ierr)
905
906 disp = int(sizeof(file_tot_part) + 2*sizeof(file_time) + sizeof(file_num_procs) &
907 & + file_num_procs*sizeof(proc_bubble_counts(1)), mpi_offset_kind)
908 call mpi_file_set_view(ifile, disp, mpi_p, view, 'native', mpi_info_null, ierr)
909
910 call mpi_file_read_all(ifile, lag_io_null, 0, mpi_p, status, ierr)
911
912 call mpi_file_close(ifile, ierr)
913 call mpi_type_free(view, ierr)
914 end if
915#endif
916
917 end subroutine s_write_lag_bubbles_results_to_text
918
919 !> Read Lagrangian bubble restart data and write bubble positions and scalar fields to the Silo database.
920 impure subroutine s_write_lag_bubbles_to_formatted_database_file(t_step)
921
922 integer, intent(in) :: t_step
923 character(len=len_trim(case_dir) + 3*name_len) :: file_loc
924 integer :: id
925
926#ifdef MFC_MPI
927 real(wp) :: time_real
928 integer, dimension(MPI_STATUS_SIZE) :: status
929 integer(KIND=MPI_OFFSET_KIND) :: disp
930 integer :: view
931 logical :: file_exist
932 integer, dimension(2) :: gsizes, lsizes, start_idx_part
933 integer :: ifile, ierr, nbub
934 real(wp) :: file_time, file_dt
935 integer :: file_num_procs, file_tot_part
936 integer, dimension(:), allocatable :: proc_bubble_counts
937 real(wp), dimension(1:1,1:lag_io_vars) :: dummy
938 character(LEN=4*name_len), dimension(num_procs) :: meshnames
939 integer, dimension(num_procs) :: meshtypes
940 real(wp) :: dummy_data
941 integer :: i
942 real(wp), dimension(:), allocatable :: bub_id
943 real(wp), dimension(:), allocatable :: px, py, pz, ppx, ppy, ppz, vx, vy, vz
944 real(wp), dimension(:), allocatable :: radius, rvel, rnot, rmax, rmin, dphidt
945 real(wp), dimension(:), allocatable :: pressure, mv, mg, betat, betac
946
947 dummy = 0._wp
948 dummy_data = 0._wp
949
950 ! Construct file path
951 write (file_loc, '(A,I0,A)') 'lag_bubbles_', t_step, '.dat'
952 file_loc = trim(case_dir) // '/restart_data' // trim(mpiiofs) // trim(file_loc)
953
954 ! Check if file exists
955 inquire (file=trim(file_loc), exist=file_exist)
956 if (.not. file_exist) then
957 call s_mpi_abort('Restart file ' // trim(file_loc) // ' does not exist!')
958 end if
959
960 if (.not. parallel_io) return
961
962 if (proc_rank == 0) then
963 call mpi_file_open(mpi_comm_self, file_loc, mpi_mode_rdonly, mpi_info_int, ifile, ierr)
964
965 call mpi_file_read(ifile, file_tot_part, 1, mpi_integer, status, ierr)
966 call mpi_file_read(ifile, file_time, 1, mpi_p, status, ierr)
967 call mpi_file_read(ifile, file_dt, 1, mpi_p, status, ierr)
968 call mpi_file_read(ifile, file_num_procs, 1, mpi_integer, status, ierr)
969
970 call mpi_file_close(ifile, ierr)
971 end if
972
973 call mpi_bcast(file_tot_part, 1, mpi_integer, 0, mpi_comm_world, ierr)
974 call mpi_bcast(file_time, 1, mpi_p, 0, mpi_comm_world, ierr)
975 call mpi_bcast(file_dt, 1, mpi_p, 0, mpi_comm_world, ierr)
976 call mpi_bcast(file_num_procs, 1, mpi_integer, 0, mpi_comm_world, ierr)
977 time_real = file_time
978
979 allocate (proc_bubble_counts(file_num_procs))
980
981 if (proc_rank == 0) then
982 call mpi_file_open(mpi_comm_self, file_loc, mpi_mode_rdonly, mpi_info_int, ifile, ierr)
983
984 ! Skip to processor counts position
985 disp = int(sizeof(file_tot_part) + 2*sizeof(file_time) + sizeof(file_num_procs), mpi_offset_kind)
986 call mpi_file_seek(ifile, disp, mpi_seek_set, ierr)
987 call mpi_file_read(ifile, proc_bubble_counts, file_num_procs, mpi_integer, status, ierr)
988
989 call mpi_file_close(ifile, ierr)
990 end if
991
992 call mpi_bcast(proc_bubble_counts, file_num_procs, mpi_integer, 0, mpi_comm_world, ierr)
993
994 ! Set time variables from file
995
996 nbub = proc_bubble_counts(proc_rank + 1)
997
998 start_idx_part(1) = 0
999 do i = 1, proc_rank
1000 start_idx_part(1) = start_idx_part(1) + proc_bubble_counts(i)
1001 end do
1002
1003 start_idx_part(2) = 0
1004 lsizes(1) = nbub
1005 lsizes(2) = lag_io_vars
1006
1007 gsizes(1) = file_tot_part
1008 gsizes(2) = lag_io_vars
1009
1010 if (nbub > 0) then
1011# 944 "/home/runner/work/MFC/MFC/src/post_process/m_data_output.fpp"
1012 allocate (bub_id(nbub))
1013# 944 "/home/runner/work/MFC/MFC/src/post_process/m_data_output.fpp"
1014 allocate (px(nbub))
1015# 944 "/home/runner/work/MFC/MFC/src/post_process/m_data_output.fpp"
1016 allocate (py(nbub))
1017# 944 "/home/runner/work/MFC/MFC/src/post_process/m_data_output.fpp"
1018 allocate (pz(nbub))
1019# 944 "/home/runner/work/MFC/MFC/src/post_process/m_data_output.fpp"
1020 allocate (ppx(nbub))
1021# 944 "/home/runner/work/MFC/MFC/src/post_process/m_data_output.fpp"
1022 allocate (ppy(nbub))
1023# 944 "/home/runner/work/MFC/MFC/src/post_process/m_data_output.fpp"
1024 allocate (ppz(nbub))
1025# 944 "/home/runner/work/MFC/MFC/src/post_process/m_data_output.fpp"
1026 allocate (vx(nbub))
1027# 944 "/home/runner/work/MFC/MFC/src/post_process/m_data_output.fpp"
1028 allocate (vy(nbub))
1029# 944 "/home/runner/work/MFC/MFC/src/post_process/m_data_output.fpp"
1030 allocate (vz(nbub))
1031# 944 "/home/runner/work/MFC/MFC/src/post_process/m_data_output.fpp"
1032 allocate (radius(nbub))
1033# 944 "/home/runner/work/MFC/MFC/src/post_process/m_data_output.fpp"
1034 allocate (rvel(nbub))
1035# 944 "/home/runner/work/MFC/MFC/src/post_process/m_data_output.fpp"
1036 allocate (rnot(nbub))
1037# 944 "/home/runner/work/MFC/MFC/src/post_process/m_data_output.fpp"
1038 allocate (rmax(nbub))
1039# 944 "/home/runner/work/MFC/MFC/src/post_process/m_data_output.fpp"
1040 allocate (rmin(nbub))
1041# 944 "/home/runner/work/MFC/MFC/src/post_process/m_data_output.fpp"
1042 allocate (dphidt(nbub))
1043# 944 "/home/runner/work/MFC/MFC/src/post_process/m_data_output.fpp"
1044 allocate (pressure(nbub))
1045# 944 "/home/runner/work/MFC/MFC/src/post_process/m_data_output.fpp"
1046 allocate (mv(nbub))
1047# 944 "/home/runner/work/MFC/MFC/src/post_process/m_data_output.fpp"
1048 allocate (mg(nbub))
1049# 944 "/home/runner/work/MFC/MFC/src/post_process/m_data_output.fpp"
1050 allocate (betat(nbub))
1051# 944 "/home/runner/work/MFC/MFC/src/post_process/m_data_output.fpp"
1052 allocate (betac(nbub))
1053# 946 "/home/runner/work/MFC/MFC/src/post_process/m_data_output.fpp"
1054 allocate (mpi_io_data_lg_bubbles(nbub,1:lag_io_vars))
1055
1056 call mpi_type_create_subarray(2, gsizes, lsizes, start_idx_part, mpi_order_fortran, mpi_p, view, ierr)
1057 call mpi_type_commit(view, ierr)
1058
1059 call mpi_file_open(mpi_comm_world, file_loc, mpi_mode_rdonly, mpi_info_int, ifile, ierr)
1060
1061 ! Skip extended header
1062 disp = int(sizeof(file_tot_part) + 2*sizeof(file_time) + sizeof(file_num_procs) &
1063 & + file_num_procs*sizeof(proc_bubble_counts(1)), mpi_offset_kind)
1064 call mpi_file_set_view(ifile, disp, mpi_p, view, 'native', mpi_info_int, ierr)
1065
1066 call mpi_file_read_all(ifile, mpi_io_data_lg_bubbles, lag_io_vars*nbub, mpi_p, status, ierr)
1067
1068 call mpi_file_close(ifile, ierr)
1069 call mpi_type_free(view, ierr)
1070
1071 ! Extract data from MPI_IO_DATA_lg_bubbles array Adjust these indices based on your actual data layout
1072# 968 "/home/runner/work/MFC/MFC/src/post_process/m_data_output.fpp"
1073 bub_id(:) = mpi_io_data_lg_bubbles(:,1)
1074# 968 "/home/runner/work/MFC/MFC/src/post_process/m_data_output.fpp"
1075 px(:) = mpi_io_data_lg_bubbles(:,2)
1076# 968 "/home/runner/work/MFC/MFC/src/post_process/m_data_output.fpp"
1077 py(:) = mpi_io_data_lg_bubbles(:,3)
1078# 968 "/home/runner/work/MFC/MFC/src/post_process/m_data_output.fpp"
1079 pz(:) = mpi_io_data_lg_bubbles(:,4)
1080# 968 "/home/runner/work/MFC/MFC/src/post_process/m_data_output.fpp"
1081 ppx(:) = mpi_io_data_lg_bubbles(:,5)
1082# 968 "/home/runner/work/MFC/MFC/src/post_process/m_data_output.fpp"
1083 ppy(:) = mpi_io_data_lg_bubbles(:,6)
1084# 968 "/home/runner/work/MFC/MFC/src/post_process/m_data_output.fpp"
1085 ppz(:) = mpi_io_data_lg_bubbles(:,7)
1086# 968 "/home/runner/work/MFC/MFC/src/post_process/m_data_output.fpp"
1087 vx(:) = mpi_io_data_lg_bubbles(:,8)
1088# 968 "/home/runner/work/MFC/MFC/src/post_process/m_data_output.fpp"
1089 vy(:) = mpi_io_data_lg_bubbles(:,9)
1090# 968 "/home/runner/work/MFC/MFC/src/post_process/m_data_output.fpp"
1091 vz(:) = mpi_io_data_lg_bubbles(:,10)
1092# 968 "/home/runner/work/MFC/MFC/src/post_process/m_data_output.fpp"
1093 radius(:) = mpi_io_data_lg_bubbles(:,11)
1094# 968 "/home/runner/work/MFC/MFC/src/post_process/m_data_output.fpp"
1095 rvel(:) = mpi_io_data_lg_bubbles(:,12)
1096# 968 "/home/runner/work/MFC/MFC/src/post_process/m_data_output.fpp"
1097 rnot(:) = mpi_io_data_lg_bubbles(:,13)
1098# 968 "/home/runner/work/MFC/MFC/src/post_process/m_data_output.fpp"
1099 rmax(:) = mpi_io_data_lg_bubbles(:,14)
1100# 968 "/home/runner/work/MFC/MFC/src/post_process/m_data_output.fpp"
1101 rmin(:) = mpi_io_data_lg_bubbles(:,15)
1102# 968 "/home/runner/work/MFC/MFC/src/post_process/m_data_output.fpp"
1103 dphidt(:) = mpi_io_data_lg_bubbles(:,16)
1104# 968 "/home/runner/work/MFC/MFC/src/post_process/m_data_output.fpp"
1105 pressure(:) = mpi_io_data_lg_bubbles(:,17)
1106# 968 "/home/runner/work/MFC/MFC/src/post_process/m_data_output.fpp"
1107 mv(:) = mpi_io_data_lg_bubbles(:,18)
1108# 968 "/home/runner/work/MFC/MFC/src/post_process/m_data_output.fpp"
1109 mg(:) = mpi_io_data_lg_bubbles(:,19)
1110# 968 "/home/runner/work/MFC/MFC/src/post_process/m_data_output.fpp"
1111 betat(:) = mpi_io_data_lg_bubbles(:,20)
1112# 968 "/home/runner/work/MFC/MFC/src/post_process/m_data_output.fpp"
1113 betac(:) = mpi_io_data_lg_bubbles(:,21)
1114# 970 "/home/runner/work/MFC/MFC/src/post_process/m_data_output.fpp"
1115
1116 ! Next, the root processor proceeds to record all of the spatial extents in the formatted database master file. In
1117 ! addition, it also records a sub-domain connectivity map so that the entire grid may be reassembled by looking at the
1118 ! master file.
1119 if (proc_rank == 0) then
1120 do i = 1, num_procs
1121 write (meshnames(i), '(A,I0,A,I0,A)') '../p', i - 1, '/', t_step, '.silo:lag_bubbles'
1122 meshtypes(i) = db_pointmesh
1123 end do
1124 err = dbset2dstrlen(len(meshnames(1)))
1125 err = dbputmmesh(out%dbroot, 'lag_bubbles', 16, num_procs, meshnames, len_trim(meshnames), meshtypes, db_f77null, &
1126 & ierr)
1127 end if
1128
1129 err = dbputpm(out%dbfile, 'lag_bubbles', 11, 3, px, py, pz, nbub, db_double, db_f77null, ierr)
1130
1131 if (lag_id_wrt) call s_write_lag_variable_to_formatted_database_file('part_id', t_step, bub_id, nbub)
1132 if (lag_vel_wrt) then
1133 call s_write_lag_variable_to_formatted_database_file('part_vel1', t_step, vx, nbub)
1134 call s_write_lag_variable_to_formatted_database_file('part_vel2', t_step, vy, nbub)
1135 if (p > 0) call s_write_lag_variable_to_formatted_database_file('part_vel3', t_step, vz, nbub)
1136 end if
1137 if (lag_rad_wrt) call s_write_lag_variable_to_formatted_database_file('part_radius', t_step, radius, nbub)
1138 if (lag_rvel_wrt) call s_write_lag_variable_to_formatted_database_file('part_rdot', t_step, rvel, nbub)
1139 if (lag_r0_wrt) call s_write_lag_variable_to_formatted_database_file('part_r0', t_step, rnot, nbub)
1140 if (lag_rmax_wrt) call s_write_lag_variable_to_formatted_database_file('part_rmax', t_step, rmax, nbub)
1141 if (lag_rmin_wrt) call s_write_lag_variable_to_formatted_database_file('part_rmin', t_step, rmin, nbub)
1142 if (lag_dphidt_wrt) call s_write_lag_variable_to_formatted_database_file('part_dphidt', t_step, dphidt, nbub)
1143 if (lag_pres_wrt) call s_write_lag_variable_to_formatted_database_file('part_pressure', t_step, pressure, nbub)
1144 if (lag_mv_wrt) call s_write_lag_variable_to_formatted_database_file('part_mv', t_step, mv, nbub)
1145 if (lag_mg_wrt) call s_write_lag_variable_to_formatted_database_file('part_mg', t_step, mg, nbub)
1146 if (lag_betat_wrt) call s_write_lag_variable_to_formatted_database_file('part_betaT', t_step, betat, nbub)
1147 if (lag_betac_wrt) call s_write_lag_variable_to_formatted_database_file('part_betaC', t_step, betac, nbub)
1148
1149 deallocate (bub_id, px, py, pz, ppx, ppy, ppz, vx, vy, vz, radius, rvel, rnot, rmax, rmin, dphidt, pressure, mv, mg, &
1150 & betat, betac)
1151 deallocate (mpi_io_data_lg_bubbles)
1152 else
1153 call mpi_type_contiguous(0, mpi_p, view, ierr)
1154 call mpi_type_commit(view, ierr)
1155
1156 call mpi_file_open(mpi_comm_world, file_loc, mpi_mode_rdonly, mpi_info_int, ifile, ierr)
1157
1158 ! Skip extended header
1159 disp = int(sizeof(file_tot_part) + 2*sizeof(file_time) + sizeof(file_num_procs) &
1160 & + file_num_procs*sizeof(proc_bubble_counts(1)), mpi_offset_kind)
1161 call mpi_file_set_view(ifile, disp, mpi_p, view, 'native', mpi_info_int, ierr)
1162
1163 call mpi_file_read_all(ifile, dummy, 0, mpi_p, status, ierr)
1164
1165 call mpi_file_close(ifile, ierr)
1166 call mpi_type_free(view, ierr)
1167
1168 if (proc_rank == 0) then
1169 do i = 1, num_procs
1170 write (meshnames(i), '(A,I0,A,I0,A)') '../p', i - 1, '/', t_step, '.silo:lag_bubbles'
1171 meshtypes(i) = db_pointmesh
1172 end do
1173 err = dbset2dstrlen(len(meshnames(1)))
1174 err = dbputmmesh(out%dbroot, 'lag_bubbles', 16, num_procs, meshnames, len_trim(meshnames), meshtypes, db_f77null, &
1175 & ierr)
1176 end if
1177
1178 err = dbsetemptyok(1)
1179 err = dbputpm(out%dbfile, 'lag_bubbles', 11, 3, dummy_data, dummy_data, dummy_data, 0, db_double, db_f77null, ierr)
1180
1181 if (lag_id_wrt) call s_write_lag_variable_to_formatted_database_file('part_id', t_step)
1182 if (lag_vel_wrt) then
1183 call s_write_lag_variable_to_formatted_database_file('part_vel1', t_step)
1184 call s_write_lag_variable_to_formatted_database_file('part_vel2', t_step)
1185 if (p > 0) call s_write_lag_variable_to_formatted_database_file('part_vel3', t_step)
1186 end if
1187 if (lag_rad_wrt) call s_write_lag_variable_to_formatted_database_file('part_radius', t_step)
1188 if (lag_rvel_wrt) call s_write_lag_variable_to_formatted_database_file('part_rdot', t_step)
1189 if (lag_r0_wrt) call s_write_lag_variable_to_formatted_database_file('part_r0', t_step)
1190 if (lag_rmax_wrt) call s_write_lag_variable_to_formatted_database_file('part_rmax', t_step)
1191 if (lag_rmin_wrt) call s_write_lag_variable_to_formatted_database_file('part_rmin', t_step)
1192 if (lag_dphidt_wrt) call s_write_lag_variable_to_formatted_database_file('part_dphidt', t_step)
1193 if (lag_pres_wrt) call s_write_lag_variable_to_formatted_database_file('part_pressure', t_step)
1194 if (lag_mv_wrt) call s_write_lag_variable_to_formatted_database_file('part_mv', t_step)
1195 if (lag_mg_wrt) call s_write_lag_variable_to_formatted_database_file('part_mg', t_step)
1196 if (lag_betat_wrt) call s_write_lag_variable_to_formatted_database_file('part_betaT', t_step)
1197 if (lag_betac_wrt) call s_write_lag_variable_to_formatted_database_file('part_betaC', t_step)
1198 end if
1199#endif
1200
1201 end subroutine s_write_lag_bubbles_to_formatted_database_file
1202
1203 !> Write a single Lagrangian bubble point-variable to the Silo database slave and master files.
1204 subroutine s_write_lag_variable_to_formatted_database_file(varname, t_step, data, nBubs)
1205
1206 character(len=*), intent(in) :: varname
1207 integer, intent(in) :: t_step
1208 real(wp), dimension(1:), intent(in), optional :: data
1209 integer, intent(in), optional :: nBubs
1210 character(len=64), dimension(num_procs) :: var_names
1211 integer, dimension(num_procs) :: var_types
1212 real(wp) :: dummy_data
1213 integer :: ierr
1214 integer :: i
1215
1216 dummy_data = 0._wp
1217
1218 if (present(nbubs) .and. present(data)) then
1219 if (proc_rank == 0) then
1220 do i = 1, num_procs
1221 write (var_names(i), '(A,I0,A,I0,A)') '../p', i - 1, '/', t_step, '.silo:' // trim(varname)
1222 var_types(i) = db_pointvar
1223 end do
1224 err = dbset2dstrlen(len(var_names(1)))
1225 err = dbputmvar(out%dbroot, trim(varname), len_trim(varname), num_procs, var_names, len_trim(var_names), &
1226 & var_types, db_f77null, ierr)
1227 end if
1228
1229 err = dbputpv1(out%dbfile, trim(varname), len_trim(varname), 'lag_bubbles', 11, data, nbubs, db_double, db_f77null, &
1230 & ierr)
1231 else
1232 if (proc_rank == 0) then
1233 do i = 1, num_procs
1234 write (var_names(i), '(A,I0,A,I0,A)') '../p', i - 1, '/', t_step, '.silo:' // trim(varname)
1235 var_types(i) = db_pointvar
1236 end do
1237 err = dbset2dstrlen(len(var_names(1)))
1238 err = dbsetemptyok(1)
1239 err = dbputmvar(out%dbroot, trim(varname), len_trim(varname), num_procs, var_names, len_trim(var_names), &
1240 & var_types, db_f77null, ierr)
1241 end if
1242
1243 err = dbsetemptyok(1)
1244 err = dbputpv1(out%dbfile, trim(varname), len_trim(varname), 'lag_bubbles', 11, dummy_data, 0, db_double, db_f77null, &
1245 & ierr)
1246 end if
1247
1248 end subroutine s_write_lag_variable_to_formatted_database_file
1249
1250 !> Convert the binary immersed-boundary state file to per-body formatted text files
1251 impure subroutine s_write_ib_state_files()
1252
1253 character(len=len_trim(case_dir) + 4*name_len) :: in_file, out_file, file_loc
1254 integer :: iu_in, ios, i, rec_id
1255 integer, allocatable, dimension(:) :: iu_out
1256 real(wp) :: rec_time
1257 real(wp), dimension(3) :: rec_force, rec_torque
1258 real(wp), dimension(3) :: rec_vel, rec_angular_vel
1259 real(wp), dimension(3) :: rec_angles, rec_centroid
1260
1261 file_loc = trim(case_dir) // '/D'
1262
1263 in_file = trim(file_loc) // '/ib_state.dat'
1264 open (newunit=iu_in, file=trim(in_file), form='unformatted', access='stream', status='old', action='read', iostat=ios)
1265 if (ios /= 0) then
1266 call s_mpi_abort('Cannot open IB state input file: ' // trim(in_file))
1267 end if
1268
1269 allocate (iu_out(num_ibs))
1270 do i = 1, num_ibs
1271 write (out_file, '(A,I0,A)') trim(file_loc) // '/ib_', i, '.txt'
1272 open (newunit=iu_out(i), file=trim(out_file), form='formatted', status='replace', action='write', iostat=ios)
1273 if (ios /= 0) then
1274 call s_mpi_abort('Cannot open IB state output file: ' // trim(out_file))
1275 end if
1276 write (iu_out(i), &
1277 & '(A)') 'mytime fx fy fz Tau_x Tau_y Tau_z vx vy vz omega_x omega_y omega_z angle_x angle_y angle_z x_c y_c z_c'
1278 end do
1279
1280 do
1281 read (iu_in, iostat=ios) rec_time, rec_id, rec_force, rec_torque, rec_vel, rec_angular_vel, rec_angles, &
1282 & rec_centroid(1), rec_centroid(2), rec_centroid(3)
1283 if (ios /= 0) exit
1284
1285 if (rec_id >= 1 .and. rec_id <= num_ibs) then
1286 write (iu_out(rec_id), '(19(ES24.16E3,1X))') rec_time, rec_force(1), rec_force(2), rec_force(3), rec_torque(1), &
1287 & rec_torque(2), rec_torque(3), rec_vel(1), rec_vel(2), rec_vel(3), rec_angular_vel(1), &
1288 & rec_angular_vel(2), rec_angular_vel(3), rec_angles(1), rec_angles(2), rec_angles(3), rec_centroid(1), &
1289 & rec_centroid(2), rec_centroid(3)
1290 end if
1291 end do
1292
1293 close (iu_in)
1294 do i = 1, num_ibs
1295 close (iu_out(i))
1296 end do
1297 deallocate (iu_out)
1298
1299 end subroutine s_write_ib_state_files
1300
1301 !> Extract the volume-fraction interface contour from primitive fields and write the coordinates to the interface data file.
1302 impure subroutine s_write_intf_data_file(q_prim_vf)
1303
1304 type(scalar_field), dimension(sys_size), intent(in) :: q_prim_vf
1305 integer :: i, j, k, l, cent
1306 integer :: counter, root !< number of data points extracted to fit shape to SH perturbations
1307 real(wp), allocatable :: x_td(:), y_td(:), x_d1(:), y_d1(:), y_d(:), x_d(:)
1308 real(wp) :: axp, axm, ayp, aym, tgp, euc_d, thres, maxalph_loc, maxalph_glb
1309
1310 allocate (x_d1(m*n))
1311 allocate (y_d1(m*n))
1312 counter = 0
1313 maxalph_loc = 0._wp
1314 do k = 0, p
1315 do j = 0, n
1316 do i = 0, m
1317 if (q_prim_vf(eqn_idx%E + 2)%sf(i, j, k) > maxalph_loc) then
1318 maxalph_loc = q_prim_vf(eqn_idx%E + 2)%sf(i, j, k)
1319 end if
1320 end do
1321 end do
1322 end do
1323
1324 call s_mpi_allreduce_max(maxalph_loc, maxalph_glb)
1325 if (p > 0) then
1326 do l = 0, p
1327 if (z_cc(l) < dz(l) .and. z_cc(l) > 0) then
1328 cent = l
1329 end if
1330 end do
1331 else
1332 cent = 0
1333 end if
1334
1335 thres = 0.9_wp*maxalph_glb
1336 do k = 0, n
1337 do j = 0, m
1338 axp = q_prim_vf(eqn_idx%E + 2)%sf(j + 1, k, cent)
1339 axm = q_prim_vf(eqn_idx%E + 2)%sf(j, k, cent)
1340 ayp = q_prim_vf(eqn_idx%E + 2)%sf(j, k + 1, cent)
1341 aym = q_prim_vf(eqn_idx%E + 2)%sf(j, k, cent)
1342 if ((axp > thres .and. axm < thres) .or. (axp < thres .and. axm > thres) .or. (ayp > thres .and. aym < thres) &
1343 & .or. (ayp < thres .and. aym > thres)) then
1344 if (counter == 0) then
1345 counter = counter + 1
1346 x_d1(counter) = x_cc(j)
1347 y_d1(counter) = y_cc(k)
1348 else
1349 tgp = sqrt(dx(j)**2 + dy(k)**2)
1350 do i = 1, counter
1351 euc_d = sqrt((x_cc(j) - x_d1(i))**2 + (y_cc(k) - y_d1(i))**2)
1352 if (euc_d < tgp) then
1353 exit
1354 else if (i == counter) then
1355 counter = counter + 1
1356 x_d1(counter) = x_cc(j)
1357 y_d1(counter) = y_cc(k)
1358 end if
1359 end do
1360 end if
1361 end if
1362 end do
1363 end do
1364
1365 allocate (x_d(counter), y_d(counter))
1366
1367 do i = 1, counter
1368 y_d(i) = y_d1(i)
1369 x_d(i) = x_d1(i)
1370 end do
1371 root = 0
1372
1373 call s_mpi_gather_data(x_d, counter, x_td, root)
1374 call s_mpi_gather_data(y_d, counter, y_td, root)
1375 if (proc_rank == 0) then
1376 do i = 1, size(x_td)
1377 if (i == size(x_td)) then
1378 write (211, '(F12.9,1X,F12.9,1X,I4)') x_td(i), y_td(i), size(x_td)
1379 else
1380 write (211, '(F12.9,1X,F12.9,1X,F3.1)') x_td(i), y_td(i), 0._wp
1381 end if
1382 end do
1383 end if
1384
1385 end subroutine s_write_intf_data_file
1386
1387 !> Compute volume-integrated kinetic, potential, and internal energies and write the energy budget to the energy data file.
1388 impure subroutine s_write_energy_data_file(q_prim_vf, q_cons_vf)
1389
1390 type(scalar_field), dimension(sys_size), intent(in) :: q_prim_vf, q_cons_vf
1391 real(wp) :: elk, egk, elp, egint, eg_phase, vb, vl, pres_av, et
1392 real(wp) :: rho, pres, dv, tmp, gamma, pi_inf, qv, maxma, maxma_glb, maxvel, c, ma
1393 real(wp), dimension(num_vels) :: vel
1394 real(wp), dimension(num_fluids) :: adv, alpha_rho
1395 integer :: i, j, k, l, s !< looping indices
1396
1397 egk = 0._wp
1398 elp = 0._wp
1399 egint = 0._wp
1400 vb = 0._wp
1401 maxvel = 0._wp
1402 maxma = 0._wp
1403 vl = 0._wp
1404 elk = 0._wp
1405 et = 0._wp
1406 vb = 0._wp
1407 dv = 0._wp
1408 pres_av = 0._wp
1409 pres = 0._wp
1410 c = 0._wp
1411
1412 do k = 0, p
1413 do j = 0, n
1414 do i = 0, m
1415 dv = dx(i)*dy(j)*dz(k)
1416 pres = q_prim_vf(eqn_idx%E)%sf(i, j, k)
1417 do s = 1, num_vels
1418 vel(s) = q_prim_vf(num_fluids + s)%sf(i, j, k)
1419 egk = egk + 0.5_wp*q_prim_vf(eqn_idx%E + 2)%sf(i, j, k)*q_prim_vf(2)%sf(i, j, k)*vel(s)*vel(s)*dv
1420 elk = elk + 0.5_wp*q_prim_vf(eqn_idx%E + 1)%sf(i, j, k)*q_prim_vf(1)%sf(i, j, k)*vel(s)*vel(s)*dv
1421 if (abs(vel(s)) > maxvel) then
1422 maxvel = abs(vel(s))
1423 end if
1424 end do
1425 do l = 1, num_fluids
1426 adv(l) = q_prim_vf(eqn_idx%E + l)%sf(i, j, k)
1427 alpha_rho(l) = q_prim_vf(l)%sf(i, j, k)
1428 end do
1429
1430 call s_phase_internal_energy(pres, adv(2), alpha_rho(2), 2, eg_phase)
1431 egint = egint + eg_phase*dv
1432
1433 call s_compute_mixture_coefficients(alpha_rho, adv, rho, gamma, pi_inf, qv)
1434
1435 call s_compute_speed_of_sound(pres, rho, gamma, pi_inf, adv, c, alpha_rho)
1436
1437 ma = maxvel/c
1438 if (ma > maxma .and. (adv(1) > (1.0_wp - 1.0e-10_wp))) then
1439 maxma = ma
1440 end if
1441 vl = vl + adv(1)*dv
1442 vb = vb + adv(2)*dv
1443 pres_av = pres_av + adv(1)*pres*dv
1444 et = et + q_cons_vf(eqn_idx%E)%sf(i, j, k)*dv
1445 end do
1446 end do
1447 end do
1448
1449 tmp = pres_av
1450 call s_mpi_allreduce_sum(tmp, pres_av)
1451 tmp = vl
1452 call s_mpi_allreduce_sum(tmp, vl)
1453
1454 call s_mpi_allreduce_max(maxma, maxma_glb)
1455 tmp = elk
1456 call s_mpi_allreduce_sum(tmp, elk)
1457 tmp = egint
1458 call s_mpi_allreduce_sum(tmp, egint)
1459 tmp = egk
1460 call s_mpi_allreduce_sum(tmp, egk)
1461 tmp = vb
1462 call s_mpi_allreduce_sum(tmp, vb)
1463 tmp = et
1464 call s_mpi_allreduce_sum(tmp, et)
1465
1466 elp = pres_av/vl*vb
1467 if (proc_rank == 0) then
1468 write (251, '(10X, 8F24.8)') elp, egint, elk, egk, et, vb, vl, maxma_glb
1469 end if
1470
1471 end subroutine s_write_energy_data_file
1472
1473 !> Read IB state and write a Silo point mesh with per-body scalar fields.
1474 impure subroutine s_write_ib_bodies_to_formatted_database_file(t_step)
1475
1476 integer, intent(in) :: t_step
1477 character(len=len_trim(case_dir) + 3*name_len) :: file_loc
1478
1479#ifdef MFC_MPI
1480 integer, parameter :: nfields_per_ib = 20
1481 real(wp) :: ib_buf(nfields_per_ib)
1482 real(wp), dimension(:,:), allocatable :: ib_data
1483 logical :: file_exist
1484 character(LEN=4*name_len), dimension(num_procs) :: meshnames
1485 integer, dimension(num_procs) :: meshtypes
1486 integer :: i, ios, file_unit
1487 integer :: ierr, nbodies
1488 integer :: r, nlocal, gbl_id
1489 character(len=10) :: t_step_string
1490 real(wp), dimension(:), allocatable :: px, py, pz
1491 real(wp), dimension(:), allocatable :: force_x, force_y, force_z
1492 real(wp), dimension(:), allocatable :: torque_x, torque_y, torque_z
1493 real(wp), dimension(:), allocatable :: vel_x, vel_y, vel_z
1494 real(wp), dimension(:), allocatable :: omega_x, omega_y, omega_z
1495 real(wp), dimension(:), allocatable :: angle_x, angle_y, angle_z
1496 real(wp), dimension(:), allocatable :: ib_diameter
1497
1498 if (proc_rank == 0) then
1499 nbodies = num_ibs
1500
1501 if (nbodies > 0) then
1502 allocate (ib_data(nbodies, nfields_per_ib))
1503 allocate (px(nbodies), py(nbodies), pz(nbodies))
1504 allocate (force_x(nbodies), force_y(nbodies), force_z(nbodies))
1505 allocate (torque_x(nbodies), torque_y(nbodies), torque_z(nbodies))
1506 allocate (vel_x(nbodies), vel_y(nbodies), vel_z(nbodies))
1507 allocate (omega_x(nbodies), omega_y(nbodies), omega_z(nbodies))
1508 allocate (angle_x(nbodies), angle_y(nbodies), angle_z(nbodies))
1509 allocate (ib_diameter(nbodies))
1510
1511 if (file_per_process) then
1512 call s_int_to_str(t_step, t_step_string)
1513 ib_data = 0._wp
1514 do r = 0, num_procs - 1
1515 write (file_loc, '(A,I0,A,i7.7,A)') 'ib_state_', t_step, '_', r, '.dat'
1516 file_loc = trim(case_dir) // '/restart_data/lustre_' // trim(t_step_string) // '/' // trim(file_loc)
1517
1518 inquire (file=trim(file_loc), exist=file_exist)
1519 if (.not. file_exist) call s_mpi_abort('Restart file ' // trim(file_loc) // ' does not exist!')
1520
1521 open (newunit=file_unit, file=trim(file_loc), form='unformatted', access='stream', status='old', iostat=ios)
1522 if (ios /= 0) call s_mpi_abort('Cannot open IB state file: ' // trim(file_loc))
1523
1524 read (file_unit, iostat=ios) nlocal
1525 if (ios /= 0) call s_mpi_abort('Error reading IB state file header: ' // trim(file_loc))
1526
1527 do i = 1, nlocal
1528 read (file_unit, iostat=ios) gbl_id
1529 if (ios /= 0) call s_mpi_abort('Error reading IB patch ID: ' // trim(file_loc))
1530 read (file_unit, iostat=ios) ib_buf
1531 if (ios /= 0) call s_mpi_abort('Error reading IB state data: ' // trim(file_loc))
1532 ib_data(gbl_id,:) = ib_buf(:)
1533 end do
1534
1535 close (file_unit)
1536 end do
1537 else
1538 ! Build path to per-timestep IB state file
1539 write (file_loc, '(A,I0,A)') '/restart_data/ib_state_', t_step, '.dat'
1540 file_loc = trim(case_dir) // trim(file_loc)
1541
1542 inquire (file=trim(file_loc), exist=file_exist)
1543 if (.not. file_exist) call s_mpi_abort('Restart file ' // trim(file_loc) // ' does not exist!')
1544
1545 open (newunit=file_unit, file=trim(file_loc), form='unformatted', access='stream', status='old', iostat=ios)
1546 if (ios /= 0) call s_mpi_abort('Cannot open IB state file: ' // trim(file_loc))
1547
1548 do i = 1, nbodies
1549 read (file_unit, iostat=ios) ib_buf
1550 if (ios /= 0) call s_mpi_abort('Error reading IB state file')
1551 ib_data(i,:) = ib_buf(:)
1552 end do
1553
1554 close (file_unit)
1555 end if
1556
1557 do i = 1, nbodies
1558 force_x(i) = ib_data(i, 2); force_y(i) = ib_data(i, 3); force_z(i) = ib_data(i, 4)
1559 torque_x(i) = ib_data(i, 5); torque_y(i) = ib_data(i, 6); torque_z(i) = ib_data(i, 7)
1560 vel_x(i) = ib_data(i, 8); vel_y(i) = ib_data(i, 9); vel_z(i) = ib_data(i, 10)
1561 omega_x(i) = ib_data(i, 11); omega_y(i) = ib_data(i, 12); omega_z(i) = ib_data(i, 13)
1562 angle_x(i) = ib_data(i, 14); angle_y(i) = ib_data(i, 15); angle_z(i) = ib_data(i, 16)
1563 px(i) = ib_data(i, 17); py(i) = ib_data(i, 18); pz(i) = ib_data(i, 19)
1564 ib_diameter(i) = ib_data(i, 20)*2.0_wp
1565 end do
1566
1567 write (meshnames(1), '(A,I0,A)') '../p0/', t_step, '.silo:ib_bodies'
1568 meshtypes(1) = db_pointmesh
1569 err = dbset2dstrlen(len(meshnames(1)))
1570 err = dbputmmesh(out%dbroot, 'ib_bodies', 16, 1, meshnames, len_trim(meshnames), meshtypes, db_f77null, ierr)
1571
1572 err = dbputpm(out%dbfile, 'ib_bodies', 9, 3, px, py, pz, nbodies, db_double, db_f77null, ierr)
1573
1574 call s_write_ib_variable('ib_force_x', t_step, force_x, nbodies)
1575 call s_write_ib_variable('ib_force_y', t_step, force_y, nbodies)
1576 call s_write_ib_variable('ib_force_z', t_step, force_z, nbodies)
1577 call s_write_ib_variable('ib_torque_x', t_step, torque_x, nbodies)
1578 call s_write_ib_variable('ib_torque_y', t_step, torque_y, nbodies)
1579 call s_write_ib_variable('ib_torque_z', t_step, torque_z, nbodies)
1580 call s_write_ib_variable('ib_vel_x', t_step, vel_x, nbodies)
1581 call s_write_ib_variable('ib_vel_y', t_step, vel_y, nbodies)
1582 call s_write_ib_variable('ib_vel_z', t_step, vel_z, nbodies)
1583 call s_write_ib_variable('ib_omega_x', t_step, omega_x, nbodies)
1584 call s_write_ib_variable('ib_omega_y', t_step, omega_y, nbodies)
1585 call s_write_ib_variable('ib_omega_z', t_step, omega_z, nbodies)
1586 call s_write_ib_variable('ib_angle_x', t_step, angle_x, nbodies)
1587 call s_write_ib_variable('ib_angle_y', t_step, angle_y, nbodies)
1588 call s_write_ib_variable('ib_angle_z', t_step, angle_z, nbodies)
1589 call s_write_ib_variable('ib_diameter', t_step, ib_diameter, nbodies)
1590
1591 deallocate (ib_data, px, py, pz, force_x, force_y, force_z)
1592 deallocate (torque_x, torque_y, torque_z, vel_x, vel_y, vel_z)
1593 deallocate (omega_x, omega_y, omega_z, angle_x, angle_y, angle_z)
1594 deallocate (ib_diameter)
1595 end if
1596 end if
1597#endif
1598
1599 end subroutine s_write_ib_bodies_to_formatted_database_file
1600
1601 !> Write a single IB point-variable to the Silo database slave and master files.
1602 subroutine s_write_ib_variable(varname, t_step, data, nBodies)
1603
1604 character(len=*), intent(in) :: varname
1605 integer, intent(in) :: t_step
1606 real(wp), dimension(:), intent(in) :: data
1607 integer, intent(in) :: nBodies
1608 character(len=4*name_len) :: var_name_entry
1609 integer :: var_type_entry, ierr
1610
1611 write (var_name_entry, '(A,I0,A)') '../p0/', t_step, '.silo:' // trim(varname)
1612 var_type_entry = db_pointvar
1613 err = dbset2dstrlen(len(var_name_entry))
1614 err = dbputmvar(out%dbroot, trim(varname), len_trim(varname), 1, var_name_entry, len_trim(var_name_entry), &
1615 & var_type_entry, db_f77null, ierr)
1616
1617 err = dbputpv1(out%dbfile, trim(varname), len_trim(varname), 'ib_bodies', 9, data, nbodies, db_double, db_f77null, ierr)
1618
1619 end subroutine s_write_ib_variable
1620
1621 !> Close the formatted database slave file and, for the root process, the master file.
1622 impure subroutine s_close_formatted_database_file()
1623
1624 integer :: ierr
1625
1626 if (format == format_silo) then
1627 ierr = dbclose(out%dbfile)
1628 if (proc_rank == 0) ierr = dbclose(out%dbroot)
1629 else
1630 close (out%dbfile)
1631 if (n == 0 .and. proc_rank == 0) close (out%dbroot)
1632 end if
1633
1634 end subroutine s_close_formatted_database_file
1635
1636 !> Close the interface data file.
1637 impure subroutine s_close_intf_data_file()
1638
1639 close (211)
1640
1641 end subroutine s_close_intf_data_file
1642
1643 !> Close the energy data file.
1644 impure subroutine s_close_energy_data_file()
1645
1646 close (251)
1647
1648 end subroutine s_close_energy_data_file
1649
1650 !> Deallocate module arrays and release all data-output resources.
1651 impure subroutine s_finalize_data_output_module()
1652
1653 deallocate (out%q_sf)
1654 if (n == 0) deallocate (out%q_root_sf)
1655 if (grid_geometry == 3) then
1656 deallocate (out%cyl_q_sf)
1657 end if
1658
1659 ! Deallocating spatial and data extents and also the variables for the offsets and the one bookkeeping the number of
1660 ! cell-boundaries in each active coordinate direction. Note that all these variables were only needed by Silo-HDF5 format
1661 ! for multidimensional data.
1662 if (format == format_silo) then
1663 deallocate (out%spatial_extents)
1664 deallocate (out%data_extents)
1665 deallocate (out%lo_offset)
1666 deallocate (out%hi_offset)
1667 deallocate (out%dims)
1668 end if
1669
1670 end subroutine s_finalize_data_output_module
1671
1672end module m_data_output
q_cons_vf
type(scalar_field), dimension(sys_size), intent(inout) q_cons_vf
Definition m_phase_change.fpp.f90:1006
k
integer, intent(in) k
Definition m_phase_change.fpp.f90:1008
j
integer, intent(in) j
Definition m_phase_change.fpp.f90:1008
l
integer, intent(in) l
Definition m_phase_change.fpp.f90:1008
m_compile_specific
Platform-specific file and directory operations: create, delete, inquire, getcwd, and basename.
Definition m_compile_specific.f90:6
m_compile_specific::my_inquire
impure subroutine my_inquire(fileloc, dircheck)
Inquire on the existence of a directory or file.
Definition m_compile_specific.f90:56
m_compile_specific::s_create_directory
impure subroutine s_create_directory(dir_name)
Create a directory and all its parents if it does not exist.
Definition m_compile_specific.f90:17
m_constants
Compile-time constant parameters: default values, tolerances, and physical constants.
Definition m_constants.fpp.f90:7
m_constants::lag_io_vars
integer, parameter lag_io_vars
Number of variables per particle for MPI_IO.
Definition m_constants.fpp.f90:70
m_constants::model_eqns_5eq
integer, parameter model_eqns_5eq
Definition m_constants.fpp.f90:148
m_constants::format_silo
integer, parameter format_silo
Definition m_constants.fpp.f90:142
m_constants::format_binary
integer, parameter format_binary
Definition m_constants.fpp.f90:143
m_constants::precision_single
integer, parameter precision_single
Definition m_constants.fpp.f90:158
m_constants::model_eqns_6eq
integer, parameter model_eqns_6eq
Definition m_constants.fpp.f90:149
m_constants::model_eqns_gamma_law
integer, parameter model_eqns_gamma_law
Definition m_constants.fpp.f90:147
m_data_output
Writes post-processed grid and flow-variable data to Silo-HDF5 or binary database files.
Definition m_data_output.fpp.f90:7
m_data_output::s_write_grid_to_formatted_database_file
impure subroutine, public s_write_grid_to_formatted_database_file(t_step)
Write the computational grid (cell-boundary coordinates) to the formatted database slave and master f...
Definition m_data_output.fpp.f90:456
m_data_output::s_write_variable_to_formatted_database_file
impure subroutine, public s_write_variable_to_formatted_database_file(varname, t_step)
Write a single flow variable field to the formatted database slave and master files for a given time ...
Definition m_data_output.fpp.f90:597
m_data_output::s_open_energy_data_file
impure subroutine, public s_open_energy_data_file()
Open the energy data file for appending volume-integrated energy budget quantities.
Definition m_data_output.fpp.f90:444
m_data_output::s_open_intf_data_file
impure subroutine, public s_open_intf_data_file()
Open the interface data file for appending extracted interface coordinates.
Definition m_data_output.fpp.f90:432
m_data_output::s_write_ib_state_files
impure subroutine, public s_write_ib_state_files()
Convert the binary immersed-boundary state file to per-body formatted text files.
Definition m_data_output.fpp.f90:1252
m_data_output::s_write_energy_data_file
impure subroutine, public s_write_energy_data_file(q_prim_vf, q_cons_vf)
Compute volume-integrated kinetic, potential, and internal energies and write the energy budget to th...
Definition m_data_output.fpp.f90:1389
m_data_output::s_write_lag_bubbles_to_formatted_database_file
impure subroutine, public s_write_lag_bubbles_to_formatted_database_file(t_step)
Read Lagrangian bubble restart data and write bubble positions and scalar fields to the Silo database...
Definition m_data_output.fpp.f90:921
m_data_output::s_write_intf_data_file
impure subroutine, public s_write_intf_data_file(q_prim_vf)
Extract the volume-fraction interface contour from primitive fields and write the coordinates to the ...
Definition m_data_output.fpp.f90:1303
m_data_output::s_write_ib_variable
subroutine s_write_ib_variable(varname, t_step, data, nbodies)
Write a single IB point-variable to the Silo database slave and master files.
Definition m_data_output.fpp.f90:1603
m_data_output::s_initialize_data_output_module
impure subroutine, public s_initialize_data_output_module()
Allocate storage arrays, configure output directories, and count flow variables for formatted databas...
Definition m_data_output.fpp.f90:41
m_data_output::s_close_energy_data_file
impure subroutine, public s_close_energy_data_file()
Close the energy data file.
Definition m_data_output.fpp.f90:1645
m_data_output::s_close_formatted_database_file
impure subroutine, public s_close_formatted_database_file()
Close the formatted database slave file and, for the root process, the master file.
Definition m_data_output.fpp.f90:1623
m_data_output::s_open_formatted_database_file
impure subroutine, public s_open_formatted_database_file(t_step)
Open (or create) the Silo-HDF5 or Binary formatted database slave and master files for a given time s...
Definition m_data_output.fpp.f90:367
m_data_output::s_close_intf_data_file
impure subroutine, public s_close_intf_data_file()
Close the interface data file.
Definition m_data_output.fpp.f90:1638
m_data_output::s_write_lag_bubbles_results_to_text
impure subroutine, public s_write_lag_bubbles_results_to_text(t_step)
Write the post-processed results in the folder 'lag_bubbles_data'.
Definition m_data_output.fpp.f90:752
m_data_output::s_finalize_data_output_module
impure subroutine, public s_finalize_data_output_module()
Deallocate module arrays and release all data-output resources.
Definition m_data_output.fpp.f90:1652
m_data_output::err
integer, private err
Definition m_data_output.fpp.f90:35
m_data_output::s_define_output_region
impure subroutine, public s_define_output_region
Compute the cell-index bounds for the user-specified partial output domain in each coordinate directi...
Definition m_data_output.fpp.f90:282
m_data_output::s_write_ib_bodies_to_formatted_database_file
impure subroutine, public s_write_ib_bodies_to_formatted_database_file(t_step)
Read IB state and write a Silo point mesh with per-body scalar fields.
Definition m_data_output.fpp.f90:1475
m_data_output::out
type(output_context), public out
Output workspace: flow variable buffers, VisIt extents/offsets, directory paths, file handles,...
Definition m_data_output.fpp.f90:32
m_data_output::s_write_lag_variable_to_formatted_database_file
subroutine s_write_lag_variable_to_formatted_database_file(varname, t_step, data, nbubs)
Write a single Lagrangian bubble point-variable to the Silo database slave and master files.
Definition m_data_output.fpp.f90:1205
m_derived_types
Shared derived types for field data, patch geometry, bubble dynamics, and MPI I/O structures.
Definition m_derived_types.fpp.f90:335
m_derived_variables
Computes derived flow quantities (sound speed, vorticity, Schlieren, etc.) from conservative and prim...
Definition m_derived_variables.fpp.f90:8
m_eos
Equations of state in Gamma/Pi form, rho e = Gamma(rho) p + Pi(rho).
Definition m_eos.fpp.f90:381
m_eos::s_phase_internal_energy
subroutine, public s_phase_internal_energy(pres, alpha, alpha_rho, i, e_phase)
Internal energy per unit volume of phase i at pressure pres: alpha (Gamma p + Pi) + alpha_rho qv,...
Definition m_eos.fpp.f90:1677
m_eos::s_compute_mixture_coefficients
subroutine, public s_compute_mixture_coefficients(alpha_rho_k, alpha_k, rho_k, gamma_k, pi_inf_k, qv_k)
Mixture coefficients of one state. Under bubbles_euler with num_fluids == 1 the sole advection slot a...
Definition m_eos.fpp.f90:1919
m_eos::s_compute_speed_of_sound
subroutine, public s_compute_speed_of_sound(pres, rho, gamma, pi_inf, adv, c, alpha_rho)
Speed of sound of a thermodynamic state. Enthalpy is not an argument: for a real state H,...
Definition m_eos.fpp.f90:2088
m_global_parameters
Global parameters for the post-process: domain geometry, equation of state, and output database setti...
Definition m_global_parameters.fpp.f90:19
m_global_parameters::grid_geometry
integer grid_geometry
Definition m_global_parameters.fpp.f90:53
m_global_parameters::offset_y
type(int_bounds_info) offset_y
Definition m_global_parameters.fpp.f90:137
m_global_parameters::mpi_io_data_lg_bubbles
real(wp), dimension(:,:), allocatable, public mpi_io_data_lg_bubbles
Definition m_global_parameters.fpp.f90:124
m_global_parameters::y_cc
real(wp), dimension(:), allocatable y_cc
Definition m_global_parameters.fpp.f90:69
m_global_parameters::proc_rank
integer proc_rank
Rank of the local processor.
Definition m_global_parameters.fpp.f90:41
m_global_parameters::adv
real(wp), dimension(:), allocatable adv
Advection variables.
Definition m_global_parameters.fpp.f90:128
m_global_parameters::z_output_idx
type(int_bounds_info) z_output_idx
Indices of domain to output for post-processing.
Definition m_global_parameters.fpp.f90:132
m_global_parameters::y_cb
real(wp), dimension(:), allocatable y_cb
Definition m_global_parameters.fpp.f90:64
m_global_parameters::x_output
type(bounds_info) x_output
Definition m_global_parameters.fpp.f90:131
m_global_parameters::dz
real(wp), dimension(:), allocatable dz
Definition m_global_parameters.fpp.f90:75
m_global_parameters::x_root_cb
real(wp), dimension(:), allocatable x_root_cb
Definition m_global_parameters.fpp.f90:64
m_global_parameters::z_cb
real(wp), dimension(:), allocatable z_cb
Definition m_global_parameters.fpp.f90:64
m_global_parameters::z_output
type(bounds_info) z_output
Portion of domain to output for post-processing.
Definition m_global_parameters.fpp.f90:131
m_global_parameters::x_output_idx
type(int_bounds_info) x_output_idx
Definition m_global_parameters.fpp.f90:132
m_global_parameters::x_cc
real(wp), dimension(:), allocatable x_cc
Definition m_global_parameters.fpp.f90:69
m_global_parameters::x_cb
real(wp), dimension(:), allocatable x_cb
Definition m_global_parameters.fpp.f90:64
m_global_parameters::dy
real(wp), dimension(:), allocatable dy
Definition m_global_parameters.fpp.f90:75
m_global_parameters::offset_x
type(int_bounds_info) offset_x
Definition m_global_parameters.fpp.f90:137
m_global_parameters::z_cc
real(wp), dimension(:), allocatable z_cc
Definition m_global_parameters.fpp.f90:69
m_global_parameters::num_procs
integer num_procs
Number of processors.
Definition m_global_parameters.fpp.f90:36
m_global_parameters::y_output_idx
type(int_bounds_info) y_output_idx
Definition m_global_parameters.fpp.f90:132
m_global_parameters::offset_z
type(int_bounds_info) offset_z
Definition m_global_parameters.fpp.f90:137
m_global_parameters::dx
real(wp), dimension(:), allocatable dx
Cell-width distributions in the x-, y- and z-coordinate directions.
Definition m_global_parameters.fpp.f90:75
m_global_parameters::m_root
integer m_root
Definition m_global_parameters.fpp.f90:44
m_global_parameters::y_output
type(bounds_info) y_output
Definition m_global_parameters.fpp.f90:131
m_helper
Utility routines for bubble model setup, coordinate transforms, array sampling, and special functions...
Definition m_helper.fpp.f90:346
m_helper::s_int_to_str
elemental subroutine, public s_int_to_str(i, res)
Convert an integer to its trimmed string representation.
Definition m_helper.fpp.f90:750
m_mpi_proxy
MPI gather and scatter operations for distributing post-process grid and flow-variable data.
Definition m_mpi_proxy.fpp.f90:7
m_mpi_proxy::s_mpi_defragment_1d_grid_variable
impure subroutine s_mpi_defragment_1d_grid_variable
Collect the sub-domain cell-boundary or cell-center location data from all processors and put back to...
Definition m_mpi_proxy.fpp.f90:482
m_mpi_proxy::s_mpi_defragment_1d_flow_variable
impure subroutine s_mpi_defragment_1d_flow_variable(q_sf, q_root_sf)
Gather the sub-domain flow variable data from all processors and reassemble it for the entire computa...
Definition m_mpi_proxy.fpp.f90:538
m_mpi_proxy::s_mpi_gather_spatial_extents
impure subroutine s_mpi_gather_spatial_extents(spatial_extents)
Gather spatial extents from all ranks for Silo database metadata.
Definition m_mpi_proxy.fpp.f90:390
m_mpi_proxy::s_mpi_gather_data_extents
impure subroutine s_mpi_gather_data_extents(q_sf, data_extents)
Gather the Silo database metadata for the flow variable's extents to boost performance of the multidi...
Definition m_mpi_proxy.fpp.f90:503
m_variables_conversion
Conservative-to-primitive variable conversion, mixture property evaluation, and pressure computation.
Definition m_variables_conversion.fpp.f90:346
m_derived_types::output_context
Output workspace for post_process: flow variable buffers, VisIt extents/offsets, directory paths,...
Definition m_derived_types.fpp.f90:513
m_derived_types::scalar_field
Derived type annexing a scalar field (SF).
Definition m_derived_types.fpp.f90:349