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 mesh_oct_m
  use messages_oct_m
  use mpi_oct_m,       only: mpi_comm
  use profiling_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
    final     :: centroids_finalize
  end type centroids_t
 
contains
 
  subroutine centroids_init(this, mesh, positions)
    class(centroids_t),  intent(out) :: this
    class(mesh_t),       intent(in ) :: mesh
    real(real64),        intent(in ) :: positions(:, :)
 
    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)
 
    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
 
 
  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)
    type(centroids_t),   intent(inout) :: this
    class(mesh_t),       intent(in   ) :: mesh
 
    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))
 
    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
      ! Domain-specific indices
      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
 
end module centroids_oct_m
 
!! Local Variables:
!! mode: f90
!! coding: utf-8
!! End: