Neko  0.9.0
A portable framework for high-order spectral element flow simulations
local_interpolation.f90
Go to the documentation of this file.
1 ! Copyright (c) 2021-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 !
37  use space, only: space_t, gl, gll
38  use num_types, only: rp
39  use point, only: point_t
40  use math, only: abscmp
42  use utils, only: neko_error
43  use field, only: field_t
44  use field_list, only: field_list_t
45  use device
46  use device_math, only: device_rzero
48  implicit none
49 
57  type(space_t), pointer :: xh => null()
59  integer :: n_points
61  real(kind=rp), allocatable :: weights_r(:,:)
62  real(kind=rp), allocatable :: weights_s(:,:)
63  real(kind=rp), allocatable :: weights_t(:,:)
64  type(c_ptr) :: weights_r_d = c_null_ptr
65  type(c_ptr) :: weights_s_d = c_null_ptr
66  type(c_ptr) :: weights_t_d = c_null_ptr
67  contains
69  procedure, pass(this) :: init => local_interpolator_init
71  procedure, pass(this) :: free => local_interpolator_free
73  procedure, pass(this) :: evaluate => local_interpolator_evaluate
75  procedure, pass(this) :: compute_weights => local_interpolator_compute_weights
76 
77  end type local_interpolator_t
78 
79 contains
80 
83  subroutine local_interpolator_init(this, Xh, r, s, t, n_points)
84  class(local_interpolator_t), intent(inout), target :: this
85  type(space_t), intent(in), target :: Xh
86  integer, intent(in) :: n_points
87  real(kind=rp) :: r(n_points), s(n_points), t(n_points)
88  integer :: size_weights
89  call this%free()
90  if ((xh%t .eq. gl) .or. (xh%t .eq. gll)) then
91  else
92  call neko_error('Unsupported interpolation')
93  end if
94 
95  this%Xh => xh
96  this%n_points = n_points
97  allocate(this%weights_r(xh%lx,n_points))
98  allocate(this%weights_s(xh%ly,n_points))
99  allocate(this%weights_t(xh%lz,n_points))
100  call this%compute_weights(r, s, t)
101  size_weights = xh%lx * n_points
102 
103  if (neko_bcknd_device .eq. 1) then
104  call device_map(this%weights_r, this%weights_r_d, size_weights)
105  call device_map(this%weights_s, this%weights_s_d, size_weights)
106  call device_map(this%weights_t, this%weights_t_d, size_weights)
107  call device_memcpy(this%weights_r, this%weights_r_d,&
108  size_weights, host_to_device, sync = .true.)
109  call device_memcpy(this%weights_s, this%weights_s_d,&
110  size_weights, host_to_device, sync = .true.)
111  call device_memcpy(this%weights_t, this%weights_t_d,&
112  size_weights, host_to_device, sync = .true.)
113  end if
114 
115  end subroutine local_interpolator_init
116 
118  subroutine local_interpolator_free(this)
119  class(local_interpolator_t), intent(inout) :: this
120 
121  if (associated(this%Xh)) this%Xh => null()
122 
123  if(allocated(this%weights_r)) deallocate(this%weights_r)
124  if(allocated(this%weights_s)) deallocate(this%weights_s)
125  if(allocated(this%weights_t)) deallocate(this%weights_t)
126  if (c_associated(this%weights_r_d)) then
127  call device_free(this%weights_r_d)
128  end if
129  if (c_associated(this%weights_s_d)) then
130  call device_free(this%weights_s_d)
131  end if
132  if (c_associated(this%weights_t_d)) then
133  call device_free(this%weights_t_d)
134  end if
135  end subroutine local_interpolator_free
136 
147  subroutine local_interpolator_compute_weights(this, r, s, t)
148  class(local_interpolator_t), intent(inout) :: this
149  real(kind=rp), intent(in) :: r(:), s(:), t(:)
150 
151  integer :: N, i, lx
152  lx = this%Xh%lx
153  n = size(r)
154 
155  do i = 1, n
156  call fd_weights_full(r(i), this%Xh%zg(:,1), lx-1, 0, this%weights_r(:,i))
157  call fd_weights_full(s(i), this%Xh%zg(:,2), lx-1, 0, this%weights_s(:,i))
158  call fd_weights_full(t(i), this%Xh%zg(:,3), lx-1, 0, this%weights_t(:,i))
159  end do
160 
162 
175  subroutine local_interpolator_evaluate(this, interp_values, el_list, field,nel)
176  class(local_interpolator_t), intent(inout) :: this
177  integer, intent(in) :: el_list(this%n_points)
178  integer, intent(in) :: nel
179  real(kind=rp), intent(inout) :: interp_values(this%n_points)
180  real(kind=rp), intent(inout) :: field(this%Xh%lxyz, nel)
181 
182 
183  call tnsr3d_el_list(interp_values, 1, field, this%Xh%lx, &
184  this%weights_r, this%weights_s, this%weights_t, el_list, this%n_points)
185 
186  end subroutine local_interpolator_evaluate
187 
188 
189 end module local_interpolation
Map a Fortran array to a device (allocate and associate)
Definition: device.F90:57
Copy data between host and device (or device and device)
Definition: device.F90:51
subroutine, public device_rzero(a_d, n)
Zero a real vector.
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:185
Fast diagonalization methods from NEKTON.
Definition: fast3d.f90:61
subroutine, public fd_weights_full(xi, x, n, m, c)
Compute finite-difference stencil weights for evaluating derivatives up to order at a point.
Definition: fast3d.f90:105
subroutine, public setup_intp(jh, jht, z_to, z_from, n_to, n_from, derivative)
Compute interpolation weights for points z_to using values at points z_from.
Definition: fast3d.f90:243
Defines a field.
Definition: field.f90:34
Routines to obtain interpolated values on a set of points with known rst coordinates in elements loca...
subroutine local_interpolator_compute_weights(this, r, s, t)
Computes interpolation weights for a list of points.
subroutine local_interpolator_evaluate(this, interp_values, el_list, field, nel)
Interpolates a list of fields based on a set of element ids.
subroutine local_interpolator_free(this)
Free pointers.
subroutine local_interpolator_init(this, Xh, r, s, t, n_points)
Initialization of point interpolation.
Definition: math.f90:60
Build configurations.
Definition: neko_config.f90:34
integer, parameter neko_bcknd_device
Definition: neko_config.f90:44
integer, parameter, public rp
Global precision used in computations.
Definition: num_types.f90:12
Implements a point.
Definition: point.f90:35
Defines a function space.
Definition: space.f90:34
integer, parameter, public gll
Definition: space.f90:48
integer, parameter, public gl
Definition: space.f90:48
Tensor operations.
Definition: tensor.f90:61
subroutine, public tnsr3d_el_list(v, nv, u, nu, A, Bt, Ct, el_list, n_pt)
Tensor product performed on a subset of the elements.
Definition: tensor.f90:188
Utilities.
Definition: utils.f90:35
field_list_t, To be able to group fields together
Definition: field_list.f90:13
Interpolation on a set of points with known rst coordinates in elements local to this process....
A point in with coordinates .
Definition: point.f90:43
The function space for the SEM solution fields.
Definition: space.f90:62