Neko  0.9.99
A portable framework for high-order spectral element flow simulations
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 !
34 module 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
43  use interpolation, only: interpolator_t
45  use field_series, only: field_series_t
47  use device, only: device_map, device_get_ptr
48  use, intrinsic :: iso_c_binding, only: c_ptr, c_null_ptr
49  implicit none
50  private
51 
52  !! Type encapsulating operator-integration-factor splitting advection
53  !! routines with dealiasing applied
54  !! Literature:
55  !! https://www.mcs.anl.gov/~fischer/nek5000/oifs.pdf
56  !! https://publications.anl.gov/anlpubs/2017/12/140626.pdf
57  !! https://dl.acm.org/doi/abs/10.1007/BF01063118
58  type, public, extends(advection_t) :: adv_oifs_t
60  integer :: ntaubd
62  type(coef_t) :: coef_gl
64  type(coef_t), pointer :: coef_gll => null()
66  type(interpolator_t) :: gll_to_gl
68  type(space_t) :: xh_gl
70  type(space_t), pointer :: xh_gll => null()
72  type(time_interpolator_t) :: dtime
74  type(field_series_t), pointer :: ulag, vlag, wlag, slag => null()
76  real(kind=rp), pointer :: ctlag(:) => null()
78  real(kind=rp), pointer :: dctlag(:) => null()
80  type(time_scheme_controller_t), pointer :: oifs_scheme => null()
82  real(kind=rp), allocatable :: cx(:), cy(:), cz(:)
84  real(kind=rp), allocatable :: c(:,:)
85  contains
88  procedure, pass(this) :: compute => adv_oifs_compute
91  procedure, pass(this) :: compute_scalar => adv_oifs_compute_scalar
93  procedure, pass(this) :: init => adv_oifs_init
95  procedure, pass(this) :: set_conv_velocity_fst
97  procedure, pass(this) :: free => adv_oifs_free
98  end type adv_oifs_t
99 
100 contains
101 
113  subroutine adv_oifs_init(this, lxd, coef, ctarget, ulag, vlag, wlag, &
114  dtlag, tlag, time_scheme, slag)
115  implicit none
116  class(adv_oifs_t) :: this
117  integer, intent(in) :: lxd
118  type(coef_t), target :: coef
119  real(kind=rp), intent(in) :: ctarget
120  type(field_series_t), target, intent(in) :: ulag, vlag, wlag
121  real(kind=rp), target, intent(in) :: dtlag(10)
122  real(kind=rp), target, intent(in) :: tlag(10)
123  type(time_scheme_controller_t), target, intent(in) :: time_scheme
124  type(field_series_t), target, optional :: slag
125  integer :: nel, n_GL, n_GL_t, n, idx, idy, idz
126  real(kind=rp) :: max_cfl_rk4
127 
128  ! stability limit for RK4 including safety factor
129  max_cfl_rk4 = 2.0
130  this%ntaubd = max(int(ctarget/max_cfl_rk4),1)
131 
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)
136 
137  call this%coef_GL%init(this%Xh_GL, coef%msh)
138 
139  nel = coef%msh%nelv
140  n_gl = nel*this%Xh_GL%lxyz
141  n = nel*coef%Xh%lxyz
142  n_gl_t = 3 * n_gl
143 
144  idx = 1
145  idy = idx + n_gl
146  idz = idy + n_gl
147 
148 
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)
158 
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))
163 
164  call this%dtime%init(1)
165  this%ulag => ulag
166  this%vlag => vlag
167  this%wlag => wlag
168  this%ctlag => tlag
169  this%dctlag => dtlag
170  this%oifs_scheme => time_scheme
171 
172  ! Initializing the convecting fields
173  ! Map the velocity fields from GLL space to GL space
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)
177 
178  ! Set the convecting field in the rst format
179  call set_convect_rst(this%c(idx,1), this%c(idy,1), this%c(idz,1), &
180  this%cx, this%cy, this%cz, this%Xh_GL, this%coef_GL)
181  ! Repeat for previous time-steps
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)
185 
186  call set_convect_rst(this%c(idx,2), this%c(idy,2), this%c(idz,2), &
187  this%cx, this%cy, this%cz, this%Xh_GL, this%coef_GL)
188 
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)
192 
193  call set_convect_rst(this%c(idx,3), this%c(idy,3), this%c(idz,3), &
194  this%cx, this%cy, this%cz, this%Xh_GL, this%coef_GL)
195 
196  ! Initilize the lagged scalar field, if present.
197  if (present(slag)) then
198  this%slag => slag
199  end if
200 
201  end subroutine adv_oifs_init
202 
204  subroutine adv_oifs_free(this)
205  class(adv_oifs_t), intent(inout) :: this
206 
207  call this%coef_GL%free()
208 
209  nullify(this%coef_GLL)
210 
211  call this%GLL_to_GL%free()
212 
213  call this%Xh_GL%free()
214 
215  nullify(this%Xh_GLL)
216 
217  call this%dtime%free()
218 
219  nullify(this%ulag)
220  nullify(this%vlag)
221  nullify(this%wlag)
222  nullify(this%slag)
223  nullify(this%ctlag)
224  nullify(this%dctlag)
225  nullify(this%oifs_scheme)
226 
227  if (allocated(this%cx)) then
228  deallocate(this%cx)
229  end if
230  if (allocated(this%cy)) then
231  deallocate(this%cy)
232  end if
233  if (allocated(this%cz)) then
234  deallocate(this%cz)
235  end if
236  if (allocated(this%c)) then
237  deallocate(this%c)
238  end if
239 
240  end subroutine adv_oifs_free
241 
247  subroutine set_conv_velocity_fst(this, u, v, w)
248  implicit none
249  class(adv_oifs_t), intent(inout) :: this
250  type(field_t), intent(inout) :: u, v, w
251  integer :: i, nel, n_GL, n_GL_t, idx, idy, idz
252 
253  nel = this%coef_GLL%msh%nelv
254  n_gl = nel*this%Xh_GL%lxyz
255  n_gl_t = 3 * n_gl
256  call copy(this%c(:,3), this%c(:,2), n_gl_t)
257  call copy(this%c(:,2), this%c(:,1), n_gl_t)
258 
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)
262 
263  idx = 1
264  idy = idx + n_gl
265  idz = idy + n_gl
266 
267  call set_convect_rst(this%c(idx,1), this%c(idy,1), this%c(idz,1), &
268  this%cx, this%cy, this%cz, this%Xh_GL, this%coef_GL)
269 
270  end subroutine set_conv_velocity_fst
271 
272 
285  subroutine adv_oifs_compute(this, vx, vy, vz, fx, fy, fz, Xh, coef, n, dt)
286  implicit none
287  class(adv_oifs_t), intent(inout) :: this
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
294 
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
301 
302  nel = coef%msh%nelv
303  n_gl = nel * this%Xh_GL%lxyz
304  n_gl_t = 3 * n_gl
305 
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)
309 
310  call dtime%init(oifs_scheme%ndiff)
311 
312  tau = ctlag(oifs_scheme%ndiff)
313 
314  call rzero(fx%x,n)
315  call rzero(fy%x,n)
316  call rzero(fz%x,n)
317 
318  call this%set_conv_velocity_fst(vx, vy, vz)
319 
320  do ilag = oifs_scheme%ndiff, 1, -1
321 
322  if (ilag .eq. 1) then
323  do i = 1, n
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)
333  end do
334  else
335  do i = 1, n
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) &
339  * coef%B(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) &
343  * coef%B(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) &
347  * coef%B(i,1,1,1)
348  end do
349  end if
350 
351  if (dctlag(ilag) .lt. eps) then
352  dtau = dt/real(ntaubd)
353  else
354  dtau = dctlag(ilag)/real(ntaubd)
355  end if
356  do itau = 1, ntaubd
357  th = tau + dtau/2.
358  tau1 = tau + dtau
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)
371  tau = tau1
372  end do
373  end do
374  end associate
375 
376  end subroutine adv_oifs_compute
389  subroutine adv_oifs_compute_scalar(this, vx, vy, vz, s, fs, Xh, coef, n, dt)
390  implicit none
391  class(adv_oifs_t), intent(inout) :: this
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
405  nel = coef%msh%nelv
406  n_gl = nel * this%Xh_GL%lxyz
407  n_gl_t = 3 * n_gl
408 
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)
412 
413  call dtime%init(oifs_scheme%ndiff)
414 
415  tau = ctlag(oifs_scheme%ndiff)
416 
417  call rzero(fs%x,n)
418 
419  call this%set_conv_velocity_fst(vx, vy, vz)
420 
421  do ilag = oifs_scheme%ndiff, 1, -1
422  if (ilag .eq. 1) then
423  do i = 1, n
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)
427  end do
428  else
429  do i = 1, n
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)
433  end do
434  end if
435 
436  if (dctlag(ilag) .lt. eps) then
437  dtau = dt/real(ntaubd)
438  else
439  dtau = dctlag(ilag)/real(ntaubd)
440  end if
441  do itau = 1, ntaubd
442  th = tau + dtau/2.
443  tau1 = tau + dtau
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)
450  tau = tau1
451  end do
452  end do
453 
454  end associate
455 
456  end subroutine adv_oifs_compute_scalar
457 
458 end module adv_oifs
double real
Definition: device_config.h:12
Return the device pointer for an associated Fortran array.
Definition: device.F90:81
Map a Fortran array to a device (allocate and associate)
Definition: device.F90:57
Subroutines to add advection terms to the RHS of a transport equation.
Definition: adv_oifs.f90:34
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:286
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:248
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:390
subroutine adv_oifs_init(this, lxd, coef, ctarget, ulag, vlag, wlag, dtlag, tlag, time_scheme, slag)
Constructor.
Definition: adv_oifs.f90:115
subroutine adv_oifs_free(this)
Destructor.
Definition: adv_oifs.f90:205
Subroutines to add advection terms to the RHS of a transport equation.
Definition: advection.f90:34
Coefficients.
Definition: coef.f90:34
Device abstraction, common interface for various accelerators.
Definition: device.F90:34
Stores a series fields.
Defines a field.
Definition: field.f90:34
Routines to interpolate between different spaces.
Definition: math.f90:60
subroutine, public copy(a, b, n)
Copy a vector .
Definition: math.f90:239
subroutine, public rzero(a, n)
Zero a real vector.
Definition: math.f90:195
Build configurations.
Definition: neko_config.f90:34
integer, parameter neko_bcknd_sx
Definition: neko_config.f90:39
integer, parameter neko_bcknd_xsmm
Definition: neko_config.f90:40
integer, parameter, public rp
Global precision used in computations.
Definition: num_types.f90:12
Operators.
Definition: operators.f90:34
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.
Definition: operators.f90:562
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:528
Defines a function space.
Definition: space.f90:34
integer, parameter, public gl
Definition: space.f90:48
Implements type time_interpolator_t.
Compound scheme for the advection and diffusion operators in a transport equation.
Base class for time integration schemes.
Definition: time_scheme.f90:61
Base abstract type for computing the advection operator.
Definition: advection.f90:46
Coefficients defined on a given (mesh, ) tuple. Arrays use indices (i,j,k,e): element e,...
Definition: coef.f90:55
Interpolation between two space::space_t.
The function space for the SEM solution fields.
Definition: space.f90:62
Provides a tool to perform interpolation in time.
Implements the logic to compute the time coefficients for the advection and diffusion operators in a ...
#define max(a, b)
Definition: tensor.cu:40