Neko 0.9.99
A portable framework for high-order spectral element flow simulations
Loading...
Searching...
No Matches
global_interpolation.F90
Go to the documentation of this file.
1! Copyright (c) 2020-2023, 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!
40 use num_types, only : rp
41 use space, only : space_t
42 use dofmap, only : dofmap_t
43 use mesh, only : mesh_t
44 use logger, only : neko_log, log_size
45 use utils, only : neko_error, neko_warning
47 use device
48 use comm, only : neko_comm, pe_size, pe_rank, mpi_gather, &
49 mpi_double_precision, mpi_integer, mpi_sum, mpi_reduce, mpi_gatherv
50 use math, only : copy
52 use structs, only : array_ptr_t
53 use, intrinsic :: iso_c_binding
54 implicit none
55 private
56
58 type, public :: global_interpolation_t
60 type(array_ptr_t) :: x
62 type(array_ptr_t) :: y
64 type(array_ptr_t) :: z
66 integer :: gdim
68 integer :: nelv
70 type(space_t), pointer :: xh
72 type(local_interpolator_t) :: local_interp
74 logical :: all_points_local = .false.
75 !! Gslib handle.
76 !! @note: Remove when we remove gslib.
77 integer :: gs_handle
78 logical :: gs_init = .false.
81 integer :: n_points
84 real(kind=rp), allocatable :: xyz(:,:)
86 integer, allocatable :: proc_owner(:)
88 integer, allocatable :: el_owner(:)
89 type(c_ptr) :: el_owner_d = c_null_ptr
92 real(kind=rp), allocatable :: rst(:,:)
95 real(kind=rp), allocatable :: dist2(:)
97 integer, allocatable :: error_code(:)
99 real(kind=rp) :: tol = 5d-13
100 contains
102 procedure, pass(this) :: init_xyz => global_interpolation_init_xyz
103 procedure, pass(this) :: init_dof => global_interpolation_init_dof
105 procedure, pass(this) :: free => global_interpolation_free
107 procedure, pass(this) :: free_points => global_interpolation_free_points
108 procedure, pass(this) :: find_points_and_redist => &
113 procedure, pass(this) :: find_points_coords => &
115 procedure, pass(this) :: find_points_xyz => global_interpolation_find_xyz
116 generic :: find_points => find_points_xyz, find_points_coords
118 procedure, pass(this) :: evaluate => global_interpolation_evaluate
119
121 generic :: init => init_dof, init_xyz
122
124
125contains
126
130 subroutine global_interpolation_init_dof(this, dof, tol)
131 class(global_interpolation_t), intent(inout) :: this
132 type(dofmap_t), target :: dof
133 real(kind=rp), optional :: tol
134
135 ! NOTE: Passing dof%x(:,1,1,1), etc in init_xyz passes down the entire
136 ! dof%x array and not a slice. It is done this way for
137 ! this%x%ptr to point to dof%x (see global_interpolation_init_xyz).
138 call this%init_xyz(dof%x(:,1,1,1), dof%y(:,1,1,1), dof%z(:,1,1,1), &
139 dof%msh%gdim, dof%msh%nelv, dof%Xh, tol = tol)
140
141 end subroutine global_interpolation_init_dof
142
152 subroutine global_interpolation_init_xyz(this, x, y, z, gdim, nelv, Xh, tol)
153 class(global_interpolation_t), intent(inout) :: this
154 real(kind=rp), intent(in), target :: x(:)
155 real(kind=rp), intent(in), target :: y(:)
156 real(kind=rp), intent(in), target :: z(:)
157 integer, intent(in) :: gdim
158 integer, intent(in) :: nelv
159 type(space_t), intent(in), target :: Xh
160 real(kind=rp), intent(in), optional :: tol
161 integer :: lx, ly, lz, max_pts_per_iter
162
163#ifdef HAVE_GSLIB
164
165 this%x%ptr => x
166 this%y%ptr => y
167 this%z%ptr => z
168 this%gdim = gdim
169 this%nelv = nelv
170 this%Xh => xh
171 if (present(tol)) this%tol = tol
172
173 ! Number of points to iterate on simultaneosuly
174 max_pts_per_iter = 128
175 lx = xh%lx
176 ly = xh%ly
177 lz = xh%lz
178
179 call fgslib_findpts_setup(this%gs_handle, &
181 this%gdim, &
182 this%x%ptr, this%y%ptr, this%z%ptr, & ! Physical nodal values
183 lx, ly, lz, nelv, & ! Mesh dimensions
184 2*lx, 2*ly, 2*lz, & ! Mesh size for bounding box computation
185 0.01, & ! relative size to expand bounding boxes by
186 lx*ly*lz*nelv, lx*ly*lz*nelv, & ! local/global hash mesh sizes
187 max_pts_per_iter, this%tol)
188 this%gs_init = .true.
189#else
190 call neko_error('Neko needs to be built with GSLIB support')
191#endif
192
193 end subroutine global_interpolation_init_xyz
194
195
198 class(global_interpolation_t), intent(inout) :: this
199
200 nullify(this%x%ptr)
201 nullify(this%y%ptr)
202 nullify(this%z%ptr)
203 nullify(this%Xh)
204
205 this%nelv = 0
206 this%gdim = 0
207
208 call this%free_points()
209
210#ifdef HAVE_GSLIB
211 if (this%gs_init) then
212 call fgslib_findpts_free(this%gs_handle)
213 this%gs_init = .false.
214 end if
215#endif
216
217 end subroutine global_interpolation_free
218
221 class(global_interpolation_t), intent(inout) :: this
222
223 this%n_points = 0
224 this%all_points_local = .false.
225
226 if (allocated(this%xyz)) deallocate(this%xyz)
227 if (allocated(this%rst)) deallocate(this%rst)
228 if (allocated(this%proc_owner)) deallocate(this%proc_owner)
229 if (allocated(this%el_owner)) deallocate(this%el_owner)
230 if (allocated(this%dist2)) deallocate(this%dist2)
231 if (allocated(this%error_code)) deallocate(this%error_code)
232
233 if (c_associated(this%el_owner_d)) then
234 call device_free(this%el_owner_d)
235 end if
236
238
241 class(global_interpolation_t), intent(inout) :: this
242 !!Perhaps this should be kind dp
243 real(kind=rp) :: xdiff, ydiff, zdiff
244 character(len=8000) :: log_buf
245 real(kind=rp), allocatable :: x_check(:)
246 real(kind=rp), allocatable :: y_check(:)
247 real(kind=rp), allocatable :: z_check(:)
248 logical :: isdiff
249 integer :: i
250
251
252#ifdef HAVE_GSLIB
253
254 ! gslib find points, which element they belong, to process etc.
255 call fgslib_findpts(this%gs_handle, &
256 this%error_code, 1, &
257 this%proc_owner, 1, &
258 this%el_owner, 1, &
259 this%rst, this%gdim, &
260 this%dist2, 1, &
261 this%xyz(1,1), this%gdim, &
262 this%xyz(2,1), this%gdim, &
263 this%xyz(3,1), this%gdim, this%n_points)
264
265 do i = 1 , this%n_points
266
267 !
268 ! Check validity of points
269 !
270 if (this%error_code(i) .eq. 1) then
271 if (this%dist2(i) .gt. this%tol) then
272 write(*,*) 'Point with coords: ',&
273 this%xyz(1,i),&
274 this%xyz(2,i),&
275 this%xyz(3,i),&
276 'Did not converge to tol. Absolute differences squared: ',&
277 this%dist2(i), 'PE rank', pe_rank
278 end if
279 end if
280
281 if (this%error_code(i) .eq. 2) &
282 write(*,*) 'Point with coords: ',&
283 this%xyz(1,i), this%xyz(2,i), this%xyz(3,i),&
284 'Outside the mesh!',&
285 ' Interpolation on these points will return 0.0. dist2: ', &
286 this%dist2(i),&
287 'el_owner, rst coords, pe: ',&
288 this%el_owner(i), this%rst(1,i), this%rst(2,i), &
289 this%rst(3,i), pe_rank
290
291 end do
292
293 allocate(x_check(this%n_points))
294 allocate(y_check(this%n_points))
295 allocate(z_check(this%n_points))
296
297 call fgslib_findpts_eval(this%gs_handle, x_check, &
298 1, this%error_code, 1, &
299 this%proc_owner, 1, this%el_owner, 1, &
300 this%rst, this%gdim, &
301 this%n_points, this%x%ptr)
302
303 call fgslib_findpts_eval(this%gs_handle, y_check, &
304 1, this%error_code, 1, &
305 this%proc_owner, 1, this%el_owner, 1, &
306 this%rst, this%gdim, &
307 this%n_points, this%y%ptr)
308
309 call fgslib_findpts_eval(this%gs_handle, z_check, &
310 1, this%error_code, 1, &
311 this%proc_owner, 1, this%el_owner, 1, &
312 this%rst, this%gdim, &
313 this%n_points, this%z%ptr)
314
315 do i = 1 , this%n_points
316
317 !
318 ! Check validity of points
319 !
320 isdiff = .false.
321 xdiff = (x_check(i)-this%xyz(1,i))**2
322 if ( xdiff .gt. this%tol) isdiff = .true.
323 ydiff = (y_check(i)-this%xyz(2,i))**2
324 if ( ydiff .gt. this%tol) isdiff = .true.
325 zdiff = (z_check(i)-this%xyz(3,i))**2
326 if ( zdiff .gt. this%tol) isdiff = .true.
327
328 if (isdiff) then
329 write(*,*) 'Points with coords: ', &
330 this%xyz(1, i), this%xyz(2, i), this%xyz(3, i), &
331 'Differ from interpolated coords: ', &
332 x_check(i), y_check(i), z_check(i), &
333 'Distance squared: ', &
334 xdiff, ydiff, zdiff
335 end if
336
337 end do
338
339 deallocate(x_check)
340 deallocate(y_check)
341 deallocate(z_check)
342
343 if (neko_bcknd_device .eq. 1) then
344 call device_memcpy(this%el_owner, this%el_owner_d, &
345 this%n_points, host_to_device, sync = .true.)
346 end if
347#else
348 call neko_error('Neko needs to be built with GSLIB support')
349#endif
351
364 subroutine global_interpolation_find_coords(this, x, y, z, n_points)
365 class(global_interpolation_t), intent(inout) :: this
366 integer :: n_points
367 !!Perhaps this should be kind dp
368 !!this is to get around that x,y,z often is 4 dimensional...
369 !!Should maybe add interface for 1d aswell
370 real(kind=rp) :: x(n_points,1,1,1)
371 real(kind=rp) :: y(n_points,1,1,1)
372 real(kind=rp) :: z(n_points,1,1,1)
373 integer :: i
374
375 call this%free_points()
376
377 this%n_points = n_points
378
380
381 do i = 1, n_points
382 this%xyz(1, i) = x(i,1,1,1)
383 this%xyz(2, i) = y(i,1,1,1)
384 this%xyz(3, i) = z(i,1,1,1)
385 end do
386
388
390
392 class(global_interpolation_t) :: this
393
394 allocate(this%xyz(3, this%n_points))
395 allocate(this%rst(3, this%n_points))
396 allocate(this%proc_owner(this%n_points))
397 allocate(this%el_owner(this%n_points))
398 allocate(this%dist2(this%n_points))
399 allocate(this%error_code(this%n_points))
400
401 if (neko_bcknd_device .eq. 1) &
402 call device_map(this%el_owner, this%el_owner_d, this%n_points)
403
405
416 subroutine global_interpolation_find_xyz(this, xyz, n_points)
417 class(global_interpolation_t), intent(inout) :: this
418 integer, intent(in) :: n_points
419 !!Perhaps this should be kind dp
420 real(kind=rp), intent(inout) :: xyz(3, n_points)
421
422
423 call this%free_points()
424
425 this%n_points = n_points
426
428
430 call copy(this%xyz, xyz, 3 * n_points)
431
433
434 end subroutine global_interpolation_find_xyz
435
447 subroutine global_interpolation_find_and_redist(this, xyz, n_points)
448 class(global_interpolation_t), intent(inout) :: this
449 integer, intent(inout) :: n_points
450 !!Perhaps this should be kind dp
451 real(kind=rp), allocatable, intent(inout) :: xyz(:,:)
452 integer :: i
453
454
455 call this%free_points()
456
457 this%n_points = n_points
459
461 call copy(this%xyz, xyz, 3 * n_points)
462
467
468 do i = 1, this%n_points
469 if (this%proc_owner(i) .ne. pe_rank) then
470 write(*,*) 'Redistribution failed on rank: ', pe_rank, &
471 'for point with coord: ', &
472 this%xyz(1, i), this%xyz(2, i), this%xyz(3, i)
473 exit
474 end if
475 end do
476
477 n_points = this%n_points
478 if (allocated(xyz)) then
479 deallocate(xyz)
480 end if
481 allocate(xyz(3, n_points))
482 call copy(xyz, this%xyz, 3*n_points)
483
484 call this%local_interp%init(this%Xh, this%rst(1,:),&
485 this%rst(2,:), this%rst(3,:), n_points)
486 this%all_points_local = .true.
487
488
490
492 class(global_interpolation_t), intent(inout) :: this
493 integer, allocatable :: n_points_per_pe(:)
494 integer, allocatable :: n_points_from_pe(:)
495 integer, allocatable :: n_point_offset_from_pe(:)
496 real(kind=rp), allocatable :: xyz_send_to_pe(:,:)
497 real(kind=rp), allocatable :: new_xyz(:,:)
498 integer :: i, j, k, ierr, n_new_points, max_n_points_to_send
499
500 n_new_points = 0
501
502 allocate(n_points_per_pe(0:(pe_size-1)))
503 allocate(n_points_from_pe(0:(pe_size-1)))
504 n_points_per_pe = 0
505 n_points_from_pe = 0
507 do i = 1, this%n_points
508 n_points_per_pe(this%proc_owner(i)) = &
509 n_points_per_pe(this%proc_owner(i)) + 1
510 end do
513 do i = 0, (pe_size - 1)
514 call mpi_reduce(n_points_per_pe(i), n_new_points, 1, mpi_integer, &
515 mpi_sum, i, neko_comm, ierr)
516 call mpi_gather(n_points_per_pe(i), 1, mpi_integer,&
517 n_points_from_pe, 1, mpi_integer, i, neko_comm, ierr)
518 end do
519
520 allocate(n_point_offset_from_pe(0:(pe_size-1)))
521 n_point_offset_from_pe(0) = 0
522 do i = 1, (pe_size - 1)
523 n_point_offset_from_pe(i) = n_points_from_pe(i-1)&
524 + n_point_offset_from_pe(i-1)
525 end do
526
527 allocate(new_xyz(3, n_new_points))
528 max_n_points_to_send = maxval(n_points_per_pe)
529 allocate(xyz_send_to_pe(3, max_n_points_to_send))
530 do i = 0, (pe_size - 1)
532 k = 0
533 do j = 1, this%n_points
534 if (this%proc_owner(j) .eq. i) then
535 k = k + 1
536 xyz_send_to_pe(:,k) = this%xyz(:,j)
537 end if
538 end do
539 if (k .ne. n_points_per_pe(i)) then
540 write(*,*) 'PE: ', pe_rank, ' has k= ', k, &
541 'points for PE:', i,' but should have: ', &
542 n_points_per_pe(i)
543 call neko_error('Error in redistribution of points')
544 end if
545 call mpi_gatherv(xyz_send_to_pe,3*n_points_per_pe(i), &
546 mpi_double_precision, new_xyz,3*n_points_from_pe, &
547 3*n_point_offset_from_pe, &
548 mpi_double_precision, i, neko_comm, ierr)
549
550 end do
551
552 call this%free_points()
553
554 this%n_points = n_new_points
556 call copy(this%xyz, new_xyz, 3 * n_new_points)
557
558 deallocate(n_point_offset_from_pe)
559 deallocate(n_points_from_pe)
560 deallocate(n_points_per_pe)
561 deallocate(xyz_send_to_pe)
562
563 end subroutine global_interpolation_redist
564
565
566
570 subroutine global_interpolation_evaluate(this, interp_values, field)
571 class(global_interpolation_t), intent(inout) :: this
572 real(kind=rp), intent(inout) :: interp_values(this%n_points)
573 real(kind=rp), intent(inout) :: field(this%nelv*this%Xh%lxyz)
574
575
576#ifdef HAVE_GSLIB
577 if (.not. this%all_points_local) then
578 call fgslib_findpts_eval(this%gs_handle, interp_values, &
579 1, this%error_code, 1, &
580 this%proc_owner, 1, this%el_owner, 1, &
581 this%rst, this%gdim, &
582 this%n_points, field)
583 else
584 if (this%n_points .gt. 0) &
585 call this%local_interp%evaluate(interp_values, this%el_owner, &
586 field, this%nelv)
587 end if
588#else
589 call neko_error('Neko needs to be built with GSLIB support')
590#endif
591
592 end subroutine global_interpolation_evaluate
593
594end module global_interpolation
Map a Fortran array to a device (allocate and associate)
Definition device.F90:68
Copy data between host and device (or device and device)
Definition device.F90:62
Definition comm.F90:1
integer pe_rank
MPI rank.
Definition comm.F90:28
type(mpi_comm) neko_comm
MPI communicator.
Definition comm.F90:16
integer pe_size
MPI size of communicator.
Definition comm.F90:31
Device abstraction, common interface for various accelerators.
Definition device.F90:34
integer, parameter, public host_to_device
Definition device.F90:46
subroutine, public device_free(x_d)
Deallocate memory on the device.
Definition device.F90:197
Defines a mapping of the degrees of freedom.
Definition dofmap.f90:35
Defines a field.
Definition field.f90:34
Implements global_interpolation given a dofmap.
subroutine global_interpolation_free(this)
Destructor.
subroutine global_interpolation_free_points(this)
Destructor for point arrays.
subroutine global_interpolation_init_point_arrays(this)
subroutine global_interpolation_init_xyz(this, x, y, z, gdim, nelv, xh, tol)
Initialize the global interpolation object on a set of coordinates.
subroutine global_interpolation_find_coords(this, x, y, z, n_points)
Finds the corresponding r,s,t coordinates in the correct global element as well as which process that...
subroutine global_interpolation_evaluate(this, interp_values, field)
Evalute the interpolated value in the points given a field on the dofmap.
subroutine global_interpolation_init_dof(this, dof, tol)
Initialize the global interpolation object on a dofmap.
subroutine global_interpolation_find_xyz(this, xyz, n_points)
Finds the corresponding r,s,t coordinates in the correct global element as well as which process that...
subroutine global_interpolation_find_and_redist(this, xyz, n_points)
Finds the corresponding r,s,t coordinates and redistributes the points to the owning rank in the corr...
subroutine global_interpolation_find_common(this)
Common routine for finding the points.
subroutine global_interpolation_redist(this)
Routines to obtain interpolated values on a set of points with known rst coordinates in elements loca...
Logging routines.
Definition log.f90:34
type(log_t), public neko_log
Global log stream.
Definition log.f90:65
integer, parameter, public log_size
Definition log.f90:42
Definition math.f90:60
subroutine, public copy(a, b, n)
Copy a vector .
Definition math.f90:238
Defines a mesh.
Definition mesh.f90:34
Build configurations.
integer, parameter neko_bcknd_device
MPI derived types.
Definition mpi_types.f90:34
integer, parameter, public rp
Global precision used in computations.
Definition num_types.f90:12
Defines a function space.
Definition space.f90:34
Defines structs that are used... Dont know if we should keep it though.
Definition structs.f90:2
Utilities.
Definition utils.f90:35
subroutine, public neko_warning(warning_msg)
Reports a warning to standard output.
Definition utils.f90:266
Implements global interpolation for arbitrary points in the domain.
Interpolation on a set of points with known rst coordinates in elements local to this process....
The function space for the SEM solution fields.
Definition space.f90:62
Pointer to array.
Definition structs.f90:14