Neko  0.9.0
A portable framework for high-order spectral element flow simulations
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 !
35  use neko_config
36  use num_types, only : rp
37  use device
38  use fast3d
39  use tensor, only : tnsr3d
40  use tensor_cpu, only : tnsr3d_cpu
41  use space, only : space_t, operator(.eq.), gl, gll
42  implicit none
43  private
44 
51  type, public :: interpolator_t
53  type(space_t), pointer :: xh
55  type(space_t), pointer :: yh
57  real(kind=rp), allocatable :: xh_to_yh(:,:), xh_to_yht(:,:)
59  real(kind=rp), allocatable :: yh_to_xh(:,:), yh_to_xht(:,:)
61  type(c_ptr) :: xh_yh_d = c_null_ptr
63  type(c_ptr) :: xh_yht_d = c_null_ptr
65  type(c_ptr) :: yh_xh_d = c_null_ptr
67  type(c_ptr) :: yh_xht_d = c_null_ptr
68  contains
69 
71  procedure, pass(this) :: init => interpolator_init
73  procedure, pass(this) :: free => interpolator_free
75  procedure, pass(this) :: map => interpolator_map
77  procedure, pass(this) :: map_host => interpolator_map_host
78 
79  end type interpolator_t
80 
81 contains
82 
86  subroutine interpolator_init(this, Xh, Yh)
87  class(interpolator_t), intent(inout), target :: this
88  type(space_t), intent(inout), target :: Xh
89  type(space_t), intent(inout), target :: Yh
90  integer :: deg_derivate
91 
92  call this%free()
93 
94  allocate(this%Xh_to_Yh(yh%lx,xh%lx))
95  allocate(this%Xh_to_YhT(xh%lx,yh%lx))
96  allocate(this%Yh_to_Xh(xh%lx,yh%lx))
97  allocate(this%Yh_to_XhT(yh%lx,xh%lx))
98 
99  if (xh%t .eq. gll .and. yh%t .eq. gll) then
100  else if ((xh%t .eq. gl .and. yh%t .eq. gll) .or. &
101  (yh%t .eq. gl .and. xh%t .eq. gll)) then
102  else
103  call neko_error('Unsupported interpolation')
104  end if
105  deg_derivate = 0
106  call setup_intp(this%Xh_to_Yh, this%Xh_to_YhT, &
107  yh%zg, xh%zg, yh%lx, xh%lx, deg_derivate)
108  call setup_intp(this%Yh_to_Xh, this%Yh_to_XhT, &
109  xh%zg, yh%zg, xh%lx, yh%lx, deg_derivate)
110 
111  this%Xh => xh
112  this%Yh => yh
113  if (neko_bcknd_device .eq. 1) then
114  call device_map(this%Xh_to_Yh, this%Xh_Yh_d, yh%lx*xh%lx)
115  call device_map(this%Xh_to_YhT, this%Xh_YhT_d, yh%lx*xh%lx)
116  call device_map(this%Yh_to_Xh, this%Yh_Xh_d, yh%lx*xh%lx)
117  call device_map(this%Yh_to_XhT, this%Yh_XhT_d, yh%lx*xh%lx)
118  call device_memcpy(this%Xh_to_Yh, this%Xh_Yh_d, yh%lx*xh%lx, &
119  host_to_device, sync=.false.)
120  call device_memcpy(this%Xh_to_YhT, this%Xh_YhT_d, yh%lx*xh%lx, &
121  host_to_device, sync=.false.)
122  call device_memcpy(this%Yh_to_Xh, this%Yh_Xh_d, yh%lx*xh%lx, &
123  host_to_device, sync=.false.)
124  call device_memcpy(this%Yh_to_XhT, this%Yh_XhT_d, yh%lx*xh%lx, &
125  host_to_device, sync=.false.)
126  end if
127 
128  end subroutine interpolator_init
129 
130  subroutine interpolator_free(this)
131  class(interpolator_t), intent(inout) :: this
132 
133  if (allocated(this%Xh_to_Yh)) then
134  deallocate(this%Xh_to_Yh)
135  end if
136  if (allocated(this%Xh_to_YhT)) then
137  deallocate(this%Xh_to_YhT)
138  end if
139  if (allocated(this%Yh_to_Xh)) then
140  deallocate(this%Yh_to_Xh)
141  end if
142  if (allocated(this%Yh_to_XhT)) then
143  deallocate(this%Yh_to_XhT)
144  end if
145  if (c_associated(this%Yh_Xh_d)) then
146  call device_free(this%Yh_Xh_d)
147  end if
148  if (c_associated(this%Yh_XhT_d)) then
149  call device_free(this%Yh_XhT_d)
150  end if
151  if (c_associated(this%Xh_Yh_d)) then
152  call device_free(this%Xh_Yh_d)
153  end if
154  if (c_associated(this%Xh_YhT_d)) then
155  call device_free(this%Xh_YhT_d)
156  end if
157 
158  end subroutine interpolator_free
159 
165  subroutine interpolator_map(this, y, x, nel, to_space)
166  class(interpolator_t), intent(inout) :: this
167  integer :: nel
168  type(space_t) :: to_space
169  real(kind=rp), intent(inout) :: x(this%Xh%lx, this%Xh%lx, this%Xh%lx, nel)
170  real(kind=rp), intent(inout) :: y(this%Yh%lx, this%Yh%lx, this%Yh%lx, nel)
171  if (to_space .eq. this%Yh) then
172  call tnsr3d(y, this%Yh%lx, x, &
173  this%Xh%lx,this%Yh_to_XhT, &
174  this%Yh_to_Xh, this%Yh_to_Xh, nel)
175  else if (to_space .eq. this%Xh) then
176  call tnsr3d(y, this%Xh%lx, x, &
177  this%Yh%lx,this%Yh_to_Xh, &
178  this%Yh_to_XhT, this%Yh_to_XhT, nel)
179  else
180  call neko_error('Invalid interpolation')
181  end if
182  end subroutine interpolator_map
183 
184 
191  subroutine interpolator_map_host(this, y, x, nel, to_space)
192  class(interpolator_t), intent(inout) :: this
193  integer :: nel
194  type(space_t) :: to_space
195  real(kind=rp), intent(inout) :: x(this%Xh%lx, this%Xh%lx, this%Xh%lx, nel)
196  real(kind=rp), intent(inout) :: y(this%Yh%lx, this%Yh%lx, this%Yh%lx, nel)
197  if (to_space .eq. this%Yh) then
198  call tnsr3d_cpu(y, this%Yh%lx, x, &
199  this%Xh%lx,this%Yh_to_XhT, &
200  this%Yh_to_Xh, this%Yh_to_Xh, nel)
201  else if (to_space .eq. this%Xh) then
202  call tnsr3d_cpu(y, this%Xh%lx, x, &
203  this%Yh%lx,this%Yh_to_Xh, &
204  this%Yh_to_XhT, this%Yh_to_XhT, nel)
205  else
206  call neko_error('Invalid interpolation')
207  end if
208  end subroutine interpolator_map_host
209 
210 
211 end module 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
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 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
Routines to interpolate between different spaces.
subroutine interpolator_map_host(this, y, x, nel, to_space)
Interpolates an array to one of Xh or Yh on host.
subroutine interpolator_map(this, y, x, nel, to_space)
Interpolates an array to one of Xh or Yh.
subroutine interpolator_init(this, Xh, Yh)
Constructor to initialize with two different spaces.
subroutine interpolator_free(this)
NEKTON map.
Definition: map.f90:3
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
Defines a function space.
Definition: space.f90:34
integer, parameter, public gll
Definition: space.f90:48
integer, parameter, public gl
Definition: space.f90:48
subroutine, public tnsr3d_cpu(v, nv, u, nu, A, Bt, Ct, nelv)
Definition: tensor_cpu.f90:878
Tensor operations.
Definition: tensor.f90:61
subroutine, public tnsr3d(v, nv, u, nu, A, Bt, Ct, nelv)
Tensor product performed on nelv elements.
Definition: tensor.f90:223
Interpolation between two space::space_t.
The function space for the SEM solution fields.
Definition: space.f90:62