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/simulation/m_data_output.fpp"
2!>
3!! @file
4!! @brief Contains module m_data_output
5
6# 1 "/home/runner/work/MFC/MFC/src/common/include/macros.fpp" 1
7# 1 "/home/runner/work/MFC/MFC/src/common/include/parallel_macros.fpp" 1
8# 1 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp" 1
9# 2 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
10# 3 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
11# 4 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
12# 5 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
13# 6 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
14
15# 8 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
16# 9 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
17# 10 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
18
19# 17 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
20
21# 46 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
22
23# 58 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
24
25# 68 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
26
27# 98 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
28
29# 110 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
30
31# 120 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
32
33# 167 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
34! New line at end of file is required for FYPP
35# 2 "/home/runner/work/MFC/MFC/src/common/include/parallel_macros.fpp" 2
36# 1 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp" 1
37# 1 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp" 1
38# 2 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
39# 3 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
40# 4 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
41# 5 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
42# 6 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
43
44# 8 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
45# 9 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
46# 10 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
47
48# 17 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
49
50# 46 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
51
52# 58 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
53
54# 68 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
55
56# 98 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
57
58# 110 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
59
60# 120 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
61
62# 167 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
63! New line at end of file is required for FYPP
64# 2 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp" 2
65
66# 4 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
67# 5 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
68# 6 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
69# 7 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
70# 8 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
71
72# 20 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
73
74# 43 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
75
76# 48 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
77
78# 53 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
79
80# 58 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
81
82# 63 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
83
84# 68 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
85
86# 81 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
87
88# 86 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
89
90# 91 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
91
92# 96 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
93
94# 101 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
95
96# 106 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
97
98# 111 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
99
100# 116 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
101
102# 121 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
103
104# 126 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
105
106# 156 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
107
108# 197 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
109
110# 211 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
111
112# 236 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
113
114# 247 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
115
116# 249 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
117# 260 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
118
119# 310 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
120
121# 320 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
122
123# 330 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
124
125# 339 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
126
127# 356 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
128
129# 366 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
130
131# 373 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
132
133# 379 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
134
135# 385 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
136
137# 391 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
138
139# 397 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
140
141# 403 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
142! New line at end of file is required for FYPP
143# 3 "/home/runner/work/MFC/MFC/src/common/include/parallel_macros.fpp" 2
144# 1 "/home/runner/work/MFC/MFC/src/common/include/acc_macros.fpp" 1
145# 1 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp" 1
146# 2 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
147# 3 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
148# 4 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
149# 5 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
150# 6 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
151
152# 8 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
153# 9 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
154# 10 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
155
156# 17 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
157
158# 46 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
159
160# 58 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
161
162# 68 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
163
164# 98 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
165
166# 110 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
167
168# 120 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
169
170# 167 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
171! New line at end of file is required for FYPP
172# 2 "/home/runner/work/MFC/MFC/src/common/include/acc_macros.fpp" 2
173
174# 7 "/home/runner/work/MFC/MFC/src/common/include/acc_macros.fpp"
175
176# 17 "/home/runner/work/MFC/MFC/src/common/include/acc_macros.fpp"
177
178# 22 "/home/runner/work/MFC/MFC/src/common/include/acc_macros.fpp"
179
180# 27 "/home/runner/work/MFC/MFC/src/common/include/acc_macros.fpp"
181
182# 32 "/home/runner/work/MFC/MFC/src/common/include/acc_macros.fpp"
183
184# 37 "/home/runner/work/MFC/MFC/src/common/include/acc_macros.fpp"
185
186# 42 "/home/runner/work/MFC/MFC/src/common/include/acc_macros.fpp"
187
188# 47 "/home/runner/work/MFC/MFC/src/common/include/acc_macros.fpp"
189
190# 52 "/home/runner/work/MFC/MFC/src/common/include/acc_macros.fpp"
191
192# 57 "/home/runner/work/MFC/MFC/src/common/include/acc_macros.fpp"
193
194# 62 "/home/runner/work/MFC/MFC/src/common/include/acc_macros.fpp"
195
196# 73 "/home/runner/work/MFC/MFC/src/common/include/acc_macros.fpp"
197
198# 78 "/home/runner/work/MFC/MFC/src/common/include/acc_macros.fpp"
199
200# 83 "/home/runner/work/MFC/MFC/src/common/include/acc_macros.fpp"
201
202# 88 "/home/runner/work/MFC/MFC/src/common/include/acc_macros.fpp"
203
204# 103 "/home/runner/work/MFC/MFC/src/common/include/acc_macros.fpp"
205
206# 131 "/home/runner/work/MFC/MFC/src/common/include/acc_macros.fpp"
207
208# 160 "/home/runner/work/MFC/MFC/src/common/include/acc_macros.fpp"
209
210# 175 "/home/runner/work/MFC/MFC/src/common/include/acc_macros.fpp"
211
212# 193 "/home/runner/work/MFC/MFC/src/common/include/acc_macros.fpp"
213
214# 215 "/home/runner/work/MFC/MFC/src/common/include/acc_macros.fpp"
215
216# 244 "/home/runner/work/MFC/MFC/src/common/include/acc_macros.fpp"
217
218# 259 "/home/runner/work/MFC/MFC/src/common/include/acc_macros.fpp"
219
220# 269 "/home/runner/work/MFC/MFC/src/common/include/acc_macros.fpp"
221
222# 278 "/home/runner/work/MFC/MFC/src/common/include/acc_macros.fpp"
223
224# 294 "/home/runner/work/MFC/MFC/src/common/include/acc_macros.fpp"
225
226# 304 "/home/runner/work/MFC/MFC/src/common/include/acc_macros.fpp"
227
228# 311 "/home/runner/work/MFC/MFC/src/common/include/acc_macros.fpp"
229! New line at end of file is required for FYPP
230# 4 "/home/runner/work/MFC/MFC/src/common/include/parallel_macros.fpp" 2
231
232! GPU parallel region (scalar reductions, maxval/minval)
233# 23 "/home/runner/work/MFC/MFC/src/common/include/parallel_macros.fpp"
234
235! GPU parallel loop over threads (most common GPU macro)
236# 43 "/home/runner/work/MFC/MFC/src/common/include/parallel_macros.fpp"
237
238! Required closing for GPU_PARALLEL_LOOP
239# 55 "/home/runner/work/MFC/MFC/src/common/include/parallel_macros.fpp"
240
241! Mark routine for device compilation
242# 112 "/home/runner/work/MFC/MFC/src/common/include/parallel_macros.fpp"
243
244! Declare device-resident data
245# 130 "/home/runner/work/MFC/MFC/src/common/include/parallel_macros.fpp"
246
247! Inner loop within a GPU parallel region
248# 145 "/home/runner/work/MFC/MFC/src/common/include/parallel_macros.fpp"
249
250! Scoped GPU data region
251# 164 "/home/runner/work/MFC/MFC/src/common/include/parallel_macros.fpp"
252
253! Host code with device pointers (for MPI with GPU buffers)
254# 193 "/home/runner/work/MFC/MFC/src/common/include/parallel_macros.fpp"
255
256! Allocate device memory (unscoped)
257# 207 "/home/runner/work/MFC/MFC/src/common/include/parallel_macros.fpp"
258
259! Free device memory
260# 219 "/home/runner/work/MFC/MFC/src/common/include/parallel_macros.fpp"
261
262! Atomic operation on device
263# 231 "/home/runner/work/MFC/MFC/src/common/include/parallel_macros.fpp"
264
265! End atomic capture block
266# 242 "/home/runner/work/MFC/MFC/src/common/include/parallel_macros.fpp"
267
268! Copy data between host and device
269# 254 "/home/runner/work/MFC/MFC/src/common/include/parallel_macros.fpp"
270
271! Synchronization barrier
272# 266 "/home/runner/work/MFC/MFC/src/common/include/parallel_macros.fpp"
273
274! Import GPU library module (openacc or omp_lib)
275# 275 "/home/runner/work/MFC/MFC/src/common/include/parallel_macros.fpp"
276
277! Emit code only for AMD compiler
278# 282 "/home/runner/work/MFC/MFC/src/common/include/parallel_macros.fpp"
279
280! Emit code for non-Cray compilers
281# 289 "/home/runner/work/MFC/MFC/src/common/include/parallel_macros.fpp"
282
283! Emit code only for Cray compiler
284# 296 "/home/runner/work/MFC/MFC/src/common/include/parallel_macros.fpp"
285
286! Emit code for non-NVIDIA compilers
287# 303 "/home/runner/work/MFC/MFC/src/common/include/parallel_macros.fpp"
288
289# 305 "/home/runner/work/MFC/MFC/src/common/include/parallel_macros.fpp"
290# 306 "/home/runner/work/MFC/MFC/src/common/include/parallel_macros.fpp"
291! New line at end of file is required for FYPP
292# 2 "/home/runner/work/MFC/MFC/src/common/include/macros.fpp" 2
293
294# 14 "/home/runner/work/MFC/MFC/src/common/include/macros.fpp"
295
296! Caution: This macro requires the use of a binding script to set CUDA_VISIBLE_DEVICES, such that we have one GPU device per MPI
297! rank. That's because for both cudaMemAdvise (preferred location) and cudaMemPrefetchAsync we use location = device_id = 0. For an
298! example see misc/nvidia_uvm/bind.sh.
299# 55 "/home/runner/work/MFC/MFC/src/common/include/macros.fpp"
300
301! Allocate and create GPU device memory
302# 75 "/home/runner/work/MFC/MFC/src/common/include/macros.fpp"
303
304! Free GPU device memory and deallocate
305# 83 "/home/runner/work/MFC/MFC/src/common/include/macros.fpp"
306
307! Cray-specific GPU pointer setup for vector fields
308# 107 "/home/runner/work/MFC/MFC/src/common/include/macros.fpp"
309
310! Cray-specific GPU pointer setup for scalar fields
311# 123 "/home/runner/work/MFC/MFC/src/common/include/macros.fpp"
312
313! Cray-specific GPU pointer setup for acoustic source spatials
314# 148 "/home/runner/work/MFC/MFC/src/common/include/macros.fpp"
315
316# 154 "/home/runner/work/MFC/MFC/src/common/include/macros.fpp"
317
318# 161 "/home/runner/work/MFC/MFC/src/common/include/macros.fpp"
319! New line at end of file is required for FYPP
320# 6 "/home/runner/work/MFC/MFC/src/simulation/m_data_output.fpp" 2
321# 1 "/home/runner/work/MFC/MFC/src/common/include/case.fpp" 1
322! This file exists so that Fypp can be run without generating case.fpp files for
323! each target. This is useful when generating documentation, for example. This
324! should also let MFC be built with CMake directly, without invoking mfc.sh.
325
326! For pre-process.
327# 8 "/home/runner/work/MFC/MFC/src/common/include/case.fpp"
328
329! For moving immersed boundaries in simulation
330# 12 "/home/runner/work/MFC/MFC/src/common/include/case.fpp"
331# 7 "/home/runner/work/MFC/MFC/src/simulation/m_data_output.fpp" 2
332
333!> @brief Writes solution data, run-time stability diagnostics (ICFL, VCFL, CCFL, Rc), and probe/center-of-mass files
335
338 use m_mpi_proxy
341 use m_helper
343 use m_sim_helpers
345 use m_ibm
348
349 implicit none
350
351 private
356
357 real(wp), public, allocatable, dimension(:,:) :: c_mass
358
359# 33 "/home/runner/work/MFC/MFC/src/simulation/m_data_output.fpp"
360#if defined(MFC_OpenACC)
361# 33 "/home/runner/work/MFC/MFC/src/simulation/m_data_output.fpp"
362!$acc declare create(c_mass)
363# 33 "/home/runner/work/MFC/MFC/src/simulation/m_data_output.fpp"
364#elif defined(MFC_OpenMP)
365# 33 "/home/runner/work/MFC/MFC/src/simulation/m_data_output.fpp"
366!$omp declare target (c_mass)
367# 33 "/home/runner/work/MFC/MFC/src/simulation/m_data_output.fpp"
368#endif
369
370 !> @name ICFL, VCFL, CCFL, and Rc stability criteria extrema over all the time-steps
371 !> @{
372 real(wp) :: icfl_max !< ICFL criterion maximum
373 real(wp) :: vcfl_max !< VCFL criterion maximum
374 real(wp) :: ccfl_max !< CCFL criterion maximum
375 real(wp) :: rc_min !< Rc criterion maximum
376 !> @}
377
378 type(scalar_field), allocatable, dimension(:) :: q_cons_temp_ds
379
380contains
381
382 !> Write data files. Dispatch subroutine that replaces procedure pointer.
383 impure subroutine s_write_data_files(q_cons_vf, q_T_sf, q_prim_vf, t_step, bc_type, beta)
384
385 type(scalar_field), dimension(sys_size), intent(inout) :: q_cons_vf
386 type(scalar_field), intent(inout) :: q_t_sf
387 type(scalar_field), dimension(sys_size), intent(inout) :: q_prim_vf
388 integer, intent(in) :: t_step
389 type(scalar_field), intent(inout), optional :: beta
390 type(integer_field), dimension(1:num_dims,-1:1), intent(in) :: bc_type
391
392 if (.not. parallel_io) then
393 call s_write_serial_data_files(q_cons_vf, q_t_sf, q_prim_vf, t_step, bc_type, beta)
394 else
395 call s_write_parallel_data_files(q_cons_vf, t_step, bc_type, beta)
396 end if
397
398 end subroutine s_write_data_files
399
400 !> Open the run-time information file and write the stability criteria table header
402
403 character(LEN=name_len), parameter :: file_name = 'run_time.inf' !< Name of the run-time information file
404 character(LEN=path_len + name_len) :: file_path !< Relative path to a file in the case directory
405 character(LEN=8) :: file_date !< Creation date of the run-time information file
406
407 file_path = trim(case_dir) // '/' // trim(file_name)
408
409 open (3, file=trim(file_path), form='formatted', status='replace')
410
411 write (3, '(A)') 'Description: Stability information at ' // 'each time-step of the simulation. This'
412 write (3, '(13X,A)') 'data is composed of the inviscid ' // 'Courant-Friedrichs-Lewy (ICFL)'
413 write (3, '(13X,A)') 'number, the viscous CFL (VCFL) number, ' // 'the capillary CFL (CCFL)'
414 write (3, '(13X,A)') 'number and the cell Reynolds (Rc) ' // 'number. Please note that only'
415 write (3, '(13X,A)') 'those stability conditions pertinent ' // 'to the physics included in'
416 write (3, '(13X,A)') 'the current computation are displayed.'
417 if (hypoelasticity) then
418 write (3, '(13X,A)') 'NOTE: the reported ICFL uses the acoustic ' // 'sound speed only; it may'
419 write (3, '(13X,A)') 'underestimate the elastic characteristic ' // 'speeds.'
420 end if
421
422 call date_and_time(date=file_date)
423
424 write (3, '(A)') 'Date: ' // file_date(5:6) // '/' // file_date(7:8) // '/' // file_date(3:4)
425
426 write (3, '(A)') ''; write (3, '(A)') ''
427
428 write (3, '(13X,A9,13X,A10,13X,A10,13X,A10)', advance="no") trim('Time-step'), trim('dt'), trim('Time'), trim('ICFL Max')
429
430 if (surface_tension) then
431 write (3, '(13X,A10)', advance="no") trim('CCFL Max')
432 end if
433
434 if (viscous) then
435 write (3, '(13X,A10,13X,A16)', advance="no") trim('VCFL Max'), trim('Rc Min')
436 end if
437
438 if (bubbles_lagrange) then
439 write (3, '(13X,A10)', advance="no") trim('N Bubbles')
440 end if
441
442 write (3, *) ! new line
443
445
446 !> Open center-of-mass data files for writing
447 impure subroutine s_open_com_files()
448
449 character(len=path_len + 3*name_len) :: file_path !< Relative path to the CoM file in the case directory
450 integer :: i !< Generic loop iterator
451
452 do i = 1, num_fluids
453 write (file_path, '(A,I0,A)') '/fluid', i, '_com.dat'
454 file_path = trim(case_dir) // trim(file_path)
455 open (i + 120, file=trim(file_path), form='formatted', position='append', status='unknown')
456 if (n == 0) then
457 write (i + 120, '(A)') ' Non-Dimensional Time ' // ' Total Mass ' // ' x-loc ' // ' Total Volume '
458 else if (p == 0) then
459 write (i + 120, &
460 & '(A)') ' Non-Dimensional Time ' // ' Total Mass ' // ' x-loc ' // ' y-loc ' &
461 & // ' Total Volume '
462 else
463 write (i + 120, &
464 & '(A)') ' Non-Dimensional Time ' // ' Total Mass ' // ' x-loc ' // ' y-loc ' // ' z-loc ' &
465 & // ' Total Volume '
466 end if
467 end do
468
469 end subroutine s_open_com_files
470
471 !> Open flow probe data files for writing
472 impure subroutine s_open_probe_files
473
474 character(LEN=path_len + 3*name_len) :: file_path !< Relative path to the probe data file in the case directory
475 integer :: i !< Generic loop iterator
476 logical :: file_exist
477
478 do i = 1, num_probes
479 write (file_path, '(A,I0,A)') '/D/probe', i, '_prim.dat'
480 file_path = trim(case_dir) // trim(file_path)
481
482 inquire (file=trim(file_path), exist=file_exist)
483
484 if (file_exist) then
485 open (i + 30, file=trim(file_path), form='formatted', status='old', position='append')
486 else
487 open (i + 30, file=trim(file_path), form='formatted', status='unknown')
488 end if
489 end do
490
491 end subroutine s_open_probe_files
492
493 !> Write stability criteria extrema to the run-time information file at the given time step
494 impure subroutine s_write_run_time_information(q_prim_vf, t_step)
495
496 type(scalar_field), dimension(sys_size), intent(in) :: q_prim_vf
497 integer, intent(in) :: t_step
498 real(wp) :: rho !< Cell-avg. density
499
500# 169 "/home/runner/work/MFC/MFC/src/simulation/m_data_output.fpp"
501 real(wp), dimension(num_fluids) :: alpha !< Cell-avg. volume fraction
502 real(wp), dimension(num_vels) :: vel !< Cell-avg. velocity
503# 172 "/home/runner/work/MFC/MFC/src/simulation/m_data_output.fpp"
504 real(wp) :: vel_sum !< Cell-avg. velocity sum
505 real(wp) :: pres !< Cell-avg. pressure
506 real(wp) :: gamma !< Cell-avg. sp. heat ratio
507 real(wp) :: pi_inf !< Cell-avg. liquid stiffness function
508 real(wp) :: qv !< Cell-avg. internal energy reference value
509 real(wp) :: c !< Cell-avg. sound speed
510 real(wp) :: h !< Cell-avg. enthalpy
511 real(wp), dimension(2) :: re !< Cell-avg. Reynolds numbers
512 integer :: j, k, l
513 real(wp) :: icfl_max_loc, icfl_max_glb !< ICFL stability extrema on local and global grids
514 real(wp) :: vcfl_max_loc, vcfl_max_glb !< VCFL stability extrema on local and global grids
515 real(wp) :: ccfl_max_loc, ccfl_max_glb !< CCFL stability extrema on local and global grids
516 real(wp) :: rc_min_loc, rc_min_glb !< Rc stability extrema on local and global grids
517 real(wp) :: icfl, vcfl, ccfl, rc
518 integer :: fl !< Fluid loop iterator
519
520 icfl_max_loc = 0._wp
521 vcfl_max_loc = 0._wp
522 ccfl_max_loc = 0._wp
523 rc_min_loc = huge(1.0_wp)
524 ! Computing Stability Criteria at Current Time-step
525
526# 193 "/home/runner/work/MFC/MFC/src/simulation/m_data_output.fpp"
527
528# 193 "/home/runner/work/MFC/MFC/src/simulation/m_data_output.fpp"
529#if defined(MFC_OpenACC)
530# 193 "/home/runner/work/MFC/MFC/src/simulation/m_data_output.fpp"
531!$acc parallel loop collapse(3) gang vector default(present) private(j, k, l, vel, alpha, Re, rho, vel_sum, pres, gamma, pi_inf, c, H, qv, icfl, vcfl, Rc, ccfl, fl) &
532# 193 "/home/runner/work/MFC/MFC/src/simulation/m_data_output.fpp"
533!$acc& reduction(max:icfl_max_loc, vcfl_max_loc, ccfl_max_loc) reduction(min:Rc_min_loc)
534# 193 "/home/runner/work/MFC/MFC/src/simulation/m_data_output.fpp"
535#elif defined(MFC_OpenMP)
536# 193 "/home/runner/work/MFC/MFC/src/simulation/m_data_output.fpp"
537
538# 193 "/home/runner/work/MFC/MFC/src/simulation/m_data_output.fpp"
539
540# 193 "/home/runner/work/MFC/MFC/src/simulation/m_data_output.fpp"
541
542# 193 "/home/runner/work/MFC/MFC/src/simulation/m_data_output.fpp"
543!$omp target teams loop defaultmap(firstprivate:scalar) bind(teams,parallel) collapse(3) defaultmap(tofrom:aggregate) defaultmap(tofrom:allocatable) defaultmap(tofrom:pointer) &
544# 193 "/home/runner/work/MFC/MFC/src/simulation/m_data_output.fpp"
545!$omp& private(j, k, l, vel, alpha, Re, rho, vel_sum, pres, gamma, pi_inf, c, H, qv, icfl, vcfl, Rc, ccfl, fl) reduction(max:icfl_max_loc, vcfl_max_loc, ccfl_max_loc) reduction(min:Rc_min_loc)
546# 193 "/home/runner/work/MFC/MFC/src/simulation/m_data_output.fpp"
547#endif
548# 196 "/home/runner/work/MFC/MFC/src/simulation/m_data_output.fpp"
549 do l = 0, p
550 do k = 0, n
551 do j = 0, m
552 call s_compute_enthalpy(q_prim_vf, pres, rho, gamma, pi_inf, re, h, alpha, vel, vel_sum, qv, j, k, l)
553
554 call s_compute_speed_of_sound(pres, rho, gamma, pi_inf, h, alpha, vel_sum, 0._wp, c, qv)
555
556 if (any_non_newtonian) then
557 re(1) = 0._wp
558 do fl = 1, num_fluids
559 if (is_non_newtonian(fl)) then
560 re(1) = re(1) + alpha(fl)*hb_mu_max(fl)
561 else
562 re(1) = re(1) + alpha(fl)*fluid_inv_re(fl)
563 end if
564 end do
565 re(1) = 1._wp/max(re(1), sgm_eps)
566 end if
567
568 call s_compute_stability_from_dt(vel, c, rho, re, j, k, l, icfl, vcfl, rc, ccfl)
569
570 icfl_max_loc = max(icfl_max_loc, icfl)
571 vcfl_max_loc = max(vcfl_max_loc, merge(vcfl, 0.0_wp, viscous))
572 ccfl_max_loc = max(ccfl_max_loc, merge(ccfl, 0.0_wp, surface_tension))
573 rc_min_loc = min(rc_min_loc, merge(rc, huge(1.0_wp), viscous))
574 end do
575 end do
576 end do
577
578# 224 "/home/runner/work/MFC/MFC/src/simulation/m_data_output.fpp"
579#if defined(MFC_OpenACC)
580# 224 "/home/runner/work/MFC/MFC/src/simulation/m_data_output.fpp"
581!$acc end parallel loop
582# 224 "/home/runner/work/MFC/MFC/src/simulation/m_data_output.fpp"
583#elif defined(MFC_OpenMP)
584# 224 "/home/runner/work/MFC/MFC/src/simulation/m_data_output.fpp"
585
586# 224 "/home/runner/work/MFC/MFC/src/simulation/m_data_output.fpp"
587!$omp end target teams loop
588# 224 "/home/runner/work/MFC/MFC/src/simulation/m_data_output.fpp"
589#endif
590 ! end: Computing Stability Criteria at Current Time-step
591
592 if (num_procs > 1) then
593 call s_mpi_reduce_stability_criteria_extrema(icfl_max_loc, vcfl_max_loc, rc_min_loc, n_el_bubs_loc, icfl_max_glb, &
594 & vcfl_max_glb, rc_min_glb, n_el_bubs_glb, ccfl_max_loc, ccfl_max_glb)
595 else
596 icfl_max_glb = icfl_max_loc
597 if (viscous) vcfl_max_glb = vcfl_max_loc
598 if (viscous) rc_min_glb = rc_min_loc
599 if (surface_tension) ccfl_max_glb = ccfl_max_loc
600 if (bubbles_lagrange) n_el_bubs_glb = n_el_bubs_loc
601 end if
602
603 if (icfl_max_glb > icfl_max) icfl_max = icfl_max_glb
604
605 if (surface_tension) then
606 if (ccfl_max_glb > ccfl_max) ccfl_max = ccfl_max_glb
607 end if
608
609 if (viscous) then
610 if (vcfl_max_glb > vcfl_max) vcfl_max = vcfl_max_glb
611 if (rc_min_glb < rc_min) rc_min = rc_min_glb
612 end if
613
614 if (proc_rank == 0) then
615 write (3, '(13X,I9,13X,F10.6,13X,F10.6,13X,F10.6)', advance="no") t_step, dt, mytime, icfl_max_glb
616
617 if (surface_tension) then
618 write (3, '(13X,F10.6)', advance="no") ccfl_max_glb
619 end if
620
621 if (viscous) then
622 write (3, '(13X,F10.6,13X,ES16.6)', advance="no") vcfl_max_glb, rc_min_glb
623 end if
624
625 if (bubbles_lagrange) then
626 write (3, '(13X,I10)', advance="no") n_el_bubs_glb
627 end if
628
629 write (3, *) ! new line
630
631 if (.not. f_approx_equal(icfl_max_glb, icfl_max_glb)) then
632 call s_mpi_abort('ICFL is NaN. Exiting.')
633 else if (icfl_max_glb > 1._wp) then
634 print *, 'icfl', icfl_max_glb
635 call s_mpi_abort('ICFL is greater than 1.0. Exiting.')
636 end if
637
638 if (viscous) then
639 if (.not. f_approx_equal(vcfl_max_glb, vcfl_max_glb)) then
640 call s_mpi_abort('VCFL is NaN. Exiting.')
641 else if (vcfl_max_glb > 1._wp) then
642 print *, 'vcfl', vcfl_max_glb
643 call s_mpi_abort('VCFL is greater than 1.0. Exiting.')
644 end if
645 end if
646
647 if (bubbles_lagrange) then
648 if (n_el_bubs_glb == 0) then
649 call s_mpi_abort('No Lagrangian bubbles remain in the domain. Exiting.')
650 end if
651 end if
652 end if
653
654 call s_mpi_barrier()
655
656 end subroutine s_write_run_time_information
657
658 !> Write grid and conservative variable data files in serial format
659 impure subroutine s_write_serial_data_files(q_cons_vf, q_T_sf, q_prim_vf, t_step, bc_type, beta)
660
661 type(scalar_field), dimension(sys_size), intent(inout) :: q_cons_vf
662 type(scalar_field), intent(inout) :: q_t_sf
663 type(scalar_field), dimension(sys_size), intent(inout) :: q_prim_vf
664 integer, intent(in) :: t_step
665 type(scalar_field), intent(inout), optional :: beta
666 type(integer_field), dimension(1:num_dims,-1:1), intent(in) :: bc_type
667 character(LEN=path_len + 2*name_len) :: t_step_dir !< Relative path to the current time-step directory
668 character(LEN=path_len + 3*name_len) :: file_path !< Relative path to the grid and conservative variables data files
669 logical :: file_exist !< Logical used to check existence of current time-step directory
670 character(LEN=15) :: fmt
671 integer :: i, j, k, l, r
672 real(wp) :: gamma, lit_gamma, pi_inf, qv !< Temporary EOS params
673
674 write (t_step_dir, '(A,I0,A,I0)') trim(case_dir) // '/p_all'
675 write (t_step_dir, '(a,i0,a,i0)') trim(case_dir) // '/p_all/p', proc_rank, '/', t_step
676
677 file_path = trim(t_step_dir) // '/.'
678 call my_inquire(file_path, file_exist)
679 if (file_exist) call s_delete_directory(trim(t_step_dir))
680 call s_create_directory(trim(t_step_dir))
681
682 file_path = trim(t_step_dir) // '/x_cb.dat'
683
684 open (2, file=trim(file_path), form='unformatted', status='new')
685 write (2) x_cb(-1:m); close (2)
686
687 if (n > 0) then
688 file_path = trim(t_step_dir) // '/y_cb.dat'
689
690 open (2, file=trim(file_path), form='unformatted', status='new')
691 write (2) y_cb(-1:n); close (2)
692
693 if (p > 0) then
694 file_path = trim(t_step_dir) // '/z_cb.dat'
695
696 open (2, file=trim(file_path), form='unformatted', status='new')
697 write (2) z_cb(-1:p); close (2)
698 end if
699 end if
700
701 do i = 1, sys_size
702 write (file_path, '(A,I0,A)') trim(t_step_dir) // '/q_cons_vf', i, '.dat'
703
704 open (2, file=trim(file_path), form='unformatted', status='new')
705
706 write (2) q_cons_vf(i)%sf(0:m,0:n,0:p); close (2)
707 end do
708
709 ! Lagrangian beta (void fraction) written as q_cons_vf(sys_size+1) to match the parallel I/O path and allow post_process to
710 ! read it.
711 if (bubbles_lagrange) then
712 write (file_path, '(A,I0,A)') trim(t_step_dir) // '/q_cons_vf', sys_size + 1, '.dat'
713
714 open (2, file=trim(file_path), form='unformatted', status='new')
715
716 write (2) beta%sf(0:m,0:n,0:p); close (2)
717 end if
718
719 if (qbmm .and. .not. polytropic) then
720 do i = 1, nb
721 do r = 1, nnode
722 write (file_path, '(A,I0,A)') trim(t_step_dir) // '/pb', sys_size + (i - 1)*nnode + r, '.dat'
723
724 open (2, file=trim(file_path), form='unformatted', status='new')
725
726 write (2) pb_ts(1)%sf(0:m,0:n,0:p,r, i); close (2)
727 end do
728 end do
729
730 do i = 1, nb
731 do r = 1, nnode
732 write (file_path, '(A,I0,A)') trim(t_step_dir) // '/mv', sys_size + (i - 1)*nnode + r, '.dat'
733
734 open (2, file=trim(file_path), form='unformatted', status='new')
735
736 write (2) mv_ts(1)%sf(0:m,0:n,0:p,r, i); close (2)
737 end do
738 end do
739 end if
740
741 ! Writing the IB markers
742 if (ib) then
743 call s_write_serial_ib_data(t_step)
744 end if
745
746 gamma = gammas(1)
747 lit_gamma = gs_min(1)
748 pi_inf = pi_infs(1)
749 qv = qvs(1)
750
751 if (precision == precision_single) then
752 fmt = "(2F30.3)"
753 else
754 fmt = "(2F40.14)"
755 end if
756
757 write (t_step_dir, '(A,I0,A,I0)') trim(case_dir) // '/D'
758 file_path = trim(t_step_dir) // '/.'
759
760 inquire (file=trim(file_path), exist=file_exist)
761
762 if (.not. file_exist) call s_create_directory(trim(t_step_dir))
763
764 if ((prim_vars_wrt .or. (n == 0 .and. p == 0)) .and. (.not. igr)) then
766 do i = 1, sys_size
767
768# 402 "/home/runner/work/MFC/MFC/src/simulation/m_data_output.fpp"
769#if defined(MFC_OpenACC)
770# 402 "/home/runner/work/MFC/MFC/src/simulation/m_data_output.fpp"
771!$acc update host(q_prim_vf(i)%sf(:, :, :))
772# 402 "/home/runner/work/MFC/MFC/src/simulation/m_data_output.fpp"
773#elif defined(MFC_OpenMP)
774# 402 "/home/runner/work/MFC/MFC/src/simulation/m_data_output.fpp"
775!$omp target update from(q_prim_vf(i)%sf(:, :, :))
776# 402 "/home/runner/work/MFC/MFC/src/simulation/m_data_output.fpp"
777#endif
778 end do
779 ! q_prim_vf(eqn_idx%bub%beg) stores the value of nb needed in riemann solvers, so replace with true primitive value
780 ! (=1._wp)
781 if (qbmm) then
782 q_prim_vf(eqn_idx%bub%beg)%sf = 1._wp
783 end if
784 end if
785
786 if (n == 0 .and. p == 0) then
787 if (model_eqns == model_eqns_5eq .and. (.not. igr)) then
788 do i = 1, sys_size
789 write (file_path, '(A,I0,A,I2.2,A,I6.6,A)') trim(t_step_dir) // '/prim.', i, '.', proc_rank, '.', t_step, '.dat'
790
791 open (2, file=trim(file_path))
792 do j = 0, m
793 ! todo: revisit change here
794 if (((i >= eqn_idx%adv%beg) .and. (i <= eqn_idx%adv%end))) then
795 write (2, fmt) x_cb(j), q_cons_vf(i)%sf(j, 0, 0)
796 else
797 write (2, fmt) x_cb(j), q_prim_vf(i)%sf(j, 0, 0)
798 end if
799 end do
800 close (2)
801 end do
802 end if
803
804 do i = 1, sys_size
805 write (file_path, '(A,I0,A,I2.2,A,I6.6,A)') trim(t_step_dir) // '/cons.', i, '.', proc_rank, '.', t_step, '.dat'
806
807 open (2, file=trim(file_path))
808 do j = 0, m
809 write (2, fmt) x_cb(j), q_cons_vf(i)%sf(j, 0, 0)
810 end do
811 close (2)
812 end do
813
814 if (qbmm .and. .not. polytropic) then
815 do i = 1, nb
816 do r = 1, nnode
817 write (file_path, '(A,I0,A,I0,A,I2.2,A,I6.6,A)') trim(t_step_dir) // '/pres.', i, '.', r, '.', proc_rank, &
818 & '.', t_step, '.dat'
819
820 open (2, file=trim(file_path))
821 do j = 0, m
822 write (2, fmt) x_cb(j), pb_ts(1)%sf(j, 0, 0, r, i)
823 end do
824 close (2)
825 end do
826 end do
827 do i = 1, nb
828 do r = 1, nnode
829 write (file_path, '(A,I0,A,I0,A,I2.2,A,I6.6,A)') trim(t_step_dir) // '/mv.', i, '.', r, '.', proc_rank, &
830 & '.', t_step, '.dat'
831
832 open (2, file=trim(file_path))
833 do j = 0, m
834 write (2, fmt) x_cb(j), mv_ts(1)%sf(j, 0, 0, r, i)
835 end do
836 close (2)
837 end do
838 end do
839 end if
840 end if
841
842 if (precision == precision_single) then
843 fmt = "(3F30.7)"
844 else
845 fmt = "(3F40.14)"
846 end if
847
848 if ((n > 0) .and. (p == 0)) then
849 do i = 1, sys_size
850 write (file_path, '(A,I0,A,I2.2,A,I6.6,A)') trim(t_step_dir) // '/cons.', i, '.', proc_rank, '.', t_step, '.dat'
851 open (2, file=trim(file_path))
852 do j = 0, m
853 do k = 0, n
854 write (2, fmt) x_cb(j), y_cb(k), q_cons_vf(i)%sf(j, k, 0)
855 end do
856 write (2, *)
857 end do
858 close (2)
859 end do
860
861 if (present(beta)) then
862 write (file_path, '(A,I0,A,I2.2,A,I6.6,A)') trim(t_step_dir) // '/beta.', i, '.', proc_rank, '.', t_step, '.dat'
863 open (2, file=trim(file_path))
864 do j = 0, m
865 do k = 0, n
866 write (2, fmt) x_cb(j), y_cb(k), beta%sf(j, k, 0)
867 end do
868 write (2, *)
869 end do
870 close (2)
871 end if
872
873 if (qbmm .and. .not. polytropic) then
874 do i = 1, nb
875 do r = 1, nnode
876 write (file_path, '(A,I0,A,I0,A,I2.2,A,I6.6,A)') trim(t_step_dir) // '/pres.', i, '.', r, '.', proc_rank, &
877 & '.', t_step, '.dat'
878
879 open (2, file=trim(file_path))
880 do j = 0, m
881 do k = 0, n
882 write (2, fmt) x_cb(j), y_cb(k), pb_ts(1)%sf(j, k, 0, r, i)
883 end do
884 end do
885 close (2)
886 end do
887 end do
888 do i = 1, nb
889 do r = 1, nnode
890 write (file_path, '(A,I0,A,I0,A,I2.2,A,I6.6,A)') trim(t_step_dir) // '/mv.', i, '.', r, '.', proc_rank, &
891 & '.', t_step, '.dat'
892
893 open (2, file=trim(file_path))
894 do j = 0, m
895 do k = 0, n
896 write (2, fmt) x_cb(j), y_cb(k), mv_ts(1)%sf(j, k, 0, r, i)
897 end do
898 end do
899 close (2)
900 end do
901 end do
902 end if
903
904 if (prim_vars_wrt .and. (.not. igr)) then
905 do i = 1, sys_size
906 write (file_path, '(A,I0,A,I2.2,A,I6.6,A)') trim(t_step_dir) // '/prim.', i, '.', proc_rank, '.', t_step, '.dat'
907
908 open (2, file=trim(file_path))
909
910 do j = 0, m
911 do k = 0, n
912 if (((i >= eqn_idx%cont%beg) .and. (i <= eqn_idx%cont%end)) .or. ((i >= eqn_idx%adv%beg) &
913 & .and. (i <= eqn_idx%adv%end))) then
914 write (2, fmt) x_cb(j), y_cb(k), q_cons_vf(i)%sf(j, k, 0)
915 else
916 write (2, fmt) x_cb(j), y_cb(k), q_prim_vf(i)%sf(j, k, 0)
917 end if
918 end do
919 write (2, *)
920 end do
921 close (2)
922 end do
923 end if
924 end if
925
926 if (precision == precision_single) then
927 fmt = "(4F30.7)"
928 else
929 fmt = "(4F40.14)"
930 end if
931
932 if (p > 0) then
933 do i = 1, sys_size
934 write (file_path, '(A,I0,A,I2.2,A,I6.6,A)') trim(t_step_dir) // '/cons.', i, '.', proc_rank, '.', t_step, '.dat'
935 open (2, file=trim(file_path))
936 do j = 0, m
937 do k = 0, n
938 do l = 0, p
939 write (2, fmt) x_cb(j), y_cb(k), z_cb(l), q_cons_vf(i)%sf(j, k, l)
940 end do
941 write (2, *)
942 end do
943 write (2, *)
944 end do
945 close (2)
946 end do
947
948 if (present(beta)) then
949 write (file_path, '(A,I0,A,I2.2,A,I6.6,A)') trim(t_step_dir) // '/beta.', i, '.', proc_rank, '.', t_step, '.dat'
950 open (2, file=trim(file_path))
951 do j = 0, m
952 do k = 0, n
953 do l = 0, p
954 write (2, fmt) x_cb(j), y_cb(k), z_cb(l), beta%sf(j, k, l)
955 end do
956 write (2, *)
957 end do
958 write (2, *)
959 end do
960 close (2)
961 end if
962
963 if (qbmm .and. .not. polytropic) then
964 do i = 1, nb
965 do r = 1, nnode
966 write (file_path, '(A,I0,A,I0,A,I2.2,A,I6.6,A)') trim(t_step_dir) // '/pres.', i, '.', r, '.', proc_rank, &
967 & '.', t_step, '.dat'
968
969 open (2, file=trim(file_path))
970 do j = 0, m
971 do k = 0, n
972 do l = 0, p
973 write (2, fmt) x_cb(j), y_cb(k), z_cb(l), pb_ts(1)%sf(j, k, l, r, i)
974 end do
975 end do
976 end do
977 close (2)
978 end do
979 end do
980 do i = 1, nb
981 do r = 1, nnode
982 write (file_path, '(A,I0,A,I0,A,I2.2,A,I6.6,A)') trim(t_step_dir) // '/mv.', i, '.', r, '.', proc_rank, &
983 & '.', t_step, '.dat'
984
985 open (2, file=trim(file_path))
986 do j = 0, m
987 do k = 0, n
988 do l = 0, p
989 write (2, fmt) x_cb(j), y_cb(k), z_cb(l), mv_ts(1)%sf(j, k, l, r, i)
990 end do
991 end do
992 end do
993 close (2)
994 end do
995 end do
996 end if
997
998 if (prim_vars_wrt .and. (.not. igr)) then
999 do i = 1, sys_size
1000 write (file_path, '(A,I0,A,I2.2,A,I6.6,A)') trim(t_step_dir) // '/prim.', i, '.', proc_rank, '.', t_step, '.dat'
1001
1002 open (2, file=trim(file_path))
1003
1004 do j = 0, m
1005 do k = 0, n
1006 do l = 0, p
1007 if (((i >= eqn_idx%cont%beg) .and. (i <= eqn_idx%cont%end)) .or. ((i >= eqn_idx%adv%beg) &
1008 & .and. (i <= eqn_idx%adv%end)) .or. ((i >= eqn_idx%species%beg) &
1009 & .and. (i <= eqn_idx%species%end))) then
1010 write (2, fmt) x_cb(j), y_cb(k), z_cb(l), q_cons_vf(i)%sf(j, k, l)
1011 else
1012 write (2, fmt) x_cb(j), y_cb(k), z_cb(l), q_prim_vf(i)%sf(j, k, l)
1013 end if
1014 end do
1015 write (2, *)
1016 end do
1017 write (2, *)
1018 end do
1019 close (2)
1020 end do
1021 end if
1022 end if
1023
1024 end subroutine s_write_serial_data_files
1025
1026 !> Write grid and conservative variable data files in parallel via MPI I/O
1027 impure subroutine s_write_parallel_data_files(q_cons_vf, t_step, bc_type, beta, q_T_sf)
1028
1029 type(scalar_field), dimension(sys_size), intent(inout) :: q_cons_vf
1030 integer, intent(in) :: t_step
1031 type(scalar_field), intent(inout), optional :: beta
1032 type(integer_field), dimension(1:num_dims,-1:1), intent(in) :: bc_type
1033 type(scalar_field), intent(inout), optional :: q_t_sf
1034
1035#ifdef MFC_MPI
1036 integer :: ifile, ierr, data_size
1037 integer, dimension(MPI_STATUS_SIZE) :: status
1038 integer(kind=MPI_OFFSET_kind) :: disp
1039 integer(kind=MPI_OFFSET_kind) :: m_mok, n_mok, p_mok
1040 integer(kind=MPI_OFFSET_kind) :: wp_mok, var_mok, str_mok
1041 integer(kind=MPI_OFFSET_kind) :: nvars_mok
1042 integer(kind=MPI_OFFSET_kind) :: mok
1043 character(LEN=path_len + 2*name_len) :: file_loc
1044 logical :: file_exist, dir_check
1045 character(len=10) :: t_step_string
1046 integer :: i !< Generic loop iterator
1047 integer :: alt_sys !< Altered system size for the lagrangian subgrid bubble model
1048 ! Down sampling variables
1049 integer :: m_ds, n_ds, p_ds
1050 integer :: m_glb_ds, n_glb_ds, p_glb_ds
1051 integer :: m_glb_save, n_glb_save, p_glb_save !< Global save size
1052
1053 if (down_sample) then
1054 call s_downsample_data(q_cons_vf, q_cons_temp_ds, m_ds, n_ds, p_ds, m_glb_ds, n_glb_ds, p_glb_ds)
1055 end if
1056
1057 if (present(beta)) then
1058 alt_sys = sys_size + 1
1059 else
1060 alt_sys = sys_size
1061 end if
1062
1063 if (file_per_process) then
1064 call s_int_to_str(t_step, t_step_string)
1065
1066 if (down_sample) then
1067 call s_initialize_mpi_data_ds(m_ds, n_ds, p_ds)
1068 else
1069 if (ib) then
1070 call s_initialize_mpi_data(q_cons_vf, ib_markers=ib_markers, ib_mpi_data=mpi_io_ib_data, qbmm_pb=pb_ts(1), &
1071 & qbmm_mv=mv_ts(1))
1072 else
1073 call s_initialize_mpi_data(q_cons_vf, qbmm_pb=pb_ts(1), qbmm_mv=mv_ts(1))
1074 end if
1075 end if
1076
1077 if (proc_rank == 0) then
1078 file_loc = trim(case_dir) // '/restart_data/lustre_' // trim(t_step_string)
1079 call my_inquire(file_loc, dir_check)
1080 if (dir_check .neqv. .true.) then
1081 call s_create_directory(trim(file_loc))
1082 end if
1083 call s_create_directory(trim(file_loc))
1084 end if
1085 call s_mpi_barrier()
1087
1088 call s_initialize_mpi_data(q_cons_vf, qbmm_pb=pb_ts(1), qbmm_mv=mv_ts(1))
1089
1090 write (file_loc, '(I0,A,i7.7,A)') t_step, '_', proc_rank, '.dat'
1091 file_loc = trim(case_dir) // '/restart_data/lustre_' // trim(t_step_string) // trim(mpiiofs) // trim(file_loc)
1092 inquire (file=trim(file_loc), exist=file_exist)
1093 if (file_exist .and. proc_rank == 0) then
1094 call mpi_file_delete(file_loc, mpi_info_int, ierr)
1095 end if
1096 call mpi_file_open(mpi_comm_self, file_loc, ior(mpi_mode_wronly, mpi_mode_create), mpi_info_int, ifile, ierr)
1097
1098 if (down_sample) then
1099 data_size = (m_ds + 3)*(n_ds + 3)*(p_ds + 3)
1100 m_glb_save = m_glb_ds + 1
1101 n_glb_save = n_glb_ds + 1
1102 p_glb_save = p_glb_ds + 1
1103 else
1104 data_size = (m + 1)*(n + 1)*(p + 1)
1105 m_glb_save = m_glb + 1
1106 n_glb_save = n_glb + 1
1107 p_glb_save = p_glb + 1
1108 end if
1109
1110 m_mok = int(m_glb_save + 1, mpi_offset_kind)
1111 n_mok = int(n_glb_save + 1, mpi_offset_kind)
1112 p_mok = int(p_glb_save + 1, mpi_offset_kind)
1113 wp_mok = int(storage_size(0._stp)/8, mpi_offset_kind)
1114 mok = int(1._wp, mpi_offset_kind)
1115 str_mok = int(name_len, mpi_offset_kind)
1116 nvars_mok = int(sys_size, mpi_offset_kind)
1117
1118 if (bubbles_euler) then
1119 do i = 1, sys_size
1120 var_mok = int(i, mpi_offset_kind)
1121
1122 call mpi_file_write_all(ifile, mpi_io_data%var(i)%sf, data_size*mpi_io_type, mpi_io_p, status, ierr)
1123 end do
1124 if (qbmm .and. .not. polytropic) then
1125 do i = sys_size + 1, sys_size + 2*nb*nnode
1126 var_mok = int(i, mpi_offset_kind)
1127
1128 call mpi_file_write_all(ifile, mpi_io_data%var(i)%sf, data_size*mpi_io_type, mpi_io_p, status, ierr)
1129 end do
1130 end if
1131 else
1132 if (down_sample) then
1133 do i = 1, sys_size ! TODO: check if sys_size is correct
1134 var_mok = int(i, mpi_offset_kind)
1135
1136 call mpi_file_write_all(ifile, q_cons_temp_ds(i)%sf, data_size*mpi_io_type, mpi_io_p, status, ierr)
1137 end do
1138 else
1139 do i = 1, sys_size ! TODO: check if sys_size is correct
1140 var_mok = int(i, mpi_offset_kind)
1141
1142 call mpi_file_write_all(ifile, mpi_io_data%var(i)%sf, data_size*mpi_io_type, mpi_io_p, status, ierr)
1143 end do
1144 end if
1145 end if
1146
1147 call mpi_file_close(ifile, ierr)
1148
1149 if (ib) then
1150 call s_write_parallel_ib_data(t_step)
1151 end if
1152 else
1153 if (ib) then
1154 call s_initialize_mpi_data(q_cons_vf, ib_markers=ib_markers, ib_mpi_data=mpi_io_ib_data, qbmm_pb=pb_ts(1), &
1155 & qbmm_mv=mv_ts(1))
1156 else if (present(beta)) then
1157 call s_initialize_mpi_data(q_cons_vf, beta=beta, qbmm_pb=pb_ts(1), qbmm_mv=mv_ts(1))
1158 else
1159 call s_initialize_mpi_data(q_cons_vf, qbmm_pb=pb_ts(1), qbmm_mv=mv_ts(1))
1160 end if
1161
1162 write (file_loc, '(I0,A)') t_step, '.dat'
1163 file_loc = trim(case_dir) // '/restart_data' // trim(mpiiofs) // trim(file_loc)
1164 inquire (file=trim(file_loc), exist=file_exist)
1165 if (file_exist .and. proc_rank == 0) then
1166 call mpi_file_delete(file_loc, mpi_info_int, ierr)
1167 end if
1168 call mpi_file_open(mpi_comm_world, file_loc, ior(mpi_mode_wronly, mpi_mode_create), mpi_info_int, ifile, ierr)
1169
1170 data_size = (m + 1)*(n + 1)*(p + 1)
1171
1172 m_mok = int(m_glb + 1, mpi_offset_kind)
1173 n_mok = int(n_glb + 1, mpi_offset_kind)
1174 p_mok = int(p_glb + 1, mpi_offset_kind)
1175 wp_mok = int(storage_size(0._stp)/8, mpi_offset_kind)
1176 mok = int(1._wp, mpi_offset_kind)
1177 str_mok = int(name_len, mpi_offset_kind)
1178 nvars_mok = int(alt_sys, mpi_offset_kind)
1179
1180 if (bubbles_euler) then
1181 do i = 1, sys_size
1182 var_mok = int(i, mpi_offset_kind)
1183
1184 disp = m_mok*max(mok, n_mok)*max(mok, p_mok)*wp_mok*(var_mok - 1)
1185
1186 call mpi_file_set_view(ifile, disp, mpi_p, mpi_io_data%view(i), 'native', mpi_info_int, ierr)
1187 call mpi_file_write_all(ifile, mpi_io_data%var(i)%sf, data_size*mpi_io_type, mpi_io_p, status, ierr)
1188 end do
1189 if (qbmm .and. .not. polytropic) then
1190 do i = sys_size + 1, sys_size + 2*nb*nnode
1191 var_mok = int(i, mpi_offset_kind)
1192
1193 disp = m_mok*max(mok, n_mok)*max(mok, p_mok)*wp_mok*(var_mok - 1)
1194
1195 call mpi_file_set_view(ifile, disp, mpi_p, mpi_io_data%view(i), 'native', mpi_info_int, ierr)
1196 call mpi_file_write_all(ifile, mpi_io_data%var(i)%sf, data_size*mpi_io_type, mpi_io_p, status, ierr)
1197 end do
1198 end if
1199 else
1200 do i = 1, sys_size ! TODO: check if sys_size is correct
1201 var_mok = int(i, mpi_offset_kind)
1202
1203 disp = m_mok*max(mok, n_mok)*max(mok, p_mok)*wp_mok*(var_mok - 1)
1204
1205 call mpi_file_set_view(ifile, disp, mpi_p, mpi_io_data%view(i), 'native', mpi_info_int, ierr)
1206 call mpi_file_write_all(ifile, mpi_io_data%var(i)%sf, data_size*mpi_io_type, mpi_io_p, status, ierr)
1207 end do
1208 end if
1209
1210 if (present(beta)) then
1211 var_mok = int(sys_size + 1, mpi_offset_kind)
1212
1213 disp = m_mok*max(mok, n_mok)*max(mok, p_mok)*wp_mok*(var_mok - 1)
1214
1215 call mpi_file_set_view(ifile, disp, mpi_p, mpi_io_data%view(sys_size + 1), 'native', mpi_info_int, ierr)
1216 call mpi_file_write_all(ifile, mpi_io_data%var(sys_size + 1)%sf, data_size*mpi_io_type, mpi_io_p, status, ierr)
1217 end if
1218
1219 call mpi_file_close(ifile, ierr)
1220
1221 if (ib) then
1222 call s_write_parallel_ib_data(t_step)
1223 end if
1224 end if
1225#endif
1226
1227 end subroutine s_write_parallel_data_files
1228
1229 !> Write immersed boundary marker data to a serial (per-processor) unformatted file
1230 subroutine s_write_serial_ib_data(time_step)
1231
1232 integer, intent(in) :: time_step
1233 character(LEN=path_len + 2*name_len) :: file_path
1234 character(LEN=path_len + 2*name_len) :: t_step_dir
1235
1236 write (t_step_dir, '(A,I0,A,I0)') trim(case_dir) // '/p_all'
1237 write (t_step_dir, '(a,i0,a,i0)') trim(case_dir) // '/p_all/p', proc_rank, '/', time_step
1238 write (file_path, '(A,I0,A)') trim(t_step_dir) // '/ib_data.dat'
1239
1240 open (2, file=trim(file_path), form='unformatted', status='new')
1241
1242
1243# 867 "/home/runner/work/MFC/MFC/src/simulation/m_data_output.fpp"
1244#if defined(MFC_OpenACC)
1245# 867 "/home/runner/work/MFC/MFC/src/simulation/m_data_output.fpp"
1246!$acc update host(ib_markers%sf)
1247# 867 "/home/runner/work/MFC/MFC/src/simulation/m_data_output.fpp"
1248#elif defined(MFC_OpenMP)
1249# 867 "/home/runner/work/MFC/MFC/src/simulation/m_data_output.fpp"
1250!$omp target update from(ib_markers%sf)
1251# 867 "/home/runner/work/MFC/MFC/src/simulation/m_data_output.fpp"
1252#endif
1253 write (2) ib_markers%sf(0:m,0:n,0:p); close (2)
1254
1255 end subroutine s_write_serial_ib_data
1256
1257 !> Write immersed boundary marker data in parallel using MPI I/O
1258 subroutine s_write_parallel_ib_data(time_step)
1259
1260 integer, intent(in) :: time_step
1261
1262#ifdef MFC_MPI
1263 character(LEN=path_len + 2*name_len) :: file_loc
1264 integer(kind=MPI_OFFSET_kind) :: disp
1265 integer(kind=MPI_OFFSET_kind) :: m_MOK, n_MOK, p_MOK
1266 integer(kind=MPI_OFFSET_kind) :: WP_MOK, var_MOK, MOK
1267 integer :: ifile, ierr, data_size
1268 integer, dimension(MPI_STATUS_SIZE) :: status
1269
1270
1271# 885 "/home/runner/work/MFC/MFC/src/simulation/m_data_output.fpp"
1272#if defined(MFC_OpenACC)
1273# 885 "/home/runner/work/MFC/MFC/src/simulation/m_data_output.fpp"
1274!$acc update host(ib_markers%sf)
1275# 885 "/home/runner/work/MFC/MFC/src/simulation/m_data_output.fpp"
1276#elif defined(MFC_OpenMP)
1277# 885 "/home/runner/work/MFC/MFC/src/simulation/m_data_output.fpp"
1278!$omp target update from(ib_markers%sf)
1279# 885 "/home/runner/work/MFC/MFC/src/simulation/m_data_output.fpp"
1280#endif
1281
1282 data_size = (m + 1)*(n + 1)*(p + 1)
1283 m_mok = int(m_glb + 1, mpi_offset_kind)
1284 n_mok = int(n_glb + 1, mpi_offset_kind)
1285 p_mok = int(p_glb + 1, mpi_offset_kind)
1286 wp_mok = int(storage_size(0._stp)/8, mpi_offset_kind)
1287 mok = int(1._wp, mpi_offset_kind)
1288
1289 write (file_loc, '(A)') 'ib.dat'
1290 file_loc = trim(case_dir) // '/restart_data' // trim(mpiiofs) // trim(file_loc)
1291
1292 call s_mpi_barrier()
1294
1295 call mpi_file_open(mpi_comm_world, file_loc, ior(mpi_mode_wronly, mpi_mode_create), mpi_info_int, ifile, ierr)
1296
1297 var_mok = int(sys_size + 1, mpi_offset_kind)
1298 disp = m_mok*max(mok, n_mok)*max(mok, p_mok)*wp_mok*(var_mok - 1 + int(time_step/t_step_save))
1299 if (time_step == 0) disp = 0
1300
1301 call mpi_file_set_view(ifile, disp, mpi_integer, mpi_io_ib_data%view, 'native', mpi_info_int, ierr)
1302 call mpi_file_write_all(ifile, mpi_io_ib_data%var%sf, data_size, mpi_integer, status, ierr)
1303 call mpi_file_close(ifile, ierr)
1304#endif
1305
1306 end subroutine s_write_parallel_ib_data
1307
1308 !> Dispatch immersed boundary data output to the serial or parallel writer
1309 subroutine s_write_ib_data_file(time_step)
1310
1311 integer, intent(in) :: time_step
1312
1313
1314# 918 "/home/runner/work/MFC/MFC/src/simulation/m_data_output.fpp"
1315#if defined(MFC_OpenACC)
1316# 918 "/home/runner/work/MFC/MFC/src/simulation/m_data_output.fpp"
1317!$acc update host(patch_ib(1:num_ibs))
1318# 918 "/home/runner/work/MFC/MFC/src/simulation/m_data_output.fpp"
1319#elif defined(MFC_OpenMP)
1320# 918 "/home/runner/work/MFC/MFC/src/simulation/m_data_output.fpp"
1321!$omp target update from(patch_ib(1:num_ibs))
1322# 918 "/home/runner/work/MFC/MFC/src/simulation/m_data_output.fpp"
1323#endif
1324
1325 if (parallel_io) then
1326 call s_write_parallel_ib_data(time_step)
1327 else
1328 call s_write_serial_ib_data(time_step)
1329 end if
1330
1331 end subroutine s_write_ib_data_file
1332
1333 !> Writes the IB state information out to file
1334 subroutine s_write_parallel_ib_state(t_step)
1335
1336 integer, intent(in) :: t_step
1337
1338#ifdef MFC_MPI
1339 character(LEN=path_len + 2*name_len) :: file_loc
1340 integer(kind=MPI_OFFSET_KIND) :: disp
1341 integer(kind=MPI_OFFSET_KIND) :: WP_MOK
1342 integer :: ifile, ierr
1343 integer, dimension(MPI_STATUS_SIZE) :: status
1344 logical :: file_exist, dir_check
1345 integer :: i, ib_idx
1346 integer, parameter :: NFIELDS_PER_IB = 20
1347 real(wp) :: ib_buf(NFIELDS_PER_IB)
1348 integer :: file_unit
1349 character(len=10) :: t_step_string
1350
1351 ! Partition IBs across ranks round-robin style
1352 integer :: ib_start, ib_end, nibs_per_rank, remainder
1353
1354 wp_mok = int(storage_size(0._wp)/8, mpi_offset_kind)
1355
1356 if (file_per_process) then
1357 call s_int_to_str(t_step, t_step_string)
1358
1359 if (proc_rank == 0) then
1360 file_loc = trim(case_dir) // '/restart_data/lustre_' // trim(t_step_string)
1361 call s_create_directory(trim(file_loc))
1362 end if
1363 call s_mpi_barrier()
1365
1366 write (file_loc, '(A,I0,A,i7.7,A)') 'ib_state_', t_step, '_', proc_rank, '.dat'
1367 file_loc = trim(case_dir) // '/restart_data/lustre_' // trim(t_step_string) // '/' // trim(file_loc)
1368
1369 inquire (file=trim(file_loc), exist=file_exist)
1370 if (file_exist) then
1371 open (newunit=file_unit, file=trim(file_loc), form='unformatted', access='stream', status='replace')
1372 else
1373 open (newunit=file_unit, file=trim(file_loc), form='unformatted', access='stream', status='new')
1374 end if
1375
1376 write (file_unit) num_local_ibs
1377 do i = 1, num_local_ibs
1378 ib_idx = local_ib_patch_ids(i)
1379 ib_buf(1) = mytime
1380 ib_buf(2:4) = patch_ib(ib_idx)%force(1:3)
1381 ib_buf(5:7) = patch_ib(ib_idx)%torque(1:3)
1382 ib_buf(8:10) = patch_ib(ib_idx)%vel(1:3)
1383 ib_buf(11:13) = patch_ib(ib_idx)%angular_vel(1:3)
1384 ib_buf(14:16) = patch_ib(ib_idx)%angles(1:3)
1385 ib_buf(17) = patch_ib(ib_idx)%x_centroid
1386 ib_buf(18) = patch_ib(ib_idx)%y_centroid
1387 ib_buf(19) = patch_ib(ib_idx)%z_centroid
1388 ib_buf(20) = patch_ib(ib_idx)%radius
1389
1390 write (file_unit) patch_ib(ib_idx)%gbl_patch_id
1391 write (file_unit) ib_buf
1392 end do
1393
1394 close (file_unit)
1395 else
1396 if (proc_rank == 0) then
1397 call s_create_directory(trim(case_dir) // '/restart_data')
1398 end if
1399 call s_mpi_barrier()
1400
1401 write (file_loc, '(A,I0,A)') '/restart_data/ib_state_', t_step, '.dat'
1402 file_loc = trim(case_dir) // trim(file_loc)
1403
1404 inquire (file=trim(file_loc), exist=file_exist)
1405 if (file_exist .and. proc_rank == 0) then
1406 call mpi_file_delete(file_loc, mpi_info_int, ierr)
1407 end if
1408 call s_mpi_barrier()
1409
1410 call mpi_file_open(mpi_comm_world, file_loc, ior(mpi_mode_wronly, mpi_mode_create), mpi_info_int, ifile, ierr)
1411
1412 do i = 1, num_local_ibs
1413 ib_idx = local_ib_patch_ids(i)
1414 ib_buf(1) = mytime
1415 ib_buf(2:4) = patch_ib(ib_idx)%force(1:3)
1416 ib_buf(5:7) = patch_ib(ib_idx)%torque(1:3)
1417 ib_buf(8:10) = patch_ib(ib_idx)%vel(1:3)
1418 ib_buf(11:13) = patch_ib(ib_idx)%angular_vel(1:3)
1419 ib_buf(14:16) = patch_ib(ib_idx)%angles(1:3)
1420 ib_buf(17) = patch_ib(ib_idx)%x_centroid
1421 ib_buf(18) = patch_ib(ib_idx)%y_centroid
1422 ib_buf(19) = patch_ib(ib_idx)%z_centroid
1423 ib_buf(20) = patch_ib(ib_idx)%radius
1424
1425 ! Global IB index determines position in file
1426 disp = int(patch_ib(ib_idx)%gbl_patch_id - 1, mpi_offset_kind)*int(nfields_per_ib, mpi_offset_kind)*wp_mok
1427
1428 call mpi_file_write_at(ifile, disp, ib_buf, nfields_per_ib, mpi_p, status, ierr)
1429 end do
1430
1431 call mpi_file_close(ifile, ierr)
1432 end if
1433#endif
1434
1435 end subroutine s_write_parallel_ib_state
1436
1437 !> Write IB state data to a per-timestep serial (unformatted) file
1438 subroutine s_write_serial_ib_state(t_step)
1439
1440 integer, intent(in) :: t_step
1441 character(LEN=path_len + 2*name_len) :: file_loc
1442 integer :: i, ios, file_unit
1443 integer, parameter :: NFIELDS_PER_IB = 20
1444 real(wp) :: ib_buf(NFIELDS_PER_IB)
1445
1446 call s_create_directory(trim(case_dir) // '/restart_data')
1447
1448 write (file_loc, '(A,I0,A)') '/restart_data/ib_state_', t_step, '.dat'
1449 file_loc = trim(case_dir) // trim(file_loc)
1450
1451 open (newunit=file_unit, file=trim(file_loc), form='unformatted', access='stream', status='replace', iostat=ios)
1452 if (ios /= 0) call s_mpi_abort('Cannot open IB state output file: ' // trim(file_loc))
1453
1454 do i = 1, num_ibs
1455 ib_buf(1) = mytime
1456 ib_buf(2:4) = patch_ib(i)%force(1:3)
1457 ib_buf(5:7) = patch_ib(i)%torque(1:3)
1458 ib_buf(8:10) = patch_ib(i)%vel(1:3)
1459 ib_buf(11:13) = patch_ib(i)%angular_vel(1:3)
1460 ib_buf(14:16) = patch_ib(i)%angles(1:3)
1461 ib_buf(17) = patch_ib(i)%x_centroid
1462 ib_buf(18) = patch_ib(i)%y_centroid
1463 ib_buf(19) = patch_ib(i)%z_centroid
1464 ib_buf(20) = patch_ib(i)%radius
1465
1466 write (file_unit) ib_buf
1467 end do
1468
1469 close (file_unit)
1470
1471 end subroutine s_write_serial_ib_state
1472
1473 !> @brief Writes IB state records to restart_data/ib_state.dat. Must be called only on rank 0.
1474 impure subroutine s_write_ib_state_file(time_step)
1475
1476 integer, intent(in) :: time_step
1477
1478
1479# 1073 "/home/runner/work/MFC/MFC/src/simulation/m_data_output.fpp"
1480#if defined(MFC_OpenACC)
1481# 1073 "/home/runner/work/MFC/MFC/src/simulation/m_data_output.fpp"
1482!$acc update host(patch_ib(1:num_ibs))
1483# 1073 "/home/runner/work/MFC/MFC/src/simulation/m_data_output.fpp"
1484#elif defined(MFC_OpenMP)
1485# 1073 "/home/runner/work/MFC/MFC/src/simulation/m_data_output.fpp"
1486!$omp target update from(patch_ib(1:num_ibs))
1487# 1073 "/home/runner/work/MFC/MFC/src/simulation/m_data_output.fpp"
1488#endif
1489
1490 if (parallel_io) then
1491 call s_write_parallel_ib_state(time_step)
1492 else
1493 call s_write_serial_ib_state(time_step)
1494 end if
1495
1496 end subroutine s_write_ib_state_file
1497
1498 !> Write center-of-mass data at the current time step
1499 impure subroutine s_write_com_files(t_step, c_mass_in)
1500
1501 integer, intent(in) :: t_step
1502 real(wp), dimension(num_fluids, 5), intent(in) :: c_mass_in
1503 integer :: i !< Generic loop iterator
1504 real(wp) :: nondim_time !< Non-dimensional time
1505
1506 if (t_step_old /= dflt_int) then
1507 nondim_time = real(t_step + t_step_old, wp)*dt
1508 else
1509 nondim_time = real(t_step, wp)*dt
1510 end if
1511
1512 if (proc_rank == 0) then
1513 if (n == 0) then
1514 do i = 1, num_fluids
1515 write (i + 120, '(6X,4F24.12)') nondim_time, c_mass_in(i, 1), c_mass_in(i, 2), c_mass_in(i, 5)
1516 end do
1517 else if (p == 0) then
1518 do i = 1, num_fluids
1519 write (i + 120, '(6X,5F24.12)') nondim_time, c_mass_in(i, 1), c_mass_in(i, 2), c_mass_in(i, 3), c_mass_in(i, 5)
1520 end do
1521 else
1522 do i = 1, num_fluids
1523 write (i + 120, '(6X,6F24.12)') nondim_time, c_mass_in(i, 1), c_mass_in(i, 2), c_mass_in(i, 3), c_mass_in(i, &
1524 & 4), c_mass_in(i, 5)
1525 end do
1526 end if
1527 end if
1528
1529 end subroutine s_write_com_files
1530
1531 !> Write flow probe data at the current time step
1532 impure subroutine s_write_probe_files(t_step, q_cons_vf, accel_mag)
1533
1534 integer, intent(in) :: t_step
1535 type(scalar_field), dimension(sys_size), intent(in) :: q_cons_vf
1536 real(wp), dimension(0:m,0:n,0:p), intent(in) :: accel_mag
1537 real(wp), dimension(-1:m) :: distx
1538 real(wp), dimension(-1:n) :: disty
1539 real(wp), dimension(-1:p) :: distz
1540
1541 ! The cell-averaged partial densities, density, velocity, pressure, volume fractions, specific heat ratio function, liquid
1542 ! stiffness function, and sound speed.
1543 real(wp) :: lit_gamma, nbub
1544 real(wp) :: rho
1545 real(wp), dimension(num_vels) :: vel
1546 real(wp) :: pres
1547 real(wp) :: ptilde
1548 real(wp) :: ptot
1549 real(wp) :: alf
1550 real(wp) :: alfgr
1551 real(wp), dimension(num_fluids) :: alpha
1552 real(wp) :: gamma
1553 real(wp) :: pi_inf
1554 real(wp) :: qv
1555 real(wp) :: c
1556 real(wp) :: m00, m10, m01, m20, m02
1557 real(wp) :: varr, varv
1558 real(wp), dimension(Nb) :: nr, r, nrdot, rdot
1559 real(wp) :: nr3
1560 real(wp) :: accel
1561 real(wp) :: int_pres
1562 real(wp) :: max_pres
1563 real(wp), dimension(2) :: re
1564 real(wp), dimension(6) :: tau_e
1565 real(wp) :: g_local
1566 real(wp) :: dyn_p, t
1567 real(wp) :: damage_state
1568 integer :: i, j, k, l, s, d !< Generic loop iterator
1569 real(wp) :: nondim_time !< Non-dimensional time
1570 real(wp) :: tmp !< Temporary variable to store quantity for mpi_allreduce
1571 real(wp) :: rhoyks(1:num_species)
1572
1573 t = dflt_t_guess
1574
1575 if (time_stepper == 23) then
1576 nondim_time = mytime
1577 else
1578 if (t_step_old /= dflt_int) then
1579 nondim_time = real(t_step + t_step_old, wp)*dt
1580 else
1581 nondim_time = real(t_step, wp)*dt
1582 end if
1583 end if
1584
1585 do i = 1, num_probes
1586 rho = 0._wp
1587 do s = 1, num_vels
1588 vel(s) = 0._wp
1589 end do
1590 pres = 0._wp
1591 gamma = 0._wp
1592 pi_inf = 0._wp
1593 qv = 0._wp
1594 c = 0._wp
1595 accel = 0._wp
1596 nr = 0._wp; r = 0._wp
1597 nrdot = 0._wp; rdot = 0._wp
1598 nbub = 0._wp
1599 m00 = 0._wp
1600 m10 = 0._wp
1601 m01 = 0._wp
1602 m20 = 0._wp
1603 m02 = 0._wp
1604 varr = 0._wp; varv = 0._wp
1605 alf = 0._wp
1606 do s = 1, (num_dims*(num_dims + 1))/2
1607 tau_e(s) = 0._wp
1608 end do
1609 damage_state = 0._wp
1610
1611 if (n == 0) then
1612 if ((probe(i)%x >= x_cb(-1)) .and. (probe(i)%x <= x_cb(m))) then
1613 do s = -1, m
1614 distx(s) = x_cb(s) - probe(i)%x
1615 if (distx(s) < 0._wp) distx(s) = 1000._wp
1616 end do
1617 j = minloc(distx, 1)
1618 if (j == 1) j = 2 ! Pick first point if probe is at edge
1619 k = 0
1620 l = 0
1621
1622 if (chemistry) then
1623 do d = 1, num_species
1624 rhoyks(d) = q_cons_vf(eqn_idx%species%beg + d - 1)%sf(j - 2, k, l)
1625 end do
1626 end if
1627
1628 ! Computing/Sharing necessary state variables
1629 if (hypoelasticity) then
1630 call s_convert_to_mixture_variables(q_cons_vf, j - 2, k, l, rho, gamma, pi_inf, qv, re, g_local, &
1631 & fluid_pp(:)%G)
1632 else
1633 call s_convert_to_mixture_variables(q_cons_vf, j - 2, k, l, rho, gamma, pi_inf, qv)
1634 end if
1635 do s = 1, num_vels
1636 vel(s) = q_cons_vf(eqn_idx%cont%end + s)%sf(j - 2, k, l)/rho
1637 end do
1638 do s = 1, num_fluids
1639 alpha(s) = q_cons_vf(eqn_idx%adv%beg + s - 1)%sf(j - 2, k, l)
1640 end do
1641
1642 dyn_p = 0.5_wp*rho*dot_product(vel, vel)
1643
1644 if (hypoelasticity) then
1645 if (cont_damage) then
1646 damage_state = q_cons_vf(eqn_idx%damage)%sf(j - 2, k, l)
1647 g_local = g_local*max((1._wp - damage_state), 0._wp)
1648 end if
1649
1650 call s_compute_pressure(q_cons_vf(eqn_idx%E)%sf(j - 2, k, l), q_cons_vf(eqn_idx%alf)%sf(j - 2, k, l), &
1651 & dyn_p, pi_inf, gamma, rho, qv, rhoyks(:), pres, t, &
1652 & q_cons_vf(eqn_idx%stress%beg)%sf(j - 2, k, l), &
1653 & q_cons_vf(eqn_idx%mom%beg)%sf(j - 2, k, l), g_local)
1654 else
1655 call s_compute_pressure(q_cons_vf(eqn_idx%E)%sf(j - 2, k, l), q_cons_vf(eqn_idx%alf)%sf(j - 2, k, l), &
1656 & dyn_p, pi_inf, gamma, rho, qv, rhoyks, pres, t)
1657 end if
1658
1659 if (hypoelasticity) then
1660 tau_e(1) = q_cons_vf(eqn_idx%stress%end)%sf(j - 2, k, l)/rho
1661 end if
1662
1663 if (bubbles_euler) then
1664 alf = q_cons_vf(eqn_idx%alf)%sf(j - 2, k, l)
1665 if (num_fluids == 3) then
1666 alfgr = q_cons_vf(eqn_idx%alf - 1)%sf(j - 2, k, l)
1667 end if
1668 do s = 1, nb
1669 nr(s) = q_cons_vf(qbmm_idx%rs(s))%sf(j - 2, k, l)
1670 nrdot(s) = q_cons_vf(qbmm_idx%vs(s))%sf(j - 2, k, l)
1671 end do
1672
1673 if (adv_n) then
1674 nbub = q_cons_vf(eqn_idx%n)%sf(j - 2, k, l)
1675 else
1676 nr3 = 0._wp
1677 do s = 1, nb
1678 nr3 = nr3 + weight(s)*(nr(s)**3._wp)
1679 end do
1680
1681 nbub = sqrt((4._wp*pi/3._wp)*nr3/alf)
1682 end if
1683#ifdef MFC_DEBUG
1684 print *, 'In probe, nbub: ', nbub
1685#endif
1686 if (qbmm) then
1687 m00 = q_cons_vf(qbmm_idx%moms(1, 1))%sf(j - 2, k, l)/nbub
1688 m10 = q_cons_vf(qbmm_idx%moms(1, 2))%sf(j - 2, k, l)/nbub
1689 m01 = q_cons_vf(qbmm_idx%moms(1, 3))%sf(j - 2, k, l)/nbub
1690 m20 = q_cons_vf(qbmm_idx%moms(1, 4))%sf(j - 2, k, l)/nbub
1691 m02 = q_cons_vf(qbmm_idx%moms(1, 6))%sf(j - 2, k, l)/nbub
1692
1693 m10 = m10/m00
1694 m01 = m01/m00
1695 m20 = m20/m00
1696 m02 = m02/m00
1697
1698 varr = m20 - m10**2._wp
1699 varv = m02 - m01**2._wp
1700 end if
1701 r(:) = nr(:)/nbub
1702 rdot(:) = nrdot(:)/nbub
1703
1704 ptilde = ptil(j - 2, k, l)
1705 ptot = pres - ptilde
1706 end if
1707
1708 ! Compute mixture sound Speed
1709 call s_compute_speed_of_sound(pres, rho, gamma, pi_inf, ((gamma + 1._wp)*pres + pi_inf)/rho, alpha, 0._wp, &
1710 & 0._wp, c, qv)
1711 if (hypoelasticity) c = sqrt(c*c + (4._wp/3._wp)*g_local/rho)
1712
1713 accel = accel_mag(j - 2, k, l)
1714 end if
1715 else if (p == 0) then
1716 if (chemistry) then
1717 do d = 1, num_species
1718 rhoyks(d) = q_cons_vf(eqn_idx%species%beg + d - 1)%sf(j - 2, k - 2, l)
1719 end do
1720 end if
1721
1722 if ((probe(i)%x >= x_cb(-1)) .and. (probe(i)%x <= x_cb(m))) then
1723 if ((probe(i)%y >= y_cb(-1)) .and. (probe(i)%y <= y_cb(n))) then
1724 do s = -1, m
1725 distx(s) = x_cb(s) - probe(i)%x
1726 if (distx(s) < 0._wp) distx(s) = 1000._wp
1727 end do
1728 do s = -1, n
1729 disty(s) = y_cb(s) - probe(i)%y
1730 if (disty(s) < 0._wp) disty(s) = 1000._wp
1731 end do
1732 j = minloc(distx, 1)
1733 k = minloc(disty, 1)
1734 if (j == 1) j = 2 ! Pick first point if probe is at edge
1735 if (k == 1) k = 2 ! Pick first point if probe is at edge
1736 l = 0
1737
1738 ! Computing/Sharing necessary state variables
1739 call s_convert_to_mixture_variables(q_cons_vf, j - 2, k - 2, l, rho, gamma, pi_inf, qv, re, g_local, &
1740 & fluid_pp(:)%G)
1741 do s = 1, num_vels
1742 vel(s) = q_cons_vf(eqn_idx%cont%end + s)%sf(j - 2, k - 2, l)/rho
1743 end do
1744 do s = 1, num_fluids
1745 alpha(s) = q_cons_vf(eqn_idx%adv%beg + s - 1)%sf(j - 2, k - 2, l)
1746 end do
1747
1748 dyn_p = 0.5_wp*rho*dot_product(vel, vel)
1749
1750 if (hypoelasticity) then
1751 if (cont_damage) then
1752 damage_state = q_cons_vf(eqn_idx%damage)%sf(j - 2, k - 2, l)
1753 g_local = g_local*max((1._wp - damage_state), 0._wp)
1754 end if
1755
1756 call s_compute_pressure(q_cons_vf(eqn_idx%E)%sf(j - 2, k - 2, l), q_cons_vf(eqn_idx%alf)%sf(j - 2, &
1757 & k - 2, l), dyn_p, pi_inf, gamma, rho, qv, rhoyks, pres, t, &
1758 & q_cons_vf(eqn_idx%stress%beg)%sf(j - 2, k - 2, l), &
1759 & q_cons_vf(eqn_idx%mom%beg)%sf(j - 2, k - 2, l), g_local)
1760 else
1761 call s_compute_pressure(q_cons_vf(eqn_idx%E)%sf(j - 2, k - 2, l), q_cons_vf(eqn_idx%alf)%sf(j - 2, &
1762 & k - 2, l), dyn_p, pi_inf, gamma, rho, qv, rhoyks, pres, t)
1763 end if
1764
1765 if (hypoelasticity) then
1766 do s = 1, 3
1767 tau_e(s) = q_cons_vf(eqn_idx%stress%beg + s - 1)%sf(j - 2, k - 2, l)/rho
1768 end do
1769 end if
1770
1771 if (bubbles_euler) then
1772 alf = q_cons_vf(eqn_idx%alf)%sf(j - 2, k - 2, l)
1773 do s = 1, nb
1774 nr(s) = q_cons_vf(qbmm_idx%rs(s))%sf(j - 2, k - 2, l)
1775 nrdot(s) = q_cons_vf(qbmm_idx%vs(s))%sf(j - 2, k - 2, l)
1776 end do
1777
1778 if (adv_n) then
1779 nbub = q_cons_vf(eqn_idx%n)%sf(j - 2, k - 2, l)
1780 else
1781 nr3 = 0._wp
1782 do s = 1, nb
1783 nr3 = nr3 + weight(s)*(nr(s)**3._wp)
1784 end do
1785
1786 nbub = sqrt((4._wp*pi/3._wp)*nr3/alf)
1787 end if
1788
1789 r(:) = nr(:)/nbub
1790 rdot(:) = nrdot(:)/nbub
1791 end if
1792 ! Compute mixture sound speed
1793 call s_compute_speed_of_sound(pres, rho, gamma, pi_inf, ((gamma + 1._wp)*pres + pi_inf)/rho, alpha, &
1794 & 0._wp, 0._wp, c, qv)
1795 if (hypoelasticity) c = sqrt(c*c + (4._wp/3._wp)*g_local/rho)
1796 end if
1797 end if
1798 else
1799 if ((probe(i)%x >= x_cb(-1)) .and. (probe(i)%x <= x_cb(m))) then
1800 if ((probe(i)%y >= y_cb(-1)) .and. (probe(i)%y <= y_cb(n))) then
1801 if ((probe(i)%z >= z_cb(-1)) .and. (probe(i)%z <= z_cb(p))) then
1802 do s = -1, m
1803 distx(s) = x_cb(s) - probe(i)%x
1804 if (distx(s) < 0._wp) distx(s) = 1000._wp
1805 end do
1806 do s = -1, n
1807 disty(s) = y_cb(s) - probe(i)%y
1808 if (disty(s) < 0._wp) disty(s) = 1000._wp
1809 end do
1810 do s = -1, p
1811 distz(s) = z_cb(s) - probe(i)%z
1812 if (distz(s) < 0._wp) distz(s) = 1000._wp
1813 end do
1814 j = minloc(distx, 1)
1815 k = minloc(disty, 1)
1816 l = minloc(distz, 1)
1817 if (j == 1) j = 2 ! Pick first point if probe is at edge
1818 if (k == 1) k = 2 ! Pick first point if probe is at edge
1819 if (l == 1) l = 2 ! Pick first point if probe is at edge
1820
1821 ! Computing/Sharing necessary state variables
1822 call s_convert_to_mixture_variables(q_cons_vf, j - 2, k - 2, l - 2, rho, gamma, pi_inf, qv, re, &
1823 & g_local, fluid_pp(:)%G)
1824 do s = 1, num_vels
1825 vel(s) = q_cons_vf(eqn_idx%cont%end + s)%sf(j - 2, k - 2, l - 2)/rho
1826 end do
1827 do s = 1, num_fluids
1828 alpha(s) = q_cons_vf(eqn_idx%adv%beg + s - 1)%sf(j - 2, k - 2, l - 2)
1829 end do
1830
1831 dyn_p = 0.5_wp*rho*dot_product(vel, vel)
1832
1833 if (chemistry) then
1834 do d = 1, num_species
1835 rhoyks(d) = q_cons_vf(eqn_idx%species%beg + d - 1)%sf(j - 2, k - 2, l - 2)
1836 end do
1837 end if
1838
1839 if (hypoelasticity) then
1840 if (cont_damage) then
1841 damage_state = q_cons_vf(eqn_idx%damage)%sf(j - 2, k - 2, l - 2)
1842 g_local = g_local*max((1._wp - damage_state), 0._wp)
1843 end if
1844
1845 call s_compute_pressure(q_cons_vf(eqn_idx%E)%sf(j - 2, k - 2, l - 2), &
1846 & q_cons_vf(eqn_idx%alf)%sf(j - 2, k - 2, l - 2), dyn_p, pi_inf, gamma, &
1847 & rho, qv, rhoyks, pres, t, q_cons_vf(eqn_idx%stress%beg)%sf(j - 2, &
1848 & k - 2, l - 2), q_cons_vf(eqn_idx%mom%beg)%sf(j - 2, k - 2, l - 2), &
1849 & g_local)
1850 else
1851 call s_compute_pressure(q_cons_vf(eqn_idx%E)%sf(j - 2, k - 2, l - 2), &
1852 & q_cons_vf(eqn_idx%alf)%sf(j - 2, k - 2, l - 2), dyn_p, pi_inf, gamma, &
1853 & rho, qv, rhoyks, pres, t)
1854 end if
1855
1856 if (hypoelasticity) then
1857 do s = 1, 6
1858 tau_e(s) = q_cons_vf(eqn_idx%stress%beg + s - 1)%sf(j - 2, k - 2, l - 2)/rho
1859 end do
1860 end if
1861
1862 ! Compute mixture sound speed
1863 call s_compute_speed_of_sound(pres, rho, gamma, pi_inf, ((gamma + 1._wp)*pres + pi_inf)/rho, alpha, &
1864 & 0._wp, 0._wp, c, qv)
1865 if (hypoelasticity) c = sqrt(c*c + (4._wp/3._wp)*g_local/rho)
1866
1867 accel = accel_mag(j - 2, k - 2, l - 2)
1868 end if
1869 end if
1870 end if
1871 end if
1872 if (num_procs > 1) then
1873# 1459 "/home/runner/work/MFC/MFC/src/simulation/m_data_output.fpp"
1874 tmp = rho
1875 call s_mpi_allreduce_sum(tmp, rho)
1876# 1459 "/home/runner/work/MFC/MFC/src/simulation/m_data_output.fpp"
1877 tmp = pres
1878 call s_mpi_allreduce_sum(tmp, pres)
1879# 1459 "/home/runner/work/MFC/MFC/src/simulation/m_data_output.fpp"
1880 tmp = gamma
1881 call s_mpi_allreduce_sum(tmp, gamma)
1882# 1459 "/home/runner/work/MFC/MFC/src/simulation/m_data_output.fpp"
1883 tmp = pi_inf
1884 call s_mpi_allreduce_sum(tmp, pi_inf)
1885# 1459 "/home/runner/work/MFC/MFC/src/simulation/m_data_output.fpp"
1886 tmp = qv
1887 call s_mpi_allreduce_sum(tmp, qv)
1888# 1459 "/home/runner/work/MFC/MFC/src/simulation/m_data_output.fpp"
1889 tmp = c
1890 call s_mpi_allreduce_sum(tmp, c)
1891# 1459 "/home/runner/work/MFC/MFC/src/simulation/m_data_output.fpp"
1892 tmp = accel
1893 call s_mpi_allreduce_sum(tmp, accel)
1894# 1462 "/home/runner/work/MFC/MFC/src/simulation/m_data_output.fpp"
1895
1896 do s = 1, num_vels
1897 tmp = vel(s)
1898 call s_mpi_allreduce_sum(tmp, vel(s))
1899 end do
1900
1901 if (bubbles_euler) then
1902# 1470 "/home/runner/work/MFC/MFC/src/simulation/m_data_output.fpp"
1903 tmp = alf
1904 call s_mpi_allreduce_sum(tmp, alf)
1905# 1470 "/home/runner/work/MFC/MFC/src/simulation/m_data_output.fpp"
1906 tmp = alfgr
1907 call s_mpi_allreduce_sum(tmp, alfgr)
1908# 1470 "/home/runner/work/MFC/MFC/src/simulation/m_data_output.fpp"
1909 tmp = nbub
1910 call s_mpi_allreduce_sum(tmp, nbub)
1911# 1470 "/home/runner/work/MFC/MFC/src/simulation/m_data_output.fpp"
1912 tmp = nr(1)
1913 call s_mpi_allreduce_sum(tmp, nr(1))
1914# 1470 "/home/runner/work/MFC/MFC/src/simulation/m_data_output.fpp"
1915 tmp = nrdot(1)
1916 call s_mpi_allreduce_sum(tmp, nrdot(1))
1917# 1470 "/home/runner/work/MFC/MFC/src/simulation/m_data_output.fpp"
1918 tmp = m00
1919 call s_mpi_allreduce_sum(tmp, m00)
1920# 1470 "/home/runner/work/MFC/MFC/src/simulation/m_data_output.fpp"
1921 tmp = r(1)
1922 call s_mpi_allreduce_sum(tmp, r(1))
1923# 1470 "/home/runner/work/MFC/MFC/src/simulation/m_data_output.fpp"
1924 tmp = rdot(1)
1925 call s_mpi_allreduce_sum(tmp, rdot(1))
1926# 1470 "/home/runner/work/MFC/MFC/src/simulation/m_data_output.fpp"
1927 tmp = ptilde
1928 call s_mpi_allreduce_sum(tmp, ptilde)
1929# 1470 "/home/runner/work/MFC/MFC/src/simulation/m_data_output.fpp"
1930 tmp = ptot
1931 call s_mpi_allreduce_sum(tmp, ptot)
1932# 1473 "/home/runner/work/MFC/MFC/src/simulation/m_data_output.fpp"
1933
1934 if (qbmm) then
1935# 1476 "/home/runner/work/MFC/MFC/src/simulation/m_data_output.fpp"
1936 tmp = varr
1937 call s_mpi_allreduce_sum(tmp, varr)
1938# 1476 "/home/runner/work/MFC/MFC/src/simulation/m_data_output.fpp"
1939 tmp = varv
1940 call s_mpi_allreduce_sum(tmp, varv)
1941# 1476 "/home/runner/work/MFC/MFC/src/simulation/m_data_output.fpp"
1942 tmp = m10
1943 call s_mpi_allreduce_sum(tmp, m10)
1944# 1476 "/home/runner/work/MFC/MFC/src/simulation/m_data_output.fpp"
1945 tmp = m01
1946 call s_mpi_allreduce_sum(tmp, m01)
1947# 1476 "/home/runner/work/MFC/MFC/src/simulation/m_data_output.fpp"
1948 tmp = m20
1949 call s_mpi_allreduce_sum(tmp, m20)
1950# 1476 "/home/runner/work/MFC/MFC/src/simulation/m_data_output.fpp"
1951 tmp = m02
1952 call s_mpi_allreduce_sum(tmp, m02)
1953# 1479 "/home/runner/work/MFC/MFC/src/simulation/m_data_output.fpp"
1954 end if
1955 end if
1956
1957 if (hypoelasticity) then
1958 do s = 1, (num_dims*(num_dims + 1))/2
1959 tmp = tau_e(s)
1960 call s_mpi_allreduce_sum(tmp, tau_e(s))
1961 end do
1962 end if
1963
1964 if (cont_damage) then
1965 tmp = damage_state
1966 call s_mpi_allreduce_sum(tmp, damage_state)
1967 end if
1968 end if
1969 if (proc_rank == 0) then
1970 if (n == 0) then
1971 if (bubbles_euler .and. (num_fluids <= 2)) then
1972 if (qbmm) then
1973 write (i + 30, '(6x,f12.6,14f28.16)') nondim_time, rho, vel(1), pres, alf, r(1), rdot(1), nr(1), &
1974 & nrdot(1), varr, varv, m10, m01, m20, m02
1975 else
1976 write (i + 30, '(6x,f12.6,8f24.8)') nondim_time, rho, vel(1), pres, alf, r(1), rdot(1), nr(1), nrdot(1)
1977 ! ptilde, & ptot
1978 end if
1979 else if (bubbles_euler .and. (num_fluids == 3)) then
1980 write (i + 30, &
1981 & '(6x,f12.6,f24.8,f24.8,f24.8,f24.8,f24.8,' // 'f24.8,f24.8,f24.8,f24.8,f24.8, f24.8)') &
1982 & nondim_time, rho, vel(1), pres, alf, alfgr, nr(1), nrdot(1), r(1), rdot(1), ptilde, ptot
1983 else if (bubbles_euler .and. num_fluids == 4) then
1984 write (i + 30, &
1985 & '(6x,f12.6,f24.8,f24.8,f24.8,f24.8,' // 'f24.8,f24.8,f24.8,f24.8,f24.8,f24.8,f24.8,f24.8,f24.8)') &
1986 & nondim_time, q_cons_vf(1)%sf(j - 2, 0, 0), q_cons_vf(2)%sf(j - 2, 0, 0), q_cons_vf(3)%sf(j - 2, &
1987 & 0, 0), q_cons_vf(4)%sf(j - 2, 0, 0), q_cons_vf(5)%sf(j - 2, 0, 0), q_cons_vf(6)%sf(j - 2, 0, 0), &
1988 & q_cons_vf(7)%sf(j - 2, 0, 0), q_cons_vf(8)%sf(j - 2, 0, 0), q_cons_vf(9)%sf(j - 2, 0, 0), &
1989 & q_cons_vf(10)%sf(j - 2, 0, 0), nbub, r(1), rdot(1)
1990 else if (hypoelasticity) then
1991 write (i + 30, '(6X,F12.6,F24.8,F24.8,F24.8,F24.8)') nondim_time, rho, vel(1), pres, tau_e(1)
1992 else
1993 write (i + 30, '(6X,F12.6,F24.8,F24.8,F24.8)') nondim_time, rho, vel(1), pres
1994 end if
1995 else if (p == 0) then
1996 if (bubbles_euler) then
1997# 1523 "/home/runner/work/MFC/MFC/src/simulation/m_data_output.fpp"
1998 write (i + 30, '(6X,10F24.8)') nondim_time, rho, vel(1), vel(2), pres, alf, nr(1), nrdot(1), r(1), &
1999 & rdot(1)
2000# 1526 "/home/runner/work/MFC/MFC/src/simulation/m_data_output.fpp"
2001 else if (hypoelasticity) then
2002# 1528 "/home/runner/work/MFC/MFC/src/simulation/m_data_output.fpp"
2003 write (i + 30, '(6X,F12.6,F24.8,F24.8,F24.8,F24.8,' // 'F24.8,F24.8,F24.8)') nondim_time, rho, &
2004 & vel(1), vel(2), pres, tau_e(1), tau_e(2), tau_e(3)
2005# 1531 "/home/runner/work/MFC/MFC/src/simulation/m_data_output.fpp"
2006 else
2007 write (i + 30, '(6X,F12.6,F24.8,F24.8,F24.8)') nondim_time, rho, vel(1), pres
2008 print *, 'time =', nondim_time, 'rho =', rho, 'pres =', pres
2009 end if
2010 else
2011# 1537 "/home/runner/work/MFC/MFC/src/simulation/m_data_output.fpp"
2012 if (hypoelasticity) then
2013 write (i + 30, '(6X,F12.6,16F24.8)') nondim_time, rho, vel(1), vel(2), vel(3), pres, gamma, pi_inf, &
2014 & qv, c, accel, tau_e(1), tau_e(2), tau_e(3), tau_e(4), tau_e(5), tau_e(6)
2015 else
2016 write (i + 30, &
2017 & '(6X,F12.6,F24.8,F24.8,F24.8,F24.8,' // 'F24.8,F24.8,F24.8,F24.8,F24.8,' // 'F24.8)') &
2018 & nondim_time, rho, vel(1), vel(2), vel(3), pres, gamma, pi_inf, qv, c, accel
2019 end if
2020# 1546 "/home/runner/work/MFC/MFC/src/simulation/m_data_output.fpp"
2021 end if
2022 end if
2023 end do
2024
2025 end subroutine s_write_probe_files
2026
2027 !> Write footer with stability criteria extrema and run-time to the information file, then close it
2029
2030 real(wp) :: run_time !< Run-time of the simulation
2031
2032 write (3, '(A)') ' '
2033 write (3, '(A)') ''
2034
2035 write (3, '(A,F9.6)') 'ICFL Max: ', icfl_max
2036 if (surface_tension) write (3, '(A,F9.6)') 'CCFL Max: ', ccfl_max
2037 if (viscous) write (3, '(A,F9.6)') 'VCFL Max: ', vcfl_max
2038 if (viscous) write (3, '(A,ES16.6)') 'Rc Min: ', rc_min
2039
2040 call cpu_time(run_time)
2041
2042 write (3, '(A)') ''
2043 write (3, '(A,I0,A)') 'Run-time: ', int(anint(run_time)), 's'
2044 write (3, '(A)') ' '
2045 close (3)
2046
2048
2049 !> Closes communication files
2050 impure subroutine s_close_com_files()
2051
2052 integer :: i !< Generic loop iterator
2053
2054 do i = 1, num_fluids
2055 close (i + 120)
2056 end do
2057
2058 end subroutine s_close_com_files
2059
2060 !> Closes probe files
2061 impure subroutine s_close_probe_files
2062
2063 integer :: i !< Generic loop iterator
2064
2065 do i = 1, num_probes
2066 close (i + 30)
2067 end do
2068
2069 end subroutine s_close_probe_files
2070
2071 !> Initialize the data output module
2073
2074 integer :: i, m_ds, n_ds, p_ds
2075
2076 if (run_time_info) then
2077 icfl_max = 0._wp
2078 if (surface_tension) then
2079 ccfl_max = 0._wp
2080 end if
2081 if (viscous) then
2082 vcfl_max = 0._wp
2083 rc_min = 1.e12_wp
2084 end if
2085 end if
2086
2087 if (probe_wrt) then
2088#ifdef MFC_DEBUG
2089# 1613 "/home/runner/work/MFC/MFC/src/simulation/m_data_output.fpp"
2090 block
2091# 1613 "/home/runner/work/MFC/MFC/src/simulation/m_data_output.fpp"
2092 use iso_fortran_env, only: output_unit
2093# 1613 "/home/runner/work/MFC/MFC/src/simulation/m_data_output.fpp"
2094
2095# 1613 "/home/runner/work/MFC/MFC/src/simulation/m_data_output.fpp"
2096 print *, 'm_data_output.fpp:1613: ', '@:ALLOCATE(c_mass(num_fluids,5))'
2097# 1613 "/home/runner/work/MFC/MFC/src/simulation/m_data_output.fpp"
2098
2099# 1613 "/home/runner/work/MFC/MFC/src/simulation/m_data_output.fpp"
2100 call flush (output_unit)
2101# 1613 "/home/runner/work/MFC/MFC/src/simulation/m_data_output.fpp"
2102 end block
2103# 1613 "/home/runner/work/MFC/MFC/src/simulation/m_data_output.fpp"
2104#endif
2105# 1613 "/home/runner/work/MFC/MFC/src/simulation/m_data_output.fpp"
2106 allocate (c_mass(num_fluids,5))
2107# 1613 "/home/runner/work/MFC/MFC/src/simulation/m_data_output.fpp"
2108
2109# 1613 "/home/runner/work/MFC/MFC/src/simulation/m_data_output.fpp"
2110
2111# 1613 "/home/runner/work/MFC/MFC/src/simulation/m_data_output.fpp"
2112#if defined(MFC_OpenACC)
2113# 1613 "/home/runner/work/MFC/MFC/src/simulation/m_data_output.fpp"
2114!$acc enter data create(c_mass)
2115# 1613 "/home/runner/work/MFC/MFC/src/simulation/m_data_output.fpp"
2116#elif defined(MFC_OpenMP)
2117# 1613 "/home/runner/work/MFC/MFC/src/simulation/m_data_output.fpp"
2118!$omp target enter data map(always,alloc:c_mass)
2119# 1613 "/home/runner/work/MFC/MFC/src/simulation/m_data_output.fpp"
2120#endif
2121 end if
2122
2123 if (down_sample) then
2124 m_ds = int((m + 1)/3) - 1
2125 n_ds = int((n + 1)/3) - 1
2126 p_ds = int((p + 1)/3) - 1
2127
2128 allocate (q_cons_temp_ds(1:sys_size))
2129 do i = 1, sys_size
2130 allocate (q_cons_temp_ds(i)%sf(-1:m_ds + 1,-1:n_ds + 1,-1:p_ds + 1))
2131 end do
2132 end if
2133
2134 end subroutine s_initialize_data_output_module
2135
2136 !> Module deallocation and/or disassociation procedures
2138
2139 integer :: i
2140
2141 if (probe_wrt) then
2142#ifdef MFC_DEBUG
2143# 1635 "/home/runner/work/MFC/MFC/src/simulation/m_data_output.fpp"
2144 block
2145# 1635 "/home/runner/work/MFC/MFC/src/simulation/m_data_output.fpp"
2146 use iso_fortran_env, only: output_unit
2147# 1635 "/home/runner/work/MFC/MFC/src/simulation/m_data_output.fpp"
2148
2149# 1635 "/home/runner/work/MFC/MFC/src/simulation/m_data_output.fpp"
2150 print *, 'm_data_output.fpp:1635: ', '@:DEALLOCATE(c_mass)'
2151# 1635 "/home/runner/work/MFC/MFC/src/simulation/m_data_output.fpp"
2152
2153# 1635 "/home/runner/work/MFC/MFC/src/simulation/m_data_output.fpp"
2154 call flush (output_unit)
2155# 1635 "/home/runner/work/MFC/MFC/src/simulation/m_data_output.fpp"
2156 end block
2157# 1635 "/home/runner/work/MFC/MFC/src/simulation/m_data_output.fpp"
2158#endif
2159# 1635 "/home/runner/work/MFC/MFC/src/simulation/m_data_output.fpp"
2160
2161# 1635 "/home/runner/work/MFC/MFC/src/simulation/m_data_output.fpp"
2162#if defined(MFC_OpenACC)
2163# 1635 "/home/runner/work/MFC/MFC/src/simulation/m_data_output.fpp"
2164!$acc exit data delete(c_mass)
2165# 1635 "/home/runner/work/MFC/MFC/src/simulation/m_data_output.fpp"
2166#elif defined(MFC_OpenMP)
2167# 1635 "/home/runner/work/MFC/MFC/src/simulation/m_data_output.fpp"
2168!$omp target exit data map(release:c_mass)
2169# 1635 "/home/runner/work/MFC/MFC/src/simulation/m_data_output.fpp"
2170#endif
2171# 1635 "/home/runner/work/MFC/MFC/src/simulation/m_data_output.fpp"
2172 deallocate (c_mass)
2173 end if
2174
2175 if (down_sample) then
2176 do i = 1, sys_size
2177 deallocate (q_cons_temp_ds(i)%sf)
2178 end do
2179 deallocate (q_cons_temp_ds)
2180 end if
2181
2182 end subroutine s_finalize_data_output_module
2183
2184end module m_data_output
type(scalar_field), dimension(sys_size), intent(inout) q_cons_vf
integer, intent(in) k
integer, intent(in) j
integer, intent(in) l
Noncharacteristic and processor boundary condition application for ghost cells and buffer regions.
Platform-specific file and directory operations: create, delete, inquire, getcwd, and basename.
impure subroutine s_delete_directory(dir_name)
Recursively delete a directory using a platform-specific system command.
impure subroutine my_inquire(fileloc, dircheck)
Inquire on the existence of a directory or file.
impure subroutine s_create_directory(dir_name)
Create a directory and all its parents if it does not exist.
Compile-time constant parameters: default values, tolerances, and physical constants.
integer, parameter model_eqns_5eq
integer, parameter name_len
Maximum name length.
real(wp), parameter dflt_t_guess
Default guess for temperature (when a previous value is not available).
integer, parameter dflt_int
Default integer value.
real(wp), parameter sgm_eps
Segmentation tolerance.
integer, parameter nnode
Number of QBMM nodes.
integer, parameter precision_single
real(wp), parameter pi
Pi.
Writes solution data, run-time stability diagnostics (ICFL, VCFL, CCFL, Rc), and probe/center-of-mass...
real(wp), dimension(:,:), allocatable, public c_mass
impure subroutine, public s_open_probe_files
Open flow probe data files for writing.
real(wp) rc_min
Rc criterion maximum.
real(wp) vcfl_max
VCFL criterion maximum.
impure subroutine, public s_write_probe_files(t_step, q_cons_vf, accel_mag)
Write flow probe data at the current time step.
subroutine s_write_parallel_ib_state(t_step)
Writes the IB state information out to file.
impure subroutine, public s_write_com_files(t_step, c_mass_in)
Write center-of-mass data at the current time step.
impure subroutine, public s_initialize_data_output_module
Initialize the data output module.
impure subroutine, public s_finalize_data_output_module
Module deallocation and/or disassociation procedures.
subroutine s_write_serial_ib_state(t_step)
Write IB state data to a per-timestep serial (unformatted) file.
subroutine s_write_serial_ib_data(time_step)
Write immersed boundary marker data to a serial (per-processor) unformatted file.
impure subroutine, public s_close_run_time_information_file
Write footer with stability criteria extrema and run-time to the information file,...
impure subroutine, public s_write_data_files(q_cons_vf, q_t_sf, q_prim_vf, t_step, bc_type, beta)
Write data files. Dispatch subroutine that replaces procedure pointer.
impure subroutine, public s_write_serial_data_files(q_cons_vf, q_t_sf, q_prim_vf, t_step, bc_type, beta)
Write grid and conservative variable data files in serial format.
impure subroutine, public s_close_probe_files
Closes probe files.
impure subroutine, public s_close_com_files()
Closes communication files.
impure subroutine, public s_write_run_time_information(q_prim_vf, t_step)
Write stability criteria extrema to the run-time information file at the given time step.
impure subroutine, public s_open_run_time_information_file
Open the run-time information file and write the stability criteria table header.
real(wp) icfl_max
ICFL criterion maximum.
impure subroutine, public s_write_parallel_data_files(q_cons_vf, t_step, bc_type, beta, q_t_sf)
Write grid and conservative variable data files in parallel via MPI I/O.
impure subroutine, public s_write_ib_state_file(time_step)
Writes IB state records to restart_data/ib_state.dat. Must be called only on rank 0.
type(scalar_field), dimension(:), allocatable q_cons_temp_ds
subroutine s_write_parallel_ib_data(time_step)
Write immersed boundary marker data in parallel using MPI I/O.
impure subroutine, public s_open_com_files()
Open center-of-mass data files for writing.
subroutine, public s_write_ib_data_file(time_step)
Dispatch immersed boundary data output to the serial or parallel writer.
real(wp) ccfl_max
CCFL criterion maximum.
Rank-staggered file access delays to prevent I/O contention on parallel file systems.
impure subroutine, public s_delay_file_access(process_rank)
Introduce a rank-dependent busy-wait delay to stagger parallel file access and reduce I/O contention.
Shared derived types for field data, patch geometry, bubble dynamics, and MPI I/O structures.
Global parameters for the computational domain, fluid properties, and simulation algorithm configurat...
real(wp) mytime
Current simulation time.
real(wp), dimension(:), allocatable fluid_inv_re
per-fluid Newtonian inverse-Re
type(int_bounds_info), dimension(1:3) idwint
real(wp), dimension(:), allocatable, target z_cb
logical any_non_newtonian
.true. if any fluid is non-Newtonian
type(qbmm_idx_info) qbmm_idx
QBMM moment index mappings (allocatable; GPU-managed separately).
integer proc_rank
Rank of the local processor.
type(mpi_io_ib_var), public mpi_io_ib_data
real(wp), dimension(:), allocatable weight
Simpson quadrature weights.
integer, dimension(num_local_ibs_max) local_ib_patch_ids
lookup table of IBs in the local compute domain
type(pres_field), dimension(:), allocatable pb_ts
integer n_el_bubs_glb
Number of Lagrangian bubbles (local and global).
type(pres_field), dimension(:), allocatable mv_ts
integer num_procs
Number of processors.
real(wp), dimension(:), allocatable, target y_cb
real(wp), dimension(:,:,:), allocatable ptil
Pressure modification.
type(mpi_io_var), public mpi_io_data
real(wp), dimension(:), allocatable hb_mu_max
logical, dimension(:), allocatable is_non_newtonian
per-fluid NN flag
real(wp), dimension(:), allocatable, target x_cb
Basic floating-point utilities: approximate equality, default detection, and coordinate bounds.
logical elemental function, public f_approx_equal(a, b, tol_input)
Check if two floating point numbers of wp are within tolerance.
Utility routines for bubble model setup, coordinate transforms, array sampling, and special functions...
subroutine, public s_downsample_data(q_cons_vf, q_cons_temp, m_ds, n_ds, p_ds, m_glb_ds, n_glb_ds, p_glb_ds)
Downsample conservative variable fields by a factor of 3 in each direction using volume averaging.
elemental subroutine, public s_int_to_str(i, res)
Convert an integer to its trimmed string representation.
Ghost-node immersed boundary method: locates ghost/image points, computes interpolation coefficients,...
type(integer_field), public ib_markers
MPI halo exchange, domain decomposition, and buffer packing/unpacking for the simulation solver.
Simulation helper routines for enthalpy computation, CFL calculation, and stability checks.
subroutine, public s_compute_enthalpy(q_prim_vf, pres, rho, gamma, pi_inf, re, h, alpha, vel, vel_sum, qv, j, k, l)
Computes enthalpy.
subroutine, public s_compute_stability_from_dt(vel, c, rho, re_l, j, k, l, icfl, vcfl, rc, ccfl)
Computes stability criterion for a specified dt.
Conservative-to-primitive variable conversion, mixture property evaluation, and pressure computation.
subroutine, public s_compute_pressure(energy, alf, dyn_p, pi_inf, gamma, rho, qv, rhoyks, pres, t, stress, mom, g, pres_mag)
Compute the pressure from the appropriate equation of state.
subroutine, public s_compute_speed_of_sound(pres, rho, gamma, pi_inf, h, adv, vel_sum, c_c, c, qv)
Compute the speed of sound from thermodynamic state variables, supporting multiple equation-of-state ...
subroutine, public s_convert_conservative_to_primitive_variables(qk_cons_vf, q_t_sf, qk_prim_vf, ibounds)
Convert conserved variables (rho*alpha, rho*u, E, alpha) to primitives (rho, u, p,...
subroutine, public s_convert_to_mixture_variables(q_vf, i, j, k, rho, gamma, pi_inf, qv, re_k, g_k, g)
Dispatch to the s_convert_mixture_to_mixture_variables and s_convert_species_to_mixture_variables sub...
Derived type annexing an integer scalar field (SF).
Derived type annexing a scalar field (SF).