poisson_multigrid.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 poisson_multigrid_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 multicomm_oct_m
  use multigrid_oct_m
  use multigrid_solver_oct_m
  use namespace_oct_m
  use poisson_corrections_oct_m
  use profiling_oct_m
  use space_oct_m
  use stencil_oct_m
 
  implicit none
 
  private
  public ::                   &
    poisson_multigrid_solver, &
    poisson_multigrid_init,   &
    poisson_multigrid_end,    &
    poisson_mg_solver_t
 
  type poisson_mg_solver_t
    private
 
    type(poisson_corr_t) :: corrector
    type(multigrid_t) :: mgrid  ! multigrid object
    type(mg_solver_t) :: mg_solver
 
  end type poisson_mg_solver_t
 
contains
 
  ! ---------------------------------------------------------
  subroutine poisson_multigrid_init(this, namespace, space, mesh, der, stencil, mc, ml, thr)
    type(poisson_mg_solver_t), intent(out)   :: this
    type(namespace_t),         intent(in)    :: namespace
    class(space_t),            intent(in)    :: space
    type(mesh_t),              intent(inout) :: mesh
    type(derivatives_t),       intent(in)    :: der
    type(stencil_t),           intent(in)    :: stencil
    type(multicomm_t),         intent(in)    :: mc
    integer,                   intent(in)    :: ml
    real(real64),              intent(in)    :: thr
 
    integer :: i
 
    push_sub(poisson_multigrid_init)
 
    call poisson_corrections_init(this%corrector, namespace, space, ml, mesh)
 
    call multigrid_init(this%mgrid, namespace, space, mesh, der, stencil, mc)
 
    ! For the multigrid solver to work, we need to set a pointer from one operator
    ! to the corresponding one on the coarser grid.
    do i = 1, this%mgrid%n_levels
      this%mgrid%level(i - 1)%der%lapl%coarser => this%mgrid%level(i)%der%lapl
    end do
 
    call multigrid_solver_init(this%mg_solver, namespace, space, mesh, thr)
 
    pop_sub(poisson_multigrid_init)
  end subroutine poisson_multigrid_init
 
 
  ! ---------------------------------------------------------
  subroutine poisson_multigrid_end(this)
    type(poisson_mg_solver_t), intent(inout) :: this
 
    push_sub(poisson_multigrid_end)
 
    call poisson_corrections_end(this%corrector)
 
    call multigrid_end(this%mgrid)
 
    pop_sub(poisson_multigrid_end)
  end subroutine poisson_multigrid_end
 
 
  ! ---------------------------------------------------------
  subroutine poisson_multigrid_solver(this, namespace, der, pot, rho)
    type(poisson_mg_solver_t),   intent(in)    :: this
    type(namespace_t),           intent(in)    :: namespace
    type(derivatives_t),         intent(in)    :: der
    real(real64),                intent(inout) :: pot(:)
    real(real64), contiguous,    intent(in)    :: rho(:)
 
    integer :: iter, ip
    real(real64) :: resnorm
    real(real64), allocatable :: vh_correction(:), res(:), cor(:), err(:)
 
    push_sub(poisson_multigrid_solver)
 
    ! correction for treating boundaries
    safe_allocate(vh_correction(1:der%mesh%np_part))
    safe_allocate(res(1:der%mesh%np))
    safe_allocate(cor(1:der%mesh%np_part))
    safe_allocate(err(1:der%mesh%np))
 
    call correct_rho(this%corrector, der, rho, res, vh_correction)
    call lalg_scal(der%mesh%np, -m_four*m_pi, res)
 
    do ip = 1, der%mesh%np
      cor(ip) = pot(ip) - vh_correction(ip)
    end do
 
    do iter = 1, this%mg_solver%maxcycles
 
      call multigrid_solver_cycle(this%mg_solver, this%mgrid%level(0)%der, this%mgrid%level(0)%der%lapl, cor, res)
      call dderivatives_lapl(der, cor, err)
      do ip = 1, der%mesh%np
        err(ip) = res(ip) - err(ip)
      end do
      resnorm =  dmf_nrm2(der%mesh, err)
 
      if (resnorm < this%mg_solver%threshold) exit
 
      if (debug%info) then
        write(message(1), '(a,i5,a,e13.6)') "Multigrid: base level: iter ", iter, " res ", resnorm
        call messages_info(1, namespace=namespace)
      end if
 
    end do
 
    if (resnorm >= this%mg_solver%threshold) then
      message(1) = 'Multigrid Poisson solver did not converge.'
      write(message(2), '(a,e14.6)') '  Norm of the residue = ', resnorm
      call messages_warning(2, namespace=namespace)
    else
      if (debug%info) then
        write(message(1), '(a,i4,a)') "Multigrid Poisson solver converged in ", iter, " iterations."
        write(message(2), '(a,e14.6)') '  Norm of the residue = ', resnorm
        call messages_info(2, namespace=namespace)
      end if
    end if
 
    do ip = 1, der%mesh%np
      pot(ip) = cor(ip) + vh_correction(ip)
    end do
 
    safe_deallocate_a(vh_correction)
    safe_deallocate_a(res)
    safe_deallocate_a(cor)
    safe_deallocate_a(err)
 
    pop_sub(poisson_multigrid_solver)
  end subroutine poisson_multigrid_solver
 
end module poisson_multigrid_oct_m
 
!! Local Variables:
!! mode: f90
!! coding: utf-8
!! End: