Skip to content

Commit 6cc7acc

Browse files
danieljvickersdanieljvickers
committed
STLs ran on GPU in 3D!
1 parent 5e58655 commit 6cc7acc

2 files changed

Lines changed: 244 additions & 226 deletions

File tree

‎src/common/m_model.fpp‎

Lines changed: 3 additions & 219 deletions
Original file line numberDiff line numberDiff line change
@@ -17,9 +17,9 @@ module m_model
1717

1818
private
1919

20-
public :: s_instantiate_STL_models, f_model_read, s_model_write, s_model_free, f_model_is_inside, models, gpu_ntrs, &
20+
public :: f_model_read, s_model_write, s_model_free, f_model_is_inside, models, gpu_ntrs, &
2121
gpu_trs_v, gpu_trs_n, gpu_boundary_v, gpu_interpolated_boundary_v, gpu_interpolate, gpu_boundary_edge_count, &
22-
gpu_total_vertices
22+
gpu_total_vertices, stl_bounding_boxes
2323

2424
! Subroutines for STL immersed boundaries
2525
public :: f_check_boundary, f_register_edge, f_check_interpolation_2D, &
@@ -38,226 +38,10 @@ module m_model
3838
integer, allocatable :: gpu_interpolate(:)
3939
integer, allocatable :: gpu_boundary_edge_count(:)
4040
integer, allocatable :: gpu_total_vertices(:)
41+
real(wp), allocatable :: stl_bounding_boxes(:, :, :)
4142

4243
contains
4344

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

0 commit comments

Comments
 (0)