multigrid_solver.F90

Go to the documentation of this file.

!! Copyright (C) 2005-2006 M. Marques, A. Castro, A. Rubio, G. Bertsch, X. Andrade
!!
!! 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_solver_oct_m
  use boundaries_oct_m
  use debug_oct_m
  use derivatives_oct_m
  use global_oct_m
  use, intrinsic :: iso_fortran_env
  use lalg_basic_oct_m
  use mesh_oct_m
  use mesh_function_oct_m
  use messages_oct_m
  use multigrid_oct_m
  use namespace_oct_m
  use nl_operator_oct_m
  use operate_f_oct_m
  use parser_oct_m
  use par_vec_oct_m
  use profiling_oct_m
  use space_oct_m
  use varinfo_oct_m
 
  implicit none
 
  integer, parameter ::       &
    GAUSS_SEIDEL        = 1,  &
    gauss_jacobi        = 2
 
  private
  public ::                   &
    multigrid_solver_cycle, &
    multigrid_solver_init,   &
    mg_solver_t
 
  type mg_solver_t
    private
 
    real(real64), public :: threshold
    real(real64) :: relax_factor
 
    integer, public :: maxcycles = 0
    integer :: presteps
    integer :: poststeps
    integer :: restriction_method
    integer :: relaxation_method
  end type mg_solver_t
 
