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_grp_t, mpi_integer, mpi_integer8
  use namespace_oct_m
  use profiling_oct_m
  use sort_oct_m
 
  implicit none
  private
 
  public :: centroids_t, &
    centroids_read_all_indices
 
  type centroids_t
    private
    integer                     :: ndim
    real(real64),   allocatable :: positions(:, :)
    integer                     :: n_centroids_global
    integer                     :: n_centroids
    integer(int64), allocatable :: icg_to_ipg(:)
    !                                                    Defined in full on every process.
    !                                                    Kept private, but one could add a getter if required
    integer,        allocatable :: ic_to_icg(:)
    integer,        allocatable :: ic_to_ip(:)
 
  contains
    ! Initialize
    procedure :: init => centroids_init
    procedure :: init_from_file => centroids_init_from_file
    ! Getters
    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
    ! Finalize
    procedure :: end => centroids_end
    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 convert_coordinates_to_global_indices(this, mesh)
    call exact_centroids_init_index_maps(this, mesh, optional_default(check_duplicates, .false.))
 
    pop_sub(centroids_init)
 
  end subroutine centroids_init
 
 
  subroutine centroids_init_from_file(this, mesh, fname, check_duplicates)
    class(centroids_t),  intent(out) :: this
    class(mesh_t),       intent(in ) :: mesh
    character(len=*),    intent(in ) :: fname
    logical,  optional,  intent(in ) :: check_duplicates
 
    push_sub(centroids_init_from_file)
 
    this%ndim = mesh%box%dim
    call centroids_read_all_indices(mesh%mpi_grp, fname, this%n_centroids_global, this%icg_to_ipg)
    ! This step might not strictly be necessary, but it populates this%positions
    call convert_global_indices_to_coordinates(this, mesh)
    call exact_centroids_init_index_maps(this, mesh, optional_default(check_duplicates, .false.))
 
    pop_sub(centroids_init_from_file)
 
  end subroutine centroids_init_from_file
 
 
  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
 
 
  subroutine centroids_output_all_indices(this, namespace, grp, fname)
    class(centroids_t), intent(in) :: this
    type(namespace_t),  intent(in) :: namespace
    type(mpi_grp_t),    intent(in) :: grp
    character(len=*),   intent(in) :: fname
 
    integer :: unit, icg
 
    push_sub(centroids_output_all_indices)
 
    write(message(1),'(a)') 'Info: Writing interpolation points to file: ' // trim(adjustl(fname))
    call messages_info(1, namespace=namespace)
 
    unit = io_open(trim(adjustl(fname)), action='write', grp=grp)
 
    if (grp%is_root()) then
      write(unit, *) 'N centroids', this%n_centroids_global
      do icg = 1, this%n_centroids_global
        write(unit, *) this%icg_to_ipg(icg)
      enddo
    endif
 
    call io_close(unit, grp=grp)
 
    pop_sub(centroids_output_all_indices)
 
  end subroutine centroids_output_all_indices
 
 
  subroutine centroids_read_all_indices(grp, fname, n_centroids_global, icg_to_ipg)
    type(mpi_grp_t),             intent(in ) :: grp
    character(len=*),            intent(in ) :: fname
    integer,                     intent(out) :: n_centroids_global
    !                                                               over the whole mesh/system
    integer(int64), allocatable, intent(out) :: icg_to_ipg(:)
    !                                                         to global mesh indices, defined for all processed
 
    integer           :: icg, iunit
    character(len=11) :: dummy
 
    push_sub(centroids_read_all_indices)
 
    iunit = io_open(trim(adjustl(fname)), action='read', grp=grp)
 
    if (grp%is_root()) then
      read(iunit, *) dummy, n_centroids_global
      safe_allocate(icg_to_ipg(1:n_centroids_global))
      do icg = 1, n_centroids_global
        read(iunit, *) icg_to_ipg(icg)
      enddo
    endif
 
    call io_close(iunit, grp=grp)
 
    call grp%bcast(n_centroids_global, 1, mpi_integer, 0)
    if (.not. grp%is_root()) safe_allocate(icg_to_ipg(1:n_centroids_global))
    call grp%bcast(icg_to_ipg, n_centroids_global, mpi_integer8, 0)
 
    pop_sub(centroids_read_all_indices)
 
  end subroutine centroids_read_all_indices
 
 
  subroutine centroids_end(this)
    class(centroids_t),  intent(inout) :: this
 
    push_sub(centroids_end)
 
    safe_deallocate_a(this%positions)
    safe_deallocate_a(this%icg_to_ipg)
    safe_deallocate_a(this%ic_to_icg)
    safe_deallocate_a(this%ic_to_ip)
 
    pop_sub(centroids_end)
 
  end subroutine centroids_end
 
 
  subroutine centroids_finalize(this)
    type(centroids_t),  intent(inout) :: this
 
    push_sub(centroids_finalize)
    call this%end()
    pop_sub(centroids_finalize)
 
  end subroutine centroids_finalize
 
 
  ! Helper functions for initialisation, hence not type-bound or exposed
 
  subroutine convert_coordinates_to_global_indices(this, mesh)
    type(centroids_t),   intent(inout) :: this
    class(mesh_t),       intent(in   ) :: mesh
 
    real(real64)                       :: chi(this%ndim)
    integer                            :: ix(this%ndim)
    integer                            :: icg
    integer(int64)                     :: ipg
 
    push_sub(convert_coordinates_to_global_indices)
 
    ! Convert coordinate to global grid index
    safe_allocate(this%icg_to_ipg(1:this%n_centroids_global))
    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
 
    pop_sub(convert_coordinates_to_global_indices)
 
  end subroutine convert_coordinates_to_global_indices
 
 
  subroutine convert_global_indices_to_coordinates(this, mesh)
    type(centroids_t),   intent(inout) :: this
    class(mesh_t),       intent(in   ) :: mesh
 
    integer                            :: icg
    integer(int64)                     :: ipg
 
    push_sub(convert_global_indices_to_coordinates)
 
    safe_allocate(this%positions(1:this%ndim, 1:this%n_centroids_global))
 
    do icg = 1, this%n_centroids_global
      ipg = this%icg_to_ipg(icg)
      this%positions(1:this%ndim, icg) = mesh_x_global(mesh, ipg)
    enddo
 
    pop_sub(convert_global_indices_to_coordinates)
 
  end subroutine convert_global_indices_to_coordinates
 
 
  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
 
    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)
 
    ! 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
 
    ! 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))
 
    ! 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: