Neko 1.99.7
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, vz_GL, s_GL), &
427 !$omp& private (f_GL, temp, dsdx, dsdy, dsdz)
428 do e = 1, coef%msh%nelv
429 ! Map advecting velocity onto the higher-order space
430 call this%GLL_to_GL%map(vx_gl, vx%x(1,1,1,e), 1, this%Xh_GL)
431 call this%GLL_to_GL%map(vy_gl, vy%x(1,1,1,e), 1, this%Xh_GL)
432 call this%GLL_to_GL%map(vz_gl, vz%x(1,1,1,e), 1, this%Xh_GL)
433
434 ! Map scalar onto the higher-order space
435 call this%GLL_to_GL%map(s_gl, s%x(1,1,1,e), 1, this%Xh_GL)
436
437 ! Gradient of s in the higher-order space
438 call opgrad(dsdx, dsdy, dsdz, s_gl, c_gl, e, e)
439
440 ! vx * ds/dx + vy * ds/dy + vz * ds/dz for each point in the element
441 do i = 1, this%Xh_GL%lxyz
442 f_gl(i) = vx_gl(i)*dsdx(i) + vy_gl(i)*dsdy(i) + vz_gl(i)*dsdz(i)
443 end do
444
445 ! Map back the contructed operator to the original space
446 call this%GLL_to_GL%map(temp, f_gl, 1, this%Xh_GLL)
447
448 idx = (e-1)*this%Xh_GLL%lxyz + 1
449
450 call sub2(fs%x(idx, 1, 1, 1), temp, this%Xh_GLL%lxyz)
451 end do
452 !$omp end parallel do
453 end if
454 end associate
455
456 end subroutine compute_scalar_advection_dealias
457
458
459 !!> Add the advection term in ALE framework using dealiasing.
460 !! @param this The object.
461 !! @param vx The x component of velocity.
462 !! @param vy The y component of velocity.
463 !! @param vz The z component of velocity.
464 !! @param wm_x The x component of mesh velocity.
465 !! @param wm_y The y component of mesh velocity.
466 !! @param wm_z The z component of mesh velocity.
467 !! @param fx The x component of source term.
468 !! @param fy The y component of source term.
469 !! @param fz The z component of source term.
470 !! @param Xh The function space.
471 !! @param coef The coefficients of the (Xh, mesh) pair.
472 !! @param n Typically the size of the mesh.
473 !! @param dt Current time-step, not required for this method.
474 !! Here, we compute: - div ( u_i * wm ).
475 !! Based on Ho, L.W. A Legendre spectral element method for
476 !! simulation of incompressible unsteady viscous free-surface flows.
477 !! Ph.D. thesis, Massachusetts Institute of Technology, 1989.
478 !! Note: In Nek5000, dealiasing is not done for this term.
479 subroutine compute_ale_advection_dealias(this, vx, vy, vz, &
480 wm_x, wm_y, wm_z, fx, fy, fz, Xh, coef, n, dt)
481 class(adv_dealias_t), intent(inout) :: this
482 type(field_t), intent(inout) :: vx, vy, vz
483 type(field_t), intent(inout) :: wm_x, wm_y, wm_z
484 type(field_t), intent(inout) :: fx, fy, fz
485 type(space_t), intent(in) :: Xh
486 type(coef_t), intent(in) :: coef
487 integer, intent(in) :: n
488 real(kind=rp), intent(in), optional :: dt
489 real(kind=rp), dimension(this%Xh_GL%lxyz) :: vx_gl, vy_gl, vz_gl
490 real(kind=rp), dimension(this%Xh_GL%lxyz) :: wm_x_gl, wm_y_gl, wm_z_gl
491 real(kind=rp), dimension(this%Xh_GL%lxyz) :: flux_gl
492 real(kind=rp), dimension(this%Xh_GL%lxyz) :: grad_x, grad_y, grad_z
493 real(kind=rp), dimension(this%Xh_GL%lxyz) :: total_div_gl
494 integer :: e, i, idx, nel, n_GL
495 real(kind=rp), dimension(this%Xh_GLL%lxyz) :: temp_x, temp_y, temp_z
496
497 nel = coef%msh%nelv
498 n_gl = nel * this%Xh_GL%lxyz
499
500 associate(c_gl => this%coef_GL)
501 if (neko_bcknd_device .eq. 1) then
502
503 ! Map mesh velocity (wm) to the GL space
504 call this%GLL_to_GL%map(this%vr, wm_x%x, nel, this%Xh_GL)
505 call this%GLL_to_GL%map(this%vs, wm_y%x, nel, this%Xh_GL)
506 call this%GLL_to_GL%map(this%vt, wm_z%x, nel, this%Xh_GL)
507
508 ! --------------------- X-Momentum
509 ! vx * wm_x
510 call this%GLL_to_GL%map(this%temp, vx%x, nel, this%Xh_GL)
511 call device_col3(this%temp_d, this%temp_d, this%vr_d, n_gl)
512 call opgrad(this%tz, this%tx, this%ty, this%temp, c_gl)
513
514 ! vx * wm_y
515 call this%GLL_to_GL%map(this%temp, vx%x, nel, this%Xh_GL)
516 call device_col3(this%temp_d, this%temp_d, this%vs_d, n_gl)
517 call opgrad(this%tx, this%tbf, this%ty, this%temp, c_gl)
518 call device_add2(this%tz_d, this%tbf_d, n_gl)
519
520 ! vx * wm_z
521 call this%GLL_to_GL%map(this%temp, vx%x, nel, this%Xh_GL)
522 call device_col3(this%temp_d, this%temp_d, this%vt_d, n_gl)
523 call opgrad(this%tx, this%ty, this%tbf, this%temp, c_gl)
524 call device_add2(this%tz_d, this%tbf_d, n_gl)
525
526 ! Map divergence back to GLL space and add to RHS
527 call this%GLL_to_GL%map(this%temp, this%tz, nel, this%Xh_GLL)
528 call device_add2(fx%x_d, this%temp_d, n)
529
530 ! --------------------- Y-Momentum
531 ! vy * wm_x
532 call this%GLL_to_GL%map(this%temp, vy%x, nel, this%Xh_GL)
533 call device_col3(this%temp_d, this%temp_d, this%vr_d, n_gl)
534 call opgrad(this%tz, this%tx, this%ty, this%temp, c_gl)
535
536 ! vy * wm_y
537 call this%GLL_to_GL%map(this%temp, vy%x, nel, this%Xh_GL)
538 call device_col3(this%temp_d, this%temp_d, this%vs_d, n_gl)
539 call opgrad(this%tx, this%tbf, this%ty, this%temp, c_gl)
540 call device_add2(this%tz_d, this%tbf_d, n_gl)
541
542 ! vy * wm_z
543 call this%GLL_to_GL%map(this%temp, vy%x, nel, this%Xh_GL)
544 call device_col3(this%temp_d, this%temp_d, this%vt_d, n_gl)
545 call opgrad(this%tx, this%ty, this%tbf, this%temp, c_gl)
546 call device_add2(this%tz_d, this%tbf_d, n_gl)
547
548 call this%GLL_to_GL%map(this%temp, this%tz, nel, this%Xh_GLL)
549 call device_add2(fy%x_d, this%temp_d, n)
550
551 ! --------------------- Z-Momentum
552 ! wz * wm_x
553 call this%GLL_to_GL%map(this%temp, vz%x, nel, this%Xh_GL)
554 call device_col3(this%temp_d, this%temp_d, this%vr_d, n_gl)
555 call opgrad(this%tz, this%tx, this%ty, this%temp, c_gl)
556
557 ! wz * wm_y
558 call this%GLL_to_GL%map(this%temp, vz%x, nel, this%Xh_GL)
559 call device_col3(this%temp_d, this%temp_d, this%vs_d, n_gl)
560 call opgrad(this%tx, this%tbf, this%ty, this%temp, c_gl)
561 call device_add2(this%tz_d, this%tbf_d, n_gl)
562
563 ! wz * wm_z
564 call this%GLL_to_GL%map(this%temp, vz%x, nel, this%Xh_GL)
565 call device_col3(this%temp_d, this%temp_d, this%vt_d, n_gl)
566 call opgrad(this%tx, this%ty, this%tbf, this%temp, c_gl)
567 call device_add2(this%tz_d, this%tbf_d, n_gl)
568
569 call this%GLL_to_GL%map(this%temp, this%tz, nel, this%Xh_GLL)
570 call device_add2(fz%x_d, this%temp_d, n)
571
572 else
573 !$omp parallel do private (e, i, idx, vx_GL, vy_GL, vz_GL), &
574 !$omp& private (wm_x_GL, wm_y_GL, wm_z_GL, flux_GL, total_div_GL), &
575 !$omp& private (grad_x, grad_y, grad_z, temp_x, temp_y, temp_z)
576 do e = 1, coef%msh%nelv
577 ! Map advecting velocity and mesh velocity onto the
578 ! higher-order space
579 call this%GLL_to_GL%map(vx_gl, vx%x(1,1,1,e), 1, this%Xh_GL)
580 call this%GLL_to_GL%map(vy_gl, vy%x(1,1,1,e), 1, this%Xh_GL)
581 call this%GLL_to_GL%map(vz_gl, vz%x(1,1,1,e), 1, this%Xh_GL)
582 call this%GLL_to_GL%map(wm_x_gl, wm_x%x(1,1,1,e), 1, this%Xh_GL)
583 call this%GLL_to_GL%map(wm_y_gl, wm_y%x(1,1,1,e), 1, this%Xh_GL)
584 call this%GLL_to_GL%map(wm_z_gl, wm_z%x(1,1,1,e), 1, this%Xh_GL)
585
586 ! I think below can be written more efficiently. Will fix it later.
587 ! This works for now.
588 ! --------------------- X-Momentum
589 total_div_gl = 0.0_rp
590 ! div(u * wm_*) = d/dx (u * wm_x) + d/dy (u * wm_y) +
591 ! d/dz (u * wm_z)
592
593 flux_gl = vx_gl * wm_x_gl
594 call opgrad(grad_x, grad_y, grad_z, flux_gl, c_gl, e, e)
595 total_div_gl = total_div_gl + grad_x
596 flux_gl = vx_gl * wm_y_gl
597 call opgrad(grad_x, grad_y, grad_z, flux_gl, c_gl, e, e)
598 total_div_gl = total_div_gl + grad_y
599 flux_gl = vx_gl * wm_z_gl
600 call opgrad(grad_x, grad_y, grad_z, flux_gl, c_gl, e, e)
601 total_div_gl = total_div_gl + grad_z
602
603 ! Map back the constructed operator to the original space
604 call this%GLL_to_GL%map(temp_x, total_div_gl, 1, this%Xh_GLL)
605
606 ! --------------------- Y-Momentum
607 total_div_gl = 0.0_rp
608 ! div(v * wm_*) = d/dx (v * wm_x) + d/dy (v * wm_y) +
609 ! d/dz (v * wm_z)
610
611 flux_gl = vy_gl * wm_x_gl
612 call opgrad(grad_x, grad_y, grad_z, flux_gl, c_gl, e, e)
613 total_div_gl = total_div_gl + grad_x
614 flux_gl = vy_gl * wm_y_gl
615 call opgrad(grad_x, grad_y, grad_z, flux_gl, c_gl, e, e)
616 total_div_gl = total_div_gl + grad_y
617 flux_gl = vy_gl * wm_z_gl
618 call opgrad(grad_x, grad_y, grad_z, flux_gl, c_gl, e, e)
619 total_div_gl = total_div_gl + grad_z
620
621 ! Map back the constructed operator to the original space
622 call this%GLL_to_GL%map(temp_y, total_div_gl, 1, this%Xh_GLL)
623
624 ! --------------------- Z-Momentum
625 total_div_gl = 0.0_rp
626 ! div(w * wm_*) = d/dx (w * wm_x) + d/dy (w * wm_y) +
627 ! d/dz (w * wm_z)
628
629 flux_gl = vz_gl * wm_x_gl
630 call opgrad(grad_x, grad_y, grad_z, flux_gl, c_gl, e, e)
631 total_div_gl = total_div_gl + grad_x
632 flux_gl = vz_gl * wm_y_gl
633 call opgrad(grad_x, grad_y, grad_z, flux_gl, c_gl, e, e)
634 total_div_gl = total_div_gl + grad_y
635 flux_gl = vz_gl * wm_z_gl
636 call opgrad(grad_x, grad_y, grad_z, flux_gl, c_gl, e, e)
637 total_div_gl = total_div_gl + grad_z
638
639 ! Map back the constructed operator to the original space
640 call this%GLL_to_GL%map(temp_z, total_div_gl, 1, this%Xh_GLL)
641
642 ! Note we add (+) here since the ALE advection term is
643 ! - div(u * wm) on the LHS. So on the RHS it will be + div(u * wm)
644
645 idx = (e-1)*this%Xh_GLL%lxyz+1
646 do concurrent(i = 0:this%Xh_GLL%lxyz-1)
647 fx%x(i+idx,1,1,1) = fx%x(i+idx,1,1,1) + temp_x(i+1)
648 fy%x(i+idx,1,1,1) = fy%x(i+idx,1,1,1) + temp_y(i+1)
649 fz%x(i+idx,1,1,1) = fz%x(i+idx,1,1,1) + temp_z(i+1)
650 end do
651 end do
652 !$omp end parallel do
653 end if
654 end associate
655 end subroutine compute_ale_advection_dealias
656
657 subroutine recompute_metrics_dealias(this, coef, moving_boundary)
658 class(adv_dealias_t), intent(inout) :: this
659 type(coef_t), intent(in) :: coef
660 logical, intent(in) :: moving_boundary
661 integer :: nel
662
663 if (.not. moving_boundary) return
664
665 nel = coef%msh%nelv
666 call this%GLL_to_GL%map(this%coef_GL%drdx, coef%drdx, nel, this%Xh_GL)
667 call this%GLL_to_GL%map(this%coef_GL%dsdx, coef%dsdx, nel, this%Xh_GL)
668 call this%GLL_to_GL%map(this%coef_GL%dtdx, coef%dtdx, nel, this%Xh_GL)
669
670 call this%GLL_to_GL%map(this%coef_GL%drdy, coef%drdy, nel, this%Xh_GL)
671 call this%GLL_to_GL%map(this%coef_GL%dsdy, coef%dsdy, nel, this%Xh_GL)
672 call this%GLL_to_GL%map(this%coef_GL%dtdy, coef%dtdy, nel, this%Xh_GL)
673
674 call this%GLL_to_GL%map(this%coef_GL%drdz, coef%drdz, nel, this%Xh_GL)
675 call this%GLL_to_GL%map(this%coef_GL%dsdz, coef%dsdz, nel, this%Xh_GL)
676 call this%GLL_to_GL%map(this%coef_GL%dtdz, coef%dtdz, nel, this%Xh_GL)
677
678 end subroutine recompute_metrics_dealias
679
680end module adv_dealias
device::device_get_ptr
Return the device pointer for an associated Fortran array.
Definition device.F90:113
device::device_map
Map a Fortran array to a device (allocate and associate)
Definition device.F90:83
device::device_unmap
Unmap a Fortran array from a device (deassociate and free)
Definition device.F90:89
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:481
adv_dealias::recompute_metrics_dealias
subroutine recompute_metrics_dealias(this, coef, moving_boundary)
Definition adv_dealias.f90:658
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:703
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:1267
device_math::device_sub2
subroutine, public device_sub2(a_d, b_d, n, strm)
Vector substraction .
Definition device_math.F90:1125
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:1069
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:50
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:63
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:64