Neko 1.99.5
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-2026, 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, device_col3, device_add2
42 use neko_config, only : neko_bcknd_device, neko_bcknd_sx, neko_bcknd_xsmm, &
43 neko_bcknd_opencl, neko_bcknd_cuda, neko_bcknd_hip, neko_bcknd_metal
44 use utils, only : neko_error
45 use operators, only : opgrad
46 use interpolation, only : interpolator_t
47 use device, only : device_map, device_get_ptr, device_unmap
48 use, intrinsic :: iso_c_binding, only : c_ptr, c_null_ptr
49 implicit none
50 private
51
53 type, public, extends(advection_t) :: adv_dealias_t
55 type(coef_t) :: coef_gl
57 type(coef_t), pointer :: coef_gll
59 type(interpolator_t) :: gll_to_gl
61 type(space_t) :: xh_gl
63 type(space_t), pointer :: xh_gll
64 real(kind=rp), allocatable :: temp(:), tbf(:)
66 real(kind=rp), allocatable :: tx(:), ty(:), tz(:)
67 real(kind=rp), allocatable :: vr(:), vs(:), vt(:)
69 type(c_ptr) :: temp_d = c_null_ptr
71 type(c_ptr) :: tbf_d = c_null_ptr
73 type(c_ptr) :: tx_d = c_null_ptr
75 type(c_ptr) :: ty_d = c_null_ptr
77 type(c_ptr) :: tz_d = c_null_ptr
79 type(c_ptr) :: vr_d = c_null_ptr
81 type(c_ptr) :: vs_d = c_null_ptr
83 type(c_ptr) :: vt_d = c_null_ptr
84
85 contains
88 procedure, pass(this) :: compute => compute_advection_dealias
91 procedure, pass(this) :: compute_scalar => compute_scalar_advection_dealias
93 procedure, pass(this) :: compute_ale => compute_ale_advection_dealias
95 procedure, pass(this) :: recompute_metrics => recompute_metrics_dealias
97 procedure, pass(this) :: init => init_dealias
99 procedure, pass(this) :: free => free_dealias
100 end type adv_dealias_t
101
102contains
103
107 subroutine init_dealias(this, lxd, coef)
108 class(adv_dealias_t), target, intent(inout) :: this
109 integer, intent(in) :: lxd
110 type(coef_t), intent(inout), target :: coef
111 integer :: nel, n_GL, n
112
113 call this%Xh_GL%init(gl, lxd, lxd, lxd)
114 this%Xh_GLL => coef%Xh
115 this%coef_GLL => coef
116 call this%GLL_to_GL%init(this%Xh_GL, this%Xh_GLL)
117
118 call this%coef_GL%init(this%Xh_GL, coef%msh)
119
120 nel = coef%msh%nelv
121 n_gl = nel*this%Xh_GL%lxyz
122 n = nel*coef%Xh%lxyz
123 call this%GLL_to_GL%map(this%coef_GL%drdx, coef%drdx, nel, this%Xh_GL)
124 call this%GLL_to_GL%map(this%coef_GL%dsdx, coef%dsdx, nel, this%Xh_GL)
125 call this%GLL_to_GL%map(this%coef_GL%dtdx, coef%dtdx, nel, this%Xh_GL)
126 call this%GLL_to_GL%map(this%coef_GL%drdy, coef%drdy, nel, this%Xh_GL)
127 call this%GLL_to_GL%map(this%coef_GL%dsdy, coef%dsdy, nel, this%Xh_GL)
128 call this%GLL_to_GL%map(this%coef_GL%dtdy, coef%dtdy, nel, this%Xh_GL)
129 call this%GLL_to_GL%map(this%coef_GL%drdz, coef%drdz, nel, this%Xh_GL)
130 call this%GLL_to_GL%map(this%coef_GL%dsdz, coef%dsdz, nel, this%Xh_GL)
131 call this%GLL_to_GL%map(this%coef_GL%dtdz, coef%dtdz, nel, this%Xh_GL)
132 if ((neko_bcknd_hip .eq. 1) .or. (neko_bcknd_cuda .eq. 1) .or. &
133 (neko_bcknd_opencl .eq. 1) .or. (neko_bcknd_metal .eq. 1) .or. &
134 (neko_bcknd_sx .eq. 1) .or. (neko_bcknd_xsmm .eq. 1)) then
135 allocate(this%temp(n_gl))
136 allocate(this%tbf(n_gl))
137 allocate(this%tx(n_gl))
138 allocate(this%ty(n_gl))
139 allocate(this%tz(n_gl))
140 allocate(this%vr(n_gl))
141 allocate(this%vs(n_gl))
142 allocate(this%vt(n_gl))
143 end if
144
145 if (neko_bcknd_device .eq. 1) then
146 call device_map(this%temp, this%temp_d, n_gl)
147 call device_map(this%tbf, this%tbf_d, n_gl)
148 call device_map(this%tx, this%tx_d, n_gl)
149 call device_map(this%ty, this%ty_d, n_gl)
150 call device_map(this%tz, this%tz_d, n_gl)
151 call device_map(this%vr, this%vr_d, n_gl)
152 call device_map(this%vs, this%vs_d, n_gl)
153 call device_map(this%vt, this%vt_d, n_gl)
154 end if
155
156 end subroutine init_dealias
157
159 subroutine free_dealias(this)
160 class(adv_dealias_t), intent(inout) :: this
161
162 if (allocated(this%temp)) then
163 if (neko_bcknd_device .eq. 1) then
164 call device_unmap(this%temp, this%temp_d)
165 end if
166 deallocate(this%temp)
167 end if
168
169 if (allocated(this%tbf)) then
170 if (neko_bcknd_device .eq. 1) then
171 call device_unmap(this%tbf, this%tbf_d)
172 end if
173 deallocate(this%tbf)
174 end if
175 if (allocated(this%tx)) then
176 if (neko_bcknd_device .eq. 1) then
177 call device_unmap(this%tx, this%tx_d)
178 end if
179 deallocate(this%tx)
180 end if
181 if (allocated(this%ty)) then
182 if (neko_bcknd_device .eq. 1) then
183 call device_unmap(this%ty, this%ty_d)
184 end if
185 deallocate(this%ty)
186 end if
187 if (allocated(this%tz)) then
188 if (neko_bcknd_device .eq. 1) then
189 call device_unmap(this%tz, this%tz_d)
190 end if
191 deallocate(this%tz)
192 end if
193 if (allocated(this%vr)) then
194 if (neko_bcknd_device .eq. 1) then
195 call device_unmap(this%vr, this%vr_d)
196 end if
197 deallocate(this%vr)
198 end if
199 if (allocated(this%vs)) then
200 if (neko_bcknd_device .eq. 1) then
201 call device_unmap(this%vs, this%vs_d)
202 end if
203 deallocate(this%vs)
204 end if
205 if (allocated(this%vt)) then
206 if (neko_bcknd_device .eq. 1) then
207 call device_unmap(this%vt, this%vt_d)
208 end if
209 deallocate(this%vt)
210 end if
211
212 call this%coef_GL%free()
213 call this%GLL_to_GL%free()
214 call this%Xh_GL%free()
215
216 nullify(this%Xh_GLL)
217 nullify(this%coef_GLL)
218
219 end subroutine free_dealias
220
221
234 subroutine compute_advection_dealias(this, vx, vy, vz, fx, fy, fz, Xh, &
235 coef, n, dt)
236 class(adv_dealias_t), intent(inout) :: this
237 type(space_t), intent(in) :: Xh
238 type(coef_t), intent(in) :: coef
239 type(field_t), intent(inout) :: vx, vy, vz
240 type(field_t), intent(inout) :: fx, fy, fz
241 integer, intent(in) :: n
242 real(kind=rp), intent(in), optional :: dt
243
244 real(kind=rp), dimension(this%Xh_GL%lxyz) :: tx, ty, tz
245 real(kind=rp), dimension(this%Xh_GL%lxyz) :: tfx, tfy, tfz
246 real(kind=rp), dimension(this%Xh_GL%lxyz) :: vr, vs, vt
247 real(kind=rp), dimension(this%Xh_GLL%lxyz) :: tempx, tempy, tempz
248 integer :: e, i, idx, nel, n_GL
249
250 nel = coef%msh%nelv
251 n_gl = nel * this%Xh_GL%lxyz
252
253 !This is extremely primitive and unoptimized on the device //Karp
254 associate(c_gl => this%coef_GL)
255 if (neko_bcknd_device .eq. 1) then
256 call this%GLL_to_GL%map(this%tx, vx%x, nel, this%Xh_GL)
257 call this%GLL_to_GL%map(this%ty, vy%x, nel, this%Xh_GL)
258 call this%GLL_to_GL%map(this%tz, vz%x, nel, this%Xh_GL)
259
260 call opgrad(this%vr, this%vs, this%vt, this%tx, c_gl)
261 call device_vdot3(this%tbf_d, this%vr_d, this%vs_d, this%vt_d, &
262 this%tx_d, this%ty_d, this%tz_d, n_gl)
263 call this%GLL_to_GL%map(this%temp, this%tbf, nel, this%Xh_GLL)
264 call device_sub2(fx%x_d, this%temp_d, n)
265
266
267 call opgrad(this%vr, this%vs, this%vt, this%ty, c_gl)
268 call device_vdot3(this%tbf_d, this%vr_d, this%vs_d, this%vt_d, &
269 this%tx_d, this%ty_d, this%tz_d, n_gl)
270 call this%GLL_to_GL%map(this%temp, this%tbf, nel, this%Xh_GLL)
271 call device_sub2(fy%x_d, this%temp_d, n)
272
273 call opgrad(this%vr, this%vs, this%vt, this%tz, c_gl)
274 call device_vdot3(this%tbf_d, this%vr_d, this%vs_d, this%vt_d, &
275 this%tx_d, this%ty_d, this%tz_d, n_gl)
276 call this%GLL_to_GL%map(this%temp, this%tbf, nel, this%Xh_GLL)
277 call device_sub2(fz%x_d, this%temp_d, n)
278
279 else if ((neko_bcknd_sx .eq. 1) .or. (neko_bcknd_xsmm .eq. 1)) then
280
281 call this%GLL_to_GL%map(this%tx, vx%x, nel, this%Xh_GL)
282 call this%GLL_to_GL%map(this%ty, vy%x, nel, this%Xh_GL)
283 call this%GLL_to_GL%map(this%tz, vz%x, nel, this%Xh_GL)
284
285 call opgrad(this%vr, this%vs, this%vt, this%tx, c_gl)
286 call vdot3(this%tbf, this%vr, this%vs, this%vt, &
287 this%tx, this%ty, this%tz, n_gl)
288 call this%GLL_to_GL%map(this%temp, this%tbf, nel, this%Xh_GLL)
289 call sub2(fx%x, this%temp, n)
290
291
292 call opgrad(this%vr, this%vs, this%vt, this%ty, c_gl)
293 call vdot3(this%tbf, this%vr, this%vs, this%vt, &
294 this%tx, this%ty, this%tz, n_gl)
295 call this%GLL_to_GL%map(this%temp, this%tbf, nel, this%Xh_GLL)
296 call sub2(fy%x, this%temp, n)
297
298 call opgrad(this%vr, this%vs, this%vt, this%tz, c_gl)
299 call vdot3(this%tbf, this%vr, this%vs, this%vt, &
300 this%tx, this%ty, this%tz, n_gl)
301 call this%GLL_to_GL%map(this%temp, this%tbf, nel, this%Xh_GLL)
302 call sub2(fz%x, this%temp, n)
303
304 else
305 !$omp parallel do private(e, i, idx, tempx, tempy, tempz), &
306 !$omp& private(tx, ty, tz, vr, vs, vt, tfx, tfy, tfz)
307 do e = 1, coef%msh%nelv
308 call this%GLL_to_GL%map(tx, vx%x(1,1,1,e), 1, this%Xh_GL)
309 call this%GLL_to_GL%map(ty, vy%x(1,1,1,e), 1, this%Xh_GL)
310 call this%GLL_to_GL%map(tz, vz%x(1,1,1,e), 1, this%Xh_GL)
311
312 call opgrad(vr, vs, vt, tx, c_gl, e, e)
313 do i = 1, this%Xh_GL%lxyz
314 tfx(i) = tx(i)*vr(i) + ty(i)*vs(i) + tz(i)*vt(i)
315 end do
316
317 call opgrad(vr, vs, vt, ty, c_gl, e, e)
318 do i = 1, this%Xh_GL%lxyz
319 tfy(i) = tx(i)*vr(i) + ty(i)*vs(i) + tz(i)*vt(i)
320 end do
321
322 call opgrad(vr, vs, vt, tz, c_gl, e, e)
323 do i = 1, this%Xh_GL%lxyz
324 tfz(i) = tx(i)*vr(i) + ty(i)*vs(i) + tz(i)*vt(i)
325 end do
326
327 call this%GLL_to_GL%map(tempx, tfx, 1, this%Xh_GLL)
328 call this%GLL_to_GL%map(tempy, tfy, 1, this%Xh_GLL)
329 call this%GLL_to_GL%map(tempz, tfz, 1, this%Xh_GLL)
330
331 idx = (e-1)*this%Xh_GLL%lxyz+1
332 do concurrent(i = 0:this%Xh_GLL%lxyz-1)
333 fx%x(i+idx,1,1,1) = fx%x(i+idx,1,1,1) - tempx(i+1)
334 fy%x(i+idx,1,1,1) = fy%x(i+idx,1,1,1) - tempy(i+1)
335 fz%x(i+idx,1,1,1) = fz%x(i+idx,1,1,1) - tempz(i+1)
336 end do
337 end do
338 !$omp end parallel do
339 end if
340 end associate
341
342 end subroutine compute_advection_dealias
343
356 subroutine compute_scalar_advection_dealias(this, vx, vy, vz, s, fs, Xh, &
357 coef, n, dt)
358 class(adv_dealias_t), intent(inout) :: this
359 type(field_t), intent(inout) :: vx, vy, vz
360 type(field_t), intent(inout) :: s
361 type(field_t), intent(inout) :: fs
362 type(space_t), intent(in) :: Xh
363 type(coef_t), intent(in) :: coef
364 integer, intent(in) :: n
365 real(kind=rp), intent(in), optional :: dt
366
367 real(kind=rp), dimension(this%Xh_GL%lxyz) :: vx_gl, vy_gl, vz_gl, s_gl
368 real(kind=rp), dimension(this%Xh_GL%lxyz) :: dsdx, dsdy, dsdz
369 real(kind=rp), dimension(this%Xh_GL%lxyz) :: f_gl
370 integer :: e, i, idx, nel, n_GL
371 real(kind=rp), dimension(this%Xh_GLL%lxyz) :: temp
372
373 nel = coef%msh%nelv
374 n_gl = nel * this%Xh_GL%lxyz
375
376 associate(c_gl => this%coef_GL)
377 if (neko_bcknd_device .eq. 1) then
378
379 ! Map advecting velocity onto the higher-order space
380 call this%GLL_to_GL%map(this%tx, vx%x, nel, this%Xh_GL)
381 call this%GLL_to_GL%map(this%ty, vy%x, nel, this%Xh_GL)
382 call this%GLL_to_GL%map(this%tz, vz%x, nel, this%Xh_GL)
383
384 ! Map the scalar onto the high-order space
385 call this%GLL_to_GL%map(this%temp, s%x, nel, this%Xh_GL)
386
387 ! Compute the scalar gradient in the high-order space
388 call opgrad(this%vr, this%vs, this%vt, this%temp, c_gl)
389
390 ! Compute the convective term, i.e dot the velocity with the
391 ! scalar grad
392 call device_vdot3(this%tbf_d, this%vr_d, this%vs_d, this%vt_d, &
393 this%tx_d, this%ty_d, this%tz_d, n_gl)
394
395 ! Map back to the original space (we reuse this%temp)
396 call this%GLL_to_GL%map(this%temp, this%tbf, nel, this%Xh_GLL)
397
398 ! Update the source term
399 call device_sub2(fs%x_d, this%temp_d, n)
400
401 else if ((neko_bcknd_sx .eq. 1) .or. (neko_bcknd_xsmm .eq. 1)) then
402
403 ! Map advecting velocity onto the higher-order space
404 call this%GLL_to_GL%map(this%tx, vx%x, nel, this%Xh_GL)
405 call this%GLL_to_GL%map(this%ty, vy%x, nel, this%Xh_GL)
406 call this%GLL_to_GL%map(this%tz, vz%x, nel, this%Xh_GL)
407
408 ! Map the scalar onto the high-order space
409 call this%GLL_to_GL%map(this%temp, s%x, nel, this%Xh_GL)
410
411 ! Compute the scalar gradient in the high-order space
412 call opgrad(this%vr, this%vs, this%vt, this%temp, c_gl)
413
414 ! Compute the convective term, i.e dot the velocity with the
415 ! scalar grad
416 call vdot3(this%tbf, this%vr, this%vs, this%vt, &
417 this%tx, this%ty, this%tz, n_gl)
418
419 ! Map back to the original space (we reuse this%temp)
420 call this%GLL_to_GL%map(this%temp, this%tbf, nel, this%Xh_GLL)
421
422 ! Update the source term
423 call sub2(fs%x, this%temp, n)
424
425 else
426 !$omp parallel do private (e, i, idx, vx_GL, vy_GL, s_GL, f_GL, temp)
427 do e = 1, coef%msh%nelv
428 ! Map advecting velocity onto the higher-order space
429 call this%GLL_to_GL%map(vx_gl, vx%x(1,1,1,e), 1, this%Xh_GL)
430 call this%GLL_to_GL%map(vy_gl, vy%x(1,1,1,e), 1, this%Xh_GL)
431 call this%GLL_to_GL%map(vz_gl, vz%x(1,1,1,e), 1, this%Xh_GL)
432
433 ! Map scalar onto the higher-order space
434 call this%GLL_to_GL%map(s_gl, s%x(1,1,1,e), 1, this%Xh_GL)
435
436 ! Gradient of s in the higher-order space
437 call opgrad(dsdx, dsdy, dsdz, s_gl, c_gl, e, e)
438
439 ! vx * ds/dx + vy * ds/dy + vz * ds/dz for each point in the element
440 do i = 1, this%Xh_GL%lxyz
441 f_gl(i) = vx_gl(i)*dsdx(i) + vy_gl(i)*dsdy(i) + vz_gl(i)*dsdz(i)
442 end do
443
444 ! Map back the contructed operator to the original space
445 call this%GLL_to_GL%map(temp, f_gl, 1, this%Xh_GLL)
446
447 idx = (e-1)*this%Xh_GLL%lxyz + 1
448
449 call sub2(fs%x(idx, 1, 1, 1), temp, this%Xh_GLL%lxyz)
450 end do
451 !$omp end parallel do
452 end if
453 end associate
454
455 end subroutine compute_scalar_advection_dealias
456
457
458 !!> Add the advection term in ALE framework using dealiasing.
459 !! @param this The object.
460 !! @param vx The x component of velocity.
461 !! @param vy The y component of velocity.
462 !! @param vz The z component of velocity.
463 !! @param wm_x The x component of mesh velocity.
464 !! @param wm_y The y component of mesh velocity.
465 !! @param wm_z The z component of mesh velocity.
466 !! @param fx The x component of source term.
467 !! @param fy The y component of source term.
468 !! @param fz The z component of source term.
469 !! @param Xh The function space.
470 !! @param coef The coefficients of the (Xh, mesh) pair.
471 !! @param n Typically the size of the mesh.
472 !! @param dt Current time-step, not required for this method.
473 !! Here, we compute: - div ( u_i * wm ).
474 !! Based on Ho, L.W. A Legendre spectral element method for
475 !! simulation of incompressible unsteady viscous free-surface flows.
476 !! Ph.D. thesis, Massachusetts Institute of Technology, 1989.
477 !! Note: In Nek5000, dealiasing is not done for this term.
478 subroutine compute_ale_advection_dealias(this, vx, vy, vz, &
479 wm_x, wm_y, wm_z, fx, fy, fz, Xh, coef, n, dt)
480 class(adv_dealias_t), intent(inout) :: this
481 type(field_t), intent(inout) :: vx, vy, vz
482 type(field_t), intent(inout) :: wm_x, wm_y, wm_z
483 type(field_t), intent(inout) :: fx, fy, fz
484 type(space_t), intent(in) :: Xh
485 type(coef_t), intent(in) :: coef
486 integer, intent(in) :: n
487 real(kind=rp), intent(in), optional :: dt
488 real(kind=rp), dimension(this%Xh_GL%lxyz) :: vx_gl, vy_gl, vz_gl
489 real(kind=rp), dimension(this%Xh_GL%lxyz) :: wm_x_gl, wm_y_gl, wm_z_gl
490 real(kind=rp), dimension(this%Xh_GL%lxyz) :: flux_gl
491 real(kind=rp), dimension(this%Xh_GL%lxyz) :: grad_x, grad_y, grad_z
492 real(kind=rp), dimension(this%Xh_GL%lxyz) :: total_div_gl
493 integer :: e, i, idx, nel, n_GL
494 real(kind=rp), dimension(this%Xh_GLL%lxyz) :: temp_x, temp_y, temp_z
495
496 nel = coef%msh%nelv
497 n_gl = nel * this%Xh_GL%lxyz
498
499 associate(c_gl => this%coef_GL)
500 if (neko_bcknd_device .eq. 1) then
501
502 ! Map mesh velocity (wm) to the GL space
503 call this%GLL_to_GL%map(this%vr, wm_x%x, nel, this%Xh_GL)
504 call this%GLL_to_GL%map(this%vs, wm_y%x, nel, this%Xh_GL)
505 call this%GLL_to_GL%map(this%vt, wm_z%x, nel, this%Xh_GL)
506
507 ! --------------------- X-Momentum
508 ! vx * wm_x
509 call this%GLL_to_GL%map(this%temp, vx%x, nel, this%Xh_GL)
510 call device_col3(this%temp_d, this%temp_d, this%vr_d, n_gl)
511 call opgrad(this%tz, this%tx, this%ty, this%temp, c_gl)
512
513 ! vx * wm_y
514 call this%GLL_to_GL%map(this%temp, vx%x, nel, this%Xh_GL)
515 call device_col3(this%temp_d, this%temp_d, this%vs_d, n_gl)
516 call opgrad(this%tx, this%tbf, this%ty, this%temp, c_gl)
517 call device_add2(this%tz_d, this%tbf_d, n_gl)
518
519 ! vx * wm_z
520 call this%GLL_to_GL%map(this%temp, vx%x, nel, this%Xh_GL)
521 call device_col3(this%temp_d, this%temp_d, this%vt_d, n_gl)
522 call opgrad(this%tx, this%ty, this%tbf, this%temp, c_gl)
523 call device_add2(this%tz_d, this%tbf_d, n_gl)
524
525 ! Map divergence back to GLL space and add to RHS
526 call this%GLL_to_GL%map(this%temp, this%tz, nel, this%Xh_GLL)
527 call device_add2(fx%x_d, this%temp_d, n)
528
529 ! --------------------- Y-Momentum
530 ! vy * wm_x
531 call this%GLL_to_GL%map(this%temp, vy%x, nel, this%Xh_GL)
532 call device_col3(this%temp_d, this%temp_d, this%vr_d, n_gl)
533 call opgrad(this%tz, this%tx, this%ty, this%temp, c_gl)
534
535 ! vy * wm_y
536 call this%GLL_to_GL%map(this%temp, vy%x, nel, this%Xh_GL)
537 call device_col3(this%temp_d, this%temp_d, this%vs_d, n_gl)
538 call opgrad(this%tx, this%tbf, this%ty, this%temp, c_gl)
539 call device_add2(this%tz_d, this%tbf_d, n_gl)
540
541 ! vy * wm_z
542 call this%GLL_to_GL%map(this%temp, vy%x, nel, this%Xh_GL)
543 call device_col3(this%temp_d, this%temp_d, this%vt_d, n_gl)
544 call opgrad(this%tx, this%ty, this%tbf, this%temp, c_gl)
545 call device_add2(this%tz_d, this%tbf_d, n_gl)
546
547 call this%GLL_to_GL%map(this%temp, this%tz, nel, this%Xh_GLL)
548 call device_add2(fy%x_d, this%temp_d, n)
549
550 ! --------------------- Z-Momentum
551 ! wz * wm_x
552 call this%GLL_to_GL%map(this%temp, vz%x, nel, this%Xh_GL)
553 call device_col3(this%temp_d, this%temp_d, this%vr_d, n_gl)
554 call opgrad(this%tz, this%tx, this%ty, this%temp, c_gl)
555
556 ! wz * wm_y
557 call this%GLL_to_GL%map(this%temp, vz%x, nel, this%Xh_GL)
558 call device_col3(this%temp_d, this%temp_d, this%vs_d, n_gl)
559 call opgrad(this%tx, this%tbf, this%ty, this%temp, c_gl)
560 call device_add2(this%tz_d, this%tbf_d, n_gl)
561
562 ! wz * wm_z
563 call this%GLL_to_GL%map(this%temp, vz%x, nel, this%Xh_GL)
564 call device_col3(this%temp_d, this%temp_d, this%vt_d, n_gl)
565 call opgrad(this%tx, this%ty, this%tbf, this%temp, c_gl)
566 call device_add2(this%tz_d, this%tbf_d, n_gl)
567
568 call this%GLL_to_GL%map(this%temp, this%tz, nel, this%Xh_GLL)
569 call device_add2(fz%x_d, this%temp_d, n)
570
571 else
572 !$omp parallel do private(e, i, flux_GL, total_div_GL, idx)
573 do e = 1, coef%msh%nelv
574 ! Map advecting velocity and mesh velocity onto the
575 ! higher-order space
576 call this%GLL_to_GL%map(vx_gl, vx%x(1,1,1,e), 1, this%Xh_GL)
577 call this%GLL_to_GL%map(vy_gl, vy%x(1,1,1,e), 1, this%Xh_GL)
578 call this%GLL_to_GL%map(vz_gl, vz%x(1,1,1,e), 1, this%Xh_GL)
579 call this%GLL_to_GL%map(wm_x_gl, wm_x%x(1,1,1,e), 1, this%Xh_GL)
580 call this%GLL_to_GL%map(wm_y_gl, wm_y%x(1,1,1,e), 1, this%Xh_GL)
581 call this%GLL_to_GL%map(wm_z_gl, wm_z%x(1,1,1,e), 1, this%Xh_GL)
582
583 ! I think below can be written more efficiently. Will fix it later.
584 ! This works for now.
585 ! --------------------- X-Momentum
586 total_div_gl = 0.0_rp
587 ! div(u * wm_*) = d/dx (u * wm_x) + d/dy (u * wm_y) +
588 ! d/dz (u * wm_z)
589
590 flux_gl = vx_gl * wm_x_gl
591 call opgrad(grad_x, grad_y, grad_z, flux_gl, c_gl, e, e)
592 total_div_gl = total_div_gl + grad_x
593 flux_gl = vx_gl * wm_y_gl
594 call opgrad(grad_x, grad_y, grad_z, flux_gl, c_gl, e, e)
595 total_div_gl = total_div_gl + grad_y
596 flux_gl = vx_gl * wm_z_gl
597 call opgrad(grad_x, grad_y, grad_z, flux_gl, c_gl, e, e)
598 total_div_gl = total_div_gl + grad_z
599
600 ! Map back the constructed operator to the original space
601 call this%GLL_to_GL%map(temp_x, total_div_gl, 1, this%Xh_GLL)
602
603 ! --------------------- Y-Momentum
604 total_div_gl = 0.0_rp
605 ! div(v * wm_*) = d/dx (v * wm_x) + d/dy (v * wm_y) +
606 ! d/dz (v * wm_z)
607
608 flux_gl = vy_gl * wm_x_gl
609 call opgrad(grad_x, grad_y, grad_z, flux_gl, c_gl, e, e)
610 total_div_gl = total_div_gl + grad_x
611 flux_gl = vy_gl * wm_y_gl
612 call opgrad(grad_x, grad_y, grad_z, flux_gl, c_gl, e, e)
613 total_div_gl = total_div_gl + grad_y
614 flux_gl = vy_gl * wm_z_gl
615 call opgrad(grad_x, grad_y, grad_z, flux_gl, c_gl, e, e)
616 total_div_gl = total_div_gl + grad_z
617
618 ! Map back the constructed operator to the original space
619 call this%GLL_to_GL%map(temp_y, total_div_gl, 1, this%Xh_GLL)
620
621 ! --------------------- Z-Momentum
622 total_div_gl = 0.0_rp
623 ! div(w * wm_*) = d/dx (w * wm_x) + d/dy (w * wm_y) +
624 ! d/dz (w * wm_z)
625
626 flux_gl = vz_gl * wm_x_gl
627 call opgrad(grad_x, grad_y, grad_z, flux_gl, c_gl, e, e)
628 total_div_gl = total_div_gl + grad_x
629 flux_gl = vz_gl * wm_y_gl
630 call opgrad(grad_x, grad_y, grad_z, flux_gl, c_gl, e, e)
631 total_div_gl = total_div_gl + grad_y
632 flux_gl = vz_gl * wm_z_gl
633 call opgrad(grad_x, grad_y, grad_z, flux_gl, c_gl, e, e)
634 total_div_gl = total_div_gl + grad_z
635
636 ! Map back the constructed operator to the original space
637 call this%GLL_to_GL%map(temp_z, total_div_gl, 1, this%Xh_GLL)
638
639 ! Note we add (+) here since the ALE advection term is
640 ! - div(u * wm) on the LHS. So on the RHS it will be + div(u * wm)
641
642 idx = (e-1)*this%Xh_GLL%lxyz+1
643 do concurrent(i = 0:this%Xh_GLL%lxyz-1)
644 fx%x(i+idx,1,1,1) = fx%x(i+idx,1,1,1) + temp_x(i+1)
645 fy%x(i+idx,1,1,1) = fy%x(i+idx,1,1,1) + temp_y(i+1)
646 fz%x(i+idx,1,1,1) = fz%x(i+idx,1,1,1) + temp_z(i+1)
647 end do
648 end do
649 !$omp end parallel do
650 end if
651 end associate
652 end subroutine compute_ale_advection_dealias
653
654 subroutine recompute_metrics_dealias(this, coef, moving_boundary)
655 class(adv_dealias_t), intent(inout) :: this
656 type(coef_t), intent(in) :: coef
657 logical, intent(in) :: moving_boundary
658 integer :: nel
659
660 if (.not. moving_boundary) return
661
662 nel = coef%msh%nelv
663 call this%GLL_to_GL%map(this%coef_GL%drdx, coef%drdx, nel, this%Xh_GL)
664 call this%GLL_to_GL%map(this%coef_GL%dsdx, coef%dsdx, nel, this%Xh_GL)
665 call this%GLL_to_GL%map(this%coef_GL%dtdx, coef%dtdx, nel, this%Xh_GL)
666
667 call this%GLL_to_GL%map(this%coef_GL%drdy, coef%drdy, nel, this%Xh_GL)
668 call this%GLL_to_GL%map(this%coef_GL%dsdy, coef%dsdy, nel, this%Xh_GL)
669 call this%GLL_to_GL%map(this%coef_GL%dtdy, coef%dtdy, nel, this%Xh_GL)
670
671 call this%GLL_to_GL%map(this%coef_GL%drdz, coef%drdz, nel, this%Xh_GL)
672 call this%GLL_to_GL%map(this%coef_GL%dsdz, coef%dsdz, nel, this%Xh_GL)
673 call this%GLL_to_GL%map(this%coef_GL%dtdz, coef%dtdz, nel, this%Xh_GL)
674
675 end subroutine recompute_metrics_dealias
676
677end module adv_dealias
device::device_get_ptr
Return the device pointer for an associated Fortran array.
Definition device.F90:108
device::device_map
Map a Fortran array to a device (allocate and associate)
Definition device.F90:78
device::device_unmap
Unmap a Fortran array from a device (deassociate and free)
Definition device.F90:84
utils::neko_error
Definition utils.f90:45
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:108
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:358
adv_dealias::free_dealias
subroutine free_dealias(this)
Destructor.
Definition adv_dealias.f90:160
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:236
adv_dealias::compute_ale_advection_dealias
subroutine compute_ale_advection_dealias(this, vx, vy, vz, wm_x, wm_y, wm_z, fx, fy, fz, xh, coef, n, dt)
Definition adv_dealias.f90:480
adv_dealias::recompute_metrics_dealias
subroutine recompute_metrics_dealias(this, coef, moving_boundary)
Definition adv_dealias.f90:655
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_add2
subroutine, public device_add2(a_d, b_d, n, strm)
Vector addition .
Definition device_math.F90:647
device_math::device_vdot3
subroutine, public device_vdot3(dot_d, u1_d, u2_d, u3_d, v1_d, v2_d, v3_d, n, strm)
Compute a dot product (3-d version) assuming vector components etc.
Definition device_math.F90:1212
device_math::device_sub2
subroutine, public device_sub2(a_d, b_d, n, strm)
Vector substraction .
Definition device_math.F90:1070
device_math::device_col3
subroutine, public device_col3(a_d, b_d, c_d, n, strm)
Vector multiplication with 3 vectors .
Definition device_math.F90:1014
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:852
math::sub2
subroutine, public sub2(a, b, n)
Vector substraction .
Definition math.f90:948
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:45
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_metal
integer, parameter neko_bcknd_metal
Definition neko_config.f90:44
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:308
space
Defines a function space.
Definition space.f90:34
space::gl
integer, parameter, public gl
Definition space.f90:49
utils
Utilities.
Definition utils.f90:35
adv_dealias::adv_dealias_t
Type encapsulating advection routines with dealiasing.
Definition adv_dealias.f90:53
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:62
field::field_t
Definition field.f90:48
interpolation::interpolator_t
Interpolation between two space::space_t.
Definition interpolation.f90:53
space::space_t
The function space for the SEM solution fields.
Definition space.f90:63