Neko  0.9.99
A portable framework for high-order spectral element flow simulations
adv_dealias.f90
Go to the documentation of this file.
1 ! Copyright (c) 2021-2024, 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 !
34 module adv_dealias
35  use advection, only: advection_t
36  use num_types, only: rp
37  use math, only: vdot3, sub2
38  use space, only: space_t, gl
39  use field, only: field_t
40  use coefs, only: coef_t
41  use device_math, only: device_vdot3, device_sub2
42  use neko_config, only: neko_bcknd_device, neko_bcknd_sx, neko_bcknd_xsmm, &
43  neko_bcknd_opencl, neko_bcknd_cuda, neko_bcknd_hip
44  use operators, only: opgrad
45  use interpolation, only: interpolator_t
46  use device, only: device_map, device_get_ptr
47  use, intrinsic :: iso_c_binding, only: c_ptr, c_null_ptr
48  implicit none
49  private
50 
52  type, public, extends(advection_t) :: adv_dealias_t
54  type(coef_t) :: coef_gl
56  type(coef_t), pointer :: coef_gll
58  type(interpolator_t) :: gll_to_gl
60  type(space_t) :: xh_gl
62  type(space_t), pointer :: xh_gll
63  real(kind=rp), allocatable :: temp(:), tbf(:)
65  real(kind=rp), allocatable :: tx(:), ty(:), tz(:)
66  real(kind=rp), allocatable :: vr(:), vs(:), vt(:)
68  type(c_ptr) :: temp_d = c_null_ptr
70  type(c_ptr) :: tbf_d = c_null_ptr
72  type(c_ptr) :: tx_d = c_null_ptr
74  type(c_ptr) :: ty_d = c_null_ptr
76  type(c_ptr) :: tz_d = c_null_ptr
78  type(c_ptr) :: vr_d = c_null_ptr
80  type(c_ptr) :: vs_d = c_null_ptr
82  type(c_ptr) :: vt_d = c_null_ptr
83 
84  contains
87  procedure, pass(this) :: compute => compute_advection_dealias
90  procedure, pass(this) :: compute_scalar => compute_scalar_advection_dealias
92  procedure, pass(this) :: init => init_dealias
94  procedure, pass(this) :: free => free_dealias
95  end type adv_dealias_t
96 
97 contains
98 
102  subroutine init_dealias(this, lxd, coef)
103  class(adv_dealias_t), target, intent(inout) :: this
104  integer, intent(in) :: lxd
105  type(coef_t), intent(inout), target :: coef
106  integer :: nel, n_GL, n
107 
108  call this%Xh_GL%init(gl, lxd, lxd, lxd)
109  this%Xh_GLL => coef%Xh
110  this%coef_GLL => coef
111  call this%GLL_to_GL%init(this%Xh_GL, this%Xh_GLL)
112 
113  call this%coef_GL%init(this%Xh_GL, coef%msh)
114 
115  nel = coef%msh%nelv
116  n_gl = nel*this%Xh_GL%lxyz
117  n = nel*coef%Xh%lxyz
118  call this%GLL_to_GL%map(this%coef_GL%drdx, coef%drdx, nel, this%Xh_GL)
119  call this%GLL_to_GL%map(this%coef_GL%dsdx, coef%dsdx, nel, this%Xh_GL)
120  call this%GLL_to_GL%map(this%coef_GL%dtdx, coef%dtdx, nel, this%Xh_GL)
121  call this%GLL_to_GL%map(this%coef_GL%drdy, coef%drdy, nel, this%Xh_GL)
122  call this%GLL_to_GL%map(this%coef_GL%dsdy, coef%dsdy, nel, this%Xh_GL)
123  call this%GLL_to_GL%map(this%coef_GL%dtdy, coef%dtdy, nel, this%Xh_GL)
124  call this%GLL_to_GL%map(this%coef_GL%drdz, coef%drdz, nel, this%Xh_GL)
125  call this%GLL_to_GL%map(this%coef_GL%dsdz, coef%dsdz, nel, this%Xh_GL)
126  call this%GLL_to_GL%map(this%coef_GL%dtdz, coef%dtdz, nel, this%Xh_GL)
127  if ((neko_bcknd_hip .eq. 1) .or. (neko_bcknd_cuda .eq. 1) .or. &
128  (neko_bcknd_opencl .eq. 1) .or. (neko_bcknd_sx .eq. 1) .or. &
129  (neko_bcknd_xsmm .eq. 1)) then
130  allocate(this%temp(n_gl))
131  allocate(this%tbf(n_gl))
132  allocate(this%tx(n_gl))
133  allocate(this%ty(n_gl))
134  allocate(this%tz(n_gl))
135  allocate(this%vr(n_gl))
136  allocate(this%vs(n_gl))
137  allocate(this%vt(n_gl))
138  end if
139 
140  if (neko_bcknd_device .eq. 1) then
141  call device_map(this%temp, this%temp_d, n_gl)
142  call device_map(this%tbf, this%tbf_d, n_gl)
143  call device_map(this%tx, this%tx_d, n_gl)
144  call device_map(this%ty, this%ty_d, n_gl)
145  call device_map(this%tz, this%tz_d, n_gl)
146  call device_map(this%vr, this%vr_d, n_gl)
147  call device_map(this%vs, this%vs_d, n_gl)
148  call device_map(this%vt, this%vt_d, n_gl)
149  end if
150 
151  end subroutine init_dealias
152 
154  subroutine free_dealias(this)
155  class(adv_dealias_t), intent(inout) :: this
156  end subroutine free_dealias
157 
158 
171  subroutine compute_advection_dealias(this, vx, vy, vz, fx, fy, fz, Xh, &
172  coef, n, dt)
173  class(adv_dealias_t), intent(inout) :: this
174  type(space_t), intent(inout) :: Xh
175  type(coef_t), intent(inout) :: coef
176  type(field_t), intent(inout) :: vx, vy, vz
177  type(field_t), intent(inout) :: fx, fy, fz
178  integer, intent(in) :: n
179  real(kind=rp), intent(in), optional :: dt
180 
181  real(kind=rp), dimension(this%Xh_GL%lxyz) :: tx, ty, tz
182  real(kind=rp), dimension(this%Xh_GL%lxyz) :: tfx, tfy, tfz
183  real(kind=rp), dimension(this%Xh_GL%lxyz) :: vr, vs, vt
184  real(kind=rp), dimension(this%Xh_GLL%lxyz) :: tempx, tempy, tempz
185  integer :: e, i, idx, nel, n_GL
186 
187  nel = coef%msh%nelv
188  n_gl = nel * this%Xh_GL%lxyz
189 
190  !This is extremely primitive and unoptimized on the device //Karp
191  associate(c_gl => this%coef_GL)
192  if (neko_bcknd_device .eq. 1) then
193  call this%GLL_to_GL%map(this%tx, vx%x, nel, this%Xh_GL)
194  call this%GLL_to_GL%map(this%ty, vy%x, nel, this%Xh_GL)
195  call this%GLL_to_GL%map(this%tz, vz%x, nel, this%Xh_GL)
196 
197  call opgrad(this%vr, this%vs, this%vt, this%tx, c_gl)
198  call device_vdot3(this%tbf_d, this%vr_d, this%vs_d, this%vt_d, &
199  this%tx_d, this%ty_d, this%tz_d, n_gl)
200  call this%GLL_to_GL%map(this%temp, this%tbf, nel, this%Xh_GLL)
201  call device_sub2(fx%x_d, this%temp_d, n)
202 
203 
204  call opgrad(this%vr, this%vs, this%vt, this%ty, c_gl)
205  call device_vdot3(this%tbf_d, this%vr_d, this%vs_d, this%vt_d, &
206  this%tx_d, this%ty_d, this%tz_d, n_gl)
207  call this%GLL_to_GL%map(this%temp, this%tbf, nel, this%Xh_GLL)
208  call device_sub2(fy%x_d, this%temp_d, n)
209 
210  call opgrad(this%vr, this%vs, this%vt, this%tz, c_gl)
211  call device_vdot3(this%tbf_d, this%vr_d, this%vs_d, this%vt_d, &
212  this%tx_d, this%ty_d, this%tz_d, n_gl)
213  call this%GLL_to_GL%map(this%temp, this%tbf, nel, this%Xh_GLL)
214  call device_sub2(fz%x_d, this%temp_d, n)
215 
216  else if ((neko_bcknd_sx .eq. 1) .or. (neko_bcknd_xsmm .eq. 1)) then
217 
218  call this%GLL_to_GL%map(this%tx, vx%x, nel, this%Xh_GL)
219  call this%GLL_to_GL%map(this%ty, vy%x, nel, this%Xh_GL)
220  call this%GLL_to_GL%map(this%tz, vz%x, nel, this%Xh_GL)
221 
222  call opgrad(this%vr, this%vs, this%vt, this%tx, c_gl)
223  call vdot3(this%tbf, this%vr, this%vs, this%vt, &
224  this%tx, this%ty, this%tz, n_gl)
225  call this%GLL_to_GL%map(this%temp, this%tbf, nel, this%Xh_GLL)
226  call sub2(fx%x, this%temp, n)
227 
228 
229  call opgrad(this%vr, this%vs, this%vt, this%ty, c_gl)
230  call vdot3(this%tbf, this%vr, this%vs, this%vt, &
231  this%tx, this%ty, this%tz, n_gl)
232  call this%GLL_to_GL%map(this%temp, this%tbf, nel, this%Xh_GLL)
233  call sub2(fy%x, this%temp, n)
234 
235  call opgrad(this%vr, this%vs, this%vt, this%tz, c_gl)
236  call vdot3(this%tbf, this%vr, this%vs, this%vt, &
237  this%tx, this%ty, this%tz, n_gl)
238  call this%GLL_to_GL%map(this%temp, this%tbf, nel, this%Xh_GLL)
239  call sub2(fz%x, this%temp, n)
240 
241  else
242 
243  do e = 1, coef%msh%nelv
244  call this%GLL_to_GL%map(tx, vx%x(1,1,1,e), 1, this%Xh_GL)
245  call this%GLL_to_GL%map(ty, vy%x(1,1,1,e), 1, this%Xh_GL)
246  call this%GLL_to_GL%map(tz, vz%x(1,1,1,e), 1, this%Xh_GL)
247 
248  call opgrad(vr, vs, vt, tx, c_gl, e, e)
249  do i = 1, this%Xh_GL%lxyz
250  tfx(i) = tx(i)*vr(i) + ty(i)*vs(i) + tz(i)*vt(i)
251  end do
252 
253  call opgrad(vr, vs, vt, ty, c_gl, e, e)
254  do i = 1, this%Xh_GL%lxyz
255  tfy(i) = tx(i)*vr(i) + ty(i)*vs(i) + tz(i)*vt(i)
256  end do
257 
258  call opgrad(vr, vs, vt, tz, c_gl, e, e)
259  do i = 1, this%Xh_GL%lxyz
260  tfz(i) = tx(i)*vr(i) + ty(i)*vs(i) + tz(i)*vt(i)
261  end do
262 
263  call this%GLL_to_GL%map(tempx, tfx, 1, this%Xh_GLL)
264  call this%GLL_to_GL%map(tempy, tfy, 1, this%Xh_GLL)
265  call this%GLL_to_GL%map(tempz, tfz, 1, this%Xh_GLL)
266 
267  idx = (e-1)*this%Xh_GLL%lxyz+1
268  do concurrent(i = 0:this%Xh_GLL%lxyz-1)
269  fx%x(i+idx,1,1,1) = fx%x(i+idx,1,1,1) - tempx(i+1)
270  fy%x(i+idx,1,1,1) = fy%x(i+idx,1,1,1) - tempy(i+1)
271  fz%x(i+idx,1,1,1) = fz%x(i+idx,1,1,1) - tempz(i+1)
272  end do
273  end do
274  end if
275  end associate
276 
277  end subroutine compute_advection_dealias
278 
291  subroutine compute_scalar_advection_dealias(this, vx, vy, vz, s, fs, Xh, &
292  coef, n, dt)
293  class(adv_dealias_t), intent(inout) :: this
294  type(field_t), intent(inout) :: vx, vy, vz
295  type(field_t), intent(inout) :: s
296  type(field_t), intent(inout) :: fs
297  type(space_t), intent(inout) :: Xh
298  type(coef_t), intent(inout) :: coef
299  integer, intent(in) :: n
300  real(kind=rp), intent(in), optional :: dt
301 
302  real(kind=rp), dimension(this%Xh_GL%lxyz) :: vx_gl, vy_gl, vz_gl, s_gl
303  real(kind=rp), dimension(this%Xh_GL%lxyz) :: dsdx, dsdy, dsdz
304  real(kind=rp), dimension(this%Xh_GL%lxyz) :: f_gl
305  integer :: e, i, idx, nel, n_GL
306  real(kind=rp), dimension(this%Xh_GLL%lxyz) :: temp
307 
308  nel = coef%msh%nelv
309  n_gl = nel * this%Xh_GL%lxyz
310 
311  associate(c_gl => this%coef_GL)
312  if (neko_bcknd_device .eq. 1) then
313 
314  ! Map advecting velocity onto the higher-order space
315  call this%GLL_to_GL%map(this%tx, vx%x, nel, this%Xh_GL)
316  call this%GLL_to_GL%map(this%ty, vy%x, nel, this%Xh_GL)
317  call this%GLL_to_GL%map(this%tz, vz%x, nel, this%Xh_GL)
318 
319  ! Map the scalar onto the high-order space
320  call this%GLL_to_GL%map(this%temp, s%x, nel, this%Xh_GL)
321 
322  ! Compute the scalar gradient in the high-order space
323  call opgrad(this%vr, this%vs, this%vt, this%temp, c_gl)
324 
325  ! Compute the convective term, i.e dot the velocity with the scalar grad
326  call device_vdot3(this%tbf_d, this%vr_d, this%vs_d, this%vt_d, &
327  this%tx_d, this%ty_d, this%tz_d, n_gl)
328 
329  ! Map back to the original space (we reuse this%temp)
330  call this%GLL_to_GL%map(this%temp, this%tbf, nel, this%Xh_GLL)
331 
332  ! Update the source term
333  call device_sub2(fs%x_d, this%temp_d, n)
334 
335  else if ((neko_bcknd_sx .eq. 1) .or. (neko_bcknd_xsmm .eq. 1)) then
336 
337  ! Map advecting velocity onto the higher-order space
338  call this%GLL_to_GL%map(this%tx, vx%x, nel, this%Xh_GL)
339  call this%GLL_to_GL%map(this%ty, vy%x, nel, this%Xh_GL)
340  call this%GLL_to_GL%map(this%tz, vz%x, nel, this%Xh_GL)
341 
342  ! Map the scalar onto the high-order space
343  call this%GLL_to_GL%map(this%temp, s%x, nel, this%Xh_GL)
344 
345  ! Compute the scalar gradient in the high-order space
346  call opgrad(this%vr, this%vs, this%vt, this%temp, c_gl)
347 
348  ! Compute the convective term, i.e dot the velocity with the scalar grad
349  call vdot3(this%tbf, this%vr, this%vs, this%vt, &
350  this%tx, this%ty, this%tz, n_gl)
351 
352  ! Map back to the original space (we reuse this%temp)
353  call this%GLL_to_GL%map(this%temp, this%tbf, nel, this%Xh_GLL)
354 
355  ! Update the source term
356  call sub2(fs%x, this%temp, n)
357 
358  else
359  do e = 1, coef%msh%nelv
360  ! Map advecting velocity onto the higher-order space
361  call this%GLL_to_GL%map(vx_gl, vx%x(1,1,1,e), 1, this%Xh_GL)
362  call this%GLL_to_GL%map(vy_gl, vy%x(1,1,1,e), 1, this%Xh_GL)
363  call this%GLL_to_GL%map(vz_gl, vz%x(1,1,1,e), 1, this%Xh_GL)
364 
365  ! Map scalar onto the higher-order space
366  call this%GLL_to_GL%map(s_gl, s%x(1,1,1,e), 1, this%Xh_GL)
367 
368  ! Gradient of s in the higher-order space
369  call opgrad(dsdx, dsdy, dsdz, s_gl, c_gl, e, e)
370 
371  ! vx * ds/dx + vy * ds/dy + vz * ds/dz for each point in the element
372  do i = 1, this%Xh_GL%lxyz
373  f_gl(i) = vx_gl(i)*dsdx(i) + vy_gl(i)*dsdy(i) + vz_gl(i)*dsdz(i)
374  end do
375 
376  ! Map back the contructed operator to the original space
377  call this%GLL_to_GL%map(temp, f_gl, 1, this%Xh_GLL)
378 
379  idx = (e-1)*this%Xh_GLL%lxyz + 1
380 
381  call sub2(fs%x(idx, 1, 1, 1), temp, this%Xh_GLL%lxyz)
382  end do
383  end if
384  end associate
385 
386  end subroutine compute_scalar_advection_dealias
387 
388 end module adv_dealias
device::device_get_ptr
Return the device pointer for an associated Fortran array.
Definition: device.F90:81
device::device_map
Map a Fortran array to a device (allocate and associate)
Definition: device.F90:57
adv_dealias
Subroutines to add advection terms to the RHS of a transport equation.
Definition: adv_dealias.f90:34
adv_dealias::compute_advection_dealias
subroutine compute_advection_dealias(this, vx, vy, vz, fx, fy, fz, Xh, coef, n, dt)
Add the advection term for the fluid, i.e. , to the RHS.
Definition: adv_dealias.f90:173
adv_dealias::init_dealias
subroutine init_dealias(this, lxd, coef)
Constructor.
Definition: adv_dealias.f90:103
adv_dealias::compute_scalar_advection_dealias
subroutine compute_scalar_advection_dealias(this, vx, vy, vz, s, fs, Xh, coef, n, dt)
Add the advection term for a scalar, i.e. , to the RHS.
Definition: adv_dealias.f90:293
adv_dealias::free_dealias
subroutine free_dealias(this)
Destructor.
Definition: adv_dealias.f90:155
advection
Subroutines to add advection terms to the RHS of a transport equation.
Definition: advection.f90:34
coefs
Coefficients.
Definition: coef.f90:34
device_math
Definition: device_math.F90:33
device_math::device_vdot3
subroutine, public device_vdot3(dot_d, u1_d, u2_d, u3_d, v1_d, v2_d, v3_d, n)
Compute a dot product (3-d version) assuming vector components etc.
Definition: device_math.F90:491
device_math::device_sub2
subroutine, public device_sub2(a_d, b_d, n)
Vector substraction .
Definition: device_math.F90:430
device
Device abstraction, common interface for various accelerators.
Definition: device.F90:34
field
Defines a field.
Definition: field.f90:34
interpolation
Routines to interpolate between different spaces.
Definition: interpolation.f90:34
math
Definition: math.f90:60
math::vdot3
subroutine, public vdot3(dot, u1, u2, u3, v1, v2, v3, n)
Compute a dot product (3-d version) assuming vector components etc.
Definition: math.f90:545
math::sub2
subroutine, public sub2(a, b, n)
Vector substraction .
Definition: math.f90:629
neko_config
Build configurations.
Definition: neko_config.f90:34
neko_config::neko_bcknd_sx
integer, parameter neko_bcknd_sx
Definition: neko_config.f90:39
neko_config::neko_bcknd_hip
integer, parameter neko_bcknd_hip
Definition: neko_config.f90:42
neko_config::neko_bcknd_device
integer, parameter neko_bcknd_device
Definition: neko_config.f90:44
neko_config::neko_bcknd_opencl
integer, parameter neko_bcknd_opencl
Definition: neko_config.f90:43
neko_config::neko_bcknd_cuda
integer, parameter neko_bcknd_cuda
Definition: neko_config.f90:41
neko_config::neko_bcknd_xsmm
integer, parameter neko_bcknd_xsmm
Definition: neko_config.f90:40
num_types
Definition: num_types.f90:1
num_types::rp
integer, parameter, public rp
Global precision used in computations.
Definition: num_types.f90:12
operators
Operators.
Definition: operators.f90:34
operators::opgrad
subroutine, public opgrad(ux, uy, uz, u, coef, es, ee)
Compute the weak gradient of a scalar field, i.e. the gradient multiplied by the mass matrix.
Definition: operators.f90:171
space
Defines a function space.
Definition: space.f90:34
space::gl
integer, parameter, public gl
Definition: space.f90:48
adv_dealias::adv_dealias_t
Type encapsulating advection routines with dealiasing.
Definition: adv_dealias.f90:52
advection::advection_t
Base abstract type for computing the advection operator.
Definition: advection.f90:46
coefs::coef_t
Coefficients defined on a given (mesh, ) tuple. Arrays use indices (i,j,k,e): element e,...
Definition: coef.f90:55
field::field_t
Definition: field.f90:46
interpolation::interpolator_t
Interpolation between two space::space_t.
Definition: interpolation.f90:51
space::space_t
The function space for the SEM solution fields.
Definition: space.f90:62