mesh_partition.F90

Go to the documentation of this file.

!! Copyright (C) 2002-2013 M. Marques, A. Castro, A. Rubio, G. Bertsch, M. Oliveira
!! Copyright (C) 2021 S. Ohlmann
!!
!! 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 mesh_partition_oct_m
  use debug_oct_m
  use global_oct_m
  use index_oct_m
  use io_oct_m
  use mesh_oct_m
  use metis_oct_m
  use messages_oct_m
  use mpi_oct_m
  use namespace_oct_m
  use parser_oct_m
  use partition_oct_m
  use profiling_oct_m
  use restart_oct_m
  use space_oct_m
  use stencil_oct_m
  use stencil_star_oct_m
 
  implicit none
 
  private
  public ::                      &
    mesh_partition,              &
    mesh_partition_from_parent,  &
    mesh_partition_dump,         &
    mesh_partition_load,         &
    mesh_partition_write_info,   &
    mesh_partition_messages_debug
 
  integer, parameter :: &
    METIS    = 1,       &
    parmetis = 2,       &
    hilbert  = 3
 
  integer, parameter :: &
    RCB    = 1,         &
    graph  = 2
 
contains
 
  ! ---------------------------------------------------------------
  ! ---------------------------------------------------------------
  subroutine mesh_partition(mesh, namespace, space, lapl_stencil, vsize)
    type(mesh_t),      intent(inout)  :: mesh
    type(namespace_t), intent(in)     :: namespace
    class(space_t),    intent(in)     :: space
    type(stencil_t),   intent(in)     :: lapl_stencil
    integer,           intent(in)     :: vsize
 
    integer              :: jp, jpart
    integer(imetis)      :: iv
    integer(int64)          :: ivtmp, inb
    integer              :: ix(space%dim), jx(space%dim)
 
    integer              :: npart
    integer              :: ipart
 
    integer(imetis)              :: ne_global
    integer(imetis)              :: nv_global
    integer(imetis)              :: ne
    integer                      :: nv
 
    integer(imetis), allocatable :: xadj_global(:)
    integer(imetis), allocatable :: adjncy_global(:)
    integer(imetis), allocatable :: adjncy(:)
    integer(imetis), allocatable :: xadj(:)
 
    integer(imetis), allocatable :: options(:)
#ifdef HAVE_METIS
    integer(imetis)              :: edgecut
#endif
    real(real32),        allocatable :: tpwgts(:)
 
    integer                      :: iunit
    integer, allocatable :: rcounts(:)
    integer(imetis)      :: tmp
    integer(imetis), allocatable :: rdispls(:)
    integer(imetis), allocatable :: part(:)
    integer(imetis), allocatable :: part_global(:)
    integer(imetis), allocatable :: vtxdist(:)
 
    type(stencil_t) :: stencil
    integer :: stencil_to_use, default_method, method
    integer :: library
    integer, parameter   :: STAR = 1, laplacian = 2
    integer(imetis), allocatable :: istart(:)
    integer, allocatable :: lsize(:)
 
    integer :: default, ierr
 
    push_sub_with_profile(mesh_partition)
 
    if (mesh%np_global == 0) then
      message(1) = 'The mesh is empty and cannot be partitioned.'
      call messages_fatal(1, namespace=namespace)
    end if
 
    !%Variable MeshPartitionPackage
    !%Type integer
    !%Section Execution::Parallelization
    !%Description
    !% Decides which library to use to perform the mesh partition.
    !% By default ParMETIS is used when available, otherwise METIS is used.
    !%Option metis 1
    !% METIS library.
    !%Option parmetis 2
    !% (Experimental) Use ParMETIS library to perform the mesh partition.
    !% Only available if the code was compiled with ParMETIS support.
    !%Option part_hilbert 3
    !% Use the ordering along the Hilbert curve for partitioning.
    !%End
    default = metis
#ifdef HAVE_PARMETIS
    default = parmetis
#endif
    call parse_variable(namespace, 'MeshPartitionPackage', default, library)
 
#if !defined(HAVE_METIS)
    if (library == metis) then
      message(1) = 'METIS was requested, but Octopus was compiled without it.'
      call messages_fatal(1, only_root_writes = .true.)
    end if
#endif
#if !defined(HAVE_PARMETIS)
    if (library == parmetis) then
      message(1) = 'PARMETIS was requested, but Octopus was compiled without it.'
      call messages_fatal(1, only_root_writes = .true.)
    end if
