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.
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, &
487 & wp)*(glb_bounds(1)%end - glb_bounds(1)%beg))
488 dist_vec(2) = y_cc(j) - (patch_ib(ib_patch_id)%y_centroid + real(gp%y_periodicity, &
489 & wp)*(glb_bounds(2)%end - glb_bounds(2)%beg))
490 dist_vec(3) = 0._wp
491 dist = sqrt(sum(dist_vec**2))
492
493 gp%levelset = dist - radius
494 if (f_approx_equal(dist, 0._wp)) then
495 gp%levelset_norm = 0._wp
496 else
497 gp%levelset_norm = dist_vec(:)/dist
498 end if
499
500 end subroutine s_circle_levelset
501
502 !> Compute the signed distance and outward normal from a ghost point to a 2D NACA airfoil surface
503 subroutine s_airfoil_levelset(gp)
504
505
506# 112 "/home/runner/work/MFC/MFC/src/simulation/m_compute_levelset.fpp"
507#if MFC_OpenACC
508# 112 "/home/runner/work/MFC/MFC/src/simulation/m_compute_levelset.fpp"
509!$acc routine seq
510# 112 "/home/runner/work/MFC/MFC/src/simulation/m_compute_levelset.fpp"
511#elif MFC_OpenMP
512# 112 "/home/runner/work/MFC/MFC/src/simulation/m_compute_levelset.fpp"
513
514# 112 "/home/runner/work/MFC/MFC/src/simulation/m_compute_levelset.fpp"
515
516# 112 "/home/runner/work/MFC/MFC/src/simulation/m_compute_levelset.fpp"
517!$omp declare target device_type(any)
518# 112 "/home/runner/work/MFC/MFC/src/simulation/m_compute_levelset.fpp"
519#endif
520
521 type(ghost_point), intent(inout) :: gp
522 real(wp) :: dist, global_dist
523 integer :: global_id, airfoil_id, Np_local
524 real(wp), dimension(3) :: dist_vec
525 real(wp), dimension(1:3) :: xy_local, offset !< x and y coordinates in local IB frame
526 real(wp), dimension(1:2) :: center
527 real(wp), dimension(1:3,1:3) :: rotation, inverse_rotation
528 integer :: i, j, k, ib_patch_id !< Loop index variables
529 ib_patch_id = gp%ib_patch_id
530 i = gp%loc(1)
531 j = gp%loc(2)
532
533 airfoil_id = patch_ib(ib_patch_id)%airfoil_id
534 np_local = ib_airfoil_grids(airfoil_id)%Np
535 center(1) = patch_ib(ib_patch_id)%x_centroid + real(gp%x_periodicity, wp)*(glb_bounds(1)%end - glb_bounds(1)%beg)
536 center(2) = patch_ib(ib_patch_id)%y_centroid + real(gp%y_periodicity, wp)*(glb_bounds(2)%end - glb_bounds(2)%beg)
537
538 inverse_rotation(:,:) = patch_ib(ib_patch_id)%rotation_matrix_inverse(:,:)
539 rotation(:,:) = patch_ib(ib_patch_id)%rotation_matrix(:,:)
540 offset(:) = patch_ib(ib_patch_id)%centroid_offset(:)
541
542 xy_local = [x_cc(i) - center(1), y_cc(j) - center(2), 0._wp] ! get coordinate frame centered on IB
543 xy_local = matmul(inverse_rotation, xy_local) ! rotate the frame into the IB's coordinate
544 xy_local = xy_local - offset ! airfoils are a patch that require a centroid offset
545
546 if (xy_local(2) >= 0._wp) then
547 ! finds the location on the airfoil grid with the minimum distance (closest)
548 do k = 1, np_local
549 dist_vec(1) = ib_airfoil_grids(airfoil_id)%upper(k)%x - xy_local(1)
550 dist_vec(2) = ib_airfoil_grids(airfoil_id)%upper(k)%y - xy_local(2)
551 dist_vec(3) = 0._wp
552 dist = sqrt(sum(dist_vec**2))
553 if (k == 1) then
554 global_dist = dist
555 global_id = k
556 else
557 if (dist < global_dist) then
558 global_dist = dist
559 global_id = k
560 end if
561 end if
562 end do
563 dist_vec(1) = ib_airfoil_grids(airfoil_id)%upper(global_id)%x - xy_local(1)
564 dist_vec(2) = ib_airfoil_grids(airfoil_id)%upper(global_id)%y - xy_local(2)
565 dist_vec(3) = 0
566 dist = global_dist
567 else
568 do k = 1, np_local
569 dist_vec(1) = ib_airfoil_grids(airfoil_id)%lower(k)%x - xy_local(1)
570 dist_vec(2) = ib_airfoil_grids(airfoil_id)%lower(k)%y - xy_local(2)
571 dist_vec(3) = 0
572 dist = sqrt(sum(dist_vec**2))
573 if (k == 1) then
574 global_dist = dist
575 global_id = k
576 else
577 if (dist < global_dist) then
578 global_dist = dist
579 global_id = k
580 end if
581 end if
582 end do
583 dist_vec(1) = ib_airfoil_grids(airfoil_id)%lower(global_id)%x - xy_local(1)
584 dist_vec(2) = ib_airfoil_grids(airfoil_id)%lower(global_id)%y - xy_local(2)
585 dist_vec(3) = 0._wp
586 dist = global_dist
587 end if
588
589 gp%levelset = dist
590 if (f_approx_equal(dist, 0._wp)) then
591 gp%levelset_norm = 0._wp
592 else
593 gp%levelset_norm = matmul(rotation, dist_vec(:))/dist ! convert the normal vector back to global grid coordinates
594 end if
595
596 end subroutine s_airfoil_levelset
597
598 !> Compute the signed distance and outward normal from a ghost point to a 3D extruded airfoil surface
599 subroutine s_3d_airfoil_levelset(gp)
600
601
602# 194 "/home/runner/work/MFC/MFC/src/simulation/m_compute_levelset.fpp"
603#if MFC_OpenACC
604# 194 "/home/runner/work/MFC/MFC/src/simulation/m_compute_levelset.fpp"
605!$acc routine seq
606# 194 "/home/runner/work/MFC/MFC/src/simulation/m_compute_levelset.fpp"
607#elif MFC_OpenMP
608# 194 "/home/runner/work/MFC/MFC/src/simulation/m_compute_levelset.fpp"
609
610# 194 "/home/runner/work/MFC/MFC/src/simulation/m_compute_levelset.fpp"
611
612# 194 "/home/runner/work/MFC/MFC/src/simulation/m_compute_levelset.fpp"
613!$omp declare target device_type(any)
614# 194 "/home/runner/work/MFC/MFC/src/simulation/m_compute_levelset.fpp"
615#endif
616
617 type(ghost_point), intent(inout) :: gp
618 real(wp) :: dist_surf, dist_side, global_dist
619 integer :: global_id, airfoil_id, Np_local
620 real(wp) :: lz, z_max, z_min
621 real(wp), dimension(3) :: dist_vec
622 real(wp), dimension(1:3) :: xyz_local, center, offset, normal !< x, y, z coordinates in local IB frame
623 real(wp), dimension(1:3,1:3) :: rotation, inverse_rotation
624 integer :: i, j, k, l, ib_patch_id !< Loop index variables
625 ib_patch_id = gp%ib_patch_id
626 i = gp%loc(1)
627 j = gp%loc(2)
628 l = gp%loc(3)
629
630 airfoil_id = patch_ib(ib_patch_id)%airfoil_id
631 np_local = ib_airfoil_grids(airfoil_id)%Np
632 center(1) = patch_ib(ib_patch_id)%x_centroid + real(gp%x_periodicity, wp)*(glb_bounds(1)%end - glb_bounds(1)%beg)
633 center(2) = patch_ib(ib_patch_id)%y_centroid + real(gp%y_periodicity, wp)*(glb_bounds(2)%end - glb_bounds(2)%beg)
634 center(3) = patch_ib(ib_patch_id)%z_centroid + real(gp%z_periodicity, wp)*(glb_bounds(3)%end - glb_bounds(3)%beg)
635
636 lz = patch_ib(ib_patch_id)%length_z
637 inverse_rotation(:,:) = patch_ib(ib_patch_id)%rotation_matrix_inverse(:,:)
638 rotation(:,:) = patch_ib(ib_patch_id)%rotation_matrix(:,:)
639 offset(:) = patch_ib(ib_patch_id)%centroid_offset(:)
640
641 z_max = lz/2
642 z_min = -lz/2
643
644 xyz_local = [x_cc(i), y_cc(j), z_cc(l)] - center
645 xyz_local = matmul(inverse_rotation, xyz_local) ! rotate the frame into the IB's coordinates
646 xyz_local = xyz_local - offset ! airfoils are a patch that require a centroid offset
647
648 if (xyz_local(2) >= 0._wp) then
649 do k = 1, np_local
650 dist_vec(1) = xyz_local(1) - ib_airfoil_grids(airfoil_id)%upper(k)%x
651 dist_vec(2) = xyz_local(2) - ib_airfoil_grids(airfoil_id)%upper(k)%y
652 dist_vec(3) = 0._wp
653 dist_surf = sqrt(sum(dist_vec**2))
654 if (k == 1) then
655 global_dist = dist_surf
656 global_id = k
657 else
658 if (dist_surf < global_dist) then
659 global_dist = dist_surf
660 global_id = k
661 end if
662 end if
663 end do
664 dist_vec(1) = ib_airfoil_grids(airfoil_id)%upper(global_id)%x - xyz_local(1)
665 dist_vec(2) = ib_airfoil_grids(airfoil_id)%upper(global_id)%y - xyz_local(2)
666 dist_vec(3) = 0._wp
667 dist_surf = global_dist
668 else
669 do k = 1, np_local
670 dist_vec(1) = ib_airfoil_grids(airfoil_id)%lower(k)%x - xyz_local(1)
671 dist_vec(2) = ib_airfoil_grids(airfoil_id)%lower(k)%y - xyz_local(2)
672 dist_vec(3) = 0
673 dist_surf = sqrt(sum(dist_vec**2))
674 if (k == 1) then
675 global_dist = dist_surf
676 global_id = k
677 else
678 if (dist_surf < global_dist) then
679 global_dist = dist_surf
680 global_id = k
681 end if
682 end if
683 end do
684 dist_vec(1) = ib_airfoil_grids(airfoil_id)%lower(global_id)%x - xyz_local(1)
685 dist_vec(2) = ib_airfoil_grids(airfoil_id)%lower(global_id)%y - xyz_local(2)
686 dist_vec(3) = 0._wp
687 dist_surf = global_dist
688 end if
689
690 dist_side = min(abs(xyz_local(3) - z_min), abs(z_max - xyz_local(3)))
691
692 if (dist_side < dist_surf) then
693 gp%levelset = dist_side
694 normal = 0._wp
695 if (f_approx_equal(dist_side, abs(xyz_local(3) - z_min))) then
696 normal(3) = -1._wp
697 else
698 normal(3) = 1._wp
699 end if
700 gp%levelset_norm = matmul(rotation, normal)
701 else
702 gp%levelset = dist_surf
703 if (f_approx_equal(dist_surf, 0._wp)) then
704 gp%levelset_norm = 0._wp
705 else
706 gp%levelset_norm = matmul(rotation, dist_vec(:)/dist_surf)
707 end if
708 end if
709
710 end subroutine s_3d_airfoil_levelset
711
712 !> Compute the signed distance and outward normal from a ghost point to a 2D rectangle
713 subroutine s_rectangle_levelset(gp)
714
715
716# 294 "/home/runner/work/MFC/MFC/src/simulation/m_compute_levelset.fpp"
717#if MFC_OpenACC
718# 294 "/home/runner/work/MFC/MFC/src/simulation/m_compute_levelset.fpp"
719!$acc routine seq
720# 294 "/home/runner/work/MFC/MFC/src/simulation/m_compute_levelset.fpp"
721#elif MFC_OpenMP
722# 294 "/home/runner/work/MFC/MFC/src/simulation/m_compute_levelset.fpp"
723
724# 294 "/home/runner/work/MFC/MFC/src/simulation/m_compute_levelset.fpp"
725
726# 294 "/home/runner/work/MFC/MFC/src/simulation/m_compute_levelset.fpp"
727!$omp declare target device_type(any)
728# 294 "/home/runner/work/MFC/MFC/src/simulation/m_compute_levelset.fpp"
729#endif
730
731 type(ghost_point), intent(inout) :: gp
732 real(wp) :: top_right(2), bottom_left(2)
733 real(wp) :: min_dist
734 real(wp) :: side_dists(4)
735 real(wp) :: length_x, length_y
736 real(wp), dimension(1:3) :: xy_local, dist_vec !< x and y coordinates in local IB frame
737 real(wp), dimension(2) :: center !< x and y coordinates in local IB frame
738 real(wp), dimension(1:3,1:3) :: rotation, inverse_rotation
739 integer :: i, j, k !< Loop index variables
740 integer :: idx !< Shortest path direction indicator
741 integer :: ib_patch_id !< patch ID
742 ib_patch_id = gp%ib_patch_id
743 i = gp%loc(1)
744 j = gp%loc(2)
745
746 length_x = patch_ib(ib_patch_id)%length_x
747 length_y = patch_ib(ib_patch_id)%length_y
748 center(1) = patch_ib(ib_patch_id)%x_centroid + real(gp%x_periodicity, wp)*(glb_bounds(1)%end - glb_bounds(1)%beg)
749 center(2) = patch_ib(ib_patch_id)%y_centroid + real(gp%y_periodicity, wp)*(glb_bounds(2)%end - glb_bounds(2)%beg)
750 inverse_rotation(:,:) = patch_ib(ib_patch_id)%rotation_matrix_inverse(:,:)
751 rotation(:,:) = patch_ib(ib_patch_id)%rotation_matrix(:,:)
752
753 top_right(1) = length_x/2
754 top_right(2) = length_y/2
755 bottom_left(1) = -length_x/2
756 bottom_left(2) = -length_y/2
757
758 ! convert grid to local coordinates
759 xy_local = [x_cc(i) - center(1), y_cc(j) - center(2), 0._wp]
760 xy_local = matmul(inverse_rotation, xy_local)
761
762 side_dists(1) = bottom_left(1) - xy_local(1)
763 side_dists(2) = top_right(1) - xy_local(1)
764 side_dists(3) = bottom_left(2) - xy_local(2)
765 side_dists(4) = top_right(2) - xy_local(2)
766 min_dist = side_dists(1)
767 idx = 1
768
769 do k = 2, 4
770 if (abs(side_dists(k)) < abs(min_dist)) then
771 idx = k
772 min_dist = side_dists(idx)
773 end if
774 end do
775
776 gp%levelset = side_dists(idx)
777 dist_vec = 0._wp
778 if (.not. f_approx_equal(side_dists(idx), 0._wp)) then
779 if (idx == 1 .or. idx == 2) then
780 ! vector points along the x axis
781 dist_vec(1) = side_dists(idx)/abs(side_dists(idx))
782 else
783 ! vector points along the y axis
784 dist_vec(2) = side_dists(idx)/abs(side_dists(idx))
785 end if
786 ! convert the normal vector back into the global coordinate system
787 gp%levelset_norm = matmul(rotation, dist_vec)
788 else
789 gp%levelset_norm = 0._wp
790 end if
791
792 end subroutine s_rectangle_levelset
793
794 !> Compute the signed distance and outward normal from a ghost point to an elliptical immersed boundary
795 subroutine s_ellipse_levelset(gp)
796
797
798# 362 "/home/runner/work/MFC/MFC/src/simulation/m_compute_levelset.fpp"
799#if MFC_OpenACC
800# 362 "/home/runner/work/MFC/MFC/src/simulation/m_compute_levelset.fpp"
801!$acc routine seq
802# 362 "/home/runner/work/MFC/MFC/src/simulation/m_compute_levelset.fpp"
803#elif MFC_OpenMP
804# 362 "/home/runner/work/MFC/MFC/src/simulation/m_compute_levelset.fpp"
805
806# 362 "/home/runner/work/MFC/MFC/src/simulation/m_compute_levelset.fpp"
807
808# 362 "/home/runner/work/MFC/MFC/src/simulation/m_compute_levelset.fpp"
809!$omp declare target device_type(any)
810# 362 "/home/runner/work/MFC/MFC/src/simulation/m_compute_levelset.fpp"
811#endif
812
813 type(ghost_point), intent(inout) :: gp
814 real(wp) :: ellipse_coeffs(2) !< a and b in the ellipse equation
815 real(wp) :: quadratic_coeffs(3) !< A, B, C in the quadratic equation to compute levelset
816 real(wp) :: length_x, length_y
817 real(wp), dimension(1:3) :: xy_local, normal_vector !< x and y coordinates in local IB frame
818 real(wp), dimension(2) :: center !< x and y coordinates in local IB frame
819 real(wp), dimension(1:3,1:3) :: rotation, inverse_rotation
820 integer :: i, j !< Loop index variables
821 integer :: ib_patch_id !< patch ID
822 ib_patch_id = gp%ib_patch_id
823 i = gp%loc(1)
824 j = gp%loc(2)
825
826 length_x = patch_ib(ib_patch_id)%length_x
827 length_y = patch_ib(ib_patch_id)%length_y
828 center(1) = patch_ib(ib_patch_id)%x_centroid + real(gp%x_periodicity, wp)*(glb_bounds(1)%end - glb_bounds(1)%beg)
829 center(2) = patch_ib(ib_patch_id)%y_centroid + real(gp%y_periodicity, wp)*(glb_bounds(2)%end - glb_bounds(2)%beg)
830 inverse_rotation(:,:) = patch_ib(ib_patch_id)%rotation_matrix_inverse(:,:)
831 rotation(:,:) = patch_ib(ib_patch_id)%rotation_matrix(:,:)
832
833 ellipse_coeffs(1) = 0.5_wp*length_x
834 ellipse_coeffs(2) = 0.5_wp*length_y
835
836 xy_local = [x_cc(i) - center(1), y_cc(j) - center(2), 0._wp]
837 xy_local = matmul(inverse_rotation, xy_local)
838
839 normal_vector = xy_local
840 ! get the normal direction via the coordinate transformation method
841 normal_vector(2) = normal_vector(2)*(ellipse_coeffs(1)/ellipse_coeffs(2))**2._wp
842 normal_vector = normal_vector/sqrt(dot_product(normal_vector, normal_vector)) ! normalize the vector
843 gp%levelset_norm = matmul(rotation, normal_vector) ! save after rotating the vector to the global frame
844
845 ! use the normal vector to set up the quadratic equation for the levelset, using A, B, and C in indices 1, 2, and 3
846 quadratic_coeffs(1) = (normal_vector(1)/ellipse_coeffs(1))**2 + (normal_vector(2)/ellipse_coeffs(2))**2
847 quadratic_coeffs(2) = 2._wp*((xy_local(1)*normal_vector(1)/(ellipse_coeffs(1)**2)) + (xy_local(2)*normal_vector(2) &
848 & /(ellipse_coeffs(2)**2)))
849 quadratic_coeffs(3) = (xy_local(1)/ellipse_coeffs(1))**2._wp + (xy_local(2)/ellipse_coeffs(2))**2._wp - 1._wp
850
851 ! compute the levelset with the quadratic equation [ -B + sqrt(B^2 - 4AC) ] / 2A
852 gp%levelset = -0.5_wp*(-quadratic_coeffs(2) + sqrt(quadratic_coeffs(2)**2._wp - 4._wp*quadratic_coeffs(1) &
853 & *quadratic_coeffs(3)))/quadratic_coeffs(1)
854
855 end subroutine s_ellipse_levelset
856
857 !> Compute the signed distance and outward normal from a ghost point to a cuboid immersed boundary
858 subroutine s_cuboid_levelset(gp)
859
860
861# 411 "/home/runner/work/MFC/MFC/src/simulation/m_compute_levelset.fpp"
862#if MFC_OpenACC
863# 411 "/home/runner/work/MFC/MFC/src/simulation/m_compute_levelset.fpp"
864!$acc routine seq
865# 411 "/home/runner/work/MFC/MFC/src/simulation/m_compute_levelset.fpp"
866#elif MFC_OpenMP
867# 411 "/home/runner/work/MFC/MFC/src/simulation/m_compute_levelset.fpp"
868
869# 411 "/home/runner/work/MFC/MFC/src/simulation/m_compute_levelset.fpp"
870
871# 411 "/home/runner/work/MFC/MFC/src/simulation/m_compute_levelset.fpp"
872!$omp declare target device_type(any)
873# 411 "/home/runner/work/MFC/MFC/src/simulation/m_compute_levelset.fpp"
874#endif
875
876 type(ghost_point), intent(inout) :: gp
877 real(wp) :: Right, Left, Bottom, Top, Front, Back
878 real(wp) :: min_dist
879 real(wp) :: dist_left, dist_right, dist_bottom, dist_top, dist_back, dist_front
880 real(wp), dimension(3) :: center
881 real(wp) :: length_x, length_y, length_z
882 real(wp), dimension(1:3) :: xyz_local, dist_vec !< x and y coordinates in local IB frame
883 real(wp), dimension(1:3,1:3) :: rotation, inverse_rotation
884 integer :: i, j, k !< Loop index variables
885 integer :: ib_patch_id !< patch ID
886 ib_patch_id = gp%ib_patch_id
887 i = gp%loc(1)
888 j = gp%loc(2)
889 k = gp%loc(3)
890
891 length_x = patch_ib(ib_patch_id)%length_x
892 length_y = patch_ib(ib_patch_id)%length_y
893 length_z = patch_ib(ib_patch_id)%length_z
894
895 center(1) = patch_ib(ib_patch_id)%x_centroid + real(gp%x_periodicity, wp)*(glb_bounds(1)%end - glb_bounds(1)%beg)
896 center(2) = patch_ib(ib_patch_id)%y_centroid + real(gp%y_periodicity, wp)*(glb_bounds(2)%end - glb_bounds(2)%beg)
897 center(3) = patch_ib(ib_patch_id)%z_centroid + real(gp%z_periodicity, wp)*(glb_bounds(3)%end - glb_bounds(3)%beg)
898
899 inverse_rotation(:,:) = patch_ib(ib_patch_id)%rotation_matrix_inverse(:,:)
900 rotation(:,:) = patch_ib(ib_patch_id)%rotation_matrix(:,:)
901
902 right = length_x/2
903 left = -length_x/2
904 top = length_y/2
905 bottom = -length_y/2
906 front = length_z/2
907 back = -length_z/2
908
909 xyz_local = [x_cc(i), y_cc(j), z_cc(k)] - center ! get coordinate frame centered on IB
910 xyz_local = matmul(inverse_rotation, xyz_local) ! rotate the frame into the IB's coordinate
911
912 dist_left = left - xyz_local(1)
913 dist_right = xyz_local(1) - right
914 dist_bottom = bottom - xyz_local(2)
915 dist_top = xyz_local(2) - top
916 dist_back = back - xyz_local(3)
917 dist_front = xyz_local(3) - front
918
919 min_dist = min(abs(dist_left), abs(dist_right), abs(dist_bottom), abs(dist_top), abs(dist_back), abs(dist_front))
920 dist_vec = 0._wp
921
922 if (f_approx_equal(min_dist, abs(dist_left))) then
923 gp%levelset = dist_left
924 if (.not. f_approx_equal(dist_left, 0._wp)) then
925 dist_vec(1) = dist_left/abs(dist_left)
926 end if
927 else if (f_approx_equal(min_dist, abs(dist_right))) then
928 gp%levelset = dist_right
929 if (.not. f_approx_equal(dist_right, 0._wp)) then
930 dist_vec(1) = -dist_right/abs(dist_right)
931 end if
932 else if (f_approx_equal(min_dist, abs(dist_bottom))) then
933 gp%levelset = dist_bottom
934 if (.not. f_approx_equal(dist_bottom, 0._wp)) then
935 dist_vec(2) = dist_bottom/abs(dist_bottom)
936 end if
937 else if (f_approx_equal(min_dist, abs(dist_top))) then
938 gp%levelset = dist_top
939 if (.not. f_approx_equal(dist_top, 0._wp)) then
940 dist_vec(2) = -dist_top/abs(dist_top)
941 end if
942 else if (f_approx_equal(min_dist, abs(dist_back))) then
943 gp%levelset = dist_back
944 if (.not. f_approx_equal(dist_back, 0._wp)) then
945 dist_vec(3) = dist_back/abs(dist_back)
946 end if
947 else if (f_approx_equal(min_dist, abs(dist_front))) then
948 gp%levelset = dist_front
949 if (.not. f_approx_equal(dist_front, 0._wp)) then
950 dist_vec(3) = -dist_front/abs(dist_front)
951 end if
952 end if
953
954 gp%levelset_norm = matmul(rotation, dist_vec)
955
956 end subroutine s_cuboid_levelset
957
958 !> Compute the signed distance and outward normal from a ghost point to a spherical immersed boundary
959 subroutine s_sphere_levelset(gp)
960
961
962# 498 "/home/runner/work/MFC/MFC/src/simulation/m_compute_levelset.fpp"
963#if MFC_OpenACC
964# 498 "/home/runner/work/MFC/MFC/src/simulation/m_compute_levelset.fpp"
965!$acc routine seq
966# 498 "/home/runner/work/MFC/MFC/src/simulation/m_compute_levelset.fpp"
967#elif MFC_OpenMP
968# 498 "/home/runner/work/MFC/MFC/src/simulation/m_compute_levelset.fpp"
969
970# 498 "/home/runner/work/MFC/MFC/src/simulation/m_compute_levelset.fpp"
971
972# 498 "/home/runner/work/MFC/MFC/src/simulation/m_compute_levelset.fpp"
973!$omp declare target device_type(any)
974# 498 "/home/runner/work/MFC/MFC/src/simulation/m_compute_levelset.fpp"
975#endif
976
977 type(ghost_point), intent(inout) :: gp
978 real(wp) :: radius, dist
979 real(wp), dimension(3) :: dist_vec, center, periodicity
980 integer :: i, j, k, ib_patch_id !< Loop index variables
981 ib_patch_id = gp%ib_patch_id
982 i = gp%loc(1)
983 j = gp%loc(2)
984 k = gp%loc(3)
985
986 radius = patch_ib(ib_patch_id)%radius
987 periodicity(1) = real(gp%x_periodicity, wp)*(glb_bounds(1)%end - glb_bounds(1)%beg)
988 periodicity(2) = real(gp%y_periodicity, wp)*(glb_bounds(2)%end - glb_bounds(2)%beg)
989 periodicity(3) = real(gp%z_periodicity, wp)*(glb_bounds(3)%end - glb_bounds(3)%beg)
990 center(1) = patch_ib(ib_patch_id)%x_centroid
991 center(2) = patch_ib(ib_patch_id)%y_centroid
992 center(3) = patch_ib(ib_patch_id)%z_centroid
993 center = center + periodicity
994
995 dist_vec(1) = x_cc(i) - center(1)
996 dist_vec(2) = y_cc(j) - center(2)
997 dist_vec(3) = z_cc(k) - center(3)
998 dist = sqrt(sum(dist_vec**2))
999 gp%levelset = dist - radius
1000 if (f_approx_equal(dist, 0._wp)) then
1001 gp%levelset_norm = (/1, 0, 0/)
1002 else
1003 gp%levelset_norm = dist_vec(:)/dist
1004 end if
1005
1006 end subroutine s_sphere_levelset
1007
1008 !> Compute the signed distance and outward normal from a ghost point to a cylindrical immersed boundary
1009 subroutine s_cylinder_levelset(gp)
1010
1011
1012# 534 "/home/runner/work/MFC/MFC/src/simulation/m_compute_levelset.fpp"
1013#if MFC_OpenACC
1014# 534 "/home/runner/work/MFC/MFC/src/simulation/m_compute_levelset.fpp"
1015!$acc routine seq
1016# 534 "/home/runner/work/MFC/MFC/src/simulation/m_compute_levelset.fpp"
1017#elif MFC_OpenMP
1018# 534 "/home/runner/work/MFC/MFC/src/simulation/m_compute_levelset.fpp"
1019
1020# 534 "/home/runner/work/MFC/MFC/src/simulation/m_compute_levelset.fpp"
1021
1022# 534 "/home/runner/work/MFC/MFC/src/simulation/m_compute_levelset.fpp"
1023!$omp declare target device_type(any)
1024# 534 "/home/runner/work/MFC/MFC/src/simulation/m_compute_levelset.fpp"
1025#endif
1026
1027 type(ghost_point), intent(inout) :: gp
1028 real(wp) :: radius
1029 real(wp), dimension(3) :: dist_sides_vec, dist_surface_vec, length
1030 real(wp), dimension(2) :: boundary
1031 real(wp) :: dist_side, dist_surface, side_pos
1032 integer :: i, j, k !< Loop index variables
1033 integer :: ib_patch_id !< patch ID
1034 real(wp), dimension(1:3) :: xyz_local, center !< x and y coordinates in local IB frame
1035 real(wp), dimension(1:3,1:3) :: rotation, inverse_rotation
1036
1037 ib_patch_id = gp%ib_patch_id
1038 i = gp%loc(1)
1039 j = gp%loc(2)
1040 k = gp%loc(3)
1041
1042 radius = patch_ib(ib_patch_id)%radius
1043 center(1) = patch_ib(ib_patch_id)%x_centroid + real(gp%x_periodicity, wp)*(glb_bounds(1)%end - glb_bounds(1)%beg)
1044 center(2) = patch_ib(ib_patch_id)%y_centroid + real(gp%y_periodicity, wp)*(glb_bounds(2)%end - glb_bounds(2)%beg)
1045 center(3) = patch_ib(ib_patch_id)%z_centroid + real(gp%z_periodicity, wp)*(glb_bounds(3)%end - glb_bounds(3)%beg)
1046 length(1) = patch_ib(ib_patch_id)%length_x
1047 length(2) = patch_ib(ib_patch_id)%length_y
1048 length(3) = patch_ib(ib_patch_id)%length_z
1049
1050 inverse_rotation(:,:) = patch_ib(ib_patch_id)%rotation_matrix_inverse(:,:)
1051 rotation(:,:) = patch_ib(ib_patch_id)%rotation_matrix(:,:)
1052
1053 if (.not. f_approx_equal(length(1), 0._wp)) then
1054 boundary(1) = -0.5_wp*length(1)
1055 boundary(2) = 0.5_wp*length(1)
1056 dist_sides_vec = (/1, 0, 0/)
1057 dist_surface_vec = (/0, 1, 1/)
1058 else if (.not. f_approx_equal(length(2), 0._wp)) then
1059 boundary(1) = -0.5_wp*length(2)
1060 boundary(2) = 0.5_wp*length(2)
1061 dist_sides_vec = (/0, 1, 0/)
1062 dist_surface_vec = (/1, 0, 1/)
1063 else if (.not. f_approx_equal(length(3), 0._wp)) then
1064 boundary(1) = -0.5_wp*length(3)
1065 boundary(2) = 0.5_wp*length(3)
1066 dist_sides_vec = (/0, 0, 1/)
1067 dist_surface_vec = (/1, 1, 0/)
1068 end if
1069
1070 xyz_local = [x_cc(i), y_cc(j), z_cc(k)] - center ! get coordinate frame centered on IB
1071 xyz_local = matmul(inverse_rotation, xyz_local) ! rotate the frame into the IB's coordinates
1072
1073 ! get distance to flat edge of cylinder
1074 side_pos = dot_product(xyz_local, dist_sides_vec)
1075 dist_side = min(abs(side_pos - boundary(1)), abs(boundary(2) - side_pos))
1076 ! get distance to curved side of cylinder
1077 dist_surface = norm2(xyz_local*dist_surface_vec) - radius
1078
1079 if (dist_side < abs(dist_surface)) then
1080 ! if the closest edge is flat
1081 gp%levelset = -dist_side
1082 if (f_approx_equal(dist_side, abs(side_pos - boundary(1)))) then
1083 gp%levelset_norm = matmul(rotation, -dist_sides_vec)
1084 else
1085 gp%levelset_norm = matmul(rotation, dist_sides_vec)
1086 end if
1087 else
1088 gp%levelset = dist_surface
1089 xyz_local = xyz_local*dist_surface_vec
1090 xyz_local = xyz_local/max(norm2(xyz_local), sgm_eps)
1091 gp%levelset_norm = matmul(rotation, xyz_local)
1092 end if
1093
1094 end subroutine s_cylinder_levelset
1095
1096 !> The STL patch is a 2/3D geometry that is imported from an STL file.
1097 subroutine s_model_levelset(gp)
1098
1099
1100# 608 "/home/runner/work/MFC/MFC/src/simulation/m_compute_levelset.fpp"
1101#if MFC_OpenACC
1102# 608 "/home/runner/work/MFC/MFC/src/simulation/m_compute_levelset.fpp"
1103!$acc routine seq
1104# 608 "/home/runner/work/MFC/MFC/src/simulation/m_compute_levelset.fpp"
1105#elif MFC_OpenMP
1106# 608 "/home/runner/work/MFC/MFC/src/simulation/m_compute_levelset.fpp"
1107
1108# 608 "/home/runner/work/MFC/MFC/src/simulation/m_compute_levelset.fpp"
1109
1110# 608 "/home/runner/work/MFC/MFC/src/simulation/m_compute_levelset.fpp"
1111!$omp declare target device_type(any)
1112# 608 "/home/runner/work/MFC/MFC/src/simulation/m_compute_levelset.fpp"
1113#endif
1114
1115 type(ghost_point), intent(inout) :: gp
1116 integer :: i, j, k, patch_id, boundary_edge_count
1117 real(wp), dimension(1:3) :: center, xyz_local
1118 real(wp) :: normals(1:3) !< Boundary normal buffer
1119 real(wp) :: distance
1120 real(wp), dimension(1:3,1:3) :: inverse_rotation, rotation
1121
1122 patch_id = gp%ib_patch_id
1123 i = gp%loc(1)
1124 j = gp%loc(2)
1125 k = gp%loc(3)
1126
1127 ! load in model values via the stl model index
1128 boundary_edge_count = gpu_boundary_edge_count(patch_ib(patch_id)%model_id)
1129
1130 center = 0._wp
1131 if (.not. f_is_default(patch_ib(patch_id)%x_centroid)) center(1) = patch_ib(patch_id)%x_centroid + real(gp%x_periodicity, &
1132 & wp)*(glb_bounds(1)%end - glb_bounds(1)%beg)
1133 if (.not. f_is_default(patch_ib(patch_id)%y_centroid)) center(2) = patch_ib(patch_id)%y_centroid + real(gp%y_periodicity, &
1134 & wp)*(glb_bounds(2)%end - glb_bounds(2)%beg)
1135 if (p > 0) then
1136 if (.not. f_is_default(patch_ib(patch_id)%z_centroid)) center(3) = patch_ib(patch_id)%z_centroid &
1137 & + real(gp%z_periodicity, wp)*(glb_bounds(3)%end - glb_bounds(3)%beg)
1138 end if
1139
1140 inverse_rotation(:,:) = patch_ib(patch_id)%rotation_matrix_inverse(:,:)
1141 rotation(:,:) = patch_ib(patch_id)%rotation_matrix(:,:)
1142
1143 ! determine where we are located in space
1144 xyz_local = (/x_cc(i) - center(1), y_cc(j) - center(2), 0._wp/)
1145 if (p > 0) then
1146 xyz_local(3) = z_cc(k) - center(3)
1147 end if
1148 xyz_local = matmul(inverse_rotation, xyz_local)
1149
1150 ! 3D models
1151 if (p > 0) then
1152 ! Get the boundary normals and shortest distance between the cell center and the model boundary
1153 call s_distance_normals_3d(gpu_ntrs(patch_ib(patch_id)%model_id), patch_ib(patch_id)%model_id, xyz_local, normals, &
1154 & distance)
1155
1156 ! Get the shortest distance between the cell center and the model boundary
1157 gp%levelset = distance
1158 gp%levelset = -abs(gp%levelset)
1159
1160 ! Assign the levelset_norm
1161 gp%levelset_norm = matmul(rotation, normals(1:3))
1162 else
1163 ! 2D models
1164 call s_distance_normals_2d(patch_ib(patch_id)%model_id, boundary_edge_count, xyz_local, normals, distance)
1165 gp%levelset = -abs(distance)
1166 gp%levelset_norm = matmul(rotation, normals(1:3))
1167 end if
1168
1169 end subroutine s_model_levelset
1170
1171end 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:504
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:960
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:796
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:1010
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:600
m_compute_levelset::s_rectangle_levelset
subroutine s_rectangle_levelset(gp)
Compute the signed distance and outward normal from a ghost point to a 2D rectangle.
Definition m_compute_levelset.fpp.f90:714
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:1098
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:859
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::glb_bounds
type(bounds_info), dimension(3) glb_bounds
Definition m_global_parameters.fpp.f90:395
m_global_parameters::y_cc
real(wp), dimension(:), allocatable, target y_cc
Definition m_global_parameters.fpp.f90:400
m_global_parameters::z_cc
real(wp), dimension(:), allocatable, target z_cc
Definition m_global_parameters.fpp.f90:400
m_global_parameters::x_cc
real(wp), dimension(:), allocatable, target x_cc
Definition m_global_parameters.fpp.f90:400
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:826
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)
Check 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:802