wall_model.f90

Go to the documentation of this file.

! Copyright (c) 2024-2026, 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 wall_sampler, only : wall_sampler_t
  use wall_sampler_fctry, only : wall_sampler_factory
  use user_intf, only : user_t
  use, intrinsic :: iso_c_binding, only : c_ptr, c_null_ptr
  use device, only : device_get_ptr
  use wall_model_device, only : wall_model_compute_mag_field_device
  use json_utils, only : json_get, json_get_or_lookup
  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
     class(wall_sampler_t), allocatable :: sampler
     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) :: setup_geometry => wall_model_setup_geometry
     procedure, pass(this) :: validate_single_sample => &
          wall_model_validate_single_sample
  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, &
          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(:)
       type(json_file), intent(inout) :: json
     end subroutine wall_model_init
  end interface
 
  abstract interface
 
     subroutine wall_model_partial_init(this, coef, scheme_name, json)
       import wall_model_t, json_file, dofmap_t, coef_t, rp, user_t
       class(wall_model_t), intent(inout) :: this
       type(coef_t), intent(in) :: coef
       character(len=*), intent(in) :: scheme_name
       type(json_file), intent(inout) :: json
     end subroutine wall_model_partial_init
  end interface
 
  abstract interface
 
     subroutine wall_model_finalize(this, msk, facet, bc_name, user)
       import wall_model_t, user_t
       class(wall_model_t), intent(inout) :: this
       integer, intent(in) :: msk(:)
       integer, intent(in) :: facet(:)
       character(len=*), optional, intent(in) :: bc_name
       type(user_t), target, optional, intent(in) :: user
     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, sampler)
    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
    class(wall_sampler_t), allocatable, intent(inout) :: sampler
 
    call this%free_base
 
    this%coef => coef
    this%dof => coef%dof
    if (.not. allocated(sampler)) then
       call neko_error('Wall model sampler has not been initialized')
    end if
    call move_alloc(sampler, this%sampler)
    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, scheme_name, json)
    class(wall_model_t), intent(inout) :: this
    type(coef_t), target, intent(in) :: coef
    character(len=*), intent(in) :: scheme_name
    type(json_file), intent(inout) :: json
 
    call this%free_base()
 
    this%coef => coef
    this%dof => coef%dof
    call wall_sampler_factory(this%sampler, json)
    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")
  end subroutine wall_model_partial_init_base
 
  subroutine wall_model_finalize_base(this, msk, facet, bc_name, user)
    class(wall_model_t), intent(inout) :: this
    integer, target, intent(in) :: msk(0:)
    integer, target, intent(in) :: facet(:)
    character(len=*), optional, intent(in) :: bc_name
    type(user_t), target, optional, intent(in) :: user
 
    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))
 
    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%setup_geometry()
    if (.not. allocated(this%sampler)) then
       call neko_error('Wall model sampler has not been initialized')
    end if
    if (present(bc_name)) then
       if (present(user)) then
          call this%sampler%finalize(this%coef, this%msk, this%facet, &
               this%n_x, this%n_y, this%n_z, bc_name, user)
       else
          call this%sampler%finalize(this%coef, this%msk, this%facet, &
               this%n_x, this%n_y, this%n_z, bc_name)
       end if
    else
       if (this%sampler%user_values) then
          call neko_error('Wall model user sampler has no boundary context')
       end if
       call this%sampler%finalize(this%coef, this%msk, this%facet, &
            this%n_x, this%n_y, this%n_z)
    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)
    this%msk_d = c_null_ptr
 
    call this%tau_x%free()
    call this%tau_y%free()
    call this%tau_z%free()
 
    if (allocated(this%sampler)) then
       call this%sampler%free()
       deallocate(this%sampler)
    end if
 
    if (allocated(this%scheme_name)) then
       deallocate(this%scheme_name)
    end if
 
    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_setup_geometry(this)
    class(wall_model_t), intent(inout) :: this
    integer :: n_nodes, fid, idx(4), i, linear
    real(kind=rp) :: normal(3)
 
    call this%coef%require_facets('wall_model')
 
    n_nodes = this%msk(0)
    this%n_nodes = n_nodes
 
    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)
 
    end do
 
    if (neko_bcknd_device .eq. 1) then
       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
 
  end subroutine wall_model_setup_geometry
 
  subroutine wall_model_validate_single_sample(this)
    class(wall_model_t), intent(inout) :: this
 
    if (this%sampler%n_samples /= 1) then
       call neko_error('This wall model requires exactly one sample per node')
    end if
  end subroutine wall_model_validate_single_sample
 
  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