Skip to content

Commit a1769d0

Browse files
danieljvickersdanieljvickers
committed
Moved mdoel-specific code to the model file for cleanliness
1 parent bc972ca commit a1769d0

4 files changed

Lines changed: 205 additions & 195 deletions

File tree

‎src/common/m_model.fpp‎

Lines changed: 181 additions & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -17,16 +17,196 @@ module m_model
1717

1818
private
1919

20-
public :: f_model_read, s_model_write, s_model_free, f_model_is_inside
20+
public :: s_instantiate_STL_models, f_model_read, s_model_write, s_model_free, f_model_is_inside, models, gpu_ntrs, gpu_trs_v, gpu_trs_n
2121

2222
! Subroutines for STL immersed boundaries
2323
public :: f_check_boundary, f_register_edge, f_check_interpolation_2D, &
2424
f_check_interpolation_3D, f_interpolate_2D, f_interpolate_3D, &
2525
f_interpolated_distance, f_normals, f_distance, f_distance_normals_3D, f_tri_area, s_pack_model_for_gpu, &
2626
f_model_is_inside_flat
2727

28+
!! array of STL models that can be allocated and then used in IB marker and levelset compute
29+
type(t_model_array), allocatable, target :: models(:)
30+
!! GPU-friendly flat arrays for STL model data
31+
integer, allocatable :: gpu_ntrs(:)
32+
real(wp), allocatable, dimension(:, :, :, :) :: gpu_trs_v
33+
real(wp), allocatable, dimension(:, :, :) :: gpu_trs_n
34+
real(wp), allocatable, dimension(:, :, :) :: gpu_boundary_v
35+
real(wp), allocatable, dimension(:, :, :) :: gpu_interpolated_boundary_v
36+
2837
contains
2938

