48 use,
intrinsic :: iso_c_binding, only: c_ptr, c_null_ptr
64 type(
coef_t),
pointer :: coef_gll => null()
70 type(
space_t),
pointer :: xh_gll => null()
76 real(kind=
rp),
pointer :: ctlag(:) => null()
78 real(kind=
rp),
pointer :: dctlag(:) => null()
82 real(kind=
rp),
allocatable :: cx(:), cy(:), cz(:)
84 real(kind=
rp),
allocatable :: c(:,:)
114 dtlag, tlag, time_scheme, slag)
117 integer,
intent(in) :: lxd
118 type(
coef_t),
target :: coef
119 real(kind=
rp),
intent(in) :: ctarget
121 real(kind=
rp),
target,
intent(in) :: dtlag(10)
122 real(kind=
rp),
target,
intent(in) :: tlag(10)
125 integer :: nel, n_GL, n_GL_t, n, idx, idy, idz
126 real(kind=
rp) :: max_cfl_rk4
130 this%ntaubd =
max(int(ctarget/max_cfl_rk4),1)
132 call this%Xh_GL%init(
gl, lxd, lxd, lxd)
133 this%Xh_GLL => coef%Xh
134 this%coef_GLL => coef
135 call this%GLL_to_GL%init(this%Xh_GL, this%Xh_GLL)
137 call this%coef_GL%init(this%Xh_GL, coef%msh)
140 n_gl = nel*this%Xh_GL%lxyz
149 call this%GLL_to_GL%map(this%coef_GL%drdx, coef%drdx, nel, this%Xh_GL)
150 call this%GLL_to_GL%map(this%coef_GL%dsdx, coef%dsdx, nel, this%Xh_GL)
151 call this%GLL_to_GL%map(this%coef_GL%dtdx, coef%dtdx, nel, this%Xh_GL)
152 call this%GLL_to_GL%map(this%coef_GL%drdy, coef%drdy, nel, this%Xh_GL)
153 call this%GLL_to_GL%map(this%coef_GL%dsdy, coef%dsdy, nel, this%Xh_GL)
154 call this%GLL_to_GL%map(this%coef_GL%dtdy, coef%dtdy, nel, this%Xh_GL)
155 call this%GLL_to_GL%map(this%coef_GL%drdz, coef%drdz, nel, this%Xh_GL)
156 call this%GLL_to_GL%map(this%coef_GL%dsdz, coef%dsdz, nel, this%Xh_GL)
157 call this%GLL_to_GL%map(this%coef_GL%dtdz, coef%dtdz, nel, this%Xh_GL)
159 allocate(this%cx(n_gl))
160 allocate(this%cy(n_gl))
161 allocate(this%cz(n_gl))
162 allocate(this%c(n_gl_t, 3))
164 call this%dtime%init(1)
174 call this%GLL_to_GL%map(this%cx, this%ulag%f%x, nel, this%Xh_GL)
175 call this%GLL_to_GL%map(this%cy, this%vlag%f%x, nel, this%Xh_GL)
176 call this%GLL_to_GL%map(this%cz, this%wlag%f%x, nel, this%Xh_GL)
180 this%cx, this%cy, this%cz, this%Xh_GL, this%coef_GL)
182 call this%GLL_to_GL%map(this%cx, this%ulag%lf(1)%x, nel, this%Xh_GL)
183 call this%GLL_to_GL%map(this%cy, this%vlag%lf(1)%x, nel, this%Xh_GL)
184 call this%GLL_to_GL%map(this%cz, this%wlag%lf(1)%x, nel, this%Xh_GL)
187 this%cx, this%cy, this%cz, this%Xh_GL, this%coef_GL)
189 call this%GLL_to_GL%map(this%cx, this%ulag%lf(2)%x, nel, this%Xh_GL)
190 call this%GLL_to_GL%map(this%cy, this%vlag%lf(2)%x, nel, this%Xh_GL)
191 call this%GLL_to_GL%map(this%cz, this%wlag%lf(2)%x, nel, this%Xh_GL)
194 this%cx, this%cy, this%cz, this%Xh_GL, this%coef_GL)
197 if (
present(slag))
then
207 call this%coef_GL%free()
209 nullify(this%coef_GLL)
211 call this%GLL_to_GL%free()
213 call this%Xh_GL%free()
217 call this%dtime%free()
225 nullify(this%oifs_scheme)
227 if (
allocated(this%cx))
then
230 if (
allocated(this%cy))
then
233 if (
allocated(this%cz))
then
236 if (
allocated(this%c))
then
250 type(
field_t),
intent(inout) :: u, v, w
251 integer :: i, nel, n_GL, n_GL_t, idx, idy, idz
253 nel = this%coef_GLL%msh%nelv
254 n_gl = nel*this%Xh_GL%lxyz
256 call copy(this%c(:,3), this%c(:,2), n_gl_t)
257 call copy(this%c(:,2), this%c(:,1), n_gl_t)
259 call this%GLL_to_GL%map(this%cx, u%x, nel, this%Xh_GL)
260 call this%GLL_to_GL%map(this%cy, v%x, nel, this%Xh_GL)
261 call this%GLL_to_GL%map(this%cz, w%x, nel, this%Xh_GL)
268 this%cx, this%cy, this%cz, this%Xh_GL, this%coef_GL)
285 subroutine adv_oifs_compute(this, vx, vy, vz, fx, fy, fz, Xh, coef, n, dt)
288 type(
field_t),
intent(inout) :: vx, vy, vz
289 type(
field_t),
intent(inout) :: fx, fy, fz
290 type(
space_t),
intent(inout) :: Xh
291 type(
coef_t),
intent(inout) :: coef
292 integer,
intent(in) :: n
293 real(kind=
rp),
intent(in),
optional :: dt
295 real(kind=
rp) :: tau, tau1, th, dtau
296 integer :: i, ilag, itau, nel, n_GL, n_GL_t
297 real(kind=
rp),
dimension(3 * coef%msh%nelv * this%Xh_GL%lxyz) :: c_r1
298 real(kind=
rp),
dimension(3 * coef%msh%nelv * this%Xh_GL%lxyz) :: c_r23
299 real(kind=
rp),
dimension(3 * coef%msh%nelv * this%Xh_GL%lxyz) :: c_r4
300 real(kind=
rp),
parameter :: eps = 1e-10
303 n_gl = nel * this%Xh_GL%lxyz
306 associate(ulag => this%ulag, vlag => this%vlag, wlag => this%wlag, &
307 ctlag => this%ctlag, dctlag => this%dctlag, dtime => this%dtime, &
308 oifs_scheme => this%oifs_scheme, ntaubd => this%ntaubd, c => this%c)
310 call dtime%init(oifs_scheme%ndiff)
312 tau = ctlag(oifs_scheme%ndiff)
318 call this%set_conv_velocity_fst(vx, vy, vz)
320 do ilag = oifs_scheme%ndiff, 1, -1
322 if (ilag .eq. 1)
then
324 fx%x(i,1,1,1) = fx%x(i,1,1,1) + &
325 oifs_scheme%diffusion_coeffs(2) &
326 * vx%x(i,1,1,1) * coef%B(i,1,1,1)
327 fy%x(i,1,1,1) = fy%x(i,1,1,1) + &
328 oifs_scheme%diffusion_coeffs(2) &
329 * vy%x(i,1,1,1) * coef%B(i,1,1,1)
330 fz%x(i,1,1,1) = fz%x(i,1,1,1) + &
331 oifs_scheme%diffusion_coeffs(2) &
332 * vz%x(i,1,1,1) * coef%B(i,1,1,1)
336 fx%x(i,1,1,1) = fx%x(i,1,1,1) + &
337 oifs_scheme%diffusion_coeffs(ilag+1) &
338 * ulag%lf(ilag-1)%x(i,1,1,1) &
340 fy%x(i,1,1,1) = fy%x(i,1,1,1) + &
341 oifs_scheme%diffusion_coeffs(ilag+1) &
342 * vlag%lf(ilag-1)%x(i,1,1,1) &
344 fz%x(i,1,1,1) = fz%x(i,1,1,1) + &
345 oifs_scheme%diffusion_coeffs(ilag+1) &
346 * wlag%lf(ilag-1)%x(i,1,1,1) &
351 if (dctlag(ilag) .lt. eps)
then
352 dtau = dt/
real(ntaubd)
354 dtau = dctlag(ilag)/
real(ntaubd)
359 call dtime%interpolate_scalar(tau, c_r1, c, ctlag, n_gl_t)
360 call dtime%interpolate_scalar(th, c_r23, c, ctlag, n_gl_t)
361 call dtime%interpolate_scalar(tau1, c_r4, c, ctlag, n_gl_t)
362 call runge_kutta(fx%x, c_r1, c_r23, c_r4, xh, this%Xh_GL, coef, &
363 this%coef_GL, this%GLL_to_GL, &
364 tau, dtau, n, nel, n_gl)
365 call runge_kutta(fy%x, c_r1, c_r23, c_r4, xh, this%Xh_GL, coef, &
366 this%coef_GL, this%GLL_to_GL, &
367 tau, dtau, n, nel, n_gl)
368 call runge_kutta(fz%x, c_r1, c_r23, c_r4, xh, this%Xh_GL, coef, &
369 this%coef_GL, this%GLL_to_GL, &
370 tau, dtau, n, nel, n_gl)
392 type(field_t),
intent(inout) :: vx, vy, vz
393 type(field_t),
intent(inout) :: fs
394 type(field_t),
intent(inout) :: s
395 type(space_t),
intent(inout) :: Xh
396 type(coef_t),
intent(inout) :: coef
397 integer,
intent(in) :: n
398 real(kind=rp),
intent(in),
optional :: dt
399 real(kind=rp) :: tau, tau1, th, dtau
400 integer :: i, ilag, itau, nel, n_GL, n_GL_t
401 real(kind=rp),
dimension(3 * coef%msh%nelv * this%Xh_GL%lxyz) :: c_r1
402 real(kind=rp),
dimension(3 * coef%msh%nelv * this%Xh_GL%lxyz) :: c_r23
403 real(kind=rp),
dimension(3 * coef%msh%nelv * this%Xh_GL%lxyz) :: c_r4
404 real(kind=rp),
parameter :: eps = 1e-10
406 n_gl = nel * this%Xh_GL%lxyz
409 associate(slag => this%slag, ctlag => this%ctlag, &
410 dctlag => this%dctlag, dtime => this%dtime, &
411 oifs_scheme => this%oifs_scheme, ntaubd => this%ntaubd, c => this%c)
413 call dtime%init(oifs_scheme%ndiff)
415 tau = ctlag(oifs_scheme%ndiff)
419 call this%set_conv_velocity_fst(vx, vy, vz)
421 do ilag = oifs_scheme%ndiff, 1, -1
422 if (ilag .eq. 1)
then
424 fs%x(i,1,1,1) = fs%x(i,1,1,1) + &
425 oifs_scheme%diffusion_coeffs(2) &
426 * s%x(i,1,1,1) * coef%B(i,1,1,1)
430 fs%x(i,1,1,1) = fs%x(i,1,1,1) + &
431 oifs_scheme%diffusion_coeffs(ilag+1) &
432 * slag%lf(ilag-1)%x(i,1,1,1) * coef%B(i,1,1,1)
436 if (dctlag(ilag) .lt. eps)
then
437 dtau = dt/
real(ntaubd)
439 dtau = dctlag(ilag)/
real(ntaubd)
444 call dtime%interpolate_scalar(tau, c_r1, c, ctlag, n_gl_t)
445 call dtime%interpolate_scalar(th, c_r23, c, ctlag, n_gl_t)
446 call dtime%interpolate_scalar(tau1, c_r4, c, ctlag, n_gl_t)
447 call runge_kutta(fs%x, c_r1, c_r23, c_r4, xh, this%Xh_GL, coef, &
448 this%coef_GL, this%GLL_to_GL, &
449 tau, dtau, n, nel, n_gl)
Return the device pointer for an associated Fortran array.
Map a Fortran array to a device (allocate and associate)
Subroutines to add advection terms to the RHS of a transport equation.
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.
subroutine set_conv_velocity_fst(this, u, v, w)
Mapping the velocity fields to GL space and transforming them to the rst format.
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.
subroutine adv_oifs_init(this, lxd, coef, ctarget, ulag, vlag, wlag, dtlag, tlag, time_scheme, slag)
Constructor.
subroutine adv_oifs_free(this)
Destructor.
Subroutines to add advection terms to the RHS of a transport equation.
Device abstraction, common interface for various accelerators.
Routines to interpolate between different spaces.
subroutine, public copy(a, b, n)
Copy a vector .
subroutine, public rzero(a, n)
Zero a real vector.
integer, parameter neko_bcknd_sx
integer, parameter neko_bcknd_xsmm
integer, parameter, public rp
Global precision used in computations.
subroutine, public runge_kutta(phi, c_r1, c_r23, c_r4, 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.
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.
Defines a function space.
integer, parameter, public gl
Implements type time_interpolator_t.
Compound scheme for the advection and diffusion operators in a transport equation.
Base class for time integration schemes.
Base abstract type for computing the advection operator.
Coefficients defined on a given (mesh, ) tuple. Arrays use indices (i,j,k,e): element e,...
Interpolation between two space::space_t.
The function space for the SEM solution fields.
Provides a tool to perform interpolation in time.
Implements the logic to compute the time coefficients for the advection and diffusion operators in a ...