Neko 1.99.2
A portable framework for high-order spectral element flow simulations
Loading...
Searching...
No Matches
adv_oifs.f90
Go to the documentation of this file.
1! Copyright (c) 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_oifs
35 use advection, only: advection_t
36 use num_types, only: rp
37 use space, only: space_t, gl
38 use field, only: field_t
39 use coefs, only: coef_t
40 use math, only: copy, rzero
41 use operators, only: runge_kutta, set_convect_rst
42 use neko_config, only: neko_bcknd_device, neko_bcknd_sx, neko_bcknd_xsmm
43 use interpolation, only: interpolator_t
44 use time_interpolator, only: time_interpolator_t
45 use field_series, only: field_series_t
46 use field_list, only: field_list_t
47 use time_scheme_controller, only: time_scheme_controller_t
48 use device, only: device_map, device_free
49 use device_math, only: device_addcol3s2, device_rzero
50 use, intrinsic :: iso_c_binding, only: c_ptr, c_null_ptr, c_associated
51 implicit none
52 private
53
54 !! Type encapsulating operator-integration-factor splitting advection
55 !! routines with dealiasing applied
56 !! Literature:
57 !! https://www.mcs.anl.gov/~fischer/nek5000/oifs.pdf
58 !! https://publications.anl.gov/anlpubs/2017/12/140626.pdf
59 !! https://dl.acm.org/doi/abs/10.1007/BF01063118
60 type, public, extends(advection_t) :: adv_oifs_t
62 integer :: ntaubd
64 type(coef_t) :: coef_gl
66 type(coef_t), pointer :: coef_gll => null()
68 type(interpolator_t) :: gll_to_gl
70 type(space_t) :: xh_gl
72 type(space_t), pointer :: xh_gll => null()
74 type(time_interpolator_t) :: dtime
76 type(field_series_t), pointer :: ulag, vlag, wlag, slag => null()
78 real(kind=rp), pointer :: ctlag(:) => null()
80 real(kind=rp), pointer :: dctlag(:) => null()
82 type(time_scheme_controller_t), pointer :: oifs_scheme => null()
84 type(field_t) :: cr_gl, cs_gl, ct_gl
86 type(field_series_t) :: convr_gl, convs_gl, convt_gl
88 type(field_t), pointer :: cr_k1, cs_k1, ct_k1
89 type(field_t), pointer :: cr_k23, cs_k23, ct_k23
90 type(field_t), pointer :: cr_k4, cs_k4, ct_k4
92 type(field_list_t) :: conv_k1, conv_k23, conv_k4
94 real(kind=rp), allocatable :: cx(:), cy(:), cz(:)
96 type(c_ptr) :: cx_d = c_null_ptr, cy_d = c_null_ptr, cz_d = c_null_ptr
97
98 contains
101 procedure, pass(this) :: compute => adv_oifs_compute
104 procedure, pass(this) :: compute_scalar => adv_oifs_compute_scalar
106 procedure, pass(this) :: init => adv_oifs_init
108 procedure, pass(this) :: set_conv_velocity_fst
110 procedure, pass(this) :: free => adv_oifs_free
111 end type adv_oifs_t
112
113contains
114
126 subroutine adv_oifs_init(this, lxd, coef, ctarget, ulag, vlag, wlag, &
127 dtlag, tlag, time_scheme, slag)
128 implicit none
129 class(adv_oifs_t) :: this
130 integer, intent(in) :: lxd
131 type(coef_t), target :: coef
132 real(kind=rp), intent(in) :: ctarget
133 type(field_series_t), target, intent(in) :: ulag, vlag, wlag
134 real(kind=rp), target, intent(in) :: dtlag(10)
135 real(kind=rp), target, intent(in) :: tlag(10)
136 type(time_scheme_controller_t), target, intent(in) :: time_scheme
137 type(field_series_t), target, optional :: slag
138 integer :: nel, n_GL, n, idx, idy, idz
139 real(kind=rp) :: max_cfl_rk4
140
141 ! stability limit for RK4 including safety factor
142 max_cfl_rk4 = 2.0
143 this%ntaubd = max(int(ctarget/max_cfl_rk4),1)
144
145 call this%Xh_GL%init(gl, lxd, lxd, lxd)
146 this%Xh_GLL => coef%Xh
147 this%coef_GLL => coef
148 call this%GLL_to_GL%init(this%Xh_GL, this%Xh_GLL)
149
150 call this%coef_GL%init(this%Xh_GL, coef%msh)
151
152 call this%cr_GL%init(coef%msh, this%Xh_GL)
153 call this%cs_GL%init(coef%msh, this%Xh_GL)
154 call this%ct_GL%init(coef%msh, this%Xh_GL)
155
156 nel = coef%msh%nelv
157 n_gl = nel*this%Xh_GL%lxyz
158 n = nel*coef%Xh%lxyz
159
160 call this%GLL_to_GL%map(this%coef_GL%drdx, coef%drdx, nel, this%Xh_GL)
161 call this%GLL_to_GL%map(this%coef_GL%dsdx, coef%dsdx, nel, this%Xh_GL)
162 call this%GLL_to_GL%map(this%coef_GL%dtdx, coef%dtdx, nel, this%Xh_GL)
163 call this%GLL_to_GL%map(this%coef_GL%drdy, coef%drdy, nel, this%Xh_GL)
164 call this%GLL_to_GL%map(this%coef_GL%dsdy, coef%dsdy, nel, this%Xh_GL)
165 call this%GLL_to_GL%map(this%coef_GL%dtdy, coef%dtdy, nel, this%Xh_GL)
166 call this%GLL_to_GL%map(this%coef_GL%drdz, coef%drdz, nel, this%Xh_GL)
167 call this%GLL_to_GL%map(this%coef_GL%dsdz, coef%dsdz, nel, this%Xh_GL)
168 call this%GLL_to_GL%map(this%coef_GL%dtdz, coef%dtdz, nel, this%Xh_GL)
169
170
171 allocate(this%cx(n_gl))
172 allocate(this%cy(n_gl))
173 allocate(this%cz(n_gl))
174
175 allocate(this%cr_k1)
176 allocate(this%cs_k1)
177 allocate(this%ct_k1)
178 allocate(this%cr_k23)
179 allocate(this%cs_k23)
180 allocate(this%ct_k23)
181 allocate(this%cr_k4)
182 allocate(this%cs_k4)
183 allocate(this%ct_k4)
184
185
186 call this%cr_k1%init(coef%msh, this%Xh_GL)
187 call this%cs_k1%init(coef%msh, this%Xh_GL)
188 call this%ct_k1%init(coef%msh, this%Xh_GL)
189
190 call this%cr_k23%init(coef%msh, this%Xh_GL)
191 call this%cs_k23%init(coef%msh, this%Xh_GL)
192 call this%ct_k23%init(coef%msh, this%Xh_GL)
193
194 call this%cr_k4%init(coef%msh, this%Xh_GL)
195 call this%cs_k4%init(coef%msh, this%Xh_GL)
196 call this%ct_k4%init(coef%msh, this%Xh_GL)
197
198 call this%conv_k1%init(3)
199 call this%conv_k23%init(3)
200 call this%conv_k4%init(3)
201
202 call this%conv_k1%assign(1, this%cr_k1)
203 call this%conv_k1%assign(2, this%cs_k1)
204 call this%conv_k1%assign(3, this%ct_k1)
205
206 call this%conv_k23%assign(1, this%cr_k23)
207 call this%conv_k23%assign(2, this%cs_k23)
208 call this%conv_k23%assign(3, this%ct_k23)
209
210 call this%conv_k4%assign(1, this%cr_k4)
211 call this%conv_k4%assign(2, this%cs_k4)
212 call this%conv_k4%assign(3, this%ct_k4)
213
214 call this%dtime%init(1)
215 this%ulag => ulag
216 this%vlag => vlag
217 this%wlag => wlag
218 this%ctlag => tlag
219 this%dctlag => dtlag
220 this%oifs_scheme => time_scheme
221
222 if (neko_bcknd_device .eq. 1) then
223 call device_map(this%cx, this%cx_d, n_gl)
224 call device_map(this%cy, this%cy_d, n_gl)
225 call device_map(this%cz, this%cz_d, n_gl)
226 end if
227
228 ! Initializing the convecting fields
229 ! Map the velocity fields from GLL space to GL space
230 call this%GLL_to_GL%map(this%cx, this%ulag%f%x, nel, this%Xh_GL)
231 call this%GLL_to_GL%map(this%cy, this%vlag%f%x, nel, this%Xh_GL)
232 call this%GLL_to_GL%map(this%cz, this%wlag%f%x, nel, this%Xh_GL)
233
234 ! Set the convecting field in the rst format
235 call set_convect_rst(this%cr_GL, this%cs_GL, this%ct_GL, &
236 this%cx, this%cy, this%cz, this%Xh_GL, this%coef_GL)
237
238 ! Set the convecting field series
239 call this%convr_GL%init(this%cr_GL, 3)
240 call this%convs_GL%init(this%cs_GL, 3)
241 call this%convt_GL%init(this%ct_GL, 3)
242
243 ! Repeat for previous time-steps
244 call this%GLL_to_GL%map(this%cx, this%ulag%lf(1)%x, nel, this%Xh_GL)
245 call this%GLL_to_GL%map(this%cy, this%vlag%lf(1)%x, nel, this%Xh_GL)
246 call this%GLL_to_GL%map(this%cz, this%wlag%lf(1)%x, nel, this%Xh_GL)
247
248 call set_convect_rst(this%cr_GL, this%cs_GL, this%ct_GL, &
249 this%cx, this%cy, this%cz, this%Xh_GL, this%coef_GL)
250
251 this%convr_GL%lf(1) = this%cr_GL
252 this%convs_GL%lf(1) = this%cs_GL
253 this%convt_GL%lf(1) = this%ct_GL
254
255 call this%GLL_to_GL%map(this%cx, this%ulag%lf(2)%x, nel, this%Xh_GL)
256 call this%GLL_to_GL%map(this%cy, this%vlag%lf(2)%x, nel, this%Xh_GL)
257 call this%GLL_to_GL%map(this%cz, this%wlag%lf(2)%x, nel, this%Xh_GL)
258
259 call set_convect_rst(this%cr_GL, this%cs_GL, this%ct_GL, &
260 this%cx, this%cy, this%cz, this%Xh_GL, this%coef_GL)
261
262 this%convr_GL%lf(2) = this%cr_GL
263 this%convs_GL%lf(2) = this%cs_GL
264 this%convt_GL%lf(2) = this%ct_GL
265
266 ! Initilize the lagged scalar field, if present.
267 if (present(slag)) then
268 this%slag => slag
269 end if
270
271 end subroutine adv_oifs_init
272
274 subroutine adv_oifs_free(this)
275 class(adv_oifs_t), intent(inout) :: this
276
277 call this%coef_GL%free()
278
279 nullify(this%coef_GLL)
280
281 call this%GLL_to_GL%free()
282
283 call this%Xh_GL%free()
284
285 nullify(this%Xh_GLL)
286
287 call this%dtime%free()
288
289 call this%cr_GL%free()
290 call this%cs_GL%free()
291 call this%ct_GL%free()
292
293 call this%convr_GL%free()
294 call this%convs_GL%free()
295 call this%convt_GL%free()
296
297 call this%conv_k1%free()
298 call this%conv_k23%free()
299 call this%conv_k4%free()
300
301 if (associated(this%cr_k1)) then
302 deallocate(this%cr_k1)
303 end if
304 if (associated(this%cs_k1)) then
305 deallocate(this%cs_k1)
306 end if
307 if (associated(this%ct_k1)) then
308 deallocate(this%ct_k1)
309 end if
310 if (associated(this%cr_k23)) then
311 deallocate(this%cr_k23)
312 end if
313 if (associated(this%cs_k23)) then
314 deallocate(this%cs_k23)
315 end if
316 if (associated(this%ct_k23)) then
317 deallocate(this%ct_k23)
318 end if
319 if (associated(this%cr_k4)) then
320 deallocate(this%cr_k4)
321 end if
322 if (associated(this%cs_k4)) then
323 deallocate(this%cs_k4)
324 end if
325 if (associated(this%ct_k4)) then
326 deallocate(this%ct_k4)
327 end if
328
329 nullify(this%ulag)
330 nullify(this%vlag)
331 nullify(this%wlag)
332 nullify(this%slag)
333 nullify(this%ctlag)
334 nullify(this%dctlag)
335 nullify(this%oifs_scheme)
336 nullify(this%cr_k1)
337 nullify(this%cs_k1)
338 nullify(this%ct_k1)
339 nullify(this%cr_k23)
340 nullify(this%cs_k23)
341 nullify(this%ct_k23)
342 nullify(this%cr_k4)
343 nullify(this%cs_k4)
344 nullify(this%ct_k4)
345
346 if (allocated(this%cx)) then
347 deallocate(this%cx)
348 end if
349 if (allocated(this%cy)) then
350 deallocate(this%cy)
351 end if
352 if (allocated(this%cz)) then
353 deallocate(this%cz)
354 end if
355 if (c_associated(this%cx_d)) then
356 call device_free(this%cx_d)
357 end if
358 if (c_associated(this%cy_d)) then
359 call device_free(this%cy_d)
360 end if
361 if (c_associated(this%cz_d)) then
362 call device_free(this%cz_d)
363 end if
364
365 end subroutine adv_oifs_free
366
372 subroutine set_conv_velocity_fst(this, u, v, w)
373 implicit none
374 class(adv_oifs_t), intent(inout) :: this
375 type(field_t), intent(inout) :: u, v, w
376 integer :: i, nel, n_GL, idx, idy, idz
377
378 nel = this%coef_GLL%msh%nelv
379 n_gl = nel*this%Xh_GL%lxyz
380
381 call this%convr_GL%update()
382 call this%convs_GL%update()
383 call this%convt_GL%update()
384
385 call this%GLL_to_GL%map(this%cx, u%x, nel, this%Xh_GL)
386 call this%GLL_to_GL%map(this%cy, v%x, nel, this%Xh_GL)
387 call this%GLL_to_GL%map(this%cz, w%x, nel, this%Xh_GL)
388
389 call set_convect_rst(this%cr_GL, this%cs_GL, this%ct_GL, &
390 this%cx, this%cy, this%cz, this%Xh_GL, this%coef_GL)
391
392 this%convr_GL%f = this%cr_GL
393 this%convs_GL%f = this%cs_GL
394 this%convt_GL%f = this%ct_GL
395
396 end subroutine set_conv_velocity_fst
397
398
411 subroutine adv_oifs_compute(this, vx, vy, vz, fx, fy, fz, Xh, coef, n, dt)
412 implicit none
413 class(adv_oifs_t), intent(inout) :: this
414 type(field_t), intent(inout) :: vx, vy, vz
415 type(field_t), intent(inout) :: fx, fy, fz
416 type(space_t), intent(in) :: Xh
417 type(coef_t), intent(in) :: coef
418 integer, intent(in) :: n
419 real(kind=rp), intent(in), optional :: dt
420 real(kind=rp) :: tau, tau1, th, dtau
421 integer :: i, ilag, itau, nel, n_GL
422
423 nel = coef%msh%nelv
424 n_gl = nel * this%Xh_GL%lxyz
425
426 associate(ulag => this%ulag, vlag => this%vlag, wlag => this%wlag, &
427 ctlag => this%ctlag, dctlag => this%dctlag, dtime => this%dtime, &
428 xh_gl => this%Xh_GL, coef_gl => this%coef_GL, ntaubd => this%ntaubd, &
429 gll_to_gl => this%GLL_to_GL, oifs_scheme => this%oifs_scheme, &
430 cr_k1 => this%cr_K1, cs_k1 => this%cs_K1, ct_k1 => this%ct_K1, &
431 cr_k23 => this%cr_K23, cs_k23 => this%cs_K23, ct_k23 => this%ct_K23, &
432 cr_k4 => this%cr_K4, cs_k4 => this%cs_K4, ct_k4 => this%ct_K4, &
433 convr_gl => this%convr_GL, convs_gl => this%convs_GL, &
434 convt_gl => this%convt_GL, conv_k1 => this%conv_k1, &
435 conv_k23 => this%conv_k23, conv_k4 => this%conv_k4)
436
437 call dtime%init(oifs_scheme%ndiff)
438
439 tau = ctlag(oifs_scheme%ndiff)
440
441 call this%set_conv_velocity_fst(vx, vy, vz)
442
443 if (neko_bcknd_device .eq. 1) then
444 call device_rzero(fx%x_d,n)
445 call device_rzero(fy%x_d,n)
446 call device_rzero(fz%x_d,n)
447 else
448 call rzero(fx%x,n)
449 call rzero(fy%x,n)
450 call rzero(fz%x,n)
451 end if
452
453 do ilag = oifs_scheme%ndiff, 1, -1
454 if (neko_bcknd_device .eq. 1) then
455 if (ilag .eq. 1) then
456 call device_addcol3s2(fx%x_d, vx%x_d, coef%B_d, &
457 oifs_scheme%diffusion_coeffs(2), n)
458 call device_addcol3s2(fy%x_d, vy%x_d, coef%B_d, &
459 oifs_scheme%diffusion_coeffs(2), n)
460 call device_addcol3s2(fz%x_d, vz%x_d, coef%B_d, &
461 oifs_scheme%diffusion_coeffs(2), n)
462 else
463 call device_addcol3s2(fx%x_d, ulag%lf(ilag-1)%x_d, coef%B_d, &
464 oifs_scheme%diffusion_coeffs(ilag+1), n)
465 call device_addcol3s2(fy%x_d, vlag%lf(ilag-1)%x_d, coef%B_d, &
466 oifs_scheme%diffusion_coeffs(ilag+1), n)
467 call device_addcol3s2(fz%x_d, wlag%lf(ilag-1)%x_d, coef%B_d, &
468 oifs_scheme%diffusion_coeffs(ilag+1), n)
469 end if
470 else
471 if (ilag .eq. 1) then
472 do i = 1, n
473 fx%x(i,1,1,1) = fx%x(i,1,1,1) + &
474 oifs_scheme%diffusion_coeffs(2) &
475 * vx%x(i,1,1,1) * coef%B(i,1,1,1)
476 fy%x(i,1,1,1) = fy%x(i,1,1,1) + &
477 oifs_scheme%diffusion_coeffs(2) &
478 * vy%x(i,1,1,1) * coef%B(i,1,1,1)
479 fz%x(i,1,1,1) = fz%x(i,1,1,1) + &
480 oifs_scheme%diffusion_coeffs(2) &
481 * vz%x(i,1,1,1) * coef%B(i,1,1,1)
482 end do
483 else
484 do i = 1, n
485 fx%x(i,1,1,1) = fx%x(i,1,1,1) + &
486 oifs_scheme%diffusion_coeffs(ilag+1) &
487 * ulag%lf(ilag-1)%x(i,1,1,1) &
488 * coef%B(i,1,1,1)
489 fy%x(i,1,1,1) = fy%x(i,1,1,1) + &
490 oifs_scheme%diffusion_coeffs(ilag+1) &
491 * vlag%lf(ilag-1)%x(i,1,1,1) &
492 * coef%B(i,1,1,1)
493 fz%x(i,1,1,1) = fz%x(i,1,1,1) + &
494 oifs_scheme%diffusion_coeffs(ilag+1) &
495 * wlag%lf(ilag-1)%x(i,1,1,1) &
496 * coef%B(i,1,1,1)
497 end do
498 end if
499 end if
500 dtau = dctlag(ilag)/real(ntaubd)
501 do itau = 1, ntaubd
502 th = tau + dtau/2.
503 tau1 = tau + dtau
504 call dtime%interpolate_scalar(tau, cr_k1, convr_gl, ctlag, n_gl)
505 call dtime%interpolate_scalar(tau, cs_k1, convs_gl, ctlag, n_gl)
506 call dtime%interpolate_scalar(tau, ct_k1, convt_gl, ctlag, n_gl)
507 call dtime%interpolate_scalar(th, cr_k23, convr_gl, ctlag, n_gl)
508 call dtime%interpolate_scalar(th, cs_k23, convs_gl, ctlag, n_gl)
509 call dtime%interpolate_scalar(th, ct_k23, convt_gl, ctlag, n_gl)
510 call dtime%interpolate_scalar(tau1, cr_k4, convr_gl, ctlag, n_gl)
511 call dtime%interpolate_scalar(tau1, cs_k4, convs_gl, ctlag, n_gl)
512 call dtime%interpolate_scalar(tau1, ct_k4, convt_gl, ctlag, n_gl)
513 call runge_kutta(fx, conv_k1, conv_k23, conv_k4, xh, xh_gl, &
514 coef, coef_gl, gll_to_gl, tau, dtau, &
515 n, nel, n_gl)
516 call runge_kutta(fy, conv_k1, conv_k23, conv_k4, xh, xh_gl, &
517 coef, coef_gl, gll_to_gl, tau, dtau, &
518 n, nel, n_gl)
519 call runge_kutta(fz, conv_k1, conv_k23, conv_k4, xh, xh_gl, &
520 coef, coef_gl, gll_to_gl, tau, dtau, &
521 n, nel, n_gl)
522 tau = tau1
523 end do
524 end do
525
526 end associate
527
528 end subroutine adv_oifs_compute
541 subroutine adv_oifs_compute_scalar(this, vx, vy, vz, s, fs, Xh, coef, n, dt)
542 implicit none
543 class(adv_oifs_t), intent(inout) :: this
544 type(field_t), intent(inout) :: vx, vy, vz
545 type(field_t), intent(inout) :: fs
546 type(field_t), intent(inout) :: s
547 type(space_t), intent(in) :: Xh
548 type(coef_t), intent(in) :: coef
549 integer, intent(in) :: n
550 real(kind=rp), intent(in), optional :: dt
551 real(kind=rp) :: tau, tau1, th, dtau
552 integer :: i, ilag, itau, nel, n_GL
553 nel = coef%msh%nelv
554 n_gl = nel * this%Xh_GL%lxyz
555
556 associate(slag => this%slag, ctlag => this%ctlag, dctlag => this%dctlag, &
557 dtime => this%dtime, xh_gl => this%Xh_GL, coef_gl => this%coef_GL, &
558 ntaubd => this%ntaubd, gll_to_gl => this%GLL_to_GL, &
559 oifs_scheme => this%oifs_scheme, cr_k1 => this%cr_K1, &
560 cs_k1 => this%cs_K1, ct_k1 => this%ct_K1, cr_k23 => this%cr_K23, &
561 cs_k23 => this%cs_K23, ct_k23 => this%ct_K23, cr_k4 => this%cr_K4, &
562 cs_k4 => this%cs_K4, ct_k4 => this%ct_K4, &
563 convr_gl => this%convr_GL, convs_gl => this%convs_GL, &
564 convt_gl => this%convt_GL, conv_k1 => this%conv_k1, &
565 conv_k23 => this%conv_k23, conv_k4 => this%conv_k4)
566
567 call dtime%init(oifs_scheme%ndiff)
568
569 tau = ctlag(oifs_scheme%ndiff)
570
571 call this%set_conv_velocity_fst(vx, vy, vz)
572
573 if (neko_bcknd_device .eq. 1) then
574 call device_rzero(fs%x_d,n)
575 else
576 call rzero(fs%x,n)
577 end if
578
579 do ilag = oifs_scheme%ndiff, 1, -1
580 if (neko_bcknd_device .eq. 1) then
581 if (ilag .eq. 1) then
582 call device_addcol3s2(fs%x_d, s%x_d, coef%B_d, &
583 oifs_scheme%diffusion_coeffs(2), n)
584 else
585 call device_addcol3s2(fs%x_d, slag%lf(ilag-1)%x_d, coef%B_d, &
586 oifs_scheme%diffusion_coeffs(ilag+1), n)
587 end if
588 else
589 if (ilag .eq. 1) then
590 do i = 1, n
591 fs%x(i,1,1,1) = fs%x(i,1,1,1) + &
592 oifs_scheme%diffusion_coeffs(2) &
593 * s%x(i,1,1,1) * coef%B(i,1,1,1)
594 end do
595 else
596 do i = 1, n
597 fs%x(i,1,1,1) = fs%x(i,1,1,1) + &
598 oifs_scheme%diffusion_coeffs(ilag+1) &
599 * slag%lf(ilag-1)%x(i,1,1,1) * coef%B(i,1,1,1)
600 end do
601 end if
602 end if
603 dtau = dctlag(ilag)/real(ntaubd)
604 do itau = 1, ntaubd
605 th = tau + dtau/2.
606 tau1 = tau + dtau
607 call dtime%interpolate_scalar(tau, cr_k1, convr_gl, ctlag, n_gl)
608 call dtime%interpolate_scalar(tau, cs_k1, convs_gl, ctlag, n_gl)
609 call dtime%interpolate_scalar(tau, ct_k1, convt_gl, ctlag, n_gl)
610 call dtime%interpolate_scalar(th, cr_k23, convr_gl, ctlag, n_gl)
611 call dtime%interpolate_scalar(th, cs_k23, convs_gl, ctlag, n_gl)
612 call dtime%interpolate_scalar(th, ct_k23, convt_gl, ctlag, n_gl)
613 call dtime%interpolate_scalar(tau1, cr_k4, convr_gl, ctlag, n_gl)
614 call dtime%interpolate_scalar(tau1, cs_k4, convs_gl, ctlag, n_gl)
615 call dtime%interpolate_scalar(tau1, ct_k4, convt_gl, ctlag, n_gl)
616 call runge_kutta(fs, conv_k1, conv_k23, conv_k4, xh, xh_gl, &
617 coef, coef_gl, gll_to_gl, tau, dtau, &
618 n, nel, n_gl)
619 tau = tau1
620 end do
621 end do
622
623 end associate
624
625 end subroutine adv_oifs_compute_scalar
626
627end module adv_oifs
real
double real
Definition device_config.h:12
device::device_map
Map a Fortran array to a device (allocate and associate)
Definition device.F90:77
adv_oifs
Subroutines to add advection terms to the RHS of a transport equation.
Definition adv_oifs.f90:34
adv_oifs::adv_oifs_compute_scalar
subroutine adv_oifs_compute_scalar(this, vx, vy, vz, s, fs, xh, coef, n, dt)
Add the advection term for a scalar, i.e. , to the RHS.
Definition adv_oifs.f90:542
adv_oifs::adv_oifs_compute
subroutine adv_oifs_compute(this, vx, vy, vz, fx, fy, fz, xh, coef, n, dt)
Add the advection term for the fluid, i.e. , to the RHS using the OIFS method.
Definition adv_oifs.f90:412
adv_oifs::set_conv_velocity_fst
subroutine set_conv_velocity_fst(this, u, v, w)
Mapping the velocity fields to GL space and transforming them to the rst format.
Definition adv_oifs.f90:373
adv_oifs::adv_oifs_init
subroutine adv_oifs_init(this, lxd, coef, ctarget, ulag, vlag, wlag, dtlag, tlag, time_scheme, slag)
Constructor.
Definition adv_oifs.f90:128
adv_oifs::adv_oifs_free
subroutine adv_oifs_free(this)
Destructor.
Definition adv_oifs.f90:275
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_rzero
subroutine, public device_rzero(a_d, n, strm)
Zero a real vector.
Definition device_math.F90:232
device_math::device_addcol3s2
subroutine, public device_addcol3s2(a_d, b_d, c_d, s, n, strm)
Returns .
Definition device_math.F90:964
device
Device abstraction, common interface for various accelerators.
Definition device.F90:34
device::device_free
subroutine, public device_free(x_d)
Deallocate memory on the device.
Definition device.F90:219
field_list
Definition field_list.f90:1
field_series
Contains the field_serties_t type.
Definition field_series.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::copy
subroutine, public copy(a, b, n)
Copy a vector .
Definition math.f90:249
math::rzero
subroutine, public rzero(a, n)
Zero a real vector.
Definition math.f90:205
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_device
integer, parameter neko_bcknd_device
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::set_convect_rst
subroutine, public set_convect_rst(cr, cs, ct, cx, cy, cz, xh, coef)
Transforms the convecting velocity field to the rst form of the GL space.
Definition operators.f90:595
operators::runge_kutta
subroutine, public runge_kutta(phi, conv_k1, conv_k23, conv_k4, xh_gll, xh_gl, coef, coef_gl, gll_to_gl, tau, dtau, n, nel, n_gl)
Compute one step of Runge Kutta time interpolation for OIFS scheme.
Definition operators.f90:635
space
Defines a function space.
Definition space.f90:34
space::gl
integer, parameter, public gl
Definition space.f90:49
time_interpolator
Implements type time_interpolator_t.
Definition time_interpolator.f90:34
time_scheme_controller
Compound scheme for the advection and diffusion operators in a transport equation.
Definition time_scheme_controller.f90:35
time_scheme
Base class for time integration schemes.
Definition time_scheme.f90:61
adv_oifs::adv_oifs_t
Definition adv_oifs.f90:60
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:56
field::field_t
Definition field.f90:47
field_list::field_list_t
field_list_t, To be able to group fields together
Definition field_list.f90:15
field_series::field_series_t
Stores a series (sequence) of fields, logically connected to a base field, and arranged according to ...
Definition field_series.f90:42
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
time_interpolator::time_interpolator_t
Provides a tool to perform interpolation in time.
Definition time_interpolator.f90:48
time_scheme_controller::time_scheme_controller_t
Implements the logic to compute the time coefficients for the advection and diffusion operators in a ...
Definition time_scheme_controller.f90:71
max
#define max(a, b)
Definition tensor.cu:40