39+
subroutine s_instantiate_STL_models()
40+
41+
! Variables for IBM+STL
42+
real(wp) :: normals(1:3) !< Boundary normal buffer
43+
integer :: boundary_vertex_count, boundary_edge_count, total_vertices !< Boundary vertex
44+
real(wp), allocatable, dimension(:, :, :) :: boundary_v !< Boundary vertex buffer
45+
real(wp), allocatable, dimension(:, :) :: interpolated_boundary_v !< Interpolated vertex buffer
46+
real(wp) :: dx_local, dy_local, dz_local !< Levelset distance buffer
47+
logical :: interpolate !< Logical variable to determine whether or not the model should be interpolated
48+
49+
integer :: i, j, k !< Generic loop iterators
50+
integer :: patch_id
51+
52+
type(t_bbox) :: bbox, bbox_old
53+
type(t_model) :: model
54+
type(ic_model_parameters) :: params
55+
56+
real(wp) :: eta
57+
real(wp), dimension(1:3) :: point, model_center
58+
real(wp) :: grid_mm(1:3, 1:2)
59+
60+
real(wp), dimension(1:4, 1:4) :: transform, transform_n
61+
62+
#ifdef MFC_PRE_PROCESS
63+
dx_local = dx; dy_local = dy
64+
if (p /= 0) dz_local = dz
65+
#else
66+
dx_local = minval(dx); dy_local = minval(dy)
67+
if (p /= 0) dz_local = minval(dz)
68+
#endif
69+
70+
do patch_id = 1, num_ibs
71+
if (patch_ib(patch_id)%geometry == 5 .or. patch_ib(patch_id)%geometry == 12) then
72+
allocate (models(patch_id)%model)
73+
print *, " * Reading model: "//trim(patch_ib(patch_id)%model_filepath)
74+
75+
model = f_model_read(patch_ib(patch_id)%model_filepath)
76+
params%scale(:) = patch_ib(patch_id)%model_scale(:)
77+
params%translate(:) = patch_ib(patch_id)%model_translate(:)
78+
params%rotate(:) = patch_ib(patch_id)%model_rotate(:)
79+
params%spc = patch_ib(patch_id)%model_spc
80+
params%threshold = patch_ib(patch_id)%model_threshold
81+
82+
if (f_approx_equal(dot_product(params%scale, params%scale), 0._wp)) then
83+
params%scale(:) = 1._wp
84+
end if
85+
86+
if (proc_rank == 0) then
87+
print *, " * Transforming model."
88+
end if
89+
90+
! Get the model center before transforming the model
91+
bbox_old = f_create_bbox(model)
92+
model_center(1:3) = (bbox_old%min(1:3) + bbox_old%max(1:3))/2._wp
93+
94+
! Compute the transform matrices for vertices and normals
95+
transform = f_create_transform_matrix(params, model_center)
96+
transform_n = f_create_transform_matrix(params)
97+
98+
call s_transform_model(model, transform, transform_n)
99+
100+
! Recreate the bounding box after transformation
101+
bbox = f_create_bbox(model)
102+
103+
! Show the number of vertices in the original STL model
104+
if (proc_rank == 0) then
105+
print *, ' * Number of input model vertices:', 3*model%ntrs
106+
end if
107+
108+
call f_check_boundary(model, boundary_v, boundary_vertex_count, boundary_edge_count)
109+
110+
! Check if the model needs interpolation
111+
if (p > 0) then
112+
call f_check_interpolation_3D(model, (/dx_local, dy_local, dz_local/), interpolate)
113+
else
114+
call f_check_interpolation_2D(boundary_v, boundary_edge_count, (/dx_local, dy_local, 0._wp/), interpolate)
115+
end if
116+
117+
! Show the number of edges and boundary edges in 2D STL models
118+
if (proc_rank == 0 .and. p == 0) then
119+
print *, ' * Number of 2D model boundary edges:', boundary_edge_count
120+
end if
121+
122+
! Interpolate the STL model along the edges (2D) and on triangle facets (3D)
123+
if (interpolate) then
124+
if (proc_rank == 0) then
125+
print *, ' * Interpolating STL vertices.'
126+
end if
127+
128+
if (p > 0) then
129+
call f_interpolate_3D(model, (/dx, dy, dz/), interpolated_boundary_v, total_vertices)
130+
else
131+
call f_interpolate_2D(boundary_v, boundary_edge_count, (/dx, dy, dz/), interpolated_boundary_v, total_vertices)
132+
end if
133+
134+
if (proc_rank == 0) then
135+
print *, ' * Total number of interpolated boundary vertices:', total_vertices
136+
end if
137+
end if
138+
139+
if (proc_rank == 0) then
140+
write (*, "(A, 3(2X, F20.10))") " > Model: Min:", bbox%min(1:3)
141+
write (*, "(A, 3(2X, F20.10))") " > Cen:", (bbox%min(1:3) + bbox%max(1:3))/2._wp
142+
write (*, "(A, 3(2X, F20.10))") " > Max:", bbox%max(1:3)
143+
144+
grid_mm(1, :) = (/minval(x_cc(0:m)) - 0.e5_wp*dx_local, maxval(x_cc(0:m)) + 0.e5_wp*dx_local/)
145+
grid_mm(2, :) = (/minval(y_cc(0:n)) - 0.e5_wp*dy_local, maxval(y_cc(0:n)) + 0.e5_wp*dy_local/)
146+
147+
if (p > 0) then
148+
grid_mm(3, :) = (/minval(z_cc(0:p)) - 0.e5_wp*dz_local, maxval(z_cc(0:p)) + 0.e5_wp*dz_local/)
149+
else
150+
grid_mm(3, :) = (/0._wp, 0._wp/)
151+
end if
152+
153+
write (*, "(A, 3(2X, F20.10))") " > Domain: Min:", grid_mm(:, 1)
154+
write (*, "(A, 3(2X, F20.10))") " > Cen:", (grid_mm(:, 1) + grid_mm(:, 2))/2._wp
155+
write (*, "(A, 3(2X, F20.10))") " > Max:", grid_mm(:, 2)
156+
end if
157+
158+
models(patch_id)%model = model
159+
models(patch_id)%boundary_v = boundary_v
160+
models(patch_id)%boundary_edge_count = boundary_edge_count
161+
if (interpolate) then
162+
models(patch_id)%interpolate = 1
163+
else
164+
models(patch_id)%interpolate = 0
165+
end if
166+
if (interpolate) then
167+
models(patch_id)%interpolated_boundary_v = interpolated_boundary_v
168+
models(patch_id)%total_vertices = total_vertices
169+
end if
170+
171+
end if
172+
end do
173+
174+
! Pack and upload flat arrays for GPU (AFTER the loop)
175+
block
176+
integer :: pid, max_ntrs
177+
178+
max_ntrs = 0
179+
do pid = 1, num_ibs
180+
if (allocated(models(pid)%model)) then
181+
call s_pack_model_for_gpu(models(pid))
182+
max_ntrs = max(max_ntrs, models(pid)%ntrs)
183+
end if
184+
end do
185+
186+
if (max_ntrs > 0) then
187+
allocate (gpu_ntrs(1:num_ibs))
188+
allocate (gpu_trs_v(1:3, 1:3, 1:max_ntrs, 1:num_ibs))
189+
allocate (gpu_trs_n(1:3, 1:max_ntrs, 1:num_ibs))
190+
191+
gpu_ntrs = 0
192+
gpu_trs_v = 0._wp
193+
gpu_trs_n = 0._wp
194+
195+
do pid = 1, num_ibs
196+
if (allocated(models(pid)%model)) then
197+
gpu_ntrs(pid) = models(pid)%ntrs
198+
gpu_trs_v(:, :, 1:models(pid)%ntrs, pid) = models(pid)%trs_v
199+
gpu_trs_n(:, 1:models(pid)%ntrs, pid) = models(pid)%trs_n
200+
end if
201+
end do
202+
203+
$:GPU_ENTER_DATA(copyin='[gpu_ntrs, gpu_trs_v, gpu_trs_n]')
204+
205+
end if
206+
end block
207+
208+
end subroutine s_instantiate_STL_models
209+
30210
!> This procedure reads a binary STL file.
31211
!! @param filepath Path to the STL file.
32212
!! @param model The binary of the STL file.

0 commit comments

Comments
 (0)