multigrid.F90

Go to the documentation of this file.

!! Copyright (C) 2002-2006 M. Marques, A. Castro, A. Rubio, G. Bertsch, X. Andrade
!! 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 multigrid_oct_m
  use batch_oct_m
  use batch_ops_oct_m
  use boundaries_oct_m
  use debug_oct_m
  use derivatives_oct_m
  use global_oct_m
  use index_oct_m
  use math_oct_m
  use mesh_oct_m
  use mesh_init_oct_m
  use messages_oct_m
  use mpi_oct_m
  use multicomm_oct_m
  use namespace_oct_m
  use parser_oct_m
  use par_vec_oct_m
  use space_oct_m
  use stencil_oct_m
  use transfer_table_oct_m
  use profiling_oct_m
 
  implicit none
 
  private
  public ::                         &
    multigrid_level_t,              &
    multigrid_t,                    &
    multigrid_init,                 &
    multigrid_end,                  &
    multigrid_mesh_half,            &
    multigrid_mesh_double,          &
    multigrid_number_of_levels,     &
    dmultigrid_fine2coarse,         &
    zmultigrid_fine2coarse,         &
    dmultigrid_coarse2fine,         &
    zmultigrid_coarse2fine,         &
    dmultigrid_fine2coarse_batch,   &
    zmultigrid_fine2coarse_batch,   &
    dmultigrid_coarse2fine_batch,   &
    zmultigrid_coarse2fine_batch,   &
    multigrid_get_transfer_tables
 
  integer, parameter, public :: &
    INJECTION  = 1,             &
    fullweight = 2
 
  type multigrid_level_t
    ! Components are public by default
    type(transfer_table_t)          :: tt
    type(mesh_t),          pointer  :: mesh => null()
    type(derivatives_t),   pointer  :: der  => null()
  end type multigrid_level_t
 
  type multigrid_t
    private
    integer, public                              :: n_levels
    type(multigrid_level_t), allocatable, public :: level(:)
 
    integer              :: tp
    integer, allocatable :: sp(:)
    integer, allocatable :: ep(:)
    integer, allocatable :: ep_part(:)
  end type multigrid_t
 
 
contains
 
  ! ---------------------------------------------------------
  subroutine multigrid_init(mgrid, namespace, space, mesh, der, stencil, mc, nlevels)
    type(multigrid_t),     target, intent(out) :: mgrid
    type(namespace_t),             intent(in)  :: namespace
    class(space_t),                intent(in)  :: space
    class(mesh_t),         target, intent(in)  :: mesh
    type(derivatives_t),   target, intent(in)  :: der
    type(stencil_t),               intent(in)  :: stencil
    type(multicomm_t),             intent(in)  :: mc
    integer, optional,             intent(in)  :: nlevels
 
    integer :: i, n_levels, np, order
 
    push_sub(multigrid_init)
 
    if (.not. present(nlevels)) then
      !%Variable MultigridLevels
      !%Type integer
      !%Default 4
      !%Section Mesh
      !%Description
      !% Number of levels in the grid hierarchy used for multigrid. Positive
      !% numbers indicate an absolute number of levels, negative
      !% numbers are subtracted from the maximum number of levels possible.
      !%Option max_levels 0
      !% Calculate the optimal number of levels for the grid.
      !%End
 
      call parse_variable(namespace, 'MultigridLevels', 4, n_levels)
 
      ! default:
      order = der%order
    else
      n_levels = nlevels
      write(message(1), '(a,i1,a)') "Set number of multigrid levels to ", n_levels, ". This ignores the value of MultigridLevels."
      call messages_info(1, namespace=namespace)
 
      !%Variable MultigridDerivativesOrder
      !%Type integer
      !%Default 1
      !%Section Mesh::Derivatives
      !%Description
      !% This variable gives the discretization order for the approximation of
      !% the differential operators on the different levels of the multigrid.
      !% For more details, see the variable DerivativesOrder.
      !% For star-general stencils, the minimum is set to 2.
      !%End
      call parse_variable(namespace, 'MultigridDerivativesOrder', 1, order)
      ! set order to a minimum of 2 for general star stencil, fails otherwise
      if (der%stencil_type == der_stargeneral) then
        order = max(2, order)
      end if
    end if
 
    ! require at least one cell per direction on the coarsest grid and adapt n_levels if needed
    if (any(mesh%idx%ll < 2**n_levels)) then
      message(1) = "Grid too small for the multigrid solver. Reducing the number of levels."
      call messages_warning(1)
      do i = n_levels, 1, -1
        if (all(mesh%idx%ll >= 2**i)) then
          n_levels = i
          exit
        end if
      end do
    end if
 
    if (n_levels <= 0) then
      n_levels = n_levels - 3
      np = mesh%np
 
      do while (np > 1)
        np = np / 8
        n_levels = n_levels + 1
      end do
    else
      n_levels = n_levels - 1
    end if
 
    mgrid%n_levels = n_levels
 
    safe_allocate(mgrid%level(0:n_levels))
 
    mgrid%level(0)%mesh => mesh
    mgrid%level(0)%der => der
 
    mgrid%level(0)%tt%n_fine = mesh%np
    safe_allocate(mgrid%level(0)%tt%fine_i(1:mesh%np))
 
    write(message(1), '(a,i3)') "Multigrid levels:", n_levels + 1
    call messages_info(1, namespace=namespace)
 
    do i = 1, mgrid%n_levels
      safe_allocate(mgrid%level(i)%mesh)
      safe_allocate(mgrid%level(i)%der)
 
      call multigrid_mesh_half(space, namespace, mgrid%level(i-1)%mesh, mgrid%level(i)%mesh, stencil, mc%base_grp)
 
      call derivatives_init(mgrid%level(i)%der, namespace, space, mesh%coord_system, order=order)
 
      call mesh_init_stage_3(mgrid%level(i)%mesh, namespace, space, stencil, mc, parent = mgrid%level(i - 1)%mesh)
 
      call multigrid_get_transfer_tables(mgrid%level(i)%tt, space, mgrid%level(i-1)%mesh, mgrid%level(i)%mesh)
 
      call derivatives_build(mgrid%level(i)%der, namespace, space, mgrid%level(i)%mesh)
 
      call mesh_write_info(mgrid%level(i)%mesh, namespace=namespace)
 
      mgrid%level(i)%der%finer => mgrid%level(i - 1)%der
      mgrid%level(i - 1)%der%coarser => mgrid%level(i)%der
      mgrid%level(i)%der%to_finer => mgrid%level(i)%tt
      mgrid%level(i - 1)%der%to_coarser => mgrid%level(i)%tt
    end do
 
    safe_allocate(mgrid%sp(0:mgrid%n_levels))
    safe_allocate(mgrid%ep(0:mgrid%n_levels))
    safe_allocate(mgrid%ep_part(0:mgrid%n_levels))
 
    mgrid%tp = 0
    do i = 0, mgrid%n_levels
      mgrid%sp(i) = 1 + mgrid%tp
      mgrid%ep(i) = mgrid%tp + mgrid%level(i)%mesh%np
      mgrid%tp = mgrid%tp + mgrid%level(i)%mesh%np_part
      mgrid%ep_part(i) = mgrid%tp
    end do
 
    pop_sub(multigrid_init)
  end subroutine multigrid_init
 
  ! ---------------------------------------------------------
  subroutine multigrid_get_transfer_tables(tt, space, fine, coarse)
    type(transfer_table_t), intent(inout) :: tt
    class(space_t),         intent(in)    :: space
    type(mesh_t),           intent(in)    :: fine, coarse
 
    integer :: i, i1, i2, i4, i8, pt
    integer :: x(space%dim), mod2(space%dim), idx(space%dim)
 
    push_sub(multigrid_get_transfer_tables)
 
    ! Currently this routine only works for 3D
    assert(space%dim == 3)
 
    tt%n_coarse = coarse%np
    safe_allocate(tt%to_coarse(1:tt%n_coarse))
 
    ! GENERATE THE TABLE TO MAP FROM THE FINE TO THE COARSE GRID
    do i = 1, tt%n_coarse
      ! locate the equivalent fine grid point
      call mesh_local_index_to_coords(coarse, i, idx)
      tt%to_coarse(i) = mesh_local_index_from_coords(fine, 2*idx)
    end do
 
    ! count
    tt%n_fine = fine%np
    safe_allocate(tt%fine_i(1:tt%n_fine))
 
    tt%n_fine1 = 0
    tt%n_fine2 = 0
    tt%n_fine4 = 0
    tt%n_fine8 = 0
    do i = 1, tt%n_fine
      call mesh_local_index_to_coords(fine, i, idx)
      mod2 = mod(idx, 2)
 
      pt = sum(abs(mod2))
 
      select case (pt)
      case (0)
        tt%n_fine1 = tt%n_fine1 + 1
        tt%fine_i(i) = 1
      case (1)
        tt%n_fine2 = tt%n_fine2 + 1
        tt%fine_i(i) = 2
      case (2)
        tt%n_fine4 = tt%n_fine4 + 1
        tt%fine_i(i) = 4
      case (3)
        tt%n_fine8 = tt%n_fine8 + 1
        tt%fine_i(i) = 8
      end select
    end do
 
    assert(tt%n_fine1 + tt%n_fine2 + tt%n_fine4 + tt%n_fine8 == tt%n_fine)
 
    safe_allocate(tt%to_fine1(1, 1:tt%n_fine1))
    safe_allocate(tt%to_fine2(1:2, 1:tt%n_fine2))
    safe_allocate(tt%to_fine4(1:4, 1:tt%n_fine4))
    safe_allocate(tt%to_fine8(1:8, 1:tt%n_fine8))
 
    ! and now build the tables
    i1 = 0
    i2 = 0
    i4 = 0
    i8 = 0
    do i = 1, fine%np
      call mesh_local_index_to_coords(fine, i, idx)
      x    = idx/2
      mod2 = mod(idx, 2)
 
      pt = sum(abs(mod2))
 
      select case (pt)
      case (0)
        i1 = i1 + 1
        tt%to_fine1(1, i1) = mesh_local_index_from_coords(coarse, x)
 
      case (1)
        i2 = i2 + 1
        tt%to_fine2(1, i2) = mesh_local_index_from_coords(coarse, x)
        tt%to_fine2(2, i2) = mesh_local_index_from_coords(coarse, x + mod2)
 
      case (2)
        i4 = i4 + 1
        tt%to_fine4(1, i4) = mesh_local_index_from_coords(coarse, [x(1)          , x(2) + mod2(2), x(3) + mod2(3)])
        tt%to_fine4(2, i4) = mesh_local_index_from_coords(coarse, [x(1) + mod2(1), x(2)          , x(3) + mod2(3)])
        tt%to_fine4(3, i4) = mesh_local_index_from_coords(coarse, [x(1) + mod2(1), x(2) + mod2(2), x(3)          ])
        tt%to_fine4(4, i4) = mesh_local_index_from_coords(coarse, x + mod2)
 
      case (3)
        i8 = i8 + 1
        tt%to_fine8(1, i8) = mesh_local_index_from_coords(coarse, x)
        tt%to_fine8(2, i8) = mesh_local_index_from_coords(coarse, [x(1) + mod2(1), x(2)          , x(3)          ])
        tt%to_fine8(3, i8) = mesh_local_index_from_coords(coarse, [x(1)          , x(2) + mod2(2), x(3)          ])
        tt%to_fine8(4, i8) = mesh_local_index_from_coords(coarse, [x(1)          , x(2)          , x(3) + mod2(3)])
        tt%to_fine8(5, i8) = mesh_local_index_from_coords(coarse, [x(1)          , x(2) + mod2(2), x(3) + mod2(3)])
        tt%to_fine8(6, i8) = mesh_local_index_from_coords(coarse, [x(1) + mod2(1), x(2)          , x(3) + mod2(3)])
        tt%to_fine8(7, i8) = mesh_local_index_from_coords(coarse, [x(1) + mod2(1), x(2) + mod2(2), x(3)          ])
        tt%to_fine8(8, i8) = mesh_local_index_from_coords(coarse, x + mod2)
 
      end select
 
    end do
 
    assert(i1 == tt%n_fine1 .and. i2 == tt%n_fine2 .and. i4 == tt%n_fine4 .and. i8 == tt%n_fine8)
 
 
    pop_sub(multigrid_get_transfer_tables)
  end subroutine multigrid_get_transfer_tables
 
  !---------------------------------------------------------------------------------
  !---------------------------------------------------------------------------------
  subroutine multigrid_mesh_half(space, namespace, mesh_in, mesh_out, stencil, grp)
    class(space_t),             intent(in)    :: space
    type(namespace_t),          intent(in)    :: namespace
    type(mesh_t),       target, intent(in)    :: mesh_in
    type(mesh_t),               intent(inout) :: mesh_out
    type(stencil_t),            intent(in)    :: stencil
    type(mpi_grp_t),            intent(in)    :: grp
 
    integer :: idim
 
    push_sub(multigrid_mesh_half)
 
    mesh_out%box              => mesh_in%box
    mesh_out%idx%dim          =  mesh_in%idx%dim
    mesh_out%use_curvilinear  =  mesh_in%use_curvilinear
    mesh_out%use_vol_pp       =  mesh_in%use_vol_pp
    mesh_out%masked_periodic_boundaries = mesh_in%masked_periodic_boundaries
    mesh_out%coord_system     => mesh_in%coord_system
 
    safe_allocate(mesh_out%spacing(1:space%dim))
    mesh_out%spacing(:) = 2*mesh_in%spacing(:)
 
    call index_init(mesh_out%idx, space%dim)
    mesh_out%idx%enlarge(:) = mesh_in%idx%enlarge(:)
    mesh_out%idx%nr(1,:) = (mesh_in%idx%nr(1,:)+mesh_in%idx%enlarge(:))/2
    mesh_out%idx%nr(2,:) = (mesh_in%idx%nr(2,:)-mesh_in%idx%enlarge(:))/2
    mesh_out%idx%ll(:) = mesh_out%idx%nr(2, :) - mesh_out%idx%nr(1, :) + 1
 
    mesh_out%idx%stride(1) = 1
    do idim = 2, mesh_out%idx%dim + 1
      mesh_out%idx%stride(idim) = mesh_out%idx%stride(idim-1) *  &
        (mesh_out%idx%ll(idim-1) + 2*mesh_out%idx%enlarge(idim-1))
    end do
 
    call mesh_init_stage_2(mesh_out, namespace, space, mesh_out%box, stencil, grp)
 
    pop_sub(multigrid_mesh_half)
  end subroutine multigrid_mesh_half
 
  !---------------------------------------------------------------------------------
  subroutine multigrid_mesh_double(space, namespace, mesh_in, mesh_out, stencil, grp)
    class(space_t),             intent(in)    :: space
    type(namespace_t),          intent(in)    :: namespace
    type(mesh_t),       target, intent(in)    :: mesh_in
    type(mesh_t),               intent(inout) :: mesh_out
    type(stencil_t),            intent(in)    :: stencil
    type(mpi_grp_t),            intent(in)    :: grp
 
    integer :: idim
    push_sub(multigrid_mesh_double)
 
    mesh_out%box              => mesh_in%box
    mesh_out%idx%dim          =  mesh_in%idx%dim
    mesh_out%use_curvilinear =  mesh_in%use_curvilinear
    mesh_out%use_vol_pp      =  mesh_in%use_vol_pp
    mesh_out%masked_periodic_boundaries = mesh_in%masked_periodic_boundaries
    mesh_out%coord_system    => mesh_in%coord_system
 
    safe_allocate(mesh_out%spacing(1:space%dim))
    mesh_out%spacing(:) = m_half*mesh_in%spacing(:)
 
    call index_init(mesh_out%idx, space%dim)
    mesh_out%idx%enlarge(:) = mesh_in%idx%enlarge(:)
    mesh_out%idx%nr(1,:) = (mesh_in%idx%nr(1,:)+mesh_in%idx%enlarge(:))*2
    mesh_out%idx%nr(2,:) = (mesh_in%idx%nr(2,:)-mesh_in%idx%enlarge(:))*2
    ! We need to make the possible number of grid points larger by one for
    ! the non-periodic dimensions because the spacing is only half of the
    ! original mesh and thus we could get points in the new boundary
    ! that are still inside the simulation box.
    ! For the periodic dimensions, we are anyway commensurate with the size
    ! of the box, so we are still commensurate when taking twice the number
    ! of points.
    do idim = space%periodic_dim + 1, space%dim
      mesh_out%idx%nr(1, idim) = mesh_out%idx%nr(1, idim) - 1
      mesh_out%idx%nr(2, idim) = mesh_out%idx%nr(2, idim) + 1
    end do
    mesh_out%idx%ll(:) = mesh_out%idx%nr(2, :) - mesh_out%idx%nr(1, :) + 1
 
    mesh_out%idx%stride(1) = 1
    do idim = 2, space%dim+1
      mesh_out%idx%stride(idim) = mesh_out%idx%stride(idim-1) *  &
        (mesh_out%idx%ll(idim-1) + 2*mesh_out%idx%enlarge(idim-1))
    end do
 
    call mesh_init_stage_2(mesh_out, namespace, space, mesh_out%box, stencil, grp)
 
    pop_sub(multigrid_mesh_double)
  end subroutine multigrid_mesh_double
 
  ! ---------------------------------------------------------
  subroutine multigrid_end(mgrid)
    type(multigrid_t), target, intent(inout) :: mgrid
 
    integer :: i
    type(multigrid_level_t), pointer :: level
 
    push_sub(multigrid_end)
 
    safe_deallocate_a(mgrid%sp)
    safe_deallocate_a(mgrid%ep)
    safe_deallocate_a(mgrid%ep_part)
 
    safe_deallocate_a(mgrid%level(0)%tt%fine_i)
 
    do i = 1, mgrid%n_levels
      level => mgrid%level(i)
 
      call derivatives_end(level%der)
      call mesh_end(level%mesh)
      safe_deallocate_p(level%mesh)
      safe_deallocate_p(level%der)
 
      safe_deallocate_a(level%tt%to_coarse)
      safe_deallocate_a(level%tt%to_fine1)
      safe_deallocate_a(level%tt%to_fine2)
      safe_deallocate_a(level%tt%to_fine4)
      safe_deallocate_a(level%tt%to_fine8)
      safe_deallocate_a(level%tt%fine_i)
    end do
 
    safe_deallocate_a(mgrid%level)
 
    pop_sub(multigrid_end)
  end subroutine multigrid_end
 
  !---------------------------------------------------------------------------------
  integer function multigrid_number_of_levels(base_der) result(number)
    type(derivatives_t), target, intent(in)  :: base_der
 
    type(derivatives_t), pointer :: next_der
 
    next_der => base_der%coarser
 
    number = 0
    next_der => base_der%coarser
    do while (associated(next_der))
      number = number + 1
      next_der => next_der%coarser
    end do
 
  end function multigrid_number_of_levels
 
  !---------------------------------------------------------------------------------
  subroutine multigrid_build_stencil(dim, weight, shift)
    integer, intent(in) :: dim
    real(real64), contiguous, intent(inout) :: weight(:)
    integer, contiguous, intent(inout) :: shift(:,:)
 
    integer :: nn, di, dj, dk, dd
 
    push_sub(multigrid_build_stencil)
 
    assert(ubound(weight, dim=1) == 3**dim)
    assert(ubound(shift, dim=1) == dim)
    assert(ubound(shift, dim=2) == 3**dim)
 
    nn = 0
    select case(dim)
    case(1)
      do di = -1, 1
        dd = abs(di)
        nn = nn + 1
        weight(nn) = m_half**dd
        shift(1,nn) = di
      end do
    case(2)
      do di = -1, 1
        do dj = -1, 1
          dd = abs(di)+abs(dj)
          nn = nn + 1
          weight(nn) = m_half**dd
          shift(1,nn) = di
          shift(2,nn) = dj
        end do
      end do
    case(3)
      do di = -1, 1
        do dj = -1, 1
          do dk = -1, 1
            dd = abs(di)+abs(dj)+abs(dk)
            nn = nn + 1
            weight(nn) = m_half**dd
            shift(1,nn) = di
            shift(2,nn) = dj
            shift(3,nn) = dk
          end do
        end do
      end do
    end select
 
    pop_sub(multigrid_build_stencil)
  end subroutine multigrid_build_stencil
 
 
#include "undef.F90"
#include "real.F90"
#include "multigrid_inc.F90"
 
#include "undef.F90"
#include "complex.F90"
#include "multigrid_inc.F90"
 
end module multigrid_oct_m
 
!! Local Variables:
!! mode: f90
!! coding: utf-8
!! End: