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