map_1d.f90

Go to the documentation of this file.

! Copyright (c) 2023-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 map_1d
  use neko_config, only : neko_bcknd_device
  use num_types, only : rp
  use space, only : space_t
  use dofmap, only : dofmap_t
  use gather_scatter, only : gs_op_add
  use mesh, only : mesh_t
  use device, only : device_memcpy, device_map, device_unmap, &
       host_to_device, device_to_host
  use comm, only : pe_size, pe_rank, neko_comm, mpi_real_precision
  use coefs, only : coef_t
  use field_list, only : field_list_t
  use matrix, only : matrix_t
  use vector, only : vector_ptr_t
  use utils, only : neko_error, neko_warning
  use math, only : glmax, glmin, glimax, glimin, cmult, add2s1, col2
  use mpi_f08, only : mpi_allreduce, mpi_sum, mpi_barrier, mpi_in_place
  use, intrinsic :: iso_c_binding
  implicit none
  private
 
  type, public :: map_1d_t
     integer, allocatable :: dir_el(:)
     integer, allocatable :: el_lvl(:)
     integer, allocatable :: pt_lvl(:, :, :, :)
     integer :: n_el_lvls
     integer :: n_gll_lvls
     type(dofmap_t), pointer :: dof => null()
     type(coef_t), pointer :: coef => null()
     type(mesh_t), pointer :: msh => null()
     integer :: dir
     real(kind=rp) :: tol = 1e-7
     real(kind=rp), allocatable :: volume_per_gll_lvl(:)
   contains
     procedure, pass(this) :: init_int => map_1d_init
     procedure, pass(this) :: init_char => map_1d_init_char
     generic :: init => init_int, init_char
     procedure, pass(this) :: free => map_1d_free
     procedure, pass(this) :: average_planes_fld_lst => &
          map_1d_average_field_list
     procedure, pass(this) :: average_planes_vec_ptr => &
          map_1d_average_vector_ptr
     generic :: average_planes => average_planes_fld_lst, average_planes_vec_ptr
  end type map_1d_t
 
 
contains
 
  subroutine map_1d_init(this, coef, dir, tol)
    class(map_1d_t) :: this
    type(coef_t), intent(inout), target :: coef
    integer, intent(in) :: dir
    real(kind=rp), intent(in) :: tol
    integer :: nelv, lx, n, i, e, lvl, ierr
    real(kind=rp), contiguous, pointer :: line(:, :, :, :)
    real(kind=rp), allocatable :: min_vals(:, :, :, :)
    real(kind=rp), allocatable :: min_temp(:, :, :, :)
    type(c_ptr) :: min_vals_d = c_null_ptr
    real(kind=rp) :: el_dim(3, 3), glb_min, glb_max, el_min, abs_tol
 
    call this%free()
 
    if (neko_bcknd_device .eq. 1) then
       if (pe_rank .eq. 0) then
          call neko_warning('map_1d does not copy indices to device, ' // &
               ' but ok if used on cpu and for io')
       end if
    end if
 
    this%dir = dir
    this%tol = tol
    this%dof => coef%dof
    this%coef => coef
    this%msh => coef%msh
    nelv = this%msh%nelv
    lx = this%dof%Xh%lx
    n = this%dof%size()
 
    if (dir .eq. 1) then
       line => this%dof%x%x
    else if (dir .eq. 2) then
       line => this%dof%y%x
    else if (dir .eq. 3) then
       line => this%dof%z%x
    else
       call neko_error('Invalid dir for geopmetric comm')
    end if
    allocate(this%dir_el(nelv))
    allocate(this%el_lvl(nelv))
    allocate(this%pt_lvl(lx, lx, lx, nelv))
    allocate(min_vals(lx, lx, lx, nelv))
    allocate(min_temp(lx, lx, lx, nelv))
    call mpi_barrier(neko_comm)
    if (neko_bcknd_device .eq. 1) then
 
       call device_map(min_vals, min_vals_d, n)
    end if
    call mpi_barrier(neko_comm)
 
    do i = 1, nelv
       ! Store which direction r,s,t corresponds to specified direction, x,y,z
       ! we assume elements are stacked on each other...
       ! Check which one of the normalized vectors are closest to dir
       ! If we want to incorporate other directions, we should look here
 
       ! This follows the point ordering defined in hex.f90
       el_dim(1, :) = abs(this%msh%elements(i)%e%pts(1)%p%x - &
            this%msh%elements(i)%e%pts(2)%p%x)
       el_dim(1, :) = el_dim(1, :)/norm2(el_dim(1, :))
       el_dim(2, :) = abs(this%msh%elements(i)%e%pts(1)%p%x - &
            this%msh%elements(i)%e%pts(3)%p%x)
       el_dim(2, :) = el_dim(2, :)/norm2(el_dim(2, :))
       el_dim(3, :) = abs(this%msh%elements(i)%e%pts(1)%p%x - &
            this%msh%elements(i)%e%pts(5)%p%x)
       el_dim(3, :) = el_dim(3, :)/norm2(el_dim(3, :))
       ! Checks which directions in rst the xyz corresponds to
       ! 1 corresponds to r, 2 to s, 3 to t and are stored in dir_el
       this%dir_el(i) = maxloc(el_dim(:, this%dir), dim = 1)
    end do
 
    glb_min = glmin(line, n)
    glb_max = glmax(line, n)
 
    ! The propagation below moves the element minima through the mesh with
    ! gather-scatter sums scaled by coef%mult. That arithmetic is not exact:
    ! coef%mult is 1/3 at some vertices of unstructured meshes. With the
    ! default floating point model of ifx even 1/1 is rounded. A propagated
    ! minimum can therefore differ from glb_min by a few roundoff errors of
    ! the coordinate magnitude. A tolerance relative to glb_min would demand
    ! exact equality when glb_min is zero, so compare with an absolute
    ! tolerance instead. The user tolerance is relative to the extent of the
    ! coordinates, so that it does not depend on where the origin is, and
    ! the roundoff floor is relative to their magnitude, which is what the
    ! propagation error scales with.
    abs_tol = max(this%tol * (glb_max - glb_min), &
         8.0_rp * epsilon(1.0_rp) * max(abs(glb_min), abs(glb_max)))
 
    i = 1
    this%el_lvl = -1
    ! Check what the minimum value in each element and put in min_vals
    do e = 1, nelv
       el_min = minval(line(:, :, :, e))
       ! Check if this element is on the bottom,
       ! in this case assign el_lvl = i = 1
       if (abs(el_min - glb_min) .le. abs_tol) then
          el_min = glb_min
          if (this%el_lvl(e) .eq. -1) this%el_lvl(e) = i
       end if
       min_vals(:, :, :, e) = el_min
    end do
 
    ! While loop where at each iteration the global maximum value
    ! propagates down one level.
    ! When the minimum value has propagated to the highest level this stops.
    ! Only works when the bottom plate of the domain is flat.
    do while (abs(glmax(min_vals, n) - glb_min) .gt. abs_tol)
 
       ! This is the assigned level
       i = i + 1
       ! There cannot be more levels than elements, so a longer propagation
       ! means the minimum never reaches some elements within abs_tol. Leave
       ! the loop and report that below.
       if (i .gt. this%msh%glb_nelv) exit
 
       do e = 1, nelv
          !Sets the value at the bottom of each element to glb_max
          if (this%dir_el(e) .eq. 1) then
             if (line(1, 1, 1, e) .gt. line(lx, 1, 1, e)) then
                min_vals(lx, :, :, e) = glb_max
             else
                min_vals(1, :, :, e) = glb_max
             end if
          end if
          if (this%dir_el(e) .eq. 2) then
             if (line(1, 1, 1, e) .gt. line(1, lx, 1, e)) then
                min_vals(:, lx, :, e) = glb_max
             else
                min_vals(:, 1, :, e) = glb_max
             end if
          end if
          if (this%dir_el(e) .eq. 3) then
             if (line(1, 1, 1, e) .gt. line(1, 1, lx, e)) then
                min_vals(:, :, lx, e) = glb_max
             else
                min_vals(:, :, 1, e) = glb_max
             end if
          end if
       end do
 
       !Make sketchy min as GS_OP_MIN is not supported with device mpi
       min_temp = min_vals
 
       if (neko_bcknd_device .eq. 1) then
          call device_memcpy(min_vals, min_vals_d, n, host_to_device, &
               sync = .false.)
       end if
 
       !Propagates the minimum value along the element boundary.
       call coef%gs_h%op(min_vals, n, gs_op_add)
 
       if (neko_bcknd_device .eq. 1) then
          call device_memcpy(min_vals, min_vals_d, n, device_to_host, &
               sync = .true.)
       end if
 
       !Obtain average along boundary
       call col2(min_vals, coef%mult, n)
       call cmult(min_temp, -1.0_rp, n)
       call add2s1(min_vals, min_temp, 2.0_rp, n)
 
       !Checks the new minimum value on each element
       !Assign this value to all points in this element in min_val
       !If the element has not already been assigned a level,
       !and it has obtained the minval, set el_lvl = i
       do e = 1, nelv
          el_min = minval(min_vals(:, :, :, e))
          if (abs(el_min - glb_min) .le. abs_tol) then
             ! Store exactly glb_min so that the roundoff of this propagation
             ! step does not accumulate into the next level.
             el_min = glb_min
             if (this%el_lvl(e) .eq. -1) this%el_lvl(e) = i
          end if
          min_vals(:, :, :, e) = el_min
       end do
    end do
 
    if (allocated(min_vals)) then
       if (c_associated(min_vals_d)) call device_unmap(min_vals, min_vals_d)
       deallocate(min_vals)
    end if
 
    if (allocated(min_temp)) deallocate(min_temp)
 
    ! An unassigned element would give negative point levels, which index
    ! volume_per_gll_lvl and the averaged output out of bounds. Return after
    ! the error, since a unit test with errors redirected to pFUnit
    ! exceptions continues past neko_error.
    if (glimin(this%el_lvl, nelv) .lt. 1) then
       call neko_error('map_1d: some elements were not assigned a level, ' // &
            'check that the mesh is layered in the requested direction')
       return
    end if
    this%n_el_lvls = glimax(this%el_lvl, nelv)
    this%n_gll_lvls = this%n_el_lvls*lx
 
    !Numbers the points in each element based on the element level
    !and its orientation
    do e = 1, nelv
       do i = 1, lx
          lvl = lx * (this%el_lvl(e) - 1) + i
          if (this%dir_el(e) .eq. 1) then
             if (line(1, 1, 1, e) .gt. line(lx, 1, 1, e)) then
                this%pt_lvl(lx-i+1, :, :, e) = lvl
             else
                this%pt_lvl(i, :, :, e) = lvl
             end if
          end if
          if (this%dir_el(e) .eq. 2) then
             if (line(1, 1, 1, e) .gt. line(1, lx, 1, e)) then
                this%pt_lvl(:, lx-i+1, :, e) = lvl
             else
                this%pt_lvl(:, i, :, e) = lvl
             end if
          end if
          if (this%dir_el(e) .eq. 3) then
             if (line(1, 1, 1, e) .gt. line(1, 1, lx, e)) then
                this%pt_lvl(:, :, lx-i+1, e) = lvl
             else
                this%pt_lvl(:, :, i, e) = lvl
             end if
          end if
       end do
    end do
 
    allocate(this%volume_per_gll_lvl(this%n_gll_lvls))
 
    this%volume_per_gll_lvl = 0.0_rp
 
    do i = 1, n
       this%volume_per_gll_lvl(this%pt_lvl(i, 1, 1, 1)) = &
            this%volume_per_gll_lvl(this%pt_lvl(i, 1, 1, 1)) + &
            coef%B(i, 1, 1, 1)
    end do
 
    call mpi_allreduce(mpi_in_place, this%volume_per_gll_lvl, this%n_gll_lvls, &
         mpi_real_precision, mpi_sum, neko_comm, ierr)
 
  end subroutine map_1d_init
 
  subroutine map_1d_init_char(this, coef, dir, tol)
    class(map_1d_t) :: this
    type(coef_t), intent(inout), target :: coef
    character(len=*), intent(in) :: dir
    real(kind=rp), intent(in) :: tol
    integer :: idir
 
    if (trim(dir) .eq. 'yz' .or. trim(dir) .eq. 'zy') then
       idir = 1
    else if (trim(dir) .eq. 'xz' .or. trim(dir) .eq. 'zx') then
       idir = 2
    else if (trim(dir) .eq. 'xy' .or. trim(dir) .eq. 'yx') then
       idir = 3
    else
       call neko_error('homogenous direction not supported')
    end if
 
    call this%init(coef, idir, tol)
 
  end subroutine map_1d_init_char
 
  subroutine map_1d_free(this)
    class(map_1d_t) :: this
 
    if (allocated(this%dir_el)) deallocate(this%dir_el)
    if (allocated(this%el_lvl)) deallocate(this%el_lvl)
    if (allocated(this%pt_lvl)) deallocate(this%pt_lvl)
    if (associated(this%dof)) nullify(this%dof)
    if (associated(this%msh)) nullify(this%msh)
    if (associated(this%coef)) nullify(this%coef)
    if (allocated(this%volume_per_gll_lvl)) deallocate(this%volume_per_gll_lvl)
    this%dir = 0
    this%n_el_lvls = 0
    this%n_gll_lvls = 0
 
  end subroutine map_1d_free
 
 
  subroutine map_1d_average_field_list(this, avg_planes, field_list)
    class(map_1d_t), intent(inout) :: this
    type(field_list_t), intent(inout) :: field_list
    type(matrix_t), intent(inout) :: avg_planes
    integer :: n, ierr, j, i
    real(kind=rp) :: coord
    call avg_planes%free()
    call avg_planes%init(this%n_gll_lvls, field_list%size() + 1)
    avg_planes = 0.0_rp
    !ugly way of getting coordinates, computes average
    n = this%dof%size()
    do i = 1, n
       if (this%dir .eq. 1) coord = this%dof%x%x(i,1,1,1)
       if (this%dir .eq. 2) coord = this%dof%y%x(i,1,1,1)
       if (this%dir .eq. 3) coord = this%dof%z%x(i,1,1,1)
       avg_planes%x(this%pt_lvl(i,1,1,1), 1) = &
            avg_planes%x(this%pt_lvl(i,1,1,1), 1) + &
            coord * this%coef%B(i,1,1,1) / &
            this%volume_per_gll_lvl(this%pt_lvl(i,1,1,1))
    end do
    do j = 2, field_list%size() + 1
       do i = 1, n
          avg_planes%x(this%pt_lvl(i,1,1,1), j) = &
               avg_planes%x(this%pt_lvl(i,1,1,1), j) + &
               field_list%items(j-1)%ptr%x(i,1,1,1) * this%coef%B(i,1,1,1) &
               /this%volume_per_gll_lvl(this%pt_lvl(i,1,1,1))
       end do
    end do
    if (pe_size .gt. 1) then
       call mpi_allreduce(mpi_in_place, avg_planes%x, &
            (field_list%size() + 1) * this%n_gll_lvls, mpi_real_precision, &
            mpi_sum, neko_comm, ierr)
    end if
 
  end subroutine map_1d_average_field_list
 
  subroutine map_1d_average_vector_ptr(this, avg_planes, vector_ptr)
    class(map_1d_t), intent(inout) :: this
    !Observe is an array...
    type(vector_ptr_t), intent(inout) :: vector_ptr(:)
    type(matrix_t), intent(inout) :: avg_planes
    integer :: n, ierr, j, i
    real(kind=rp) :: coord
 
    call avg_planes%free()
    call avg_planes%init(this%n_gll_lvls, size(vector_ptr) + 1)
    !ugly way of getting coordinates, computes average
    avg_planes = 0.0_rp
 
    n = this%dof%size()
    do i = 1, n
       if (this%dir .eq. 1) coord = this%dof%x%x(i,1,1,1)
       if (this%dir .eq. 2) coord = this%dof%y%x(i,1,1,1)
       if (this%dir .eq. 3) coord = this%dof%z%x(i,1,1,1)
       avg_planes%x(this%pt_lvl(i,1,1,1), 1) = &
            avg_planes%x(this%pt_lvl(i,1,1,1), 1) + &
            coord * this%coef%B(i,1,1,1) / &
            this%volume_per_gll_lvl(this%pt_lvl(i,1,1,1))
    end do
    do j = 2, size(vector_ptr) + 1
       do i = 1, n
          avg_planes%x(this%pt_lvl(i,1,1,1), j) = &
               avg_planes%x(this%pt_lvl(i,1,1,1), j) + &
               vector_ptr(j-1)%ptr%x(i)*this%coef%B(i,1,1,1) &
               /this%volume_per_gll_lvl(this%pt_lvl(i,1,1,1))
       end do
    end do
    call mpi_allreduce(mpi_in_place, avg_planes%x, &
         (size(vector_ptr) + 1) * this%n_gll_lvls, mpi_real_precision, &
         mpi_sum, neko_comm, ierr)
 
 
  end subroutine map_1d_average_vector_ptr
 
end module map_1d