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# 52 "/home/runner/work/MFC/MFC/src/common/include/macros.fpp"
300
301! Allocate and create GPU device memory
302# 72 "/home/runner/work/MFC/MFC/src/common/include/macros.fpp"
303
304! Free GPU device memory and deallocate
305# 80 "/home/runner/work/MFC/MFC/src/common/include/macros.fpp"
306
307! Cray-specific GPU pointer setup for vector fields
308# 104 "/home/runner/work/MFC/MFC/src/common/include/macros.fpp"
309
310! Cray-specific GPU pointer setup for scalar fields
311# 120 "/home/runner/work/MFC/MFC/src/common/include/macros.fpp"
312
313! Cray-specific GPU pointer setup for acoustic source spatials
314# 145 "/home/runner/work/MFC/MFC/src/common/include/macros.fpp"
315
316# 151 "/home/runner/work/MFC/MFC/src/common/include/macros.fpp"
317
318# 158 "/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), dimension(2) :: re !< Cell-avg. Reynolds numbers
511 integer :: j, k, l
512 real(wp) :: icfl_max_loc, icfl_max_glb !< ICFL stability extrema on local and global grids
513 real(wp) :: vcfl_max_loc, vcfl_max_glb !< VCFL stability extrema on local and global grids
514 real(wp) :: ccfl_max_loc, ccfl_max_glb !< CCFL stability extrema on local and global grids
515 real(wp) :: rc_min_loc, rc_min_glb !< Rc stability extrema on local and global grids
516 real(wp) :: icfl, vcfl, ccfl, rc
517 integer :: fl !< Fluid loop iterator
518
519 icfl_max_loc = 0._wp
520 vcfl_max_loc = 0._wp
521 ccfl_max_loc = 0._wp
522 rc_min_loc = huge(1.0_wp)
523 ! Computing Stability Criteria at Current Time-step
524
525# 192 "/home/runner/work/MFC/MFC/src/simulation/m_data_output.fpp"
526
527# 192 "/home/runner/work/MFC/MFC/src/simulation/m_data_output.fpp"
528#if defined(MFC_OpenACC)
529# 192 "/home/runner/work/MFC/MFC/src/simulation/m_data_output.fpp"
530!$acc parallel loop collapse(3) gang vector default(present) private(j, k, l, vel, alpha, Re, rho, vel_sum, pres, gamma, pi_inf, c, qv, icfl, vcfl, Rc, ccfl, fl) &
531# 192 "/home/runner/work/MFC/MFC/src/simulation/m_data_output.fpp"
532!$acc& reduction(max:icfl_max_loc, vcfl_max_loc, ccfl_max_loc) reduction(min:Rc_min_loc)
533# 192 "/home/runner/work/MFC/MFC/src/simulation/m_data_output.fpp"
534#elif defined(MFC_OpenMP)
535# 192 "/home/runner/work/MFC/MFC/src/simulation/m_data_output.fpp"
536
537# 192 "/home/runner/work/MFC/MFC/src/simulation/m_data_output.fpp"
538
539# 192 "/home/runner/work/MFC/MFC/src/simulation/m_data_output.fpp"
540
541# 192 "/home/runner/work/MFC/MFC/src/simulation/m_data_output.fpp"
542!$omp target teams loop defaultmap(firstprivate:scalar) bind(teams,parallel) collapse(3) defaultmap(tofrom:aggregate) defaultmap(tofrom:allocatable) defaultmap(tofrom:pointer) &
543# 192 "/home/runner/work/MFC/MFC/src/simulation/m_data_output.fpp"
544!$omp& private(j, k, l, vel, alpha, Re, rho, vel_sum, pres, gamma, pi_inf, c, qv, icfl, vcfl, Rc, ccfl, fl) reduction(max:icfl_max_loc, vcfl_max_loc, ccfl_max_loc) reduction(min:Rc_min_loc)
545# 192 "/home/runner/work/MFC/MFC/src/simulation/m_data_output.fpp"
546#endif
547# 195 "/home/runner/work/MFC/MFC/src/simulation/m_data_output.fpp"
548 do l = 0, p
549 do k = 0, n
550 do j = 0, m
551 call s_compute_cell_state(q_prim_vf, pres, rho, gamma, pi_inf, re, alpha, vel, vel_sum, qv, j, k, l)
552
553 call s_compute_speed_of_sound(pres, rho, gamma, pi_inf, alpha, c)
554
555 if (any_non_newtonian) then
556 re(1) = 0._wp
557 do fl = 1, num_fluids
558 if (is_non_newtonian(fl)) then
559 re(1) = re(1) + alpha(fl)*hb_mu_max(fl)
560 else
561 re(1) = re(1) + alpha(fl)*fluid_inv_re(fl)
562 end if
563 end do
564 re(1) = 1._wp/max(re(1), sgm_eps)
565 end if
566
567 call s_compute_stability_from_dt(vel, c, rho, re, j, k, l, icfl, vcfl, rc, ccfl)
568
569 icfl_max_loc = max(icfl_max_loc, icfl)
570 vcfl_max_loc = max(vcfl_max_loc, merge(vcfl, 0.0_wp, viscous))
571 ccfl_max_loc = max(ccfl_max_loc, merge(ccfl, 0.0_wp, surface_tension))
572 rc_min_loc = min(rc_min_loc, merge(rc, huge(1.0_wp), viscous))
573 end do
574 end do
575 end do
576
577# 223 "/home/runner/work/MFC/MFC/src/simulation/m_data_output.fpp"
578#if defined(MFC_OpenACC)
579# 223 "/home/runner/work/MFC/MFC/src/simulation/m_data_output.fpp"
580!$acc end parallel loop
581# 223 "/home/runner/work/MFC/MFC/src/simulation/m_data_output.fpp"
582#elif defined(MFC_OpenMP)
583# 223 "/home/runner/work/MFC/MFC/src/simulation/m_data_output.fpp"
584
585# 223 "/home/runner/work/MFC/MFC/src/simulation/m_data_output.fpp"
586!$omp end target teams loop
587# 223 "/home/runner/work/MFC/MFC/src/simulation/m_data_output.fpp"
588#endif
589 ! end: Computing Stability Criteria at Current Time-step
590
591 if (num_procs > 1) then
592 call s_mpi_reduce_stability_criteria_extrema(icfl_max_loc, vcfl_max_loc, rc_min_loc, n_el_bubs_loc, icfl_max_glb, &
593 & vcfl_max_glb, rc_min_glb, n_el_bubs_glb, ccfl_max_loc, ccfl_max_glb)
594 else
595 icfl_max_glb = icfl_max_loc
596 if (viscous) vcfl_max_glb = vcfl_max_loc
597 if (viscous) rc_min_glb = rc_min_loc
598 if (surface_tension) ccfl_max_glb = ccfl_max_loc
599 if (bubbles_lagrange) n_el_bubs_glb = n_el_bubs_loc
600 end if
601
602 if (icfl_max_glb > icfl_max) icfl_max = icfl_max_glb
603
604 if (surface_tension) then
605 if (ccfl_max_glb > ccfl_max) ccfl_max = ccfl_max_glb
606 end if
607
608 if (viscous) then
609 if (vcfl_max_glb > vcfl_max) vcfl_max = vcfl_max_glb
610 if (rc_min_glb < rc_min) rc_min = rc_min_glb
611 end if
612
613 if (proc_rank == 0) then
614 write (3, '(13X,I9,13X,F10.6,13X,F10.6,13X,F10.6)', advance="no") t_step, dt, mytime, icfl_max_glb
615
616 if (surface_tension) then
617 write (3, '(13X,F10.6)', advance="no") ccfl_max_glb
618 end if
619
620 if (viscous) then
621 write (3, '(13X,F10.6,13X,ES16.6)', advance="no") vcfl_max_glb, rc_min_glb
622 end if
623
624 if (bubbles_lagrange) then
625 write (3, '(13X,I10)', advance="no") n_el_bubs_glb
626 end if
627
628 write (3, *) ! new line
629
630 if (.not. f_approx_equal(icfl_max_glb, icfl_max_glb)) then
631 call s_mpi_abort('ICFL is NaN. Exiting.')
632 else if (icfl_max_glb > 1._wp) then
633 print *, 'icfl', icfl_max_glb
634 call s_mpi_abort('ICFL is greater than 1.0. Exiting.')
635 end if
636
637 if (viscous) then
638 if (.not. f_approx_equal(vcfl_max_glb, vcfl_max_glb)) then
639 call s_mpi_abort('VCFL is NaN. Exiting.')
640 else if (vcfl_max_glb > 1._wp) then
641 print *, 'vcfl', vcfl_max_glb
642 call s_mpi_abort('VCFL is greater than 1.0. Exiting.')
643 end if
644 end if
645
646 if (bubbles_lagrange) then
647 if (n_el_bubs_glb == 0) then
648 call s_mpi_abort('No Lagrangian bubbles remain in the domain. Exiting.')
649 end if
650 end if
651 end if
652
653 call s_mpi_barrier()
654
655 end subroutine s_write_run_time_information
656
657 !> Write grid and conservative variable data files in serial format
658 impure subroutine s_write_serial_data_files(q_cons_vf, q_T_sf, q_prim_vf, t_step, bc_type, beta)
659
660 type(scalar_field), dimension(sys_size), intent(inout) :: q_cons_vf
661 type(scalar_field), intent(inout) :: q_t_sf
662 type(scalar_field), dimension(sys_size), intent(inout) :: q_prim_vf
663 integer, intent(in) :: t_step
664 type(scalar_field), intent(inout), optional :: beta
665 type(integer_field), dimension(1:num_dims,-1:1), intent(in) :: bc_type
666 character(LEN=path_len + 2*name_len) :: t_step_dir !< Relative path to the current time-step directory
667 character(LEN=path_len + 3*name_len) :: file_path !< Relative path to the grid and conservative variables data files
668 logical :: file_exist !< Logical used to check existence of current time-step directory
669 character(LEN=15) :: fmt
670 integer :: i, j, k, l, r
671 real(wp) :: gamma, lit_gamma, pi_inf, qv !< Temporary EOS params
672
673 write (t_step_dir, '(A,I0,A,I0)') trim(case_dir) // '/p_all'
674 write (t_step_dir, '(a,i0,a,i0)') trim(case_dir) // '/p_all/p', proc_rank, '/', t_step
675
676 file_path = trim(t_step_dir) // '/.'
677 call my_inquire(file_path, file_exist)
678 if (file_exist) call s_delete_directory(trim(t_step_dir))
679 call s_create_directory(trim(t_step_dir))
680
681 file_path = trim(t_step_dir) // '/x_cb.dat'
682
683 open (2, file=trim(file_path), form='unformatted', status='new')
684 write (2) x_cb(-1:m); close (2)
685
686 if (n > 0) then
687 file_path = trim(t_step_dir) // '/y_cb.dat'
688
689 open (2, file=trim(file_path), form='unformatted', status='new')
690 write (2) y_cb(-1:n); close (2)
691
692 if (p > 0) then
693 file_path = trim(t_step_dir) // '/z_cb.dat'
694
695 open (2, file=trim(file_path), form='unformatted', status='new')
696 write (2) z_cb(-1:p); close (2)
697 end if
698 end if
699
700 do i = 1, sys_size
701 write (file_path, '(A,I0,A)') trim(t_step_dir) // '/q_cons_vf', i, '.dat'
702
703 open (2, file=trim(file_path), form='unformatted', status='new')
704
705 write (2) q_cons_vf(i)%sf(0:m,0:n,0:p); close (2)
706 end do
707
708 ! Lagrangian beta (void fraction) written as q_cons_vf(sys_size+1) to match the parallel I/O path and allow post_process to
709 ! read it.
710 if (bubbles_lagrange) then
711 write (file_path, '(A,I0,A)') trim(t_step_dir) // '/q_cons_vf', sys_size + 1, '.dat'
712
713 open (2, file=trim(file_path), form='unformatted', status='new')
714
715 write (2) beta%sf(0:m,0:n,0:p); close (2)
716 end if
717
718 if (qbmm .and. .not. polytropic) then
719 do i = 1, nb
720 do r = 1, nnode
721 write (file_path, '(A,I0,A)') trim(t_step_dir) // '/pb', sys_size + (i - 1)*nnode + r, '.dat'
722
723 open (2, file=trim(file_path), form='unformatted', status='new')
724
725 write (2) pb_ts(1)%sf(0:m,0:n,0:p,r, i); close (2)
726 end do
727 end do
728
729 do i = 1, nb
730 do r = 1, nnode
731 write (file_path, '(A,I0,A)') trim(t_step_dir) // '/mv', sys_size + (i - 1)*nnode + r, '.dat'
732
733 open (2, file=trim(file_path), form='unformatted', status='new')
734
735 write (2) mv_ts(1)%sf(0:m,0:n,0:p,r, i); close (2)
736 end do
737 end do
738 end if
739
740 ! Writing the IB markers
741 if (ib) then
742 call s_write_serial_ib_data(t_step)
743 end if
744
745 gamma = gammas(1)
746 lit_gamma = isentrope_n(1)
747 pi_inf = pi_infs(1)
748 qv = qvs(1)
749
750 if (precision == precision_single) then
751 fmt = "(2F30.3)"
752 else
753 fmt = "(2F40.14)"
754 end if
755
756 write (t_step_dir, '(A,I0,A,I0)') trim(case_dir) // '/D'
757 file_path = trim(t_step_dir) // '/.'
758
759 inquire (file=trim(file_path), exist=file_exist)
760
761 if (.not. file_exist) call s_create_directory(trim(t_step_dir))
762
763 if ((prim_vars_wrt .or. (n == 0 .and. p == 0)) .and. (.not. igr)) then
765 do i = 1, sys_size
766
767# 401 "/home/runner/work/MFC/MFC/src/simulation/m_data_output.fpp"
768#if defined(MFC_OpenACC)
769# 401 "/home/runner/work/MFC/MFC/src/simulation/m_data_output.fpp"
770!$acc update host(q_prim_vf(i)%sf(:, :, :))
771# 401 "/home/runner/work/MFC/MFC/src/simulation/m_data_output.fpp"
772#elif defined(MFC_OpenMP)
773# 401 "/home/runner/work/MFC/MFC/src/simulation/m_data_output.fpp"
774!$omp target update from(q_prim_vf(i)%sf(:, :, :))
775# 401 "/home/runner/work/MFC/MFC/src/simulation/m_data_output.fpp"
776#endif
777 end do
778 ! q_prim_vf(eqn_idx%bub%beg) stores the value of nb needed in riemann solvers, so replace with true primitive value
779 ! (=1._wp)
780 if (qbmm) then
781 q_prim_vf(eqn_idx%bub%beg)%sf = 1._wp
782 end if
783 end if
784
785 if (n == 0 .and. p == 0) then
786 if (model_eqns == model_eqns_5eq .and. (.not. igr)) then
787 do i = 1, sys_size
788 write (file_path, '(A,I0,A,I2.2,A,I6.6,A)') trim(t_step_dir) // '/prim.', i, '.', proc_rank, '.', t_step, '.dat'
789
790 open (2, file=trim(file_path))
791 do j = 0, m
792 ! todo: revisit change here
793 if (((i >= eqn_idx%adv%beg) .and. (i <= eqn_idx%adv%end))) then
794 write (2, fmt) x_cb(j), q_cons_vf(i)%sf(j, 0, 0)
795 else
796 write (2, fmt) x_cb(j), q_prim_vf(i)%sf(j, 0, 0)
797 end if
798 end do
799 close (2)
800 end do
801 end if
802
803 do i = 1, sys_size
804 write (file_path, '(A,I0,A,I2.2,A,I6.6,A)') trim(t_step_dir) // '/cons.', i, '.', proc_rank, '.', t_step, '.dat'
805
806 open (2, file=trim(file_path))
807 do j = 0, m
808 write (2, fmt) x_cb(j), q_cons_vf(i)%sf(j, 0, 0)
809 end do
810 close (2)
811 end do
812
813 if (qbmm .and. .not. polytropic) then
814 do i = 1, nb
815 do r = 1, nnode
816 write (file_path, '(A,I0,A,I0,A,I2.2,A,I6.6,A)') trim(t_step_dir) // '/pres.', i, '.', r, '.', proc_rank, &
817 & '.', t_step, '.dat'
818
819 open (2, file=trim(file_path))
820 do j = 0, m
821 write (2, fmt) x_cb(j), pb_ts(1)%sf(j, 0, 0, r, i)
822 end do
823 close (2)
824 end do
825 end do
826 do i = 1, nb
827 do r = 1, nnode
828 write (file_path, '(A,I0,A,I0,A,I2.2,A,I6.6,A)') trim(t_step_dir) // '/mv.', i, '.', r, '.', proc_rank, &
829 & '.', t_step, '.dat'
830
831 open (2, file=trim(file_path))
832 do j = 0, m
833 write (2, fmt) x_cb(j), mv_ts(1)%sf(j, 0, 0, r, i)
834 end do
835 close (2)
836 end do
837 end do
838 end if
839 end if
840
841 if (precision == precision_single) then
842 fmt = "(3F30.7)"
843 else
844 fmt = "(3F40.14)"
845 end if
846
847 if ((n > 0) .and. (p == 0)) then
848 do i = 1, sys_size
849 write (file_path, '(A,I0,A,I2.2,A,I6.6,A)') trim(t_step_dir) // '/cons.', i, '.', proc_rank, '.', t_step, '.dat'
850 open (2, file=trim(file_path))
851 do j = 0, m
852 do k = 0, n
853 write (2, fmt) x_cb(j), y_cb(k), q_cons_vf(i)%sf(j, k, 0)
854 end do
855 write (2, *)
856 end do
857 close (2)
858 end do
859
860 if (present(beta)) then
861 write (file_path, '(A,I0,A,I2.2,A,I6.6,A)') trim(t_step_dir) // '/beta.', i, '.', proc_rank, '.', t_step, '.dat'
862 open (2, file=trim(file_path))
863 do j = 0, m
864 do k = 0, n
865 write (2, fmt) x_cb(j), y_cb(k), beta%sf(j, k, 0)
866 end do
867 write (2, *)
868 end do
869 close (2)
870 end if
871
872 if (qbmm .and. .not. polytropic) then
873 do i = 1, nb
874 do r = 1, nnode
875 write (file_path, '(A,I0,A,I0,A,I2.2,A,I6.6,A)') trim(t_step_dir) // '/pres.', i, '.', r, '.', proc_rank, &
876 & '.', t_step, '.dat'
877
878 open (2, file=trim(file_path))
879 do j = 0, m
880 do k = 0, n
881 write (2, fmt) x_cb(j), y_cb(k), pb_ts(1)%sf(j, k, 0, r, i)
882 end do
883 end do
884 close (2)
885 end do
886 end do
887 do i = 1, nb
888 do r = 1, nnode
889 write (file_path, '(A,I0,A,I0,A,I2.2,A,I6.6,A)') trim(t_step_dir) // '/mv.', i, '.', r, '.', proc_rank, &
890 & '.', t_step, '.dat'
891
892 open (2, file=trim(file_path))
893 do j = 0, m
894 do k = 0, n
895 write (2, fmt) x_cb(j), y_cb(k), mv_ts(1)%sf(j, k, 0, r, i)
896 end do
897 end do
898 close (2)
899 end do
900 end do
901 end if
902
903 if (prim_vars_wrt .and. (.not. igr)) then
904 do i = 1, sys_size
905 write (file_path, '(A,I0,A,I2.2,A,I6.6,A)') trim(t_step_dir) // '/prim.', i, '.', proc_rank, '.', t_step, '.dat'
906
907 open (2, file=trim(file_path))
908
909 do j = 0, m
910 do k = 0, n
911 if (((i >= eqn_idx%cont%beg) .and. (i <= eqn_idx%cont%end)) .or. ((i >= eqn_idx%adv%beg) &
912 & .and. (i <= eqn_idx%adv%end))) then
913 write (2, fmt) x_cb(j), y_cb(k), q_cons_vf(i)%sf(j, k, 0)
914 else
915 write (2, fmt) x_cb(j), y_cb(k), q_prim_vf(i)%sf(j, k, 0)
916 end if
917 end do
918 write (2, *)
919 end do
920 close (2)
921 end do
922 end if
923 end if
924
925 if (precision == precision_single) then
926 fmt = "(4F30.7)"
927 else
928 fmt = "(4F40.14)"
929 end if
930
931 if (p > 0) then
932 do i = 1, sys_size
933 write (file_path, '(A,I0,A,I2.2,A,I6.6,A)') trim(t_step_dir) // '/cons.', i, '.', proc_rank, '.', t_step, '.dat'
934 open (2, file=trim(file_path))
935 do j = 0, m
936 do k = 0, n
937 do l = 0, p
938 write (2, fmt) x_cb(j), y_cb(k), z_cb(l), q_cons_vf(i)%sf(j, k, l)
939 end do
940 write (2, *)
941 end do
942 write (2, *)
943 end do
944 close (2)
945 end do
946
947 if (present(beta)) then
948 write (file_path, '(A,I0,A,I2.2,A,I6.6,A)') trim(t_step_dir) // '/beta.', i, '.', proc_rank, '.', t_step, '.dat'
949 open (2, file=trim(file_path))
950 do j = 0, m
951 do k = 0, n
952 do l = 0, p
953 write (2, fmt) x_cb(j), y_cb(k), z_cb(l), beta%sf(j, k, l)
954 end do
955 write (2, *)
956 end do
957 write (2, *)
958 end do
959 close (2)
960 end if
961
962 if (qbmm .and. .not. polytropic) then
963 do i = 1, nb
964 do r = 1, nnode
965 write (file_path, '(A,I0,A,I0,A,I2.2,A,I6.6,A)') trim(t_step_dir) // '/pres.', i, '.', r, '.', proc_rank, &
966 & '.', t_step, '.dat'
967
968 open (2, file=trim(file_path))
969 do j = 0, m
970 do k = 0, n
971 do l = 0, p
972 write (2, fmt) x_cb(j), y_cb(k), z_cb(l), pb_ts(1)%sf(j, k, l, r, i)
973 end do
974 end do
975 end do
976 close (2)
977 end do
978 end do
979 do i = 1, nb
980 do r = 1, nnode
981 write (file_path, '(A,I0,A,I0,A,I2.2,A,I6.6,A)') trim(t_step_dir) // '/mv.', i, '.', r, '.', proc_rank, &
982 & '.', t_step, '.dat'
983
984 open (2, file=trim(file_path))
985 do j = 0, m
986 do k = 0, n
987 do l = 0, p
988 write (2, fmt) x_cb(j), y_cb(k), z_cb(l), mv_ts(1)%sf(j, k, l, r, i)
989 end do
990 end do
991 end do
992 close (2)
993 end do
994 end do
995 end if
996
997 if (prim_vars_wrt .and. (.not. igr)) then
998 do i = 1, sys_size
999 write (file_path, '(A,I0,A,I2.2,A,I6.6,A)') trim(t_step_dir) // '/prim.', i, '.', proc_rank, '.', t_step, '.dat'
1000
1001 open (2, file=trim(file_path))
1002
1003 do j = 0, m
1004 do k = 0, n
1005 do l = 0, p
1006 if (((i >= eqn_idx%cont%beg) .and. (i <= eqn_idx%cont%end)) .or. ((i >= eqn_idx%adv%beg) &
1007 & .and. (i <= eqn_idx%adv%end)) .or. ((i >= eqn_idx%species%beg) &
1008 & .and. (i <= eqn_idx%species%end))) then
1009 write (2, fmt) x_cb(j), y_cb(k), z_cb(l), q_cons_vf(i)%sf(j, k, l)
1010 else
1011 write (2, fmt) x_cb(j), y_cb(k), z_cb(l), q_prim_vf(i)%sf(j, k, l)
1012 end if
1013 end do
1014 write (2, *)
1015 end do
1016 write (2, *)
1017 end do
1018 close (2)
1019 end do
1020 end if
1021 end if
1022
1023 end subroutine s_write_serial_data_files
1024
1025 !> Write grid and conservative variable data files in parallel via MPI I/O
1026 impure subroutine s_write_parallel_data_files(q_cons_vf, t_step, bc_type, beta, q_T_sf)
1027
1028 type(scalar_field), dimension(sys_size), intent(inout) :: q_cons_vf
1029 integer, intent(in) :: t_step
1030 type(scalar_field), intent(inout), optional :: beta
1031 type(integer_field), dimension(1:num_dims,-1:1), intent(in) :: bc_type
1032 type(scalar_field), intent(inout), optional :: q_t_sf
1033
1034#ifdef MFC_MPI
1035 integer :: ifile, ierr, data_size
1036 integer, dimension(MPI_STATUS_SIZE) :: status
1037 integer(kind=MPI_OFFSET_kind) :: disp
1038 integer(kind=MPI_OFFSET_kind) :: m_mok, n_mok, p_mok
1039 integer(kind=MPI_OFFSET_kind) :: wp_mok, var_mok, str_mok
1040 integer(kind=MPI_OFFSET_kind) :: nvars_mok
1041 integer(kind=MPI_OFFSET_kind) :: mok
1042 character(LEN=path_len + 2*name_len) :: file_loc
1043 logical :: file_exist, dir_check
1044 character(len=10) :: t_step_string
1045 integer :: i !< Generic loop iterator
1046 integer :: alt_sys !< Altered system size for the lagrangian subgrid bubble model
1047 ! Down sampling variables
1048 integer :: m_ds, n_ds, p_ds
1049 integer :: m_glb_ds, n_glb_ds, p_glb_ds
1050 integer :: m_glb_save, n_glb_save, p_glb_save !< Global save size
1051
1052 if (down_sample) then
1053 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)
1054 end if
1055
1056 if (present(beta)) then
1057 alt_sys = sys_size + 1
1058 else
1059 alt_sys = sys_size
1060 end if
1061
1062 if (file_per_process) then
1063 call s_int_to_str(t_step, t_step_string)
1064
1065 if (down_sample) then
1066 call s_initialize_mpi_data_ds(m_ds, n_ds, p_ds)
1067 else
1068 if (ib) then
1069 call s_initialize_mpi_data(q_cons_vf, ib_markers=ib_markers, ib_mpi_data=mpi_io_ib_data, qbmm_pb=pb_ts(1), &
1070 & qbmm_mv=mv_ts(1))
1071 else
1072 call s_initialize_mpi_data(q_cons_vf, qbmm_pb=pb_ts(1), qbmm_mv=mv_ts(1))
1073 end if
1074 end if
1075
1076 if (proc_rank == 0) then
1077 file_loc = trim(case_dir) // '/restart_data/lustre_' // trim(t_step_string)
1078 call my_inquire(file_loc, dir_check)
1079 if (dir_check .neqv. .true.) then
1080 call s_create_directory(trim(file_loc))
1081 end if
1082 call s_create_directory(trim(file_loc))
1083 end if
1084 call s_mpi_barrier()
1086
1087 call s_initialize_mpi_data(q_cons_vf, qbmm_pb=pb_ts(1), qbmm_mv=mv_ts(1))
1088
1089 write (file_loc, '(I0,A,i7.7,A)') t_step, '_', proc_rank, '.dat'
1090 file_loc = trim(case_dir) // '/restart_data/lustre_' // trim(t_step_string) // trim(mpiiofs) // trim(file_loc)
1091 inquire (file=trim(file_loc), exist=file_exist)
1092 if (file_exist .and. proc_rank == 0) then
1093 call mpi_file_delete(file_loc, mpi_info_int, ierr)
1094 end if
1095 call mpi_file_open(mpi_comm_self, file_loc, ior(mpi_mode_wronly, mpi_mode_create), mpi_info_int, ifile, ierr)
1096
1097 if (down_sample) then
1098 data_size = (m_ds + 3)*(n_ds + 3)*(p_ds + 3)
1099 m_glb_save = m_glb_ds + 1
1100 n_glb_save = n_glb_ds + 1
1101 p_glb_save = p_glb_ds + 1
1102 else
1103 data_size = (m + 1)*(n + 1)*(p + 1)
1104 m_glb_save = m_glb + 1
1105 n_glb_save = n_glb + 1
1106 p_glb_save = p_glb + 1
1107 end if
1108
1109 m_mok = int(m_glb_save + 1, mpi_offset_kind)
1110 n_mok = int(n_glb_save + 1, mpi_offset_kind)
1111 p_mok = int(p_glb_save + 1, mpi_offset_kind)
1112 wp_mok = int(storage_size(0._stp)/8, mpi_offset_kind)
1113 mok = int(1._wp, mpi_offset_kind)
1114 str_mok = int(name_len, mpi_offset_kind)
1115 nvars_mok = int(sys_size, mpi_offset_kind)
1116
1117 if (bubbles_euler) then
1118 do i = 1, sys_size
1119 var_mok = int(i, mpi_offset_kind)
1120
1121 call mpi_file_write_all(ifile, mpi_io_data%var(i)%sf, data_size*mpi_io_type, mpi_io_p, status, ierr)
1122 end do
1123 if (qbmm .and. .not. polytropic) then
1124 do i = sys_size + 1, sys_size + 2*nb*nnode
1125 var_mok = int(i, mpi_offset_kind)
1126
1127 call mpi_file_write_all(ifile, mpi_io_data%var(i)%sf, data_size*mpi_io_type, mpi_io_p, status, ierr)
1128 end do
1129 end if
1130 else
1131 if (down_sample) then
1132 do i = 1, sys_size ! TODO: check if sys_size is correct
1133 var_mok = int(i, mpi_offset_kind)
1134
1135 call mpi_file_write_all(ifile, q_cons_temp_ds(i)%sf, data_size*mpi_io_type, mpi_io_p, status, ierr)
1136 end do
1137 else
1138 do i = 1, sys_size ! TODO: check if sys_size is correct
1139 var_mok = int(i, mpi_offset_kind)
1140
1141 call mpi_file_write_all(ifile, mpi_io_data%var(i)%sf, data_size*mpi_io_type, mpi_io_p, status, ierr)
1142 end do
1143 end if
1144 end if
1145
1146 call mpi_file_close(ifile, ierr)
1147
1148 if (ib) then
1149 call s_write_parallel_ib_data(t_step)
1150 end if
1151 else
1152 if (ib) then
1153 call s_initialize_mpi_data(q_cons_vf, ib_markers=ib_markers, ib_mpi_data=mpi_io_ib_data, qbmm_pb=pb_ts(1), &
1154 & qbmm_mv=mv_ts(1))
1155 else if (present(beta)) then
1156 call s_initialize_mpi_data(q_cons_vf, beta=beta, qbmm_pb=pb_ts(1), qbmm_mv=mv_ts(1))
1157 else
1158 call s_initialize_mpi_data(q_cons_vf, qbmm_pb=pb_ts(1), qbmm_mv=mv_ts(1))
1159 end if
1160
1161 write (file_loc, '(I0,A)') t_step, '.dat'
1162 file_loc = trim(case_dir) // '/restart_data' // trim(mpiiofs) // trim(file_loc)
1163 inquire (file=trim(file_loc), exist=file_exist)
1164 if (file_exist .and. proc_rank == 0) then
1165 call mpi_file_delete(file_loc, mpi_info_int, ierr)
1166 end if
1167 call mpi_file_open(mpi_comm_world, file_loc, ior(mpi_mode_wronly, mpi_mode_create), mpi_info_int, ifile, ierr)
1168
1169 data_size = (m + 1)*(n + 1)*(p + 1)
1170
1171 m_mok = int(m_glb + 1, mpi_offset_kind)
1172 n_mok = int(n_glb + 1, mpi_offset_kind)
1173 p_mok = int(p_glb + 1, mpi_offset_kind)
1174 wp_mok = int(storage_size(0._stp)/8, mpi_offset_kind)
1175 mok = int(1._wp, mpi_offset_kind)
1176 str_mok = int(name_len, mpi_offset_kind)
1177 nvars_mok = int(alt_sys, mpi_offset_kind)
1178
1179 if (bubbles_euler) then
1180 do i = 1, sys_size
1181 var_mok = int(i, mpi_offset_kind)
1182
1183 disp = m_mok*max(mok, n_mok)*max(mok, p_mok)*wp_mok*(var_mok - 1)
1184
1185 call mpi_file_set_view(ifile, disp, mpi_p, mpi_io_data%view(i), 'native', mpi_info_int, ierr)
1186 call mpi_file_write_all(ifile, mpi_io_data%var(i)%sf, data_size*mpi_io_type, mpi_io_p, status, ierr)
1187 end do
1188 if (qbmm .and. .not. polytropic) then
1189 do i = sys_size + 1, sys_size + 2*nb*nnode
1190 var_mok = int(i, mpi_offset_kind)
1191
1192 disp = m_mok*max(mok, n_mok)*max(mok, p_mok)*wp_mok*(var_mok - 1)
1193
1194 call mpi_file_set_view(ifile, disp, mpi_p, mpi_io_data%view(i), 'native', mpi_info_int, ierr)
1195 call mpi_file_write_all(ifile, mpi_io_data%var(i)%sf, data_size*mpi_io_type, mpi_io_p, status, ierr)
1196 end do
1197 end if
1198 else
1199 do i = 1, sys_size ! TODO: check if sys_size is correct
1200 var_mok = int(i, mpi_offset_kind)
1201
1202 disp = m_mok*max(mok, n_mok)*max(mok, p_mok)*wp_mok*(var_mok - 1)
1203
1204 call mpi_file_set_view(ifile, disp, mpi_p, mpi_io_data%view(i), 'native', mpi_info_int, ierr)
1205 call mpi_file_write_all(ifile, mpi_io_data%var(i)%sf, data_size*mpi_io_type, mpi_io_p, status, ierr)
1206 end do
1207 end if
1208
1209 if (present(beta)) then
1210 var_mok = int(sys_size + 1, mpi_offset_kind)
1211
1212 disp = m_mok*max(mok, n_mok)*max(mok, p_mok)*wp_mok*(var_mok - 1)
1213
1214 call mpi_file_set_view(ifile, disp, mpi_p, mpi_io_data%view(sys_size + 1), 'native', mpi_info_int, ierr)
1215 call mpi_file_write_all(ifile, mpi_io_data%var(sys_size + 1)%sf, data_size*mpi_io_type, mpi_io_p, status, ierr)
1216 end if
1217
1218 call mpi_file_close(ifile, ierr)
1219
1220 if (ib) then
1221 call s_write_parallel_ib_data(t_step)
1222 end if
1223 end if
1224#endif
1225
1226 end subroutine s_write_parallel_data_files
1227
1228 !> Write immersed boundary marker data to a serial (per-processor) unformatted file
1229 subroutine s_write_serial_ib_data(time_step)
1230
1231 integer, intent(in) :: time_step
1232 character(LEN=path_len + 2*name_len) :: file_path
1233 character(LEN=path_len + 2*name_len) :: t_step_dir
1234
1235 write (t_step_dir, '(A,I0,A,I0)') trim(case_dir) // '/p_all'
1236 write (t_step_dir, '(a,i0,a,i0)') trim(case_dir) // '/p_all/p', proc_rank, '/', time_step
1237 write (file_path, '(A,I0,A)') trim(t_step_dir) // '/ib_data.dat'
1238
1239 open (2, file=trim(file_path), form='unformatted', status='new')
1240
1241
1242# 866 "/home/runner/work/MFC/MFC/src/simulation/m_data_output.fpp"
1243#if defined(MFC_OpenACC)
1244# 866 "/home/runner/work/MFC/MFC/src/simulation/m_data_output.fpp"
1245!$acc update host(ib_markers%sf)
1246# 866 "/home/runner/work/MFC/MFC/src/simulation/m_data_output.fpp"
1247#elif defined(MFC_OpenMP)
1248# 866 "/home/runner/work/MFC/MFC/src/simulation/m_data_output.fpp"
1249!$omp target update from(ib_markers%sf)
1250# 866 "/home/runner/work/MFC/MFC/src/simulation/m_data_output.fpp"
1251#endif
1252 write (2) ib_markers%sf(0:m,0:n,0:p); close (2)
1253
1254 end subroutine s_write_serial_ib_data
1255
1256 !> Write immersed boundary marker data in parallel using MPI I/O
1257 subroutine s_write_parallel_ib_data(time_step)
1258
1259 integer, intent(in) :: time_step
1260
1261#ifdef MFC_MPI
1262 character(LEN=path_len + 2*name_len) :: file_loc
1263 integer(kind=MPI_OFFSET_kind) :: disp
1264 integer(kind=MPI_OFFSET_kind) :: m_MOK, n_MOK, p_MOK
1265 integer(kind=MPI_OFFSET_kind) :: WP_MOK, var_MOK, MOK
1266 integer :: ifile, ierr, data_size
1267 integer, dimension(MPI_STATUS_SIZE) :: status
1268
1269
1270# 884 "/home/runner/work/MFC/MFC/src/simulation/m_data_output.fpp"
1271#if defined(MFC_OpenACC)
1272# 884 "/home/runner/work/MFC/MFC/src/simulation/m_data_output.fpp"
1273!$acc update host(ib_markers%sf)
1274# 884 "/home/runner/work/MFC/MFC/src/simulation/m_data_output.fpp"
1275#elif defined(MFC_OpenMP)
1276# 884 "/home/runner/work/MFC/MFC/src/simulation/m_data_output.fpp"
1277!$omp target update from(ib_markers%sf)
1278# 884 "/home/runner/work/MFC/MFC/src/simulation/m_data_output.fpp"
1279#endif
1280
1281 data_size = (m + 1)*(n + 1)*(p + 1)
1282 m_mok = int(m_glb + 1, mpi_offset_kind)
1283 n_mok = int(n_glb + 1, mpi_offset_kind)
1284 p_mok = int(p_glb + 1, mpi_offset_kind)
1285 wp_mok = int(storage_size(0._stp)/8, mpi_offset_kind)
1286 mok = int(1._wp, mpi_offset_kind)
1287
1288 write (file_loc, '(A)') 'ib.dat'
1289 file_loc = trim(case_dir) // '/restart_data' // trim(mpiiofs) // trim(file_loc)
1290
1291 call s_mpi_barrier()
1293
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 real(wp) :: rhoyks(1:num_species)
1571
1572 t = dflt_t_guess
1573
1574 if (time_stepper == 23) then
1575 nondim_time = mytime
1576 else
1577 if (t_step_old /= dflt_int) then
1578 nondim_time = real(t_step + t_step_old, wp)*dt
1579 else
1580 nondim_time = real(t_step, wp)*dt
1581 end if
1582 end if
1583
1584 do i = 1, num_probes
1585 rho = 0._wp
1586 do s = 1, num_vels
1587 vel(s) = 0._wp
1588 end do
1589 pres = 0._wp
1590 gamma = 0._wp
1591 pi_inf = 0._wp
1592 qv = 0._wp
1593 c = 0._wp
1594 accel = 0._wp
1595 nr = 0._wp; r = 0._wp
1596 nrdot = 0._wp; rdot = 0._wp
1597 nbub = 0._wp
1598 m00 = 0._wp
1599 m10 = 0._wp
1600 m01 = 0._wp
1601 m20 = 0._wp
1602 m02 = 0._wp
1603 varr = 0._wp; varv = 0._wp
1604 alf = 0._wp
1605 do s = 1, (num_dims*(num_dims + 1))/2
1606 tau_e(s) = 0._wp
1607 end do
1608 damage_state = 0._wp
1609
1610 if (n == 0) then
1611 if ((probe(i)%x >= x_cb(-1)) .and. (probe(i)%x <= x_cb(m))) then
1612 do s = -1, m
1613 distx(s) = x_cb(s) - probe(i)%x
1614 if (distx(s) < 0._wp) distx(s) = 1000._wp
1615 end do
1616 j = minloc(distx, 1)
1617 if (j == 1) j = 2 ! Pick first point if probe is at edge
1618 k = 0
1619 l = 0
1620
1621 if (chemistry) then
1622 do d = 1, num_species
1623 rhoyks(d) = q_cons_vf(eqn_idx%species%beg + d - 1)%sf(j - 2, k, l)
1624 end do
1625 end if
1626
1627 ! Computing/Sharing necessary state variables
1628 if (hypoelasticity) then
1629 call s_convert_to_mixture_variables(q_cons_vf, j - 2, k, l, rho, gamma, pi_inf, qv, re, g_local, &
1630 & fluid_pp(:)%G)
1631 else
1632 call s_convert_to_mixture_variables(q_cons_vf, j - 2, k, l, rho, gamma, pi_inf, qv)
1633 end if
1634 do s = 1, num_vels
1635 vel(s) = q_cons_vf(eqn_idx%cont%end + s)%sf(j - 2, k, l)/rho
1636 end do
1637 do s = 1, num_fluids
1638 alpha(s) = q_cons_vf(eqn_idx%adv%beg + s - 1)%sf(j - 2, k, l)
1639 end do
1640
1641 dyn_p = 0.5_wp*rho*dot_product(vel, vel)
1642
1643 if (hypoelasticity) 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(eqn_idx%E)%sf(j - 2, k, l), q_cons_vf(eqn_idx%alf)%sf(j - 2, k, l), &
1650 & dyn_p, pi_inf, gamma, rho, qv, rhoyks(:), pres, t, &
1651 & f_hypoelastic_energy(q_cons_vf, j - 2, k, l, rho, g_local))
1652 else
1653 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), &
1654 & dyn_p, pi_inf, gamma, rho, qv, rhoyks, pres, t)
1655 end if
1656
1657 if (hypoelasticity) then
1658 tau_e(1) = q_cons_vf(eqn_idx%stress%end)%sf(j - 2, k, l)/rho
1659 end if
1660
1661 if (bubbles_euler) then
1662 alf = q_cons_vf(eqn_idx%alf)%sf(j - 2, k, l)
1663 if (num_fluids == 3) then
1664 alfgr = q_cons_vf(eqn_idx%alf - 1)%sf(j - 2, k, l)
1665 end if
1666 do s = 1, nb
1667 nr(s) = q_cons_vf(qbmm_idx%rs(s))%sf(j - 2, k, l)
1668 nrdot(s) = q_cons_vf(qbmm_idx%vs(s))%sf(j - 2, k, l)
1669 end do
1670
1671 if (adv_n) then
1672 nbub = q_cons_vf(eqn_idx%n)%sf(j - 2, k, l)
1673 else
1674 nr3 = 0._wp
1675 do s = 1, nb
1676 nr3 = nr3 + weight(s)*(nr(s)**3._wp)
1677 end do
1678
1679 nbub = sqrt((4._wp*pi/3._wp)*nr3/alf)
1680 end if
1681#ifdef MFC_DEBUG
1682 print *, 'In probe, nbub: ', nbub
1683#endif
1684 if (qbmm) then
1685 m00 = q_cons_vf(qbmm_idx%moms(1, 1))%sf(j - 2, k, l)/nbub
1686 m10 = q_cons_vf(qbmm_idx%moms(1, 2))%sf(j - 2, k, l)/nbub
1687 m01 = q_cons_vf(qbmm_idx%moms(1, 3))%sf(j - 2, k, l)/nbub
1688 m20 = q_cons_vf(qbmm_idx%moms(1, 4))%sf(j - 2, k, l)/nbub
1689 m02 = q_cons_vf(qbmm_idx%moms(1, 6))%sf(j - 2, k, l)/nbub
1690
1691 m10 = m10/m00
1692 m01 = m01/m00
1693 m20 = m20/m00
1694 m02 = m02/m00
1695
1696 varr = m20 - m10**2._wp
1697 varv = m02 - m01**2._wp
1698 end if
1699 r(:) = nr(:)/nbub
1700 rdot(:) = nrdot(:)/nbub
1701
1702 ptilde = ptil(j - 2, k, l)
1703 ptot = pres - ptilde
1704 end if
1705
1706 ! Compute mixture sound Speed
1707 call s_compute_speed_of_sound(pres, rho, gamma, pi_inf, alpha, c)
1708 if (hypoelasticity) c = sqrt(c*c + (4._wp/3._wp)*g_local/rho)
1709
1710 accel = accel_mag(j - 2, k, l)
1711 end if
1712 else if (p == 0) then
1713 if (chemistry) then
1714 do d = 1, num_species
1715 rhoyks(d) = q_cons_vf(eqn_idx%species%beg + d - 1)%sf(j - 2, k - 2, l)
1716 end do
1717 end if
1718
1719 if ((probe(i)%x >= x_cb(-1)) .and. (probe(i)%x <= x_cb(m))) then
1720 if ((probe(i)%y >= y_cb(-1)) .and. (probe(i)%y <= y_cb(n))) then
1721 do s = -1, m
1722 distx(s) = x_cb(s) - probe(i)%x
1723 if (distx(s) < 0._wp) distx(s) = 1000._wp
1724 end do
1725 do s = -1, n
1726 disty(s) = y_cb(s) - probe(i)%y
1727 if (disty(s) < 0._wp) disty(s) = 1000._wp
1728 end do
1729 j = minloc(distx, 1)
1730 k = minloc(disty, 1)
1731 if (j == 1) j = 2 ! Pick first point if probe is at edge
1732 if (k == 1) k = 2 ! Pick first point if probe is at edge
1733 l = 0
1734
1735 ! Computing/Sharing necessary state variables
1736 call s_convert_to_mixture_variables(q_cons_vf, j - 2, k - 2, l, rho, gamma, pi_inf, qv, re, g_local, &
1737 & fluid_pp(:)%G)
1738 do s = 1, num_vels
1739 vel(s) = q_cons_vf(eqn_idx%cont%end + s)%sf(j - 2, k - 2, l)/rho
1740 end do
1741 do s = 1, num_fluids
1742 alpha(s) = q_cons_vf(eqn_idx%adv%beg + s - 1)%sf(j - 2, k - 2, l)
1743 end do
1744
1745 dyn_p = 0.5_wp*rho*dot_product(vel, vel)
1746
1747 if (hypoelasticity) 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(eqn_idx%E)%sf(j - 2, k - 2, l), q_cons_vf(eqn_idx%alf)%sf(j - 2, &
1754 & k - 2, l), dyn_p, pi_inf, gamma, rho, qv, rhoyks, pres, t, &
1755 & f_hypoelastic_energy(q_cons_vf, j - 2, k - 2, l, rho, g_local))
1756 else
1757 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, &
1758 & k - 2, l), dyn_p, pi_inf, gamma, rho, qv, rhoyks, pres, t)
1759 end if
1760
1761 if (hypoelasticity) then
1762 do s = 1, 3
1763 tau_e(s) = q_cons_vf(eqn_idx%stress%beg + s - 1)%sf(j - 2, k - 2, l)/rho
1764 end do
1765 end if
1766
1767 if (bubbles_euler) then
1768 alf = q_cons_vf(eqn_idx%alf)%sf(j - 2, k - 2, l)
1769 do s = 1, nb
1770 nr(s) = q_cons_vf(qbmm_idx%rs(s))%sf(j - 2, k - 2, l)
1771 nrdot(s) = q_cons_vf(qbmm_idx%vs(s))%sf(j - 2, k - 2, l)
1772 end do
1773
1774 if (adv_n) then
1775 nbub = q_cons_vf(eqn_idx%n)%sf(j - 2, k - 2, l)
1776 else
1777 nr3 = 0._wp
1778 do s = 1, nb
1779 nr3 = nr3 + weight(s)*(nr(s)**3._wp)
1780 end do
1781
1782 nbub = sqrt((4._wp*pi/3._wp)*nr3/alf)
1783 end if
1784
1785 r(:) = nr(:)/nbub
1786 rdot(:) = nrdot(:)/nbub
1787 end if
1788 ! Compute mixture sound speed
1789 call s_compute_speed_of_sound(pres, rho, gamma, pi_inf, alpha, c)
1790 if (hypoelasticity) c = sqrt(c*c + (4._wp/3._wp)*g_local/rho)
1791 end if
1792 end if
1793 else
1794 if ((probe(i)%x >= x_cb(-1)) .and. (probe(i)%x <= x_cb(m))) then
1795 if ((probe(i)%y >= y_cb(-1)) .and. (probe(i)%y <= y_cb(n))) then
1796 if ((probe(i)%z >= z_cb(-1)) .and. (probe(i)%z <= z_cb(p))) then
1797 do s = -1, m
1798 distx(s) = x_cb(s) - probe(i)%x
1799 if (distx(s) < 0._wp) distx(s) = 1000._wp
1800 end do
1801 do s = -1, n
1802 disty(s) = y_cb(s) - probe(i)%y
1803 if (disty(s) < 0._wp) disty(s) = 1000._wp
1804 end do
1805 do s = -1, p
1806 distz(s) = z_cb(s) - probe(i)%z
1807 if (distz(s) < 0._wp) distz(s) = 1000._wp
1808 end do
1809 j = minloc(distx, 1)
1810 k = minloc(disty, 1)
1811 l = minloc(distz, 1)
1812 if (j == 1) j = 2 ! Pick first point if probe is at edge
1813 if (k == 1) k = 2 ! Pick first point if probe is at edge
1814 if (l == 1) l = 2 ! Pick first point if probe is at edge
1815
1816 ! Computing/Sharing necessary state variables
1817 call s_convert_to_mixture_variables(q_cons_vf, j - 2, k - 2, l - 2, rho, gamma, pi_inf, qv, re, &
1818 & g_local, fluid_pp(:)%G)
1819 do s = 1, num_vels
1820 vel(s) = q_cons_vf(eqn_idx%cont%end + s)%sf(j - 2, k - 2, l - 2)/rho
1821 end do
1822 do s = 1, num_fluids
1823 alpha(s) = q_cons_vf(eqn_idx%adv%beg + s - 1)%sf(j - 2, k - 2, l - 2)
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 (hypoelasticity) 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(eqn_idx%E)%sf(j - 2, k - 2, l - 2), &
1841 & q_cons_vf(eqn_idx%alf)%sf(j - 2, k - 2, l - 2), dyn_p, pi_inf, gamma, &
1842 & rho, qv, rhoyks, pres, t, f_hypoelastic_energy(q_cons_vf, j - 2, k - 2, &
1843 & l - 2, rho, 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 if (hypoelasticity) then
1851 do s = 1, 6
1852 tau_e(s) = q_cons_vf(eqn_idx%stress%beg + s - 1)%sf(j - 2, k - 2, l - 2)/rho
1853 end do
1854 end if
1855
1856 ! Compute mixture sound speed
1857 call s_compute_speed_of_sound(pres, rho, gamma, pi_inf, alpha, c)
1858 if (hypoelasticity) c = sqrt(c*c + (4._wp/3._wp)*g_local/rho)
1859
1860 accel = accel_mag(j - 2, k - 2, l - 2)
1861 end if
1862 end if
1863 end if
1864 end if
1865 if (num_procs > 1) then
1866# 1452 "/home/runner/work/MFC/MFC/src/simulation/m_data_output.fpp"
1867 tmp = rho
1868 call s_mpi_allreduce_sum(tmp, rho)
1869# 1452 "/home/runner/work/MFC/MFC/src/simulation/m_data_output.fpp"
1870 tmp = pres
1871 call s_mpi_allreduce_sum(tmp, pres)
1872# 1452 "/home/runner/work/MFC/MFC/src/simulation/m_data_output.fpp"
1873 tmp = gamma
1874 call s_mpi_allreduce_sum(tmp, gamma)
1875# 1452 "/home/runner/work/MFC/MFC/src/simulation/m_data_output.fpp"
1876 tmp = pi_inf
1877 call s_mpi_allreduce_sum(tmp, pi_inf)
1878# 1452 "/home/runner/work/MFC/MFC/src/simulation/m_data_output.fpp"
1879 tmp = qv
1880 call s_mpi_allreduce_sum(tmp, qv)
1881# 1452 "/home/runner/work/MFC/MFC/src/simulation/m_data_output.fpp"
1882 tmp = c
1883 call s_mpi_allreduce_sum(tmp, c)
1884# 1452 "/home/runner/work/MFC/MFC/src/simulation/m_data_output.fpp"
1885 tmp = accel
1886 call s_mpi_allreduce_sum(tmp, accel)
1887# 1455 "/home/runner/work/MFC/MFC/src/simulation/m_data_output.fpp"
1888
1889 do s = 1, num_vels
1890 tmp = vel(s)
1891 call s_mpi_allreduce_sum(tmp, vel(s))
1892 end do
1893
1894 if (bubbles_euler) then
1895# 1463 "/home/runner/work/MFC/MFC/src/simulation/m_data_output.fpp"
1896 tmp = alf
1897 call s_mpi_allreduce_sum(tmp, alf)
1898# 1463 "/home/runner/work/MFC/MFC/src/simulation/m_data_output.fpp"
1899 tmp = alfgr
1900 call s_mpi_allreduce_sum(tmp, alfgr)
1901# 1463 "/home/runner/work/MFC/MFC/src/simulation/m_data_output.fpp"
1902 tmp = nbub
1903 call s_mpi_allreduce_sum(tmp, nbub)
1904# 1463 "/home/runner/work/MFC/MFC/src/simulation/m_data_output.fpp"
1905 tmp = nr(1)
1906 call s_mpi_allreduce_sum(tmp, nr(1))
1907# 1463 "/home/runner/work/MFC/MFC/src/simulation/m_data_output.fpp"
1908 tmp = nrdot(1)
1909 call s_mpi_allreduce_sum(tmp, nrdot(1))
1910# 1463 "/home/runner/work/MFC/MFC/src/simulation/m_data_output.fpp"
1911 tmp = m00
1912 call s_mpi_allreduce_sum(tmp, m00)
1913# 1463 "/home/runner/work/MFC/MFC/src/simulation/m_data_output.fpp"
1914 tmp = r(1)
1915 call s_mpi_allreduce_sum(tmp, r(1))
1916# 1463 "/home/runner/work/MFC/MFC/src/simulation/m_data_output.fpp"
1917 tmp = rdot(1)
1918 call s_mpi_allreduce_sum(tmp, rdot(1))
1919# 1463 "/home/runner/work/MFC/MFC/src/simulation/m_data_output.fpp"
1920 tmp = ptilde
1921 call s_mpi_allreduce_sum(tmp, ptilde)
1922# 1463 "/home/runner/work/MFC/MFC/src/simulation/m_data_output.fpp"
1923 tmp = ptot
1924 call s_mpi_allreduce_sum(tmp, ptot)
1925# 1466 "/home/runner/work/MFC/MFC/src/simulation/m_data_output.fpp"
1926
1927 if (qbmm) then
1928# 1469 "/home/runner/work/MFC/MFC/src/simulation/m_data_output.fpp"
1929 tmp = varr
1930 call s_mpi_allreduce_sum(tmp, varr)
1931# 1469 "/home/runner/work/MFC/MFC/src/simulation/m_data_output.fpp"
1932 tmp = varv
1933 call s_mpi_allreduce_sum(tmp, varv)
1934# 1469 "/home/runner/work/MFC/MFC/src/simulation/m_data_output.fpp"
1935 tmp = m10
1936 call s_mpi_allreduce_sum(tmp, m10)
1937# 1469 "/home/runner/work/MFC/MFC/src/simulation/m_data_output.fpp"
1938 tmp = m01
1939 call s_mpi_allreduce_sum(tmp, m01)
1940# 1469 "/home/runner/work/MFC/MFC/src/simulation/m_data_output.fpp"
1941 tmp = m20
1942 call s_mpi_allreduce_sum(tmp, m20)
1943# 1469 "/home/runner/work/MFC/MFC/src/simulation/m_data_output.fpp"
1944 tmp = m02
1945 call s_mpi_allreduce_sum(tmp, m02)
1946# 1472 "/home/runner/work/MFC/MFC/src/simulation/m_data_output.fpp"
1947 end if
1948 end if
1949
1950 if (hypoelasticity) then
1951 do s = 1, (num_dims*(num_dims + 1))/2
1952 tmp = tau_e(s)
1953 call s_mpi_allreduce_sum(tmp, tau_e(s))
1954 end do
1955 end if
1956
1957 if (cont_damage) then
1958 tmp = damage_state
1959 call s_mpi_allreduce_sum(tmp, damage_state)
1960 end if
1961 end if
1962 if (proc_rank == 0) then
1963 if (n == 0) then
1964 if (bubbles_euler .and. (num_fluids <= 2)) then
1965 if (qbmm) then
1966 write (i + 30, '(6x,f12.6,14f28.16)') nondim_time, rho, vel(1), pres, alf, r(1), rdot(1), nr(1), &
1967 & nrdot(1), varr, varv, m10, m01, m20, m02
1968 else
1969 write (i + 30, '(6x,f12.6,8f24.8)') nondim_time, rho, vel(1), pres, alf, r(1), rdot(1), nr(1), nrdot(1)
1970 ! ptilde, & ptot
1971 end if
1972 else if (bubbles_euler .and. (num_fluids == 3)) then
1973 write (i + 30, &
1974 & '(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)') &
1975 & nondim_time, rho, vel(1), pres, alf, alfgr, nr(1), nrdot(1), r(1), rdot(1), ptilde, ptot
1976 else if (bubbles_euler .and. num_fluids == 4) then
1977 write (i + 30, &
1978 & '(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)') &
1979 & 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, &
1980 & 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), &
1981 & 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), &
1982 & q_cons_vf(10)%sf(j - 2, 0, 0), nbub, r(1), rdot(1)
1983 else if (hypoelasticity) then
1984 write (i + 30, '(6X,F12.6,F24.8,F24.8,F24.8,F24.8)') nondim_time, rho, vel(1), pres, tau_e(1)
1985 else
1986 write (i + 30, '(6X,F12.6,F24.8,F24.8,F24.8)') nondim_time, rho, vel(1), pres
1987 end if
1988 else if (p == 0) then
1989 if (bubbles_euler) then
1990# 1516 "/home/runner/work/MFC/MFC/src/simulation/m_data_output.fpp"
1991 write (i + 30, '(6X,10F24.8)') nondim_time, rho, vel(1), vel(2), pres, alf, nr(1), nrdot(1), r(1), &
1992 & rdot(1)
1993# 1519 "/home/runner/work/MFC/MFC/src/simulation/m_data_output.fpp"
1994 else if (hypoelasticity) then
1995# 1521 "/home/runner/work/MFC/MFC/src/simulation/m_data_output.fpp"
1996 write (i + 30, '(6X,F12.6,F24.8,F24.8,F24.8,F24.8,' // 'F24.8,F24.8,F24.8)') nondim_time, rho, &
1997 & vel(1), vel(2), pres, tau_e(1), tau_e(2), tau_e(3)
1998# 1524 "/home/runner/work/MFC/MFC/src/simulation/m_data_output.fpp"
1999 else
2000 write (i + 30, '(6X,F12.6,F24.8,F24.8,F24.8)') nondim_time, rho, vel(1), pres
2001 print *, 'time =', nondim_time, 'rho =', rho, 'pres =', pres
2002 end if
2003 else
2004# 1530 "/home/runner/work/MFC/MFC/src/simulation/m_data_output.fpp"
2005 if (hypoelasticity) then
2006 write (i + 30, '(6X,F12.6,16F24.8)') nondim_time, rho, vel(1), vel(2), vel(3), pres, gamma, pi_inf, &
2007 & qv, c, accel, tau_e(1), tau_e(2), tau_e(3), tau_e(4), tau_e(5), tau_e(6)
2008 else
2009 write (i + 30, &
2010 & '(6X,F12.6,F24.8,F24.8,F24.8,F24.8,' // 'F24.8,F24.8,F24.8,F24.8,F24.8,' // 'F24.8)') &
2011 & nondim_time, rho, vel(1), vel(2), vel(3), pres, gamma, pi_inf, qv, c, accel
2012 end if
2013# 1539 "/home/runner/work/MFC/MFC/src/simulation/m_data_output.fpp"
2014 end if
2015 end if
2016 end do
2017
2018 end subroutine s_write_probe_files
2019
2020 !> Write footer with stability criteria extrema and run-time to the information file, then close it
2022
2023 real(wp) :: run_time !< Run-time of the simulation
2024
2025 write (3, '(A)') ' '
2026 write (3, '(A)') ''
2027
2028 write (3, '(A,F9.6)') 'ICFL Max: ', icfl_max
2029 if (surface_tension) write (3, '(A,F9.6)') 'CCFL Max: ', ccfl_max
2030 if (viscous) write (3, '(A,F9.6)') 'VCFL Max: ', vcfl_max
2031 if (viscous) write (3, '(A,ES16.6)') 'Rc Min: ', rc_min
2032
2033 call cpu_time(run_time)
2034
2035 write (3, '(A)') ''
2036 write (3, '(A,I0,A)') 'Run-time: ', int(anint(run_time)), 's'
2037 write (3, '(A)') ' '
2038 close (3)
2039
2041
2042 !> Closes communication files
2043 impure subroutine s_close_com_files()
2044
2045 integer :: i !< Generic loop iterator
2046
2047 do i = 1, num_fluids
2048 close (i + 120)
2049 end do
2050
2051 end subroutine s_close_com_files
2052
2053 !> Closes probe files
2054 impure subroutine s_close_probe_files
2055
2056 integer :: i !< Generic loop iterator
2057
2058 do i = 1, num_probes
2059 close (i + 30)
2060 end do
2061
2062 end subroutine s_close_probe_files
2063
2064 !> Initialize the data output module
2066
2067 integer :: i, m_ds, n_ds, p_ds
2068
2069 if (run_time_info) then
2070 icfl_max = 0._wp
2071 if (surface_tension) then
2072 ccfl_max = 0._wp
2073 end if
2074 if (viscous) then
2075 vcfl_max = 0._wp
2076 rc_min = 1.e12_wp
2077 end if
2078 end if
2079
2080 if (probe_wrt) then
2081#ifdef MFC_DEBUG
2082# 1606 "/home/runner/work/MFC/MFC/src/simulation/m_data_output.fpp"
2083 block
2084# 1606 "/home/runner/work/MFC/MFC/src/simulation/m_data_output.fpp"
2085 use iso_fortran_env, only: output_unit
2086# 1606 "/home/runner/work/MFC/MFC/src/simulation/m_data_output.fpp"
2087
2088# 1606 "/home/runner/work/MFC/MFC/src/simulation/m_data_output.fpp"
2089 print *, 'm_data_output.fpp:1606: ', '@:ALLOCATE(c_mass(num_fluids,5))'
2090# 1606 "/home/runner/work/MFC/MFC/src/simulation/m_data_output.fpp"
2091
2092# 1606 "/home/runner/work/MFC/MFC/src/simulation/m_data_output.fpp"
2093 call flush (output_unit)
2094# 1606 "/home/runner/work/MFC/MFC/src/simulation/m_data_output.fpp"
2095 end block
2096# 1606 "/home/runner/work/MFC/MFC/src/simulation/m_data_output.fpp"
2097#endif
2098# 1606 "/home/runner/work/MFC/MFC/src/simulation/m_data_output.fpp"
2099 allocate (c_mass(num_fluids,5))
2100# 1606 "/home/runner/work/MFC/MFC/src/simulation/m_data_output.fpp"
2101
2102# 1606 "/home/runner/work/MFC/MFC/src/simulation/m_data_output.fpp"
2103
2104# 1606 "/home/runner/work/MFC/MFC/src/simulation/m_data_output.fpp"
2105#if defined(MFC_OpenACC)
2106# 1606 "/home/runner/work/MFC/MFC/src/simulation/m_data_output.fpp"
2107!$acc enter data create(c_mass)
2108# 1606 "/home/runner/work/MFC/MFC/src/simulation/m_data_output.fpp"
2109#elif defined(MFC_OpenMP)
2110# 1606 "/home/runner/work/MFC/MFC/src/simulation/m_data_output.fpp"
2111!$omp target enter data map(always,alloc:c_mass)
2112# 1606 "/home/runner/work/MFC/MFC/src/simulation/m_data_output.fpp"
2113#endif
2114 end if
2115
2116 if (down_sample) then
2117 m_ds = int((m + 1)/3) - 1
2118 n_ds = int((n + 1)/3) - 1
2119 p_ds = int((p + 1)/3) - 1
2120
2121 allocate (q_cons_temp_ds(1:sys_size))
2122 do i = 1, sys_size
2123 allocate (q_cons_temp_ds(i)%sf(-1:m_ds + 1,-1:n_ds + 1,-1:p_ds + 1))
2124 end do
2125 end if
2126
2127 end subroutine s_initialize_data_output_module
2128
2129 !> Module deallocation and/or disassociation procedures
2131
2132 integer :: i
2133
2134 if (probe_wrt) then
2135#ifdef MFC_DEBUG
2136# 1628 "/home/runner/work/MFC/MFC/src/simulation/m_data_output.fpp"
2137 block
2138# 1628 "/home/runner/work/MFC/MFC/src/simulation/m_data_output.fpp"
2139 use iso_fortran_env, only: output_unit
2140# 1628 "/home/runner/work/MFC/MFC/src/simulation/m_data_output.fpp"
2141
2142# 1628 "/home/runner/work/MFC/MFC/src/simulation/m_data_output.fpp"
2143 print *, 'm_data_output.fpp:1628: ', '@:DEALLOCATE(c_mass)'
2144# 1628 "/home/runner/work/MFC/MFC/src/simulation/m_data_output.fpp"
2145
2146# 1628 "/home/runner/work/MFC/MFC/src/simulation/m_data_output.fpp"
2147 call flush (output_unit)
2148# 1628 "/home/runner/work/MFC/MFC/src/simulation/m_data_output.fpp"
2149 end block
2150# 1628 "/home/runner/work/MFC/MFC/src/simulation/m_data_output.fpp"
2151#endif
2152# 1628 "/home/runner/work/MFC/MFC/src/simulation/m_data_output.fpp"
2153
2154# 1628 "/home/runner/work/MFC/MFC/src/simulation/m_data_output.fpp"
2155#if defined(MFC_OpenACC)
2156# 1628 "/home/runner/work/MFC/MFC/src/simulation/m_data_output.fpp"
2157!$acc exit data delete(c_mass)
2158# 1628 "/home/runner/work/MFC/MFC/src/simulation/m_data_output.fpp"
2159#elif defined(MFC_OpenMP)
2160# 1628 "/home/runner/work/MFC/MFC/src/simulation/m_data_output.fpp"
2161!$omp target exit data map(release:c_mass)
2162# 1628 "/home/runner/work/MFC/MFC/src/simulation/m_data_output.fpp"
2163#endif
2164# 1628 "/home/runner/work/MFC/MFC/src/simulation/m_data_output.fpp"
2165 deallocate (c_mass)
2166 end if
2167
2168 if (down_sample) then
2169 do i = 1, sys_size
2170 deallocate (q_cons_temp_ds(i)%sf)
2171 end do
2172 deallocate (q_cons_temp_ds)
2173 end if
2174
2175 end subroutine s_finalize_data_output_module
2176
2177end 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 cell state, CFL calculation, and stability checks.
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.
subroutine, public s_compute_cell_state(q_prim_vf, pres, rho, gamma, pi_inf, re, alpha, vel, vel_sum, qv, j, k, l)
Computes the mixture coefficients, velocity and pressure of one cell.
Conservative-to-primitive variable conversion, mixture property evaluation, and pressure computation.
real(wp) function, public f_hypoelastic_energy(q_cons_vf, j, k, l, rho, g)
Hypoelastic strain energy at one cell, summed over the stress components.
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_compute_pressure(energy, alf, dyn_p, pi_inf, gamma, rho, qv, rhoyks, pres, t, e_e_in, pres_mag)
Compute the pressure from the appropriate equation of state.
subroutine, public s_compute_speed_of_sound(pres, rho, gamma, pi_inf, adv, c)
Speed of sound of a thermodynamic state. Enthalpy is not an argument: for a real state H,...
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).