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