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