Neko 1.99.1
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
49 use mpi_f08, only : mpi_gather, mpi_double_precision, mpi_integer, mpi_sum, &
50 mpi_reduce, mpi_gatherv
51 use math, only : copy
53 use structs, only : array_ptr_t
54 use, intrinsic :: iso_c_binding
55 implicit none
56 private
57
59 type, public :: global_interpolation_t
61 type(array_ptr_t) :: x
63 type(array_ptr_t) :: y
65 type(array_ptr_t) :: z
67 integer :: gdim
69 integer :: nelv
71 type(space_t), pointer :: xh
73 type(local_interpolator_t) :: local_interp
75 logical :: all_points_local = .false.
76 !! Gslib handle.
77 !! @note: Remove when we remove gslib.
78 integer :: gs_handle
79 logical :: gs_init = .false.
82 integer :: n_points
85 real(kind=rp), allocatable :: xyz(:,:)
87 integer, allocatable :: proc_owner(:)
89 integer, allocatable :: el_owner(:)
90 type(c_ptr) :: el_owner_d = c_null_ptr
93 real(kind=rp), allocatable :: rst(:,:)
96 real(kind=rp), allocatable :: dist2(:)
98 integer, allocatable :: error_code(:)
100 real(kind=rp) :: tol = 5d-13
101 contains
103 procedure, pass(this) :: init_xyz => global_interpolation_init_xyz
104 procedure, pass(this) :: init_dof => global_interpolation_init_dof
106 procedure, pass(this) :: free => global_interpolation_free
108 procedure, pass(this) :: free_points => global_interpolation_free_points
109 procedure, pass(this) :: find_points_and_redist => &
114 procedure, pass(this) :: find_points_coords => &
116 procedure, pass(this) :: find_points_xyz => global_interpolation_find_xyz
117 generic :: find_points => find_points_xyz, find_points_coords
119 procedure, pass(this) :: evaluate => global_interpolation_evaluate
120
122 generic :: init => init_dof, init_xyz
123
125
126contains
127
131 subroutine global_interpolation_init_dof(this, dof, tol)
132 class(global_interpolation_t), intent(inout) :: this
133 type(dofmap_t), target :: dof
134 real(kind=rp), optional :: tol
135
136 ! NOTE: Passing dof%x(:,1,1,1), etc in init_xyz passes down the entire
137 ! dof%x array and not a slice. It is done this way for
138 ! this%x%ptr to point to dof%x (see global_interpolation_init_xyz).
139 call this%init_xyz(dof%x(:,1,1,1), dof%y(:,1,1,1), dof%z(:,1,1,1), &
140 dof%msh%gdim, dof%msh%nelv, dof%Xh, tol = tol)
141
142 end subroutine global_interpolation_init_dof
143
153 subroutine global_interpolation_init_xyz(this, x, y, z, gdim, nelv, Xh, tol)
154 class(global_interpolation_t), intent(inout) :: this
155 real(kind=rp), intent(in), target :: x(:)
156 real(kind=rp), intent(in), target :: y(:)
157 real(kind=rp), intent(in), target :: z(:)
158 integer, intent(in) :: gdim
159 integer, intent(in) :: nelv
160 type(space_t), intent(in), target :: Xh
161 real(kind=rp), intent(in), optional :: tol
162 integer :: lx, ly, lz, max_pts_per_iter
163
164#ifdef HAVE_GSLIB
165
166 this%x%ptr => x
167 this%y%ptr => y
168 this%z%ptr => z
169 this%gdim = gdim
170 this%nelv = nelv
171 this%Xh => xh
172 if (present(tol)) this%tol = tol
173
174 ! Number of points to iterate on simultaneosuly
175 max_pts_per_iter = 128
176 lx = xh%lx
177 ly = xh%ly
178 lz = xh%lz
179
180 call fgslib_findpts_setup(this%gs_handle, &
182 this%gdim, &
183 this%x%ptr, this%y%ptr, this%z%ptr, & ! Physical nodal values
184 lx, ly, lz, nelv, & ! Mesh dimensions
185 2*lx, 2*ly, 2*lz, & ! Mesh size for bounding box computation
186 0.01, & ! relative size to expand bounding boxes by
187 lx*ly*lz*nelv, lx*ly*lz*nelv, & ! local/global hash mesh sizes
188 max_pts_per_iter, this%tol)
189 this%gs_init = .true.
190#else
191 call neko_error('Neko needs to be built with GSLIB support')
192#endif
193
194 end subroutine global_interpolation_init_xyz
195
196
199 class(global_interpolation_t), intent(inout) :: this
200
201 nullify(this%x%ptr)
202 nullify(this%y%ptr)
203 nullify(this%z%ptr)
204 nullify(this%Xh)
205
206 this%nelv = 0
207 this%gdim = 0
208
209 call this%free_points()
210
211#ifdef HAVE_GSLIB
212 if (this%gs_init) then
213 call fgslib_findpts_free(this%gs_handle)
214 this%gs_init = .false.
215 end if
216#endif
217
218 end subroutine global_interpolation_free
219
222 class(global_interpolation_t), intent(inout) :: this
223
224 this%n_points = 0
225 this%all_points_local = .false.
226
227 if (allocated(this%xyz)) deallocate(this%xyz)
228 if (allocated(this%rst)) deallocate(this%rst)
229 if (allocated(this%proc_owner)) deallocate(this%proc_owner)
230 if (allocated(this%el_owner)) deallocate(this%el_owner)
231 if (allocated(this%dist2)) deallocate(this%dist2)
232 if (allocated(this%error_code)) deallocate(this%error_code)
233
234 if (c_associated(this%el_owner_d)) then
235 call device_free(this%el_owner_d)
236 end if
237
239
242 class(global_interpolation_t), intent(inout) :: this
243 !!Perhaps this should be kind dp
244 real(kind=rp) :: xdiff, ydiff, zdiff
245 character(len=8000) :: log_buf
246 real(kind=rp), allocatable :: x_check(:)
247 real(kind=rp), allocatable :: y_check(:)
248 real(kind=rp), allocatable :: z_check(:)
249 logical :: isdiff
250 integer :: i
251
252
253#ifdef HAVE_GSLIB
254
255 ! gslib find points, which element they belong, to process etc.
256 call fgslib_findpts(this%gs_handle, &
257 this%error_code, 1, &
258 this%proc_owner, 1, &
259 this%el_owner, 1, &
260 this%rst, this%gdim, &
261 this%dist2, 1, &
262 this%xyz(1,1), this%gdim, &
263 this%xyz(2,1), this%gdim, &
264 this%xyz(3,1), this%gdim, this%n_points)
265
266 do i = 1 , this%n_points
267
268 !
269 ! Check validity of points
270 !
271 if (this%error_code(i) .eq. 1) then
272 if (this%dist2(i) .gt. this%tol) then
273 write(*,*) 'Point with coords: ',&
274 this%xyz(1,i),&
275 this%xyz(2,i),&
276 this%xyz(3,i),&
277 'Did not converge to tol. Absolute differences squared: ',&
278 this%dist2(i), 'PE rank', pe_rank
279 end if
280 end if
281
282 if (this%error_code(i) .eq. 2) &
283 write(*,*) 'Point with coords: ',&
284 this%xyz(1,i), this%xyz(2,i), this%xyz(3,i),&
285 'Outside the mesh!',&
286 ' Interpolation on these points will return 0.0. dist2: ', &
287 this%dist2(i),&
288 'el_owner, rst coords, pe: ',&
289 this%el_owner(i), this%rst(1,i), this%rst(2,i), &
290 this%rst(3,i), pe_rank
291
292 end do
293
294 allocate(x_check(this%n_points))
295 allocate(y_check(this%n_points))
296 allocate(z_check(this%n_points))
297
298 call fgslib_findpts_eval(this%gs_handle, x_check, &
299 1, this%error_code, 1, &
300 this%proc_owner, 1, this%el_owner, 1, &
301 this%rst, this%gdim, &
302 this%n_points, this%x%ptr)
303
304 call fgslib_findpts_eval(this%gs_handle, y_check, &
305 1, this%error_code, 1, &
306 this%proc_owner, 1, this%el_owner, 1, &
307 this%rst, this%gdim, &
308 this%n_points, this%y%ptr)
309
310 call fgslib_findpts_eval(this%gs_handle, z_check, &
311 1, this%error_code, 1, &
312 this%proc_owner, 1, this%el_owner, 1, &
313 this%rst, this%gdim, &
314 this%n_points, this%z%ptr)
315
316 do i = 1 , this%n_points
317
318 !
319 ! Check validity of points
320 !
321 isdiff = .false.
322 xdiff = (x_check(i)-this%xyz(1,i))**2
323 if ( xdiff .gt. this%tol) isdiff = .true.
324 ydiff = (y_check(i)-this%xyz(2,i))**2
325 if ( ydiff .gt. this%tol) isdiff = .true.
326 zdiff = (z_check(i)-this%xyz(3,i))**2
327 if ( zdiff .gt. this%tol) isdiff = .true.
328
329 if (isdiff) then
330 write(*,*) 'Points with coords: ', &
331 this%xyz(1, i), this%xyz(2, i), this%xyz(3, i), &
332 'Differ from interpolated coords: ', &
333 x_check(i), y_check(i), z_check(i), &
334 'Distance squared: ', &
335 xdiff, ydiff, zdiff
336 end if
337
338 end do
339
340 deallocate(x_check)
341 deallocate(y_check)
342 deallocate(z_check)
343
344 if (neko_bcknd_device .eq. 1) then
345 call device_memcpy(this%el_owner, this%el_owner_d, &
346 this%n_points, host_to_device, sync = .true.)
347 end if
348#else
349 call neko_error('Neko needs to be built with GSLIB support')
350#endif
352
365 subroutine global_interpolation_find_coords(this, x, y, z, n_points)
366 class(global_interpolation_t), intent(inout) :: this
367 integer :: n_points
368 !!Perhaps this should be kind dp
369 !!this is to get around that x,y,z often is 4 dimensional...
370 !!Should maybe add interface for 1d aswell
371 real(kind=rp) :: x(n_points,1,1,1)
372 real(kind=rp) :: y(n_points,1,1,1)
373 real(kind=rp) :: z(n_points,1,1,1)
374 integer :: i
375
376 call this%free_points()
377
378 this%n_points = n_points
379
381
382 do i = 1, n_points
383 this%xyz(1, i) = x(i,1,1,1)
384 this%xyz(2, i) = y(i,1,1,1)
385 this%xyz(3, i) = z(i,1,1,1)
386 end do
387
389
391
393 class(global_interpolation_t) :: this
394
395 allocate(this%xyz(3, this%n_points))
396 allocate(this%rst(3, this%n_points))
397 allocate(this%proc_owner(this%n_points))
398 allocate(this%el_owner(this%n_points))
399 allocate(this%dist2(this%n_points))
400 allocate(this%error_code(this%n_points))
401
402 if (neko_bcknd_device .eq. 1) &
403 call device_map(this%el_owner, this%el_owner_d, this%n_points)
404
406
417 subroutine global_interpolation_find_xyz(this, xyz, n_points)
418 class(global_interpolation_t), intent(inout) :: this
419 integer, intent(in) :: n_points
420 !!Perhaps this should be kind dp
421 real(kind=rp), intent(inout) :: xyz(3, n_points)
422
423
424 call this%free_points()
425
426 this%n_points = n_points
427
429
431 call copy(this%xyz, xyz, 3 * n_points)
432
434
435 end subroutine global_interpolation_find_xyz
436
448 subroutine global_interpolation_find_and_redist(this, xyz, n_points)
449 class(global_interpolation_t), intent(inout) :: this
450 integer, intent(inout) :: n_points
451 !!Perhaps this should be kind dp
452 real(kind=rp), allocatable, intent(inout) :: xyz(:,:)
453 integer :: i
454
455
456 call this%free_points()
457
458 this%n_points = n_points
460
462 call copy(this%xyz, xyz, 3 * n_points)
463
468
469 do i = 1, this%n_points
470 if (this%proc_owner(i) .ne. pe_rank) then
471 write(*,*) 'Redistribution failed on rank: ', pe_rank, &
472 'for point with coord: ', &
473 this%xyz(1, i), this%xyz(2, i), this%xyz(3, i)
474 exit
475 end if
476 end do
477
478 n_points = this%n_points
479 if (allocated(xyz)) then
480 deallocate(xyz)
481 end if
482 allocate(xyz(3, n_points))
483 call copy(xyz, this%xyz, 3*n_points)
484
485 call this%local_interp%init(this%Xh, this%rst(1,:),&
486 this%rst(2,:), this%rst(3,:), n_points)
487 this%all_points_local = .true.
488
489
491
493 class(global_interpolation_t), intent(inout) :: this
494 integer, allocatable :: n_points_per_pe(:)
495 integer, allocatable :: n_points_from_pe(:)
496 integer, allocatable :: n_point_offset_from_pe(:)
497 real(kind=rp), allocatable :: xyz_send_to_pe(:,:)
498 real(kind=rp), allocatable :: new_xyz(:,:)
499 integer :: i, j, k, ierr, n_new_points, max_n_points_to_send
500
501 n_new_points = 0
502
503 allocate(n_points_per_pe(0:(pe_size-1)))
504 allocate(n_points_from_pe(0:(pe_size-1)))
505 n_points_per_pe = 0
506 n_points_from_pe = 0
508 do i = 1, this%n_points
509 n_points_per_pe(this%proc_owner(i)) = &
510 n_points_per_pe(this%proc_owner(i)) + 1
511 end do
514 do i = 0, (pe_size - 1)
515 call mpi_reduce(n_points_per_pe(i), n_new_points, 1, mpi_integer, &
516 mpi_sum, i, neko_comm, ierr)
517 call mpi_gather(n_points_per_pe(i), 1, mpi_integer,&
518 n_points_from_pe, 1, mpi_integer, i, neko_comm, ierr)
519 end do
520
521 allocate(n_point_offset_from_pe(0:(pe_size-1)))
522 n_point_offset_from_pe(0) = 0
523 do i = 1, (pe_size - 1)
524 n_point_offset_from_pe(i) = n_points_from_pe(i-1)&
525 + n_point_offset_from_pe(i-1)
526 end do
527
528 allocate(new_xyz(3, n_new_points))
529 max_n_points_to_send = maxval(n_points_per_pe)
530 allocate(xyz_send_to_pe(3, max_n_points_to_send))
531 do i = 0, (pe_size - 1)
533 k = 0
534 do j = 1, this%n_points
535 if (this%proc_owner(j) .eq. i) then
536 k = k + 1
537 xyz_send_to_pe(:,k) = this%xyz(:,j)
538 end if
539 end do
540 if (k .ne. n_points_per_pe(i)) then
541 write(*,*) 'PE: ', pe_rank, ' has k= ', k, &
542 'points for PE:', i,' but should have: ', &
543 n_points_per_pe(i)
544 call neko_error('Error in redistribution of points')
545 end if
546 call mpi_gatherv(xyz_send_to_pe,3*n_points_per_pe(i), &
547 mpi_double_precision, new_xyz,3*n_points_from_pe, &
548 3*n_point_offset_from_pe, &
549 mpi_double_precision, i, neko_comm, ierr)
550
551 end do
552
553 call this%free_points()
554
555 this%n_points = n_new_points
557 call copy(this%xyz, new_xyz, 3 * n_new_points)
558
559 deallocate(n_point_offset_from_pe)
560 deallocate(n_points_from_pe)
561 deallocate(n_points_per_pe)
562 deallocate(xyz_send_to_pe)
563
564 end subroutine global_interpolation_redist
565
566
567
571 subroutine global_interpolation_evaluate(this, interp_values, field)
572 class(global_interpolation_t), intent(inout) :: this
573 real(kind=rp), intent(inout) :: interp_values(this%n_points)
574 real(kind=rp), intent(inout) :: field(this%nelv*this%Xh%lxyz)
575
576
577#ifdef HAVE_GSLIB
578 if (.not. this%all_points_local) then
579 call fgslib_findpts_eval(this%gs_handle, interp_values, &
580 1, this%error_code, 1, &
581 this%proc_owner, 1, this%el_owner, 1, &
582 this%rst, this%gdim, &
583 this%n_points, field)
584 else
585 if (this%n_points .gt. 0) &
586 call this%local_interp%evaluate(interp_values, this%el_owner, &
587 field, this%nelv)
588 end if
589#else
590 call neko_error('Neko needs to be built with GSLIB support')
591#endif
592
593 end subroutine global_interpolation_evaluate
594
595end module global_interpolation
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
Definition comm.F90:1
integer, public pe_size
MPI size of communicator.
Definition comm.F90:58
integer, public pe_rank
MPI rank.
Definition comm.F90:55
type(mpi_comm), public neko_comm
MPI communicator.
Definition comm.F90:42
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
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:70
integer, parameter, public log_size
Definition log.f90:46
Definition math.f90:60
subroutine, public copy(a, b, n)
Copy a vector .
Definition math.f90:255
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:284
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