Neko 0.9.99
A portable framework for high-order spectral element flow simulations
Loading...
Searching...
No Matches
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!
34module 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
97contains
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(in) :: Xh
175 type(coef_t), intent(in) :: 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(in) :: Xh
298 type(coef_t), intent(in) :: 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
388end 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::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
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
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:544
math::sub2
subroutine, public sub2(a, b, n)
Vector substraction .
Definition math.f90:628
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