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