contains
 
  ! ---------------------------------------------------------
  subroutine multigrid_solver_init(this, namespace, space, mesh, thr)
    type(mg_solver_t), intent(out)   :: this
    type(namespace_t), intent(in)    :: namespace
    class(space_t),    intent(in)    :: space
    type(mesh_t),      intent(inout) :: mesh
    real(real64),      intent(in)    :: thr
 
    push_sub(multigrid_solver_init)
 
    this%threshold = thr
 
    !%Variable MultigridPresmoothingSteps
    !%Type integer
    !%Default 1
    !%Section Hamiltonian::Poisson::Multigrid
    !%Description
    !% Number of Gauss-Seidel smoothing steps before coarse-level
    !% correction in the multigrid solver.
    !%End
    call parse_variable(namespace, 'MultigridPresmoothingSteps', 1, this%presteps)
 
    !%Variable MultigridPostsmoothingSteps
    !%Type integer
    !%Default 4
    !%Section Hamiltonian::Poisson::Multigrid
    !%Description
    !% Number of Gauss-Seidel smoothing steps after coarse-level
    !% correction in the multigrid solver.
    !%End
    call parse_variable(namespace, 'MultigridPostsmoothingSteps', 4, this%poststeps)
 
    !%Variable MultigridMaxCycles
    !%Type integer
    !%Default 50
    !%Section Hamiltonian::Poisson::Multigrid
    !%Description
    !% Maximum number of multigrid cycles that are performed if
    !% convergence is not achieved.
    !%End
    call parse_variable(namespace, 'MultigridMaxCycles', 50, this%maxcycles)
 
    !%Variable MultigridRestrictionMethod
    !%Type integer
    !%Default fullweight
    !%Section Hamiltonian::Poisson::Multigrid
    !%Description
    !% Method used from fine-to-coarse grid transfer.
    !%Option injection 1
    !% Injection
    !%Option fullweight 2
    !% Fullweight restriction
    !%End
    call parse_variable(namespace, 'MultigridRestrictionMethod', 2, this%restriction_method)
    if (.not. varinfo_valid_option('MultigridRestrictionMethod', this%restriction_method)) then
      call messages_input_error(namespace, 'MultigridRestrictionMethod')
    end if
    call messages_print_var_option("MultigridRestrictionMethod", this%restriction_method, namespace=namespace)
 
    !%Variable MultigridRelaxationMethod
    !%Type integer
    !%Section Hamiltonian::Poisson::Multigrid
    !%Description
    !% Method used to solve the linear system approximately in each grid for the
    !% multigrid procedure that solves a linear equation like the Poisson equation. Default is <tt>gauss_seidel</tt>,
    !% unless curvilinear coordinates are used, in which case the default is <tt>gauss_jacobi</tt>.
    !%Option gauss_seidel 1
    !% Gauss-Seidel.
    !%Option gauss_jacobi 2
    !% Gauss-Jacobi.
    !%End
    if (mesh%use_curvilinear) then
      call parse_variable(namespace, 'MultigridRelaxationMethod', gauss_jacobi, this%relaxation_method)
    else
      call parse_variable(namespace, 'MultigridRelaxationMethod', gauss_seidel, this%relaxation_method)
    end if
 
    if (.not. varinfo_valid_option('MultigridRelaxationMethod', this%relaxation_method)) then
      call messages_input_error(namespace, 'MultigridRelaxationMethod')
    end if
    call messages_print_var_option("MultigridRelaxationMethod", this%relaxation_method, namespace=namespace)
 
    !%Variable MultigridRelaxationFactor
    !%Type float
    !%Section Hamiltonian::Poisson::Multigrid
    !%Description
    !% Relaxation factor of the relaxation operator used for the
    !% multigrid method. This is mainly for debugging,
    !% since overrelaxation does not help in a multigrid scheme.
    !% The default is 1.0, except 0.6666 for the <tt>gauss_jacobi</tt> method.
    !%End
    if (this%relaxation_method == gauss_jacobi) then
      call parse_variable(namespace, 'MultigridRelaxationFactor', 0.6666_real64, this%relax_factor)
    else
      call parse_variable(namespace, 'MultigridRelaxationFactor', m_one, this%relax_factor)
    end if
 
    pop_sub(multigrid_solver_init)
  end subroutine multigrid_solver_init
 
 
  ! ---------------------------------------------------------
  recursive subroutine multigrid_solver_cycle(this, der, op, sol, rhs)
    type(mg_solver_t),           intent(in)    :: this
    type(derivatives_t),         intent(in)    :: der
    type(nl_operator_t),        intent(in)    :: op
    real(real64), contiguous,   intent(inout) :: sol(:)
    real(real64), contiguous,   intent(in)    :: rhs(:)
 
    integer :: ip, iter
    real(real64)   :: resnorm
    real(real64), allocatable :: residue(:), coarse_residue(:), correction(:), coarse_correction(:)
 
    push_sub(multigrid_solver_cycle)
 
    safe_allocate(residue(1:der%mesh%np_part))
 
    if (associated(der%coarser)) then
      assert(associated(op%coarser))
 
      safe_allocate(correction(1:der%mesh%np_part))
      safe_allocate(coarse_residue(1:der%coarser%mesh%np_part))
      safe_allocate(coarse_correction(1:der%coarser%mesh%np_part))
 
      if (debug%info) then
        message(1) = "Multigrid restriction"
        call messages_info(1)
      end if
      call multigrid_relax(this, der%mesh, der, op, sol, rhs, this%presteps)
 
      call get_residual()
 
      call dmultigrid_fine2coarse(der%to_coarser, der, der%coarser%mesh, residue, coarse_residue, this%restriction_method)
 
      coarse_correction = m_zero
      call multigrid_solver_cycle(this, der%coarser, op%coarser, coarse_correction, coarse_residue)
 
      if (debug%info) then
        message(1) = "Multigrid prolongation"
        call messages_info(1)
      end if
      correction = m_zero
      call dmultigrid_coarse2fine(der%to_coarser, der%coarser, der%mesh, coarse_correction, correction)
 
      call lalg_axpy(der%mesh%np, m_one, correction, sol)
 
      safe_deallocate_a(correction)
      safe_deallocate_a(coarse_residue)
      safe_deallocate_a(coarse_correction)
 
      call multigrid_relax(this, der%mesh, der, op, sol, rhs, this%poststeps)
 
    else ! Coarsest grid
 
 
      do iter = 1, this%maxcycles
 
        call multigrid_relax(this, der%mesh, der, op, sol, rhs, this%presteps + this%poststeps)
 
        call get_residual()
 
        resnorm = dmf_nrm2(der%mesh, residue)
        if (resnorm < this%threshold) exit
 
      end do
 
      if (debug%info) then
        write(message(1), '(a,i4,a)') "Multigrid coarsest grid solver converged in ", iter, " iterations."
        write(message(2), '(a,es18.6)') "Residue norm is ", resnorm
        call messages_info(2)
      end if
 
    end if
 
    safe_deallocate_a(residue)
 
    pop_sub(multigrid_solver_cycle)
 
  contains
 
    subroutine get_residual()
      ! Compute the residue
      call dderivatives_perform(op, der, sol, residue)
      !$omp parallel do
      do ip = 1, der%mesh%np
        residue(ip) = rhs(ip) - residue(ip)
      end do
    end subroutine get_residual
 
  end subroutine multigrid_solver_cycle
 
  ! ---------------------------------------------------------
  subroutine multigrid_relax(this, mesh, der, op, sol, rhs, steps)
    type(mg_solver_t),          intent(in)    :: this
    type(mesh_t),               intent(in)    :: mesh
    type(derivatives_t),        intent(in)    :: der
    type(nl_operator_t),        intent(in)    :: op
    real(real64), contiguous,   intent(inout) :: sol(:)
    real(real64), contiguous,   intent(in)    :: rhs(:)
    integer,                    intent(in)    :: steps
 
    integer :: istep
    integer :: ip, nn
    real(real64)   :: point, factor
    real(real64), allocatable :: op_sol(:), diag(:)
 
    push_sub(multigrid_relax)
    call profiling_in("MG_GAUSS_SEIDEL")
 
    select case (this%relaxation_method)
 
    case (gauss_seidel)
 
      do istep = 1, steps
 
        call boundaries_set(der%boundaries, der%mesh, sol)
 
        if (mesh%parallel_in_domains) then
          call dpar_vec_ghost_update(mesh%pv, sol)
        end if
 
        nn = op%stencil%size
 
        if (op%const_w) then
          factor = -m_one/op%w(op%stencil%center, 1)*this%relax_factor
          call dgauss_seidel(op%stencil%size, op%w(1, 1), op%nri, &
            op%ri(1, 1), op%rimap_inv(1), op%rimap_inv(2), factor, sol(1), rhs(1))
        else
          !$omp parallel do private(point)
          do ip = 1, mesh%np
            point = sum(op%w(1:nn, ip)*sol(op%index(1:nn, ip)))
            sol(ip) = sol(ip) - this%relax_factor/op%w(op%stencil%center, ip)*(point-rhs(ip))
          end do
        end if
 
      end do
      call profiling_count_operations(mesh%np*(steps + 1)*(2*nn + 3))
 
    case (gauss_jacobi)
 
      safe_allocate(op_sol(1:mesh%np))
      safe_allocate(diag(1:mesh%np))
 
      call dnl_operator_operate_diag(op, diag)
      !$omp parallel do
      do ip = 1, mesh%np
        diag(ip) = this%relax_factor/diag(ip)
      end do
 
      do istep = 1, steps
        call dderivatives_perform(op, der, sol, op_sol)
        !$omp parallel do
        do ip = 1, mesh%np
          sol(ip) = sol(ip) - diag(ip)*(op_sol(ip) - rhs(ip))
        end do
      end do
 
      safe_deallocate_a(diag)
      safe_deallocate_a(op_sol)
 
    end select
 
    call profiling_out("MG_GAUSS_SEIDEL")
    pop_sub(multigrid_relax)
 
  end subroutine multigrid_relax
 
end module multigrid_solver_oct_m
 
!! Local Variables:
!! mode: f90
!! coding: utf-8
!! End: