49 mpi_double_precision, mpi_integer, mpi_sum, mpi_reduce, mpi_gatherv
53 use,
intrinsic :: iso_c_binding
74 logical :: all_points_local = .false.
78 logical :: gs_init = .false.
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
108 procedure, pass(this) :: find_points_and_redist => &
113 procedure, pass(this) :: find_points_coords => &
116 generic :: find_points => find_points_xyz, find_points_coords
121 generic :: init => init_dof, init_xyz
133 real(kind=
rp),
optional :: tol
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)
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
171 if (
present(tol)) this%tol = tol
174 max_pts_per_iter = 128
179 call fgslib_findpts_setup(this%gs_handle, &
182 this%x%ptr, this%y%ptr, this%z%ptr, &
186 lx*ly*lz*nelv, lx*ly*lz*nelv, &
187 max_pts_per_iter, this%tol)
188 this%gs_init = .true.
190 call neko_error(
'Neko needs to be built with GSLIB support')
208 call this%free_points()
211 if (this%gs_init)
then
212 call fgslib_findpts_free(this%gs_handle)
213 this%gs_init = .false.
224 this%all_points_local = .false.
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)
233 if (c_associated(this%el_owner_d))
then
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(:)
255 call fgslib_findpts(this%gs_handle, &
256 this%error_code, 1, &
257 this%proc_owner, 1, &
259 this%rst, this%gdim, &
261 this%xyz(1,1), this%gdim, &
262 this%xyz(2,1), this%gdim, &
263 this%xyz(3,1), this%gdim, this%n_points)
265 do i = 1 , this%n_points
270 if (this%error_code(i) .eq. 1)
then
271 if (this%dist2(i) .gt. this%tol)
then
272 write(*,*)
'Point with coords: ',&
276 'Did not converge to tol. Absolute differences squared: ',&
277 this%dist2(i),
'PE rank',
pe_rank
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: ', &
287 'el_owner, rst coords, pe: ',&
288 this%el_owner(i), this%rst(1,i), this%rst(2,i), &
293 allocate(x_check(this%n_points))
294 allocate(y_check(this%n_points))
295 allocate(z_check(this%n_points))
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)
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)
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)
315 do i = 1 , this%n_points
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.
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: ', &
348 call neko_error(
'Neko needs to be built with GSLIB support')
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)
375 call this%free_points()
377 this%n_points = 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)
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))
402 call device_map(this%el_owner, this%el_owner_d, this%n_points)
418 integer,
intent(in) :: n_points
420 real(kind=
rp),
intent(inout) :: xyz(3, n_points)
423 call this%free_points()
425 this%n_points = n_points
430 call copy(this%xyz, xyz, 3 * n_points)
449 integer,
intent(inout) :: n_points
451 real(kind=
rp),
allocatable,
intent(inout) :: xyz(:,:)
455 call this%free_points()
457 this%n_points = n_points
461 call copy(this%xyz, xyz, 3 * n_points)
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)
477 n_points = this%n_points
478 if (
allocated(xyz))
then
481 allocate(xyz(3, n_points))
482 call copy(xyz, this%xyz, 3*n_points)
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.
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
502 allocate(n_points_per_pe(0:(
pe_size-1)))
503 allocate(n_points_from_pe(0:(
pe_size-1)))
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
514 call mpi_reduce(n_points_per_pe(i), n_new_points, 1, mpi_integer, &
516 call mpi_gather(n_points_per_pe(i), 1, mpi_integer,&
517 n_points_from_pe, 1, mpi_integer, i,
neko_comm, ierr)
520 allocate(n_point_offset_from_pe(0:(
pe_size-1)))
521 n_point_offset_from_pe(0) = 0
523 n_point_offset_from_pe(i) = n_points_from_pe(i-1)&
524 + n_point_offset_from_pe(i-1)
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))
533 do j = 1, this%n_points
534 if (this%proc_owner(j) .eq. i)
then
536 xyz_send_to_pe(:,k) = this%xyz(:,j)
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: ', &
543 call neko_error(
'Error in redistribution of points')
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)
552 call this%free_points()
554 this%n_points = n_new_points
556 call copy(this%xyz, new_xyz, 3 * n_new_points)
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)
572 real(kind=
rp),
intent(inout) :: interp_values(this%n_points)
573 real(kind=
rp),
intent(inout) ::
field(this%nelv*this%Xh%lxyz)
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)
584 if (this%n_points .gt. 0) &
585 call this%local_interp%evaluate(interp_values, this%el_owner, &
589 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)
type(mpi_comm) neko_comm
MPI communicator.
integer pe_size
MPI size of 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.