#endif
 
    !%Variable MeshPartitionStencil
    !%Type integer
    !%Default stencil_star
    !%Section Execution::Parallelization
    !%Description
    !% To partition the mesh, it is necessary to calculate the connection
    !% graph connecting the points. This variable selects which stencil
    !% is used to do this.
    !%Option stencil_star 1
    !% An order-one star stencil.
    !%Option laplacian 2
    !% The stencil used for the Laplacian is used to calculate the
    !% partition. This in principle should give a better partition, but
    !% it is slower and requires more memory.
    !%End
    call parse_variable(namespace, 'MeshPartitionStencil', star, stencil_to_use)
 
    if (stencil_to_use == star) then
      call stencil_star_get_lapl(stencil, mesh%box%dim, order = 1)
    else if (stencil_to_use == laplacian) then
      call stencil_copy(lapl_stencil, stencil)
    else
      call messages_input_error(namespace, 'MeshPartitionStencil')
    end if
 
    ! Shortcut to the global number of vertices
    if (mesh%np_global > huge(nv_global)) then
      message(1) = "Error: too many grid points for this version of metis/parmetis."
      message(2) = "Please use the 64-bit version."
      call messages_fatal(2, namespace=namespace)
    end if
    nv_global = i8_to_imetis(mesh%np_global)
 
    call partition_init(mesh%partition, imetis_to_i8(nv_global), mesh%mpi_grp)
 
    ! Get number of partitions.
    npart = mesh%mpi_grp%size
    ipart = mesh%mpi_grp%rank + 1
 
    safe_allocate(istart(1:npart))
    safe_allocate(lsize(1:npart))
    safe_allocate(vtxdist(1:npart+1))
    call partition_get_local_size(mesh%partition, ivtmp, nv)
    assert(ivtmp <= huge(0_imetis))
    iv = i8_to_imetis(ivtmp)
 
    ! Allocate local output matrix
    safe_allocate(part(1:nv))
    if (library == hilbert) then
      part = ipart
    else
 
      call mesh%mpi_grp%allgather(iv, 1, mpi_metis_int, istart(1), 1, mpi_metis_int)
      call mesh%mpi_grp%allgather(nv, 1, mpi_integer, lsize(1), 1, mpi_integer)
 
      vtxdist(1:npart) = istart(1:npart)
      vtxdist(npart+1) = nv_global + 1
 
      ! Create graph with each point being represented by a
      ! vertex and edges between neighbouring points.
      safe_allocate(xadj(1:nv+1))
      safe_allocate(adjncy(1:i4_to_imetis(stencil%size - 1)*nv))
      ne = 1
      ! Iterate over number of vertices.
      do iv = 1, lsize(ipart)
        ! Get coordinates of point iv (vertex iv).
        call mesh_global_index_to_coords(mesh, imetis_to_i8(istart(ipart)-1 + iv), ix)
 
        ! Set entry in index table.
        xadj(iv) = ne
        ! Check all possible neighbours.
        do jp = 1, stencil%size
          if (jp == stencil%center) cycle
 
          ! Store coordinates of possible neighbors, they
          ! are needed several times in the check below.
          jx = ix + stencil%points(:, jp)
 
          if (all(jx >= mesh%idx%nr(1, :)) .and. all(jx <= mesh%idx%nr(2, :))) then
            ! Only points inside the mesh or its enlargement
            ! are included in the graph.
            inb = mesh_global_index_from_coords(mesh, jx)
            if (inb /= 0 .and. inb <= nv_global) then
              ! Store a new edge and increment edge counter.
              adjncy(ne) = i8_to_imetis(inb)
              ne = ne + 1
            end if
          end if
        end do
      end do
      ! We started with ne=1 for simplicity, so ne is off by one in the end.
      ne = ne - 1
 
      ! Set the total number of edges plus one as last index.
      ! (NOTE: the plus one is because we are using Fortran numbering)
      ! The reason is: neighbours of node i are stored in adjncy(xadj(i):xadj(i+1)-1).
      ! Setting the last index as mentioned makes special handling of last element
      ! unnecessary (this indexing is a METIS requirement).
      xadj(nv+1) = ne + 1
 
 
      if (debug%info .or. library == metis) then
        !Gather the global xadj and adjncy arrays
        safe_allocate(rcounts(1:npart))
        safe_allocate(rdispls(1:npart))
 
        safe_allocate(xadj_global(1:nv_global + 1))
        xadj_global(1) = 1
        rcounts(1:npart) = lsize(1:npart)
        rdispls(1:npart) = vtxdist(1:npart) - 1
        call mesh%mpi_grp%allgatherv(xadj(2:), lsize(ipart), mpi_metis_int, &
          xadj_global(2:), rcounts, rdispls, mpi_metis_int)
        do jpart = 2, npart
          do iv = 1, lsize(jpart)
            xadj_global(istart(jpart) + iv) = xadj_global(istart(jpart) + iv) + xadj_global(istart(jpart)) - 1
          end do
        end do
 
        ne_global = xadj_global(nv_global + 1)
        safe_allocate(adjncy_global(1:ne_global))
        do jpart = 1, npart
          rdispls(jpart) = xadj_global(vtxdist(jpart)) - 1
          tmp = xadj_global(vtxdist(jpart+1)) - 1 - rdispls(jpart)
          assert(tmp < huge(0_int32))
          rcounts(jpart) = imetis_to_i4(tmp)
        end do
        assert(xadj(nv+1)-1 < huge(0_int32))
        call mesh%mpi_grp%allgatherv(adjncy, imetis_to_i4(xadj(nv+1)-1), mpi_metis_int, &
          adjncy_global, rcounts, rdispls, mpi_metis_int)
 
        safe_deallocate_a(rcounts)
        safe_deallocate_a(rdispls)
      end if
 
 
      if (debug%info) then
        ! DEBUG output. Write graph to file mesh_graph.txt.
        message(1) = 'Info: Adjacency lists of the graph representing the grid'
        message(2) = 'Info: are stored in debug/mesh_partition/mesh_graph.txt.'
        message(3) = 'Info: Compatible with METIS programs pmetis and kmetis.'
        message(4) = 'Info: First line contains number of vertices and edges.'
        message(5) = 'Info: Edges are not directed and appear twice in the lists.'
        call messages_info(5)
        if (mpi_grp_is_root(mpi_world)) then
          call io_mkdir('debug/mesh_partition', namespace)
          iunit = io_open('debug/mesh_partition/mesh_graph.txt', namespace, action='write')
          write(iunit, *) nv_global, ne_global/2
          do iv = 1, nv
            write(iunit, *) adjncy_global(xadj_global(iv):xadj_global(iv+1) - 1)
          end do
          call io_close(iunit)
        end if
        message(1) = "Info: Done writing mesh_graph.txt."
        call messages_info(1, namespace=namespace)
      end if
 
      ! The sum of all tpwgts elements has to be 1 and
      ! we don`t care about the weights. So; 1/npart
      safe_allocate(tpwgts(1:vsize))
      tpwgts(1:vsize) = real(1.0, real32)/real(vsize, real32)
 
      select case (library)
      case (metis)
        safe_allocate(options(1:40))
        options = 0
#ifdef HAVE_METIS
        call oct_metis_setdefaultoptions(options(1)) ! is equal to: options = -1
#endif
        options(metis_option_numbering) = 1 ! Fortran style: start counting from 1
 
        if (vsize  <  8) then
          default_method = rcb
        else
          default_method = graph
        end if
 
        !%Variable MeshPartition
        !%Type integer
        !%Section Execution::Parallelization
        !%Description
        !% When using METIS to perform the mesh partitioning, decides which
        !% algorithm is used. By default, <tt>graph</tt> partitioning
        !% is used for 8 or more partitions, and <tt>rcb</tt> for fewer.
        !%Option rcb 1
        !% Recursive coordinate bisection partitioning.
        !%Option graph 2
        !% Graph partitioning (called 'k-way' by METIS).
        !%End
        call parse_variable(namespace, 'MeshPartition', default_method, method)
 
        safe_allocate(part_global(1:nv_global))
 
        !Now we can call METIS
        select case (method)
        case (rcb)
          message(1) = 'Info: Using METIS 5 multilevel recursive bisection to partition the mesh.'
          call messages_info(1, namespace=namespace)
#ifdef HAVE_METIS
          ierr = oct_metis_partgraphrecursive(nv_global, 1_imetis, xadj_global(1), adjncy_global(1), &
            i4_to_imetis(vsize), tpwgts(1), 1.01_4, options(1), edgecut, part_global(1))
#endif
          call metis_error_code(ierr)
        case (graph)
          message(1) = 'Info: Using METIS 5 multilevel k-way algorithm to partition the mesh.'
          call messages_info(1, namespace=namespace)
#ifdef HAVE_METIS
          ierr = oct_metis_partgraphkway(nv_global, 1_imetis, xadj_global(1), adjncy_global(1), &
            i4_to_imetis(vsize), tpwgts(1), 1.01_4, options(1), edgecut, part_global(1))
#endif
          call metis_error_code(ierr)
        case default
          message(1) = 'Selected partition method is not available in METIS 5.'
          call messages_fatal(1, namespace=namespace)
        end select
 
        part(1:nv) = part_global(istart(ipart):istart(ipart) + nv - 1)
 
        safe_deallocate_a(options)
        safe_deallocate_a(part_global)
 
      case (parmetis)
 
        safe_allocate(options(1:3))
        options = 0 ! For the moment we use default options
 
        write(message(1),'(a)') 'Info: Using ParMETIS multilevel k-way algorithm to partition the mesh.'
        call messages_info(1, namespace=namespace)
 
        ! Call to ParMETIS with no imbalance tolerance
#ifdef HAVE_PARMETIS
        call oct_parmetis_v3_partkway(vtxdist(1), xadj(1), adjncy(1), &
          1_imetis, i4_to_imetis(mesh%mpi_grp%size), tpwgts(1), 1.05_4, &
          options(1), edgecut, part(1), mesh%mpi_grp%comm%MPI_VAL)
#endif
 
        ! Deallocate matrices
        safe_deallocate_a(options)
 
      end select
 
      safe_deallocate_a(vtxdist)
 
      if (debug%info .or. library == metis) then
        safe_deallocate_a(xadj_global)
        safe_deallocate_a(adjncy_global)
      end if
      safe_deallocate_a(tpwgts)
      safe_deallocate_a(xadj)
      safe_deallocate_a(adjncy)
      safe_deallocate_a(istart)
      safe_deallocate_a(lsize)
    end if
 
    assert(all(part(1:nv) > 0))
    assert(all(part(1:nv) <= vsize))
    call partition_set(mesh%partition, imetis_to_i4(part))
 
    safe_deallocate_a(part)
 
    call stencil_end(stencil)
    pop_sub_with_profile(mesh_partition)
 
  contains
 
    subroutine metis_error_code(ierr)
      integer, intent(in) :: ierr
 
      push_sub(mesh_partition.metis_error_code)
 
      select case (ierr)
      case (metis_ok)
        !Everything OK
      case (metis_error_input)
        message(1) = "Metis: Input error."
        call messages_fatal(1, namespace=namespace)
      case (metis_error_memory)
        message(1) = "Metis: Unable to allocate required memory."
        call messages_fatal(1, namespace=namespace)
      case (metis_error)
        message(1) = "Metis: Some type of error."
        call messages_fatal(1, namespace=namespace)
      end select
 
      pop_sub(mesh_partition.metis_error_code)
    end subroutine metis_error_code
 
  end subroutine mesh_partition
 
  ! ----------------------------------------------------
  subroutine mesh_partition_from_parent(mesh, parent)
    type(mesh_t), intent(inout) :: mesh
    type(mesh_t), intent(in)    :: parent
 
    integer :: np, ip_local
    integer(int64) :: istart, ip_global
    integer(int64), allocatable :: points(:)
    integer, allocatable :: part(:)
    integer :: idx(1:3)
 
    push_sub(mesh_partition_from_parent)
 
    !Initialize partition
    call partition_init(mesh%partition, mesh%np_global, mesh%mpi_grp)
    call partition_get_local_size(mesh%partition, istart, np)
 
    !Get the partition number of each point from the parent mesh
    safe_allocate(points(1:np))
    safe_allocate(part(1:np))
    do ip_local = 1, np
      ip_global = istart + ip_local - 1
      call mesh_global_index_to_coords(mesh, ip_global, idx)
      points(ip_local) = mesh_global_index_from_coords(parent, 2*idx)
    end do
    call partition_get_partition_number(parent%partition, np, points, part)
 
    !Set the partition
    call partition_set(mesh%partition, part)
 
    !Free memory
    safe_deallocate_a(points)
    safe_deallocate_a(part)
 
    pop_sub(mesh_partition_from_parent)
  end subroutine mesh_partition_from_parent
 
  ! ----------------------------------------------------
  subroutine mesh_partition_dump(restart, mesh, vsize, ierr)
    type(restart_t), intent(in)  :: restart
    type(mesh_t),    intent(in)  :: mesh
    integer,         intent(in)  :: vsize
    integer,         intent(out) :: ierr
 
    integer :: err
    character(len=6) :: numstring
 
    push_sub_with_profile(mesh_partition_dump)
 
    ierr = 0
 
    if (restart_skip(restart)) then
      pop_sub_with_profile(mesh_partition_dump)
      return
    end if
 
    if (debug%info) then
      message(1) = "Debug: Writing mesh partition restart."
      call messages_info(1)
    end if
 
    write(numstring, '(i6.6)') vsize
 
    call partition_dump(mesh%partition, restart_dir(restart), 'inner_partition_'//trim(numstring)//'.obf', err)
    if (err /= 0) then
      message(1) = "Unable to write inner mesh partition to 'inner_partition_"//trim(numstring)//".obf'"
      call messages_warning(1)
      ierr = ierr + 1
    end if
 
    if (debug%info) then
      message(1) = "Debug: Writing mesh partition restart done."
      call messages_info(1)
    end if
 
    pop_sub_with_profile(mesh_partition_dump)
  end subroutine mesh_partition_dump
 
  ! ----------------------------------------------------
  subroutine mesh_partition_load(restart, mesh, ierr)
    type(restart_t), intent(in)    :: restart
    type(mesh_t),    intent(inout) :: mesh
    integer,         intent(out)   :: ierr
 
    integer :: err
    character(len=6) :: numstring
    character(len=MAX_PATH_LEN) :: filename, dir
 
    push_sub(mesh_partition_load)
 
    ierr = 0
 
    if (restart_skip(restart)) then
      ierr = -1
      pop_sub(mesh_partition_load)
      return
    end if
 
    if (debug%info) then
      message(1) = "Debug: Reading mesh partition restart."
      call messages_info(1)
    end if
 
    call partition_init(mesh%partition, mesh%np_global, mesh%mpi_grp)
 
    write(numstring, '(i6.6)') mesh%mpi_grp%size
    dir = restart_dir(restart)
 
    !Read inner partition
    filename = 'inner_partition_'//numstring//'.obf'
    if (.not. io_file_exists(trim(dir)//"/"//filename)) then
      message(1) = "Unable to read inner mesh partition, file '"//trim(filename)//"' does not exist."
      call messages_warning(1)
      ierr = ierr + 1
    else
      call partition_load(mesh%partition, restart_dir(restart), filename, err)
      if (err /= 0) then
        message(1) = "Unable to read inner mesh partition from '"//trim(filename)//"'."
        call messages_warning(1)
        ierr = ierr + 2
      end if
    end if
 
    ! Free the memory in case we were unable to read the partitions
    if (ierr /= 0) then
      call partition_end(mesh%partition)
    end if
 
    if (debug%info) then
      message(1) = "Debug: Reading mesh partition restart done."
      call messages_info(1)
    end if
 
    pop_sub(mesh_partition_load)
  end subroutine mesh_partition_load
 
 
  ! ----------------------------------------------------------------------
  subroutine mesh_partition_write_info(mesh, iunit, namespace)
    type(mesh_t),                intent(in) :: mesh
    integer,           optional, intent(in) :: iunit
    type(namespace_t), optional, intent(in) :: namespace
 
    integer(int64) :: npoints
    integer :: npart
    integer, allocatable :: nghost(:), nbound(:), nlocal(:), nneigh(:)
    integer :: num_ghost, num_bound, num_local, num_neigh
    real(real64) :: quality
 
    integer :: ipart
    logical, allocatable :: is_a_neigh(:)
    real(real64) :: scal
 
    push_sub_with_profile(mesh_partition_write_info)
 
    ! Build information about the mesh partition
    npart = mesh%mpi_grp%size
    npoints = mesh%np_part_global
 
    safe_allocate(nghost(1:npart))
    safe_allocate(nbound(1:npart))
    safe_allocate(nlocal(1:npart))
    safe_allocate(nneigh(1:npart))
    safe_allocate(is_a_neigh(1:npart))
 
    ipart = mesh%mpi_grp%rank + 1
 
    num_local = mesh%np
    num_ghost = mesh%pv%np_ghost
    num_bound = mesh%pv%np_bndry
 
    is_a_neigh = mesh%pv%ghost_rcounts > 0
    num_neigh = count(is_a_neigh(1:npart))
 
    safe_deallocate_a(is_a_neigh)
 
    call mesh%mpi_grp%gather(num_neigh, 1, mpi_integer, nneigh(1), 1, mpi_integer, 0)
    call mesh%mpi_grp%gather(num_local, 1, mpi_integer, nlocal(1), 1, mpi_integer, 0)
    call mesh%mpi_grp%gather(num_ghost, 1, mpi_integer, nghost(1), 1, mpi_integer, 0)
    call mesh%mpi_grp%gather(num_bound, 1, mpi_integer, nbound(1), 1, mpi_integer, 0)
 
    ! Calculate partition quality
    scal = real(npart, real64) /npoints
 
    quality = m_zero
 
    quality = quality + (real(maxval(nlocal), real64) - real(minval(nlocal), real64))**3
    quality = quality + (sum(real(nghost, real64) **2))
 
    quality = m_one/(m_one + quality)
 
 
    ! Only the root has the information (nlocal, nghost, ...)
    if (mpi_grp_is_root(mesh%mpi_grp)) then
      ! Write information about the partition
      message(1) = &
        'Info: Mesh partition:'
      message(2) = ''
      call messages_info(2, iunit=iunit, namespace=namespace)
 
      write(message(1),'(a,e16.6)') &
        '      Partition quality:', quality
      message(2) = ''
      call messages_info(2, iunit=iunit, namespace=namespace)
 
      write(message(1),'(a)') &
        '                 Neighbours         Ghost points'
      write(message(2),'(a,i5,a,i10)') &
        '      Average  :      ', nint(sum(real(nneigh, real64) )/npart), '           ', nint(sum(real(nghost, real64) )/npart)
      write(message(3),'(a,i5,a,i10)') &
        '      Minimum  :      ', minval(nneigh),    '           ', minval(nghost)
      write(message(4),'(a,i5,a,i10)') &
        '      Maximum  :      ', maxval(nneigh),    '           ', maxval(nghost)
      message(5) = ''
      call messages_info(5, iunit=iunit, namespace=namespace)
 
      do ipart = 1, npart
        write(message(1),'(a,i5)')  &
          '      Nodes in domain-group  ', ipart
        write(message(2),'(a,i10,a,i10)') &
          '        Neighbours     :', nneigh(ipart), &
          '        Local points    :', nlocal(ipart)
        write(message(3),'(a,i10,a,i10)') &
          '        Ghost points   :', nghost(ipart), &
          '        Boundary points :', nbound(ipart)
        call messages_info(3, iunit=iunit, namespace=namespace)
      end do
 
      message(1) = ''
      call messages_info(1, iunit=iunit, namespace=namespace)
    end if
 
    safe_deallocate_a(nghost)
    safe_deallocate_a(nbound)
    safe_deallocate_a(nlocal)
    safe_deallocate_a(nneigh)
 
    pop_sub_with_profile(mesh_partition_write_info)
  end subroutine mesh_partition_write_info
 
  ! ----------------------------------------------------
  subroutine mesh_partition_messages_debug(mesh, namespace)
    type(mesh_t),      intent(in)    :: mesh
    type(namespace_t), intent(in)    :: namespace
 
    integer              :: ip
    integer(int64)          :: ipg
    integer              :: iunit          ! For debug output to files.
    character(len=6)     :: filenum
 
    if (.not. debug%info) return
 
    push_sub(mesh_partition_messages_debug)
 
    call io_mkdir('debug/mesh_partition', namespace)
 
    ! Debug output. Write points of each partition in a different file.
    write(filenum, '(i6.6)') mesh%mpi_grp%rank+1
 
    ! without boundary
    iunit = io_open('debug/mesh_partition/mesh_partition.'//filenum, &
      namespace, action='write')
    do ip = 1, mesh%np
      ipg = mesh_local2global(mesh, ip)
      write(iunit, '(i8,99f18.8)') ipg, mesh_x_global(mesh, ipg)
    end do
    call io_close(iunit)
 
    ! with boundary included
    iunit = io_open('debug/mesh_partition/mesh_partition_all.'//filenum, &
      namespace, action='write')
    do ip = 1, mesh%np_part
      ipg = mesh_local2global(mesh, ip)
      write(iunit, '(i8,99f18.8)') ipg, mesh_x_global(mesh, ipg)
    end do
    call io_close(iunit)
 
    ! points from enlargement
    if (mpi_grp_is_root(mpi_world)) then
      iunit = io_open('debug/mesh_partition/mesh_partition_boundary', &
        namespace, action='write')
      do ipg = mesh%np_global+1, mesh%np_part_global
        write(iunit, '(i8,99f18.8)') ipg, mesh_x_global(mesh, ipg)
      end do
      call io_close(iunit)
    end if
 
    pop_sub(mesh_partition_messages_debug)
 
  end subroutine mesh_partition_messages_debug
 
end module mesh_partition_oct_m
 
!! Local Variables:
!! mode: f90
!! coding: utf-8
!! End: