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