centroids.F90

Go to the documentation of this file.

!! Copyright (C) 2024 A. Buccheri
!!
!! This program is free software; you can redistribute it and/or modify
!! it under the terms of the GNU General Public License as published by
!! the Free Software Foundation; either version 2, or (at your option)
!! any later version.
!!
!! This program is distributed in the hope that it will be useful,
!! but WITHOUT ANY WARRANTY; without even the implied warranty of
!! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
!! GNU General Public License for more details.
!!
!! You should have received a copy of the GNU General Public License
!! along with this program; if not, write to the Free Software
!! Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA
!! 02110-1301, USA.
 
#include "global.h"
 
module centroids_oct_m
  use, intrinsic :: iso_fortran_env
  use comm_oct_m,      only: comm_allreduce
  use debug_oct_m
  use global_oct_m
  use io_oct_m
  use mesh_oct_m
  use messages_oct_m
  use mpi_oct_m,       only: mpi_comm, mpi_world
  use profiling_oct_m
  use sort_oct_m
 
  implicit none
  private
 
  public :: centroids_t
 
  type centroids_t
    private
    integer                     :: ndim
    real(real64),   allocatable :: positions(:, :)
    integer                     :: n_centroids_global
    integer                     :: n_centroids
    integer(int64), allocatable :: icg_to_ipg(:)
    !                                                    Kept private, but one could add a getter if required
    integer,        allocatable :: ic_to_icg(:)
    integer,        allocatable :: ic_to_ip(:)
 
  contains
    procedure :: init => centroids_init
    procedure :: npoints => centroids_get_n_centroids
    procedure :: npoints_global => centroids_get_n_centroids_global
    procedure :: get_local_position => centroids_get_local_position
    procedure :: get_all_positions => centroids_get_all_positions
    ! ic -> ip
    procedure :: to_local_mesh_index => centroids_get_local_mesh_index
    procedure :: local_mesh_indices => centroids_get_local_mesh_indices
    ! ic -> icg
    procedure :: to_global_index => centroids_local_index_to_global_index
    procedure :: global_indices => centroids_local_index_to_global_indices
    ! ic -> ipg
    procedure :: to_global_mesh_index => centroids_global_centroid_to_global_mesh_index
    procedure :: global_mesh_indices => centroids_global_centroid_to_global_mesh_indices
    ! IO
    procedure :: output_all_indices => centroids_output_all_indices
    final     :: centroids_finalize
  end type centroids_t
 
contains
 
  subroutine centroids_init(this, mesh, positions, check_duplicates)
    class(centroids_t),  intent(out) :: this
    class(mesh_t),       intent(in ) :: mesh
    real(real64),        intent(in ) :: positions(:, :)
    logical,  optional,  intent(in ) :: check_duplicates
    !                                      following discretisation of cartesian coordinates
 
    push_sub(centroids_init)
 
    this%ndim = size(positions, 1)
    this%n_centroids_global = size(positions, 2)
    safe_allocate_source(this%positions(1:this%ndim, 1:this%n_centroids_global), positions)
    call exact_centroids_init_index_maps(this, mesh, optional_default(check_duplicates, .false.))
 
    pop_sub(centroids_init)
 
  end subroutine centroids_init
 
 
  pure function centroids_get_n_centroids(this) result(n_centroids)
    class(centroids_t),  intent(in) :: this
    integer                         :: n_centroids
 
    n_centroids = this%n_centroids
 
  end function centroids_get_n_centroids
 
 
  pure function centroids_get_n_centroids_global(this) result(n_centroids_global)
    class(centroids_t),  intent(in) :: this
    integer                         :: n_centroids_global
 
    n_centroids_global = this%n_centroids_global
 
  end function centroids_get_n_centroids_global
 
 
  pure function centroids_global_centroid_to_global_mesh_index(this, ic) result(ipg)
    class(centroids_t),  intent(in) :: this
    integer,             intent(in) :: ic
    integer(int64)                  :: ipg
    integer                         :: icg
 
    icg = this%ic_to_icg(ic)
    ipg = this%icg_to_ipg(icg)
 
  end function centroids_global_centroid_to_global_mesh_index
 
 
  pure function centroids_global_centroid_to_global_mesh_indices(this) result(global_mesh_indices)
    class(centroids_t),  intent(in) :: this
    integer(int64), allocatable     :: global_mesh_indices(:)
 
    integer                         :: ic, icg
 
    allocate(global_mesh_indices(this%n_centroids))
    do ic = 1, this%n_centroids
      icg = this%ic_to_icg(ic)
      global_mesh_indices(ic) = this%icg_to_ipg(icg)
    enddo
 
  end function centroids_global_centroid_to_global_mesh_indices
 
 
  pure function centroids_get_local_mesh_index(this, ic) result(ip)
    class(centroids_t),  intent(in) :: this
    integer,             intent(in) :: ic
    integer                         :: ip
 
    ip = this%ic_to_ip(ic)
 
  end function centroids_get_local_mesh_index
 
 
  pure function centroids_get_local_mesh_indices(this) result(mesh_indices)
    class(centroids_t),  intent(in) :: this
    integer, allocatable            :: mesh_indices(:)
 
    allocate(mesh_indices(this%n_centroids), source=this%ic_to_ip)
 
  end function centroids_get_local_mesh_indices
 
 
  pure function centroids_local_index_to_global_index(this, ic) result(icg)
    class(centroids_t),  intent(in) :: this
    integer,             intent(in) :: ic
    integer                         :: icg
 
    icg = this%ic_to_icg(ic)
 
  end function centroids_local_index_to_global_index
 
 
  pure function centroids_local_index_to_global_indices(this) result(global_centroid_indices)
    class(centroids_t),  intent(in) :: this
    integer, allocatable            :: global_centroid_indices(:)
 
    allocate(global_centroid_indices(this%n_centroids), source=this%ic_to_icg)
 
  end function centroids_local_index_to_global_indices
 
 
  pure function centroids_get_local_position(this, mesh, ic) result(pos)
    class(centroids_t),  intent(in) :: this
    class(mesh_t),       intent(in) :: mesh
    integer,             intent(in) :: ic
 
    real(real64)                    :: pos(this%ndim)
 
    pos(:) = mesh%x(this%ic_to_ip(ic), :)
 
  end function centroids_get_local_position
 
 
  pure function centroids_get_all_positions(this) result(positions)
    class(centroids_t),  intent(in) :: this
    real(real64), allocatable :: positions(:, :)
 
    allocate(positions(this%ndim, this%n_centroids_global), source=this%positions)
 
  end function centroids_get_all_positions
 
 
  ! TODO(Alex) Would be useful to output a domain index
  subroutine centroids_output_all_indices(this, fname)
    class(centroids_t), intent(in) :: this
    character(len=*),   intent(in) :: fname
 
    integer :: unit, icg
 
    push_sub(centroids_output_all_indices)
 
    unit = io_open(trim(adjustl(fname)), action='write')
 
    write(unit, *) 'N centroids', this%n_centroids_global
    do icg = 1, this%n_centroids_global
      write(unit, *) this%icg_to_ipg(icg)
    enddo
 
    call io_close(unit)
 
    pop_sub(centroids_output_all_indices)
 
  end subroutine centroids_output_all_indices
 
 
  subroutine centroids_finalize(this)
    type(centroids_t),  intent(inout) :: this
 
    safe_deallocate_a(this%positions)
    safe_deallocate_a(this%icg_to_ipg)
    safe_deallocate_a(this%ic_to_icg)
 
  end subroutine centroids_finalize
 
  ! Helper function for initialisation, hence not type-bound or exposed
 
  subroutine exact_centroids_init_index_maps(this, mesh, check_duplicates)
    type(centroids_t),   intent(inout) :: this
    class(mesh_t),       intent(in   ) :: mesh
    logical,             intent(in   ) :: check_duplicates
    !                                           following discretisation of cartesian coordinates
 
    real(real64)                       :: chi(this%ndim)
    integer                            :: ix(this%ndim)
    integer                            :: icg, ip
    integer(int64)                     :: ipg
    integer, allocatable               :: tmp_ic_to_icg(:), tmp_ic_to_ip(:)
#if !defined(NDEBUG)
    integer                            :: n_centroids_global
#endif
 
    push_sub(exact_centroids_init_index_maps)
 
    safe_allocate(this%icg_to_ipg(1:this%n_centroids_global))
 
    ! Local to this process
    this%n_centroids = 0
    ! Allocated with an upper bound
    safe_allocate(tmp_ic_to_icg(1:this%n_centroids_global))
    safe_allocate(tmp_ic_to_ip(1:this%n_centroids_global))
 
    ! Convert coordinate to global grid index
    do icg = 1, this%n_centroids_global
      chi = mesh%coord_system%from_cartesian(this%positions(:, icg))
      ! See the definition of mesh_x_global, for example
      ix = nint(chi / mesh%spacing)
      ipg = mesh_global_index_from_coords(mesh, ix)
      this%icg_to_ipg(icg) = ipg
    enddo
 
    ! Having discretised the coordinates to grid indices, check for and remove any duplicate points
    ! TODO(Alex) Annoyance with this is that it changes the order of the centroid indices
    ! compared to not running check_duplicates, due to use of sorting.
    ! Additionally, this may not be the desired behaviour. I.e. it might be better to select
    ! the next-best discrete point, or stop Octopus
    if (check_duplicates) then
      call remove_duplicates_int64(this%icg_to_ipg)
      this%n_centroids_global = size(this%icg_to_ipg)
    endif
 
    ! Assign grid indices local to this domain
    do icg = 1, this%n_centroids_global
      ipg = this%icg_to_ipg(icg)
      ip = mesh_global2local(mesh, ipg)
      if (ip > 0 .and. ip <= mesh%np) then
        this%n_centroids = this%n_centroids + 1
        tmp_ic_to_icg(this%n_centroids) = icg
        tmp_ic_to_ip(this%n_centroids) = ip
      endif
    enddo
 
#if !defined(NDEBUG)
    ! Total number of centroids must be conserved over all domains
    n_centroids_global = this%n_centroids
    call comm_allreduce(mesh%mpi_grp, n_centroids_global)
    assert(n_centroids_global == this%n_centroids_global)
#endif
 
    ! Reallocate with correct size
    safe_allocate_source(this%ic_to_icg(1:this%n_centroids), tmp_ic_to_icg(1:this%n_centroids))
    safe_deallocate_a(tmp_ic_to_icg)
 
    safe_allocate_source(this%ic_to_ip(1:this%n_centroids), tmp_ic_to_ip(1:this%n_centroids))
    safe_deallocate_a(tmp_ic_to_ip)
 
    pop_sub(exact_centroids_init_index_maps)
 
  end subroutine exact_centroids_init_index_maps
 
  subroutine remove_duplicates_int64(indices)
    integer(int64), intent(inout), allocatable :: indices(:)
    ! Out: Ordered set of indices, with any duplicates removed
 
    integer                        :: n_duplicates, cnt, n, i
    integer(int64),    allocatable :: unique_indices(:)
 
    push_sub(remove_duplicates_int64)
 
    n = size(indices)
    call sort(indices)
 
    allocate(unique_indices(n))
 
    n_duplicates = 0
    cnt = 1
    unique_indices(1) = indices(1)
 
    do i = 1, n - 1
      if (indices(i) == indices(i+1)) then
        n_duplicates = n_duplicates + 1
        cycle
      endif
      cnt = cnt + 1
      unique_indices(cnt) = indices(i+1)
    enddo
 
    safe_deallocate_a(indices)
    safe_allocate_source(indices(1:cnt), unique_indices(1:cnt))
    safe_deallocate_a(unique_indices)
 
    if (debug%info .and. n_duplicates > 0) then
      write(message(1), '(I0)') n_duplicates
      message(1) = trim(message(1)) // " duplicates found after discretisation of centroid positions"
      write(message(2), '(I0)') cnt
      message(2) = " Removed the duplicate/s and resized 'indices' array to " // trim(message(2))
      call messages_warning(2)
    endif
 
    pop_sub(remove_duplicates_int64)
 
  end subroutine remove_duplicates_int64
 
end module centroids_oct_m
 
!! Local Variables:
!! mode: f90
!! coding: utf-8
!! End: