develop.f90

! ============================================================================ !
! Supporting functions
! ============================================================================ !
module develop
  use case, only: case_t
  use global_interpolation, only: global_interpolation_t
  use json_utils, only: json_get_or_default
  use logger, only: neko_log, log_size
  use num_types, only: rp
  use tuple, only: tuple4_i4_t
  use utils, only: neko_error
  use hex, only: hex_t
  implicit none
 
  integer, public, parameter :: facet_type_interior = 0
  integer, public, parameter :: facet_type_outlet = 1
  integer, public, parameter :: facet_type_inlet = 2
  integer, public, parameter :: facet_type_wall = 3
  integer, public, parameter :: facet_type_sympln = 4
 
contains
 
  pure function cross(a, b) result(c)
    real(kind=rp), dimension(3), intent(in) :: a
    real(kind=rp), dimension(3), intent(in) :: b
    real(kind=rp), dimension(3) :: c
 
    c(1) = a(2) * b(3) - a(3) * b(2)
    c(2) = a(3) * b(1) - a(1) * b(3)
    c(3) = a(1) * b(2) - a(2) * b(1)
 
  end function cross
 
  pure function dot(a, b) result(d)
    real(kind=rp), dimension(3), intent(in) :: a(3)
    real(kind=rp), dimension(3), intent(in) :: b(3)
    real(kind=rp) :: d
 
    d = a(1) * b(1) + a(2) * b(2) + a(3) * b(3)
 
  end function dot
 
  pure function area(nv, vertices) result(a)
    integer, intent(in) :: nv
    real(kind=rp), dimension(3, nv), intent(in) :: vertices
    real(kind=rp) :: a
 
    ! Internal variables
    real(kind=rp), dimension(3) :: v1, e1, e2, tmp
    real(kind=rp) :: cp(3)
    integer :: v
 
    ! Initialize the origin to the first vertex, which reduce the
    ! computational cost by allowing us to skip the first and last cross
    ! products.
    v1 = vertices(:, 1)
    cp = 0.0
    do v = 2, nv - 1
       e1 = (vertices(:, v) - v1)
       e2 = (vertices(:, v + 1) - v1)
       tmp = cross(e1, e2)
       cp = cp + 0.5 * tmp
    end do
 
    a = sqrt(dot(cp, cp))
 
  end function area
 
  subroutine get_facets(C, out, facet_type)
    type(case_t), intent(in) :: C
    integer, dimension(:, :), allocatable, intent(out) :: out
    integer, intent(in), optional :: facet_type
 
    ! Local variables
    integer :: e, f, facet_id
    integer :: N_elem, N_facets
    integer :: facet_type_
 
    if (present(facet_type)) then
       facet_type_ = facet_type
    else
       facet_type_ = 1
    end if
 
    ! Total number of elements
    n_elem = c%msh%nelv
 
    ! Count number of output facets (Facet_type is 1)
    n_facets = 0
    do e = 1, n_elem
       do f = 1, 6
          if (c%msh%facet_type(f, e) .eq. facet_type_) then
             n_facets = n_facets + 1
          end if
       end do
    end do
 
    ! Allocate the outlet_elements array
    allocate (out(2, n_facets))
 
    ! Loop across the elements and store the outlet elements
    facet_id = 1
    do e = 1, n_elem
       do f = 1, 6
          if (c%msh%facet_type(f, e) .eq. facet_type_) then
             out(1, facet_id) = e
             out(2, facet_id) = f
             facet_id = facet_id + 1
          end if
       end do
    end do
 
  end subroutine get_facets
 
  subroutine get_facets_outlet(C, out)
    type(case_t), intent(in) :: C
    integer, dimension(:, :), allocatable, intent(out) :: out
 
    call get_facets(c, out, facet_type_outlet)
  end subroutine get_facets_outlet
 
  pure function sum_weighted(a, w) result(s)
    real(kind=rp), dimension(:), intent(in) :: a
    real(kind=rp), dimension(:), intent(in), optional :: w
    real(kind=rp) :: s
 
    if (present(w)) then
       s = sum(a * w)
    else
       s = sum(a)
    end if
 
  end function sum_weighted
 
  pure function average_weighted(a, w) result(m)
    real(kind=rp), dimension(:), intent(in) :: a
    real(kind=rp), dimension(:), intent(in), optional :: w
    real(kind=rp) :: m
 
    if (present(w)) then
       m = sum_weighted(a, w) / sum(w)
    else
       m = sum(a) / size(a)
    end if
 
  end function average_weighted
 
  subroutine get_facet_nodes(C, element_id, facet_id, nodes)
    type(case_t), intent(in) :: C
    integer, intent(in) :: element_id
    integer, intent(in) :: facet_id
    real(kind=rp), dimension(:, :), allocatable, intent(out) :: nodes
 
    ! Local variables
    integer :: n, v
    integer :: N_nodes
    type(tuple4_i4_t) :: t_hex
 
    ! Determine the number of nodes in the facet
    n_nodes = 0
    select type (ele => c%msh%elements(element_id)%e)
    type is (hex_t)
       n_nodes = 4
    end select
 
    if (n_nodes .eq. 0) then
       call neko_error("Error: the facet shape is not supported.")
    end if
 
    ! Allocate the nodes array
    if (allocated(nodes) .and. size(nodes, 2) .eq. n_nodes) then
       nodes = 0.0_rp
    else if (.not. allocated(nodes)) then
       allocate (nodes(3, n_nodes))
    else
       call neko_error("Error: the nodes array is not the correct size.")
    end if
 
    ! Get the nodes
    select type (ele => c%msh%elements(element_id)%e)
    type is (hex_t)
       call ele%facet_order(t_hex, facet_id)
       do n = 1, n_nodes
          v = t_hex%x(n)
          nodes(:, n) = c%msh%points(v)%x
       end do
    end select
 
  end subroutine get_facet_nodes
 
  subroutine estimate_temperature(neko_case)
    type(case_t), intent(inout) :: neko_case
    type(global_interpolation_t) :: interpolator
 
    real(kind=rp) :: target_temperature
    real(kind=rp) :: temperature_mean
    real(kind=rp), dimension(:), allocatable :: temperature_local
 
    integer, allocatable :: facet_list(:, :)
 
    real(kind=rp), dimension(:, :), allocatable :: facet_nodes
    real(kind=rp), dimension(:, :), allocatable :: facet_centers
    real(kind=rp), dimension(:), allocatable :: facet_area
 
    integer :: element_id, facet_id
    integer :: N_facets
    integer :: f, n
 
    character(len=LOG_SIZE) :: log_buf
 
    ! Read the case file for options
    call json_get_or_default(neko_case%params, 'topopt.target_temperature', &
         target_temperature, 0.5_rp)
 
    ! Initialize the global interpolation
    call interpolator%init(neko_case%scalar%dm_xh)
 
    ! Get the list of outlet facets
    call get_facets(neko_case, facet_list)
 
    ! Allocate the facet_list array
    n_facets = size(facet_list, 2)
    allocate (facet_centers(3, n_facets))
    allocate (facet_area(n_facets))
    allocate (temperature_local(n_facets))
 
    ! Loop across the outlet elements and print the scalar field
    do f = 1, n_facets
       element_id = facet_list(1, f)
       facet_id = facet_list(2, f)
 
       ! Get the nodes for the facet
       call get_facet_nodes(neko_case, element_id, facet_id, facet_nodes)
 
       ! Compute the center of the facet
       facet_centers(:, f) = 0.0
       do n = 1, size(facet_nodes, 2)
          facet_centers(:, f) = facet_centers(:, f) + facet_nodes(:, n)
       end do
       facet_centers(:, f) = facet_centers(:, f) / size(facet_nodes, 2)
 
       ! Compute the area of the facet
       facet_area(f) = area(size(facet_nodes, 2), facet_nodes)
    end do
 
    ! Find the outlet temperature at the supplied list of points
    call interpolator%find_points_xyz(facet_centers, n_facets)
    call interpolator%evaluate(temperature_local, neko_case%scalar%s%x)
 
    temperature_local = temperature_local - target_temperature
    temperature_mean = average_weighted(temperature_local, facet_area)
 
    write (log_buf, '(a,f15.7)') &
         "Outlet area-weighted average temperature deviation: ", &
         temperature_mean
    call neko_log%message(log_buf)
 
    call interpolator%free()
 
  end subroutine estimate_temperature
 
end module develop