MFC
Exascale flow solver
Loading...
Searching...
No Matches
m_thinc.fpp.f90
Go to the documentation of this file.
1# 1 "/home/runner/work/MFC/MFC/src/simulation/m_thinc.fpp"
2!>
3!! @file
4!! @brief Contains module m_thinc
5
6# 1 "/home/runner/work/MFC/MFC/src/common/include/macros.fpp" 1
7# 1 "/home/runner/work/MFC/MFC/src/common/include/parallel_macros.fpp" 1
8# 1 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp" 1
9# 2 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
10# 3 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
11# 4 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
12# 5 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
13# 6 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
14
15# 8 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
16# 9 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
17# 10 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
18
19# 17 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
20
21# 46 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
22
23# 58 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
24
25# 68 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
26
27# 98 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
28
29# 110 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
30
31# 120 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
32
33# 145 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
34! New line at end of file is required for FYPP
35# 2 "/home/runner/work/MFC/MFC/src/common/include/parallel_macros.fpp" 2
36# 1 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp" 1
37# 1 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp" 1
38# 2 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
39# 3 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
40# 4 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
41# 5 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
42# 6 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
43
44# 8 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
45# 9 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
46# 10 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
47
48# 17 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
49
50# 46 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
51
52# 58 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
53
54# 68 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
55
56# 98 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
57
58# 110 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
59
60# 120 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
61
62# 145 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
63! New line at end of file is required for FYPP
64# 2 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp" 2
65
66# 4 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
67# 5 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
68# 6 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
69# 7 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
70# 8 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
71
72# 20 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
73
74# 43 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
75
76# 48 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
77
78# 53 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
79
80# 58 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
81
82# 63 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
83
84# 68 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
85
86# 76 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
87
88# 81 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
89
90# 86 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
91
92# 91 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
93
94# 96 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
95
96# 101 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
97
98# 106 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
99
100# 111 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
101
102# 116 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
103
104# 121 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
105
106# 151 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
107
108# 192 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
109
110# 206 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
111
112# 231 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
113
114# 242 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
115
116# 244 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
117# 255 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
118
119# 284 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
120
121# 294 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
122
123# 304 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
124
125# 313 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
126
127# 330 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
128
129# 340 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
130
131# 347 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
132
133# 353 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
134
135# 359 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
136
137# 365 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
138
139# 371 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
140
141# 377 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
142! New line at end of file is required for FYPP
143# 3 "/home/runner/work/MFC/MFC/src/common/include/parallel_macros.fpp" 2
144# 1 "/home/runner/work/MFC/MFC/src/common/include/acc_macros.fpp" 1
145# 1 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp" 1
146# 2 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
147# 3 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
148# 4 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
149# 5 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
150# 6 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
151
152# 8 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
153# 9 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
154# 10 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
155
156# 17 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
157
158# 46 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
159
160# 58 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
161
162# 68 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
163
164# 98 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
165
166# 110 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
167
168# 120 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
169
170# 145 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
171! New line at end of file is required for FYPP
172# 2 "/home/runner/work/MFC/MFC/src/common/include/acc_macros.fpp" 2
173
174# 7 "/home/runner/work/MFC/MFC/src/common/include/acc_macros.fpp"
175
176# 17 "/home/runner/work/MFC/MFC/src/common/include/acc_macros.fpp"
177
178# 22 "/home/runner/work/MFC/MFC/src/common/include/acc_macros.fpp"
179
180# 27 "/home/runner/work/MFC/MFC/src/common/include/acc_macros.fpp"
181
182# 32 "/home/runner/work/MFC/MFC/src/common/include/acc_macros.fpp"
183
184# 37 "/home/runner/work/MFC/MFC/src/common/include/acc_macros.fpp"
185
186# 42 "/home/runner/work/MFC/MFC/src/common/include/acc_macros.fpp"
187
188# 47 "/home/runner/work/MFC/MFC/src/common/include/acc_macros.fpp"
189
190# 52 "/home/runner/work/MFC/MFC/src/common/include/acc_macros.fpp"
191
192# 57 "/home/runner/work/MFC/MFC/src/common/include/acc_macros.fpp"
193
194# 62 "/home/runner/work/MFC/MFC/src/common/include/acc_macros.fpp"
195
196# 73 "/home/runner/work/MFC/MFC/src/common/include/acc_macros.fpp"
197
198# 78 "/home/runner/work/MFC/MFC/src/common/include/acc_macros.fpp"
199
200# 83 "/home/runner/work/MFC/MFC/src/common/include/acc_macros.fpp"
201
202# 88 "/home/runner/work/MFC/MFC/src/common/include/acc_macros.fpp"
203
204# 103 "/home/runner/work/MFC/MFC/src/common/include/acc_macros.fpp"
205
206# 131 "/home/runner/work/MFC/MFC/src/common/include/acc_macros.fpp"
207
208# 160 "/home/runner/work/MFC/MFC/src/common/include/acc_macros.fpp"
209
210# 175 "/home/runner/work/MFC/MFC/src/common/include/acc_macros.fpp"
211
212# 193 "/home/runner/work/MFC/MFC/src/common/include/acc_macros.fpp"
213
214# 215 "/home/runner/work/MFC/MFC/src/common/include/acc_macros.fpp"
215
216# 244 "/home/runner/work/MFC/MFC/src/common/include/acc_macros.fpp"
217
218# 259 "/home/runner/work/MFC/MFC/src/common/include/acc_macros.fpp"
219
220# 269 "/home/runner/work/MFC/MFC/src/common/include/acc_macros.fpp"
221
222# 278 "/home/runner/work/MFC/MFC/src/common/include/acc_macros.fpp"
223
224# 294 "/home/runner/work/MFC/MFC/src/common/include/acc_macros.fpp"
225
226# 304 "/home/runner/work/MFC/MFC/src/common/include/acc_macros.fpp"
227
228# 311 "/home/runner/work/MFC/MFC/src/common/include/acc_macros.fpp"
229! New line at end of file is required for FYPP
230# 4 "/home/runner/work/MFC/MFC/src/common/include/parallel_macros.fpp" 2
231
232! GPU parallel region (scalar reductions, maxval/minval)
233# 23 "/home/runner/work/MFC/MFC/src/common/include/parallel_macros.fpp"
234
235! GPU parallel loop over threads (most common GPU macro)
236# 43 "/home/runner/work/MFC/MFC/src/common/include/parallel_macros.fpp"
237
238! Required closing for GPU_PARALLEL_LOOP
239# 55 "/home/runner/work/MFC/MFC/src/common/include/parallel_macros.fpp"
240
241! Mark routine for device compilation
242# 112 "/home/runner/work/MFC/MFC/src/common/include/parallel_macros.fpp"
243
244! Declare device-resident data
245# 130 "/home/runner/work/MFC/MFC/src/common/include/parallel_macros.fpp"
246
247! Inner loop within a GPU parallel region
248# 145 "/home/runner/work/MFC/MFC/src/common/include/parallel_macros.fpp"
249
250! Scoped GPU data region
251# 164 "/home/runner/work/MFC/MFC/src/common/include/parallel_macros.fpp"
252
253! Host code with device pointers (for MPI with GPU buffers)
254# 193 "/home/runner/work/MFC/MFC/src/common/include/parallel_macros.fpp"
255
256! Allocate device memory (unscoped)
257# 207 "/home/runner/work/MFC/MFC/src/common/include/parallel_macros.fpp"
258
259! Free device memory
260# 219 "/home/runner/work/MFC/MFC/src/common/include/parallel_macros.fpp"
261
262! Atomic operation on device
263# 231 "/home/runner/work/MFC/MFC/src/common/include/parallel_macros.fpp"
264
265! End atomic capture block
266# 242 "/home/runner/work/MFC/MFC/src/common/include/parallel_macros.fpp"
267
268! Copy data between host and device
269# 254 "/home/runner/work/MFC/MFC/src/common/include/parallel_macros.fpp"
270
271! Synchronization barrier
272# 266 "/home/runner/work/MFC/MFC/src/common/include/parallel_macros.fpp"
273
274! Import GPU library module (openacc or omp_lib)
275# 275 "/home/runner/work/MFC/MFC/src/common/include/parallel_macros.fpp"
276
277! Emit code only for AMD compiler
278# 282 "/home/runner/work/MFC/MFC/src/common/include/parallel_macros.fpp"
279
280! Emit code for non-Cray compilers
281# 289 "/home/runner/work/MFC/MFC/src/common/include/parallel_macros.fpp"
282
283! Emit code only for Cray compiler
284# 296 "/home/runner/work/MFC/MFC/src/common/include/parallel_macros.fpp"
285
286! Emit code for non-NVIDIA compilers
287# 303 "/home/runner/work/MFC/MFC/src/common/include/parallel_macros.fpp"
288
289# 305 "/home/runner/work/MFC/MFC/src/common/include/parallel_macros.fpp"
290# 306 "/home/runner/work/MFC/MFC/src/common/include/parallel_macros.fpp"
291! New line at end of file is required for FYPP
292# 2 "/home/runner/work/MFC/MFC/src/common/include/macros.fpp" 2
293
294# 14 "/home/runner/work/MFC/MFC/src/common/include/macros.fpp"
295
296! Caution: This macro requires the use of a binding script to set CUDA_VISIBLE_DEVICES, such that we have one GPU device per MPI
297! rank. That's because for both cudaMemAdvise (preferred location) and cudaMemPrefetchAsync we use location = device_id = 0. For an
298! example see misc/nvidia_uvm/bind.sh. NVIDIA unified memory page placement hint
299# 57 "/home/runner/work/MFC/MFC/src/common/include/macros.fpp"
300
301! Allocate and create GPU device memory
302# 77 "/home/runner/work/MFC/MFC/src/common/include/macros.fpp"
303
304! Free GPU device memory and deallocate
305# 85 "/home/runner/work/MFC/MFC/src/common/include/macros.fpp"
306
307! Cray-specific GPU pointer setup for vector fields
308# 109 "/home/runner/work/MFC/MFC/src/common/include/macros.fpp"
309
310! Cray-specific GPU pointer setup for scalar fields
311# 125 "/home/runner/work/MFC/MFC/src/common/include/macros.fpp"
312
313! Cray-specific GPU pointer setup for acoustic source spatials
314# 150 "/home/runner/work/MFC/MFC/src/common/include/macros.fpp"
315
316# 156 "/home/runner/work/MFC/MFC/src/common/include/macros.fpp"
317
318# 163 "/home/runner/work/MFC/MFC/src/common/include/macros.fpp"
319! New line at end of file is required for FYPP
320# 6 "/home/runner/work/MFC/MFC/src/simulation/m_thinc.fpp" 2
321
322!> @brief THINC and MTHINC interface compression for volume fraction sharpening. THINC (int_comp=1): 1D directional interface
323!! compression applied after MUSCL/WENO reconstruction. Uses hyperbolic tangent profile per dimension. MTHINC (int_comp=2):
324!! Multi-dimensional THINC that reconstructs a tanh profile oriented along the interface normal (computed from the gradient of
325!! alpha), then integrates that profile over each cell face using Gaussian quadrature. A Newton iteration enforces the conservation
326!! constraint (cell-averaged alpha is preserved). Reference: B. Xie and F. Xiao, "Toward efficient and accurate interface capturing
327!! on arbitrary hybrid unstructured grids: The THINC method with quadratic surface representation and Gaussian quadrature," Journal
328!! of Computational Physics, vol. 349, pp. 415-440, 2017.
330
333 use m_helper
335
336#ifdef MFC_OpenACC
337 use openacc
338#endif
339
340 implicit none
341
343
344 !> 3-point Gauss-Legendre quadrature on [-1/2, 1/2]
345 !> Node locations: +-sqrt(3/5)/2, 0
346 real(wp) :: gq3_pts(3) = [-5e-1_wp*0.7745966692414834_wp, 0._wp, 5e-1_wp*0.7745966692414834_wp]
347 !> Weights: 5/18, 8/18, 5/18
348 real(wp) :: gq3_wts(3) = [5._wp/18._wp, 8._wp/18._wp, 5._wp/18._wp]
349
350# 34 "/home/runner/work/MFC/MFC/src/simulation/m_thinc.fpp"
351#if defined(MFC_OpenACC)
352# 34 "/home/runner/work/MFC/MFC/src/simulation/m_thinc.fpp"
353!$acc declare copyin(gq3_pts, gq3_wts)
354# 34 "/home/runner/work/MFC/MFC/src/simulation/m_thinc.fpp"
355#elif defined(MFC_OpenMP)
356# 34 "/home/runner/work/MFC/MFC/src/simulation/m_thinc.fpp"
357!$omp declare target to(gq3_pts, gq3_wts)
358# 34 "/home/runner/work/MFC/MFC/src/simulation/m_thinc.fpp"
359#endif
360
361 real(wp), allocatable, dimension(:,:,:,:) :: mthinc_nhat !> Unit normal vector
362 real(wp), allocatable, dimension(:,:,:) :: mthinc_d !> Interface position parameter
363
364# 38 "/home/runner/work/MFC/MFC/src/simulation/m_thinc.fpp"
365#if defined(MFC_OpenACC)
366# 38 "/home/runner/work/MFC/MFC/src/simulation/m_thinc.fpp"
367!$acc declare create(mthinc_nhat, mthinc_d)
368# 38 "/home/runner/work/MFC/MFC/src/simulation/m_thinc.fpp"
369#elif defined(MFC_OpenMP)
370# 38 "/home/runner/work/MFC/MFC/src/simulation/m_thinc.fpp"
371!$omp declare target (mthinc_nhat, mthinc_d)
372# 38 "/home/runner/work/MFC/MFC/src/simulation/m_thinc.fpp"
373#endif
374
375contains
376
377 !> @brief Stable difference: log_cosh(a+h) - log_cosh(a-h) = 2*atanh(tanh(a)*tanh(h)). Avoids catastrophic cancellation when h
378 !! is small relative to a.
379 function f_log_cosh_diff(a, h) result(res)
380
381
382# 46 "/home/runner/work/MFC/MFC/src/simulation/m_thinc.fpp"
383#if MFC_OpenACC
384# 46 "/home/runner/work/MFC/MFC/src/simulation/m_thinc.fpp"
385!$acc routine seq
386# 46 "/home/runner/work/MFC/MFC/src/simulation/m_thinc.fpp"
387#elif MFC_OpenMP
388# 46 "/home/runner/work/MFC/MFC/src/simulation/m_thinc.fpp"
389
390# 46 "/home/runner/work/MFC/MFC/src/simulation/m_thinc.fpp"
391
392# 46 "/home/runner/work/MFC/MFC/src/simulation/m_thinc.fpp"
393!$omp declare target device_type(any)
394# 46 "/home/runner/work/MFC/MFC/src/simulation/m_thinc.fpp"
395#endif
396 real(wp), intent(in) :: a, h
397 real(wp) :: res, t
398
399 t = tanh(a)*tanh(h)
400 res = 2._wp*atanh(sign(min(abs(t), 1._wp - epsilon(1._wp)), t))
401
402 end function f_log_cosh_diff
403
404 !> @brief Analytical 1-D integral of the THINC function
405 function f_thinc_integral_1d(a, b) result(res)
406
407
408# 58 "/home/runner/work/MFC/MFC/src/simulation/m_thinc.fpp"
409#if MFC_OpenACC
410# 58 "/home/runner/work/MFC/MFC/src/simulation/m_thinc.fpp"
411!$acc routine seq
412# 58 "/home/runner/work/MFC/MFC/src/simulation/m_thinc.fpp"
413#elif MFC_OpenMP
414# 58 "/home/runner/work/MFC/MFC/src/simulation/m_thinc.fpp"
415
416# 58 "/home/runner/work/MFC/MFC/src/simulation/m_thinc.fpp"
417
418# 58 "/home/runner/work/MFC/MFC/src/simulation/m_thinc.fpp"
419!$omp declare target device_type(any)
420# 58 "/home/runner/work/MFC/MFC/src/simulation/m_thinc.fpp"
421#endif
422 real(wp), intent(in) :: a, b
423 real(wp) :: res
424
425 if (abs(b) < verysmall) then
426 res = 5e-1_wp*(1._wp + tanh(a))
427 else
428 res = 5e-1_wp + f_log_cosh_diff(a, 5e-1_wp*b)/(2._wp*b)
429 end if
430
431 end function f_thinc_integral_1d
432
433 !> @brief Volume integral of H(xi) = 0.5*(1 + tanh(beta*(n.xi + d))) over the cell [-1/2, 1/2]^ndim
434 function f_mthinc_volume_integral(n1, n2, n3, d, beta, ndim) result(res)
435
436
437# 73 "/home/runner/work/MFC/MFC/src/simulation/m_thinc.fpp"
438#if MFC_OpenACC
439# 73 "/home/runner/work/MFC/MFC/src/simulation/m_thinc.fpp"
440!$acc routine seq
441# 73 "/home/runner/work/MFC/MFC/src/simulation/m_thinc.fpp"
442#elif MFC_OpenMP
443# 73 "/home/runner/work/MFC/MFC/src/simulation/m_thinc.fpp"
444
445# 73 "/home/runner/work/MFC/MFC/src/simulation/m_thinc.fpp"
446
447# 73 "/home/runner/work/MFC/MFC/src/simulation/m_thinc.fpp"
448!$omp declare target device_type(any)
449# 73 "/home/runner/work/MFC/MFC/src/simulation/m_thinc.fpp"
450#endif
451 real(wp), intent(in) :: n1, n2, n3, d, beta
452 integer, intent(in) :: ndim
453 real(wp) :: res, a
454 integer :: q1, q2
455
456 if (ndim == 1) then
457 res = f_thinc_integral_1d(beta*d, beta*n1)
458 else if (ndim == 2) then
459 res = 0._wp
460 do q1 = 1, 3
461 a = beta*(n2*gq3_pts(q1) + d)
462 res = res + gq3_wts(q1)*f_thinc_integral_1d(a, beta*n1)
463 end do
464 else
465 res = 0._wp
466 do q1 = 1, 3
467 do q2 = 1, 3
468 a = beta*(n2*gq3_pts(q1) + n3*gq3_pts(q2) + d)
469 res = res + gq3_wts(q1)*gq3_wts(q2)*f_thinc_integral_1d(a, beta*n1)
470 end do
471 end do
472 end if
473
474 end function f_mthinc_volume_integral
475
476 !> @brief Derivative dV/dd of the volume integral (for Newton iteration)
477 function f_mthinc_volume_integral_dd(n1, n2, n3, d, beta, ndim) result(res)
478
479
480# 102 "/home/runner/work/MFC/MFC/src/simulation/m_thinc.fpp"
481#if MFC_OpenACC
482# 102 "/home/runner/work/MFC/MFC/src/simulation/m_thinc.fpp"
483!$acc routine seq
484# 102 "/home/runner/work/MFC/MFC/src/simulation/m_thinc.fpp"
485#elif MFC_OpenMP
486# 102 "/home/runner/work/MFC/MFC/src/simulation/m_thinc.fpp"
487
488# 102 "/home/runner/work/MFC/MFC/src/simulation/m_thinc.fpp"
489
490# 102 "/home/runner/work/MFC/MFC/src/simulation/m_thinc.fpp"
491!$omp declare target device_type(any)
492# 102 "/home/runner/work/MFC/MFC/src/simulation/m_thinc.fpp"
493#endif
494 real(wp), intent(in) :: n1, n2, n3, d, beta
495 integer, intent(in) :: ndim
496 real(wp) :: res, th
497 integer :: q1, q2, q3
498
499 res = 0._wp
500 if (ndim == 1) then
501 do q1 = 1, 3
502 th = tanh(beta*(n1*gq3_pts(q1) + d))
503 res = res + gq3_wts(q1)*(1._wp - th*th)
504 end do
505 else if (ndim == 2) then
506 do q1 = 1, 3
507 do q2 = 1, 3
508 th = tanh(beta*(n1*gq3_pts(q1) + n2*gq3_pts(q2) + d))
509 res = res + gq3_wts(q1)*gq3_wts(q2)*(1._wp - th*th)
510 end do
511 end do
512 else
513 do q1 = 1, 3
514 do q2 = 1, 3
515 do q3 = 1, 3
516 th = tanh(beta*(n1*gq3_pts(q1) + n2*gq3_pts(q2) + n3*gq3_pts(q3) + d))
517 res = res + gq3_wts(q1)*gq3_wts(q2)*gq3_wts(q3)*(1._wp - th*th)
518 end do
519 end do
520 end do
521 end if
522 res = 5e-1_wp*beta*res
523
524 end function f_mthinc_volume_integral_dd
525
526 !> @brief Solve for the interface-position parameter d
527 function f_mthinc_solve_d(n1, n2, n3, beta, alpha_cell, ndim) result(d)
528
529
530# 138 "/home/runner/work/MFC/MFC/src/simulation/m_thinc.fpp"
531#if MFC_OpenACC
532# 138 "/home/runner/work/MFC/MFC/src/simulation/m_thinc.fpp"
533!$acc routine seq
534# 138 "/home/runner/work/MFC/MFC/src/simulation/m_thinc.fpp"
535#elif MFC_OpenMP
536# 138 "/home/runner/work/MFC/MFC/src/simulation/m_thinc.fpp"
537
538# 138 "/home/runner/work/MFC/MFC/src/simulation/m_thinc.fpp"
539
540# 138 "/home/runner/work/MFC/MFC/src/simulation/m_thinc.fpp"
541!$omp declare target device_type(any)
542# 138 "/home/runner/work/MFC/MFC/src/simulation/m_thinc.fpp"
543#endif
544 real(wp), intent(in) :: n1, n2, n3, beta, alpha_cell
545 integer, intent(in) :: ndim
546 real(wp) :: d, v, residual, dv
547 integer :: iter
548
549 d = 0._wp
550 do iter = 1, 30
551 v = f_mthinc_volume_integral(n1, n2, n3, d, beta, ndim)
552 residual = v - alpha_cell
553 if (abs(residual) < verysmall) exit
554 dv = f_mthinc_volume_integral_dd(n1, n2, n3, d, beta, ndim)
555 if (abs(dv) < verysmall) exit
556 d = d - residual/dv
557 end do
558
559 end function f_mthinc_solve_d
560
561 !> @brief Face-averaged THINC function at a cell face. face_dir: 1=x, 2=y, 3=z. face_pos: -0.5 (low) or +0.5 (high). The face
562 !! coordinate in face_dir is fixed; remaining directions are integrated over [-1/2, 1/2] (analytically along one, Gauss for
563 !! others).
564 function f_mthinc_face_average(n1, n2, n3, d, beta, face_dir, face_pos, ndim) result(res)
565
566
567# 161 "/home/runner/work/MFC/MFC/src/simulation/m_thinc.fpp"
568#if MFC_OpenACC
569# 161 "/home/runner/work/MFC/MFC/src/simulation/m_thinc.fpp"
570!$acc routine seq
571# 161 "/home/runner/work/MFC/MFC/src/simulation/m_thinc.fpp"
572#elif MFC_OpenMP
573# 161 "/home/runner/work/MFC/MFC/src/simulation/m_thinc.fpp"
574
575# 161 "/home/runner/work/MFC/MFC/src/simulation/m_thinc.fpp"
576
577# 161 "/home/runner/work/MFC/MFC/src/simulation/m_thinc.fpp"
578!$omp declare target device_type(any)
579# 161 "/home/runner/work/MFC/MFC/src/simulation/m_thinc.fpp"
580#endif
581 real(wp), intent(in) :: n1, n2, n3, d, beta, face_pos
582 integer, intent(in) :: face_dir, ndim
583 real(wp) :: res, a, n_face, n_t0, n_t1
584 integer :: q
585
586 if (ndim == 1) then
587 res = 5e-1_wp*(1._wp + tanh(beta*(n1*face_pos + d)))
588 else if (ndim == 2) then
589 ! One transverse direction
590 if (face_dir == 1) then
591 n_face = n1; n_t0 = n2
592 else
593 n_face = n2; n_t0 = n1
594 end if
595 a = beta*(n_face*face_pos + d)
596 res = f_thinc_integral_1d(a, beta*n_t0)
597 else
598 ! Two transverse directions
599 if (face_dir == 1) then
600 n_face = n1; n_t0 = n2; n_t1 = n3
601 else if (face_dir == 2) then
602 n_face = n2; n_t0 = n1; n_t1 = n3
603 else
604 n_face = n3; n_t0 = n1; n_t1 = n2
605 end if
606 res = 0._wp
607 do q = 1, 3
608 a = beta*(n_face*face_pos + n_t0*gq3_pts(q) + d)
609 res = res + gq3_wts(q)*f_thinc_integral_1d(a, beta*n_t1)
610 end do
611 end if
612
613 end function f_mthinc_face_average
614
616
617 if (int_comp == int_comp_mthinc) then
618#ifdef MFC_DEBUG
619# 199 "/home/runner/work/MFC/MFC/src/simulation/m_thinc.fpp"
620 block
621# 199 "/home/runner/work/MFC/MFC/src/simulation/m_thinc.fpp"
622 use iso_fortran_env, only: output_unit
623# 199 "/home/runner/work/MFC/MFC/src/simulation/m_thinc.fpp"
624
625# 199 "/home/runner/work/MFC/MFC/src/simulation/m_thinc.fpp"
626 print *, 'm_thinc.fpp:199: ', '@:ALLOCATE(mthinc_nhat(1:3, idwbuff(1)%beg:idwbuff(1)%end, idwbuff(2)%beg:idwbuff(2)%end, idwbuff(3)%beg:idwbuff(3)%end))'
627# 199 "/home/runner/work/MFC/MFC/src/simulation/m_thinc.fpp"
628
629# 199 "/home/runner/work/MFC/MFC/src/simulation/m_thinc.fpp"
630 call flush (output_unit)
631# 199 "/home/runner/work/MFC/MFC/src/simulation/m_thinc.fpp"
632 end block
633# 199 "/home/runner/work/MFC/MFC/src/simulation/m_thinc.fpp"
634#endif
635# 199 "/home/runner/work/MFC/MFC/src/simulation/m_thinc.fpp"
636 allocate (mthinc_nhat(1:3, idwbuff(1)%beg:idwbuff(1)%end, idwbuff(2)%beg:idwbuff(2)%end, idwbuff(3)%beg:idwbuff(3)%end))
637# 199 "/home/runner/work/MFC/MFC/src/simulation/m_thinc.fpp"
638
639# 199 "/home/runner/work/MFC/MFC/src/simulation/m_thinc.fpp"
640
641# 199 "/home/runner/work/MFC/MFC/src/simulation/m_thinc.fpp"
642#if defined(MFC_OpenACC)
643# 199 "/home/runner/work/MFC/MFC/src/simulation/m_thinc.fpp"
644!$acc enter data create(mthinc_nhat)
645# 199 "/home/runner/work/MFC/MFC/src/simulation/m_thinc.fpp"
646#elif defined(MFC_OpenMP)
647# 199 "/home/runner/work/MFC/MFC/src/simulation/m_thinc.fpp"
648!$omp target enter data map(always,alloc:mthinc_nhat)
649# 199 "/home/runner/work/MFC/MFC/src/simulation/m_thinc.fpp"
650#endif
651# 201 "/home/runner/work/MFC/MFC/src/simulation/m_thinc.fpp"
652#ifdef MFC_DEBUG
653# 201 "/home/runner/work/MFC/MFC/src/simulation/m_thinc.fpp"
654 block
655# 201 "/home/runner/work/MFC/MFC/src/simulation/m_thinc.fpp"
656 use iso_fortran_env, only: output_unit
657# 201 "/home/runner/work/MFC/MFC/src/simulation/m_thinc.fpp"
658
659# 201 "/home/runner/work/MFC/MFC/src/simulation/m_thinc.fpp"
660 print *, 'm_thinc.fpp:201: ', '@:ALLOCATE(mthinc_d(idwbuff(1)%beg:idwbuff(1)%end, idwbuff(2)%beg:idwbuff(2)%end, idwbuff(3)%beg:idwbuff(3)%end))'
661# 201 "/home/runner/work/MFC/MFC/src/simulation/m_thinc.fpp"
662
663# 201 "/home/runner/work/MFC/MFC/src/simulation/m_thinc.fpp"
664 call flush (output_unit)
665# 201 "/home/runner/work/MFC/MFC/src/simulation/m_thinc.fpp"
666 end block
667# 201 "/home/runner/work/MFC/MFC/src/simulation/m_thinc.fpp"
668#endif
669# 201 "/home/runner/work/MFC/MFC/src/simulation/m_thinc.fpp"
670 allocate (mthinc_d(idwbuff(1)%beg:idwbuff(1)%end, idwbuff(2)%beg:idwbuff(2)%end, idwbuff(3)%beg:idwbuff(3)%end))
671# 201 "/home/runner/work/MFC/MFC/src/simulation/m_thinc.fpp"
672
673# 201 "/home/runner/work/MFC/MFC/src/simulation/m_thinc.fpp"
674
675# 201 "/home/runner/work/MFC/MFC/src/simulation/m_thinc.fpp"
676#if defined(MFC_OpenACC)
677# 201 "/home/runner/work/MFC/MFC/src/simulation/m_thinc.fpp"
678!$acc enter data create(mthinc_d)
679# 201 "/home/runner/work/MFC/MFC/src/simulation/m_thinc.fpp"
680#elif defined(MFC_OpenMP)
681# 201 "/home/runner/work/MFC/MFC/src/simulation/m_thinc.fpp"
682!$omp target enter data map(always,alloc:mthinc_d)
683# 201 "/home/runner/work/MFC/MFC/src/simulation/m_thinc.fpp"
684#endif
685 end if
686
687 end subroutine s_initialize_thinc_module
688
689 !> @brief Compute the unit normal and interface-position parameter d at each interface cell
691
692 type(scalar_field), dimension(:), intent(in) :: v_vf
693 integer :: j, k, l
694 real(wp) :: nr_x, nr_y, nr_z, nmag, nmax, ac
695 type(int_bounds_info), dimension(3) :: id_norm
696
697
698# 214 "/home/runner/work/MFC/MFC/src/simulation/m_thinc.fpp"
699
700# 214 "/home/runner/work/MFC/MFC/src/simulation/m_thinc.fpp"
701#if defined(MFC_OpenACC)
702# 214 "/home/runner/work/MFC/MFC/src/simulation/m_thinc.fpp"
703!$acc parallel loop collapse(3) gang vector default(present) private(j, k, l)
704# 214 "/home/runner/work/MFC/MFC/src/simulation/m_thinc.fpp"
705#elif defined(MFC_OpenMP)
706# 214 "/home/runner/work/MFC/MFC/src/simulation/m_thinc.fpp"
707
708# 214 "/home/runner/work/MFC/MFC/src/simulation/m_thinc.fpp"
709
710# 214 "/home/runner/work/MFC/MFC/src/simulation/m_thinc.fpp"
711
712# 214 "/home/runner/work/MFC/MFC/src/simulation/m_thinc.fpp"
713!$omp target teams loop defaultmap(firstprivate:scalar) bind(teams,parallel) collapse(3) defaultmap(tofrom:aggregate) defaultmap(tofrom:allocatable) defaultmap(tofrom:pointer) private(j, k, l)
714# 214 "/home/runner/work/MFC/MFC/src/simulation/m_thinc.fpp"
715#endif
716 do l = idwbuff(3)%beg, idwbuff(3)%end
717 do k = idwbuff(2)%beg, idwbuff(2)%end
718 do j = idwbuff(1)%beg, idwbuff(1)%end
719 mthinc_nhat(1, j, k, l) = 0._wp
720 mthinc_nhat(2, j, k, l) = 0._wp
721 mthinc_nhat(3, j, k, l) = 0._wp
722 mthinc_d(j, k, l) = 0._wp
723 end do
724 end do
725 end do
726
727# 225 "/home/runner/work/MFC/MFC/src/simulation/m_thinc.fpp"
728#if defined(MFC_OpenACC)
729# 225 "/home/runner/work/MFC/MFC/src/simulation/m_thinc.fpp"
730!$acc end parallel loop
731# 225 "/home/runner/work/MFC/MFC/src/simulation/m_thinc.fpp"
732#elif defined(MFC_OpenMP)
733# 225 "/home/runner/work/MFC/MFC/src/simulation/m_thinc.fpp"
734
735# 225 "/home/runner/work/MFC/MFC/src/simulation/m_thinc.fpp"
736!$omp end target teams loop
737# 225 "/home/runner/work/MFC/MFC/src/simulation/m_thinc.fpp"
738#endif
739
740 id_norm(1)%beg = idwbuff(1)%beg + 1; id_norm(1)%end = idwbuff(1)%end - 1
741 id_norm(2)%beg = 0; id_norm(2)%end = 0
742 id_norm(3)%beg = 0; id_norm(3)%end = 0
743 if (n > 0) then
744 id_norm(2)%beg = idwbuff(2)%beg + 1; id_norm(2)%end = idwbuff(2)%end - 1
745 end if
746 if (p > 0) then
747 id_norm(3)%beg = idwbuff(3)%beg + 1; id_norm(3)%end = idwbuff(3)%end - 1
748 end if
749
750 ! Compute unit normal and solve for d at interior
751
752# 238 "/home/runner/work/MFC/MFC/src/simulation/m_thinc.fpp"
753
754# 238 "/home/runner/work/MFC/MFC/src/simulation/m_thinc.fpp"
755#if defined(MFC_OpenACC)
756# 238 "/home/runner/work/MFC/MFC/src/simulation/m_thinc.fpp"
757!$acc parallel loop collapse(3) gang vector default(present) private(j, k, l, nr_x, nr_y, nr_z, nmag, nmax, ac) copyin(id_norm)
758# 238 "/home/runner/work/MFC/MFC/src/simulation/m_thinc.fpp"
759#elif defined(MFC_OpenMP)
760# 238 "/home/runner/work/MFC/MFC/src/simulation/m_thinc.fpp"
761
762# 238 "/home/runner/work/MFC/MFC/src/simulation/m_thinc.fpp"
763
764# 238 "/home/runner/work/MFC/MFC/src/simulation/m_thinc.fpp"
765
766# 238 "/home/runner/work/MFC/MFC/src/simulation/m_thinc.fpp"
767!$omp target teams loop defaultmap(firstprivate:scalar) bind(teams,parallel) collapse(3) defaultmap(tofrom:aggregate) defaultmap(tofrom:allocatable) defaultmap(tofrom:pointer) &
768# 238 "/home/runner/work/MFC/MFC/src/simulation/m_thinc.fpp"
769!$omp& private(j, k, l, nr_x, nr_y, nr_z, nmag, nmax, ac) map(to:id_norm)
770# 238 "/home/runner/work/MFC/MFC/src/simulation/m_thinc.fpp"
771#endif
772 do l = id_norm(3)%beg, id_norm(3)%end
773 do k = id_norm(2)%beg, id_norm(2)%end
774 do j = id_norm(1)%beg, id_norm(1)%end
775 ac = v_vf(eqn_idx%adv%beg)%sf(j, k, l)
776
777 if (ac >= ic_eps .and. ac <= 1._wp - ic_eps) then
778 nr_x = (v_vf(eqn_idx%adv%beg)%sf(j + 1, k, l) - v_vf(eqn_idx%adv%beg)%sf(j - 1, k, &
779 & l))*(x_cb(j) - x_cb(j - 1))/(x_cc(j + 1) - x_cc(j - 1))
780
781 nr_y = 0._wp
782 if (n > 0) then
783 nr_y = (v_vf(eqn_idx%adv%beg)%sf(j, k + 1, l) - v_vf(eqn_idx%adv%beg)%sf(j, k - 1, &
784 & l))*(y_cb(k) - y_cb(k - 1))/(y_cc(k + 1) - y_cc(k - 1))
785 end if
786
787 nr_z = 0._wp
788 if (p > 0) then
789 nr_z = (v_vf(eqn_idx%adv%beg)%sf(j, k, l + 1) - v_vf(eqn_idx%adv%beg)%sf(j, k, &
790 & l - 1))*(z_cb(l) - z_cb(l - 1))/(z_cc(l + 1) - z_cc(l - 1))
791 end if
792
793 nmag = sqrt(nr_x*nr_x + nr_y*nr_y + nr_z*nr_z)
794
795 if (nmag > verysmall) then
796 nr_x = nr_x/nmag
797 nr_y = nr_y/nmag
798 nr_z = nr_z/nmag
799
800 ! Snap near-grid-aligned normals to exact alignment
801 nmax = max(abs(nr_x), abs(nr_y), abs(nr_z))
802 if (abs(nr_x) < mthinc_align_tol*nmax) nr_x = 0._wp
803 if (abs(nr_y) < mthinc_align_tol*nmax) nr_y = 0._wp
804 if (abs(nr_z) < mthinc_align_tol*nmax) nr_z = 0._wp
805 nmag = sqrt(nr_x*nr_x + nr_y*nr_y + nr_z*nr_z)
806 if (nmag > verysmall) then
807 nr_x = nr_x/nmag
808 nr_y = nr_y/nmag
809 nr_z = nr_z/nmag
810
811 mthinc_nhat(1, j, k, l) = nr_x
812 mthinc_nhat(2, j, k, l) = nr_y
813 mthinc_nhat(3, j, k, l) = nr_z
814
815 mthinc_d(j, k, l) = f_mthinc_solve_d(nr_x, nr_y, nr_z, ic_beta, ac, num_dims)
816 end if
817 end if
818 end if
819 end do
820 end do
821 end do
822
823# 289 "/home/runner/work/MFC/MFC/src/simulation/m_thinc.fpp"
824#if defined(MFC_OpenACC)
825# 289 "/home/runner/work/MFC/MFC/src/simulation/m_thinc.fpp"
826!$acc end parallel loop
827# 289 "/home/runner/work/MFC/MFC/src/simulation/m_thinc.fpp"
828#elif defined(MFC_OpenMP)
829# 289 "/home/runner/work/MFC/MFC/src/simulation/m_thinc.fpp"
830
831# 289 "/home/runner/work/MFC/MFC/src/simulation/m_thinc.fpp"
832!$omp end target teams loop
833# 289 "/home/runner/work/MFC/MFC/src/simulation/m_thinc.fpp"
834#endif
835
836 end subroutine s_compute_mthinc_normals
837
838 !> @brief Applies THINC (int_comp=1) or MTHINC (int_comp=2) interface compression to sharpen volume-fraction and density
839 !! reconstructions at material interfaces. Called after WENO/MUSCL reconstruction per direction.
840 subroutine s_thinc_compression(v_rs_ws, vL_rs_vf_x, vR_rs_vf_x, recon_dir, is1_d, is2_d, is3_d)
841
842 real(wp), dimension(idwbuff(1)%beg:,idwbuff(2)%beg:,idwbuff(3)%beg:,1:), intent(in) :: v_rs_ws
843 real(wp), dimension(idwbuff(1)%beg:,idwbuff(2)%beg:,idwbuff(3)%beg:,1:), intent(inout) :: vl_rs_vf_x, vr_rs_vf_x
844 integer, intent(in) :: recon_dir
845 type(int_bounds_info), intent(in) :: is1_d, is2_d, is3_d
846 integer :: j, k, l
847 real(wp) :: acl, acr, ac, athinc, qmin, qmax, a, b, c
848 real(wp) :: sgn, moncon, beta_eff
849 real(wp) :: nh1, nh2, nh3, d_local, rho1, rho2
850 real(wp) :: rho_b, rho_e
851
852# 311 "/home/runner/work/MFC/MFC/src/simulation/m_thinc.fpp"
853# 312 "/home/runner/work/MFC/MFC/src/simulation/m_thinc.fpp"
854# 313 "/home/runner/work/MFC/MFC/src/simulation/m_thinc.fpp"
855 if (recon_dir == 1) then
856
857# 314 "/home/runner/work/MFC/MFC/src/simulation/m_thinc.fpp"
858
859# 314 "/home/runner/work/MFC/MFC/src/simulation/m_thinc.fpp"
860#if defined(MFC_OpenACC)
861# 314 "/home/runner/work/MFC/MFC/src/simulation/m_thinc.fpp"
862!$acc parallel loop collapse(3) gang vector default(present) private(j, k, l, aCL, aC, aCR, aTHINC, moncon, sgn, qmin, qmax, A, B, C, beta_eff, nh1, nh2, nh3, d_local, rho1, rho2, rho_b, rho_e)
863# 314 "/home/runner/work/MFC/MFC/src/simulation/m_thinc.fpp"
864#elif defined(MFC_OpenMP)
865# 314 "/home/runner/work/MFC/MFC/src/simulation/m_thinc.fpp"
866
867# 314 "/home/runner/work/MFC/MFC/src/simulation/m_thinc.fpp"
868
869# 314 "/home/runner/work/MFC/MFC/src/simulation/m_thinc.fpp"
870
871# 314 "/home/runner/work/MFC/MFC/src/simulation/m_thinc.fpp"
872!$omp target teams loop defaultmap(firstprivate:scalar) bind(teams,parallel) collapse(3) defaultmap(tofrom:aggregate) defaultmap(tofrom:allocatable) defaultmap(tofrom:pointer) &
873# 314 "/home/runner/work/MFC/MFC/src/simulation/m_thinc.fpp"
874!$omp& private(j, k, l, aCL, aC, aCR, aTHINC, moncon, sgn, qmin, qmax, A, B, C, beta_eff, nh1, nh2, nh3, d_local, rho1, rho2, rho_b, rho_e)
875# 314 "/home/runner/work/MFC/MFC/src/simulation/m_thinc.fpp"
876#endif
877# 316 "/home/runner/work/MFC/MFC/src/simulation/m_thinc.fpp"
878 do l = is3_d%beg, is3_d%end
879 do k = is2_d%beg, is2_d%end
880 do j = is1_d%beg, is1_d%end
881 acl = v_rs_ws(j - 1, k, l, eqn_idx%adv%beg)
882 ac = v_rs_ws(j, k, l, eqn_idx%adv%beg)
883 acr = v_rs_ws(j + 1, k, l, eqn_idx%adv%beg)
884
885 if (ac >= ic_eps .and. ac <= 1._wp - ic_eps) then
886 if (int_comp == int_comp_mthinc .and. n > 0) then ! MTHINC
887 ! Map reshaped (j,k,l) to physical (ix,iy,iz)
888
889 nh1 = mthinc_nhat(1, j, k, l)
890 nh2 = mthinc_nhat(2, j, k, l)
891 nh3 = mthinc_nhat(3, j, k, l)
892 d_local = mthinc_d(j, k, l)
893
894 ! Skip if no valid normal was computed
895 if (nh1*nh1 + nh2*nh2 + nh3*nh3 > 5e-1_wp) then
896 rho1 = v_rs_ws(j, k, l, eqn_idx%cont%beg)/ac
897 rho2 = v_rs_ws(j, k, l, eqn_idx%cont%end)/(1._wp - ac)
898
899 ! Left face
900 athinc = f_mthinc_face_average(nh1, nh2, nh3, d_local, ic_beta, 1, -5e-1_wp, &
901 & num_dims)
902 if (athinc < ic_eps) athinc = ic_eps
903 if (athinc > 1._wp - ic_eps) athinc = 1._wp - ic_eps
904 vl_rs_vf_x(j, k, l, eqn_idx%cont%beg) = rho1*athinc
905 vl_rs_vf_x(j, k, l, eqn_idx%cont%end) = rho2*(1._wp - athinc)
906 vl_rs_vf_x(j, k, l, eqn_idx%adv%beg) = athinc
907 vl_rs_vf_x(j, k, l, eqn_idx%adv%end) = 1._wp - athinc
908
909 ! Right face
910 athinc = f_mthinc_face_average(nh1, nh2, nh3, d_local, ic_beta, 1, 5e-1_wp, &
911 & num_dims)
912 if (athinc < ic_eps) athinc = ic_eps
913 if (athinc > 1._wp - ic_eps) athinc = 1._wp - ic_eps
914 vr_rs_vf_x(j, k, l, eqn_idx%cont%beg) = rho1*athinc
915 vr_rs_vf_x(j, k, l, eqn_idx%cont%end) = rho2*(1._wp - athinc)
916 vr_rs_vf_x(j, k, l, eqn_idx%adv%beg) = athinc
917 vr_rs_vf_x(j, k, l, eqn_idx%adv%end) = 1._wp - athinc
918 end if
919 else ! THINC
920 moncon = (acr - ac)*(ac - acl)
921
922 if (moncon > moncon_cutoff) then
923 if (acr - acl > 0._wp) then
924 sgn = 1._wp
925 else
926 sgn = -1._wp
927 end if
928
929 beta_eff = ic_beta
930
931 qmin = min(acr, acl)
932 qmax = max(acr, acl) - qmin
933
934 c = (ac - qmin + sgm_eps)/(qmax + sgm_eps)
935 b = exp(sgn*beta_eff*(2._wp*c - 1._wp))
936 a = (b/cosh(beta_eff) - 1._wp)/tanh(beta_eff)
937
938 rho_b = v_rs_ws(j, k, l, eqn_idx%cont%beg)/ac
939 rho_e = v_rs_ws(j, k, l, eqn_idx%cont%end)/(1._wp - ac)
940
941 ! Left face
942 athinc = qmin + 5e-1_wp*qmax*(1._wp + sgn*a)
943 if (athinc < ic_eps) athinc = ic_eps
944 if (athinc > 1._wp - ic_eps) athinc = 1._wp - ic_eps
945 vl_rs_vf_x(j, k, l, eqn_idx%cont%beg) = rho_b*athinc
946 vl_rs_vf_x(j, k, l, eqn_idx%cont%end) = rho_e*(1._wp - athinc)
947 vl_rs_vf_x(j, k, l, eqn_idx%adv%beg) = athinc
948 vl_rs_vf_x(j, k, l, eqn_idx%adv%end) = 1._wp - athinc
949
950 ! Right face
951 athinc = qmin + 5e-1_wp*qmax*(1._wp + sgn*(tanh(beta_eff) + a)/(1._wp + a*tanh(beta_eff)))
952 if (athinc < ic_eps) athinc = ic_eps
953 if (athinc > 1._wp - ic_eps) athinc = 1._wp - ic_eps
954 vr_rs_vf_x(j, k, l, eqn_idx%cont%beg) = rho_b*athinc
955 vr_rs_vf_x(j, k, l, eqn_idx%cont%end) = rho_e*(1._wp - athinc)
956 vr_rs_vf_x(j, k, l, eqn_idx%adv%beg) = athinc
957 vr_rs_vf_x(j, k, l, eqn_idx%adv%end) = 1._wp - athinc
958 end if
959 end if
960 end if
961 end do
962 end do
963 end do
964
965# 402 "/home/runner/work/MFC/MFC/src/simulation/m_thinc.fpp"
966#if defined(MFC_OpenACC)
967# 402 "/home/runner/work/MFC/MFC/src/simulation/m_thinc.fpp"
968!$acc end parallel loop
969# 402 "/home/runner/work/MFC/MFC/src/simulation/m_thinc.fpp"
970#elif defined(MFC_OpenMP)
971# 402 "/home/runner/work/MFC/MFC/src/simulation/m_thinc.fpp"
972
973# 402 "/home/runner/work/MFC/MFC/src/simulation/m_thinc.fpp"
974!$omp end target teams loop
975# 402 "/home/runner/work/MFC/MFC/src/simulation/m_thinc.fpp"
976#endif
977 end if
978# 311 "/home/runner/work/MFC/MFC/src/simulation/m_thinc.fpp"
979# 312 "/home/runner/work/MFC/MFC/src/simulation/m_thinc.fpp"
980# 313 "/home/runner/work/MFC/MFC/src/simulation/m_thinc.fpp"
981 if (recon_dir == 2) then
982
983# 314 "/home/runner/work/MFC/MFC/src/simulation/m_thinc.fpp"
984
985# 314 "/home/runner/work/MFC/MFC/src/simulation/m_thinc.fpp"
986#if defined(MFC_OpenACC)
987# 314 "/home/runner/work/MFC/MFC/src/simulation/m_thinc.fpp"
988!$acc parallel loop collapse(3) gang vector default(present) private(j, k, l, aCL, aC, aCR, aTHINC, moncon, sgn, qmin, qmax, A, B, C, beta_eff, nh1, nh2, nh3, d_local, rho1, rho2, rho_b, rho_e)
989# 314 "/home/runner/work/MFC/MFC/src/simulation/m_thinc.fpp"
990#elif defined(MFC_OpenMP)
991# 314 "/home/runner/work/MFC/MFC/src/simulation/m_thinc.fpp"
992
993# 314 "/home/runner/work/MFC/MFC/src/simulation/m_thinc.fpp"
994
995# 314 "/home/runner/work/MFC/MFC/src/simulation/m_thinc.fpp"
996
997# 314 "/home/runner/work/MFC/MFC/src/simulation/m_thinc.fpp"
998!$omp target teams loop defaultmap(firstprivate:scalar) bind(teams,parallel) collapse(3) defaultmap(tofrom:aggregate) defaultmap(tofrom:allocatable) defaultmap(tofrom:pointer) &
999# 314 "/home/runner/work/MFC/MFC/src/simulation/m_thinc.fpp"
1000!$omp& private(j, k, l, aCL, aC, aCR, aTHINC, moncon, sgn, qmin, qmax, A, B, C, beta_eff, nh1, nh2, nh3, d_local, rho1, rho2, rho_b, rho_e)
1001# 314 "/home/runner/work/MFC/MFC/src/simulation/m_thinc.fpp"
1002#endif
1003# 316 "/home/runner/work/MFC/MFC/src/simulation/m_thinc.fpp"
1004 do l = is3_d%beg, is3_d%end
1005 do k = is1_d%beg, is1_d%end
1006 do j = is2_d%beg, is2_d%end
1007 acl = v_rs_ws(j, k - 1, l, eqn_idx%adv%beg)
1008 ac = v_rs_ws(j, k, l, eqn_idx%adv%beg)
1009 acr = v_rs_ws(j, k + 1, l, eqn_idx%adv%beg)
1010
1011 if (ac >= ic_eps .and. ac <= 1._wp - ic_eps) then
1012 if (int_comp == int_comp_mthinc .and. n > 0) then ! MTHINC
1013 ! Map reshaped (j,k,l) to physical (ix,iy,iz)
1014
1015 nh1 = mthinc_nhat(1, j, k, l)
1016 nh2 = mthinc_nhat(2, j, k, l)
1017 nh3 = mthinc_nhat(3, j, k, l)
1018 d_local = mthinc_d(j, k, l)
1019
1020 ! Skip if no valid normal was computed
1021 if (nh1*nh1 + nh2*nh2 + nh3*nh3 > 5e-1_wp) then
1022 rho1 = v_rs_ws(j, k, l, eqn_idx%cont%beg)/ac
1023 rho2 = v_rs_ws(j, k, l, eqn_idx%cont%end)/(1._wp - ac)
1024
1025 ! Left face
1026 athinc = f_mthinc_face_average(nh1, nh2, nh3, d_local, ic_beta, 2, -5e-1_wp, &
1027 & num_dims)
1028 if (athinc < ic_eps) athinc = ic_eps
1029 if (athinc > 1._wp - ic_eps) athinc = 1._wp - ic_eps
1030 vl_rs_vf_x(j, k, l, eqn_idx%cont%beg) = rho1*athinc
1031 vl_rs_vf_x(j, k, l, eqn_idx%cont%end) = rho2*(1._wp - athinc)
1032 vl_rs_vf_x(j, k, l, eqn_idx%adv%beg) = athinc
1033 vl_rs_vf_x(j, k, l, eqn_idx%adv%end) = 1._wp - athinc
1034
1035 ! Right face
1036 athinc = f_mthinc_face_average(nh1, nh2, nh3, d_local, ic_beta, 2, 5e-1_wp, &
1037 & num_dims)
1038 if (athinc < ic_eps) athinc = ic_eps
1039 if (athinc > 1._wp - ic_eps) athinc = 1._wp - ic_eps
1040 vr_rs_vf_x(j, k, l, eqn_idx%cont%beg) = rho1*athinc
1041 vr_rs_vf_x(j, k, l, eqn_idx%cont%end) = rho2*(1._wp - athinc)
1042 vr_rs_vf_x(j, k, l, eqn_idx%adv%beg) = athinc
1043 vr_rs_vf_x(j, k, l, eqn_idx%adv%end) = 1._wp - athinc
1044 end if
1045 else ! THINC
1046 moncon = (acr - ac)*(ac - acl)
1047
1048 if (moncon > moncon_cutoff) then
1049 if (acr - acl > 0._wp) then
1050 sgn = 1._wp
1051 else
1052 sgn = -1._wp
1053 end if
1054
1055 beta_eff = ic_beta
1056
1057 qmin = min(acr, acl)
1058 qmax = max(acr, acl) - qmin
1059
1060 c = (ac - qmin + sgm_eps)/(qmax + sgm_eps)
1061 b = exp(sgn*beta_eff*(2._wp*c - 1._wp))
1062 a = (b/cosh(beta_eff) - 1._wp)/tanh(beta_eff)
1063
1064 rho_b = v_rs_ws(j, k, l, eqn_idx%cont%beg)/ac
1065 rho_e = v_rs_ws(j, k, l, eqn_idx%cont%end)/(1._wp - ac)
1066
1067 ! Left face
1068 athinc = qmin + 5e-1_wp*qmax*(1._wp + sgn*a)
1069 if (athinc < ic_eps) athinc = ic_eps
1070 if (athinc > 1._wp - ic_eps) athinc = 1._wp - ic_eps
1071 vl_rs_vf_x(j, k, l, eqn_idx%cont%beg) = rho_b*athinc
1072 vl_rs_vf_x(j, k, l, eqn_idx%cont%end) = rho_e*(1._wp - athinc)
1073 vl_rs_vf_x(j, k, l, eqn_idx%adv%beg) = athinc
1074 vl_rs_vf_x(j, k, l, eqn_idx%adv%end) = 1._wp - athinc
1075
1076 ! Right face
1077 athinc = qmin + 5e-1_wp*qmax*(1._wp + sgn*(tanh(beta_eff) + a)/(1._wp + a*tanh(beta_eff)))
1078 if (athinc < ic_eps) athinc = ic_eps
1079 if (athinc > 1._wp - ic_eps) athinc = 1._wp - ic_eps
1080 vr_rs_vf_x(j, k, l, eqn_idx%cont%beg) = rho_b*athinc
1081 vr_rs_vf_x(j, k, l, eqn_idx%cont%end) = rho_e*(1._wp - athinc)
1082 vr_rs_vf_x(j, k, l, eqn_idx%adv%beg) = athinc
1083 vr_rs_vf_x(j, k, l, eqn_idx%adv%end) = 1._wp - athinc
1084 end if
1085 end if
1086 end if
1087 end do
1088 end do
1089 end do
1090
1091# 402 "/home/runner/work/MFC/MFC/src/simulation/m_thinc.fpp"
1092#if defined(MFC_OpenACC)
1093# 402 "/home/runner/work/MFC/MFC/src/simulation/m_thinc.fpp"
1094!$acc end parallel loop
1095# 402 "/home/runner/work/MFC/MFC/src/simulation/m_thinc.fpp"
1096#elif defined(MFC_OpenMP)
1097# 402 "/home/runner/work/MFC/MFC/src/simulation/m_thinc.fpp"
1098
1099# 402 "/home/runner/work/MFC/MFC/src/simulation/m_thinc.fpp"
1100!$omp end target teams loop
1101# 402 "/home/runner/work/MFC/MFC/src/simulation/m_thinc.fpp"
1102#endif
1103 end if
1104# 311 "/home/runner/work/MFC/MFC/src/simulation/m_thinc.fpp"
1105# 312 "/home/runner/work/MFC/MFC/src/simulation/m_thinc.fpp"
1106# 313 "/home/runner/work/MFC/MFC/src/simulation/m_thinc.fpp"
1107 if (recon_dir == 3) then
1108
1109# 314 "/home/runner/work/MFC/MFC/src/simulation/m_thinc.fpp"
1110
1111# 314 "/home/runner/work/MFC/MFC/src/simulation/m_thinc.fpp"
1112#if defined(MFC_OpenACC)
1113# 314 "/home/runner/work/MFC/MFC/src/simulation/m_thinc.fpp"
1114!$acc parallel loop collapse(3) gang vector default(present) private(j, k, l, aCL, aC, aCR, aTHINC, moncon, sgn, qmin, qmax, A, B, C, beta_eff, nh1, nh2, nh3, d_local, rho1, rho2, rho_b, rho_e)
1115# 314 "/home/runner/work/MFC/MFC/src/simulation/m_thinc.fpp"
1116#elif defined(MFC_OpenMP)
1117# 314 "/home/runner/work/MFC/MFC/src/simulation/m_thinc.fpp"
1118
1119# 314 "/home/runner/work/MFC/MFC/src/simulation/m_thinc.fpp"
1120
1121# 314 "/home/runner/work/MFC/MFC/src/simulation/m_thinc.fpp"
1122
1123# 314 "/home/runner/work/MFC/MFC/src/simulation/m_thinc.fpp"
1124!$omp target teams loop defaultmap(firstprivate:scalar) bind(teams,parallel) collapse(3) defaultmap(tofrom:aggregate) defaultmap(tofrom:allocatable) defaultmap(tofrom:pointer) &
1125# 314 "/home/runner/work/MFC/MFC/src/simulation/m_thinc.fpp"
1126!$omp& private(j, k, l, aCL, aC, aCR, aTHINC, moncon, sgn, qmin, qmax, A, B, C, beta_eff, nh1, nh2, nh3, d_local, rho1, rho2, rho_b, rho_e)
1127# 314 "/home/runner/work/MFC/MFC/src/simulation/m_thinc.fpp"
1128#endif
1129# 316 "/home/runner/work/MFC/MFC/src/simulation/m_thinc.fpp"
1130 do l = is1_d%beg, is1_d%end
1131 do k = is2_d%beg, is2_d%end
1132 do j = is3_d%beg, is3_d%end
1133 acl = v_rs_ws(j, k, l - 1, eqn_idx%adv%beg)
1134 ac = v_rs_ws(j, k, l, eqn_idx%adv%beg)
1135 acr = v_rs_ws(j, k, l + 1, eqn_idx%adv%beg)
1136
1137 if (ac >= ic_eps .and. ac <= 1._wp - ic_eps) then
1138 if (int_comp == int_comp_mthinc .and. n > 0) then ! MTHINC
1139 ! Map reshaped (j,k,l) to physical (ix,iy,iz)
1140
1141 nh1 = mthinc_nhat(1, j, k, l)
1142 nh2 = mthinc_nhat(2, j, k, l)
1143 nh3 = mthinc_nhat(3, j, k, l)
1144 d_local = mthinc_d(j, k, l)
1145
1146 ! Skip if no valid normal was computed
1147 if (nh1*nh1 + nh2*nh2 + nh3*nh3 > 5e-1_wp) then
1148 rho1 = v_rs_ws(j, k, l, eqn_idx%cont%beg)/ac
1149 rho2 = v_rs_ws(j, k, l, eqn_idx%cont%end)/(1._wp - ac)
1150
1151 ! Left face
1152 athinc = f_mthinc_face_average(nh1, nh2, nh3, d_local, ic_beta, 3, -5e-1_wp, &
1153 & num_dims)
1154 if (athinc < ic_eps) athinc = ic_eps
1155 if (athinc > 1._wp - ic_eps) athinc = 1._wp - ic_eps
1156 vl_rs_vf_x(j, k, l, eqn_idx%cont%beg) = rho1*athinc
1157 vl_rs_vf_x(j, k, l, eqn_idx%cont%end) = rho2*(1._wp - athinc)
1158 vl_rs_vf_x(j, k, l, eqn_idx%adv%beg) = athinc
1159 vl_rs_vf_x(j, k, l, eqn_idx%adv%end) = 1._wp - athinc
1160
1161 ! Right face
1162 athinc = f_mthinc_face_average(nh1, nh2, nh3, d_local, ic_beta, 3, 5e-1_wp, &
1163 & num_dims)
1164 if (athinc < ic_eps) athinc = ic_eps
1165 if (athinc > 1._wp - ic_eps) athinc = 1._wp - ic_eps
1166 vr_rs_vf_x(j, k, l, eqn_idx%cont%beg) = rho1*athinc
1167 vr_rs_vf_x(j, k, l, eqn_idx%cont%end) = rho2*(1._wp - athinc)
1168 vr_rs_vf_x(j, k, l, eqn_idx%adv%beg) = athinc
1169 vr_rs_vf_x(j, k, l, eqn_idx%adv%end) = 1._wp - athinc
1170 end if
1171 else ! THINC
1172 moncon = (acr - ac)*(ac - acl)
1173
1174 if (moncon > moncon_cutoff) then
1175 if (acr - acl > 0._wp) then
1176 sgn = 1._wp
1177 else
1178 sgn = -1._wp
1179 end if
1180
1181 beta_eff = ic_beta
1182
1183 qmin = min(acr, acl)
1184 qmax = max(acr, acl) - qmin
1185
1186 c = (ac - qmin + sgm_eps)/(qmax + sgm_eps)
1187 b = exp(sgn*beta_eff*(2._wp*c - 1._wp))
1188 a = (b/cosh(beta_eff) - 1._wp)/tanh(beta_eff)
1189
1190 rho_b = v_rs_ws(j, k, l, eqn_idx%cont%beg)/ac
1191 rho_e = v_rs_ws(j, k, l, eqn_idx%cont%end)/(1._wp - ac)
1192
1193 ! Left face
1194 athinc = qmin + 5e-1_wp*qmax*(1._wp + sgn*a)
1195 if (athinc < ic_eps) athinc = ic_eps
1196 if (athinc > 1._wp - ic_eps) athinc = 1._wp - ic_eps
1197 vl_rs_vf_x(j, k, l, eqn_idx%cont%beg) = rho_b*athinc
1198 vl_rs_vf_x(j, k, l, eqn_idx%cont%end) = rho_e*(1._wp - athinc)
1199 vl_rs_vf_x(j, k, l, eqn_idx%adv%beg) = athinc
1200 vl_rs_vf_x(j, k, l, eqn_idx%adv%end) = 1._wp - athinc
1201
1202 ! Right face
1203 athinc = qmin + 5e-1_wp*qmax*(1._wp + sgn*(tanh(beta_eff) + a)/(1._wp + a*tanh(beta_eff)))
1204 if (athinc < ic_eps) athinc = ic_eps
1205 if (athinc > 1._wp - ic_eps) athinc = 1._wp - ic_eps
1206 vr_rs_vf_x(j, k, l, eqn_idx%cont%beg) = rho_b*athinc
1207 vr_rs_vf_x(j, k, l, eqn_idx%cont%end) = rho_e*(1._wp - athinc)
1208 vr_rs_vf_x(j, k, l, eqn_idx%adv%beg) = athinc
1209 vr_rs_vf_x(j, k, l, eqn_idx%adv%end) = 1._wp - athinc
1210 end if
1211 end if
1212 end if
1213 end do
1214 end do
1215 end do
1216
1217# 402 "/home/runner/work/MFC/MFC/src/simulation/m_thinc.fpp"
1218#if defined(MFC_OpenACC)
1219# 402 "/home/runner/work/MFC/MFC/src/simulation/m_thinc.fpp"
1220!$acc end parallel loop
1221# 402 "/home/runner/work/MFC/MFC/src/simulation/m_thinc.fpp"
1222#elif defined(MFC_OpenMP)
1223# 402 "/home/runner/work/MFC/MFC/src/simulation/m_thinc.fpp"
1224
1225# 402 "/home/runner/work/MFC/MFC/src/simulation/m_thinc.fpp"
1226!$omp end target teams loop
1227# 402 "/home/runner/work/MFC/MFC/src/simulation/m_thinc.fpp"
1228#endif
1229 end if
1230# 405 "/home/runner/work/MFC/MFC/src/simulation/m_thinc.fpp"
1231
1232 end subroutine s_thinc_compression
1233
1235
1236 if (int_comp == int_comp_mthinc) then
1237#ifdef MFC_DEBUG
1238# 411 "/home/runner/work/MFC/MFC/src/simulation/m_thinc.fpp"
1239 block
1240# 411 "/home/runner/work/MFC/MFC/src/simulation/m_thinc.fpp"
1241 use iso_fortran_env, only: output_unit
1242# 411 "/home/runner/work/MFC/MFC/src/simulation/m_thinc.fpp"
1243
1244# 411 "/home/runner/work/MFC/MFC/src/simulation/m_thinc.fpp"
1245 print *, 'm_thinc.fpp:411: ', '@:DEALLOCATE(mthinc_nhat)'
1246# 411 "/home/runner/work/MFC/MFC/src/simulation/m_thinc.fpp"
1247
1248# 411 "/home/runner/work/MFC/MFC/src/simulation/m_thinc.fpp"
1249 call flush (output_unit)
1250# 411 "/home/runner/work/MFC/MFC/src/simulation/m_thinc.fpp"
1251 end block
1252# 411 "/home/runner/work/MFC/MFC/src/simulation/m_thinc.fpp"
1253#endif
1254# 411 "/home/runner/work/MFC/MFC/src/simulation/m_thinc.fpp"
1255
1256# 411 "/home/runner/work/MFC/MFC/src/simulation/m_thinc.fpp"
1257#if defined(MFC_OpenACC)
1258# 411 "/home/runner/work/MFC/MFC/src/simulation/m_thinc.fpp"
1259!$acc exit data delete(mthinc_nhat)
1260# 411 "/home/runner/work/MFC/MFC/src/simulation/m_thinc.fpp"
1261#elif defined(MFC_OpenMP)
1262# 411 "/home/runner/work/MFC/MFC/src/simulation/m_thinc.fpp"
1263!$omp target exit data map(release:mthinc_nhat)
1264# 411 "/home/runner/work/MFC/MFC/src/simulation/m_thinc.fpp"
1265#endif
1266# 411 "/home/runner/work/MFC/MFC/src/simulation/m_thinc.fpp"
1267 deallocate (mthinc_nhat)
1268#ifdef MFC_DEBUG
1269# 412 "/home/runner/work/MFC/MFC/src/simulation/m_thinc.fpp"
1270 block
1271# 412 "/home/runner/work/MFC/MFC/src/simulation/m_thinc.fpp"
1272 use iso_fortran_env, only: output_unit
1273# 412 "/home/runner/work/MFC/MFC/src/simulation/m_thinc.fpp"
1274
1275# 412 "/home/runner/work/MFC/MFC/src/simulation/m_thinc.fpp"
1276 print *, 'm_thinc.fpp:412: ', '@:DEALLOCATE(mthinc_d)'
1277# 412 "/home/runner/work/MFC/MFC/src/simulation/m_thinc.fpp"
1278
1279# 412 "/home/runner/work/MFC/MFC/src/simulation/m_thinc.fpp"
1280 call flush (output_unit)
1281# 412 "/home/runner/work/MFC/MFC/src/simulation/m_thinc.fpp"
1282 end block
1283# 412 "/home/runner/work/MFC/MFC/src/simulation/m_thinc.fpp"
1284#endif
1285# 412 "/home/runner/work/MFC/MFC/src/simulation/m_thinc.fpp"
1286
1287# 412 "/home/runner/work/MFC/MFC/src/simulation/m_thinc.fpp"
1288#if defined(MFC_OpenACC)
1289# 412 "/home/runner/work/MFC/MFC/src/simulation/m_thinc.fpp"
1290!$acc exit data delete(mthinc_d)
1291# 412 "/home/runner/work/MFC/MFC/src/simulation/m_thinc.fpp"
1292#elif defined(MFC_OpenMP)
1293# 412 "/home/runner/work/MFC/MFC/src/simulation/m_thinc.fpp"
1294!$omp target exit data map(release:mthinc_d)
1295# 412 "/home/runner/work/MFC/MFC/src/simulation/m_thinc.fpp"
1296#endif
1297# 412 "/home/runner/work/MFC/MFC/src/simulation/m_thinc.fpp"
1298 deallocate (mthinc_d)
1299 end if
1300
1301 end subroutine s_finalize_thinc_module
1302
1303end module m_thinc
integer, intent(in) k
integer, intent(in) j
integer, intent(in) l
Compile-time constant parameters: default values, tolerances, and physical constants.
integer, parameter int_comp_mthinc
real(wp), parameter verysmall
Very small number.
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...
type(int_bounds_info), dimension(1:3) idwbuff
Utility routines for bubble model setup, coordinate transforms, array sampling, and special functions...
THINC and MTHINC interface compression for volume fraction sharpening. THINC (int_comp=1): 1D directi...
real(wp) function f_thinc_integral_1d(a, b)
Analytical 1-D integral of the THINC function.
real(wp), dimension(3) gq3_pts
3-point Gauss-Legendre quadrature on [-1/2, 1/2] Node locations: +-sqrt(3/5)/2, 0
subroutine, public s_finalize_thinc_module()
real(wp), dimension(:,:,:), allocatable mthinc_d
subroutine, public s_initialize_thinc_module()
real(wp) function f_mthinc_solve_d(n1, n2, n3, beta, alpha_cell, ndim)
Solve for the interface-position parameter d.
real(wp), dimension(:,:,:,:), allocatable mthinc_nhat
real(wp) function f_mthinc_volume_integral(n1, n2, n3, d, beta, ndim)
Volume integral of H(xi) = 0.5*(1 + tanh(beta*(n.xi + d))) over the cell [-1/2, 1/2]^ndim.
subroutine, public s_compute_mthinc_normals(v_vf)
Compute the unit normal and interface-position parameter d at each interface cell.
subroutine, public s_thinc_compression(v_rs_ws, vl_rs_vf_x, vr_rs_vf_x, recon_dir, is1_d, is2_d, is3_d)
Applies THINC (int_comp=1) or MTHINC (int_comp=2) interface compression to sharpen volume-fraction an...
real(wp), dimension(3) gq3_wts
Weights: 5/18, 8/18, 5/18.
real(wp) function f_log_cosh_diff(a, h)
Stable difference: log_cosh(a+h) - log_cosh(a-h) = 2*atanh(tanh(a)*tanh(h)). Avoids catastrophic canc...
real(wp) function f_mthinc_face_average(n1, n2, n3, d, beta, face_dir, face_pos, ndim)
Face-averaged THINC function at a cell face. face_dir: 1=x, 2=y, 3=z. face_pos: -0....
real(wp) function f_mthinc_volume_integral_dd(n1, n2, n3, d, beta, ndim)
Derivative dV/dd of the volume integral (for Newton iteration).