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 parser_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
    logical :: periodic = .false.
 
    integer :: multigrid_cycle
 
  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)
 
    this%periodic = space%is_periodic()
    if (.not. this%periodic) then
      call poisson_corrections_init(this%corrector, namespace, space, ml, mesh)
    end if
 
    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, mesh, thr)
 
    !%Variable PoissonMultigridCycle
    !%Type integer
    !%Section Hamiltonian::Poisson::Multigrid
    !%Description
    !% The flavor of multigrid cycle
    !%Option v_shape 1
    !% V-shape cycle
    !%Option w_shape 2
    !% W-shape cycle
    !%Option fmg 3
    !% Full multigrid solver
    !%End
    call parse_variable(namespace, 'PoissonMultigridCycle', mg_v_shape, this%multigrid_cycle)
 
    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)
 
    if (.not. this%periodic) 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(out) :: pot(:)
    real(real64), contiguous,    intent(in)    :: rho(:)
 
    integer :: ip
    real(real64), allocatable :: vh_correction(:), res(:), cor(:)
 
    push_sub(poisson_multigrid_solver)
 
    safe_allocate(res(1:der%mesh%np_part))
    safe_allocate(cor(1:der%mesh%np_part))
    cor = m_zero
 
    if (this%periodic) then
      res(1:der%mesh%np) = rho(1:der%mesh%np)
      call lalg_scal(der%mesh%np, -m_four*m_pi, res)
      call dmf_remove_average(der%mesh, res)
    else
      ! correction for treating boundaries
      safe_allocate(vh_correction(1:der%mesh%np_part))
      call correct_rho(this%corrector, der, rho, res, vh_correction)
      call lalg_scal(der%mesh%np, -m_four*m_pi, res)
    end if
 
    select case (this%multigrid_cycle)
    case(mg_v_shape, mg_w_shape)
      call multigrid_iterative_solver(this%mg_solver, namespace, this%mgrid%level(0)%der, this%mgrid%level(0)%der%lapl, &
        cor, res, this%multigrid_cycle)
    case(mg_fmg)
      call multigrid_fmg_solver(this%mg_solver, namespace, this%mgrid%level(0)%der, this%mgrid%level(0)%der%lapl, cor, res)
    case default
      assert(.false.)
    end select
 
    if (this%periodic) then
      call dmf_remove_average(der%mesh, cor)
      pot(1:der%mesh%np) = cor(1:der%mesh%np)
    else
      do ip = 1, der%mesh%np
        pot(ip) = cor(ip) + vh_correction(ip)
      end do
      safe_deallocate_a(vh_correction)
    end if
    safe_deallocate_a(res)
    safe_deallocate_a(cor)
 
    pop_sub(poisson_multigrid_solver)
  end subroutine poisson_multigrid_solver
 
end module poisson_multigrid_oct_m
 
!! Local Variables:
!! mode: f90
!! coding: utf-8
!! End: