poisson_cg.F90

Go to the documentation of this file.

!! Cropyright (C) 2004-2006 M. Marques, A. Castro, A. Rubio, G. Bertsch
!!
!! 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_cg_oct_m
  use debug_oct_m
  use derivatives_oct_m
  use global_oct_m
  use iso_c_binding, only: c_loc
  use, intrinsic :: iso_fortran_env
  use lalg_basic_oct_m
  use lapl_operator_oct_m
  use mesh_oct_m
  use mesh_function_oct_m
  use messages_oct_m
  use namespace_oct_m
  use poisson_corrections_oct_m
  use profiling_oct_m
  use solvers_oct_m
 
  implicit none
 
  private
  public ::           &
    poisson_cg_init,  &
    poisson_cg_end,   &
    poisson_cg1,      &
    poisson_cg2
 
  real(real64), public :: threshold
  integer, public :: maxiter
 
contains
 
 
  ! ---------------------------------------------------------
  subroutine poisson_cg_init(thr, itr)
    integer,      intent(in) :: itr
    real(real64), intent(in) :: thr
 
    push_sub(poisson_cg_init)
    threshold = thr
    maxiter = itr
    pop_sub(poisson_cg_init)
  end subroutine poisson_cg_init
 
 
  ! ---------------------------------------------------------
  subroutine poisson_cg_end
 
  end subroutine poisson_cg_end
 
 
  ! ---------------------------------------------------------
  subroutine poisson_cg1(namespace, der, corrector, pot, rho)
    type(namespace_t),            intent(in)    :: namespace
    type(derivatives_t),  target, intent(in)    :: der
    type(poisson_corr_t),         intent(in)    :: corrector
    real(real64), contiguous,     intent(inout) :: pot(:)
    real(real64),                 intent(in)    :: rho(:)
 
    integer :: iter, ip
    real(real64) :: res
    real(real64), allocatable :: wk(:), lwk(:), rhs(:), sol(:)
 
    push_sub(poisson_cg1)
 
    safe_allocate( wk(1:der%mesh%np_part))
    safe_allocate(lwk(1:der%mesh%np_part))
    safe_allocate(rhs(1:der%mesh%np_part))
    safe_allocate(sol(1:der%mesh%np_part))
 
    ! build initial guess for the potential
    wk(1:der%mesh%np) = pot(1:der%mesh%np)
    if (der%periodic_dim > 0) then
      call dderivatives_lapl(der, wk, lwk, ghost_update = .true.)
    else
      call poisson_boundary_conditions(corrector, der%mesh, rho, wk)
      call dderivatives_lapl(der, wk, lwk, ghost_update = .true., set_bc = .false.)
    end if
 
    ! Solving -\Delta = rhs, as CG require SDP operator
    !$omp parallel do
    do ip = 1, der%mesh%np
      rhs(ip) = (m_four*m_pi*rho(ip) + lwk(ip))
    end do
    !$omp end parallel do
    if (der%periodic_dim > 0) call dmf_remove_average(der%mesh, rhs)
    safe_deallocate_a(wk)
    safe_deallocate_a(lwk) ! they are no longer needed
 
    der_pointer  => der
    mesh_pointer => der%mesh
    mesh_aux => der%mesh
    call lalg_copy(der%mesh%np, rhs, sol)
    iter = maxiter
    if (der%periodic_dim > 0) then
      call dconjugate_gradients(der%mesh%np, sol, rhs, dlaplacian_op, internal_dotp, iter, res, &
        threshold, userdata=[c_loc(der)], nullspace=dmf_remove_average_aux)
    else
      call dconjugate_gradients(der%mesh%np, sol, rhs, dlaplacian_op, internal_dotp, iter, res, threshold, userdata=[c_loc(der)])
    end if
    if (res >= threshold) then
      message(1) = 'Conjugate-gradients Poisson solver did not converge.'
      write(message(2), '(a,i8)')    '  Iter = ',iter
      write(message(3), '(a,e14.6)') '  Res = ', res
      call messages_warning(3, namespace=namespace)
    end if
    nullify(der_pointer, mesh_pointer)
    call lalg_axpy(der%mesh%np, m_one, sol, pot)
 
    safe_deallocate_a(rhs)
    safe_deallocate_a(sol)
    pop_sub(poisson_cg1)
  end subroutine poisson_cg1
 
 
  ! ---------------------------------------------------------
  subroutine poisson_cg2(namespace, der, pot, rho)
    type(namespace_t),           intent(in)    :: namespace
    type(derivatives_t), target, intent(in)    :: der
    real(real64), contiguous,    intent(inout) :: pot(:)
    real(real64),                intent(in)    :: rho(:)
 
    integer :: iter, ip
    real(real64), allocatable :: potc(:), rhs(:)
    real(real64) :: res
 
    push_sub(poisson_cg2)
 
    iter = maxiter
    der_pointer  => der
    mesh_pointer => der%mesh
 
    safe_allocate(rhs(1:der%mesh%np))
    safe_allocate(potc(1:der%mesh%np_part))
 
    ! Solving -\Delta = rhs, as CG require SDP operator
    do ip = 1, der%mesh%np
      rhs(ip) = 4.0_real64*m_pi*rho(ip)
    end do
    call lalg_copy(der%mesh%np, pot, potc)
 
    call dconjugate_gradients(der%mesh%np, potc, rhs, dlaplacian_op, internal_dotp, iter, res, threshold, userdata=[c_loc(der)])
 
    if (res >= threshold) then
      message(1) = 'Conjugate-gradients Poisson solver did not converge.'
      write(message(2), '(a,i8)')    '  Iter = ', iter
      write(message(3), '(a,e14.6)') '  Res = ', res
      call messages_warning(3, namespace=namespace)
    end if
 
    call lalg_copy(der%mesh%np, potc, pot)
 
    nullify(der_pointer, mesh_pointer)
 
    safe_deallocate_a(rhs)
    safe_deallocate_a(potc)
 
    pop_sub(poisson_cg2)
  end subroutine poisson_cg2
 
end module poisson_cg_oct_m
 
!! Local Variables:
!! mode: f90
!! coding: utf-8
!! End: