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