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# 145 "/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# 145 "/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# 76 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
87
88# 81 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
89
90# 86 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
91
92# 91 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
93
94# 96 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
95
96# 101 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
97
98# 106 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
99
100# 111 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
101
102# 116 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
103
104# 121 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
105
106# 151 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
107
108# 192 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
109
110# 206 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
111
112# 231 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
113
114# 242 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
115
116# 244 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
117# 255 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
118
119# 284 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
120
121# 294 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
122
123# 304 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
124
125# 313 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
126
127# 330 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
128
129# 340 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
130
131# 347 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
132
133# 353 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
134
135# 359 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
136
137# 365 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
138
139# 371 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
140
141# 377 "/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# 145 "/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# 57 "/home/runner/work/MFC/MFC/src/common/include/macros.fpp"
300
301! Allocate and create GPU device memory
302# 77 "/home/runner/work/MFC/MFC/src/common/include/macros.fpp"
303
304! Free GPU device memory and deallocate
305# 85 "/home/runner/work/MFC/MFC/src/common/include/macros.fpp"
306
307! Cray-specific GPU pointer setup for vector fields
308# 109 "/home/runner/work/MFC/MFC/src/common/include/macros.fpp"
309
310! Cray-specific GPU pointer setup for scalar fields
311# 125 "/home/runner/work/MFC/MFC/src/common/include/macros.fpp"
312
313! Cray-specific GPU pointer setup for acoustic source spatials
314# 150 "/home/runner/work/MFC/MFC/src/common/include/macros.fpp"
315
316# 156 "/home/runner/work/MFC/MFC/src/common/include/macros.fpp"
317
318# 163 "/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
418 call date_and_time(date=file_date)
419
420 write (3, '(A)') 'Date: ' // file_date(5:6) // '/' // file_date(7:8) // '/' // file_date(3:4)
421
422 write (3, '(A)') ''; write (3, '(A)') ''
423
424 write (3, '(13X,A9,13X,A10,13X,A10,13X,A10)', advance="no") trim('Time-step'), trim('dt'), trim('Time'), trim('ICFL Max')
425
426 if (surface_tension) then
427 write (3, '(13X,A10)', advance="no") trim('CCFL Max')
428 end if
429
430 if (viscous) then
431 write (3, '(13X,A10,13X,A16)', advance="no") trim('VCFL Max'), trim('Rc Min')
432 end if
433
434 if (bubbles_lagrange) then
435 write (3, '(13X,A10)', advance="no") trim('N Bubbles')
436 end if
437
438 write (3, *) ! new line
439
441
442 !> Open center-of-mass data files for writing
443 impure subroutine s_open_com_files()
444
445 character(len=path_len + 3*name_len) :: file_path !< Relative path to the CoM file in the case directory
446 integer :: i !< Generic loop iterator
447
448 do i = 1, num_fluids
449 write (file_path, '(A,I0,A)') '/fluid', i, '_com.dat'
450 file_path = trim(case_dir) // trim(file_path)
451 open (i + 120, file=trim(file_path), form='formatted', position='append', status='unknown')
452 if (n == 0) then
453 write (i + 120, '(A)') ' Non-Dimensional Time ' // ' Total Mass ' // ' x-loc ' // ' Total Volume '
454 else if (p == 0) then
455 write (i + 120, &
456 & '(A)') ' Non-Dimensional Time ' // ' Total Mass ' // ' x-loc ' // ' y-loc ' &
457 & // ' Total Volume '
458 else
459 write (i + 120, &
460 & '(A)') ' Non-Dimensional Time ' // ' Total Mass ' // ' x-loc ' // ' y-loc ' // ' z-loc ' &
461 & // ' Total Volume '
462 end if
463 end do
464
465 end subroutine s_open_com_files
466
467 !> Open flow probe data files for writing
468 impure subroutine s_open_probe_files
469
470 character(LEN=path_len + 3*name_len) :: file_path !< Relative path to the probe data file in the case directory
471 integer :: i !< Generic loop iterator
472 logical :: file_exist
473
474 do i = 1, num_probes
475 write (file_path, '(A,I0,A)') '/D/probe', i, '_prim.dat'
476 file_path = trim(case_dir) // trim(file_path)
477
478 inquire (file=trim(file_path), exist=file_exist)
479
480 if (file_exist) then
481 open (i + 30, file=trim(file_path), form='formatted', status='old', position='append')
482 else
483 open (i + 30, file=trim(file_path), form='formatted', status='unknown')
484 end if
485 end do
486
487 if (integral_wrt) then
488 do i = 1, num_integrals
489 write (file_path, '(A,I0,A)') '/D/integral', i, '_prim.dat'
490 file_path = trim(case_dir) // trim(file_path)
491
492 open (i + 70, file=trim(file_path), form='formatted', position='append', status='unknown')
493 end do
494 end if
495
496 end subroutine s_open_probe_files
497
498 !> Write stability criteria extrema to the run-time information file at the given time step
499 impure subroutine s_write_run_time_information(q_prim_vf, t_step)
500
501 type(scalar_field), dimension(sys_size), intent(in) :: q_prim_vf
502 integer, intent(in) :: t_step
503 real(wp) :: rho !< Cell-avg. density
504
505# 174 "/home/runner/work/MFC/MFC/src/simulation/m_data_output.fpp"
506 real(wp), dimension(num_fluids) :: alpha !< Cell-avg. volume fraction
507 real(wp), dimension(num_vels) :: vel !< Cell-avg. velocity
508# 177 "/home/runner/work/MFC/MFC/src/simulation/m_data_output.fpp"
509 real(wp) :: vel_sum !< Cell-avg. velocity sum
510 real(wp) :: pres !< Cell-avg. pressure
511 real(wp) :: gamma !< Cell-avg. sp. heat ratio
512 real(wp) :: pi_inf !< Cell-avg. liquid stiffness function
513 real(wp) :: qv !< Cell-avg. internal energy reference value
514 real(wp) :: c !< Cell-avg. sound speed
515 real(wp) :: h !< Cell-avg. enthalpy
516 real(wp), dimension(2) :: re !< Cell-avg. Reynolds numbers
517 integer :: j, k, l
518 real(wp) :: icfl_max_loc, icfl_max_glb !< ICFL stability extrema on local and global grids
519 real(wp) :: vcfl_max_loc, vcfl_max_glb !< VCFL stability extrema on local and global grids
520 real(wp) :: ccfl_max_loc, ccfl_max_glb !< CCFL stability extrema on local and global grids
521 real(wp) :: rc_min_loc, rc_min_glb !< Rc stability extrema on local and global grids
522 real(wp) :: icfl, vcfl, ccfl, rc
523 integer :: fl !< Fluid loop iterator
524
525 icfl_max_loc = 0._wp
526 vcfl_max_loc = 0._wp
527 ccfl_max_loc = 0._wp
528 rc_min_loc = huge(1.0_wp)
529 ! Computing Stability Criteria at Current Time-step
530
531# 198 "/home/runner/work/MFC/MFC/src/simulation/m_data_output.fpp"
532
533# 198 "/home/runner/work/MFC/MFC/src/simulation/m_data_output.fpp"
534#if defined(MFC_OpenACC)
535# 198 "/home/runner/work/MFC/MFC/src/simulation/m_data_output.fpp"
536!$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) &
537# 198 "/home/runner/work/MFC/MFC/src/simulation/m_data_output.fpp"
538!$acc& reduction(max:icfl_max_loc, vcfl_max_loc, ccfl_max_loc) reduction(min:Rc_min_loc)
539# 198 "/home/runner/work/MFC/MFC/src/simulation/m_data_output.fpp"
540#elif defined(MFC_OpenMP)
541# 198 "/home/runner/work/MFC/MFC/src/simulation/m_data_output.fpp"
542
543# 198 "/home/runner/work/MFC/MFC/src/simulation/m_data_output.fpp"
544
545# 198 "/home/runner/work/MFC/MFC/src/simulation/m_data_output.fpp"
546
547# 198 "/home/runner/work/MFC/MFC/src/simulation/m_data_output.fpp"
548!$omp target teams loop defaultmap(firstprivate:scalar) bind(teams,parallel) collapse(3) defaultmap(tofrom:aggregate) defaultmap(tofrom:allocatable) defaultmap(tofrom:pointer) &
549# 198 "/home/runner/work/MFC/MFC/src/simulation/m_data_output.fpp"
550!$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)
551# 198 "/home/runner/work/MFC/MFC/src/simulation/m_data_output.fpp"
552#endif
553# 201 "/home/runner/work/MFC/MFC/src/simulation/m_data_output.fpp"
554 do l = 0, p
555 do k = 0, n
556 do j = 0, m
557 call s_compute_enthalpy(q_prim_vf, pres, rho, gamma, pi_inf, re, h, alpha, vel, vel_sum, qv, j, k, l)
558
559 call s_compute_speed_of_sound(pres, rho, gamma, pi_inf, h, alpha, vel_sum, 0._wp, c, qv)
560
561 if (any_non_newtonian) then
562 re(1) = 0._wp
563 do fl = 1, num_fluids
564 if (is_non_newtonian(fl)) then
565 re(1) = re(1) + alpha(fl)*hb_mu_max(fl)
566 else
567 re(1) = re(1) + alpha(fl)*fluid_inv_re(fl)
568 end if
569 end do
570 re(1) = 1._wp/max(re(1), sgm_eps)
571 end if
572
573 call s_compute_stability_from_dt(vel, c, rho, re, j, k, l, icfl, vcfl, rc, ccfl)
574
575 icfl_max_loc = max(icfl_max_loc, icfl)
576 vcfl_max_loc = max(vcfl_max_loc, merge(vcfl, 0.0_wp, viscous))
577 ccfl_max_loc = max(ccfl_max_loc, merge(ccfl, 0.0_wp, surface_tension))
578 rc_min_loc = min(rc_min_loc, merge(rc, huge(1.0_wp), viscous))
579 end do
580 end do
581 end do
582
583# 229 "/home/runner/work/MFC/MFC/src/simulation/m_data_output.fpp"
584#if defined(MFC_OpenACC)
585# 229 "/home/runner/work/MFC/MFC/src/simulation/m_data_output.fpp"
586!$acc end parallel loop
587# 229 "/home/runner/work/MFC/MFC/src/simulation/m_data_output.fpp"
588#elif defined(MFC_OpenMP)
589# 229 "/home/runner/work/MFC/MFC/src/simulation/m_data_output.fpp"
590
591# 229 "/home/runner/work/MFC/MFC/src/simulation/m_data_output.fpp"
592!$omp end target teams loop
593# 229 "/home/runner/work/MFC/MFC/src/simulation/m_data_output.fpp"
594#endif
595 ! end: Computing Stability Criteria at Current Time-step
596
597 if (num_procs > 1) then
598 call s_mpi_reduce_stability_criteria_extrema(icfl_max_loc, vcfl_max_loc, rc_min_loc, n_el_bubs_loc, icfl_max_glb, &
599 & vcfl_max_glb, rc_min_glb, n_el_bubs_glb, ccfl_max_loc, ccfl_max_glb)
600 else
601 icfl_max_glb = icfl_max_loc
602 if (viscous) vcfl_max_glb = vcfl_max_loc
603 if (viscous) rc_min_glb = rc_min_loc
604 if (surface_tension) ccfl_max_glb = ccfl_max_loc
605 if (bubbles_lagrange) n_el_bubs_glb = n_el_bubs_loc
606 end if
607
608 if (icfl_max_glb > icfl_max) icfl_max = icfl_max_glb
609
610 if (surface_tension) then
611 if (ccfl_max_glb > ccfl_max) ccfl_max = ccfl_max_glb
612 end if
613
614 if (viscous) then
615 if (vcfl_max_glb > vcfl_max) vcfl_max = vcfl_max_glb
616 if (rc_min_glb < rc_min) rc_min = rc_min_glb
617 end if
618
619 if (proc_rank == 0) then
620 write (3, '(13X,I9,13X,F10.6,13X,F10.6,13X,F10.6)', advance="no") t_step, dt, mytime, icfl_max_glb
621
622 if (surface_tension) then
623 write (3, '(13X,F10.6)', advance="no") ccfl_max_glb
624 end if
625
626 if (viscous) then
627 write (3, '(13X,F10.6,13X,ES16.6)', advance="no") vcfl_max_glb, rc_min_glb
628 end if
629
630 if (bubbles_lagrange) then
631 write (3, '(13X,I10)', advance="no") n_el_bubs_glb
632 end if
633
634 write (3, *) ! new line
635
636 if (.not. f_approx_equal(icfl_max_glb, icfl_max_glb)) then
637 call s_mpi_abort('ICFL is NaN. Exiting.')
638 else if (icfl_max_glb > 1._wp) then
639 print *, 'icfl', icfl_max_glb
640 call s_mpi_abort('ICFL is greater than 1.0. Exiting.')
641 end if
642
643 if (viscous) then
644 if (.not. f_approx_equal(vcfl_max_glb, vcfl_max_glb)) then
645 call s_mpi_abort('VCFL is NaN. Exiting.')
646 else if (vcfl_max_glb > 1._wp) then
647 print *, 'vcfl', vcfl_max_glb
648 call s_mpi_abort('VCFL is greater than 1.0. Exiting.')
649 end if
650 end if
651
652 if (bubbles_lagrange) then
653 if (n_el_bubs_glb == 0) then
654 call s_mpi_abort('No Lagrangian bubbles remain in the domain. Exiting.')
655 end if
656 end if
657 end if
658
659 call s_mpi_barrier()
660
661 end subroutine s_write_run_time_information
662
663 !> Write grid and conservative variable data files in serial format
664 impure subroutine s_write_serial_data_files(q_cons_vf, q_T_sf, q_prim_vf, t_step, bc_type, beta)
665
666 type(scalar_field), dimension(sys_size), intent(inout) :: q_cons_vf
667 type(scalar_field), intent(inout) :: q_t_sf
668 type(scalar_field), dimension(sys_size), intent(inout) :: q_prim_vf
669 integer, intent(in) :: t_step
670 type(scalar_field), intent(inout), optional :: beta
671 type(integer_field), dimension(1:num_dims,-1:1), intent(in) :: bc_type
672 character(LEN=path_len + 2*name_len) :: t_step_dir !< Relative path to the current time-step directory
673 character(LEN=path_len + 3*name_len) :: file_path !< Relative path to the grid and conservative variables data files
674 logical :: file_exist !< Logical used to check existence of current time-step directory
675 character(LEN=15) :: fmt
676 integer :: i, j, k, l, r
677 real(wp) :: gamma, lit_gamma, pi_inf, qv !< Temporary EOS params
678
679 write (t_step_dir, '(A,I0,A,I0)') trim(case_dir) // '/p_all'
680 write (t_step_dir, '(a,i0,a,i0)') trim(case_dir) // '/p_all/p', proc_rank, '/', t_step
681
682 file_path = trim(t_step_dir) // '/.'
683 call my_inquire(file_path, file_exist)
684 if (file_exist) call s_delete_directory(trim(t_step_dir))
685 call s_create_directory(trim(t_step_dir))
686
687 file_path = trim(t_step_dir) // '/x_cb.dat'
688
689 open (2, file=trim(file_path), form='unformatted', status='new')
690 write (2) x_cb(-1:m); close (2)
691
692 if (n > 0) then
693 file_path = trim(t_step_dir) // '/y_cb.dat'
694
695 open (2, file=trim(file_path), form='unformatted', status='new')
696 write (2) y_cb(-1:n); close (2)
697
698 if (p > 0) then
699 file_path = trim(t_step_dir) // '/z_cb.dat'
700
701 open (2, file=trim(file_path), form='unformatted', status='new')
702 write (2) z_cb(-1:p); close (2)
703 end if
704 end if
705
706 do i = 1, sys_size
707 write (file_path, '(A,I0,A)') trim(t_step_dir) // '/q_cons_vf', i, '.dat'
708
709 open (2, file=trim(file_path), form='unformatted', status='new')
710
711 write (2) q_cons_vf(i)%sf(0:m,0:n,0:p); close (2)
712 end do
713
714 ! Lagrangian beta (void fraction) written as q_cons_vf(sys_size+1) to match the parallel I/O path and allow post_process to
715 ! read it.
716 if (bubbles_lagrange) then
717 write (file_path, '(A,I0,A)') trim(t_step_dir) // '/q_cons_vf', sys_size + 1, '.dat'
718
719 open (2, file=trim(file_path), form='unformatted', status='new')
720
721 write (2) beta%sf(0:m,0:n,0:p); close (2)
722 end if
723
724 if (qbmm .and. .not. polytropic) then
725 do i = 1, nb
726 do r = 1, nnode
727 write (file_path, '(A,I0,A)') trim(t_step_dir) // '/pb', sys_size + (i - 1)*nnode + r, '.dat'
728
729 open (2, file=trim(file_path), form='unformatted', status='new')
730
731 write (2) pb_ts(1)%sf(0:m,0:n,0:p,r, i); close (2)
732 end do
733 end do
734
735 do i = 1, nb
736 do r = 1, nnode
737 write (file_path, '(A,I0,A)') trim(t_step_dir) // '/mv', sys_size + (i - 1)*nnode + r, '.dat'
738
739 open (2, file=trim(file_path), form='unformatted', status='new')
740
741 write (2) mv_ts(1)%sf(0:m,0:n,0:p,r, i); close (2)
742 end do
743 end do
744 end if
745
746 ! Writing the IB markers
747 if (ib) then
748 call s_write_serial_ib_data(t_step)
749 end if
750
751 gamma = gammas(1)
752 lit_gamma = gs_min(1)
753 pi_inf = pi_infs(1)
754 qv = qvs(1)
755
756 if (precision == precision_single) then
757 fmt = "(2F30.3)"
758 else
759 fmt = "(2F40.14)"
760 end if
761
762 write (t_step_dir, '(A,I0,A,I0)') trim(case_dir) // '/D'
763 file_path = trim(t_step_dir) // '/.'
764
765 inquire (file=trim(file_path), exist=file_exist)
766
767 if (.not. file_exist) call s_create_directory(trim(t_step_dir))
768
769 if ((prim_vars_wrt .or. (n == 0 .and. p == 0)) .and. (.not. igr)) then
771 do i = 1, sys_size
772
773# 407 "/home/runner/work/MFC/MFC/src/simulation/m_data_output.fpp"
774#if defined(MFC_OpenACC)
775# 407 "/home/runner/work/MFC/MFC/src/simulation/m_data_output.fpp"
776!$acc update host(q_prim_vf(i)%sf(:, :, :))
777# 407 "/home/runner/work/MFC/MFC/src/simulation/m_data_output.fpp"
778#elif defined(MFC_OpenMP)
779# 407 "/home/runner/work/MFC/MFC/src/simulation/m_data_output.fpp"
780!$omp target update from(q_prim_vf(i)%sf(:, :, :))
781# 407 "/home/runner/work/MFC/MFC/src/simulation/m_data_output.fpp"
782#endif
783 end do
784 ! q_prim_vf(eqn_idx%bub%beg) stores the value of nb needed in riemann solvers, so replace with true primitive value
785 ! (=1._wp)
786 if (qbmm) then
787 q_prim_vf(eqn_idx%bub%beg)%sf = 1._wp
788 end if
789 end if
790
791 if (n == 0 .and. p == 0) then
792 if (model_eqns == model_eqns_5eq .and. (.not. igr)) then
793 do i = 1, sys_size
794 write (file_path, '(A,I0,A,I2.2,A,I6.6,A)') trim(t_step_dir) // '/prim.', i, '.', proc_rank, '.', t_step, '.dat'
795
796 open (2, file=trim(file_path))
797 do j = 0, m
798 ! todo: revisit change here
799 if (((i >= eqn_idx%adv%beg) .and. (i <= eqn_idx%adv%end))) then
800 write (2, fmt) x_cb(j), q_cons_vf(i)%sf(j, 0, 0)
801 else
802 write (2, fmt) x_cb(j), q_prim_vf(i)%sf(j, 0, 0)
803 end if
804 end do
805 close (2)
806 end do
807 end if
808
809 do i = 1, sys_size
810 write (file_path, '(A,I0,A,I2.2,A,I6.6,A)') trim(t_step_dir) // '/cons.', i, '.', proc_rank, '.', t_step, '.dat'
811
812 open (2, file=trim(file_path))
813 do j = 0, m
814 write (2, fmt) x_cb(j), q_cons_vf(i)%sf(j, 0, 0)
815 end do
816 close (2)
817 end do
818
819 if (qbmm .and. .not. polytropic) then
820 do i = 1, nb
821 do r = 1, nnode
822 write (file_path, '(A,I0,A,I0,A,I2.2,A,I6.6,A)') trim(t_step_dir) // '/pres.', i, '.', r, '.', proc_rank, &
823 & '.', t_step, '.dat'
824
825 open (2, file=trim(file_path))
826 do j = 0, m
827 write (2, fmt) x_cb(j), pb_ts(1)%sf(j, 0, 0, r, i)
828 end do
829 close (2)
830 end do
831 end do
832 do i = 1, nb
833 do r = 1, nnode
834 write (file_path, '(A,I0,A,I0,A,I2.2,A,I6.6,A)') trim(t_step_dir) // '/mv.', i, '.', r, '.', proc_rank, &
835 & '.', t_step, '.dat'
836
837 open (2, file=trim(file_path))
838 do j = 0, m
839 write (2, fmt) x_cb(j), mv_ts(1)%sf(j, 0, 0, r, i)
840 end do
841 close (2)
842 end do
843 end do
844 end if
845 end if
846
847 if (precision == precision_single) then
848 fmt = "(3F30.7)"
849 else
850 fmt = "(3F40.14)"
851 end if
852
853 if ((n > 0) .and. (p == 0)) then
854 do i = 1, sys_size
855 write (file_path, '(A,I0,A,I2.2,A,I6.6,A)') trim(t_step_dir) // '/cons.', i, '.', proc_rank, '.', t_step, '.dat'
856 open (2, file=trim(file_path))
857 do j = 0, m
858 do k = 0, n
859 write (2, fmt) x_cb(j), y_cb(k), q_cons_vf(i)%sf(j, k, 0)
860 end do
861 write (2, *)
862 end do
863 close (2)
864 end do
865
866 if (present(beta)) then
867 write (file_path, '(A,I0,A,I2.2,A,I6.6,A)') trim(t_step_dir) // '/beta.', i, '.', proc_rank, '.', t_step, '.dat'
868 open (2, file=trim(file_path))
869 do j = 0, m
870 do k = 0, n
871 write (2, fmt) x_cb(j), y_cb(k), beta%sf(j, k, 0)
872 end do
873 write (2, *)
874 end do
875 close (2)
876 end if
877
878 if (qbmm .and. .not. polytropic) then
879 do i = 1, nb
880 do r = 1, nnode
881 write (file_path, '(A,I0,A,I0,A,I2.2,A,I6.6,A)') trim(t_step_dir) // '/pres.', i, '.', r, '.', proc_rank, &
882 & '.', t_step, '.dat'
883
884 open (2, file=trim(file_path))
885 do j = 0, m
886 do k = 0, n
887 write (2, fmt) x_cb(j), y_cb(k), pb_ts(1)%sf(j, k, 0, r, i)
888 end do
889 end do
890 close (2)
891 end do
892 end do
893 do i = 1, nb
894 do r = 1, nnode
895 write (file_path, '(A,I0,A,I0,A,I2.2,A,I6.6,A)') trim(t_step_dir) // '/mv.', i, '.', r, '.', proc_rank, &
896 & '.', t_step, '.dat'
897
898 open (2, file=trim(file_path))
899 do j = 0, m
900 do k = 0, n
901 write (2, fmt) x_cb(j), y_cb(k), mv_ts(1)%sf(j, k, 0, r, i)
902 end do
903 end do
904 close (2)
905 end do
906 end do
907 end if
908
909 if (prim_vars_wrt .and. (.not. igr)) then
910 do i = 1, sys_size
911 write (file_path, '(A,I0,A,I2.2,A,I6.6,A)') trim(t_step_dir) // '/prim.', i, '.', proc_rank, '.', t_step, '.dat'
912
913 open (2, file=trim(file_path))
914
915 do j = 0, m
916 do k = 0, n
917 if (((i >= eqn_idx%cont%beg) .and. (i <= eqn_idx%cont%end)) .or. ((i >= eqn_idx%adv%beg) &
918 & .and. (i <= eqn_idx%adv%end))) then
919 write (2, fmt) x_cb(j), y_cb(k), q_cons_vf(i)%sf(j, k, 0)
920 else
921 write (2, fmt) x_cb(j), y_cb(k), q_prim_vf(i)%sf(j, k, 0)
922 end if
923 end do
924 write (2, *)
925 end do
926 close (2)
927 end do
928 end if
929 end if
930
931 if (precision == precision_single) then
932 fmt = "(4F30.7)"
933 else
934 fmt = "(4F40.14)"
935 end if
936
937 if (p > 0) then
938 do i = 1, sys_size
939 write (file_path, '(A,I0,A,I2.2,A,I6.6,A)') trim(t_step_dir) // '/cons.', i, '.', proc_rank, '.', t_step, '.dat'
940 open (2, file=trim(file_path))
941 do j = 0, m
942 do k = 0, n
943 do l = 0, p
944 write (2, fmt) x_cb(j), y_cb(k), z_cb(l), q_cons_vf(i)%sf(j, k, l)
945 end do
946 write (2, *)
947 end do
948 write (2, *)
949 end do
950 close (2)
951 end do
952
953 if (present(beta)) then
954 write (file_path, '(A,I0,A,I2.2,A,I6.6,A)') trim(t_step_dir) // '/beta.', i, '.', proc_rank, '.', t_step, '.dat'
955 open (2, file=trim(file_path))
956 do j = 0, m
957 do k = 0, n
958 do l = 0, p
959 write (2, fmt) x_cb(j), y_cb(k), z_cb(l), beta%sf(j, k, l)
960 end do
961 write (2, *)
962 end do
963 write (2, *)
964 end do
965 close (2)
966 end if
967
968 if (qbmm .and. .not. polytropic) then
969 do i = 1, nb
970 do r = 1, nnode
971 write (file_path, '(A,I0,A,I0,A,I2.2,A,I6.6,A)') trim(t_step_dir) // '/pres.', i, '.', r, '.', proc_rank, &
972 & '.', t_step, '.dat'
973
974 open (2, file=trim(file_path))
975 do j = 0, m
976 do k = 0, n
977 do l = 0, p
978 write (2, fmt) x_cb(j), y_cb(k), z_cb(l), pb_ts(1)%sf(j, k, l, r, i)
979 end do
980 end do
981 end do
982 close (2)
983 end do
984 end do
985 do i = 1, nb
986 do r = 1, nnode
987 write (file_path, '(A,I0,A,I0,A,I2.2,A,I6.6,A)') trim(t_step_dir) // '/mv.', i, '.', r, '.', proc_rank, &
988 & '.', t_step, '.dat'
989
990 open (2, file=trim(file_path))
991 do j = 0, m
992 do k = 0, n
993 do l = 0, p
994 write (2, fmt) x_cb(j), y_cb(k), z_cb(l), mv_ts(1)%sf(j, k, l, r, i)
995 end do
996 end do
997 end do
998 close (2)
999 end do
1000 end do
1001 end if
1002
1003 if (prim_vars_wrt .and. (.not. igr)) then
1004 do i = 1, sys_size
1005 write (file_path, '(A,I0,A,I2.2,A,I6.6,A)') trim(t_step_dir) // '/prim.', i, '.', proc_rank, '.', t_step, '.dat'
1006
1007 open (2, file=trim(file_path))
1008
1009 do j = 0, m
1010 do k = 0, n
1011 do l = 0, p
1012 if (((i >= eqn_idx%cont%beg) .and. (i <= eqn_idx%cont%end)) .or. ((i >= eqn_idx%adv%beg) &
1013 & .and. (i <= eqn_idx%adv%end)) .or. ((i >= eqn_idx%species%beg) &
1014 & .and. (i <= eqn_idx%species%end))) then
1015 write (2, fmt) x_cb(j), y_cb(k), z_cb(l), q_cons_vf(i)%sf(j, k, l)
1016 else
1017 write (2, fmt) x_cb(j), y_cb(k), z_cb(l), q_prim_vf(i)%sf(j, k, l)
1018 end if
1019 end do
1020 write (2, *)
1021 end do
1022 write (2, *)
1023 end do
1024 close (2)
1025 end do
1026 end if
1027 end if
1028
1029 end subroutine s_write_serial_data_files
1030
1031 !> Write grid and conservative variable data files in parallel via MPI I/O
1032 impure subroutine s_write_parallel_data_files(q_cons_vf, t_step, bc_type, beta, q_T_sf)
1033
1034 type(scalar_field), dimension(sys_size), intent(inout) :: q_cons_vf
1035 integer, intent(in) :: t_step
1036 type(scalar_field), intent(inout), optional :: beta
1037 type(integer_field), dimension(1:num_dims,-1:1), intent(in) :: bc_type
1038 type(scalar_field), intent(inout), optional :: q_t_sf
1039
1040#ifdef MFC_MPI
1041 integer :: ifile, ierr, data_size
1042 integer, dimension(MPI_STATUS_SIZE) :: status
1043 integer(kind=MPI_OFFSET_kind) :: disp
1044 integer(kind=MPI_OFFSET_kind) :: m_mok, n_mok, p_mok
1045 integer(kind=MPI_OFFSET_kind) :: wp_mok, var_mok, str_mok
1046 integer(kind=MPI_OFFSET_kind) :: nvars_mok
1047 integer(kind=MPI_OFFSET_kind) :: mok
1048 character(LEN=path_len + 2*name_len) :: file_loc
1049 logical :: file_exist, dir_check
1050 character(len=10) :: t_step_string
1051 integer :: i !< Generic loop iterator
1052 integer :: alt_sys !< Altered system size for the lagrangian subgrid bubble model
1053 ! Down sampling variables
1054 integer :: m_ds, n_ds, p_ds
1055 integer :: m_glb_ds, n_glb_ds, p_glb_ds
1056 integer :: m_glb_save, n_glb_save, p_glb_save !< Global save size
1057
1058 if (down_sample) then
1059 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)
1060 end if
1061
1062 if (present(beta)) then
1063 alt_sys = sys_size + 1
1064 else
1065 alt_sys = sys_size
1066 end if
1067
1068 if (file_per_process) then
1069 call s_int_to_str(t_step, t_step_string)
1070
1071 if (down_sample) then
1072 call s_initialize_mpi_data_ds(q_cons_temp_ds)
1073 else
1074 if (ib) then
1075 call s_initialize_mpi_data(q_cons_vf, ib_markers)
1076 else
1077 call s_initialize_mpi_data(q_cons_vf)
1078 end if
1079 end if
1080
1081 if (proc_rank == 0) then
1082 file_loc = trim(case_dir) // '/restart_data/lustre_' // trim(t_step_string)
1083 call my_inquire(file_loc, dir_check)
1084 if (dir_check .neqv. .true.) then
1085 call s_create_directory(trim(file_loc))
1086 end if
1087 call s_create_directory(trim(file_loc))
1088 end if
1089 call s_mpi_barrier()
1091
1092 call s_initialize_mpi_data(q_cons_vf)
1093
1094 write (file_loc, '(I0,A,i7.7,A)') t_step, '_', proc_rank, '.dat'
1095 file_loc = trim(case_dir) // '/restart_data/lustre_' // trim(t_step_string) // trim(mpiiofs) // trim(file_loc)
1096 inquire (file=trim(file_loc), exist=file_exist)
1097 if (file_exist .and. proc_rank == 0) then
1098 call mpi_file_delete(file_loc, mpi_info_int, ierr)
1099 end if
1100 call mpi_file_open(mpi_comm_self, file_loc, ior(mpi_mode_wronly, mpi_mode_create), mpi_info_int, ifile, ierr)
1101
1102 if (down_sample) then
1103 data_size = (m_ds + 3)*(n_ds + 3)*(p_ds + 3)
1104 m_glb_save = m_glb_ds + 1
1105 n_glb_save = n_glb_ds + 1
1106 p_glb_save = p_glb_ds + 1
1107 else
1108 data_size = (m + 1)*(n + 1)*(p + 1)
1109 m_glb_save = m_glb + 1
1110 n_glb_save = n_glb + 1
1111 p_glb_save = p_glb + 1
1112 end if
1113
1114 m_mok = int(m_glb_save + 1, mpi_offset_kind)
1115 n_mok = int(n_glb_save + 1, mpi_offset_kind)
1116 p_mok = int(p_glb_save + 1, mpi_offset_kind)
1117 wp_mok = int(storage_size(0._stp)/8, mpi_offset_kind)
1118 mok = int(1._wp, mpi_offset_kind)
1119 str_mok = int(name_len, mpi_offset_kind)
1120 nvars_mok = int(sys_size, mpi_offset_kind)
1121
1122 if (bubbles_euler) then
1123 do i = 1, sys_size
1124 var_mok = int(i, mpi_offset_kind)
1125
1126 call mpi_file_write_all(ifile, mpi_io_data%var(i)%sf, data_size*mpi_io_type, mpi_io_p, status, ierr)
1127 end do
1128 if (qbmm .and. .not. polytropic) then
1129 do i = sys_size + 1, sys_size + 2*nb*nnode
1130 var_mok = int(i, mpi_offset_kind)
1131
1132 call mpi_file_write_all(ifile, mpi_io_data%var(i)%sf, data_size*mpi_io_type, mpi_io_p, status, ierr)
1133 end do
1134 end if
1135 else
1136 if (down_sample) then
1137 do i = 1, sys_size ! TODO: check if sys_size is correct
1138 var_mok = int(i, mpi_offset_kind)
1139
1140 call mpi_file_write_all(ifile, q_cons_temp_ds(i)%sf, data_size*mpi_io_type, mpi_io_p, status, ierr)
1141 end do
1142 else
1143 do i = 1, sys_size ! TODO: check if sys_size is correct
1144 var_mok = int(i, mpi_offset_kind)
1145
1146 call mpi_file_write_all(ifile, mpi_io_data%var(i)%sf, data_size*mpi_io_type, mpi_io_p, status, ierr)
1147 end do
1148 end if
1149 end if
1150
1151 call mpi_file_close(ifile, ierr)
1152
1153 if (ib) then
1154 call s_write_parallel_ib_data(t_step)
1155 end if
1156 else
1157 if (ib) then
1158 call s_initialize_mpi_data(q_cons_vf, ib_markers)
1159 else if (present(beta)) then
1160 call s_initialize_mpi_data(q_cons_vf, beta=beta)
1161 else
1162 call s_initialize_mpi_data(q_cons_vf)
1163 end if
1164
1165 write (file_loc, '(I0,A)') t_step, '.dat'
1166 file_loc = trim(case_dir) // '/restart_data' // trim(mpiiofs) // trim(file_loc)
1167 inquire (file=trim(file_loc), exist=file_exist)
1168 if (file_exist .and. proc_rank == 0) then
1169 call mpi_file_delete(file_loc, mpi_info_int, ierr)
1170 end if
1171 call mpi_file_open(mpi_comm_world, file_loc, ior(mpi_mode_wronly, mpi_mode_create), mpi_info_int, ifile, ierr)
1172
1173 data_size = (m + 1)*(n + 1)*(p + 1)
1174
1175 m_mok = int(m_glb + 1, mpi_offset_kind)
1176 n_mok = int(n_glb + 1, mpi_offset_kind)
1177 p_mok = int(p_glb + 1, mpi_offset_kind)
1178 wp_mok = int(storage_size(0._stp)/8, mpi_offset_kind)
1179 mok = int(1._wp, mpi_offset_kind)
1180 str_mok = int(name_len, mpi_offset_kind)
1181 nvars_mok = int(alt_sys, mpi_offset_kind)
1182
1183 if (bubbles_euler) then
1184 do i = 1, sys_size
1185 var_mok = int(i, mpi_offset_kind)
1186
1187 disp = m_mok*max(mok, n_mok)*max(mok, p_mok)*wp_mok*(var_mok - 1)
1188
1189 call mpi_file_set_view(ifile, disp, mpi_p, mpi_io_data%view(i), 'native', mpi_info_int, ierr)
1190 call mpi_file_write_all(ifile, mpi_io_data%var(i)%sf, data_size*mpi_io_type, mpi_io_p, status, ierr)
1191 end do
1192 if (qbmm .and. .not. polytropic) then
1193 do i = sys_size + 1, sys_size + 2*nb*nnode
1194 var_mok = int(i, mpi_offset_kind)
1195
1196 disp = m_mok*max(mok, n_mok)*max(mok, p_mok)*wp_mok*(var_mok - 1)
1197
1198 call mpi_file_set_view(ifile, disp, mpi_p, mpi_io_data%view(i), 'native', mpi_info_int, ierr)
1199 call mpi_file_write_all(ifile, mpi_io_data%var(i)%sf, data_size*mpi_io_type, mpi_io_p, status, ierr)
1200 end do
1201 end if
1202 else
1203 do i = 1, sys_size ! TODO: check if sys_size is correct
1204 var_mok = int(i, mpi_offset_kind)
1205
1206 disp = m_mok*max(mok, n_mok)*max(mok, p_mok)*wp_mok*(var_mok - 1)
1207
1208 call mpi_file_set_view(ifile, disp, mpi_p, mpi_io_data%view(i), 'native', mpi_info_int, ierr)
1209 call mpi_file_write_all(ifile, mpi_io_data%var(i)%sf, data_size*mpi_io_type, mpi_io_p, status, ierr)
1210 end do
1211 end if
1212
1213 if (present(beta)) then
1214 var_mok = int(sys_size + 1, mpi_offset_kind)
1215
1216 disp = m_mok*max(mok, n_mok)*max(mok, p_mok)*wp_mok*(var_mok - 1)
1217
1218 call mpi_file_set_view(ifile, disp, mpi_p, mpi_io_data%view(sys_size + 1), 'native', mpi_info_int, ierr)
1219 call mpi_file_write_all(ifile, mpi_io_data%var(sys_size + 1)%sf, data_size*mpi_io_type, mpi_io_p, status, ierr)
1220 end if
1221
1222 call mpi_file_close(ifile, ierr)
1223
1224 if (ib) then
1225 call s_write_parallel_ib_data(t_step)
1226 end if
1227 end if
1228#endif
1229
1230 end subroutine s_write_parallel_data_files
1231
1232 !> Write immersed boundary marker data to a serial (per-processor) unformatted file
1233 subroutine s_write_serial_ib_data(time_step)
1234
1235 integer, intent(in) :: time_step
1236 character(LEN=path_len + 2*name_len) :: file_path
1237 character(LEN=path_len + 2*name_len) :: t_step_dir
1238
1239 write (t_step_dir, '(A,I0,A,I0)') trim(case_dir) // '/p_all'
1240 write (t_step_dir, '(a,i0,a,i0)') trim(case_dir) // '/p_all/p', proc_rank, '/', time_step
1241 write (file_path, '(A,I0,A)') trim(t_step_dir) // '/ib_data.dat'
1242
1243 open (2, file=trim(file_path), form='unformatted', status='new')
1244
1245
1246# 870 "/home/runner/work/MFC/MFC/src/simulation/m_data_output.fpp"
1247#if defined(MFC_OpenACC)
1248# 870 "/home/runner/work/MFC/MFC/src/simulation/m_data_output.fpp"
1249!$acc update host(ib_markers%sf)
1250# 870 "/home/runner/work/MFC/MFC/src/simulation/m_data_output.fpp"
1251#elif defined(MFC_OpenMP)
1252# 870 "/home/runner/work/MFC/MFC/src/simulation/m_data_output.fpp"
1253!$omp target update from(ib_markers%sf)
1254# 870 "/home/runner/work/MFC/MFC/src/simulation/m_data_output.fpp"
1255#endif
1256 write (2) ib_markers%sf(0:m,0:n,0:p); close (2)
1257
1258 end subroutine s_write_serial_ib_data
1259
1260 !> Write immersed boundary marker data in parallel using MPI I/O
1261 subroutine s_write_parallel_ib_data(time_step)
1262
1263 integer, intent(in) :: time_step
1264
1265#ifdef MFC_MPI
1266 character(LEN=path_len + 2*name_len) :: file_loc
1267 integer(kind=MPI_OFFSET_kind) :: disp
1268 integer(kind=MPI_OFFSET_kind) :: m_MOK, n_MOK, p_MOK
1269 integer(kind=MPI_OFFSET_kind) :: WP_MOK, var_MOK, MOK
1270 integer :: ifile, ierr, data_size
1271 integer, dimension(MPI_STATUS_SIZE) :: status
1272
1273
1274# 888 "/home/runner/work/MFC/MFC/src/simulation/m_data_output.fpp"
1275#if defined(MFC_OpenACC)
1276# 888 "/home/runner/work/MFC/MFC/src/simulation/m_data_output.fpp"
1277!$acc update host(ib_markers%sf)
1278# 888 "/home/runner/work/MFC/MFC/src/simulation/m_data_output.fpp"
1279#elif defined(MFC_OpenMP)
1280# 888 "/home/runner/work/MFC/MFC/src/simulation/m_data_output.fpp"
1281!$omp target update from(ib_markers%sf)
1282# 888 "/home/runner/work/MFC/MFC/src/simulation/m_data_output.fpp"
1283#endif
1284
1285 data_size = (m + 1)*(n + 1)*(p + 1)
1286 m_mok = int(m_glb + 1, mpi_offset_kind)
1287 n_mok = int(n_glb + 1, mpi_offset_kind)
1288 p_mok = int(p_glb + 1, mpi_offset_kind)
1289 wp_mok = int(storage_size(0._stp)/8, mpi_offset_kind)
1290 mok = int(1._wp, mpi_offset_kind)
1291
1292 write (file_loc, '(A)') 'ib.dat'
1293 file_loc = trim(case_dir) // '/restart_data' // trim(mpiiofs) // trim(file_loc)
1294 call mpi_file_open(mpi_comm_world, file_loc, ior(mpi_mode_wronly, mpi_mode_create), mpi_info_int, ifile, ierr)
1295
1296 var_mok = int(sys_size + 1, mpi_offset_kind)
1297 disp = m_mok*max(mok, n_mok)*max(mok, p_mok)*wp_mok*(var_mok - 1 + int(time_step/t_step_save))
1298 if (time_step == 0) disp = 0
1299
1300 call mpi_file_set_view(ifile, disp, mpi_integer, mpi_io_ib_data%view, 'native', mpi_info_int, ierr)
1301 call mpi_file_write_all(ifile, mpi_io_ib_data%var%sf, data_size, mpi_integer, status, ierr)
1302 call mpi_file_close(ifile, ierr)
1303#endif
1304
1305 end subroutine s_write_parallel_ib_data
1306
1307 !> Dispatch immersed boundary data output to the serial or parallel writer
1308 subroutine s_write_ib_data_file(time_step)
1309
1310 integer, intent(in) :: time_step
1311
1312
1313# 917 "/home/runner/work/MFC/MFC/src/simulation/m_data_output.fpp"
1314#if defined(MFC_OpenACC)
1315# 917 "/home/runner/work/MFC/MFC/src/simulation/m_data_output.fpp"
1316!$acc update host(patch_ib(1:num_ibs))
1317# 917 "/home/runner/work/MFC/MFC/src/simulation/m_data_output.fpp"
1318#elif defined(MFC_OpenMP)
1319# 917 "/home/runner/work/MFC/MFC/src/simulation/m_data_output.fpp"
1320!$omp target update from(patch_ib(1:num_ibs))
1321# 917 "/home/runner/work/MFC/MFC/src/simulation/m_data_output.fpp"
1322#endif
1323
1324 if (parallel_io) then
1325 call s_write_parallel_ib_data(time_step)
1326 else
1327 call s_write_serial_ib_data(time_step)
1328 end if
1329
1330 end subroutine s_write_ib_data_file
1331
1332 !> Writes the IB state information out to file
1333 subroutine s_write_parallel_ib_state(t_step)
1334
1335 integer, intent(in) :: t_step
1336
1337#ifdef MFC_MPI
1338 character(LEN=path_len + 2*name_len) :: file_loc
1339 integer(kind=MPI_OFFSET_KIND) :: disp
1340 integer(kind=MPI_OFFSET_KIND) :: WP_MOK
1341 integer :: ifile, ierr
1342 integer, dimension(MPI_STATUS_SIZE) :: status
1343 logical :: file_exist, dir_check
1344 integer :: i, ib_idx
1345 integer, parameter :: NFIELDS_PER_IB = 20
1346 real(wp) :: ib_buf(NFIELDS_PER_IB)
1347 integer :: file_unit
1348 character(len=10) :: t_step_string
1349
1350 ! Partition IBs across ranks round-robin style
1351 integer :: ib_start, ib_end, nibs_per_rank, remainder
1352
1353 wp_mok = int(storage_size(0._wp)/8, mpi_offset_kind)
1354
1355 if (file_per_process) then
1356 call s_int_to_str(t_step, t_step_string)
1357
1358 if (proc_rank == 0) then
1359 file_loc = trim(case_dir) // '/restart_data/lustre_' // trim(t_step_string)
1360 call s_create_directory(trim(file_loc))
1361 end if
1362 call s_mpi_barrier()
1364
1365 write (file_loc, '(A,I0,A,i7.7,A)') 'ib_state_', t_step, '_', proc_rank, '.dat'
1366 file_loc = trim(case_dir) // '/restart_data/lustre_' // trim(t_step_string) // '/' // trim(file_loc)
1367
1368 inquire (file=trim(file_loc), exist=file_exist)
1369 if (file_exist) then
1370 open (newunit=file_unit, file=trim(file_loc), form='unformatted', access='stream', status='replace')
1371 else
1372 open (newunit=file_unit, file=trim(file_loc), form='unformatted', access='stream', status='new')
1373 end if
1374
1375 write (file_unit) num_local_ibs
1376 do i = 1, num_local_ibs
1377 ib_idx = local_ib_patch_ids(i)
1378 ib_buf(1) = mytime
1379 ib_buf(2:4) = patch_ib(ib_idx)%force(1:3)
1380 ib_buf(5:7) = patch_ib(ib_idx)%torque(1:3)
1381 ib_buf(8:10) = patch_ib(ib_idx)%vel(1:3)
1382 ib_buf(11:13) = patch_ib(ib_idx)%angular_vel(1:3)
1383 ib_buf(14:16) = patch_ib(ib_idx)%angles(1:3)
1384 ib_buf(17) = patch_ib(ib_idx)%x_centroid
1385 ib_buf(18) = patch_ib(ib_idx)%y_centroid
1386 ib_buf(19) = patch_ib(ib_idx)%z_centroid
1387 ib_buf(20) = patch_ib(ib_idx)%radius
1388
1389 write (file_unit) patch_ib(ib_idx)%gbl_patch_id
1390 write (file_unit) ib_buf
1391 end do
1392
1393 close (file_unit)
1394 else
1395 if (proc_rank == 0) then
1396 call s_create_directory(trim(case_dir) // '/restart_data')
1397 end if
1398 call s_mpi_barrier()
1399
1400 write (file_loc, '(A,I0,A)') '/restart_data/ib_state_', t_step, '.dat'
1401 file_loc = trim(case_dir) // trim(file_loc)
1402
1403 inquire (file=trim(file_loc), exist=file_exist)
1404 if (file_exist .and. proc_rank == 0) then
1405 call mpi_file_delete(file_loc, mpi_info_int, ierr)
1406 end if
1407 call s_mpi_barrier()
1408
1409 call mpi_file_open(mpi_comm_world, file_loc, ior(mpi_mode_wronly, mpi_mode_create), mpi_info_int, ifile, ierr)
1410
1411 do i = 1, num_local_ibs
1412 ib_idx = local_ib_patch_ids(i)
1413 ib_buf(1) = mytime
1414 ib_buf(2:4) = patch_ib(ib_idx)%force(1:3)
1415 ib_buf(5:7) = patch_ib(ib_idx)%torque(1:3)
1416 ib_buf(8:10) = patch_ib(ib_idx)%vel(1:3)
1417 ib_buf(11:13) = patch_ib(ib_idx)%angular_vel(1:3)
1418 ib_buf(14:16) = patch_ib(ib_idx)%angles(1:3)
1419 ib_buf(17) = patch_ib(ib_idx)%x_centroid
1420 ib_buf(18) = patch_ib(ib_idx)%y_centroid
1421 ib_buf(19) = patch_ib(ib_idx)%z_centroid
1422 ib_buf(20) = patch_ib(ib_idx)%radius
1423
1424 ! Global IB index determines position in file
1425 disp = int(patch_ib(ib_idx)%gbl_patch_id - 1, mpi_offset_kind)*int(nfields_per_ib, mpi_offset_kind)*wp_mok
1426
1427 call mpi_file_write_at(ifile, disp, ib_buf, nfields_per_ib, mpi_p, status, ierr)
1428 end do
1429
1430 call mpi_file_close(ifile, ierr)
1431 end if
1432#endif
1433
1434 end subroutine s_write_parallel_ib_state
1435
1436 !> Write IB state data to a per-timestep serial (unformatted) file
1437 subroutine s_write_serial_ib_state(t_step)
1438
1439 integer, intent(in) :: t_step
1440 character(LEN=path_len + 2*name_len) :: file_loc
1441 integer :: i, ios, file_unit
1442 integer, parameter :: NFIELDS_PER_IB = 20
1443 real(wp) :: ib_buf(NFIELDS_PER_IB)
1444
1445 call s_create_directory(trim(case_dir) // '/restart_data')
1446
1447 write (file_loc, '(A,I0,A)') '/restart_data/ib_state_', t_step, '.dat'
1448 file_loc = trim(case_dir) // trim(file_loc)
1449
1450 open (newunit=file_unit, file=trim(file_loc), form='unformatted', access='stream', status='replace', iostat=ios)
1451 if (ios /= 0) call s_mpi_abort('Cannot open IB state output file: ' // trim(file_loc))
1452
1453 do i = 1, num_ibs
1454 ib_buf(1) = mytime
1455 ib_buf(2:4) = patch_ib(i)%force(1:3)
1456 ib_buf(5:7) = patch_ib(i)%torque(1:3)
1457 ib_buf(8:10) = patch_ib(i)%vel(1:3)
1458 ib_buf(11:13) = patch_ib(i)%angular_vel(1:3)
1459 ib_buf(14:16) = patch_ib(i)%angles(1:3)
1460 ib_buf(17) = patch_ib(i)%x_centroid
1461 ib_buf(18) = patch_ib(i)%y_centroid
1462 ib_buf(19) = patch_ib(i)%z_centroid
1463 ib_buf(20) = patch_ib(i)%radius
1464
1465 write (file_unit) ib_buf
1466 end do
1467
1468 close (file_unit)
1469
1470 end subroutine s_write_serial_ib_state
1471
1472 !> @brief Writes IB state records to restart_data/ib_state.dat. Must be called only on rank 0.
1473 impure subroutine s_write_ib_state_file(time_step)
1474
1475 integer, intent(in) :: time_step
1476
1477
1478# 1072 "/home/runner/work/MFC/MFC/src/simulation/m_data_output.fpp"
1479#if defined(MFC_OpenACC)
1480# 1072 "/home/runner/work/MFC/MFC/src/simulation/m_data_output.fpp"
1481!$acc update host(patch_ib(1:num_ibs))
1482# 1072 "/home/runner/work/MFC/MFC/src/simulation/m_data_output.fpp"
1483#elif defined(MFC_OpenMP)
1484# 1072 "/home/runner/work/MFC/MFC/src/simulation/m_data_output.fpp"
1485!$omp target update from(patch_ib(1:num_ibs))
1486# 1072 "/home/runner/work/MFC/MFC/src/simulation/m_data_output.fpp"
1487#endif
1488
1489 if (parallel_io) then
1490 call s_write_parallel_ib_state(time_step)
1491 else
1492 call s_write_serial_ib_state(time_step)
1493 end if
1494
1495 end subroutine s_write_ib_state_file
1496
1497 !> Write center-of-mass data at the current time step
1498 impure subroutine s_write_com_files(t_step, c_mass_in)
1499
1500 integer, intent(in) :: t_step
1501 real(wp), dimension(num_fluids, 5), intent(in) :: c_mass_in
1502 integer :: i !< Generic loop iterator
1503 real(wp) :: nondim_time !< Non-dimensional time
1504
1505 if (t_step_old /= dflt_int) then
1506 nondim_time = real(t_step + t_step_old, wp)*dt
1507 else
1508 nondim_time = real(t_step, wp)*dt
1509 end if
1510
1511 if (proc_rank == 0) then
1512 if (n == 0) then
1513 do i = 1, num_fluids
1514 write (i + 120, '(6X,4F24.12)') nondim_time, c_mass_in(i, 1), c_mass_in(i, 2), c_mass_in(i, 5)
1515 end do
1516 else if (p == 0) then
1517 do i = 1, num_fluids
1518 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)
1519 end do
1520 else
1521 do i = 1, num_fluids
1522 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, &
1523 & 4), c_mass_in(i, 5)
1524 end do
1525 end if
1526 end if
1527
1528 end subroutine s_write_com_files
1529
1530 !> Write flow probe data at the current time step
1531 impure subroutine s_write_probe_files(t_step, q_cons_vf, accel_mag)
1532
1533 integer, intent(in) :: t_step
1534 type(scalar_field), dimension(sys_size), intent(in) :: q_cons_vf
1535 real(wp), dimension(0:m,0:n,0:p), intent(in) :: accel_mag
1536 real(wp), dimension(-1:m) :: distx
1537 real(wp), dimension(-1:n) :: disty
1538 real(wp), dimension(-1:p) :: distz
1539
1540 ! The cell-averaged partial densities, density, velocity, pressure, volume fractions, specific heat ratio function, liquid
1541 ! stiffness function, and sound speed.
1542 real(wp) :: lit_gamma, nbub
1543 real(wp) :: rho
1544 real(wp), dimension(num_vels) :: vel
1545 real(wp) :: pres
1546 real(wp) :: ptilde
1547 real(wp) :: ptot
1548 real(wp) :: alf
1549 real(wp) :: alfgr
1550 real(wp), dimension(num_fluids) :: alpha
1551 real(wp) :: gamma
1552 real(wp) :: pi_inf
1553 real(wp) :: qv
1554 real(wp) :: c
1555 real(wp) :: m00, m10, m01, m20, m02
1556 real(wp) :: varr, varv
1557 real(wp), dimension(Nb) :: nr, r, nrdot, rdot
1558 real(wp) :: nr3
1559 real(wp) :: accel
1560 real(wp) :: int_pres
1561 real(wp) :: max_pres
1562 real(wp), dimension(2) :: re
1563 real(wp), dimension(6) :: tau_e
1564 real(wp) :: g_local
1565 real(wp) :: dyn_p, t
1566 real(wp) :: damage_state
1567 integer :: i, j, k, l, s, d !< Generic loop iterator
1568 real(wp) :: nondim_time !< Non-dimensional time
1569 real(wp) :: tmp !< Temporary variable to store quantity for mpi_allreduce
1570 integer :: npts !< Number of included integral points
1571 real(wp) :: rad, thickness !< For integral quantities
1572 logical :: trigger !< For integral quantities
1573 real(wp) :: rhoyks(1:num_species)
1574
1575 t = dflt_t_guess
1576
1577 if (time_stepper == 23) then
1578 nondim_time = mytime
1579 else
1580 if (t_step_old /= dflt_int) then
1581 nondim_time = real(t_step + t_step_old, wp)*dt
1582 else
1583 nondim_time = real(t_step, wp)*dt
1584 end if
1585 end if
1586
1587 do i = 1, num_probes
1588 rho = 0._wp
1589 do s = 1, num_vels
1590 vel(s) = 0._wp
1591 end do
1592 pres = 0._wp
1593 gamma = 0._wp
1594 pi_inf = 0._wp
1595 qv = 0._wp
1596 c = 0._wp
1597 accel = 0._wp
1598 nr = 0._wp; r = 0._wp
1599 nrdot = 0._wp; rdot = 0._wp
1600 nbub = 0._wp
1601 m00 = 0._wp
1602 m10 = 0._wp
1603 m01 = 0._wp
1604 m20 = 0._wp
1605 m02 = 0._wp
1606 varr = 0._wp; varv = 0._wp
1607 alf = 0._wp
1608 do s = 1, (num_dims*(num_dims + 1))/2
1609 tau_e(s) = 0._wp
1610 end do
1611 damage_state = 0._wp
1612
1613 if (n == 0) then
1614 if ((probe(i)%x >= x_cb(-1)) .and. (probe(i)%x <= x_cb(m))) then
1615 do s = -1, m
1616 distx(s) = x_cb(s) - probe(i)%x
1617 if (distx(s) < 0._wp) distx(s) = 1000._wp
1618 end do
1619 j = minloc(distx, 1)
1620 if (j == 1) j = 2 ! Pick first point if probe is at edge
1621 k = 0
1622 l = 0
1623
1624 if (chemistry) then
1625 do d = 1, num_species
1626 rhoyks(d) = q_cons_vf(eqn_idx%species%beg + d - 1)%sf(j - 2, k, l)
1627 end do
1628 end if
1629
1630 ! Computing/Sharing necessary state variables
1631 if (elasticity) then
1632 call s_convert_to_mixture_variables(q_cons_vf, j - 2, k, l, rho, gamma, pi_inf, qv, re, g_local, &
1633 & fluid_pp(:)%G)
1634 else
1635 call s_convert_to_mixture_variables(q_cons_vf, j - 2, k, l, rho, gamma, pi_inf, qv)
1636 end if
1637 do s = 1, num_vels
1638 vel(s) = q_cons_vf(eqn_idx%cont%end + s)%sf(j - 2, k, l)/rho
1639 end do
1640
1641 dyn_p = 0.5_wp*rho*dot_product(vel, vel)
1642
1643 if (elasticity) then
1644 if (cont_damage) then
1645 damage_state = q_cons_vf(eqn_idx%damage)%sf(j - 2, k, l)
1646 g_local = g_local*max((1._wp - damage_state), 0._wp)
1647 end if
1648
1649 call s_compute_pressure(q_cons_vf(1)%sf(j - 2, k, l), q_cons_vf(eqn_idx%alf)%sf(j - 2, k, l), dyn_p, &
1650 & pi_inf, gamma, rho, qv, rhoyks(:), pres, t, &
1651 & q_cons_vf(eqn_idx%stress%beg)%sf(j - 2, k, l), &
1652 & q_cons_vf(eqn_idx%mom%beg)%sf(j - 2, k, l), g_local)
1653 else
1654 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), &
1655 & dyn_p, pi_inf, gamma, rho, qv, rhoyks, pres, t)
1656 end if
1657
1658 if (model_eqns == model_eqns_4eq) then
1659 lit_gamma = gammas(1)
1660 else if (elasticity) then
1661 tau_e(1) = q_cons_vf(eqn_idx%stress%end)%sf(j - 2, k, l)/rho
1662 end if
1663
1664 if (bubbles_euler) then
1665 alf = q_cons_vf(eqn_idx%alf)%sf(j - 2, k, l)
1666 if (num_fluids == 3) then
1667 alfgr = q_cons_vf(eqn_idx%alf - 1)%sf(j - 2, k, l)
1668 end if
1669 do s = 1, nb
1670 nr(s) = q_cons_vf(qbmm_idx%rs(s))%sf(j - 2, k, l)
1671 nrdot(s) = q_cons_vf(qbmm_idx%vs(s))%sf(j - 2, k, l)
1672 end do
1673
1674 if (adv_n) then
1675 nbub = q_cons_vf(eqn_idx%n)%sf(j - 2, k, l)
1676 else
1677 nr3 = 0._wp
1678 do s = 1, nb
1679 nr3 = nr3 + weight(s)*(nr(s)**3._wp)
1680 end do
1681
1682 nbub = sqrt((4._wp*pi/3._wp)*nr3/alf)
1683 end if
1684#ifdef MFC_DEBUG
1685 print *, 'In probe, nbub: ', nbub
1686#endif
1687 if (qbmm) then
1688 m00 = q_cons_vf(qbmm_idx%moms(1, 1))%sf(j - 2, k, l)/nbub
1689 m10 = q_cons_vf(qbmm_idx%moms(1, 2))%sf(j - 2, k, l)/nbub
1690 m01 = q_cons_vf(qbmm_idx%moms(1, 3))%sf(j - 2, k, l)/nbub
1691 m20 = q_cons_vf(qbmm_idx%moms(1, 4))%sf(j - 2, k, l)/nbub
1692 m02 = q_cons_vf(qbmm_idx%moms(1, 6))%sf(j - 2, k, l)/nbub
1693
1694 m10 = m10/m00
1695 m01 = m01/m00
1696 m20 = m20/m00
1697 m02 = m02/m00
1698
1699 varr = m20 - m10**2._wp
1700 varv = m02 - m01**2._wp
1701 end if
1702 r(:) = nr(:)/nbub
1703 rdot(:) = nrdot(:)/nbub
1704
1705 ptilde = ptil(j - 2, k, l)
1706 ptot = pres - ptilde
1707 end if
1708
1709 ! Compute mixture sound Speed
1710 call s_compute_speed_of_sound(pres, rho, gamma, pi_inf, ((gamma + 1._wp)*pres + pi_inf)/rho, alpha, 0._wp, &
1711 & 0._wp, c, qv)
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
1745 dyn_p = 0.5_wp*rho*dot_product(vel, vel)
1746
1747 if (elasticity) then
1748 if (cont_damage) then
1749 damage_state = q_cons_vf(eqn_idx%damage)%sf(j - 2, k - 2, l)
1750 g_local = g_local*max((1._wp - damage_state), 0._wp)
1751 end if
1752
1753 call s_compute_pressure(q_cons_vf(1)%sf(j - 2, k - 2, l), q_cons_vf(eqn_idx%alf)%sf(j - 2, k - 2, l), &
1754 & dyn_p, pi_inf, gamma, rho, qv, rhoyks, pres, t, &
1755 & q_cons_vf(eqn_idx%stress%beg)%sf(j - 2, k - 2, l), &
1756 & q_cons_vf(eqn_idx%mom%beg)%sf(j - 2, k - 2, l), g_local)
1757 else
1758 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, &
1759 & k - 2, l), dyn_p, pi_inf, gamma, rho, qv, rhoyks, pres, t)
1760 end if
1761
1762 if (model_eqns == model_eqns_4eq) then
1763 lit_gamma = gs_min(1)
1764 else if (elasticity) then
1765 do s = 1, 3
1766 tau_e(s) = q_cons_vf(s)%sf(j - 2, k - 2, l)/rho
1767 end do
1768 end if
1769
1770 if (bubbles_euler) then
1771 alf = q_cons_vf(eqn_idx%alf)%sf(j - 2, k - 2, l)
1772 do s = 1, nb
1773 nr(s) = q_cons_vf(qbmm_idx%rs(s))%sf(j - 2, k - 2, l)
1774 nrdot(s) = q_cons_vf(qbmm_idx%vs(s))%sf(j - 2, k - 2, l)
1775 end do
1776
1777 if (adv_n) then
1778 nbub = q_cons_vf(eqn_idx%n)%sf(j - 2, k - 2, l)
1779 else
1780 nr3 = 0._wp
1781 do s = 1, nb
1782 nr3 = nr3 + weight(s)*(nr(s)**3._wp)
1783 end do
1784
1785 nbub = sqrt((4._wp*pi/3._wp)*nr3/alf)
1786 end if
1787
1788 r(:) = nr(:)/nbub
1789 rdot(:) = nrdot(:)/nbub
1790 end if
1791 ! Compute mixture sound speed
1792 call s_compute_speed_of_sound(pres, rho, gamma, pi_inf, ((gamma + 1._wp)*pres + pi_inf)/rho, alpha, &
1793 & 0._wp, 0._wp, c, qv)
1794 end if
1795 end if
1796 else
1797 if ((probe(i)%x >= x_cb(-1)) .and. (probe(i)%x <= x_cb(m))) then
1798 if ((probe(i)%y >= y_cb(-1)) .and. (probe(i)%y <= y_cb(n))) then
1799 if ((probe(i)%z >= z_cb(-1)) .and. (probe(i)%z <= z_cb(p))) then
1800 do s = -1, m
1801 distx(s) = x_cb(s) - probe(i)%x
1802 if (distx(s) < 0._wp) distx(s) = 1000._wp
1803 end do
1804 do s = -1, n
1805 disty(s) = y_cb(s) - probe(i)%y
1806 if (disty(s) < 0._wp) disty(s) = 1000._wp
1807 end do
1808 do s = -1, p
1809 distz(s) = z_cb(s) - probe(i)%z
1810 if (distz(s) < 0._wp) distz(s) = 1000._wp
1811 end do
1812 j = minloc(distx, 1)
1813 k = minloc(disty, 1)
1814 l = minloc(distz, 1)
1815 if (j == 1) j = 2 ! Pick first point if probe is at edge
1816 if (k == 1) k = 2 ! Pick first point if probe is at edge
1817 if (l == 1) l = 2 ! Pick first point if probe is at edge
1818
1819 ! Computing/Sharing necessary state variables
1820 call s_convert_to_mixture_variables(q_cons_vf, j - 2, k - 2, l - 2, rho, gamma, pi_inf, qv, re, &
1821 & g_local, fluid_pp(:)%G)
1822 do s = 1, num_vels
1823 vel(s) = q_cons_vf(eqn_idx%cont%end + s)%sf(j - 2, k - 2, l - 2)/rho
1824 end do
1825
1826 dyn_p = 0.5_wp*rho*dot_product(vel, vel)
1827
1828 if (chemistry) then
1829 do d = 1, num_species
1830 rhoyks(d) = q_cons_vf(eqn_idx%species%beg + d - 1)%sf(j - 2, k - 2, l - 2)
1831 end do
1832 end if
1833
1834 if (elasticity) then
1835 if (cont_damage) then
1836 damage_state = q_cons_vf(eqn_idx%damage)%sf(j - 2, k - 2, l - 2)
1837 g_local = g_local*max((1._wp - damage_state), 0._wp)
1838 end if
1839
1840 call s_compute_pressure(q_cons_vf(1)%sf(j - 2, k - 2, l - 2), q_cons_vf(eqn_idx%alf)%sf(j - 2, &
1841 & k - 2, l - 2), dyn_p, pi_inf, gamma, rho, qv, rhoyks, pres, t, &
1842 & q_cons_vf(eqn_idx%stress%beg)%sf(j - 2, k - 2, l - 2), &
1843 & q_cons_vf(eqn_idx%mom%beg)%sf(j - 2, k - 2, l - 2), g_local)
1844 else
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)
1848 end if
1849
1850 ! Compute mixture sound speed
1851 call s_compute_speed_of_sound(pres, rho, gamma, pi_inf, ((gamma + 1._wp)*pres + pi_inf)/rho, alpha, &
1852 & 0._wp, 0._wp, c, qv)
1853
1854 accel = accel_mag(j - 2, k - 2, l - 2)
1855 end if
1856 end if
1857 end if
1858 end if
1859 if (num_procs > 1) then
1860# 1446 "/home/runner/work/MFC/MFC/src/simulation/m_data_output.fpp"
1861 tmp = rho
1862 call s_mpi_allreduce_sum(tmp, rho)
1863# 1446 "/home/runner/work/MFC/MFC/src/simulation/m_data_output.fpp"
1864 tmp = pres
1865 call s_mpi_allreduce_sum(tmp, pres)
1866# 1446 "/home/runner/work/MFC/MFC/src/simulation/m_data_output.fpp"
1867 tmp = gamma
1868 call s_mpi_allreduce_sum(tmp, gamma)
1869# 1446 "/home/runner/work/MFC/MFC/src/simulation/m_data_output.fpp"
1870 tmp = pi_inf
1871 call s_mpi_allreduce_sum(tmp, pi_inf)
1872# 1446 "/home/runner/work/MFC/MFC/src/simulation/m_data_output.fpp"
1873 tmp = qv
1874 call s_mpi_allreduce_sum(tmp, qv)
1875# 1446 "/home/runner/work/MFC/MFC/src/simulation/m_data_output.fpp"
1876 tmp = c
1877 call s_mpi_allreduce_sum(tmp, c)
1878# 1446 "/home/runner/work/MFC/MFC/src/simulation/m_data_output.fpp"
1879 tmp = accel
1880 call s_mpi_allreduce_sum(tmp, accel)
1881# 1449 "/home/runner/work/MFC/MFC/src/simulation/m_data_output.fpp"
1882
1883 do s = 1, num_vels
1884 tmp = vel(s)
1885 call s_mpi_allreduce_sum(tmp, vel(s))
1886 end do
1887
1888 if (bubbles_euler) then
1889# 1457 "/home/runner/work/MFC/MFC/src/simulation/m_data_output.fpp"
1890 tmp = alf
1891 call s_mpi_allreduce_sum(tmp, alf)
1892# 1457 "/home/runner/work/MFC/MFC/src/simulation/m_data_output.fpp"
1893 tmp = alfgr
1894 call s_mpi_allreduce_sum(tmp, alfgr)
1895# 1457 "/home/runner/work/MFC/MFC/src/simulation/m_data_output.fpp"
1896 tmp = nbub
1897 call s_mpi_allreduce_sum(tmp, nbub)
1898# 1457 "/home/runner/work/MFC/MFC/src/simulation/m_data_output.fpp"
1899 tmp = nr(1)
1900 call s_mpi_allreduce_sum(tmp, nr(1))
1901# 1457 "/home/runner/work/MFC/MFC/src/simulation/m_data_output.fpp"
1902 tmp = nrdot(1)
1903 call s_mpi_allreduce_sum(tmp, nrdot(1))
1904# 1457 "/home/runner/work/MFC/MFC/src/simulation/m_data_output.fpp"
1905 tmp = m00
1906 call s_mpi_allreduce_sum(tmp, m00)
1907# 1457 "/home/runner/work/MFC/MFC/src/simulation/m_data_output.fpp"
1908 tmp = r(1)
1909 call s_mpi_allreduce_sum(tmp, r(1))
1910# 1457 "/home/runner/work/MFC/MFC/src/simulation/m_data_output.fpp"
1911 tmp = rdot(1)
1912 call s_mpi_allreduce_sum(tmp, rdot(1))
1913# 1457 "/home/runner/work/MFC/MFC/src/simulation/m_data_output.fpp"
1914 tmp = ptilde
1915 call s_mpi_allreduce_sum(tmp, ptilde)
1916# 1457 "/home/runner/work/MFC/MFC/src/simulation/m_data_output.fpp"
1917 tmp = ptot
1918 call s_mpi_allreduce_sum(tmp, ptot)
1919# 1460 "/home/runner/work/MFC/MFC/src/simulation/m_data_output.fpp"
1920
1921 if (qbmm) then
1922# 1463 "/home/runner/work/MFC/MFC/src/simulation/m_data_output.fpp"
1923 tmp = varr
1924 call s_mpi_allreduce_sum(tmp, varr)
1925# 1463 "/home/runner/work/MFC/MFC/src/simulation/m_data_output.fpp"
1926 tmp = varv
1927 call s_mpi_allreduce_sum(tmp, varv)
1928# 1463 "/home/runner/work/MFC/MFC/src/simulation/m_data_output.fpp"
1929 tmp = m10
1930 call s_mpi_allreduce_sum(tmp, m10)
1931# 1463 "/home/runner/work/MFC/MFC/src/simulation/m_data_output.fpp"
1932 tmp = m01
1933 call s_mpi_allreduce_sum(tmp, m01)
1934# 1463 "/home/runner/work/MFC/MFC/src/simulation/m_data_output.fpp"
1935 tmp = m20
1936 call s_mpi_allreduce_sum(tmp, m20)
1937# 1463 "/home/runner/work/MFC/MFC/src/simulation/m_data_output.fpp"
1938 tmp = m02
1939 call s_mpi_allreduce_sum(tmp, m02)
1940# 1466 "/home/runner/work/MFC/MFC/src/simulation/m_data_output.fpp"
1941 end if
1942 end if
1943
1944 if (elasticity) then
1945 do s = 1, (num_dims*(num_dims + 1))/2
1946 tmp = tau_e(s)
1947 call s_mpi_allreduce_sum(tmp, tau_e(s))
1948 end do
1949 end if
1950
1951 if (cont_damage) then
1952 tmp = damage_state
1953 call s_mpi_allreduce_sum(tmp, damage_state)
1954 end if
1955 end if
1956 if (proc_rank == 0) then
1957 if (n == 0) then
1958 if (bubbles_euler .and. (num_fluids <= 2)) then
1959 if (qbmm) then
1960 write (i + 30, '(6x,f12.6,14f28.16)') nondim_time, rho, vel(1), pres, alf, r(1), rdot(1), nr(1), &
1961 & nrdot(1), varr, varv, m10, m01, m20, m02
1962 else
1963 write (i + 30, '(6x,f12.6,8f24.8)') nondim_time, rho, vel(1), pres, alf, r(1), rdot(1), nr(1), nrdot(1)
1964 ! ptilde, & ptot
1965 end if
1966 else if (bubbles_euler .and. (num_fluids == 3)) then
1967 write (i + 30, &
1968 & '(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)') &
1969 & nondim_time, rho, vel(1), pres, alf, alfgr, nr(1), nrdot(1), r(1), rdot(1), ptilde, ptot
1970 else if (bubbles_euler .and. num_fluids == 4) then
1971 write (i + 30, &
1972 & '(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)') &
1973 & 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, &
1974 & 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), &
1975 & 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), &
1976 & q_cons_vf(10)%sf(j - 2, 0, 0), nbub, r(1), rdot(1)
1977 else
1978 write (i + 30, '(6X,F12.6,F24.8,F24.8,F24.8)') nondim_time, rho, vel(1), pres
1979 end if
1980 else if (p == 0) then
1981 if (bubbles_euler) then
1982# 1508 "/home/runner/work/MFC/MFC/src/simulation/m_data_output.fpp"
1983 write (i + 30, '(6X,10F24.8)') nondim_time, rho, vel(1), vel(2), pres, alf, nr(1), nrdot(1), r(1), &
1984 & rdot(1)
1985# 1511 "/home/runner/work/MFC/MFC/src/simulation/m_data_output.fpp"
1986 else if (elasticity) then
1987# 1513 "/home/runner/work/MFC/MFC/src/simulation/m_data_output.fpp"
1988 write (i + 30, '(6X,F12.6,F24.8,F24.8,F24.8,F24.8,' // 'F24.8,F24.8,F24.8)') nondim_time, rho, &
1989 & vel(1), vel(2), pres, tau_e(1), tau_e(2), tau_e(3)
1990# 1516 "/home/runner/work/MFC/MFC/src/simulation/m_data_output.fpp"
1991 else
1992 write (i + 30, '(6X,F12.6,F24.8,F24.8,F24.8)') nondim_time, rho, vel(1), pres
1993 print *, 'time =', nondim_time, 'rho =', rho, 'pres =', pres
1994 end if
1995 else
1996# 1522 "/home/runner/work/MFC/MFC/src/simulation/m_data_output.fpp"
1997 write (i + 30, &
1998 & '(6X,F12.6,F24.8,F24.8,F24.8,F24.8,' // 'F24.8,F24.8,F24.8,F24.8,F24.8,' // 'F24.8)') &
1999 & nondim_time, rho, vel(1), vel(2), vel(3), pres, gamma, pi_inf, qv, c, accel
2000# 1526 "/home/runner/work/MFC/MFC/src/simulation/m_data_output.fpp"
2001 end if
2002 end if
2003 end do
2004
2005 if (integral_wrt .and. bubbles_euler) then
2006 if (n == 0) then
2007 do i = 1, num_integrals
2008 int_pres = 0._wp
2009 max_pres = 0._wp
2010 k = 0; l = 0
2011 npts = 0
2012 do j = 1, m
2013 pres = 0._wp
2014 do s = 1, num_vels
2015 vel(s) = 0._wp
2016 end do
2017 rho = 0._wp
2018 pres = 0._wp
2019 gamma = 0._wp
2020 pi_inf = 0._wp
2021 qv = 0._wp
2022
2023 if ((integral(i)%xmin <= x_cb(j)) .and. (integral(i)%xmax >= x_cb(j))) then
2024 npts = npts + 1
2025 call s_convert_to_mixture_variables(q_cons_vf, j, k, l, rho, gamma, pi_inf, qv, re)
2026 do s = 1, num_vels
2027 vel(s) = q_cons_vf(eqn_idx%cont%end + s)%sf(j, k, l)/rho
2028 end do
2029
2030 pres = ((q_cons_vf(eqn_idx%E)%sf(j, k, l) - 0.5_wp*(q_cons_vf(eqn_idx%mom%beg)%sf(j, k, &
2031 & l)**2._wp)/rho)/(1._wp - q_cons_vf(eqn_idx%alf)%sf(j, k, l)) - pi_inf - qv)/gamma
2032 int_pres = int_pres + (pres - 1._wp)**2._wp
2033 end if
2034 end do
2035 int_pres = sqrt(int_pres/(1._wp*npts))
2036
2037 if (num_procs > 1) then
2038 tmp = int_pres
2039 call s_mpi_allreduce_sum(tmp, int_pres)
2040 end if
2041
2042 if (proc_rank == 0) then
2043 if (bubbles_euler .and. (num_fluids <= 2)) then
2044 write (i + 70, '(6x,f12.6,f24.8)') nondim_time, int_pres
2045 end if
2046 end if
2047 end do
2048 else if (p == 0) then
2049 if (num_integrals /= 3) then
2050 call s_mpi_abort('Incorrect number of integrals')
2051 end if
2052
2053 rad = integral(1)%xmax
2054 thickness = integral(1)%xmin
2055
2056 do i = 1, num_integrals
2057 int_pres = 0._wp
2058 max_pres = 0._wp
2059 l = 0
2060 npts = 0
2061 do j = 1, m
2062 do k = 1, n
2063 trigger = .false.
2064 if (i == 1) then
2065 ! inner portion
2066 if (sqrt(x_cb(j)**2._wp + y_cb(k)**2._wp) < (rad - 0.5_wp*thickness)) trigger = .true.
2067 else if (i == 2) then
2068 ! net region
2069 if (sqrt(x_cb(j)**2._wp + y_cb(k)**2._wp) > (rad - 0.5_wp*thickness) .and. sqrt(x_cb(j)**2._wp &
2070 & + y_cb(k)**2._wp) < (rad + 0.5_wp*thickness)) trigger = .true.
2071 else if (i == 3) then
2072 ! everything else
2073 if (sqrt(x_cb(j)**2._wp + y_cb(k)**2._wp) > (rad + 0.5_wp*thickness)) trigger = .true.
2074 end if
2075
2076 pres = 0._wp
2077 do s = 1, num_vels
2078 vel(s) = 0._wp
2079 end do
2080 rho = 0._wp
2081 pres = 0._wp
2082 gamma = 0._wp
2083 pi_inf = 0._wp
2084 qv = 0._wp
2085
2086 if (trigger) then
2087 npts = npts + 1
2088 call s_convert_to_mixture_variables(q_cons_vf, j, k, l, rho, gamma, pi_inf, qv, re)
2089 do s = 1, num_vels
2090 vel(s) = q_cons_vf(eqn_idx%cont%end + s)%sf(j, k, l)/rho
2091 end do
2092
2093 pres = ((q_cons_vf(eqn_idx%E)%sf(j, k, l) - 0.5_wp*(q_cons_vf(eqn_idx%mom%beg)%sf(j, k, &
2094 & l)**2._wp)/rho)/(1._wp - q_cons_vf(eqn_idx%alf)%sf(j, k, l)) - pi_inf - qv)/gamma
2095 int_pres = int_pres + abs(pres - 1._wp)
2096 max_pres = max(max_pres, abs(pres - 1._wp))
2097 end if
2098 end do
2099 end do
2100
2101 if (npts > 0) then
2102 int_pres = int_pres/(1._wp*npts)
2103 else
2104 int_pres = 0._wp
2105 end if
2106
2107 if (num_procs > 1) then
2108 tmp = int_pres
2109 call s_mpi_allreduce_sum(tmp, int_pres)
2110
2111 tmp = max_pres
2112 call s_mpi_allreduce_max(tmp, max_pres)
2113 end if
2114
2115 if (proc_rank == 0) then
2116 if (bubbles_euler .and. (num_fluids <= 2)) then
2117 write (i + 70, '(6x,f12.6,f24.8,f24.8)') nondim_time, int_pres, max_pres
2118 end if
2119 end if
2120 end do
2121 end if
2122 end if
2123
2124 end subroutine s_write_probe_files
2125
2126 !> Write footer with stability criteria extrema and run-time to the information file, then close it
2128
2129 real(wp) :: run_time !< Run-time of the simulation
2130
2131 write (3, '(A)') ' '
2132 write (3, '(A)') ''
2133
2134 write (3, '(A,F9.6)') 'ICFL Max: ', icfl_max
2135 if (surface_tension) write (3, '(A,F9.6)') 'CCFL Max: ', ccfl_max
2136 if (viscous) write (3, '(A,F9.6)') 'VCFL Max: ', vcfl_max
2137 if (viscous) write (3, '(A,ES16.6)') 'Rc Min: ', rc_min
2138
2139 call cpu_time(run_time)
2140
2141 write (3, '(A)') ''
2142 write (3, '(A,I0,A)') 'Run-time: ', int(anint(run_time)), 's'
2143 write (3, '(A)') ' '
2144 close (3)
2145
2147
2148 !> Closes communication files
2149 impure subroutine s_close_com_files()
2150
2151 integer :: i !< Generic loop iterator
2152
2153 do i = 1, num_fluids
2154 close (i + 120)
2155 end do
2156
2157 end subroutine s_close_com_files
2158
2159 !> Closes probe files
2160 impure subroutine s_close_probe_files
2161
2162 integer :: i !< Generic loop iterator
2163
2164 do i = 1, num_probes
2165 close (i + 30)
2166 end do
2167
2168 end subroutine s_close_probe_files
2169
2170 !> Initialize the data output module
2172
2173 integer :: i, m_ds, n_ds, p_ds
2174
2175 if (run_time_info) then
2176 icfl_max = 0._wp
2177 if (surface_tension) then
2178 ccfl_max = 0._wp
2179 end if
2180 if (viscous) then
2181 vcfl_max = 0._wp
2182 rc_min = 1.e12_wp
2183 end if
2184 end if
2185
2186 if (probe_wrt) then
2187#ifdef MFC_DEBUG
2188# 1712 "/home/runner/work/MFC/MFC/src/simulation/m_data_output.fpp"
2189 block
2190# 1712 "/home/runner/work/MFC/MFC/src/simulation/m_data_output.fpp"
2191 use iso_fortran_env, only: output_unit
2192# 1712 "/home/runner/work/MFC/MFC/src/simulation/m_data_output.fpp"
2193
2194# 1712 "/home/runner/work/MFC/MFC/src/simulation/m_data_output.fpp"
2195 print *, 'm_data_output.fpp:1712: ', '@:ALLOCATE(c_mass(num_fluids,5))'
2196# 1712 "/home/runner/work/MFC/MFC/src/simulation/m_data_output.fpp"
2197
2198# 1712 "/home/runner/work/MFC/MFC/src/simulation/m_data_output.fpp"
2199 call flush (output_unit)
2200# 1712 "/home/runner/work/MFC/MFC/src/simulation/m_data_output.fpp"
2201 end block
2202# 1712 "/home/runner/work/MFC/MFC/src/simulation/m_data_output.fpp"
2203#endif
2204# 1712 "/home/runner/work/MFC/MFC/src/simulation/m_data_output.fpp"
2205 allocate (c_mass(num_fluids,5))
2206# 1712 "/home/runner/work/MFC/MFC/src/simulation/m_data_output.fpp"
2207
2208# 1712 "/home/runner/work/MFC/MFC/src/simulation/m_data_output.fpp"
2209
2210# 1712 "/home/runner/work/MFC/MFC/src/simulation/m_data_output.fpp"
2211#if defined(MFC_OpenACC)
2212# 1712 "/home/runner/work/MFC/MFC/src/simulation/m_data_output.fpp"
2213!$acc enter data create(c_mass)
2214# 1712 "/home/runner/work/MFC/MFC/src/simulation/m_data_output.fpp"
2215#elif defined(MFC_OpenMP)
2216# 1712 "/home/runner/work/MFC/MFC/src/simulation/m_data_output.fpp"
2217!$omp target enter data map(always,alloc:c_mass)
2218# 1712 "/home/runner/work/MFC/MFC/src/simulation/m_data_output.fpp"
2219#endif
2220 end if
2221
2222 if (down_sample) then
2223 m_ds = int((m + 1)/3) - 1
2224 n_ds = int((n + 1)/3) - 1
2225 p_ds = int((p + 1)/3) - 1
2226
2227 allocate (q_cons_temp_ds(1:sys_size))
2228 do i = 1, sys_size
2229 allocate (q_cons_temp_ds(i)%sf(-1:m_ds + 1,-1:n_ds + 1,-1:p_ds + 1))
2230 end do
2231 end if
2232
2233 end subroutine s_initialize_data_output_module
2234
2235 !> Module deallocation and/or disassociation procedures
2237
2238 integer :: i
2239
2240 if (probe_wrt) then
2241#ifdef MFC_DEBUG
2242# 1734 "/home/runner/work/MFC/MFC/src/simulation/m_data_output.fpp"
2243 block
2244# 1734 "/home/runner/work/MFC/MFC/src/simulation/m_data_output.fpp"
2245 use iso_fortran_env, only: output_unit
2246# 1734 "/home/runner/work/MFC/MFC/src/simulation/m_data_output.fpp"
2247
2248# 1734 "/home/runner/work/MFC/MFC/src/simulation/m_data_output.fpp"
2249 print *, 'm_data_output.fpp:1734: ', '@:DEALLOCATE(c_mass)'
2250# 1734 "/home/runner/work/MFC/MFC/src/simulation/m_data_output.fpp"
2251
2252# 1734 "/home/runner/work/MFC/MFC/src/simulation/m_data_output.fpp"
2253 call flush (output_unit)
2254# 1734 "/home/runner/work/MFC/MFC/src/simulation/m_data_output.fpp"
2255 end block
2256# 1734 "/home/runner/work/MFC/MFC/src/simulation/m_data_output.fpp"
2257#endif
2258# 1734 "/home/runner/work/MFC/MFC/src/simulation/m_data_output.fpp"
2259
2260# 1734 "/home/runner/work/MFC/MFC/src/simulation/m_data_output.fpp"
2261#if defined(MFC_OpenACC)
2262# 1734 "/home/runner/work/MFC/MFC/src/simulation/m_data_output.fpp"
2263!$acc exit data delete(c_mass)
2264# 1734 "/home/runner/work/MFC/MFC/src/simulation/m_data_output.fpp"
2265#elif defined(MFC_OpenMP)
2266# 1734 "/home/runner/work/MFC/MFC/src/simulation/m_data_output.fpp"
2267!$omp target exit data map(release:c_mass)
2268# 1734 "/home/runner/work/MFC/MFC/src/simulation/m_data_output.fpp"
2269#endif
2270# 1734 "/home/runner/work/MFC/MFC/src/simulation/m_data_output.fpp"
2271 deallocate (c_mass)
2272 end if
2273
2274 if (down_sample) then
2275 do i = 1, sys_size
2276 deallocate (q_cons_temp_ds(i)%sf)
2277 end do
2278 deallocate (q_cons_temp_ds)
2279 end if
2280
2281 end subroutine s_finalize_data_output_module
2282
2283end 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_4eq
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 delayfileaccess(processrank)
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
real(wp), dimension(:), allocatable qvs
real(wp), dimension(:), allocatable pi_infs
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 gammas
real(wp), dimension(:), allocatable gs_min
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 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_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_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).