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