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