Neko 1.99.1
A portable framework for high-order spectral element flow simulations
Loading...
Searching...
No Matches
facet_normal.f90
Go to the documentation of this file.
1! Copyright (c) 2020-2025, The Neko Authors
2! All rights reserved.
3!
4! Redistribution and use in source and binary forms, with or without
5! modification, are permitted provided that the following conditions
6! are met:
7!
8! * Redistributions of source code must retain the above copyright
9! notice, this list of conditions and the following disclaimer.
10!
11! * Redistributions in binary form must reproduce the above
12! copyright notice, this list of conditions and the following
13! disclaimer in the documentation and/or other materials provided
14! with the distribution.
15!
16! * Neither the name of the authors nor the names of its
17! contributors may be used to endorse or promote products derived
18! from this software without specific prior written permission.
19!
20! THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
21! "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
22! LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS
23! FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE
24! COPYRIGHT OWNER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT,
25! INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING,
26! BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES;
27! LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER
28! CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT
29! LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN
30! ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
31! POSSIBILITY OF SUCH DAMAGE.
32!
36 use num_types, only : rp
38 use math, only: cfill_mask
41 use vector, only : vector_t
42 use coefs, only : coef_t
43 use bc, only : bc_t
45 use json_module, only : json_file
46 use, intrinsic :: iso_c_binding, only : c_ptr, c_null_ptr, c_associated
47 use htable, only : htable_i4_t
50 use time_state, only : time_state_t
51 implicit none
52 private
53
55 type, public, extends(bc_t) :: facet_normal_t
56 integer, allocatable :: unique_mask(:)
57 type(c_ptr) :: unique_mask_d = c_null_ptr
58 type(vector_t) :: nx, ny, nz, work
59 contains
60 procedure, pass(this) :: apply_scalar => facet_normal_apply_scalar
61 procedure, pass(this) :: apply_scalar_dev => facet_normal_apply_scalar_dev
62 procedure, pass(this) :: apply_vector => facet_normal_apply_vector
63 procedure, pass(this) :: apply_vector_dev => facet_normal_apply_vector_dev
64 procedure, pass(this) :: apply_surfvec => facet_normal_apply_surfvec
65 procedure, pass(this) :: apply_surfvec_dev => &
68 procedure, pass(this) :: init => facet_normal_init
70 procedure, pass(this) :: init_from_components => &
73 procedure, pass(this) :: free => facet_normal_free
75 procedure, pass(this) :: finalize => facet_normal_finalize
76 end type facet_normal_t
77
78contains
79
83 subroutine facet_normal_init(this, coef, json)
84 class(facet_normal_t), intent(inout), target :: this
85 type(coef_t), target, intent(in) :: coef
86 type(json_file), intent(inout) ::json
87
88 call this%init_from_components(coef)
89 end subroutine facet_normal_init
90
93 subroutine facet_normal_init_from_components(this, coef)
94 class(facet_normal_t), intent(inout), target :: this
95 type(coef_t), target, intent(in) :: coef
96
97 call this%init_base(coef)
99
101 subroutine facet_normal_apply_scalar(this, x, n, time, strong)
102 class(facet_normal_t), intent(inout) :: this
103 integer, intent(in) :: n
104 real(kind=rp), intent(inout), dimension(n) :: x
105 type(time_state_t), intent(in), optional :: time
106 logical, intent(in), optional :: strong
107 end subroutine facet_normal_apply_scalar
108
110 subroutine facet_normal_apply_scalar_dev(this, x_d, time, strong, strm)
111 class(facet_normal_t), intent(inout), target :: this
112 type(c_ptr),intent(inout) :: x_d
113 type(time_state_t), intent(in), optional :: time
114 logical, intent(in), optional :: strong
115 type(c_ptr), intent(inout) :: strm
116
117 end subroutine facet_normal_apply_scalar_dev
118
120 subroutine facet_normal_apply_vector_dev(this, x_d, y_d, z_d, time, &
121 strong, strm)
122 class(facet_normal_t), intent(inout), target :: this
123 type(c_ptr), intent(inout) :: x_d
124 type(c_ptr), intent(inout) :: y_d
125 type(c_ptr), intent(inout) :: z_d
126 type(time_state_t), intent(in), optional :: time
127 logical, intent(in), optional :: strong
128 type(c_ptr), intent(inout) :: strm
129
130 end subroutine facet_normal_apply_vector_dev
131
133 subroutine facet_normal_apply_vector(this, x, y, z, n, time, strong)
134 class(facet_normal_t), intent(inout) :: this
135 integer, intent(in) :: n
136 real(kind=rp), intent(inout), dimension(n) :: x
137 real(kind=rp), intent(inout), dimension(n) :: y
138 real(kind=rp), intent(inout), dimension(n) :: z
139 type(time_state_t), intent(in), optional :: time
140 logical, intent(in), optional :: strong
141 end subroutine facet_normal_apply_vector
142
144 subroutine facet_normal_apply_surfvec(this, x, y, z, u, v, w, n, time)
145 class(facet_normal_t), intent(in) :: this
146 integer, intent(in) :: n
147 real(kind=rp), intent(inout), dimension(n) :: x
148 real(kind=rp), intent(inout), dimension(n) :: y
149 real(kind=rp), intent(inout), dimension(n) :: z
150 real(kind=rp), intent(inout), dimension(n) :: u
151 real(kind=rp), intent(inout), dimension(n) :: v
152 real(kind=rp), intent(inout), dimension(n) :: w
153 type(time_state_t), intent(in), optional :: time
154 integer :: i, m, k, idx(4), facet
155 real(kind=rp) :: normal(3), area
156
157 m = this%unique_mask(0)
158
159 do i = 1, m
160 k = this%unique_mask(i)
161 x(k) = u(k) * this%nx%x(i)
162 y(k) = v(k) * this%ny%x(i)
163 z(k) = w(k) * this%nz%x(i)
164 end do
165
166 end subroutine facet_normal_apply_surfvec
167
169 subroutine facet_normal_apply_surfvec_dev(this, x_d, y_d, z_d, &
170 u_d, v_d, w_d, time, strm)
171 class(facet_normal_t), intent(in), target :: this
172 type(c_ptr) :: x_d, y_d, z_d, u_d, v_d, w_d
173 type(time_state_t), intent(in), optional :: time
174 type(c_ptr), optional :: strm
175 type(c_ptr) :: strm_
176 integer :: n, m
177
178 n = this%coef%dof%size()
179 m = this%unique_mask(0)
180
181 if (present(strm)) then
182 strm_ = strm
183 else
184 strm_ = glb_cmd_queue
185 end if
186
187 if (m .gt. 0) then
188 call device_masked_gather_copy_0(this%work%x_d, u_d, this%unique_mask_d, &
189 n, m, strm_)
190 call device_col2(this%work%x_d, this%nx%x_d, m, strm_)
191 call device_masked_scatter_copy_0(x_d, this%work%x_d, &
192 this%unique_mask_d, n, m, strm_)
193 call device_masked_gather_copy_0(this%work%x_d, v_d, this%unique_mask_d, &
194 n , m, strm_)
195 call device_col2(this%work%x_d, this%ny%x_d, m, strm_)
196 call device_masked_scatter_copy_0(y_d, this%work%x_d, &
197 this%unique_mask_d, n, m, strm_)
198 call device_masked_gather_copy_0(this%work%x_d, w_d, this%unique_mask_d, &
199 n, m, strm_)
200 call device_col2(this%work%x_d, this%nz%x_d, m, strm)
201 call device_masked_scatter_copy_0(z_d, this%work%x_d, &
202 this%unique_mask_d, n, m, strm_)
203 end if
204
205 end subroutine facet_normal_apply_surfvec_dev
206
208 subroutine facet_normal_free(this)
209 class(facet_normal_t), target, intent(inout) :: this
210
211 call this%free_base()
212 if (allocated(this%unique_mask)) then
213 deallocate(this%unique_mask)
214 end if
215 if (c_associated(this%unique_mask_d)) then
216 call device_free(this%unique_mask_d)
217 end if
218
219 call this%nx%free()
220 call this%ny%free()
221 call this%nz%free()
222 call this%work%free()
223
224 end subroutine facet_normal_free
225
227 subroutine facet_normal_finalize(this, only_facets)
228 class(facet_normal_t), target, intent(inout) :: this
229 logical, optional, intent(in) :: only_facets
230 logical :: only_facets_
231 type(htable_i4_t) :: unique_point_idx
232 integer :: htable_data, rcode, i, j, idx(4), facet
233 real(kind=rp) :: area, normal(3)
234
235 if (present(only_facets)) then
236 if (only_facets .eqv. .false.) then
237 call neko_error("For facet_normal_t, only_facets has to be true.")
238 end if
239 end if
240
241 call this%finalize_base(.true.)
242 ! This part is purely needed to ensure that contributions
243 ! for all faces a point is on is properly summed up.
244 ! If one simply uses the original mask, if a point is on a corner
245 ! where both faces are on the boundary
246 ! one will only get the contribution from one face, not both
247 ! We solve this by adding up the normals of both faces for these points
248 ! and storing this sum in this%nx, this%ny, this%nz.
249 ! As both contrbutions are added already,
250 ! we also ensure that we only visit each point once
251 ! and create a new mask with only unique points (this%unique_mask).
252 if (allocated(this%unique_mask)) then
253 deallocate(this%unique_mask)
254 end if
255 if (c_associated(this%unique_mask_d)) then
256 call device_free(this%unique_mask_d)
257 end if
258
259 call unique_point_idx%init(this%msk(0), htable_data)
260 j = 0
261 do i = 1, this%msk(0)
262 if (unique_point_idx%get(this%msk(i),htable_data) .ne. 0) then
263 j = j + 1
264 htable_data = j
265 call unique_point_idx%set(this%msk(i), j)
266 end if
267 end do
268
269 ! Only allocate work vectors if size is non-zero
270 if (unique_point_idx%num_entries() .gt. 0 ) then
271 call this%nx%init(unique_point_idx%num_entries())
272 call this%ny%init(unique_point_idx%num_entries())
273 call this%nz%init(unique_point_idx%num_entries())
274 call this%work%init(unique_point_idx%num_entries())
275 end if
276 allocate(this%unique_mask(0:unique_point_idx%num_entries()))
277
278 this%unique_mask(0) = unique_point_idx%num_entries()
279 do i = 1, this%unique_mask(0)
280 this%unique_mask(i) = 0
281 end do
282
283
284 do i = 1, this%msk(0)
285 rcode = unique_point_idx%get(this%msk(i), htable_data)
286 if (rcode .ne. 0) call neko_error("Facet normal: htable get failed.")
287 this%unique_mask(htable_data) = this%msk(i)
288 facet = this%facet(i)
289
290 idx = nonlinear_index(this%msk(i), this%Xh%lx, this%Xh%lx, this%Xh%lx)
291 normal = this%coef%get_normal(idx(1), idx(2), idx(3), idx(4), facet)
292 area = this%coef%get_area(idx(1), idx(2), idx(3), idx(4), facet)
293 normal = normal * area !Scale normal by area
294 this%nx%x(htable_data) = this%nx%x(htable_data) + normal(1)
295 this%ny%x(htable_data) = this%ny%x(htable_data) + normal(2)
296 this%nz%x(htable_data) = this%nz%x(htable_data) + normal(3)
297 end do
298
299 if (neko_bcknd_device .eq. 1 .and. &
300 (unique_point_idx%num_entries() .gt. 0 )) then
301 call device_map(this%unique_mask, this%unique_mask_d, &
302 size(this%unique_mask))
303 call device_memcpy(this%unique_mask, this%unique_mask_d, &
304 size(this%unique_mask), host_to_device, sync = .true.)
305 call device_memcpy(this%nx%x, this%nx%x_d, &
306 this%nx%size(), host_to_device, sync = .true.)
307 call device_memcpy(this%ny%x, this%ny%x_d, &
308 this%ny%size(), host_to_device, sync = .true.)
309 call device_memcpy(this%nz%x, this%nz%x_d, &
310 this%nz%size(), host_to_device, sync = .true.)
311 end if
312
313 call unique_point_idx%free()
314
315 end subroutine facet_normal_finalize
316
317end module facet_normal
__inline__ __device__ void nonlinear_index(const int idx, const int lx, int *index)
Definition bc_utils.h:44
Map a Fortran array to a device (allocate and associate)
Definition device.F90:72
Copy data between host and device (or device and device)
Definition device.F90:66
Defines a boundary condition.
Definition bc.f90:34
Coefficients.
Definition coef.f90:34
subroutine, public device_masked_scatter_copy_0(a_d, b_d, mask_d, n, n_mask, strm)
Scatter a masked vector .
subroutine, public device_col2(a_d, b_d, n, strm)
Vector multiplication .
subroutine, public device_masked_gather_copy_0(a_d, b_d, mask_d, n, n_mask, strm)
Gather a masked vector .
Device abstraction, common interface for various accelerators.
Definition device.F90:34
integer, parameter, public host_to_device
Definition device.F90:47
subroutine, public device_free(x_d)
Deallocate memory on the device.
Definition device.F90:214
integer, parameter, public device_to_host
Definition device.F90:47
type(c_ptr), bind(C), public glb_cmd_queue
Global command queue.
Definition device.F90:51
Dirichlet condition applied in the facet normal direction.
subroutine facet_normal_init_from_components(this, coef)
Constructor from components.
subroutine facet_normal_apply_vector_dev(this, x_d, y_d, z_d, time, strong, strm)
No-op vector apply on device.
subroutine facet_normal_finalize(this, only_facets)
Finalize.
subroutine facet_normal_apply_scalar_dev(this, x_d, time, strong, strm)
No-op scalar apply on device.
subroutine facet_normal_init(this, coef, json)
Constructor.
subroutine facet_normal_apply_scalar(this, x, n, time, strong)
No-op scalar apply.
subroutine facet_normal_apply_surfvec_dev(this, x_d, y_d, z_d, u_d, v_d, w_d, time, strm)
Apply in facet normal direction (vector valued, device version)
subroutine facet_normal_apply_surfvec(this, x, y, z, u, v, w, n, time)
Apply in facet normal direction (vector valued)
subroutine facet_normal_apply_vector(this, x, y, z, n, time, strong)
No-op vector apply.
subroutine facet_normal_free(this)
Destructor.
Implements a hash table ADT.
Definition htable.f90:36
Definition math.f90:60
subroutine, public cfill_mask(a, c, n, mask, n_mask)
Fill a constant to a masked vector. .
Definition math.f90:403
Build configurations.
integer, parameter neko_bcknd_device
integer, parameter, public rp
Global precision used in computations.
Definition num_types.f90:12
Module with things related to the simulation time.
Utilities.
Definition utils.f90:35
Defines a vector.
Definition vector.f90:34
Base type for a boundary condition.
Definition bc.f90:61
Coefficients defined on a given (mesh, ) tuple. Arrays use indices (i,j,k,e): element e,...
Definition coef.f90:55
Dirichlet condition in facet normal direction.
Integer based hash table.
Definition htable.f90:82
A struct that contains all info about the time, expand as needed.