49 use mpi_f08,
only : mpi_gather, mpi_double_precision, mpi_integer, mpi_sum, &
50 mpi_reduce, mpi_gatherv
54 use,
intrinsic :: iso_c_binding
75 logical :: all_points_local = .false.
79 logical :: gs_init = .false.
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
109 procedure, pass(this) :: find_points_and_redist => &
114 procedure, pass(this) :: find_points_coords => &
117 generic :: find_points => find_points_xyz, find_points_coords
122 generic :: init => init_dof, init_xyz
134 real(kind=
rp),
optional :: tol
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)
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
172 if (
present(tol)) this%tol = tol
175 max_pts_per_iter = 128
180 call fgslib_findpts_setup(this%gs_handle, &
183 this%x%ptr, this%y%ptr, this%z%ptr, &
187 lx*ly*lz*nelv, lx*ly*lz*nelv, &
188 max_pts_per_iter, this%tol)
189 this%gs_init = .true.
191 call neko_error(
'Neko needs to be built with GSLIB support')
209 call this%free_points()
212 if (this%gs_init)
then
213 call fgslib_findpts_free(this%gs_handle)
214 this%gs_init = .false.
225 this%all_points_local = .false.
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)
234 if (c_associated(this%el_owner_d))
then
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(:)
256 call fgslib_findpts(this%gs_handle, &
257 this%error_code, 1, &
258 this%proc_owner, 1, &
260 this%rst, this%gdim, &
262 this%xyz(1,1), this%gdim, &
263 this%xyz(2,1), this%gdim, &
264 this%xyz(3,1), this%gdim, this%n_points)
266 do i = 1 , this%n_points
271 if (this%error_code(i) .eq. 1)
then
272 if (this%dist2(i) .gt. this%tol)
then
273 write(*,*)
'Point with coords: ',&
277 'Did not converge to tol. Absolute differences squared: ',&
278 this%dist2(i),
'PE rank',
pe_rank
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: ', &
288 'el_owner, rst coords, pe: ',&
289 this%el_owner(i), this%rst(1,i), this%rst(2,i), &
294 allocate(x_check(this%n_points))
295 allocate(y_check(this%n_points))
296 allocate(z_check(this%n_points))
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)
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)
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)
316 do i = 1 , this%n_points
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.
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: ', &
349 call neko_error(
'Neko needs to be built with GSLIB support')
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)
376 call this%free_points()
378 this%n_points = 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)
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))
403 call device_map(this%el_owner, this%el_owner_d, this%n_points)
419 integer,
intent(in) :: n_points
421 real(kind=
rp),
intent(inout) :: xyz(3, n_points)
424 call this%free_points()
426 this%n_points = n_points
431 call copy(this%xyz, xyz, 3 * n_points)
450 integer,
intent(inout) :: n_points
452 real(kind=
rp),
allocatable,
intent(inout) :: xyz(:,:)
456 call this%free_points()
458 this%n_points = n_points
462 call copy(this%xyz, xyz, 3 * n_points)
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)
478 n_points = this%n_points
479 if (
allocated(xyz))
then
482 allocate(xyz(3, n_points))
483 call copy(xyz, this%xyz, 3*n_points)
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.
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
503 allocate(n_points_per_pe(0:(
pe_size-1)))
504 allocate(n_points_from_pe(0:(
pe_size-1)))
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
515 call mpi_reduce(n_points_per_pe(i), n_new_points, 1, mpi_integer, &
517 call mpi_gather(n_points_per_pe(i), 1, mpi_integer,&
518 n_points_from_pe, 1, mpi_integer, i,
neko_comm, ierr)
521 allocate(n_point_offset_from_pe(0:(
pe_size-1)))
522 n_point_offset_from_pe(0) = 0
524 n_point_offset_from_pe(i) = n_points_from_pe(i-1)&
525 + n_point_offset_from_pe(i-1)
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))
534 do j = 1, this%n_points
535 if (this%proc_owner(j) .eq. i)
then
537 xyz_send_to_pe(:,k) = this%xyz(:,j)
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: ', &
544 call neko_error(
'Error in redistribution of points')
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)
553 call this%free_points()
555 this%n_points = n_new_points
557 call copy(this%xyz, new_xyz, 3 * n_new_points)
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)
573 real(kind=
rp),
intent(inout) :: interp_values(this%n_points)
574 real(kind=
rp),
intent(inout) ::
field(this%nelv*this%Xh%lxyz)
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)
585 if (this%n_points .gt. 0) &
586 call this%local_interp%evaluate(interp_values, this%el_owner, &
590 call neko_error(
'Neko needs to be built with GSLIB support')
Map a Fortran array to a device (allocate and associate)
Copy data between host and device (or device and device)
integer, public pe_size
MPI size of communicator.
integer, public pe_rank
MPI rank.
type(mpi_comm), public neko_comm
MPI communicator.
Device abstraction, common interface for various accelerators.
integer, parameter, public host_to_device
subroutine, public device_free(x_d)
Deallocate memory on the device.
Defines a mapping of the degrees of freedom.
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...
type(log_t), public neko_log
Global log stream.
integer, parameter, public log_size
subroutine, public copy(a, b, n)
Copy a vector .
integer, parameter neko_bcknd_device
integer, parameter, public rp
Global precision used in computations.
Defines a function space.
Defines structs that are used... Dont know if we should keep it though.
subroutine, public neko_warning(warning_msg)
Reports a warning to standard output.
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.