49 use,
intrinsic :: iso_c_binding
63 logical :: all_points_local = .false.
67 logical :: gs_init = .false.
73 real(kind=
rp),
allocatable :: xyz(:,:)
75 integer,
allocatable :: proc_owner(:)
77 integer,
allocatable :: el_owner(:)
78 type(c_ptr) :: el_owner_d = c_null_ptr
81 real(kind=
rp),
allocatable :: rst(:,:)
84 real(kind=
rp),
allocatable :: dist2(:)
86 integer,
allocatable :: error_code(:)
88 real(kind=
rp) :: tol = 5e-13
96 procedure, pass(this) :: find_points_and_redist => &
103 generic :: find_points => find_points_xyz, find_points_coords
117 real(kind=
rp),
optional :: tol
118 integer :: lx, ly, lz, nelv, max_pts_per_iter
123 if(
present(tol)) this%tol = tol
130 nelv = this%mesh%nelv
132 max_pts_per_iter = 128
134 call fgslib_findpts_setup(this%gs_handle, &
137 dof%x, dof%y, dof%z, &
141 lx*ly*lz*nelv, lx*ly*lz*nelv, &
142 max_pts_per_iter, this%tol)
143 this%gs_init = .true.
145 call neko_error(
'Neko needs to be built with GSLIB support')
159 call this%free_points()
162 if (this%gs_init)
then
163 call fgslib_findpts_free(this%gs_handle)
164 this%gs_init = .false.
175 this%all_points_local = .false.
177 if (
allocated(this%xyz))
deallocate(this%xyz)
178 if (
allocated(this%rst))
deallocate(this%rst)
179 if (
allocated(this%proc_owner))
deallocate(this%proc_owner)
180 if (
allocated(this%el_owner))
deallocate(this%el_owner)
181 if (
allocated(this%dist2))
deallocate(this%dist2)
182 if (
allocated(this%error_code))
deallocate(this%error_code)
184 if (c_associated(this%el_owner_d))
then
185 call device_free(this%el_owner_d)
194 real(kind=
rp) :: xdiff, ydiff, zdiff
195 character(len=8000) :: log_buf
196 real(kind=
rp),
allocatable :: x_check(:)
197 real(kind=
rp),
allocatable :: y_check(:)
198 real(kind=
rp),
allocatable :: z_check(:)
206 call fgslib_findpts(this%gs_handle, &
207 this%error_code, 1, &
208 this%proc_owner, 1, &
210 this%rst, this%mesh%gdim, &
212 this%xyz(1,1), this%mesh%gdim, &
213 this%xyz(2,1), this%mesh%gdim, &
214 this%xyz(3,1), this%mesh%gdim, this%n_points)
221 if (this%error_code(i) .eq. 1)
then
222 if (this%dist2(i) .gt. this%tol)
then
223 write(*,*)
'Point with coords: ',&
227 'Did not converge to tol. Absolute differences squared: ',&
228 this%dist2(i),
'PE rank',
pe_rank
232 if (this%error_code(i) .eq. 2) &
233 write(*,*)
'Point with coords: ',&
234 this%xyz(1,i), this%xyz(2,i), this%xyz(3,i),&
235 'Outside the mesh!',&
236 ' Interpolation on these points will return 0.0. dist2: ', &
238 'el_owner, rst coords, pe: ',&
239 this%el_owner(i), this%rst(1,i), this%rst(2,i), this%rst(3,i),
pe_rank
243 allocate(x_check(this%n_points))
244 allocate(y_check(this%n_points))
245 allocate(z_check(this%n_points))
247 call fgslib_findpts_eval(this%gs_handle, x_check, &
248 1, this%error_code, 1, &
249 this%proc_owner, 1, this%el_owner, 1, &
250 this%rst, this%mesh%gdim, &
251 this%n_points, this%dof%x)
253 call fgslib_findpts_eval(this%gs_handle, y_check, &
254 1, this%error_code, 1, &
255 this%proc_owner, 1, this%el_owner, 1, &
256 this%rst, this%mesh%gdim, &
257 this%n_points, this%dof%y)
259 call fgslib_findpts_eval(this%gs_handle, z_check, &
260 1, this%error_code, 1, &
261 this%proc_owner, 1, this%el_owner, 1, &
262 this%rst, this%mesh%gdim, &
263 this%n_points, this%dof%z)
272 xdiff = (x_check(i)-this%xyz(1,i))**2
273 if ( xdiff .gt. this%tol) isdiff = .true.
274 ydiff = (y_check(i)-this%xyz(2,i))**2
275 if ( ydiff .gt. this%tol) isdiff = .true.
276 zdiff = (z_check(i)-this%xyz(3,i))**2
277 if ( zdiff .gt. this%tol) isdiff = .true.
280 write(*,*)
'Points with coords: ',&
281 this%xyz(1,i),this%xyz(2,i),this%xyz(3,i), &
282 'Differ from interpolated coords: ',&
283 x_check(i), y_check(i), z_check(i),&
284 'Distance squared: ',&
294 if (neko_bcknd_device .eq. 1)
then
295 call device_memcpy(this%el_owner, this%el_owner_d, &
296 this%n_points, host_to_device, sync = .true.)
299 call neko_error(
'Neko needs to be built with GSLIB support')
321 real(kind=
rp) :: x(n_points,1,1,1)
322 real(kind=
rp) :: y(n_points,1,1,1)
323 real(kind=
rp) :: z(n_points,1,1,1)
326 call this%free_points()
328 this%n_points = n_points
333 this%xyz(1,i) = x(i,1,1,1)
334 this%xyz(2,i) = y(i,1,1,1)
335 this%xyz(3,i) = z(i,1,1,1)
345 allocate(this%xyz(3,this%n_points))
346 allocate(this%rst(3,this%n_points))
347 allocate(this%proc_owner(this%n_points))
348 allocate(this%el_owner(this%n_points))
349 allocate(this%dist2(this%n_points))
350 allocate(this%error_code(this%n_points))
352 if (neko_bcknd_device .eq. 1) &
353 call device_map(this%el_owner, this%el_owner_d,this%n_points)
369 integer,
intent(in) :: n_points
371 real(kind=
rp),
intent(inout) :: xyz(3,n_points)
374 call this%free_points()
376 this%n_points = n_points
381 call copy(this%xyz,xyz,3*n_points)
399 integer,
intent(inout) :: n_points
401 real(kind=
rp),
allocatable,
intent(inout) :: xyz(:,:)
405 call this%free_points()
407 this%n_points = n_points
411 call copy(this%xyz,xyz,3*n_points)
418 do i = 1, this%n_points
419 if (this%proc_owner(i) .ne.
pe_rank)
then
420 write(*,*)
'Redistribution failed on rank: ',
pe_rank,&
421 'for point with coord: ', &
422 this%xyz(1,i),this%xyz(2,i),this%xyz(3,i)
427 n_points = this%n_points
429 allocate(xyz(3,n_points))
430 call copy(xyz,this%xyz,3*n_points)
432 call this%local_interp%init(this%Xh, this%rst(1,:),&
433 this%rst(2,:), this%rst(3,:), n_points)
434 this%all_points_local = .true.
441 integer,
allocatable :: n_points_per_pe(:)
442 integer,
allocatable :: n_points_from_pe(:)
443 integer,
allocatable :: n_point_offset_from_pe(:)
444 real(kind=
rp),
allocatable :: xyz_send_to_pe(:,:)
445 real(kind=
rp),
allocatable :: new_xyz(:,:)
446 integer :: i, j, k, ierr, n_new_points, max_n_points_to_send
450 allocate(n_points_per_pe(0:(
pe_size-1)))
451 allocate(n_points_from_pe(0:(
pe_size-1)))
455 do i = 1,this%n_points
456 n_points_per_pe(this%proc_owner(i)) = n_points_per_pe(this%proc_owner(i)) + 1
461 call mpi_reduce(n_points_per_pe(i),n_new_points,1,mpi_integer,mpi_sum, i,
neko_comm, ierr)
462 call mpi_gather(n_points_per_pe(i), 1, mpi_integer,&
463 n_points_from_pe, 1, mpi_integer, i,
neko_comm, ierr)
466 allocate(n_point_offset_from_pe(0:(
pe_size-1)))
467 n_point_offset_from_pe(0) = 0
469 n_point_offset_from_pe(i) = n_points_from_pe(i-1)&
470 + n_point_offset_from_pe(i-1)
473 allocate(new_xyz(3,n_new_points))
474 max_n_points_to_send = maxval(n_points_per_pe)
475 allocate(xyz_send_to_pe(3,max_n_points_to_send))
479 do j = 1, this%n_points
480 if (this%proc_owner(j) .eq. i)
then
482 xyz_send_to_pe(:,k) = this%xyz(:,j)
485 if (k .ne. n_points_per_pe(i))
then
486 write(*,*)
'PE: ',
pe_rank,
' has k= ', k,&
487 'points for PE:', i,
' but should have: ',&
489 call neko_error(
'Error in redistribution of points')
491 call mpi_gatherv(xyz_send_to_pe,3*n_points_per_pe(i),&
492 mpi_double_precision, new_xyz,3*n_points_from_pe,&
493 3*n_point_offset_from_pe,&
494 mpi_double_precision, i,
neko_comm, ierr)
498 call this%free_points()
500 this%n_points = n_new_points
502 call copy(this%xyz,new_xyz,3*n_new_points)
504 deallocate(n_point_offset_from_pe)
505 deallocate(n_points_from_pe)
506 deallocate(n_points_per_pe)
507 deallocate(xyz_send_to_pe)
518 real(kind=
rp),
intent(inout) :: interp_values(this%n_points)
519 real(kind=
rp),
intent(inout) ::
field(this%dof%size())
522 if (.not. this%all_points_local)
then
523 call fgslib_findpts_eval(this%gs_handle, interp_values, &
524 1, this%error_code, 1, &
525 this%proc_owner, 1, this%el_owner, 1, &
526 this%rst, this%mesh%gdim, &
527 this%n_points,
field)
529 if (this%n_points .gt. 0) &
530 call this%local_interp%evaluate(interp_values, this%el_owner,&
531 field, this%mesh%nelv)
534 call neko_error(
'Neko needs to be built with GSLIB support')
type(mpi_comm) neko_comm
MPI communicator.
integer pe_size
MPI size of communicator.
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_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_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)
subroutine global_interpolation_init(this, dof, tol)
Evaluate only local points.
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, public rp
Global precision used in computations.
Defines a function space.
subroutine neko_warning(warning_msg)
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.