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