Neko 1.99.1
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 nullify(this%ulag)
302 nullify(this%vlag)
303 nullify(this%wlag)
304 nullify(this%slag)
305 nullify(this%ctlag)
306 nullify(this%dctlag)
307 nullify(this%oifs_scheme)
308 nullify(this%cr_k1)
309 nullify(this%cs_k1)
310 nullify(this%ct_k1)
311 nullify(this%cr_k23)
312 nullify(this%cs_k23)
313 nullify(this%ct_k23)
314 nullify(this%cr_k4)
315 nullify(this%cs_k4)
316 nullify(this%ct_k4)
317
318 if (allocated(this%cx)) then
319 deallocate(this%cx)
320 end if
321 if (allocated(this%cy)) then
322 deallocate(this%cy)
323 end if
324 if (allocated(this%cz)) then
325 deallocate(this%cz)
326 end if
327 if (c_associated(this%cx_d)) then
328 call device_free(this%cx_d)
329 end if
330 if (c_associated(this%cy_d)) then
331 call device_free(this%cy_d)
332 end if
333 if (c_associated(this%cz_d)) then
334 call device_free(this%cz_d)
335 end if
336
337 end subroutine adv_oifs_free
338
344 subroutine set_conv_velocity_fst(this, u, v, w)
345 implicit none
346 class(adv_oifs_t), intent(inout) :: this
347 type(field_t), intent(inout) :: u, v, w
348 integer :: i, nel, n_GL, idx, idy, idz
349
350 nel = this%coef_GLL%msh%nelv
351 n_gl = nel*this%Xh_GL%lxyz
352
353 call this%convr_GL%update()
354 call this%convs_GL%update()
355 call this%convt_GL%update()
356
357 call this%GLL_to_GL%map(this%cx, u%x, nel, this%Xh_GL)
358 call this%GLL_to_GL%map(this%cy, v%x, nel, this%Xh_GL)
359 call this%GLL_to_GL%map(this%cz, w%x, nel, this%Xh_GL)
360
361 call set_convect_rst(this%cr_GL, this%cs_GL, this%ct_GL, &
362 this%cx, this%cy, this%cz, this%Xh_GL, this%coef_GL)
363
364 this%convr_GL%f = this%cr_GL
365 this%convs_GL%f = this%cs_GL
366 this%convt_GL%f = this%ct_GL
367
368 end subroutine set_conv_velocity_fst
369
370
383 subroutine adv_oifs_compute(this, vx, vy, vz, fx, fy, fz, Xh, coef, n, dt)
384 implicit none
385 class(adv_oifs_t), intent(inout) :: this
386 type(field_t), intent(inout) :: vx, vy, vz
387 type(field_t), intent(inout) :: fx, fy, fz
388 type(space_t), intent(in) :: Xh
389 type(coef_t), intent(in) :: coef
390 integer, intent(in) :: n
391 real(kind=rp), intent(in), optional :: dt
392 real(kind=rp) :: tau, tau1, th, dtau
393 integer :: i, ilag, itau, nel, n_GL
394
395 nel = coef%msh%nelv
396 n_gl = nel * this%Xh_GL%lxyz
397
398 associate(ulag => this%ulag, vlag => this%vlag, wlag => this%wlag, &
399 ctlag => this%ctlag, dctlag => this%dctlag, dtime => this%dtime, &
400 xh_gl => this%Xh_GL, coef_gl => this%coef_GL, ntaubd => this%ntaubd, &
401 gll_to_gl => this%GLL_to_GL, oifs_scheme => this%oifs_scheme, &
402 cr_k1 => this%cr_K1, cs_k1 => this%cs_K1, ct_k1 => this%ct_K1, &
403 cr_k23 => this%cr_K23, cs_k23 => this%cs_K23, ct_k23 => this%ct_K23, &
404 cr_k4 => this%cr_K4, cs_k4 => this%cs_K4, ct_k4 => this%ct_K4, &
405 convr_gl => this%convr_GL, convs_gl => this%convs_GL, &
406 convt_gl => this%convt_GL, conv_k1 => this%conv_k1, &
407 conv_k23 => this%conv_k23, conv_k4 => this%conv_k4)
408
409 call dtime%init(oifs_scheme%ndiff)
410
411 tau = ctlag(oifs_scheme%ndiff)
412
413 call this%set_conv_velocity_fst(vx, vy, vz)
414
415 if (neko_bcknd_device .eq. 1) then
416 call device_rzero(fx%x_d,n)
417 call device_rzero(fy%x_d,n)
418 call device_rzero(fz%x_d,n)
419 else
420 call rzero(fx%x,n)
421 call rzero(fy%x,n)
422 call rzero(fz%x,n)
423 end if
424
425 do ilag = oifs_scheme%ndiff, 1, -1
426 if (neko_bcknd_device .eq. 1) then
427 if (ilag .eq. 1) then
428 call device_addcol3s2(fx%x_d, vx%x_d, coef%B_d, &
429 oifs_scheme%diffusion_coeffs(2), n)
430 call device_addcol3s2(fy%x_d, vy%x_d, coef%B_d, &
431 oifs_scheme%diffusion_coeffs(2), n)
432 call device_addcol3s2(fz%x_d, vz%x_d, coef%B_d, &
433 oifs_scheme%diffusion_coeffs(2), n)
434 else
435 call device_addcol3s2(fx%x_d, ulag%lf(ilag-1)%x_d, coef%B_d, &
436 oifs_scheme%diffusion_coeffs(ilag+1), n)
437 call device_addcol3s2(fy%x_d, vlag%lf(ilag-1)%x_d, coef%B_d, &
438 oifs_scheme%diffusion_coeffs(ilag+1), n)
439 call device_addcol3s2(fz%x_d, wlag%lf(ilag-1)%x_d, coef%B_d, &
440 oifs_scheme%diffusion_coeffs(ilag+1), n)
441 end if
442 else
443 if (ilag .eq. 1) then
444 do i = 1, n
445 fx%x(i,1,1,1) = fx%x(i,1,1,1) + &
446 oifs_scheme%diffusion_coeffs(2) &
447 * vx%x(i,1,1,1) * coef%B(i,1,1,1)
448 fy%x(i,1,1,1) = fy%x(i,1,1,1) + &
449 oifs_scheme%diffusion_coeffs(2) &
450 * vy%x(i,1,1,1) * coef%B(i,1,1,1)
451 fz%x(i,1,1,1) = fz%x(i,1,1,1) + &
452 oifs_scheme%diffusion_coeffs(2) &
453 * vz%x(i,1,1,1) * coef%B(i,1,1,1)
454 end do
455 else
456 do i = 1, n
457 fx%x(i,1,1,1) = fx%x(i,1,1,1) + &
458 oifs_scheme%diffusion_coeffs(ilag+1) &
459 * ulag%lf(ilag-1)%x(i,1,1,1) &
460 * coef%B(i,1,1,1)
461 fy%x(i,1,1,1) = fy%x(i,1,1,1) + &
462 oifs_scheme%diffusion_coeffs(ilag+1) &
463 * vlag%lf(ilag-1)%x(i,1,1,1) &
464 * coef%B(i,1,1,1)
465 fz%x(i,1,1,1) = fz%x(i,1,1,1) + &
466 oifs_scheme%diffusion_coeffs(ilag+1) &
467 * wlag%lf(ilag-1)%x(i,1,1,1) &
468 * coef%B(i,1,1,1)
469 end do
470 end if
471 end if
472 dtau = dctlag(ilag)/real(ntaubd)
473 do itau = 1, ntaubd
474 th = tau + dtau/2.
475 tau1 = tau + dtau
476 call dtime%interpolate_scalar(tau, cr_k1, convr_gl, ctlag, n_gl)
477 call dtime%interpolate_scalar(tau, cs_k1, convs_gl, ctlag, n_gl)
478 call dtime%interpolate_scalar(tau, ct_k1, convt_gl, ctlag, n_gl)
479 call dtime%interpolate_scalar(th, cr_k23, convr_gl, ctlag, n_gl)
480 call dtime%interpolate_scalar(th, cs_k23, convs_gl, ctlag, n_gl)
481 call dtime%interpolate_scalar(th, ct_k23, convt_gl, ctlag, n_gl)
482 call dtime%interpolate_scalar(tau1, cr_k4, convr_gl, ctlag, n_gl)
483 call dtime%interpolate_scalar(tau1, cs_k4, convs_gl, ctlag, n_gl)
484 call dtime%interpolate_scalar(tau1, ct_k4, convt_gl, ctlag, n_gl)
485 call runge_kutta(fx, conv_k1, conv_k23, conv_k4, xh, xh_gl, &
486 coef, coef_gl, gll_to_gl, tau, dtau, &
487 n, nel, n_gl)
488 call runge_kutta(fy, conv_k1, conv_k23, conv_k4, xh, xh_gl, &
489 coef, coef_gl, gll_to_gl, tau, dtau, &
490 n, nel, n_gl)
491 call runge_kutta(fz, conv_k1, conv_k23, conv_k4, xh, xh_gl, &
492 coef, coef_gl, gll_to_gl, tau, dtau, &
493 n, nel, n_gl)
494 tau = tau1
495 end do
496 end do
497
498 end associate
499
500 end subroutine adv_oifs_compute
513 subroutine adv_oifs_compute_scalar(this, vx, vy, vz, s, fs, Xh, coef, n, dt)
514 implicit none
515 class(adv_oifs_t), intent(inout) :: this
516 type(field_t), intent(inout) :: vx, vy, vz
517 type(field_t), intent(inout) :: fs
518 type(field_t), intent(inout) :: s
519 type(space_t), intent(in) :: Xh
520 type(coef_t), intent(in) :: coef
521 integer, intent(in) :: n
522 real(kind=rp), intent(in), optional :: dt
523 real(kind=rp) :: tau, tau1, th, dtau
524 integer :: i, ilag, itau, nel, n_GL
525 nel = coef%msh%nelv
526 n_gl = nel * this%Xh_GL%lxyz
527
528 associate(slag => this%slag, ctlag => this%ctlag, dctlag => this%dctlag, &
529 dtime => this%dtime, xh_gl => this%Xh_GL, coef_gl => this%coef_GL, &
530 ntaubd => this%ntaubd, gll_to_gl => this%GLL_to_GL, &
531 oifs_scheme => this%oifs_scheme, cr_k1 => this%cr_K1, &
532 cs_k1 => this%cs_K1, ct_k1 => this%ct_K1, cr_k23 => this%cr_K23, &
533 cs_k23 => this%cs_K23, ct_k23 => this%ct_K23, cr_k4 => this%cr_K4, &
534 cs_k4 => this%cs_K4, ct_k4 => this%ct_K4, &
535 convr_gl => this%convr_GL, convs_gl => this%convs_GL, &
536 convt_gl => this%convt_GL, conv_k1 => this%conv_k1, &
537 conv_k23 => this%conv_k23, conv_k4 => this%conv_k4)
538
539 call dtime%init(oifs_scheme%ndiff)
540
541 tau = ctlag(oifs_scheme%ndiff)
542
543 call this%set_conv_velocity_fst(vx, vy, vz)
544
545 if (neko_bcknd_device .eq. 1) then
546 call device_rzero(fs%x_d,n)
547 else
548 call rzero(fs%x,n)
549 end if
550
551 do ilag = oifs_scheme%ndiff, 1, -1
552 if (neko_bcknd_device .eq. 1) then
553 if (ilag .eq. 1) then
554 call device_addcol3s2(fs%x_d, s%x_d, coef%B_d, &
555 oifs_scheme%diffusion_coeffs(2), n)
556 else
557 call device_addcol3s2(fs%x_d, slag%lf(ilag-1)%x_d, coef%B_d, &
558 oifs_scheme%diffusion_coeffs(ilag+1), n)
559 end if
560 else
561 if (ilag .eq. 1) then
562 do i = 1, n
563 fs%x(i,1,1,1) = fs%x(i,1,1,1) + &
564 oifs_scheme%diffusion_coeffs(2) &
565 * s%x(i,1,1,1) * coef%B(i,1,1,1)
566 end do
567 else
568 do i = 1, n
569 fs%x(i,1,1,1) = fs%x(i,1,1,1) + &
570 oifs_scheme%diffusion_coeffs(ilag+1) &
571 * slag%lf(ilag-1)%x(i,1,1,1) * coef%B(i,1,1,1)
572 end do
573 end if
574 end if
575 dtau = dctlag(ilag)/real(ntaubd)
576 do itau = 1, ntaubd
577 th = tau + dtau/2.
578 tau1 = tau + dtau
579 call dtime%interpolate_scalar(tau, cr_k1, convr_gl, ctlag, n_gl)
580 call dtime%interpolate_scalar(tau, cs_k1, convs_gl, ctlag, n_gl)
581 call dtime%interpolate_scalar(tau, ct_k1, convt_gl, ctlag, n_gl)
582 call dtime%interpolate_scalar(th, cr_k23, convr_gl, ctlag, n_gl)
583 call dtime%interpolate_scalar(th, cs_k23, convs_gl, ctlag, n_gl)
584 call dtime%interpolate_scalar(th, ct_k23, convt_gl, ctlag, n_gl)
585 call dtime%interpolate_scalar(tau1, cr_k4, convr_gl, ctlag, n_gl)
586 call dtime%interpolate_scalar(tau1, cs_k4, convs_gl, ctlag, n_gl)
587 call dtime%interpolate_scalar(tau1, ct_k4, convt_gl, ctlag, n_gl)
588 call runge_kutta(fs, conv_k1, conv_k23, conv_k4, xh, xh_gl, &
589 coef, coef_gl, gll_to_gl, tau, dtau, &
590 n, nel, n_gl)
591 tau = tau1
592 end do
593 end do
594
595 end associate
596
597 end subroutine adv_oifs_compute_scalar
598
599end 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:514
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:384
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:345
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:587
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:627
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:55
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