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