wall_model.f90

Go to the documentation of this file.

! Copyright (c) 2024, 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 wall_model
  use num_types, only : rp
  use field, only : field_t
  use json_module, only : json_file
  use registry, only : neko_registry
  use dofmap, only : dofmap_t
  use coefs, only : coef_t
  use neko_config, only : neko_bcknd_device
  use device, only : device_memcpy, host_to_device
  use vector, only : vector_t
  use utils, only : neko_error, nonlinear_index
  use math, only : glmin, glmax
  use comm, only : pe_rank
  use logger, only : neko_log, neko_log_debug, log_size
  use file, only : file_t
  use, intrinsic :: iso_c_binding, only : c_ptr, c_null_ptr, c_associated
  use device, only : device_map, device_free, device_get_ptr
  use wall_model_device, only : wall_model_compute_mag_field_device
  use json_utils, only : json_get
  implicit none
  private
 
  type, abstract, public :: wall_model_t
     type(coef_t), pointer :: coef => null()
     type(dofmap_t), pointer :: dof => null()
     type(field_t), pointer :: mu => null()
     type(field_t), pointer :: rho => null()
     character(len=:), allocatable :: scheme_name
     integer, pointer :: msk(:) => null()
     type(c_ptr) :: msk_d = c_null_ptr
     integer, pointer :: facet(:) => null()
     type(vector_t) :: tau_x
     type(vector_t) :: tau_y
     type(vector_t) :: tau_z
     type(vector_t) :: n_x
     type(vector_t) :: n_y
     type(vector_t) :: n_z
     integer, allocatable :: ind_r(:)
     type(c_ptr) :: ind_r_d = c_null_ptr
     integer, allocatable :: ind_s(:)
     type(c_ptr) :: ind_s_d = c_null_ptr
     integer, allocatable :: ind_t(:)
     type(c_ptr) :: ind_t_d = c_null_ptr
     integer, allocatable :: ind_e(:)
     type(c_ptr) :: ind_e_d = c_null_ptr
     type(vector_t) :: h
     integer :: h_index = 0
     integer :: n_nodes = 0
     type(field_t), pointer :: tau_field => null()
   contains
     procedure, pass(this) :: init_base => wall_model_init_base
     procedure, pass(this) :: partial_init_base => wall_model_partial_init_base
     procedure, pass(this) :: finalize_base => wall_model_finalize_base
     procedure, pass(this) :: free_base => wall_model_free_base
     procedure, pass(this) :: compute_mag_field => wall_model_compute_mag_field
     procedure(wall_model_init), pass(this), deferred :: init
     procedure(wall_model_partial_init), pass(this), deferred :: partial_init
     procedure(wall_model_finalize), pass(this), deferred :: finalize
     procedure(wall_model_free), pass(this), deferred :: free
     procedure(wall_model_compute), pass(this), deferred :: compute
     procedure, pass(this) :: find_points => wall_model_find_points
  end type wall_model_t
 
  abstract interface
 
     subroutine wall_model_compute(this, t, tstep)
       import wall_model_t, rp
       class(wall_model_t), intent(inout) :: this
       real(kind=rp), intent(in) :: t
       integer, intent(in) :: tstep
     end subroutine wall_model_compute
  end interface
 
  abstract interface
 
     subroutine wall_model_init(this, scheme_name, coef, msk, facet, &
          h_index, json)
       import wall_model_t, json_file, dofmap_t, coef_t, rp
       class(wall_model_t), intent(inout) :: this
       character(len=*), intent(in) :: scheme_name
       type(coef_t), intent(in) :: coef
       integer, intent(in) :: msk(:)
       integer, intent(in) :: facet(:)
       integer, intent(in) :: h_index
       type(json_file), intent(inout) :: json
     end subroutine wall_model_init
  end interface
 
  abstract interface
 
     subroutine wall_model_partial_init(this, coef, json)
       import wall_model_t, json_file, dofmap_t, coef_t, rp
       class(wall_model_t), intent(inout) :: this
       type(coef_t), intent(in) :: coef
       type(json_file), intent(inout) :: json
     end subroutine wall_model_partial_init
  end interface
 
  abstract interface
 
     subroutine wall_model_finalize(this, msk, facet)
       import wall_model_t
       class(wall_model_t), intent(inout) :: this
       integer, intent(in) :: msk(:)
       integer, intent(in) :: facet(:)
     end subroutine wall_model_finalize
  end interface
 
  abstract interface
 
     subroutine wall_model_free(this)
       import wall_model_t
       class(wall_model_t), intent(inout) :: this
     end subroutine wall_model_free
  end interface
 
  interface
 
     module subroutine wall_model_factory(object, scheme_name, coef, msk, &
          facet, json)
       class(wall_model_t), allocatable, intent(inout) :: object
       character(len=*), intent(in) :: scheme_name
       type(coef_t), intent(in) :: coef
       integer, intent(in) :: msk(:)
       integer, intent(in) :: facet(:)
       type(json_file), intent(inout) :: json
     end subroutine wall_model_factory
  end interface
 
  interface
 
     module subroutine wall_model_allocator(object, type_name)
       class(wall_model_t), allocatable, intent(inout) :: object
       character(len=:), allocatable, intent(in) :: type_name
     end subroutine wall_model_allocator
  end interface
 
  !
  ! Machinery for injecting user-defined types
  !
 
  abstract interface
     subroutine wall_model_allocate(obj)
       import wall_model_t
       class(wall_model_t), allocatable, intent(inout) :: obj
     end subroutine wall_model_allocate
  end interface
 
  interface
 
     module subroutine register_wall_model(type_name, allocator)
       character(len=*), intent(in) :: type_name
       procedure(wall_model_allocate), pointer, intent(in) :: allocator
     end subroutine register_wall_model
  end interface
 
  ! A name-allocator pair for user-defined types. A helper type to define a
  ! registry of custom allocators.
  type allocator_entry
     character(len=20) :: type_name
     procedure(wall_model_allocate), pointer, nopass :: allocator
  end type allocator_entry
 
  type(allocator_entry), allocatable :: wall_model_registry(:)
 
  integer :: wall_model_registry_size = 0
 
  public :: wall_model_factory, wall_model_allocator, register_wall_model, &
       wall_model_allocate
 
