time_scheme_controller.f90

Go to the documentation of this file.

! Copyright (c) 2023, The Neko Authors
! All rights reserved.
!
! Redistribution and use in source and binary forms, with or without
! modification, are permitted provided that the following conditions
! are met:
!
!   * Redistributions of source code must retain the above copyright
!     notice, this list of conditions and the following disclaimer.
!
!   * Redistributions in binary form must reproduce the above
!     copyright notice, this list of conditions and the following
!     disclaimer in the documentation and/or other materials provided
!     with the distribution.
!
!   * Neither the name of the authors nor the names of its
!     contributors may be used to endorse or promote products derived
!     from this software without specific prior written permission.
!
! THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
! "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
! LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS
! FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE
! COPYRIGHT OWNER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT,
! INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING,
! BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES;
! LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER
! CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT
! LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN
! ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
! POSSIBILITY OF SUCH DAMAGE.
!
module time_scheme_controller
  use neko_config, only : neko_bcknd_device
  use num_types, only : rp
  use bdf_time_scheme, only : bdf_time_scheme_t
  use ext_time_scheme, only : ext_time_scheme_t
  use ab_time_scheme, only : ab_time_scheme_t
  use device, only : device_free, device_map, device_memcpy, device_sync, &
       host_to_device
  use vector, only : vector_t
  use, intrinsic :: iso_c_binding
  implicit none
  private
  type, public :: time_scheme_controller_t
     type(ext_time_scheme_t) :: ext
     type(ab_time_scheme_t) :: ab
     type(bdf_time_scheme_t) :: bdf
 
     type(vector_t) :: advection_coeffs
     type(vector_t) :: diffusion_coeffs
     integer :: ndiff = 0
     integer :: nadv = 0
     integer :: advection_time_order = 3
     integer :: diffusion_time_order
     type(c_ptr) :: advection_coeffs_d = c_null_ptr
     type(c_ptr) :: diffusion_coeffs_d = c_null_ptr
 
   contains
     procedure, pass(this) :: init => time_scheme_controller_init
     procedure, pass(this) :: free => time_scheme_controller_free
     procedure, pass(this) :: set_coeffs => &
          time_scheme_controller_set_coeffs
  end type time_scheme_controller_t
 
contains
 
  subroutine time_scheme_controller_init(this, torder)
    implicit none
    class(time_scheme_controller_t) :: this
    integer :: torder
 
    call this%free()
 
    this%diffusion_time_order = torder
 
    ! Force 1st order advection when diffusion is 1st order
    if (torder .eq. 1) then
       this%advection_time_order = 1
    end if
 
    call this%advection_coeffs%init(4)
    call this%diffusion_coeffs%init(4)
  end subroutine time_scheme_controller_init
 
  subroutine time_scheme_controller_free(this)
    implicit none
    class(time_scheme_controller_t) :: this
 
    call this%advection_coeffs%free()
    call this%diffusion_coeffs%free()
  end subroutine time_scheme_controller_free
 
  subroutine time_scheme_controller_set_coeffs(this, dt)
    implicit none
    class(time_scheme_controller_t) :: this
    real(kind=rp), intent(inout), dimension(10) :: dt
    real(kind=rp), dimension(4) :: adv_coeffs_old
    real(kind=rp), dimension(4) :: diff_coeffs_old
 
    associate( &
         nadv => this%nadv, &
         ndiff => this%ndiff, &
         adv_coeffs => this%advection_coeffs, &
         adv_coeffs_d => this%advection_coeffs_d, &
         diff_coeffs => this%diffusion_coeffs, &
         diff_coeffs_d => this%diffusion_coeffs_d)
 
      adv_coeffs_old = adv_coeffs%x
      diff_coeffs_old = diff_coeffs%x
 
      ! On unified memory backends the coefficient arrays may alias
      ! their device pointers, and the host-side rewrite below races
      ! in-flight device work reading them. Rewriting identical values
      ! is benign, so only drain the device when the coefficients can
      ! actually change (scheme order still ramping up, or a varying
      ! time step)
      if (neko_bcknd_device .eq. 1) then
         if ((ndiff .lt. this%diffusion_time_order) .or. &
              (nadv .lt. this%advection_time_order) .or. &
              (dt(1) .ne. dt(2))) then
            call device_sync()
         end if
      end if
 
      ! Increment the order of the scheme if below time_order
      ndiff = ndiff + 1
      ndiff = min(ndiff, this%diffusion_time_order)
      nadv = nadv + 1
      nadv = min(nadv, this%advection_time_order)
 
      call this%bdf%compute_coeffs(diff_coeffs%x, dt, ndiff)
 
      if (nadv .eq. 1) then
         ! Forward euler
         call this%ext%compute_coeffs(adv_coeffs%x, dt, nadv)
      else if (nadv .eq. 2) then
         if (ndiff .eq. 1) then
            ! 2nd order Adam-Bashforth, currently never used
            call this%ab%compute_coeffs(adv_coeffs%x, dt, nadv)
         else
            ! Linear extrapolation
            call this%ext%compute_coeffs(adv_coeffs%x, dt, nadv)
         end if
      else if (nadv .eq. 3) then
         if (ndiff .eq. 1) then
            ! 3rd order Adam-Bashforth, currently never used
            call this%ab%compute_coeffs(adv_coeffs%x, dt, nadv)
         else if (ndiff .eq. 2) then
            ! The modified EXT scheme
            call this%ext%compute_modified_coeffs(adv_coeffs%x, dt)
         else
            ! Quadratic extrapolation
            call this%ext%compute_coeffs(adv_coeffs%x, dt, nadv)
         end if
      end if
 
      if (neko_bcknd_device .eq. 1) then
         if (maxval(abs(adv_coeffs%x - adv_coeffs_old)) .gt. 1e-10_rp) then
            call adv_coeffs%copy_from(host_to_device, sync = .false.)
         end if
 
         if (maxval(abs(diff_coeffs%x - diff_coeffs_old)) .gt. 1e-10_rp) then
            call diff_coeffs%copy_from(host_to_device, sync = .false.)
         end if
      end if
    end associate
 
  end subroutine time_scheme_controller_set_coeffs
end module time_scheme_controller