MFC
Exascale flow solver
Loading...
Searching...
No Matches
m_model.fpp.f90
Go to the documentation of this file.
1# 1 "/home/runner/work/MFC/MFC/src/common/m_model.fpp"
2!>
3!! @file
4!! @author Henry Le Berre <hberre3@gatech.edu>
5!! @brief Contains module m_model
6
7# 1 "/home/runner/work/MFC/MFC/src/common/include/macros.fpp" 1
8# 1 "/home/runner/work/MFC/MFC/src/common/include/parallel_macros.fpp" 1
9# 1 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp" 1
10# 2 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
11# 3 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
12# 4 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
13# 5 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
14# 6 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
15
16# 8 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
17# 9 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
18# 10 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
19
20# 17 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
21
22# 46 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
23
24# 58 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
25
26# 68 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
27
28# 98 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
29
30# 110 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
31
32# 120 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
33! New line at end of file is required for FYPP
34# 2 "/home/runner/work/MFC/MFC/src/common/include/parallel_macros.fpp" 2
35# 1 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp" 1
36# 1 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp" 1
37# 2 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
38# 3 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
39# 4 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
40# 5 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
41# 6 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
42
43# 8 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
44# 9 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
45# 10 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
46
47# 17 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
48
49# 46 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
50
51# 58 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
52
53# 68 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
54
55# 98 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
56
57# 110 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
58
59# 120 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
60! New line at end of file is required for FYPP
61# 2 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp" 2
62
63# 4 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
64# 5 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
65# 6 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
66# 7 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
67# 8 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
68
69# 20 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
70
71# 43 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
72
73# 48 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
74
75# 53 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
76
77# 58 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
78
79# 63 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
80
81# 68 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
82
83# 76 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
84
85# 81 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
86
87# 86 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
88
89# 91 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
90
91# 96 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
92
93# 101 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
94
95# 106 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
96
97# 111 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
98
99# 116 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
100
101# 121 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
102
103# 151 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
104
105# 192 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
106
107# 206 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
108
109# 231 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
110
111# 242 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
112
113# 244 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
114# 255 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
115
116# 284 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
117
118# 294 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
119
120# 304 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
121
122# 313 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
123
124# 330 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
125
126# 340 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
127
128# 347 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
129
130# 353 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
131
132# 359 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
133
134# 365 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
135
136# 371 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
137
138# 377 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
139! New line at end of file is required for FYPP
140# 3 "/home/runner/work/MFC/MFC/src/common/include/parallel_macros.fpp" 2
141# 1 "/home/runner/work/MFC/MFC/src/common/include/acc_macros.fpp" 1
142# 1 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp" 1
143# 2 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
144# 3 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
145# 4 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
146# 5 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
147# 6 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
148
149# 8 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
150# 9 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
151# 10 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
152
153# 17 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
154
155# 46 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
156
157# 58 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
158
159# 68 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
160
161# 98 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
162
163# 110 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
164
165# 120 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
166! New line at end of file is required for FYPP
167# 2 "/home/runner/work/MFC/MFC/src/common/include/acc_macros.fpp" 2
168
169# 7 "/home/runner/work/MFC/MFC/src/common/include/acc_macros.fpp"
170
171# 17 "/home/runner/work/MFC/MFC/src/common/include/acc_macros.fpp"
172
173# 22 "/home/runner/work/MFC/MFC/src/common/include/acc_macros.fpp"
174
175# 27 "/home/runner/work/MFC/MFC/src/common/include/acc_macros.fpp"
176
177# 32 "/home/runner/work/MFC/MFC/src/common/include/acc_macros.fpp"
178
179# 37 "/home/runner/work/MFC/MFC/src/common/include/acc_macros.fpp"
180
181# 42 "/home/runner/work/MFC/MFC/src/common/include/acc_macros.fpp"
182
183# 47 "/home/runner/work/MFC/MFC/src/common/include/acc_macros.fpp"
184
185# 52 "/home/runner/work/MFC/MFC/src/common/include/acc_macros.fpp"
186
187# 57 "/home/runner/work/MFC/MFC/src/common/include/acc_macros.fpp"
188
189# 62 "/home/runner/work/MFC/MFC/src/common/include/acc_macros.fpp"
190
191# 73 "/home/runner/work/MFC/MFC/src/common/include/acc_macros.fpp"
192
193# 78 "/home/runner/work/MFC/MFC/src/common/include/acc_macros.fpp"
194
195# 83 "/home/runner/work/MFC/MFC/src/common/include/acc_macros.fpp"
196
197# 88 "/home/runner/work/MFC/MFC/src/common/include/acc_macros.fpp"
198
199# 103 "/home/runner/work/MFC/MFC/src/common/include/acc_macros.fpp"
200
201# 131 "/home/runner/work/MFC/MFC/src/common/include/acc_macros.fpp"
202
203# 160 "/home/runner/work/MFC/MFC/src/common/include/acc_macros.fpp"
204
205# 175 "/home/runner/work/MFC/MFC/src/common/include/acc_macros.fpp"
206
207# 193 "/home/runner/work/MFC/MFC/src/common/include/acc_macros.fpp"
208
209# 215 "/home/runner/work/MFC/MFC/src/common/include/acc_macros.fpp"
210
211# 244 "/home/runner/work/MFC/MFC/src/common/include/acc_macros.fpp"
212
213# 259 "/home/runner/work/MFC/MFC/src/common/include/acc_macros.fpp"
214
215# 269 "/home/runner/work/MFC/MFC/src/common/include/acc_macros.fpp"
216
217# 278 "/home/runner/work/MFC/MFC/src/common/include/acc_macros.fpp"
218
219# 294 "/home/runner/work/MFC/MFC/src/common/include/acc_macros.fpp"
220
221# 304 "/home/runner/work/MFC/MFC/src/common/include/acc_macros.fpp"
222
223# 311 "/home/runner/work/MFC/MFC/src/common/include/acc_macros.fpp"
224! New line at end of file is required for FYPP
225# 4 "/home/runner/work/MFC/MFC/src/common/include/parallel_macros.fpp" 2
226
227! GPU parallel region (scalar reductions, maxval/minval)
228# 23 "/home/runner/work/MFC/MFC/src/common/include/parallel_macros.fpp"
229
230! GPU parallel loop over threads (most common GPU macro)
231# 43 "/home/runner/work/MFC/MFC/src/common/include/parallel_macros.fpp"
232
233! Required closing for GPU_PARALLEL_LOOP
234# 55 "/home/runner/work/MFC/MFC/src/common/include/parallel_macros.fpp"
235
236! Mark routine for device compilation
237# 112 "/home/runner/work/MFC/MFC/src/common/include/parallel_macros.fpp"
238
239! Declare device-resident data
240# 130 "/home/runner/work/MFC/MFC/src/common/include/parallel_macros.fpp"
241
242! Inner loop within a GPU parallel region
243# 145 "/home/runner/work/MFC/MFC/src/common/include/parallel_macros.fpp"
244
245! Scoped GPU data region
246# 164 "/home/runner/work/MFC/MFC/src/common/include/parallel_macros.fpp"
247
248! Host code with device pointers (for MPI with GPU buffers)
249# 193 "/home/runner/work/MFC/MFC/src/common/include/parallel_macros.fpp"
250
251! Allocate device memory (unscoped)
252# 207 "/home/runner/work/MFC/MFC/src/common/include/parallel_macros.fpp"
253
254! Free device memory
255# 219 "/home/runner/work/MFC/MFC/src/common/include/parallel_macros.fpp"
256
257! Atomic operation on device
258# 231 "/home/runner/work/MFC/MFC/src/common/include/parallel_macros.fpp"
259
260! End atomic capture block
261# 242 "/home/runner/work/MFC/MFC/src/common/include/parallel_macros.fpp"
262
263! Copy data between host and device
264# 254 "/home/runner/work/MFC/MFC/src/common/include/parallel_macros.fpp"
265
266! Synchronization barrier
267# 266 "/home/runner/work/MFC/MFC/src/common/include/parallel_macros.fpp"
268
269! Import GPU library module (openacc or omp_lib)
270# 275 "/home/runner/work/MFC/MFC/src/common/include/parallel_macros.fpp"
271
272! Emit code only for AMD compiler
273# 282 "/home/runner/work/MFC/MFC/src/common/include/parallel_macros.fpp"
274
275! Emit code for non-Cray compilers
276# 289 "/home/runner/work/MFC/MFC/src/common/include/parallel_macros.fpp"
277
278! Emit code only for Cray compiler
279# 296 "/home/runner/work/MFC/MFC/src/common/include/parallel_macros.fpp"
280
281! Emit code for non-NVIDIA compilers
282# 303 "/home/runner/work/MFC/MFC/src/common/include/parallel_macros.fpp"
283
284# 305 "/home/runner/work/MFC/MFC/src/common/include/parallel_macros.fpp"
285# 306 "/home/runner/work/MFC/MFC/src/common/include/parallel_macros.fpp"
286! New line at end of file is required for FYPP
287# 2 "/home/runner/work/MFC/MFC/src/common/include/macros.fpp" 2
288
289# 14 "/home/runner/work/MFC/MFC/src/common/include/macros.fpp"
290
291! Caution: This macro requires the use of a binding script to set CUDA_VISIBLE_DEVICES, such that we have one GPU device per MPI
292! rank. That's because for both cudaMemAdvise (preferred location) and cudaMemPrefetchAsync we use location = device_id = 0. For an
293! example see misc/nvidia_uvm/bind.sh. NVIDIA unified memory page placement hint
294# 57 "/home/runner/work/MFC/MFC/src/common/include/macros.fpp"
295
296! Allocate and create GPU device memory
297# 77 "/home/runner/work/MFC/MFC/src/common/include/macros.fpp"
298
299! Free GPU device memory and deallocate
300# 85 "/home/runner/work/MFC/MFC/src/common/include/macros.fpp"
301
302! Cray-specific GPU pointer setup for vector fields
303# 109 "/home/runner/work/MFC/MFC/src/common/include/macros.fpp"
304
305! Cray-specific GPU pointer setup for scalar fields
306# 125 "/home/runner/work/MFC/MFC/src/common/include/macros.fpp"
307
308! Cray-specific GPU pointer setup for acoustic source spatials
309# 150 "/home/runner/work/MFC/MFC/src/common/include/macros.fpp"
310
311# 156 "/home/runner/work/MFC/MFC/src/common/include/macros.fpp"
312
313# 163 "/home/runner/work/MFC/MFC/src/common/include/macros.fpp"
314! New line at end of file is required for FYPP
315# 7 "/home/runner/work/MFC/MFC/src/common/m_model.fpp" 2
316
317!> @brief Binary STL file reader and processor for immersed boundary geometry
318module m_model
319
320 use m_helper
321 use m_mpi_proxy
322 use m_derived_types
323 use iso_c_binding, only: c_char, c_int32_t, c_int16_t, c_float
324
325 implicit none
326
327 private
328
329 public :: f_model_read, s_model_write, s_model_free, f_model_is_inside, models, gpu_ntrs, gpu_trs_v, gpu_trs_n, &
330 & gpu_boundary_v, gpu_boundary_edge_count, gpu_total_vertices, stl_bounding_boxes
331
332 ! Subroutines for STL immersed boundaries
333 public :: s_check_boundary, s_register_edge, f_model_is_inside_flat, s_distance_normals_3d, s_distance_normals_2d, &
334 & s_pack_model_for_gpu
335
336#ifdef MFC_SIMULATION
337 public :: s_instantiate_stl_models
338#endif
339
340 type(t_model_array), allocatable, target :: models(:) !< STL/OBJ models for IB markers and levelset
341 integer, allocatable :: gpu_ntrs(:) !< GPU-friendly flat arrays for STL model data
342 real(wp), allocatable, dimension(:,:,:,:) :: gpu_trs_v
343 real(wp), allocatable, dimension(:,:,:) :: gpu_trs_n
344 real(wp), allocatable, dimension(:,:,:,:) :: gpu_boundary_v
345 integer, allocatable :: gpu_boundary_edge_count(:)
346 integer, allocatable :: gpu_total_vertices(:)
347 real(wp), allocatable :: stl_bounding_boxes(:,:,:)
348
349# 39 "/home/runner/work/MFC/MFC/src/common/m_model.fpp"
350#if defined(MFC_OpenACC)
351# 39 "/home/runner/work/MFC/MFC/src/common/m_model.fpp"
352!$acc declare create(gpu_ntrs, gpu_trs_v, gpu_trs_n, gpu_boundary_v, gpu_boundary_edge_count, gpu_total_vertices)
353# 39 "/home/runner/work/MFC/MFC/src/common/m_model.fpp"
354#elif defined(MFC_OpenMP)
355# 39 "/home/runner/work/MFC/MFC/src/common/m_model.fpp"
356!$omp declare target (gpu_ntrs, gpu_trs_v, gpu_trs_n, gpu_boundary_v, gpu_boundary_edge_count, gpu_total_vertices)
357# 39 "/home/runner/work/MFC/MFC/src/common/m_model.fpp"
358#endif
359
360contains
361
362 !> Read a binary STL file.
363 impure subroutine s_read_stl_binary(filepath, model)
364
365 character(LEN=*), intent(in) :: filepath
366 type(t_model), intent(out) :: model
367 integer :: i, iunit, iostat
368 character(kind=c_char, len=80) :: header
369 integer(kind=c_int32_t) :: ntriangles
370 real(kind=c_float) :: normal(3), v(3, 3), v_norm
371 integer(kind=c_int16_t) :: attribute
372
373 open (newunit=iunit, file=filepath, action='READ', form='UNFORMATTED', status='OLD', iostat=iostat, access='STREAM')
374
375 if (iostat /= 0) then
376 print *, "Error: could not open Binary STL file ", filepath
377
378 call s_mpi_abort()
379 end if
380
381 read (iunit, iostat=iostat) header, ntriangles
382
383 if (iostat /= 0) then
384 print *, "Error: could not read header from Binary STL file ", filepath
385
386 call s_mpi_abort()
387 end if
388
389 model%ntrs = ntriangles
390
391 allocate (model%trs(model%ntrs))
392
393 do i = 1, model%ntrs
394 read (iunit) normal(:), v(1,:), v(2,:), v(3,:), attribute
395
396 model%trs(i)%v = v
397 model%trs(i)%n = normal
398 v_norm = sqrt(normal(1)**2 + normal(2)**2 + normal(3)**2)
399 if (v_norm > 0._wp) model%trs(i)%n = normal/v_norm
400 end do
401
402 close (iunit)
403
404 end subroutine s_read_stl_binary
405
406 !> Read an ASCII STL file.
407 impure subroutine s_read_stl_ascii(filepath, model)
408
409 character(LEN=*), intent(in) :: filepath
410 type(t_model), intent(out) :: model
411 integer :: i, j, iunit, iostat
412 character(80) :: line, buffered_line
413 logical :: is_buffered
414 real(wp) :: normal(3), v_norm
415
416 is_buffered = .false.
417
418 open (newunit=iunit, file=filepath, action='READ', form='FORMATTED', status='OLD', iostat=iostat, access='STREAM')
419
420 if (iostat /= 0) then
421 print *, "Error: could not open ASCII STL file ", filepath
422 call s_mpi_abort()
423 end if
424
425 model%ntrs = 0
426 do
427 if (is_buffered) then
428 line = buffered_line
429 is_buffered = .false.
430 else
431 if (.not. f_read_line(iunit, line)) exit
432 end if
433
434 if (line(1:6) == "facet ") then
435 model%ntrs = model%ntrs + 1
436 end if
437 end do
438
439 allocate (model%trs(model%ntrs))
440
441 rewind(iunit)
442
443 i = 1
444 do
445 if (is_buffered) then
446 line = buffered_line
447 is_buffered = .false.
448 else
449 if (.not. f_read_line(iunit, line)) exit
450 end if
451
452 if (line(1:5) == "solid") cycle
453 if (line(1:8) == "endsolid") exit
454
455 if (line(1:12) /= "facet normal") then
456 print *, "Error: expected facet normal in STL file ", filepath
457 call s_mpi_abort()
458 end if
459
460 call s_skip_ignored_lines(iunit, buffered_line, is_buffered)
461 read (line(13:), *) normal
462 v_norm = sqrt(normal(1)**2 + normal(2)**2 + normal(3)**2)
463 if (v_norm > 0._wp) model%trs(i)%n = normal/v_norm
464
465 call s_skip_ignored_lines(iunit, buffered_line, is_buffered)
466 if (is_buffered) then
467 line = buffered_line
468 is_buffered = .false.
469 else
470 read (iunit, '(A)') line
471 end if
472
473 do j = 1, 3
474 if (is_buffered) then
475 line = buffered_line
476 is_buffered = .false.
477 else
478 if (.not. f_read_line(iunit, line)) exit
479 end if
480
481 if (line(1:6) /= "vertex") then
482 print *, "Error: expected vertex in STL file ", filepath
483 call s_mpi_abort()
484 end if
485
486 call s_skip_ignored_lines(iunit, buffered_line, is_buffered)
487 read (line(7:), *) model%trs(i)%v(j,:)
488 end do
489
490 if (is_buffered) then
491 line = buffered_line
492 is_buffered = .false.
493 else
494 if (.not. f_read_line(iunit, line)) exit
495 end if
496
497 if (is_buffered) then
498 line = buffered_line
499 is_buffered = .false.
500 else
501 if (.not. f_read_line(iunit, line)) exit
502 end if
503
504 if (line(1:8) /= "endfacet") then
505 print *, "Error: expected endfacet in STL file ", filepath
506 call s_mpi_abort()
507 end if
508
509 i = i + 1
510 end do
511
512 end subroutine s_read_stl_ascii
513
514 !> Read an STL file.
515 impure subroutine s_read_stl(filepath, model)
516
517 character(LEN=*), intent(in) :: filepath
518 type(t_model), intent(out) :: model
519 integer :: iunit, iostat
520 character(80) :: line
521
522 open (newunit=iunit, file=filepath, action='READ', form='FORMATTED', status='OLD', iostat=iostat, access='STREAM')
523
524 if (iostat /= 0) then
525 print *, "Error: could not open STL file ", filepath
526
527 call s_mpi_abort()
528 end if
529
530 read (iunit, '(A)') line
531
532 close (iunit)
533
534 if (line(1:5) == "solid") then
535 call s_read_stl_ascii(filepath, model)
536 else
537 call s_read_stl_binary(filepath, model)
538 end if
539
540 end subroutine s_read_stl
541
542 !> Read an OBJ file.
543 impure subroutine s_read_obj(filepath, model)
544
545 character(LEN=*), intent(in) :: filepath
546 type(t_model), intent(out) :: model
547 integer :: i, j, k, l, iv3, iunit, iostat, nvertices
548 real(wp), dimension(1:3), allocatable :: vertices(:,:)
549 character(80) :: line
550
551 open (newunit=iunit, file=filepath, action='READ', form='FORMATTED', status='OLD', iostat=iostat, access='STREAM')
552
553 if (iostat /= 0) then
554 print *, "Error: could not open model file ", filepath
555
556 call s_mpi_abort()
557 end if
558
559 nvertices = 0
560 model%ntrs = 0
561 do
562 if (.not. f_read_line(iunit, line)) exit
563
564 select case (line(1:2))
565 case ("v ")
566 nvertices = nvertices + 1
567 case ("f ")
568 model%ntrs = model%ntrs + 1
569 end select
570 end do
571
572 rewind(iunit)
573
574 allocate (vertices(nvertices,1:3))
575 allocate (model%trs(model%ntrs))
576
577 i = 1
578 j = 1
579
580 do
581 if (.not. f_read_line(iunit, line)) exit
582
583 select case (line(1:2))
584 case ("g ")
585 case ("vn")
586 case ("vt")
587 case ("l ")
588 case ("v ")
589 read (line(3:), *) vertices(i,:)
590 i = i + 1
591 case ("f ")
592 read (line(3:), *) k, l, iv3
593 model%trs(j)%v(1,:) = vertices(k,:)
594 model%trs(j)%v(2,:) = vertices(l,:)
595 model%trs(j)%v(3,:) = vertices(iv3,:)
596 j = j + 1
597 case default
598 print *, "Error: unknown line type in OBJ file ", filepath
599 print *, "Line: ", line
600
601 call s_mpi_abort()
602 end select
603 end do
604
605 deallocate (vertices)
606
607 close (iunit)
608
609 end subroutine s_read_obj
610
611 !> Read a mesh from a file.
612 impure function f_model_read(filepath) result(model)
613
614 character(LEN=*), intent(in) :: filepath
615 type(t_model) :: model
616
617 select case (filepath(len(trim(filepath)) - 3:len(trim(filepath))))
618 case (".stl")
619 call s_read_stl(filepath, model)
620 case (".obj")
621 call s_read_obj(filepath, model)
622 case default
623 print *, "Error: unknown model file format for file ", filepath
624
625 call s_mpi_abort()
626 end select
627
628 end function f_model_read
629
630 !> Write a binary STL file.
631 impure subroutine s_write_stl(filepath, model)
632
633 character(LEN=*), intent(in) :: filepath
634 type(t_model), intent(in) :: model
635 integer :: i, j, iunit, iostat
636 character(kind=c_char, len=80), parameter :: header = "Model file written by MFC."
637 integer(kind=c_int32_t) :: ntriangles
638 real(wp) :: normal(3), v(3)
639 integer(kind=c_int16_t) :: attribute
640
641 open (newunit=iunit, file=filepath, action='WRITE', form='UNFORMATTED', iostat=iostat, access='STREAM')
642
643 if (iostat /= 0) then
644 print *, "Error: could not open STL file ", filepath
645
646 call s_mpi_abort()
647 end if
648
649 ntriangles = model%ntrs
650 write (iunit, iostat=iostat) header, ntriangles
651
652 if (iostat /= 0) then
653 print *, "Error: could not write header to STL file ", filepath
654
655 call s_mpi_abort()
656 end if
657
658 do i = 1, model%ntrs
659 normal = model%trs(i)%n
660 write (iunit) normal
661
662 do j = 1, 3
663 v = model%trs(i)%v(j,:)
664 write (iunit) v(:)
665 end do
666
667 attribute = 0
668 write (iunit) attribute
669 end do
670
671 close (iunit)
672
673 end subroutine s_write_stl
674
675 !> Write an OBJ file.
676 impure subroutine s_write_obj(filepath, model)
677
678 character(LEN=*), intent(in) :: filepath
679 type(t_model), intent(in) :: model
680 integer :: iunit, iostat
681 integer :: i, j
682
683 open (newunit=iunit, file=filepath, action='WRITE', form='FORMATTED', iostat=iostat, access='STREAM')
684
685 if (iostat /= 0) then
686 print *, "Error: could not open OBJ file ", filepath
687
688 call s_mpi_abort()
689 end if
690
691 write (iunit, '(A)') "# Model file written by MFC."
692
693 do i = 1, model%ntrs
694 do j = 1, 3
695 write (iunit, '(A, " ", (f30.20), " ", (f30.20), " ", (f30.20))') "v", model%trs(i)%v(j, 1), model%trs(i)%v(j, &
696 & 2), model%trs(i)%v(j, 3)
697 end do
698
699 write (iunit, '(A, " ", I0, " ", I0, " ", I0)') "f", i*3 - 2, i*3 - 1, i*3
700 end do
701
702 close (iunit)
703
704 end subroutine s_write_obj
705
706 !> Write a mesh to a file.
707 impure subroutine s_model_write(filepath, model)
708
709 character(LEN=*), intent(in) :: filepath
710 type(t_model), intent(in) :: model
711
712 select case (filepath(len(trim(filepath)) - 3:len(trim(filepath))))
713 case (".stl")
714 call s_write_stl(filepath, model)
715 case (".obj")
716 call s_write_obj(filepath, model)
717 case default
718 print *, "Error: unknown model file format for file ", filepath
719
720 call s_mpi_abort()
721 end select
722
723 end subroutine s_model_write
724
725 !> Free the memory allocated for an STL mesh.
726 subroutine s_model_free(model)
727
728 type(t_model), intent(inout) :: model
729
730 deallocate (model%trs)
731
732 end subroutine s_model_free
733
734 !> Read the next non-blank, non-comment line from an STL or OBJ model file.
735 impure function f_read_line(iunit, line) result(bIsLine)
736
737 integer, intent(in) :: iunit
738 character(80), intent(out) :: line
739 logical :: bisline
740 integer :: iostat
741
742 bisline = .true.
743
744 do
745 read (iunit, '(A)', iostat=iostat) line
746
747 if (iostat < 0) then
748 bisline = .false.
749 exit
750 end if
751
752 line = adjustl(trim(line))
753
754 if (len(trim(line)) == 0) cycle
755 if (line(1:5) == "solid") cycle
756 if (line(1:1) == "#") cycle
757
758 exit
759 end do
760
761 end function f_read_line
762
763 !> Read the next non-comment line from a model file, using a buffered look-ahead mechanism.
764 impure subroutine s_skip_ignored_lines(iunit, buffered_line, is_buffered)
765
766 integer, intent(in) :: iunit
767 character(80), intent(inout) :: buffered_line
768 logical, intent(inout) :: is_buffered
769 character(80) :: line
770
771 if (is_buffered) then
772 line = buffered_line
773 is_buffered = .false.
774 else
775 if (.not. f_read_line(iunit, line)) return
776 end if
777
778 buffered_line = line
779 is_buffered = .true.
780
781 end subroutine s_skip_ignored_lines
782
783 !> Generate a pseudo-random number using a seed-based xorshift, replacing the Fortran intrinsic which is incompatible with GPU
784 !! routines
785 function f_model_random_number(seed) result(rval)
786
787 ! $:GPU_ROUTINE(parallelism='[seq]')
788
789 integer, intent(inout) :: seed
790 real(wp) :: rval
791
792 seed = ieor(seed, ishft(seed, 13))
793 seed = ieor(seed, ishft(seed, -17))
794 seed = ieor(seed, ishft(seed, 5))
795
796 rval = abs(real(seed, wp))/real(huge(seed), wp)
797
798 end function f_model_random_number
799
800 !> Determine whether a point is inside a model using stochastic ray casting.
801 !! @param spc Number of samples per cell.
802 impure function f_model_is_inside(model, point, spacing, spc) result(fraction)
803
804 ! $:GPU_ROUTINE(parallelism='[seq]')
805
806 type(t_model), intent(in) :: model
807 real(wp), dimension(1:3), intent(in) :: point
808 real(wp), dimension(1:3), intent(in) :: spacing
809 integer, intent(in) :: spc
810 real(wp) :: phi, theta
811 integer :: rand_seed
812 real(wp) :: fraction
813 type(t_ray) :: ray
814 integer :: i, j, k, ninorout, nhits
815 real(wp), dimension(1:spc,1:3) :: ray_origins, ray_dirs
816
817 rand_seed = int(point(1)*73856093._wp) + int(point(2)*19349663._wp) + int(point(3)*83492791._wp)
818 if (rand_seed == 0) rand_seed = 1
819
820 ! generate our random collection or rays
821 do i = 1, spc
822 do k = 1, 3
823 ! random jitter in the origin helps us estimate volume fraction instead of only at the cell center
824 ray_origins(i, k) = point(k) + (f_model_random_number(rand_seed) - 0.5_wp)*spacing(k)
825 ! cast sample rays in all directions
826 ray_dirs(i, k) = f_model_random_number(rand_seed) - 0.5_wp
827 end do
828 ray_dirs(i,:) = ray_dirs(i,:)/sqrt(sum(ray_dirs(i,:)*ray_dirs(i,:)))
829 end do
830
831 ! ray trace
832 ninorout = 0
833 do i = 1, spc
834 ray%o = ray_origins(i,:)
835 ray%d = ray_dirs(i,:)
836
837 nhits = 0
838 do j = 1, model%ntrs
839 ! count the number of triangles this ray intersects
840 if (f_intersects_triangle(ray, model%trs(j)) == 1) then
841 nhits = nhits + 1
842 end if
843 end do
844
845 ! if the ray hits an odd number of triangles on its way out, then it must be on the inside of the model
846 ninorout = ninorout + mod(nhits, 2)
847 end do
848
849 fraction = real(ninorout)/real(spc)
850
851 end function f_model_is_inside
852
853 !> Determine if a point is inside a surface using the generalized winding number (Jacobson et al., SIGGRAPH 2013). In 3D, sums
854 !! the solid angle subtended by each triangle (Van Oosterom-Strackee formula). In 2D (p==0), sums the signed angle subtended by
855 !! each boundary edge. Returns ~1.0 inside, ~0.0 outside. Unlike ray casting, this is robust to small triangles/edges and vertex
856 !! winding order.
857 !! @return fraction Winding number (~1.0 inside, ~0.0 outside).
858 function f_model_is_inside_flat(ntrs, pid, point) result(fraction)
859
860
861# 541 "/home/runner/work/MFC/MFC/src/common/m_model.fpp"
862#if MFC_OpenACC
863# 541 "/home/runner/work/MFC/MFC/src/common/m_model.fpp"
864!$acc routine seq
865# 541 "/home/runner/work/MFC/MFC/src/common/m_model.fpp"
866#elif MFC_OpenMP
867# 541 "/home/runner/work/MFC/MFC/src/common/m_model.fpp"
868
869# 541 "/home/runner/work/MFC/MFC/src/common/m_model.fpp"
870
871# 541 "/home/runner/work/MFC/MFC/src/common/m_model.fpp"
872!$omp declare target device_type(any)
873# 541 "/home/runner/work/MFC/MFC/src/common/m_model.fpp"
874#endif
875
876 integer, intent(in) :: ntrs
877 integer, intent(in) :: pid
878 real(wp), dimension(1:3), intent(in) :: point
879 real(wp) :: fraction
880 real(wp) :: r1(3), r2(3), r3(3)
881 real(wp) :: r1_mag, r2_mag, r3_mag
882 real(wp) :: numerator, denominator
883 real(wp) :: d1(2), d2(2)
884 integer :: q
885
886 fraction = 0.0_wp
887
888 if (p == 0) then
889 ! 2D winding number: sum signed angles subtended by each boundary edge at the query point.
890 do q = 1, gpu_boundary_edge_count(pid)
891 d1(1) = gpu_boundary_v(q, 1, 1, pid) - point(1)
892 d1(2) = gpu_boundary_v(q, 1, 2, pid) - point(2)
893 d2(1) = gpu_boundary_v(q, 2, 1, pid) - point(1)
894 d2(2) = gpu_boundary_v(q, 2, 2, pid) - point(2)
895
896 ! Signed angle = atan2(d1 x d2, d1 . d2)
897 fraction = fraction + atan2(d1(1)*d2(2) - d1(2)*d2(1), d1(1)*d2(1) + d1(2)*d2(2))
898 end do
899
900 ! 2D winding number = total angle / (2*pi)
901 fraction = fraction/(2.0_wp*acos(-1.0_wp))
902 else
903 ! 3D winding number: sum solid angles via Van Oosterom-Strackee formula.
904 do q = 1, ntrs
905 r1 = gpu_trs_v(1,:,q, pid) - point
906 r2 = gpu_trs_v(2,:,q, pid) - point
907 r3 = gpu_trs_v(3,:,q, pid) - point
908
909 r1_mag = sqrt(dot_product(r1, r1))
910 r2_mag = sqrt(dot_product(r2, r2))
911 r3_mag = sqrt(dot_product(r3, r3))
912
913 ! Skip if query point is coincident with a vertex (magnitudes are zero/subnormal).
914 if (r1_mag*r2_mag*r3_mag < tiny(1.0_wp)) cycle
915
916 ! tan(Omega/2) = numerator / denominator numerator = scalar triple product r1 . (r2 x r3)
917 numerator = r1(1)*(r2(2)*r3(3) - r2(3)*r3(2)) + r1(2)*(r2(3)*r3(1) - r2(1)*r3(3)) + r1(3)*(r2(1)*r3(2) - r2(2) &
918 & *r3(1))
919
920 denominator = r1_mag*r2_mag*r3_mag + dot_product(r1, r2)*r3_mag + dot_product(r2, r3)*r1_mag + dot_product(r3, &
921 & r1)*r2_mag
922
923 fraction = fraction + atan2(numerator, denominator)
924 end do
925
926 ! Each atan2 returns Omega/2 per triangle; divide by 2*pi to get winding number = sum(Omega)/(4*pi).
927 fraction = fraction/(2.0_wp*acos(-1.0_wp))
928 end if
929
930 end function f_model_is_inside_flat
931
932 !> Check if a ray intersects a triangle using the Moller-Trumbore algorithm (barycentric coordinates). Unlike the previous
933 !! cross-product sign test, this is vertex winding-order independent.
934 !! @return 1 if the ray intersects the triangle, 0 otherwise.
935 function f_intersects_triangle(ray, triangle) result(intersects)
936
937
938# 604 "/home/runner/work/MFC/MFC/src/common/m_model.fpp"
939#if MFC_OpenACC
940# 604 "/home/runner/work/MFC/MFC/src/common/m_model.fpp"
941!$acc routine seq
942# 604 "/home/runner/work/MFC/MFC/src/common/m_model.fpp"
943#elif MFC_OpenMP
944# 604 "/home/runner/work/MFC/MFC/src/common/m_model.fpp"
945
946# 604 "/home/runner/work/MFC/MFC/src/common/m_model.fpp"
947
948# 604 "/home/runner/work/MFC/MFC/src/common/m_model.fpp"
949!$omp declare target device_type(any)
950# 604 "/home/runner/work/MFC/MFC/src/common/m_model.fpp"
951#endif
952
953 type(t_ray), intent(in) :: ray
954 type(t_triangle), intent(in) :: triangle
955 integer :: intersects
956 real(wp) :: edge1(3), edge2(3), h(3), s(3), q(3)
957 real(wp) :: a, f, u, v, t
958
959 intersects = 0
960
961 edge1 = triangle%v(2,:) - triangle%v(1,:)
962 edge2 = triangle%v(3,:) - triangle%v(1,:)
963 h = f_cross(ray%d, edge2)
964 a = dot_product(edge1, h)
965
966 ! Ray nearly parallel to triangle plane. In single precision builds epsilon(1.0) ~ 1.2e-7, so use 10*epsilon as a floor.
967 if (abs(a) < max(1e-7_wp, 10.0_wp*epsilon(1.0_wp))) return
968
969 f = 1.0_wp/a
970 s = ray%o - triangle%v(1,:)
971 u = f*dot_product(s, h)
972
973 if (u < 0.0_wp .or. u > 1.0_wp) return
974
975 q = f_cross(s, edge1)
976 v = f*dot_product(ray%d, q)
977
978 if (v < 0.0_wp .or. u + v > 1.0_wp) return
979
980 t = f*dot_product(edge2, q)
981
982 if (t > 0.0_wp) intersects = 1
983
984 end function f_intersects_triangle
985
986 !> Check and label edges shared by two or more triangle facets of the 2D STL model.
987 subroutine s_check_boundary(model, boundary_v, boundary_vertex_count, boundary_edge_count)
988
989 type(t_model), intent(in) :: model
990 real(wp), allocatable, intent(out), dimension(:,:,:) :: boundary_v !< Output boundary vertices/normals
991 integer, intent(out) :: boundary_vertex_count, boundary_edge_count !< Output boundary vertex/edge count
992 integer :: i, j !< Model index iterator
993 integer :: edge_count, edge_index, store_index !< Boundary edge index iterator
994 real(wp), dimension(1:2,1:2) :: edge !< Edge end points buffer
995 real(wp), dimension(1:2) :: boundary_edge !< Boundary edge end points buffer
996 real(wp), dimension(1:(3*model%ntrs),1:2,1:2) :: temp_boundary_v !< Temporary boundary vertex buffer
997 integer, dimension(1:(3*model%ntrs)) :: edge_occurrence !< The manifoldness of the edges
998 real(wp) :: edgetan, initial, v_norm, xnormal, ynormal !< The manifoldness of the edges
999 ! Total number of edges in 2D STL
1000
1001 edge_count = 3*model%ntrs
1002
1003 ! Initialize edge_occurrence array to zero
1004 edge_occurrence = 0
1005 edge_index = 0
1006
1007 ! Collect all edges of all triangles and store them
1008 do i = 1, model%ntrs
1009 ! First edge (v1, v2)
1010 edge(1,1:2) = model%trs(i)%v(1,1:2)
1011 edge(2,1:2) = model%trs(i)%v(2,1:2)
1012 call s_register_edge(temp_boundary_v, edge, edge_index, edge_count)
1013
1014 ! Second edge (v2, v3)
1015 edge(1,1:2) = model%trs(i)%v(2,1:2)
1016 edge(2,1:2) = model%trs(i)%v(3,1:2)
1017 call s_register_edge(temp_boundary_v, edge, edge_index, edge_count)
1018
1019 ! Third edge (v3, v1)
1020 edge(1,1:2) = model%trs(i)%v(3,1:2)
1021 edge(2,1:2) = model%trs(i)%v(1,1:2)
1022 call s_register_edge(temp_boundary_v, edge, edge_index, edge_count)
1023 end do
1024
1025 ! Check all edges and count repeated edges
1026
1027# 679 "/home/runner/work/MFC/MFC/src/common/m_model.fpp"
1028
1029# 679 "/home/runner/work/MFC/MFC/src/common/m_model.fpp"
1030#if defined(MFC_OpenACC)
1031# 679 "/home/runner/work/MFC/MFC/src/common/m_model.fpp"
1032!$acc parallel loop collapse(2) gang vector default(present) private(i, j) copy(temp_boundary_v, edge_occurrence)
1033# 679 "/home/runner/work/MFC/MFC/src/common/m_model.fpp"
1034#elif defined(MFC_OpenMP)
1035# 679 "/home/runner/work/MFC/MFC/src/common/m_model.fpp"
1036
1037# 679 "/home/runner/work/MFC/MFC/src/common/m_model.fpp"
1038
1039# 679 "/home/runner/work/MFC/MFC/src/common/m_model.fpp"
1040
1041# 679 "/home/runner/work/MFC/MFC/src/common/m_model.fpp"
1042!$omp target teams loop defaultmap(firstprivate:scalar) bind(teams,parallel) collapse(2) defaultmap(tofrom:aggregate) defaultmap(tofrom:allocatable) defaultmap(tofrom:pointer) private(i, j) map(tofrom:temp_boundary_v, edge_occurrence)
1043# 679 "/home/runner/work/MFC/MFC/src/common/m_model.fpp"
1044#endif
1045 do i = 1, edge_count
1046 do j = 1, edge_count
1047 if (i /= j) then
1048 if (((abs(temp_boundary_v(i, 1, 1) - temp_boundary_v(j, 1, &
1049 & 1)) < threshold_edge_zero) .and. (abs(temp_boundary_v(i, 1, 2) - temp_boundary_v(j, 1, &
1050 & 2)) < threshold_edge_zero) .and. (abs(temp_boundary_v(i, 2, 1) - temp_boundary_v(j, 2, &
1051 & 1)) < threshold_edge_zero) .and. (abs(temp_boundary_v(i, 2, 2) - temp_boundary_v(j, 2, &
1052 & 2)) < threshold_edge_zero)) .or. ((abs(temp_boundary_v(i, 1, 1) - temp_boundary_v(j, 2, &
1053 & 1)) < threshold_edge_zero) .and. (abs(temp_boundary_v(i, 1, 2) - temp_boundary_v(j, 2, &
1054 & 2)) < threshold_edge_zero) .and. (abs(temp_boundary_v(i, 2, 1) - temp_boundary_v(j, 1, &
1055 & 1)) < threshold_edge_zero) .and. (abs(temp_boundary_v(i, 2, 2) - temp_boundary_v(j, 1, &
1056 & 2)) < threshold_edge_zero))) then
1057
1058# 692 "/home/runner/work/MFC/MFC/src/common/m_model.fpp"
1059#if defined(MFC_OpenACC)
1060# 692 "/home/runner/work/MFC/MFC/src/common/m_model.fpp"
1061!$acc atomic update
1062# 692 "/home/runner/work/MFC/MFC/src/common/m_model.fpp"
1063#elif defined(MFC_OpenMP)
1064# 692 "/home/runner/work/MFC/MFC/src/common/m_model.fpp"
1065!$omp atomic update
1066# 692 "/home/runner/work/MFC/MFC/src/common/m_model.fpp"
1067#endif
1068 edge_occurrence(i) = edge_occurrence(i) + 1
1069 end if
1070 end if
1071 end do
1072 end do
1073
1074# 698 "/home/runner/work/MFC/MFC/src/common/m_model.fpp"
1075#if defined(MFC_OpenACC)
1076# 698 "/home/runner/work/MFC/MFC/src/common/m_model.fpp"
1077!$acc end parallel loop
1078# 698 "/home/runner/work/MFC/MFC/src/common/m_model.fpp"
1079#elif defined(MFC_OpenMP)
1080# 698 "/home/runner/work/MFC/MFC/src/common/m_model.fpp"
1081
1082# 698 "/home/runner/work/MFC/MFC/src/common/m_model.fpp"
1083!$omp end target teams loop
1084# 698 "/home/runner/work/MFC/MFC/src/common/m_model.fpp"
1085#endif
1086
1087 ! Count the number of boundary vertices/edges
1088 boundary_vertex_count = 0
1089 boundary_edge_count = 0
1090
1091 do i = 1, edge_count
1092 if (edge_occurrence(i) == 0) then
1093 boundary_vertex_count = boundary_vertex_count + 2
1094 boundary_edge_count = boundary_edge_count + 1
1095 end if
1096 end do
1097
1098 ! Allocate the boundary_v array based on the number of boundary edges
1099 allocate (boundary_v(boundary_edge_count,1:3,1:2))
1100
1101 ! Store boundary vertices
1102 store_index = 0
1103 do i = 1, edge_count
1104 if (edge_occurrence(i) == 0) then
1105 store_index = store_index + 1
1106 boundary_v(store_index, 1,1:2) = temp_boundary_v(i, 1,1:2)
1107 boundary_v(store_index, 2,1:2) = temp_boundary_v(i, 2,1:2)
1108 end if
1109 end do
1110
1111 ! Find/store the normal vector of the boundary edges
1112 do i = 1, boundary_edge_count
1113 boundary_edge(1) = boundary_v(i, 2, 1) - boundary_v(i, 1, 1)
1114 boundary_edge(2) = boundary_v(i, 2, 2) - boundary_v(i, 1, 2)
1115 edgetan = boundary_edge(1)/sign(max(sgm_eps, abs(boundary_edge(2))), boundary_edge(2))
1116
1117 if (abs(boundary_edge(2)) < threshold_vector_zero) then
1118 if (edgetan > 0._wp) then
1119 ynormal = -1
1120 xnormal = 0._wp
1121 else
1122 ynormal = 1
1123 xnormal = 0._wp
1124 end if
1125 else
1126 initial = boundary_edge(2)
1127 ynormal = -edgetan*initial
1128 xnormal = initial
1129 end if
1130
1131 v_norm = sqrt(xnormal**2 + ynormal**2)
1132 boundary_v(i, 3, 1) = xnormal/v_norm
1133 boundary_v(i, 3, 2) = ynormal/v_norm
1134 end do
1135
1136 end subroutine s_check_boundary
1137
1138 !> Append the edge end vertices to a temporary buffer.
1139 subroutine s_register_edge(temp_boundary_v, edge, edge_index, edge_count)
1140
1141 integer, intent(inout) :: edge_index !< Edge index iterator
1142 integer, intent(inout) :: edge_count !< Total number of edges
1143 real(wp), intent(in), dimension(1:2,1:2) :: edge !< Edges end points to be registered
1144 real(wp), dimension(1:edge_count,1:2,1:2), intent(inout) :: temp_boundary_v !< Temporary edge end vertex buffer
1145 ! Increment edge index and store the edge
1146
1147 edge_index = edge_index + 1
1148 temp_boundary_v(edge_index, 1,1:2) = edge(1,1:2)
1149 temp_boundary_v(edge_index, 2,1:2) = edge(2,1:2)
1150
1151 end subroutine s_register_edge
1152
1153 !> Determine the levelset distance and normals of 3D models by computing the exact closest point via projection onto triangle
1154 !! surfaces.
1155 subroutine s_distance_normals_3d(ntrs, pid, point, normals, distance)
1156
1157
1158# 770 "/home/runner/work/MFC/MFC/src/common/m_model.fpp"
1159#if MFC_OpenACC
1160# 770 "/home/runner/work/MFC/MFC/src/common/m_model.fpp"
1161!$acc routine seq
1162# 770 "/home/runner/work/MFC/MFC/src/common/m_model.fpp"
1163#elif MFC_OpenMP
1164# 770 "/home/runner/work/MFC/MFC/src/common/m_model.fpp"
1165
1166# 770 "/home/runner/work/MFC/MFC/src/common/m_model.fpp"
1167
1168# 770 "/home/runner/work/MFC/MFC/src/common/m_model.fpp"
1169!$omp declare target device_type(any)
1170# 770 "/home/runner/work/MFC/MFC/src/common/m_model.fpp"
1171#endif
1172
1173 integer, intent(in) :: ntrs
1174 integer, intent(in) :: pid
1175 real(wp), dimension(1:3), intent(in) :: point
1176 real(wp), dimension(1:3), intent(out) :: normals
1177 real(wp), intent(out) :: distance
1178 integer :: i, j, l
1179 real(wp) :: dist_min, dist_proj, dist_v, dist_e, t
1180 real(wp) :: v1(1:3), v2(1:3), v3(1:3)
1181 real(wp) :: e0(1:3), e1(1:3), pv(1:3)
1182 real(wp) :: n(1:3), proj(1:3), norm_vec(1:3)
1183 real(wp) :: d, ndot, denom, norm_mag
1184 real(wp) :: u, v_bary, w
1185 real(wp) :: l00, l01, l11, l20, l21
1186 real(wp) :: edge(1:3), pe(1:3)
1187 real(wp) :: verts(1:3,1:3)
1188
1189 dist_min = initial_distance_buffer
1190 normals = 0._wp
1191
1192 do i = 1, ntrs
1193 ! Triangle vertices
1194 v1(:) = gpu_trs_v(1,:,i, pid)
1195 v2(:) = gpu_trs_v(2,:,i, pid)
1196 v3(:) = gpu_trs_v(3,:,i, pid)
1197
1198 ! Triangle normal
1199 n(:) = gpu_trs_n(:,i, pid)
1200
1201 ! Project point onto triangle plane
1202 pv(:) = point(:) - v1(:)
1203 d = dot_product(pv, n)
1204 if (abs(d) >= dist_min) cycle ! minimum distance is not small enough, no need to check validity
1205 proj(:) = point(:) - d*n(:)
1206
1207 ! Check if projection is inside triangle using barycentric coordinates
1208 e0(:) = v2(:) - v1(:)
1209 e1(:) = v3(:) - v1(:)
1210 pv(:) = proj(:) - v1(:)
1211
1212 l00 = dot_product(e0, e0)
1213 l01 = dot_product(e0, e1)
1214 l11 = dot_product(e1, e1)
1215 l20 = dot_product(pv, e0)
1216 l21 = dot_product(pv, e1)
1217
1218 denom = l00*l11 - l01*l01
1219
1220 ! compute the barycentric coordinates of the projection in the triangle
1221 if (abs(denom) > 0._wp) then
1222 v_bary = (l11*l20 - l01*l21)/denom
1223 w = (l00*l21 - l01*l20)/denom
1224 u = 1._wp - v_bary - w
1225 else
1226 u = -1._wp
1227 v_bary = -1._wp
1228 w = -1._wp
1229 end if
1230
1231 ! If projection is inside triangle
1232 if (u >= 0._wp .and. v_bary >= 0._wp .and. w >= 0._wp) then
1233 dist_proj = sqrt((point(1) - proj(1))**2 + (point(2) - proj(2))**2 + (point(3) - proj(3))**2)
1234
1235 if (dist_proj < dist_min) then
1236 dist_min = dist_proj
1237 normals(:) = n(:)
1238 end if
1239 else
1240 ! Projection outside triangle: check edges and vertices
1241 verts(:,1) = v1(:)
1242 verts(:,2) = v2(:)
1243 verts(:,3) = v3(:)
1244
1245 ! Check three edges
1246 do j = 1, 3
1247 edge(:) = verts(:,mod(j, 3) + 1) - verts(:,j)
1248 pe(:) = point(:) - verts(:,j)
1249
1250 t = dot_product(pe, edge)/max(dot_product(edge, edge), 1.e-30_wp)
1251
1252 if (t >= 0._wp .and. t <= 1._wp) then
1253 proj(:) = verts(:,j) + t*edge(:)
1254 dist_e = sqrt((point(1) - proj(1))**2 + (point(2) - proj(2))**2 + (point(3) - proj(3))**2)
1255
1256 if (dist_e < dist_min) then
1257 dist_min = dist_e
1258 norm_vec(:) = proj(:) - point(:)
1259 if (dist_e > 0._wp) norm_vec = norm_vec/dist_e
1260 ! Snap to triangle normal if nearly parallel
1261 if (f_approx_equal(dot_product(norm_vec, n), 1._wp)) then
1262 normals(:) = n(:)
1263 else
1264 normals(:) = norm_vec(:)
1265 end if
1266 end if
1267 else if (t < 0._wp) then
1268 dist_v = sqrt((point(1) - verts(1, j))**2 + (point(2) - verts(2, j))**2 + (point(3) - verts(3, j))**2)
1269
1270 if (dist_v < dist_min) then
1271 dist_min = dist_v
1272 norm_vec(:) = verts(:,j) - point(:)
1273 norm_mag = sqrt(dot_product(norm_vec, norm_vec))
1274 if (norm_mag > 0._wp) norm_vec = norm_vec/norm_mag
1275 normals(:) = norm_vec(:)
1276 end if
1277 else
1278 dist_v = sqrt((point(1) - verts(1, mod(j, 3) + 1))**2 + (point(2) - verts(2, mod(j, &
1279 & 3) + 1))**2 + (point(3) - verts(3, mod(j, 3) + 1))**2)
1280
1281 if (dist_v < dist_min) then
1282 dist_min = dist_v
1283 norm_vec(:) = verts(:,mod(j, 3) + 1) - point(:)
1284 norm_mag = sqrt(dot_product(norm_vec, norm_vec))
1285 if (norm_mag > 0._wp) norm_vec = norm_vec/norm_mag
1286 normals(:) = norm_vec(:)
1287 end if
1288 end if
1289 end do
1290 end if
1291 end do
1292
1293 distance = dist_min
1294
1295 end subroutine s_distance_normals_3d
1296
1297 !> Determine the levelset distance and normals of 2D models by computing the exact closest point via projection onto boundary
1298 !! edges.
1299 subroutine s_distance_normals_2d(pid, boundary_edge_count, point, normals, distance)
1300
1301
1302# 900 "/home/runner/work/MFC/MFC/src/common/m_model.fpp"
1303#if MFC_OpenACC
1304# 900 "/home/runner/work/MFC/MFC/src/common/m_model.fpp"
1305!$acc routine seq
1306# 900 "/home/runner/work/MFC/MFC/src/common/m_model.fpp"
1307#elif MFC_OpenMP
1308# 900 "/home/runner/work/MFC/MFC/src/common/m_model.fpp"
1309
1310# 900 "/home/runner/work/MFC/MFC/src/common/m_model.fpp"
1311
1312# 900 "/home/runner/work/MFC/MFC/src/common/m_model.fpp"
1313!$omp declare target device_type(any)
1314# 900 "/home/runner/work/MFC/MFC/src/common/m_model.fpp"
1315#endif
1316
1317 integer, intent(in) :: pid
1318 integer, intent(in) :: boundary_edge_count
1319 real(wp), dimension(1:3), intent(in) :: point
1320 real(wp), dimension(1:3), intent(out) :: normals
1321 real(wp), intent(out) :: distance
1322 integer :: i
1323 real(wp) :: dist_min, dist, t, norm_mag
1324 real(wp) :: v1(1:2), v2(1:2), edge(1:2), pv(1:2)
1325 real(wp) :: edge_len_sq, proj(1:2), norm(1:2), c
1326
1327 dist_min = initial_distance_buffer
1328 normals = 0._wp
1329 norm = 0._wp
1330
1331 do i = 1, boundary_edge_count
1332 ! Edge endpoints
1333 v1(1) = gpu_boundary_v(i, 1, 1, pid)
1334 v1(2) = gpu_boundary_v(i, 1, 2, pid)
1335 v2(1) = gpu_boundary_v(i, 2, 1, pid)
1336 v2(2) = gpu_boundary_v(i, 2, 2, pid)
1337
1338 ! Edge vector and point-to-v1 vector
1339 edge = v2 - v1
1340 pv(1) = point(1) - v1(1)
1341 pv(2) = point(2) - v1(2)
1342 edge_len_sq = dot_product(edge, edge)
1343
1344 ! Parameter of projection onto the edge line
1345 if (edge_len_sq > 0._wp) then
1346 t = dot_product(pv, edge)/edge_len_sq
1347 else
1348 t = 0._wp
1349 end if
1350
1351 ! Check if projection falls within the segment
1352 if (t >= 0._wp .and. t <= 1._wp) then
1353 proj = v1 + t*edge
1354 dist = sqrt((point(1) - proj(1))**2 + (point(2) - proj(2))**2)
1355 norm(1) = gpu_boundary_v(i, 3, 1, pid)
1356 norm(2) = gpu_boundary_v(i, 3, 2, pid)
1357 else if (t < 0._wp) then ! negative t means that v1 is the closest point on the edge
1358 dist = sqrt((point(1) - v1(1))**2 + (point(2) - v1(2))**2)
1359 norm(1) = v1(1) - point(1)
1360 norm(2) = v1(2) - point(2)
1361 norm = norm/dist
1362 else ! t > 1 means that v2 is the closest point on the line edge
1363 dist = sqrt((point(1) - v2(1))**2 + (point(2) - v2(2))**2)
1364 norm(1) = v2(1) - point(1)
1365 norm(2) = v2(2) - point(2)
1366 norm = norm/dist
1367 end if
1368
1369 if (dist < dist_min) then
1370 dist_min = dist
1371 normals(1) = norm(1)
1372 normals(2) = norm(2)
1373 end if
1374 end do
1375
1376 distance = dist_min
1377
1378 end subroutine s_distance_normals_2d
1379
1380#ifdef MFC_SIMULATION
1381 !> Load, transform, and register STL/OBJ immersed-boundary models onto the simulation grid.
1382 subroutine s_instantiate_stl_models()
1383
1384 ! Variables for IBM+STL
1385 real(wp) :: normals(1:3) !< Boundary normal buffer
1386 integer :: boundary_vertex_count, boundary_edge_count, total_vertices !< Boundary vertex
1387 real(wp), allocatable, dimension(:,:,:) :: boundary_v !< Boundary vertex buffer
1388 real(wp) :: dx_local, dy_local, dz_local !< Levelset distance buffer
1389 integer :: i, j, k !< Generic loop iterators
1390 integer :: patch_id
1391 type(t_bbox) :: bbox, bbox_old
1392 type(t_model) :: model
1393 type(ic_model_parameters) :: params
1394 real(wp) :: eta
1395 real(wp), dimension(1:3) :: point, model_center
1396 real(wp) :: grid_mm(1:3,1:2)
1397 real(wp), dimension(1:4,1:4) :: transform, transform_n
1398
1399 dx_local = minval(dx); dy_local = minval(dy)
1400 if (p /= 0) dz_local = minval(dz)
1401
1402 allocate (stl_bounding_boxes(num_ibs,1:3,1:3))
1403
1404 do patch_id = 1, num_ibs
1405 if (patch_ib(patch_id)%geometry == 5 .or. patch_ib(patch_id)%geometry == 12) then
1406 allocate (models(patch_id)%model)
1407 print *, " * Reading model: " // trim(patch_ib(patch_id)%model_filepath)
1408
1409 model = f_model_read(patch_ib(patch_id)%model_filepath)
1410 params%scale(:) = patch_ib(patch_id)%model_scale(:)
1411 params%translate(:) = patch_ib(patch_id)%model_translate(:)
1412 params%rotate(:) = patch_ib(patch_id)%model_rotate(:)
1413 params%spc = patch_ib(patch_id)%model_spc
1414 params%threshold = patch_ib(patch_id)%model_threshold
1415
1416 if (f_approx_equal(dot_product(params%scale, params%scale), 0._wp)) then
1417 params%scale(:) = 1._wp
1418 end if
1419
1420 if (proc_rank == 0) then
1421 print *, " * Transforming model."
1422 end if
1423
1424 ! Get the model center before transforming the model
1425 bbox_old = f_create_bbox(model)
1426 model_center(1:3) = (bbox_old%min(1:3) + bbox_old%max(1:3))/2._wp
1427
1428 ! Compute the transform matrices for vertices and normals
1429 transform = f_create_transform_matrix(params, model_center)
1430 transform_n = f_create_transform_matrix(params)
1431
1432 call s_transform_model(model, transform, transform_n)
1433
1434 ! Recreate the bounding box after transformation
1435 bbox = f_create_bbox(model)
1436
1437 ! Show the number of vertices in the original STL model
1438 if (proc_rank == 0) then
1439 print *, ' * Number of input model vertices:', 3*model%ntrs
1440 end if
1441
1442 ! Need the cells that form the boundary of the flat projection in 2D
1443 if (p == 0) call s_check_boundary(model, boundary_v, boundary_vertex_count, boundary_edge_count)
1444
1445 ! Show the number of edges and boundary edges in 2D STL models
1446 if (proc_rank == 0 .and. p == 0) then
1447 print *, ' * Number of 2D model boundary edges:', boundary_edge_count
1448 end if
1449
1450 if (proc_rank == 0) then
1451 write (*, "(A, 3(2X, F20.10))") " > Model: Min:", bbox%min(1:3)
1452 write (*, "(A, 3(2X, F20.10))") " > Cen:", (bbox%min(1:3) + bbox%max(1:3))/2._wp
1453 write (*, "(A, 3(2X, F20.10))") " > Max:", bbox%max(1:3)
1454
1455 grid_mm(1,:) = (/minval(x_cc(0:m)) - 0.5_wp*dx_local, maxval(x_cc(0:m)) + 0.5_wp*dx_local/)
1456 grid_mm(2,:) = (/minval(y_cc(0:n)) - 0.5_wp*dy_local, maxval(y_cc(0:n)) + 0.5_wp*dy_local/)
1457
1458 if (p > 0) then
1459 grid_mm(3,:) = (/minval(z_cc(0:p)) - 0.5_wp*dz_local, maxval(z_cc(0:p)) + 0.5_wp*dz_local/)
1460 else
1461 grid_mm(3,:) = (/0._wp, 0._wp/)
1462 end if
1463
1464 write (*, "(A, 3(2X, F20.10))") " > Domain: Min:", grid_mm(:,1)
1465 write (*, "(A, 3(2X, F20.10))") " > Cen:", (grid_mm(:,1) + grid_mm(:,2))/2._wp
1466 write (*, "(A, 3(2X, F20.10))") " > Max:", grid_mm(:,2)
1467 end if
1468
1469 stl_bounding_boxes(patch_id, 1,1:3) = [bbox%min(1), (bbox%min(1) + bbox%max(1))/2._wp, bbox%max(1)]
1470 stl_bounding_boxes(patch_id, 2,1:3) = [bbox%min(2), (bbox%min(2) + bbox%max(2))/2._wp, bbox%max(2)]
1471 stl_bounding_boxes(patch_id, 3,1:3) = [bbox%min(3), (bbox%min(3) + bbox%max(3))/2._wp, bbox%max(3)]
1472
1473 models(patch_id)%model = model
1474 if (p == 0) then
1475 models(patch_id)%boundary_v = boundary_v
1476 models(patch_id)%boundary_edge_count = boundary_edge_count
1477 end if
1478 end if
1479 end do
1480
1481 ! Pack and upload flat arrays for GPU (AFTER the loop)
1482 block
1483 integer :: pid, max_ntrs
1484 integer :: max_bv1, max_bv2, max_bv3, max_iv1, max_iv2
1485
1486 max_ntrs = 0
1487 max_bv1 = 0; max_bv2 = 0; max_bv3 = 0
1488 max_iv1 = 0; max_iv2 = 0
1489
1490 do pid = 1, num_ibs
1491 if (allocated(models(pid)%model)) then
1492 call s_pack_model_for_gpu(models(pid))
1493 max_ntrs = max(max_ntrs, models(pid)%ntrs)
1494 end if
1495 if (allocated(models(pid)%boundary_v)) then
1496 max_bv1 = max(max_bv1, size(models(pid)%boundary_v, 1))
1497 max_bv2 = max(max_bv2, size(models(pid)%boundary_v, 2))
1498 max_bv3 = max(max_bv3, size(models(pid)%boundary_v, 3))
1499 end if
1500 end do
1501
1502 if (max_ntrs > 0) then
1503#ifdef MFC_DEBUG
1504# 1088 "/home/runner/work/MFC/MFC/src/common/m_model.fpp"
1505 block
1506# 1088 "/home/runner/work/MFC/MFC/src/common/m_model.fpp"
1507 use iso_fortran_env, only: output_unit
1508# 1088 "/home/runner/work/MFC/MFC/src/common/m_model.fpp"
1509
1510# 1088 "/home/runner/work/MFC/MFC/src/common/m_model.fpp"
1511 print *, 'm_model.fpp:1088: ', '@:ALLOCATE(gpu_ntrs(1:num_ibs))'
1512# 1088 "/home/runner/work/MFC/MFC/src/common/m_model.fpp"
1513
1514# 1088 "/home/runner/work/MFC/MFC/src/common/m_model.fpp"
1515 call flush (output_unit)
1516# 1088 "/home/runner/work/MFC/MFC/src/common/m_model.fpp"
1517 end block
1518# 1088 "/home/runner/work/MFC/MFC/src/common/m_model.fpp"
1519#endif
1520# 1088 "/home/runner/work/MFC/MFC/src/common/m_model.fpp"
1521 allocate (gpu_ntrs(1:num_ibs))
1522# 1088 "/home/runner/work/MFC/MFC/src/common/m_model.fpp"
1523
1524# 1088 "/home/runner/work/MFC/MFC/src/common/m_model.fpp"
1525
1526# 1088 "/home/runner/work/MFC/MFC/src/common/m_model.fpp"
1527#if defined(MFC_OpenACC)
1528# 1088 "/home/runner/work/MFC/MFC/src/common/m_model.fpp"
1529!$acc enter data create(gpu_ntrs)
1530# 1088 "/home/runner/work/MFC/MFC/src/common/m_model.fpp"
1531#elif defined(MFC_OpenMP)
1532# 1088 "/home/runner/work/MFC/MFC/src/common/m_model.fpp"
1533!$omp target enter data map(always,alloc:gpu_ntrs)
1534# 1088 "/home/runner/work/MFC/MFC/src/common/m_model.fpp"
1535#endif
1536#ifdef MFC_DEBUG
1537# 1089 "/home/runner/work/MFC/MFC/src/common/m_model.fpp"
1538 block
1539# 1089 "/home/runner/work/MFC/MFC/src/common/m_model.fpp"
1540 use iso_fortran_env, only: output_unit
1541# 1089 "/home/runner/work/MFC/MFC/src/common/m_model.fpp"
1542
1543# 1089 "/home/runner/work/MFC/MFC/src/common/m_model.fpp"
1544 print *, 'm_model.fpp:1089: ', '@:ALLOCATE(gpu_trs_v(1:3, 1:3, 1:max_ntrs, 1:num_ibs))'
1545# 1089 "/home/runner/work/MFC/MFC/src/common/m_model.fpp"
1546
1547# 1089 "/home/runner/work/MFC/MFC/src/common/m_model.fpp"
1548 call flush (output_unit)
1549# 1089 "/home/runner/work/MFC/MFC/src/common/m_model.fpp"
1550 end block
1551# 1089 "/home/runner/work/MFC/MFC/src/common/m_model.fpp"
1552#endif
1553# 1089 "/home/runner/work/MFC/MFC/src/common/m_model.fpp"
1554 allocate (gpu_trs_v(1:3, 1:3, 1:max_ntrs, 1:num_ibs))
1555# 1089 "/home/runner/work/MFC/MFC/src/common/m_model.fpp"
1556
1557# 1089 "/home/runner/work/MFC/MFC/src/common/m_model.fpp"
1558
1559# 1089 "/home/runner/work/MFC/MFC/src/common/m_model.fpp"
1560#if defined(MFC_OpenACC)
1561# 1089 "/home/runner/work/MFC/MFC/src/common/m_model.fpp"
1562!$acc enter data create(gpu_trs_v)
1563# 1089 "/home/runner/work/MFC/MFC/src/common/m_model.fpp"
1564#elif defined(MFC_OpenMP)
1565# 1089 "/home/runner/work/MFC/MFC/src/common/m_model.fpp"
1566!$omp target enter data map(always,alloc:gpu_trs_v)
1567# 1089 "/home/runner/work/MFC/MFC/src/common/m_model.fpp"
1568#endif
1569#ifdef MFC_DEBUG
1570# 1090 "/home/runner/work/MFC/MFC/src/common/m_model.fpp"
1571 block
1572# 1090 "/home/runner/work/MFC/MFC/src/common/m_model.fpp"
1573 use iso_fortran_env, only: output_unit
1574# 1090 "/home/runner/work/MFC/MFC/src/common/m_model.fpp"
1575
1576# 1090 "/home/runner/work/MFC/MFC/src/common/m_model.fpp"
1577 print *, 'm_model.fpp:1090: ', '@:ALLOCATE(gpu_trs_n(1:3, 1:max_ntrs, 1:num_ibs))'
1578# 1090 "/home/runner/work/MFC/MFC/src/common/m_model.fpp"
1579
1580# 1090 "/home/runner/work/MFC/MFC/src/common/m_model.fpp"
1581 call flush (output_unit)
1582# 1090 "/home/runner/work/MFC/MFC/src/common/m_model.fpp"
1583 end block
1584# 1090 "/home/runner/work/MFC/MFC/src/common/m_model.fpp"
1585#endif
1586# 1090 "/home/runner/work/MFC/MFC/src/common/m_model.fpp"
1587 allocate (gpu_trs_n(1:3, 1:max_ntrs, 1:num_ibs))
1588# 1090 "/home/runner/work/MFC/MFC/src/common/m_model.fpp"
1589
1590# 1090 "/home/runner/work/MFC/MFC/src/common/m_model.fpp"
1591
1592# 1090 "/home/runner/work/MFC/MFC/src/common/m_model.fpp"
1593#if defined(MFC_OpenACC)
1594# 1090 "/home/runner/work/MFC/MFC/src/common/m_model.fpp"
1595!$acc enter data create(gpu_trs_n)
1596# 1090 "/home/runner/work/MFC/MFC/src/common/m_model.fpp"
1597#elif defined(MFC_OpenMP)
1598# 1090 "/home/runner/work/MFC/MFC/src/common/m_model.fpp"
1599!$omp target enter data map(always,alloc:gpu_trs_n)
1600# 1090 "/home/runner/work/MFC/MFC/src/common/m_model.fpp"
1601#endif
1602#ifdef MFC_DEBUG
1603# 1091 "/home/runner/work/MFC/MFC/src/common/m_model.fpp"
1604 block
1605# 1091 "/home/runner/work/MFC/MFC/src/common/m_model.fpp"
1606 use iso_fortran_env, only: output_unit
1607# 1091 "/home/runner/work/MFC/MFC/src/common/m_model.fpp"
1608
1609# 1091 "/home/runner/work/MFC/MFC/src/common/m_model.fpp"
1610 print *, 'm_model.fpp:1091: ', '@:ALLOCATE(gpu_boundary_edge_count(1:num_ibs))'
1611# 1091 "/home/runner/work/MFC/MFC/src/common/m_model.fpp"
1612
1613# 1091 "/home/runner/work/MFC/MFC/src/common/m_model.fpp"
1614 call flush (output_unit)
1615# 1091 "/home/runner/work/MFC/MFC/src/common/m_model.fpp"
1616 end block
1617# 1091 "/home/runner/work/MFC/MFC/src/common/m_model.fpp"
1618#endif
1619# 1091 "/home/runner/work/MFC/MFC/src/common/m_model.fpp"
1620 allocate (gpu_boundary_edge_count(1:num_ibs))
1621# 1091 "/home/runner/work/MFC/MFC/src/common/m_model.fpp"
1622
1623# 1091 "/home/runner/work/MFC/MFC/src/common/m_model.fpp"
1624
1625# 1091 "/home/runner/work/MFC/MFC/src/common/m_model.fpp"
1626#if defined(MFC_OpenACC)
1627# 1091 "/home/runner/work/MFC/MFC/src/common/m_model.fpp"
1628!$acc enter data create(gpu_boundary_edge_count)
1629# 1091 "/home/runner/work/MFC/MFC/src/common/m_model.fpp"
1630#elif defined(MFC_OpenMP)
1631# 1091 "/home/runner/work/MFC/MFC/src/common/m_model.fpp"
1632!$omp target enter data map(always,alloc:gpu_boundary_edge_count)
1633# 1091 "/home/runner/work/MFC/MFC/src/common/m_model.fpp"
1634#endif
1635#ifdef MFC_DEBUG
1636# 1092 "/home/runner/work/MFC/MFC/src/common/m_model.fpp"
1637 block
1638# 1092 "/home/runner/work/MFC/MFC/src/common/m_model.fpp"
1639 use iso_fortran_env, only: output_unit
1640# 1092 "/home/runner/work/MFC/MFC/src/common/m_model.fpp"
1641
1642# 1092 "/home/runner/work/MFC/MFC/src/common/m_model.fpp"
1643 print *, 'm_model.fpp:1092: ', '@:ALLOCATE(gpu_total_vertices(1:num_ibs))'
1644# 1092 "/home/runner/work/MFC/MFC/src/common/m_model.fpp"
1645
1646# 1092 "/home/runner/work/MFC/MFC/src/common/m_model.fpp"
1647 call flush (output_unit)
1648# 1092 "/home/runner/work/MFC/MFC/src/common/m_model.fpp"
1649 end block
1650# 1092 "/home/runner/work/MFC/MFC/src/common/m_model.fpp"
1651#endif
1652# 1092 "/home/runner/work/MFC/MFC/src/common/m_model.fpp"
1653 allocate (gpu_total_vertices(1:num_ibs))
1654# 1092 "/home/runner/work/MFC/MFC/src/common/m_model.fpp"
1655
1656# 1092 "/home/runner/work/MFC/MFC/src/common/m_model.fpp"
1657
1658# 1092 "/home/runner/work/MFC/MFC/src/common/m_model.fpp"
1659#if defined(MFC_OpenACC)
1660# 1092 "/home/runner/work/MFC/MFC/src/common/m_model.fpp"
1661!$acc enter data create(gpu_total_vertices)
1662# 1092 "/home/runner/work/MFC/MFC/src/common/m_model.fpp"
1663#elif defined(MFC_OpenMP)
1664# 1092 "/home/runner/work/MFC/MFC/src/common/m_model.fpp"
1665!$omp target enter data map(always,alloc:gpu_total_vertices)
1666# 1092 "/home/runner/work/MFC/MFC/src/common/m_model.fpp"
1667#endif
1668
1669 gpu_ntrs = 0
1670 gpu_trs_v = 0._wp
1671 gpu_trs_n = 0._wp
1672 gpu_boundary_edge_count = 0
1673 gpu_total_vertices = 0
1674
1675 if (max_bv1 > 0) then
1676#ifdef MFC_DEBUG
1677# 1101 "/home/runner/work/MFC/MFC/src/common/m_model.fpp"
1678 block
1679# 1101 "/home/runner/work/MFC/MFC/src/common/m_model.fpp"
1680 use iso_fortran_env, only: output_unit
1681# 1101 "/home/runner/work/MFC/MFC/src/common/m_model.fpp"
1682
1683# 1101 "/home/runner/work/MFC/MFC/src/common/m_model.fpp"
1684 print *, 'm_model.fpp:1101: ', '@:ALLOCATE(gpu_boundary_v(1:max_bv1, 1:max_bv2, 1:max_bv3, 1:num_ibs))'
1685# 1101 "/home/runner/work/MFC/MFC/src/common/m_model.fpp"
1686
1687# 1101 "/home/runner/work/MFC/MFC/src/common/m_model.fpp"
1688 call flush (output_unit)
1689# 1101 "/home/runner/work/MFC/MFC/src/common/m_model.fpp"
1690 end block
1691# 1101 "/home/runner/work/MFC/MFC/src/common/m_model.fpp"
1692#endif
1693# 1101 "/home/runner/work/MFC/MFC/src/common/m_model.fpp"
1694 allocate (gpu_boundary_v(1:max_bv1, 1:max_bv2, 1:max_bv3, 1:num_ibs))
1695# 1101 "/home/runner/work/MFC/MFC/src/common/m_model.fpp"
1696
1697# 1101 "/home/runner/work/MFC/MFC/src/common/m_model.fpp"
1698
1699# 1101 "/home/runner/work/MFC/MFC/src/common/m_model.fpp"
1700#if defined(MFC_OpenACC)
1701# 1101 "/home/runner/work/MFC/MFC/src/common/m_model.fpp"
1702!$acc enter data create(gpu_boundary_v)
1703# 1101 "/home/runner/work/MFC/MFC/src/common/m_model.fpp"
1704#elif defined(MFC_OpenMP)
1705# 1101 "/home/runner/work/MFC/MFC/src/common/m_model.fpp"
1706!$omp target enter data map(always,alloc:gpu_boundary_v)
1707# 1101 "/home/runner/work/MFC/MFC/src/common/m_model.fpp"
1708#endif
1709 gpu_boundary_v = 0._wp
1710 end if
1711
1712 do pid = 1, num_ibs
1713 if (allocated(models(pid)%model)) then
1714 gpu_ntrs(pid) = models(pid)%ntrs
1715 gpu_trs_v(:,:,1:models(pid)%ntrs,pid) = models(pid)%trs_v
1716 gpu_trs_n(:,1:models(pid)%ntrs,pid) = models(pid)%trs_n
1717 gpu_boundary_edge_count(pid) = models(pid)%boundary_edge_count
1718 gpu_total_vertices(pid) = models(pid)%total_vertices
1719 end if
1720 if (allocated(models(pid)%boundary_v) .and. p == 0) then
1721 gpu_boundary_v(1:size(models(pid)%boundary_v, 1),1:size(models(pid)%boundary_v, 2), &
1722 & 1:size(models(pid)%boundary_v, 3),pid) = models(pid)%boundary_v
1723 end if
1724 end do
1725
1726
1727# 1119 "/home/runner/work/MFC/MFC/src/common/m_model.fpp"
1728#if defined(MFC_OpenACC)
1729# 1119 "/home/runner/work/MFC/MFC/src/common/m_model.fpp"
1730!$acc update device(gpu_ntrs, gpu_trs_v, gpu_trs_n, gpu_boundary_edge_count, gpu_total_vertices)
1731# 1119 "/home/runner/work/MFC/MFC/src/common/m_model.fpp"
1732#elif defined(MFC_OpenMP)
1733# 1119 "/home/runner/work/MFC/MFC/src/common/m_model.fpp"
1734!$omp target update to(gpu_ntrs, gpu_trs_v, gpu_trs_n, gpu_boundary_edge_count, gpu_total_vertices)
1735# 1119 "/home/runner/work/MFC/MFC/src/common/m_model.fpp"
1736#endif
1737 if (allocated(gpu_boundary_v)) then
1738
1739# 1121 "/home/runner/work/MFC/MFC/src/common/m_model.fpp"
1740#if defined(MFC_OpenACC)
1741# 1121 "/home/runner/work/MFC/MFC/src/common/m_model.fpp"
1742!$acc update device(gpu_boundary_v)
1743# 1121 "/home/runner/work/MFC/MFC/src/common/m_model.fpp"
1744#elif defined(MFC_OpenMP)
1745# 1121 "/home/runner/work/MFC/MFC/src/common/m_model.fpp"
1746!$omp target update to(gpu_boundary_v)
1747# 1121 "/home/runner/work/MFC/MFC/src/common/m_model.fpp"
1748#endif
1749 end if
1750 end if
1751 end block
1752
1753 end subroutine s_instantiate_stl_models
1754#endif
1755
1756 !> Pack triangle vertices and normals from a model into flat arrays for GPU transfer.
1757 subroutine s_pack_model_for_gpu(ma)
1758
1759 type(t_model_array), intent(inout) :: ma
1760 integer :: i
1761
1762 ma%ntrs = ma%model%ntrs
1763 allocate (ma%trs_v(1:3,1:3,1:ma%ntrs))
1764 allocate (ma%trs_n(1:3,1:ma%ntrs))
1765
1766 do i = 1, ma%ntrs
1767 ma%trs_v(:,:,i) = ma%model%trs(i)%v(:,:)
1768 ma%trs_n(:,i) = ma%model%trs(i)%n(:)
1769 end do
1770
1771 end subroutine s_pack_model_for_gpu
1772
1773end module m_model
k
integer, intent(in) k
Definition m_phase_change.fpp.f90:1064
j
integer, intent(in) j
Definition m_phase_change.fpp.f90:1064
l
integer, intent(in) l
Definition m_phase_change.fpp.f90:1064
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_helper
Utility routines for bubble model setup, coordinate transforms, array sampling, and special functions...
Definition m_helper.fpp.f90:328
m_helper::s_transform_model
subroutine, public s_transform_model(model, matrix, matrix_n)
Transform a model by a matrix, one triangle at a time.
Definition m_helper.fpp.f90:886
m_helper::f_create_bbox
type(t_bbox) function, public f_create_bbox(model)
Create a bounding box for a model.
Definition m_helper.fpp.f90:899
m_helper::f_create_transform_matrix
real(wp) function, dimension(1:4, 1:4), public f_create_transform_matrix(param, center)
Create a transformation matrix.
Definition m_helper.fpp.f90:821
m_helper::f_cross
pure real(wp) function, dimension(3), public f_cross(a, b)
Compute the cross product of two vectors.
Definition m_helper.fpp.f90:781
m_model
Binary STL file reader and processor for immersed boundary geometry.
Definition m_model.fpp.f90:318
m_model::s_model_write
impure subroutine, public s_model_write(filepath, model)
Write a mesh to a file.
Definition m_model.fpp.f90:708
m_model::s_write_stl
impure subroutine s_write_stl(filepath, model)
Write a binary STL file.
Definition m_model.fpp.f90:632
m_model::f_intersects_triangle
integer function f_intersects_triangle(ray, triangle)
Check if a ray intersects a triangle using the Moller-Trumbore algorithm (barycentric coordinates)....
Definition m_model.fpp.f90:936
m_model::f_read_line
impure logical function f_read_line(iunit, line)
Read the next non-blank, non-comment line from an STL or OBJ model file.
Definition m_model.fpp.f90:736
m_model::gpu_total_vertices
integer, dimension(:), allocatable, public gpu_total_vertices
Definition m_model.fpp.f90:346
m_model::s_read_obj
impure subroutine s_read_obj(filepath, model)
Read an OBJ file.
Definition m_model.fpp.f90:544
m_model::f_model_is_inside
impure real(wp) function, public f_model_is_inside(model, point, spacing, spc)
Determine whether a point is inside a model using stochastic ray casting.
Definition m_model.fpp.f90:803
m_model::f_model_is_inside_flat
real(wp) function, public f_model_is_inside_flat(ntrs, pid, point)
Determine if a point is inside a surface using the generalized winding number (Jacobson et al....
Definition m_model.fpp.f90:859
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::s_skip_ignored_lines
impure subroutine s_skip_ignored_lines(iunit, buffered_line, is_buffered)
Read the next non-comment line from a model file, using a buffered look-ahead mechanism.
Definition m_model.fpp.f90:765
m_model::gpu_trs_v
real(wp), dimension(:,:,:,:), allocatable, public gpu_trs_v
Definition m_model.fpp.f90:342
m_model::s_model_free
subroutine, public s_model_free(model)
Free the memory allocated for an STL mesh.
Definition m_model.fpp.f90:727
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::gpu_boundary_v
real(wp), dimension(:,:,:,:), allocatable, public gpu_boundary_v
Definition m_model.fpp.f90:344
m_model::s_register_edge
subroutine, public s_register_edge(temp_boundary_v, edge, edge_index, edge_count)
Append the edge end vertices to a temporary buffer.
Definition m_model.fpp.f90:1140
m_model::s_read_stl
impure subroutine s_read_stl(filepath, model)
Read an STL file.
Definition m_model.fpp.f90:516
m_model::f_model_read
impure type(t_model) function, public f_model_read(filepath)
Read a mesh from a file.
Definition m_model.fpp.f90:613
m_model::gpu_trs_n
real(wp), dimension(:,:,:), allocatable, public gpu_trs_n
Definition m_model.fpp.f90:343
m_model::s_check_boundary
subroutine, public s_check_boundary(model, boundary_v, boundary_vertex_count, boundary_edge_count)
Check and label edges shared by two or more triangle facets of the 2D STL model.
Definition m_model.fpp.f90:988
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::models
type(t_model_array), dimension(:), allocatable, target, public models
STL/OBJ models for IB markers and levelset.
Definition m_model.fpp.f90:340
m_model::s_write_obj
impure subroutine s_write_obj(filepath, model)
Write an OBJ file.
Definition m_model.fpp.f90:677
m_model::stl_bounding_boxes
real(wp), dimension(:,:,:), allocatable, public stl_bounding_boxes
Definition m_model.fpp.f90:347
m_model::s_read_stl_ascii
impure subroutine s_read_stl_ascii(filepath, model)
Read an ASCII STL file.
Definition m_model.fpp.f90:408
m_model::s_read_stl_binary
impure subroutine s_read_stl_binary(filepath, model)
Read a binary STL file.
Definition m_model.fpp.f90:364
m_model::s_instantiate_stl_models
subroutine, public s_instantiate_stl_models()
Load, transform, and register STL/OBJ immersed-boundary models onto the simulation grid.
Definition m_model.fpp.f90:1383
m_model::s_pack_model_for_gpu
subroutine, public s_pack_model_for_gpu(ma)
Pack triangle vertices and normals from a model into flat arrays for GPU transfer.
Definition m_model.fpp.f90:1758
m_model::gpu_boundary_edge_count
integer, dimension(:), allocatable, public gpu_boundary_edge_count
Definition m_model.fpp.f90:345
m_model::f_model_random_number
real(wp) function f_model_random_number(seed)
Generate a pseudo-random number using a seed-based xorshift, replacing the Fortran intrinsic which is...
Definition m_model.fpp.f90:786
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::ic_model_parameters
Defines parameters for a Model Patch.
Definition m_derived_types.fpp.f90:449
m_derived_types::t_bbox
Definition m_derived_types.fpp.f90:468
m_derived_types::t_model_array
Definition m_derived_types.fpp.f90:478
m_derived_types::t_model
Definition m_derived_types.fpp.f90:473
m_derived_types::t_ray
Definition m_derived_types.fpp.f90:463
m_derived_types::t_triangle
Definition m_derived_types.fpp.f90:458