contains
  subroutine wall_model_init_base(this, scheme_name, coef, msk, facet, index)
    class(wall_model_t), intent(inout) :: this
    type(coef_t), target, intent(in) :: coef
    integer, target, intent(in) :: msk(0:)
    integer, target, intent(in) :: facet(0:)
    character(len=*) :: scheme_name
    integer, intent(in) :: index
 
    call this%free_base
 
    this%coef => coef
    this%dof => coef%dof
    this%h_index = index
    this%scheme_name = trim(scheme_name)
    this%mu => neko_registry%get_field_by_name(this%scheme_name // "_mu")
    this%rho => neko_registry%get_field_by_name(this%scheme_name // &
         "_rho")
 
    call neko_registry%add_field(this%dof, "tau", &
         ignore_existing = .true.)
    this%tau_field => neko_registry%get_field("tau")
 
    call this%finalize_base(msk, facet)
  end subroutine wall_model_init_base
 
  subroutine wall_model_partial_init_base(this, coef, json)
    class(wall_model_t), intent(inout) :: this
    type(coef_t), target, intent(in) :: coef
    type(json_file), intent(inout) :: json
 
    call this%free_base()
 
    this%coef => coef
    this%dof => coef%dof
    call json_get(json, "h_index", this%h_index)
    call json_get(json, "scheme_name", this%scheme_name)
 
    this%mu => neko_registry%get_field_by_name(this%scheme_name // "_mu")
    this%rho => neko_registry%get_field_by_name(this%scheme_name // &
         "_rho")
 
    call neko_registry%add_field(this%dof, "tau", &
         ignore_existing = .true.)
    this%tau_field => neko_registry%get_field("tau")
  end subroutine wall_model_partial_init_base
 
  subroutine wall_model_finalize_base(this, msk, facet)
    class(wall_model_t), intent(inout) :: this
    integer, target, intent(in) :: msk(0:)
    integer, target, intent(in) :: facet(:)
 
    this%msk(0:msk(0)) => msk
    if (neko_bcknd_device .eq. 1) this%msk_d = device_get_ptr(msk)
    this%facet(0:msk(0)) => facet
 
    call this%tau_x%init(this%msk(0))
    call this%tau_y%init(this%msk(0))
    call this%tau_z%init(this%msk(0))
 
    allocate(this%ind_r(this%msk(0)))
    allocate(this%ind_s(this%msk(0)))
    allocate(this%ind_t(this%msk(0)))
    allocate(this%ind_e(this%msk(0)))
 
    call this%h%init(this%msk(0))
    call this%n_x%init(this%msk(0))
    call this%n_y%init(this%msk(0))
    call this%n_z%init(this%msk(0))
 
    call this%find_points()
 
    ! Initialize pointers for device
    if (neko_bcknd_device .eq. 1) then
       call device_map(this%ind_r, this%ind_r_d, this%n_nodes)
       call device_map(this%ind_s, this%ind_s_d, this%n_nodes)
       call device_map(this%ind_t, this%ind_t_d, this%n_nodes)
       call device_map(this%ind_e, this%ind_e_d, this%n_nodes)
       call device_memcpy(this%ind_r, this%ind_r_d, this%n_nodes, &
            host_to_device, sync = .false.)
       call device_memcpy(this%ind_s, this%ind_s_d, this%n_nodes, &
            host_to_device, sync = .false.)
       call device_memcpy(this%ind_t, this%ind_t_d, this%n_nodes, &
            host_to_device, sync = .false.)
       call device_memcpy(this%ind_e, this%ind_e_d, this%n_nodes, &
            host_to_device, sync = .false.)
    end if
 
  end subroutine wall_model_finalize_base
 
  subroutine wall_model_free_base(this)
    class(wall_model_t), intent(inout) :: this
 
    nullify(this%coef)
    nullify(this%msk)
    nullify(this%facet)
    nullify(this%tau_field)
    nullify(this%mu)
    nullify(this%rho)
 
    call this%tau_x%free()
    call this%tau_y%free()
    call this%tau_z%free()
 
    if (allocated(this%ind_r)) then
       deallocate(this%ind_r)
    end if
    if (allocated(this%ind_s)) then
       deallocate(this%ind_s)
    end if
    if (allocated(this%ind_t)) then
       deallocate(this%ind_t)
    end if
    if (allocated(this%ind_e)) then
       deallocate(this%ind_e)
    end if
 
    if (c_associated(this%msk_d)) then
       call device_free(this%msk_d)
    end if
    if (c_associated(this%ind_r_d)) then
       call device_free(this%ind_r_d)
    end if
    if (c_associated(this%ind_s_d)) then
       call device_free(this%ind_s_d)
    end if
    if (c_associated(this%ind_t_d)) then
       call device_free(this%ind_t_d)
    end if
    if (c_associated(this%ind_e_d)) then
       call device_free(this%ind_e_d)
    end if
 
    if (allocated(this%scheme_name)) then
       deallocate(this%scheme_name)
    end if
 
 
    call this%h%free()
    call this%n_x%free()
    call this%n_y%free()
    call this%n_z%free()
 
    nullify(this%dof)
  end subroutine wall_model_free_base
 
  subroutine wall_model_find_points(this)
    class(wall_model_t), intent(inout) :: this
    integer :: n_nodes, fid, idx(4), i, linear
    real(kind=rp) :: normal(3), p(3), x, y, z, xw, yw, zw, magp
    real(kind=rp) :: hmin, hmax
    type(field_t), pointer :: h_field
    type(file_t) :: h_file
    character(len=LOG_SIZE), allocatable :: log_msg
 
    n_nodes = this%msk(0)
    this%n_nodes = n_nodes
 
    call neko_registry%add_field(this%coef%dof, "sampling_height", &
         ignore_existing=.true.)
 
    h_field => neko_registry%get_field_by_name("sampling_height")
 
    do i = 1, n_nodes
       linear = this%msk(i)
       fid = this%facet(i)
       idx = nonlinear_index(linear, this%coef%Xh%lx, this%coef%Xh%ly,&
            this%coef%Xh%lz)
       normal = this%coef%get_normal(idx(1), idx(2), idx(3), idx(4), fid)
 
       this%n_x%x(i) = normal(1)
       this%n_y%x(i) = normal(2)
       this%n_z%x(i) = normal(3)
 
       ! inward normal
       normal = -normal
 
       select case (fid)
       case (1)
          this%ind_r(i) = idx(1) + this%h_index
          this%ind_s(i) = idx(2)
          this%ind_t(i) = idx(3)
       case (2)
          this%ind_r(i) = idx(1) - this%h_index
          this%ind_s(i) = idx(2)
          this%ind_t(i) = idx(3)
       case (3)
          this%ind_r(i) = idx(1)
          this%ind_s(i) = idx(2) + this%h_index
          this%ind_t(i) = idx(3)
       case (4)
          this%ind_r(i) = idx(1)
          this%ind_s(i) = idx(2) - this%h_index
          this%ind_t(i) = idx(3)
       case (5)
          this%ind_r(i) = idx(1)
          this%ind_s(i) = idx(2)
          this%ind_t(i) = idx(3) + this%h_index
       case (6)
          this%ind_r(i) = idx(1)
          this%ind_s(i) = idx(2)
          this%ind_t(i) = idx(3) - this%h_index
       case default
          call neko_error("The face index is not correct ")
       end select
       this%ind_e(i) = idx(4)
 
       ! Location of the wall node
       xw = this%dof%x(idx(1), idx(2), idx(3), idx(4))
       yw = this%dof%y(idx(1), idx(2), idx(3), idx(4))
       zw = this%dof%z(idx(1), idx(2), idx(3), idx(4))
 
       ! Location of the sampling point
       x = this%dof%x(this%ind_r(i), this%ind_s(i), this%ind_t(i), &
            this%ind_e(i))
       y = this%dof%y(this%ind_r(i), this%ind_s(i), this%ind_t(i), &
            this%ind_e(i))
       z = this%dof%z(this%ind_r(i), this%ind_s(i), this%ind_t(i), &
            this%ind_e(i))
 
       ! Vector from the sampling point to the wall
       p(1) = x - xw
       p(2) = y - yw
       p(3) = z - zw
 
       ! Total distance to the sampling point
       magp = sqrt(p(1)**2 + p(2)**2 + p(3)**2)
 
       ! Project on the normal direction to get h
       this%h%x(i) = p(1)*normal(1) + p(2)*normal(2) + p(3)*normal(3)
 
       h_field%x(linear,1,1,1) = this%h%x(i)
 
       ! Look at how much the total distance distance from the normal and warn
       ! if significant
       if ((this%h%x(i) - magp) / magp > 0.1) then
          write(log_msg,*) "Significant misalignment between wall normal and"
          call neko_log%message(log_msg, neko_log_debug)
          write(log_msg,*) "sampling point direction at wall node", xw, yw, zw
          call neko_log%message(log_msg, neko_log_debug)
       end if
    end do
 
!    hmin = glmin(this%h%x, n_nodes)
!    hmax = glmax(this%h%x, n_nodes)
!    if (pe_rank .eq. 0) then
!       write(*, "(A, F10.4, F10.4)") "   h min / max:", hmin, hmax
!    end if
 
    if (neko_bcknd_device .eq. 1) then
       call device_memcpy(this%h%x, this%h%x_d, n_nodes, host_to_device,&
            sync = .false.)
       call device_memcpy(this%n_x%x, this%n_x%x_d, n_nodes, host_to_device, &
            sync = .false.)
       call device_memcpy(this%n_y%x, this%n_y%x_d, n_nodes, host_to_device, &
            sync = .false.)
       call device_memcpy(this%n_z%x, this%n_z%x_d, n_nodes, host_to_device, &
            sync = .true.)
    end if
 
    ! Each wall_model bc will do a write unfortunately... But very helpful
    ! for setup debugging.
    call h_file%init("sampling_height.fld")
    call h_file%write(h_field)
  end subroutine wall_model_find_points
 
  subroutine wall_model_compute_mag_field(this)
    class(wall_model_t), intent(inout) :: this
    integer :: i, m
    real(kind=rp) :: magtau
 
    m = this%msk(0)
    if (m > 0) then
       if (neko_bcknd_device .eq. 1) then
          call wall_model_compute_mag_field_device(this%tau_x%x_d, &
               this%tau_y%x_d, &
               this%tau_z%x_d, &
               this%tau_field%x_d, &
               this%msk_d, m)
       else
          do i = 1, m
             magtau = sqrt(this%tau_x%x(i)**2 + &
                  this%tau_y%x(i)**2 + &
                  this%tau_z%x(i)**2)
             this%tau_field%x(this%msk(i),1,1,1) = magtau
          end do
       end if
    end if
 
  end subroutine wall_model_compute_mag_field
 
end module wall_model