MFC
Exascale flow solver
Loading...
Searching...
No Matches
m_compute_levelset.fpp.f90
Go to the documentation of this file.
1# 1 "/home/runner/work/MFC/MFC/src/simulation/m_compute_levelset.fpp"
2!>
3!!@file
4!! @brief Contains module m_compute_levelset
5
6# 1 "/home/runner/work/MFC/MFC/src/common/include/macros.fpp" 1
7# 1 "/home/runner/work/MFC/MFC/src/common/include/parallel_macros.fpp" 1
8# 1 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp" 1
9# 2 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
10# 3 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
11# 4 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
12# 5 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
13# 6 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
14
15# 8 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
16# 9 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
17# 10 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
18
19# 17 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
20
21# 46 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
22
23# 58 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
24
25# 68 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
26
27# 98 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
28
29# 110 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
30
31# 120 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
32
33# 145 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
34! New line at end of file is required for FYPP
35# 2 "/home/runner/work/MFC/MFC/src/common/include/parallel_macros.fpp" 2
36# 1 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp" 1
37# 1 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp" 1
38# 2 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
39# 3 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
40# 4 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
41# 5 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
42# 6 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
43
44# 8 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
45# 9 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
46# 10 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
47
48# 17 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
49
50# 46 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
51
52# 58 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
53
54# 68 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
55
56# 98 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
57
58# 110 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
59
60# 120 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
61
62# 145 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
63! New line at end of file is required for FYPP
64# 2 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp" 2
65
66# 4 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
67# 5 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
68# 6 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
69# 7 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
70# 8 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
71
72# 20 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
73
74# 43 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
75
76# 48 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
77
78# 53 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
79
80# 58 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
81
82# 63 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
83
84# 68 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
85
86# 76 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
87
88# 81 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
89
90# 86 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
91
92# 91 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
93
94# 96 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
95
96# 101 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
97
98# 106 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
99
100# 111 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
101
102# 116 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
103
104# 121 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
105
106# 151 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
107
108# 192 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
109
110# 206 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
111
112# 231 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
113
114# 242 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
115
116# 244 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
117# 255 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
118
119# 284 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
120
121# 294 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
122
123# 304 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
124
125# 313 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
126
127# 330 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
128
129# 340 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
130
131# 347 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
132
133# 353 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
134
135# 359 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
136
137# 365 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
138
139# 371 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
140
141# 377 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
142! New line at end of file is required for FYPP
143# 3 "/home/runner/work/MFC/MFC/src/common/include/parallel_macros.fpp" 2
144# 1 "/home/runner/work/MFC/MFC/src/common/include/acc_macros.fpp" 1
145# 1 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp" 1
146# 2 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
147# 3 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
148# 4 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
149# 5 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
150# 6 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
151
152# 8 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
153# 9 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
154# 10 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
155
156# 17 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
157
158# 46 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
159
160# 58 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
161
162# 68 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
163
164# 98 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
165
166# 110 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
167
168# 120 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
169
170# 145 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
171! New line at end of file is required for FYPP
172# 2 "/home/runner/work/MFC/MFC/src/common/include/acc_macros.fpp" 2
173
174# 7 "/home/runner/work/MFC/MFC/src/common/include/acc_macros.fpp"
175
176# 17 "/home/runner/work/MFC/MFC/src/common/include/acc_macros.fpp"
177
178# 22 "/home/runner/work/MFC/MFC/src/common/include/acc_macros.fpp"
179
180# 27 "/home/runner/work/MFC/MFC/src/common/include/acc_macros.fpp"
181
182# 32 "/home/runner/work/MFC/MFC/src/common/include/acc_macros.fpp"
183
184# 37 "/home/runner/work/MFC/MFC/src/common/include/acc_macros.fpp"
185
186# 42 "/home/runner/work/MFC/MFC/src/common/include/acc_macros.fpp"
187
188# 47 "/home/runner/work/MFC/MFC/src/common/include/acc_macros.fpp"
189
190# 52 "/home/runner/work/MFC/MFC/src/common/include/acc_macros.fpp"
191
192# 57 "/home/runner/work/MFC/MFC/src/common/include/acc_macros.fpp"
193
194# 62 "/home/runner/work/MFC/MFC/src/common/include/acc_macros.fpp"
195
196# 73 "/home/runner/work/MFC/MFC/src/common/include/acc_macros.fpp"
197
198# 78 "/home/runner/work/MFC/MFC/src/common/include/acc_macros.fpp"
199
200# 83 "/home/runner/work/MFC/MFC/src/common/include/acc_macros.fpp"
201
202# 88 "/home/runner/work/MFC/MFC/src/common/include/acc_macros.fpp"
203
204# 103 "/home/runner/work/MFC/MFC/src/common/include/acc_macros.fpp"
205
206# 131 "/home/runner/work/MFC/MFC/src/common/include/acc_macros.fpp"
207
208# 160 "/home/runner/work/MFC/MFC/src/common/include/acc_macros.fpp"
209
210# 175 "/home/runner/work/MFC/MFC/src/common/include/acc_macros.fpp"
211
212# 193 "/home/runner/work/MFC/MFC/src/common/include/acc_macros.fpp"
213
214# 215 "/home/runner/work/MFC/MFC/src/common/include/acc_macros.fpp"
215
216# 244 "/home/runner/work/MFC/MFC/src/common/include/acc_macros.fpp"
217
218# 259 "/home/runner/work/MFC/MFC/src/common/include/acc_macros.fpp"
219
220# 269 "/home/runner/work/MFC/MFC/src/common/include/acc_macros.fpp"
221
222# 278 "/home/runner/work/MFC/MFC/src/common/include/acc_macros.fpp"
223
224# 294 "/home/runner/work/MFC/MFC/src/common/include/acc_macros.fpp"
225
226# 304 "/home/runner/work/MFC/MFC/src/common/include/acc_macros.fpp"
227
228# 311 "/home/runner/work/MFC/MFC/src/common/include/acc_macros.fpp"
229! New line at end of file is required for FYPP
230# 4 "/home/runner/work/MFC/MFC/src/common/include/parallel_macros.fpp" 2
231
232! GPU parallel region (scalar reductions, maxval/minval)
233# 23 "/home/runner/work/MFC/MFC/src/common/include/parallel_macros.fpp"
234
235! GPU parallel loop over threads (most common GPU macro)
236# 43 "/home/runner/work/MFC/MFC/src/common/include/parallel_macros.fpp"
237
238! Required closing for GPU_PARALLEL_LOOP
239# 55 "/home/runner/work/MFC/MFC/src/common/include/parallel_macros.fpp"
240
241! Mark routine for device compilation
242# 112 "/home/runner/work/MFC/MFC/src/common/include/parallel_macros.fpp"
243
244! Declare device-resident data
245# 130 "/home/runner/work/MFC/MFC/src/common/include/parallel_macros.fpp"
246
247! Inner loop within a GPU parallel region
248# 145 "/home/runner/work/MFC/MFC/src/common/include/parallel_macros.fpp"
249
250! Scoped GPU data region
251# 164 "/home/runner/work/MFC/MFC/src/common/include/parallel_macros.fpp"
252
253! Host code with device pointers (for MPI with GPU buffers)
254# 193 "/home/runner/work/MFC/MFC/src/common/include/parallel_macros.fpp"
255
256! Allocate device memory (unscoped)
257# 207 "/home/runner/work/MFC/MFC/src/common/include/parallel_macros.fpp"
258
259! Free device memory
260# 219 "/home/runner/work/MFC/MFC/src/common/include/parallel_macros.fpp"
261
262! Atomic operation on device
263# 231 "/home/runner/work/MFC/MFC/src/common/include/parallel_macros.fpp"
264
265! End atomic capture block
266# 242 "/home/runner/work/MFC/MFC/src/common/include/parallel_macros.fpp"
267
268! Copy data between host and device
269# 254 "/home/runner/work/MFC/MFC/src/common/include/parallel_macros.fpp"
270
271! Synchronization barrier
272# 266 "/home/runner/work/MFC/MFC/src/common/include/parallel_macros.fpp"
273
274! Import GPU library module (openacc or omp_lib)
275# 275 "/home/runner/work/MFC/MFC/src/common/include/parallel_macros.fpp"
276
277! Emit code only for AMD compiler
278# 282 "/home/runner/work/MFC/MFC/src/common/include/parallel_macros.fpp"
279
280! Emit code for non-Cray compilers
281# 289 "/home/runner/work/MFC/MFC/src/common/include/parallel_macros.fpp"
282
283! Emit code only for Cray compiler
284# 296 "/home/runner/work/MFC/MFC/src/common/include/parallel_macros.fpp"
285
286! Emit code for non-NVIDIA compilers
287# 303 "/home/runner/work/MFC/MFC/src/common/include/parallel_macros.fpp"
288
289# 305 "/home/runner/work/MFC/MFC/src/common/include/parallel_macros.fpp"
290# 306 "/home/runner/work/MFC/MFC/src/common/include/parallel_macros.fpp"
291! New line at end of file is required for FYPP
292# 2 "/home/runner/work/MFC/MFC/src/common/include/macros.fpp" 2
293
294# 14 "/home/runner/work/MFC/MFC/src/common/include/macros.fpp"
295
296! Caution: This macro requires the use of a binding script to set CUDA_VISIBLE_DEVICES, such that we have one GPU device per MPI
297! rank. That's because for both cudaMemAdvise (preferred location) and cudaMemPrefetchAsync we use location = device_id = 0. For an
298! example see misc/nvidia_uvm/bind.sh. NVIDIA unified memory page placement hint
299# 57 "/home/runner/work/MFC/MFC/src/common/include/macros.fpp"
300
301! Allocate and create GPU device memory
302# 77 "/home/runner/work/MFC/MFC/src/common/include/macros.fpp"
303
304! Free GPU device memory and deallocate
305# 85 "/home/runner/work/MFC/MFC/src/common/include/macros.fpp"
306
307! Cray-specific GPU pointer setup for vector fields
308# 109 "/home/runner/work/MFC/MFC/src/common/include/macros.fpp"
309
310! Cray-specific GPU pointer setup for scalar fields
311# 125 "/home/runner/work/MFC/MFC/src/common/include/macros.fpp"
312
313! Cray-specific GPU pointer setup for acoustic source spatials
314# 150 "/home/runner/work/MFC/MFC/src/common/include/macros.fpp"
315
316# 156 "/home/runner/work/MFC/MFC/src/common/include/macros.fpp"
317
318# 163 "/home/runner/work/MFC/MFC/src/common/include/macros.fpp"
319! New line at end of file is required for FYPP
320# 6 "/home/runner/work/MFC/MFC/src/simulation/m_compute_levelset.fpp" 2
321
322!> @brief Computes signed-distance level-set fields and surface normals for immersed-boundary patch geometries
323module m_compute_levelset
324
325 use m_ib_patches
326 use m_model
327 use m_derived_types
328 use m_global_parameters
329 use m_mpi_proxy
330 use m_helper_basic
331
332 implicit none
333
334 private; public :: s_apply_levelset
335
336contains
337
338 !> Dispatch level-set distance and normal computations for all ghost points based on patch geometry type
339 impure subroutine s_apply_levelset(gps, num_gps)
340
341 type(ghost_point), dimension(:), intent(inout) :: gps
342 integer, intent(in) :: num_gps
343 integer :: i, patch_id, patch_geometry
344
345 ! 3D Patch Geometries
346
347 if (p > 0) then
348
349# 33 "/home/runner/work/MFC/MFC/src/simulation/m_compute_levelset.fpp"
350
351# 33 "/home/runner/work/MFC/MFC/src/simulation/m_compute_levelset.fpp"
352#if defined(MFC_OpenACC)
353# 33 "/home/runner/work/MFC/MFC/src/simulation/m_compute_levelset.fpp"
354!$acc parallel loop gang vector default(present) private(i, patch_id, patch_geometry) copy(gps) copyin(patch_ib(1:num_ibs) , ib_airfoil, ib_airfoil_grids)
355# 33 "/home/runner/work/MFC/MFC/src/simulation/m_compute_levelset.fpp"
356#elif defined(MFC_OpenMP)
357# 33 "/home/runner/work/MFC/MFC/src/simulation/m_compute_levelset.fpp"
358
359# 33 "/home/runner/work/MFC/MFC/src/simulation/m_compute_levelset.fpp"
360
361# 33 "/home/runner/work/MFC/MFC/src/simulation/m_compute_levelset.fpp"
362
363# 33 "/home/runner/work/MFC/MFC/src/simulation/m_compute_levelset.fpp"
364!$omp target teams loop defaultmap(firstprivate:scalar) bind(teams,parallel) defaultmap(tofrom:aggregate) defaultmap(tofrom:allocatable) defaultmap(tofrom:pointer) private(i, patch_id, patch_geometry) &
365# 33 "/home/runner/work/MFC/MFC/src/simulation/m_compute_levelset.fpp"
366!$omp& map(tofrom:gps) map(to:patch_ib(1:num_ibs) , ib_airfoil, ib_airfoil_grids)
367# 33 "/home/runner/work/MFC/MFC/src/simulation/m_compute_levelset.fpp"
368#endif
369# 35 "/home/runner/work/MFC/MFC/src/simulation/m_compute_levelset.fpp"
370 do i = 1, num_gps
371 patch_id = gps(i)%ib_patch_id
372 patch_geometry = patch_ib(patch_id)%geometry
373
374 if (patch_geometry == 8) then
375 call s_sphere_levelset(gps(i))
376 else if (patch_geometry == 9) then
377 call s_cuboid_levelset(gps(i))
378 else if (patch_geometry == 10) then
379 call s_cylinder_levelset(gps(i))
380 else if (patch_geometry == 11) then
381 call s_3d_airfoil_levelset(gps(i))
382 else if (patch_geometry == 12) then
383 call s_model_levelset(gps(i))
384 end if
385 end do
386
387# 51 "/home/runner/work/MFC/MFC/src/simulation/m_compute_levelset.fpp"
388#if defined(MFC_OpenACC)
389# 51 "/home/runner/work/MFC/MFC/src/simulation/m_compute_levelset.fpp"
390!$acc end parallel loop
391# 51 "/home/runner/work/MFC/MFC/src/simulation/m_compute_levelset.fpp"
392#elif defined(MFC_OpenMP)
393# 51 "/home/runner/work/MFC/MFC/src/simulation/m_compute_levelset.fpp"
394
395# 51 "/home/runner/work/MFC/MFC/src/simulation/m_compute_levelset.fpp"
396!$omp end target teams loop
397# 51 "/home/runner/work/MFC/MFC/src/simulation/m_compute_levelset.fpp"
398#endif
399
400 ! 2D Patch Geometries
401 else if (n > 0) then
402
403# 55 "/home/runner/work/MFC/MFC/src/simulation/m_compute_levelset.fpp"
404
405# 55 "/home/runner/work/MFC/MFC/src/simulation/m_compute_levelset.fpp"
406#if defined(MFC_OpenACC)
407# 55 "/home/runner/work/MFC/MFC/src/simulation/m_compute_levelset.fpp"
408!$acc parallel loop gang vector default(present) private(i, patch_id, patch_geometry) copy(gps) copyin(patch_ib(1:num_ibs) , ib_airfoil, ib_airfoil_grids)
409# 55 "/home/runner/work/MFC/MFC/src/simulation/m_compute_levelset.fpp"
410#elif defined(MFC_OpenMP)
411# 55 "/home/runner/work/MFC/MFC/src/simulation/m_compute_levelset.fpp"
412
413# 55 "/home/runner/work/MFC/MFC/src/simulation/m_compute_levelset.fpp"
414
415# 55 "/home/runner/work/MFC/MFC/src/simulation/m_compute_levelset.fpp"
416
417# 55 "/home/runner/work/MFC/MFC/src/simulation/m_compute_levelset.fpp"
418!$omp target teams loop defaultmap(firstprivate:scalar) bind(teams,parallel) defaultmap(tofrom:aggregate) defaultmap(tofrom:allocatable) defaultmap(tofrom:pointer) private(i, patch_id, patch_geometry) &
419# 55 "/home/runner/work/MFC/MFC/src/simulation/m_compute_levelset.fpp"
420!$omp& map(tofrom:gps) map(to:patch_ib(1:num_ibs) , ib_airfoil, ib_airfoil_grids)
421# 55 "/home/runner/work/MFC/MFC/src/simulation/m_compute_levelset.fpp"
422#endif
423# 57 "/home/runner/work/MFC/MFC/src/simulation/m_compute_levelset.fpp"
424 do i = 1, num_gps
425 patch_id = gps(i)%ib_patch_id
426 patch_geometry = patch_ib(patch_id)%geometry
427
428 if (patch_geometry == 2) then
429 call s_circle_levelset(gps(i))
430 else if (patch_geometry == 3) then
431 call s_rectangle_levelset(gps(i))
432 else if (patch_geometry == 4) then
433 call s_airfoil_levelset(gps(i))
434 else if (patch_geometry == 5) then
435 call s_model_levelset(gps(i))
436 else if (patch_geometry == 6) then
437 call s_ellipse_levelset(gps(i))
438 end if
439 end do
440
441# 73 "/home/runner/work/MFC/MFC/src/simulation/m_compute_levelset.fpp"
442#if defined(MFC_OpenACC)
443# 73 "/home/runner/work/MFC/MFC/src/simulation/m_compute_levelset.fpp"
444!$acc end parallel loop
445# 73 "/home/runner/work/MFC/MFC/src/simulation/m_compute_levelset.fpp"
446#elif defined(MFC_OpenMP)
447# 73 "/home/runner/work/MFC/MFC/src/simulation/m_compute_levelset.fpp"
448
449# 73 "/home/runner/work/MFC/MFC/src/simulation/m_compute_levelset.fpp"
450!$omp end target teams loop
451# 73 "/home/runner/work/MFC/MFC/src/simulation/m_compute_levelset.fpp"
452#endif
453 end if
454
455 end subroutine s_apply_levelset
456
457 !> Compute the signed distance and outward normal from a ghost point to a circular immersed boundary
458 subroutine s_circle_levelset(gp)
459
460
461# 81 "/home/runner/work/MFC/MFC/src/simulation/m_compute_levelset.fpp"
462#if MFC_OpenACC
463# 81 "/home/runner/work/MFC/MFC/src/simulation/m_compute_levelset.fpp"
464!$acc routine seq
465# 81 "/home/runner/work/MFC/MFC/src/simulation/m_compute_levelset.fpp"
466#elif MFC_OpenMP
467# 81 "/home/runner/work/MFC/MFC/src/simulation/m_compute_levelset.fpp"
468
469# 81 "/home/runner/work/MFC/MFC/src/simulation/m_compute_levelset.fpp"
470
471# 81 "/home/runner/work/MFC/MFC/src/simulation/m_compute_levelset.fpp"
472!$omp declare target device_type(any)
473# 81 "/home/runner/work/MFC/MFC/src/simulation/m_compute_levelset.fpp"
474#endif
475
476 type(ghost_point), intent(inout) :: gp
477 real(wp) :: radius, dist
478 real(wp), dimension(3) :: dist_vec
479 integer :: i, j, ib_patch_id !< Loop index variables
480 ib_patch_id = gp%ib_patch_id
481 i = gp%loc(1)
482 j = gp%loc(2)
483
484 radius = patch_ib(ib_patch_id)%radius
485
486 dist_vec(1) = x_cc(i) - (patch_ib(ib_patch_id)%x_centroid + real(gp%x_periodicity, wp)*(x_domain%end - x_domain%beg))
487 dist_vec(2) = y_cc(j) - (patch_ib(ib_patch_id)%y_centroid + real(gp%y_periodicity, wp)*(y_domain%end - y_domain%beg))
488 dist_vec(3) = 0._wp
489 dist = sqrt(sum(dist_vec**2))
490
491 gp%levelset = dist - radius
492 if (f_approx_equal(dist, 0._wp)) then
493 gp%levelset_norm = 0._wp
494 else
495 gp%levelset_norm = dist_vec(:)/dist
496 end if
497
498 end subroutine s_circle_levelset
499
500 !> Compute the signed distance and outward normal from a ghost point to a 2D NACA airfoil surface
501 subroutine s_airfoil_levelset(gp)
502
503
504# 110 "/home/runner/work/MFC/MFC/src/simulation/m_compute_levelset.fpp"
505#if MFC_OpenACC
506# 110 "/home/runner/work/MFC/MFC/src/simulation/m_compute_levelset.fpp"
507!$acc routine seq
508# 110 "/home/runner/work/MFC/MFC/src/simulation/m_compute_levelset.fpp"
509#elif MFC_OpenMP
510# 110 "/home/runner/work/MFC/MFC/src/simulation/m_compute_levelset.fpp"
511
512# 110 "/home/runner/work/MFC/MFC/src/simulation/m_compute_levelset.fpp"
513
514# 110 "/home/runner/work/MFC/MFC/src/simulation/m_compute_levelset.fpp"
515!$omp declare target device_type(any)
516# 110 "/home/runner/work/MFC/MFC/src/simulation/m_compute_levelset.fpp"
517#endif
518
519 type(ghost_point), intent(inout) :: gp
520 real(wp) :: dist, global_dist
521 integer :: global_id, airfoil_id, Np_local
522 real(wp), dimension(3) :: dist_vec
523 real(wp), dimension(1:3) :: xy_local, offset !< x and y coordinates in local IB frame
524 real(wp), dimension(1:2) :: center
525 real(wp), dimension(1:3,1:3) :: rotation, inverse_rotation
526 integer :: i, j, k, ib_patch_id !< Loop index variables
527 ib_patch_id = gp%ib_patch_id
528 i = gp%loc(1)
529 j = gp%loc(2)
530
531 airfoil_id = patch_ib(ib_patch_id)%airfoil_id
532 np_local = ib_airfoil_grids(airfoil_id)%Np
533 center(1) = patch_ib(ib_patch_id)%x_centroid + real(gp%x_periodicity, wp)*(x_domain%end - x_domain%beg)
534 center(2) = patch_ib(ib_patch_id)%y_centroid + real(gp%y_periodicity, wp)*(y_domain%end - y_domain%beg)
535 inverse_rotation(:,:) = patch_ib(ib_patch_id)%rotation_matrix_inverse(:,:)
536 rotation(:,:) = patch_ib(ib_patch_id)%rotation_matrix(:,:)
537 offset(:) = patch_ib(ib_patch_id)%centroid_offset(:)
538
539 xy_local = [x_cc(i) - center(1), y_cc(j) - center(2), 0._wp] ! get coordinate frame centered on IB
540 xy_local = matmul(inverse_rotation, xy_local) ! rotate the frame into the IB's coordinate
541 xy_local = xy_local - offset ! airfoils are a patch that require a centroid offset
542
543 if (xy_local(2) >= 0._wp) then
544 ! finds the location on the airfoil grid with the minimum distance (closest)
545 do k = 1, np_local
546 dist_vec(1) = ib_airfoil_grids(airfoil_id)%upper(k)%x - xy_local(1)
547 dist_vec(2) = ib_airfoil_grids(airfoil_id)%upper(k)%y - xy_local(2)
548 dist_vec(3) = 0._wp
549 dist = sqrt(sum(dist_vec**2))
550 if (k == 1) then
551 global_dist = dist
552 global_id = k
553 else
554 if (dist < global_dist) then
555 global_dist = dist
556 global_id = k
557 end if
558 end if
559 end do
560 dist_vec(1) = ib_airfoil_grids(airfoil_id)%upper(global_id)%x - xy_local(1)
561 dist_vec(2) = ib_airfoil_grids(airfoil_id)%upper(global_id)%y - xy_local(2)
562 dist_vec(3) = 0
563 dist = global_dist
564 else
565 do k = 1, np_local
566 dist_vec(1) = ib_airfoil_grids(airfoil_id)%lower(k)%x - xy_local(1)
567 dist_vec(2) = ib_airfoil_grids(airfoil_id)%lower(k)%y - xy_local(2)
568 dist_vec(3) = 0
569 dist = sqrt(sum(dist_vec**2))
570 if (k == 1) then
571 global_dist = dist
572 global_id = k
573 else
574 if (dist < global_dist) then
575 global_dist = dist
576 global_id = k
577 end if
578 end if
579 end do
580 dist_vec(1) = ib_airfoil_grids(airfoil_id)%lower(global_id)%x - xy_local(1)
581 dist_vec(2) = ib_airfoil_grids(airfoil_id)%lower(global_id)%y - xy_local(2)
582 dist_vec(3) = 0._wp
583 dist = global_dist
584 end if
585
586 gp%levelset = dist
587 if (f_approx_equal(dist, 0._wp)) then
588 gp%levelset_norm = 0._wp
589 else
590 gp%levelset_norm = matmul(rotation, dist_vec(:))/dist ! convert the normal vector back to global grid coordinates
591 end if
592
593 end subroutine s_airfoil_levelset
594
595 !> Compute the signed distance and outward normal from a ghost point to a 3D extruded airfoil surface
596 subroutine s_3d_airfoil_levelset(gp)
597
598
599# 191 "/home/runner/work/MFC/MFC/src/simulation/m_compute_levelset.fpp"
600#if MFC_OpenACC
601# 191 "/home/runner/work/MFC/MFC/src/simulation/m_compute_levelset.fpp"
602!$acc routine seq
603# 191 "/home/runner/work/MFC/MFC/src/simulation/m_compute_levelset.fpp"
604#elif MFC_OpenMP
605# 191 "/home/runner/work/MFC/MFC/src/simulation/m_compute_levelset.fpp"
606
607# 191 "/home/runner/work/MFC/MFC/src/simulation/m_compute_levelset.fpp"
608
609# 191 "/home/runner/work/MFC/MFC/src/simulation/m_compute_levelset.fpp"
610!$omp declare target device_type(any)
611# 191 "/home/runner/work/MFC/MFC/src/simulation/m_compute_levelset.fpp"
612#endif
613
614 type(ghost_point), intent(inout) :: gp
615 real(wp) :: dist_surf, dist_side, global_dist
616 integer :: global_id, airfoil_id, Np_local
617 real(wp) :: lz, z_max, z_min
618 real(wp), dimension(3) :: dist_vec
619 real(wp), dimension(1:3) :: xyz_local, center, offset, normal !< x, y, z coordinates in local IB frame
620 real(wp), dimension(1:3,1:3) :: rotation, inverse_rotation
621 integer :: i, j, k, l, ib_patch_id !< Loop index variables
622 ib_patch_id = gp%ib_patch_id
623 i = gp%loc(1)
624 j = gp%loc(2)
625 l = gp%loc(3)
626
627 airfoil_id = patch_ib(ib_patch_id)%airfoil_id
628 np_local = ib_airfoil_grids(airfoil_id)%Np
629 center(1) = patch_ib(ib_patch_id)%x_centroid + real(gp%x_periodicity, wp)*(x_domain%end - x_domain%beg)
630 center(2) = patch_ib(ib_patch_id)%y_centroid + real(gp%y_periodicity, wp)*(y_domain%end - y_domain%beg)
631 center(3) = patch_ib(ib_patch_id)%z_centroid + real(gp%z_periodicity, wp)*(z_domain%end - z_domain%beg)
632 lz = patch_ib(ib_patch_id)%length_z
633 inverse_rotation(:,:) = patch_ib(ib_patch_id)%rotation_matrix_inverse(:,:)
634 rotation(:,:) = patch_ib(ib_patch_id)%rotation_matrix(:,:)
635 offset(:) = patch_ib(ib_patch_id)%centroid_offset(:)
636
637 z_max = lz/2
638 z_min = -lz/2
639
640 xyz_local = [x_cc(i), y_cc(j), z_cc(l)] - center
641 xyz_local = matmul(inverse_rotation, xyz_local) ! rotate the frame into the IB's coordinates
642 xyz_local = xyz_local - offset ! airfoils are a patch that require a centroid offset
643
644 if (xyz_local(2) >= 0._wp) then
645 do k = 1, np_local
646 dist_vec(1) = xyz_local(1) - ib_airfoil_grids(airfoil_id)%upper(k)%x
647 dist_vec(2) = xyz_local(2) - ib_airfoil_grids(airfoil_id)%upper(k)%y
648 dist_vec(3) = 0._wp
649 dist_surf = sqrt(sum(dist_vec**2))
650 if (k == 1) then
651 global_dist = dist_surf
652 global_id = k
653 else
654 if (dist_surf < global_dist) then
655 global_dist = dist_surf
656 global_id = k
657 end if
658 end if
659 end do
660 dist_vec(1) = ib_airfoil_grids(airfoil_id)%upper(global_id)%x - xyz_local(1)
661 dist_vec(2) = ib_airfoil_grids(airfoil_id)%upper(global_id)%y - xyz_local(2)
662 dist_vec(3) = 0._wp
663 dist_surf = global_dist
664 else
665 do k = 1, np_local
666 dist_vec(1) = ib_airfoil_grids(airfoil_id)%lower(k)%x - xyz_local(1)
667 dist_vec(2) = ib_airfoil_grids(airfoil_id)%lower(k)%y - xyz_local(2)
668 dist_vec(3) = 0
669 dist_surf = sqrt(sum(dist_vec**2))
670 if (k == 1) then
671 global_dist = dist_surf
672 global_id = k
673 else
674 if (dist_surf < global_dist) then
675 global_dist = dist_surf
676 global_id = k
677 end if
678 end if
679 end do
680 dist_vec(1) = ib_airfoil_grids(airfoil_id)%lower(global_id)%x - xyz_local(1)
681 dist_vec(2) = ib_airfoil_grids(airfoil_id)%lower(global_id)%y - xyz_local(2)
682 dist_vec(3) = 0._wp
683 dist_surf = global_dist
684 end if
685
686 dist_side = min(abs(xyz_local(3) - z_min), abs(z_max - xyz_local(3)))
687
688 if (dist_side < dist_surf) then
689 gp%levelset = dist_side
690 normal = 0._wp
691 if (f_approx_equal(dist_side, abs(xyz_local(3) - z_min))) then
692 normal(3) = -1._wp
693 else
694 normal(3) = 1._wp
695 end if
696 gp%levelset_norm = matmul(rotation, normal)
697 else
698 gp%levelset = dist_surf
699 if (f_approx_equal(dist_surf, 0._wp)) then
700 gp%levelset_norm = 0._wp
701 else
702 gp%levelset_norm = matmul(rotation, dist_vec(:)/dist_surf)
703 end if
704 end if
705
706 end subroutine s_3d_airfoil_levelset
707
708 !> Subroutine for computing the levelset values at a ghost point belonging to the rectangle IB
709 subroutine s_rectangle_levelset(gp)
710
711
712# 290 "/home/runner/work/MFC/MFC/src/simulation/m_compute_levelset.fpp"
713#if MFC_OpenACC
714# 290 "/home/runner/work/MFC/MFC/src/simulation/m_compute_levelset.fpp"
715!$acc routine seq
716# 290 "/home/runner/work/MFC/MFC/src/simulation/m_compute_levelset.fpp"
717#elif MFC_OpenMP
718# 290 "/home/runner/work/MFC/MFC/src/simulation/m_compute_levelset.fpp"
719
720# 290 "/home/runner/work/MFC/MFC/src/simulation/m_compute_levelset.fpp"
721
722# 290 "/home/runner/work/MFC/MFC/src/simulation/m_compute_levelset.fpp"
723!$omp declare target device_type(any)
724# 290 "/home/runner/work/MFC/MFC/src/simulation/m_compute_levelset.fpp"
725#endif
726
727 type(ghost_point), intent(inout) :: gp
728 real(wp) :: top_right(2), bottom_left(2)
729 real(wp) :: min_dist
730 real(wp) :: side_dists(4)
731 real(wp) :: length_x, length_y
732 real(wp), dimension(1:3) :: xy_local, dist_vec !< x and y coordinates in local IB frame
733 real(wp), dimension(2) :: center !< x and y coordinates in local IB frame
734 real(wp), dimension(1:3,1:3) :: rotation, inverse_rotation
735 integer :: i, j, k !< Loop index variables
736 integer :: idx !< Shortest path direction indicator
737 integer :: ib_patch_id !< patch ID
738 ib_patch_id = gp%ib_patch_id
739 i = gp%loc(1)
740 j = gp%loc(2)
741
742 length_x = patch_ib(ib_patch_id)%length_x
743 length_y = patch_ib(ib_patch_id)%length_y
744 center(1) = patch_ib(ib_patch_id)%x_centroid + real(gp%x_periodicity, wp)*(x_domain%end - x_domain%beg)
745 center(2) = patch_ib(ib_patch_id)%y_centroid + real(gp%y_periodicity, wp)*(y_domain%end - y_domain%beg)
746 inverse_rotation(:,:) = patch_ib(ib_patch_id)%rotation_matrix_inverse(:,:)
747 rotation(:,:) = patch_ib(ib_patch_id)%rotation_matrix(:,:)
748
749 top_right(1) = length_x/2
750 top_right(2) = length_y/2
751 bottom_left(1) = -length_x/2
752 bottom_left(2) = -length_y/2
753
754 ! convert grid to local coordinates
755 xy_local = [x_cc(i) - center(1), y_cc(j) - center(2), 0._wp]
756 xy_local = matmul(inverse_rotation, xy_local)
757
758 side_dists(1) = bottom_left(1) - xy_local(1)
759 side_dists(2) = top_right(1) - xy_local(1)
760 side_dists(3) = bottom_left(2) - xy_local(2)
761 side_dists(4) = top_right(2) - xy_local(2)
762 min_dist = side_dists(1)
763 idx = 1
764
765 do k = 2, 4
766 if (abs(side_dists(k)) < abs(min_dist)) then
767 idx = k
768 min_dist = side_dists(idx)
769 end if
770 end do
771
772 gp%levelset = side_dists(idx)
773 dist_vec = 0._wp
774 if (.not. f_approx_equal(side_dists(idx), 0._wp)) then
775 if (idx == 1 .or. idx == 2) then
776 ! vector points along the x axis
777 dist_vec(1) = side_dists(idx)/abs(side_dists(idx))
778 else
779 ! vector points along the y axis
780 dist_vec(2) = side_dists(idx)/abs(side_dists(idx))
781 end if
782 ! convert the normal vector back into the global coordinate system
783 gp%levelset_norm = matmul(rotation, dist_vec)
784 else
785 gp%levelset_norm = 0._wp
786 end if
787
788 end subroutine s_rectangle_levelset
789
790 !> Compute the signed distance and outward normal from a ghost point to an elliptical immersed boundary
791 subroutine s_ellipse_levelset(gp)
792
793
794# 358 "/home/runner/work/MFC/MFC/src/simulation/m_compute_levelset.fpp"
795#if MFC_OpenACC
796# 358 "/home/runner/work/MFC/MFC/src/simulation/m_compute_levelset.fpp"
797!$acc routine seq
798# 358 "/home/runner/work/MFC/MFC/src/simulation/m_compute_levelset.fpp"
799#elif MFC_OpenMP
800# 358 "/home/runner/work/MFC/MFC/src/simulation/m_compute_levelset.fpp"
801
802# 358 "/home/runner/work/MFC/MFC/src/simulation/m_compute_levelset.fpp"
803
804# 358 "/home/runner/work/MFC/MFC/src/simulation/m_compute_levelset.fpp"
805!$omp declare target device_type(any)
806# 358 "/home/runner/work/MFC/MFC/src/simulation/m_compute_levelset.fpp"
807#endif
808
809 type(ghost_point), intent(inout) :: gp
810 real(wp) :: ellipse_coeffs(2) !< a and b in the ellipse equation
811 real(wp) :: quadratic_coeffs(3) !< A, B, C in the quadratic equation to compute levelset
812 real(wp) :: length_x, length_y
813 real(wp), dimension(1:3) :: xy_local, normal_vector !< x and y coordinates in local IB frame
814 real(wp), dimension(2) :: center !< x and y coordinates in local IB frame
815 real(wp), dimension(1:3,1:3) :: rotation, inverse_rotation
816 integer :: i, j !< Loop index variables
817 integer :: ib_patch_id !< patch ID
818 ib_patch_id = gp%ib_patch_id
819 i = gp%loc(1)
820 j = gp%loc(2)
821
822 length_x = patch_ib(ib_patch_id)%length_x
823 length_y = patch_ib(ib_patch_id)%length_y
824 center(1) = patch_ib(ib_patch_id)%x_centroid + real(gp%x_periodicity, wp)*(x_domain%end - x_domain%beg)
825 center(2) = patch_ib(ib_patch_id)%y_centroid + real(gp%y_periodicity, wp)*(y_domain%end - y_domain%beg)
826 inverse_rotation(:,:) = patch_ib(ib_patch_id)%rotation_matrix_inverse(:,:)
827 rotation(:,:) = patch_ib(ib_patch_id)%rotation_matrix(:,:)
828
829 ellipse_coeffs(1) = 0.5_wp*length_x
830 ellipse_coeffs(2) = 0.5_wp*length_y
831
832 xy_local = [x_cc(i) - center(1), y_cc(j) - center(2), 0._wp]
833 xy_local = matmul(inverse_rotation, xy_local)
834
835 normal_vector = xy_local
836 ! get the normal direction via the coordinate transformation method
837 normal_vector(2) = normal_vector(2)*(ellipse_coeffs(1)/ellipse_coeffs(2))**2._wp
838 normal_vector = normal_vector/sqrt(dot_product(normal_vector, normal_vector)) ! normalize the vector
839 gp%levelset_norm = matmul(rotation, normal_vector) ! save after rotating the vector to the global frame
840
841 ! use the normal vector to set up the quadratic equation for the levelset, using A, B, and C in indices 1, 2, and 3
842 quadratic_coeffs(1) = (normal_vector(1)/ellipse_coeffs(1))**2 + (normal_vector(2)/ellipse_coeffs(2))**2
843 quadratic_coeffs(2) = 2._wp*((xy_local(1)*normal_vector(1)/(ellipse_coeffs(1)**2)) + (xy_local(2)*normal_vector(2) &
844 & /(ellipse_coeffs(2)**2)))
845 quadratic_coeffs(3) = (xy_local(1)/ellipse_coeffs(1))**2._wp + (xy_local(2)/ellipse_coeffs(2))**2._wp - 1._wp
846
847 ! compute the levelset with the quadratic equation [ -B + sqrt(B^2 - 4AC) ] / 2A
848 gp%levelset = -0.5_wp*(-quadratic_coeffs(2) + sqrt(quadratic_coeffs(2)**2._wp - 4._wp*quadratic_coeffs(1) &
849 & *quadratic_coeffs(3)))/quadratic_coeffs(1)
850
851 end subroutine s_ellipse_levelset
852
853 !> Compute the signed distance and outward normal from a ghost point to a cuboid immersed boundary
854 subroutine s_cuboid_levelset(gp)
855
856
857# 407 "/home/runner/work/MFC/MFC/src/simulation/m_compute_levelset.fpp"
858#if MFC_OpenACC
859# 407 "/home/runner/work/MFC/MFC/src/simulation/m_compute_levelset.fpp"
860!$acc routine seq
861# 407 "/home/runner/work/MFC/MFC/src/simulation/m_compute_levelset.fpp"
862#elif MFC_OpenMP
863# 407 "/home/runner/work/MFC/MFC/src/simulation/m_compute_levelset.fpp"
864
865# 407 "/home/runner/work/MFC/MFC/src/simulation/m_compute_levelset.fpp"
866
867# 407 "/home/runner/work/MFC/MFC/src/simulation/m_compute_levelset.fpp"
868!$omp declare target device_type(any)
869# 407 "/home/runner/work/MFC/MFC/src/simulation/m_compute_levelset.fpp"
870#endif
871
872 type(ghost_point), intent(inout) :: gp
873 real(wp) :: Right, Left, Bottom, Top, Front, Back
874 real(wp) :: min_dist
875 real(wp) :: dist_left, dist_right, dist_bottom, dist_top, dist_back, dist_front
876 real(wp), dimension(3) :: center
877 real(wp) :: length_x, length_y, length_z
878 real(wp), dimension(1:3) :: xyz_local, dist_vec !< x and y coordinates in local IB frame
879 real(wp), dimension(1:3,1:3) :: rotation, inverse_rotation
880 integer :: i, j, k !< Loop index variables
881 integer :: ib_patch_id !< patch ID
882 ib_patch_id = gp%ib_patch_id
883 i = gp%loc(1)
884 j = gp%loc(2)
885 k = gp%loc(3)
886
887 length_x = patch_ib(ib_patch_id)%length_x
888 length_y = patch_ib(ib_patch_id)%length_y
889 length_z = patch_ib(ib_patch_id)%length_z
890
891 center(1) = patch_ib(ib_patch_id)%x_centroid + real(gp%x_periodicity, wp)*(x_domain%end - x_domain%beg)
892 center(2) = patch_ib(ib_patch_id)%y_centroid + real(gp%y_periodicity, wp)*(y_domain%end - y_domain%beg)
893 center(3) = patch_ib(ib_patch_id)%z_centroid + real(gp%z_periodicity, wp)*(z_domain%end - z_domain%beg)
894
895 inverse_rotation(:,:) = patch_ib(ib_patch_id)%rotation_matrix_inverse(:,:)
896 rotation(:,:) = patch_ib(ib_patch_id)%rotation_matrix(:,:)
897
898 right = length_x/2
899 left = -length_x/2
900 top = length_y/2
901 bottom = -length_y/2
902 front = length_z/2
903 back = -length_z/2
904
905 xyz_local = [x_cc(i), y_cc(j), z_cc(k)] - center ! get coordinate frame centered on IB
906 xyz_local = matmul(inverse_rotation, xyz_local) ! rotate the frame into the IB's coordinate
907
908 dist_left = left - xyz_local(1)
909 dist_right = xyz_local(1) - right
910 dist_bottom = bottom - xyz_local(2)
911 dist_top = xyz_local(2) - top
912 dist_back = back - xyz_local(3)
913 dist_front = xyz_local(3) - front
914
915 min_dist = min(abs(dist_left), abs(dist_right), abs(dist_bottom), abs(dist_top), abs(dist_back), abs(dist_front))
916 dist_vec = 0._wp
917
918 if (f_approx_equal(min_dist, abs(dist_left))) then
919 gp%levelset = dist_left
920 if (.not. f_approx_equal(dist_left, 0._wp)) then
921 dist_vec(1) = dist_left/abs(dist_left)
922 end if
923 else if (f_approx_equal(min_dist, abs(dist_right))) then
924 gp%levelset = dist_right
925 if (.not. f_approx_equal(dist_right, 0._wp)) then
926 dist_vec(1) = -dist_right/abs(dist_right)
927 end if
928 else if (f_approx_equal(min_dist, abs(dist_bottom))) then
929 gp%levelset = dist_bottom
930 if (.not. f_approx_equal(dist_bottom, 0._wp)) then
931 dist_vec(2) = dist_bottom/abs(dist_bottom)
932 end if
933 else if (f_approx_equal(min_dist, abs(dist_top))) then
934 gp%levelset = dist_top
935 if (.not. f_approx_equal(dist_top, 0._wp)) then
936 dist_vec(2) = -dist_top/abs(dist_top)
937 end if
938 else if (f_approx_equal(min_dist, abs(dist_back))) then
939 gp%levelset = dist_back
940 if (.not. f_approx_equal(dist_back, 0._wp)) then
941 dist_vec(3) = dist_back/abs(dist_back)
942 end if
943 else if (f_approx_equal(min_dist, abs(dist_front))) then
944 gp%levelset = dist_front
945 if (.not. f_approx_equal(dist_front, 0._wp)) then
946 dist_vec(3) = -dist_front/abs(dist_front)
947 end if
948 end if
949
950 gp%levelset_norm = matmul(rotation, dist_vec)
951
952 end subroutine s_cuboid_levelset
953
954 !> Compute the signed distance and outward normal from a ghost point to a spherical immersed boundary
955 subroutine s_sphere_levelset(gp)
956
957
958# 494 "/home/runner/work/MFC/MFC/src/simulation/m_compute_levelset.fpp"
959#if MFC_OpenACC
960# 494 "/home/runner/work/MFC/MFC/src/simulation/m_compute_levelset.fpp"
961!$acc routine seq
962# 494 "/home/runner/work/MFC/MFC/src/simulation/m_compute_levelset.fpp"
963#elif MFC_OpenMP
964# 494 "/home/runner/work/MFC/MFC/src/simulation/m_compute_levelset.fpp"
965
966# 494 "/home/runner/work/MFC/MFC/src/simulation/m_compute_levelset.fpp"
967
968# 494 "/home/runner/work/MFC/MFC/src/simulation/m_compute_levelset.fpp"
969!$omp declare target device_type(any)
970# 494 "/home/runner/work/MFC/MFC/src/simulation/m_compute_levelset.fpp"
971#endif
972
973 type(ghost_point), intent(inout) :: gp
974 real(wp) :: radius, dist
975 real(wp), dimension(3) :: dist_vec, center, periodicity
976 integer :: i, j, k, ib_patch_id !< Loop index variables
977 ib_patch_id = gp%ib_patch_id
978 i = gp%loc(1)
979 j = gp%loc(2)
980 k = gp%loc(3)
981
982 radius = patch_ib(ib_patch_id)%radius
983 periodicity(1) = real(gp%x_periodicity, wp)*(x_domain%end - x_domain%beg)
984 periodicity(2) = real(gp%y_periodicity, wp)*(y_domain%end - y_domain%beg)
985 periodicity(3) = real(gp%z_periodicity, wp)*(z_domain%end - z_domain%beg)
986 center(1) = patch_ib(ib_patch_id)%x_centroid
987 center(2) = patch_ib(ib_patch_id)%y_centroid
988 center(3) = patch_ib(ib_patch_id)%z_centroid
989 center = center + periodicity
990
991 dist_vec(1) = x_cc(i) - center(1)
992 dist_vec(2) = y_cc(j) - center(2)
993 dist_vec(3) = z_cc(k) - center(3)
994 dist = sqrt(sum(dist_vec**2))
995 gp%levelset = dist - radius
996 if (f_approx_equal(dist, 0._wp)) then
997 gp%levelset_norm = (/1, 0, 0/)
998 else
999 gp%levelset_norm = dist_vec(:)/dist
1000 end if
1001
1002 end subroutine s_sphere_levelset
1003
1004 !> Compute the signed distance and outward normal from a ghost point to a cylindrical immersed boundary
1005 subroutine s_cylinder_levelset(gp)
1006
1007
1008# 530 "/home/runner/work/MFC/MFC/src/simulation/m_compute_levelset.fpp"
1009#if MFC_OpenACC
1010# 530 "/home/runner/work/MFC/MFC/src/simulation/m_compute_levelset.fpp"
1011!$acc routine seq
1012# 530 "/home/runner/work/MFC/MFC/src/simulation/m_compute_levelset.fpp"
1013#elif MFC_OpenMP
1014# 530 "/home/runner/work/MFC/MFC/src/simulation/m_compute_levelset.fpp"
1015
1016# 530 "/home/runner/work/MFC/MFC/src/simulation/m_compute_levelset.fpp"
1017
1018# 530 "/home/runner/work/MFC/MFC/src/simulation/m_compute_levelset.fpp"
1019!$omp declare target device_type(any)
1020# 530 "/home/runner/work/MFC/MFC/src/simulation/m_compute_levelset.fpp"
1021#endif
1022
1023 type(ghost_point), intent(inout) :: gp
1024 real(wp) :: radius
1025 real(wp), dimension(3) :: dist_sides_vec, dist_surface_vec, length
1026 real(wp), dimension(2) :: boundary
1027 real(wp) :: dist_side, dist_surface, side_pos
1028 integer :: i, j, k !< Loop index variables
1029 integer :: ib_patch_id !< patch ID
1030 real(wp), dimension(1:3) :: xyz_local, center !< x and y coordinates in local IB frame
1031 real(wp), dimension(1:3,1:3) :: rotation, inverse_rotation
1032
1033 ib_patch_id = gp%ib_patch_id
1034 i = gp%loc(1)
1035 j = gp%loc(2)
1036 k = gp%loc(3)
1037
1038 radius = patch_ib(ib_patch_id)%radius
1039 center(1) = patch_ib(ib_patch_id)%x_centroid + real(gp%x_periodicity, wp)*(x_domain%end - x_domain%beg)
1040 center(2) = patch_ib(ib_patch_id)%y_centroid + real(gp%y_periodicity, wp)*(y_domain%end - y_domain%beg)
1041 center(3) = patch_ib(ib_patch_id)%z_centroid + real(gp%z_periodicity, wp)*(z_domain%end - z_domain%beg)
1042 length(1) = patch_ib(ib_patch_id)%length_x
1043 length(2) = patch_ib(ib_patch_id)%length_y
1044 length(3) = patch_ib(ib_patch_id)%length_z
1045
1046 inverse_rotation(:,:) = patch_ib(ib_patch_id)%rotation_matrix_inverse(:,:)
1047 rotation(:,:) = patch_ib(ib_patch_id)%rotation_matrix(:,:)
1048
1049 if (.not. f_approx_equal(length(1), 0._wp)) then
1050 boundary(1) = -0.5_wp*length(1)
1051 boundary(2) = 0.5_wp*length(1)
1052 dist_sides_vec = (/1, 0, 0/)
1053 dist_surface_vec = (/0, 1, 1/)
1054 else if (.not. f_approx_equal(length(2), 0._wp)) then
1055 boundary(1) = -0.5_wp*length(2)
1056 boundary(2) = 0.5_wp*length(2)
1057 dist_sides_vec = (/0, 1, 0/)
1058 dist_surface_vec = (/1, 0, 1/)
1059 else if (.not. f_approx_equal(length(3), 0._wp)) then
1060 boundary(1) = -0.5_wp*length(3)
1061 boundary(2) = 0.5_wp*length(3)
1062 dist_sides_vec = (/0, 0, 1/)
1063 dist_surface_vec = (/1, 1, 0/)
1064 end if
1065
1066 xyz_local = [x_cc(i), y_cc(j), z_cc(k)] - center ! get coordinate frame centered on IB
1067 xyz_local = matmul(inverse_rotation, xyz_local) ! rotate the frame into the IB's coordinates
1068
1069 ! get distance to flat edge of cylinder
1070 side_pos = dot_product(xyz_local, dist_sides_vec)
1071 dist_side = min(abs(side_pos - boundary(1)), abs(boundary(2) - side_pos))
1072 ! get distance to curved side of cylinder
1073 dist_surface = norm2(xyz_local*dist_surface_vec) - radius
1074
1075 if (dist_side < abs(dist_surface)) then
1076 ! if the closest edge is flat
1077 gp%levelset = -dist_side
1078 if (f_approx_equal(dist_side, abs(side_pos - boundary(1)))) then
1079 gp%levelset_norm = matmul(rotation, -dist_sides_vec)
1080 else
1081 gp%levelset_norm = matmul(rotation, dist_sides_vec)
1082 end if
1083 else
1084 gp%levelset = dist_surface
1085 xyz_local = xyz_local*dist_surface_vec
1086 xyz_local = xyz_local/max(norm2(xyz_local), sgm_eps)
1087 gp%levelset_norm = matmul(rotation, xyz_local)
1088 end if
1089
1090 end subroutine s_cylinder_levelset
1091
1092 !> The STL patch is a 2/3D geometry that is imported from an STL file.
1093 subroutine s_model_levelset(gp)
1094
1095
1096# 604 "/home/runner/work/MFC/MFC/src/simulation/m_compute_levelset.fpp"
1097#if MFC_OpenACC
1098# 604 "/home/runner/work/MFC/MFC/src/simulation/m_compute_levelset.fpp"
1099!$acc routine seq
1100# 604 "/home/runner/work/MFC/MFC/src/simulation/m_compute_levelset.fpp"
1101#elif MFC_OpenMP
1102# 604 "/home/runner/work/MFC/MFC/src/simulation/m_compute_levelset.fpp"
1103
1104# 604 "/home/runner/work/MFC/MFC/src/simulation/m_compute_levelset.fpp"
1105
1106# 604 "/home/runner/work/MFC/MFC/src/simulation/m_compute_levelset.fpp"
1107!$omp declare target device_type(any)
1108# 604 "/home/runner/work/MFC/MFC/src/simulation/m_compute_levelset.fpp"
1109#endif
1110
1111 type(ghost_point), intent(inout) :: gp
1112 integer :: i, j, k, patch_id, boundary_edge_count
1113 real(wp), dimension(1:3) :: center, xyz_local
1114 real(wp) :: normals(1:3) !< Boundary normal buffer
1115 real(wp) :: distance
1116 real(wp), dimension(1:3,1:3) :: inverse_rotation, rotation
1117
1118 patch_id = gp%ib_patch_id
1119 i = gp%loc(1)
1120 j = gp%loc(2)
1121 k = gp%loc(3)
1122
1123 ! load in model values via the stl model index
1124 boundary_edge_count = gpu_boundary_edge_count(patch_ib(patch_id)%model_id)
1125
1126 center = 0._wp
1127 if (.not. f_is_default(patch_ib(patch_id)%x_centroid)) center(1) = patch_ib(patch_id)%x_centroid + real(gp%x_periodicity, &
1128 & wp)*(x_domain%end - x_domain%beg)
1129 if (.not. f_is_default(patch_ib(patch_id)%y_centroid)) center(2) = patch_ib(patch_id)%y_centroid + real(gp%y_periodicity, &
1130 & wp)*(y_domain%end - y_domain%beg)
1131 if (p > 0) then
1132 if (.not. f_is_default(patch_ib(patch_id)%z_centroid)) center(3) = patch_ib(patch_id)%z_centroid &
1133 & + real(gp%z_periodicity, wp)*(z_domain%end - z_domain%beg)
1134 end if
1135
1136 inverse_rotation(:,:) = patch_ib(patch_id)%rotation_matrix_inverse(:,:)
1137 rotation(:,:) = patch_ib(patch_id)%rotation_matrix(:,:)
1138
1139 ! determine where we are located in space
1140 xyz_local = (/x_cc(i) - center(1), y_cc(j) - center(2), 0._wp/)
1141 if (p > 0) then
1142 xyz_local(3) = z_cc(k) - center(3)
1143 end if
1144 xyz_local = matmul(inverse_rotation, xyz_local)
1145
1146 ! 3D models
1147 if (p > 0) then
1148 ! Get the boundary normals and shortest distance between the cell center and the model boundary
1149 call s_distance_normals_3d(gpu_ntrs(patch_ib(patch_id)%model_id), patch_ib(patch_id)%model_id, xyz_local, normals, &
1150 & distance)
1151
1152 ! Get the shortest distance between the cell center and the model boundary
1153 gp%levelset = distance
1154 gp%levelset = -abs(gp%levelset)
1155
1156 ! Assign the levelset_norm
1157 gp%levelset_norm = matmul(rotation, normals(1:3))
1158 else
1159 ! 2D models
1160 call s_distance_normals_2d(patch_ib(patch_id)%model_id, boundary_edge_count, xyz_local, normals, distance)
1161 gp%levelset = -abs(distance)
1162 gp%levelset_norm = matmul(rotation, normals(1:3))
1163 end if
1164
1165 end subroutine s_model_levelset
1166
1167end module m_compute_levelset
m_compute_levelset
Computes signed-distance level-set fields and surface normals for immersed-boundary patch geometries.
Definition m_compute_levelset.fpp.f90:323
m_compute_levelset::s_apply_levelset
impure subroutine, public s_apply_levelset(gps, num_gps)
Dispatch level-set distance and normal computations for all ghost points based on patch geometry type...
Definition m_compute_levelset.fpp.f90:340
m_compute_levelset::s_airfoil_levelset
subroutine s_airfoil_levelset(gp)
Compute the signed distance and outward normal from a ghost point to a 2D NACA airfoil surface.
Definition m_compute_levelset.fpp.f90:502
m_compute_levelset::s_sphere_levelset
subroutine s_sphere_levelset(gp)
Compute the signed distance and outward normal from a ghost point to a spherical immersed boundary.
Definition m_compute_levelset.fpp.f90:956
m_compute_levelset::s_circle_levelset
subroutine s_circle_levelset(gp)
Compute the signed distance and outward normal from a ghost point to a circular immersed boundary.
Definition m_compute_levelset.fpp.f90:459
m_compute_levelset::s_ellipse_levelset
subroutine s_ellipse_levelset(gp)
Compute the signed distance and outward normal from a ghost point to an elliptical immersed boundary.
Definition m_compute_levelset.fpp.f90:792
m_compute_levelset::s_cylinder_levelset
subroutine s_cylinder_levelset(gp)
Compute the signed distance and outward normal from a ghost point to a cylindrical immersed boundary.
Definition m_compute_levelset.fpp.f90:1006
m_compute_levelset::s_3d_airfoil_levelset
subroutine s_3d_airfoil_levelset(gp)
Compute the signed distance and outward normal from a ghost point to a 3D extruded airfoil surface.
Definition m_compute_levelset.fpp.f90:597
m_compute_levelset::s_rectangle_levelset
subroutine s_rectangle_levelset(gp)
Subroutine for computing the levelset values at a ghost point belonging to the rectangle IB.
Definition m_compute_levelset.fpp.f90:710
m_compute_levelset::s_model_levelset
subroutine s_model_levelset(gp)
The STL patch is a 2/3D geometry that is imported from an STL file.
Definition m_compute_levelset.fpp.f90:1094
m_compute_levelset::s_cuboid_levelset
subroutine s_cuboid_levelset(gp)
Compute the signed distance and outward normal from a ghost point to a cuboid immersed boundary.
Definition m_compute_levelset.fpp.f90:855
m_derived_types
Shared derived types for field data, patch geometry, bubble dynamics, and MPI I/O structures.
Definition m_derived_types.fpp.f90:323
m_global_parameters
Global parameters for the computational domain, fluid properties, and simulation algorithm configurat...
Definition m_global_parameters.fpp.f90:334
m_global_parameters::z_domain
type(bounds_info) z_domain
Definition m_global_parameters.fpp.f90:545
m_global_parameters::x_domain
type(bounds_info) x_domain
Definition m_global_parameters.fpp.f90:545
m_global_parameters::y_cc
real(wp), dimension(:), allocatable, target y_cc
Definition m_global_parameters.fpp.f90:399
m_global_parameters::y_domain
type(bounds_info) y_domain
Definition m_global_parameters.fpp.f90:545
m_global_parameters::z_cc
real(wp), dimension(:), allocatable, target z_cc
Definition m_global_parameters.fpp.f90:399
m_global_parameters::x_cc
real(wp), dimension(:), allocatable, target x_cc
Definition m_global_parameters.fpp.f90:399
m_global_parameters::ib_airfoil_grids
type(ib_airfoil_grid), dimension(num_ib_airfoils_max) ib_airfoil_grids
Per-airfoil computed surface grids.
Definition m_global_parameters.fpp.f90:738
m_helper_basic
Basic floating-point utilities: approximate equality, default detection, and coordinate bounds.
Definition m_helper_basic.fpp.f90:323
m_helper_basic::f_approx_equal
logical elemental function, public f_approx_equal(a, b, tol_input)
Check if two floating point numbers of wp are within tolerance.
Definition m_helper_basic.fpp.f90:340
m_helper_basic::f_is_default
logical elemental function, public f_is_default(var)
Checks if a real(wp) variable is of default value.
Definition m_helper_basic.fpp.f90:428
m_ib_patches
Allocate memory and read initial condition data for IC extrusion.
Definition m_ib_patches.fpp.f90:393
m_model
Binary STL file reader and processor for immersed boundary geometry.
Definition m_model.fpp.f90:324
m_model::s_distance_normals_2d
subroutine, public s_distance_normals_2d(pid, boundary_edge_count, point, normals, distance)
Determine the levelset distance and normals of 2D models by computing the exact closest point via pro...
Definition m_model.fpp.f90:1183
m_model::gpu_ntrs
integer, dimension(:), allocatable, public gpu_ntrs
GPU-friendly flat arrays for STL model data.
Definition m_model.fpp.f90:347
m_model::s_distance_normals_3d
subroutine, public s_distance_normals_3d(ntrs, pid, point, normals, distance)
Determine the levelset distance and normals of 3D models by computing the exact closest point via pro...
Definition m_model.fpp.f90:1039
m_model::gpu_boundary_edge_count
integer, dimension(:), allocatable, public gpu_boundary_edge_count
Definition m_model.fpp.f90:351
m_mpi_proxy
MPI halo exchange, domain decomposition, and buffer packing/unpacking for the simulation solver.
Definition m_mpi_proxy.fpp.f90:334
m_derived_types::ghost_point
Ghost Point for Immersed Boundaries.
Definition m_derived_types.fpp.f90:791