poisson_corrections.F90

Go to the documentation of this file.

!! Copyright (C) 2002-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_corrections_oct_m
  use debug_oct_m
  use derivatives_oct_m
  use global_oct_m
  use iso_c_binding
  use lalg_basic_oct_m
  use loct_math_oct_m
  use math_oct_m
  use mesh_function_oct_m
  use mesh_oct_m
  use messages_oct_m
  use namespace_oct_m
  use nl_operator_oct_m
  use parser_oct_m
  use profiling_oct_m
  use space_oct_m
 
  implicit none
 
  private
  public ::                   &
    poisson_corrections_init, &
    poisson_corrections_end,  &
    correct_rho,              &
    poisson_boundary_conditions, &
    poisson_corr_t,           &
    internal_dotp,            &
    der_pointer,              &
    mesh_pointer
 
  real(real64), parameter :: alpha_ = m_five
 
  integer, parameter  ::    &
    CORR_NONE      = 0,     &
    corr_multipole = 1,     &
    corr_exact     = 3
 
  type poisson_corr_t
    private
    integer :: method = corr_none
    integer :: maxl
    real(real64), allocatable :: phi(:, :)
    real(real64), allocatable :: aux(:, :)
    real(real64), allocatable :: gaussian(:)
  end type poisson_corr_t
 
  type(derivatives_t), pointer :: der_pointer
  type(mesh_t),        pointer :: mesh_pointer
 
contains
 
  ! ---------------------------------------------------------
  subroutine poisson_corrections_init(this, namespace, space, ml, mesh)
    type(poisson_corr_t), intent(out) :: this
    type(namespace_t),    intent(in)  :: namespace
    class(space_t),       intent(in)  :: space
    integer,              intent(in)  :: ml
    type(mesh_t),         intent(in)  :: mesh
 
    real(real64) :: alpha, beta, ylm, rr, xx(space%dim)
    integer :: ip, ll, add_lm, lldfac, jj, mm
 
    push_sub(poisson_corrections_init)
 
    assert(space%dim == 3)
 
    if (space%is_periodic()) then
      call messages_not_implemented("Poisson boundary corrections for periodic systems", namespace=namespace)
    end if
 
    !%Variable PoissonSolverBoundaries
    !%Type integer
    !%Section Hamiltonian::Poisson
    !%Default multipole
    !%Description
    !% For finite systems, some Poisson solvers (<tt>multigrid</tt>,
    !% <tt>cg_corrected</tt>, and <tt>fft</tt> with <tt>PoissonFFTKernel = multipole_correction</tt>)
    !% require the calculation of the
    !% boundary conditions with an auxiliary method. This variable selects that method.
    !%Option multipole 1
    !% A multipole expansion of the density is used to approximate the potential on the boundaries.
    !%Option exact 3
    !% An exact integration of the Poisson equation is done over the boundaries. This option is
    !% experimental, and not implemented for domain parallelization.
    !%End
    call parse_variable(namespace, 'PoissonSolverBoundaries', corr_multipole, this%method)
 
    select case (this%method)
    case (corr_multipole)
      this%maxl = ml
 
      add_lm = 0
      do ll = 0, this%maxl
        do mm = -ll, ll
          add_lm = add_lm + 1
        end do
      end do
 
      safe_allocate(this%phi(1:mesh%np, 1:add_lm))
      safe_allocate(this%aux(1:mesh%np, 1:add_lm))
      safe_allocate(this%gaussian(1:mesh%np))
 
      alpha = alpha_ * mesh%spacing(1)
      do ip = 1, mesh%np
        call mesh_r(mesh, ip, rr, coords = xx)
        this%gaussian(ip) = exp(-(rr/alpha)**2)
        add_lm = 1
        do ll = 0, this%maxl
          lldfac = 1
          do jj = 1, 2*ll+1, 2
            lldfac = lldfac * jj
          end do
          beta = sqrt(m_pi)*2**(ll+3) / lldfac
          do mm = -ll, ll
            call ylmr_real(xx, ll, mm, ylm)
            if (rr > m_epsilon) then
              this%phi(ip, add_lm) = beta*isubl(ll, rr/alpha)*ylm/rr**(ll+1)
              this%aux(ip, add_lm) = rr**ll*ylm
            else
              this%phi(ip, add_lm) = beta*ylm / alpha
              if (ll == 0) then
                this%aux(ip, add_lm) = ylm
              else
                this%aux(ip, add_lm) = m_zero
              end if
            end if
            add_lm = add_lm + 1
          end do
        end do
      end do
 
    case (corr_exact)
      call messages_experimental('Exact Poisson solver boundaries', namespace=namespace)
      if (mesh%parallel_in_domains) call messages_not_implemented('Exact Poisson solver boundaries with domain parallelization')
 
    end select
 
    pop_sub(poisson_corrections_init)
 
  contains
 
    ! ---------------------------------------------------------
    real(real64) function isubl( ll, xx)
      integer, intent(in) :: ll
      real(real64),   intent(in) :: xx
 
      ! no push_sub, called too frequently
      isubl = m_half * gamma((ll + m_half)) * (m_one - loct_incomplete_gamma(ll+m_half, xx**2))
 
    end function isubl
 
  end subroutine poisson_corrections_init
 
 
  ! ---------------------------------------------------------
  subroutine poisson_corrections_end(this)
    type(poisson_corr_t), intent(inout) :: this
 
    push_sub(poisson_corrections_end)
 
    select case (this%method)
    case (corr_multipole)
      safe_deallocate_a(this%phi)
      safe_deallocate_a(this%aux)
      safe_deallocate_a(this%gaussian)
    case (corr_exact)
    end select
 
    pop_sub(poisson_corrections_end)
  end subroutine poisson_corrections_end
 
 
  ! ---------------------------------------------------------
  subroutine correct_rho(this, der, rho, rho_corrected, vh_correction)
    type(poisson_corr_t),     intent(in)  :: this
    type(derivatives_t),      intent(in)  :: der
    real(real64), contiguous, intent(in)  :: rho(:)
    real(real64), contiguous, intent(out) :: rho_corrected(:)
    real(real64), contiguous, intent(out) :: vh_correction(:)
 
    integer :: ip, add_lm, ll, mm, lldfac, jj, ip2
    real(real64)   :: alpha, vv, rr
    real(real64), allocatable :: mult(:)
    real(real64), allocatable :: betal(:)
 
    push_sub(correct_rho)
    call profiling_in("POISSON_CORRECT")
 
    assert(ubound(vh_correction, dim = 1) == der%mesh%np_part)
 
    select case (this%method)
    case (corr_multipole)
 
      safe_allocate(mult(1:(this%maxl+1)**2))
      call get_multipoles(this, der%mesh, rho, this%maxl, mult)
 
      alpha = alpha_ * der%mesh%spacing(1)
 
      safe_allocate(betal(1:(this%maxl+1)**2))
      add_lm = 1
      do ll = 0, this%maxl
        do mm = -ll, ll
          lldfac = 1
          do jj = 1, 2*ll+1, 2
            lldfac = lldfac*jj
          end do
          betal(add_lm) = (2**(ll + 2)) / (alpha**(2*ll+3) * sqrt(m_pi) * lldfac)
          add_lm = add_lm + 1
        end do
      end do
 
      call lalg_copy(der%mesh%np, rho, rho_corrected)
      vh_correction = m_zero
      add_lm = 1
      do ll = 0, this%maxl
        do mm = -ll, ll
          do ip = 1, der%mesh%np
            rho_corrected(ip) = rho_corrected(ip) - mult(add_lm)*betal(add_lm)*this%aux(ip, add_lm)*this%gaussian(ip)
          end do
          call lalg_axpy(der%mesh%np, mult(add_lm), this%phi(:, add_lm), vh_correction)
          add_lm = add_lm + 1
        end do
      end do
 
      safe_deallocate_a(mult)
      safe_deallocate_a(betal)
 
    case (corr_exact)
 
      do ip = 1, der%mesh%np
        vh_correction(ip) = m_zero
      end do
 
      do ip = der%mesh%np + 1, der%mesh%np_part
        vv = m_zero
        do ip2 = 1, der%mesh%np
          rr = norm2((der%mesh%x(1:der%dim, ip) - der%mesh%x(1:der%dim, ip2)))
          vv = vv + rho(ip2)/rr
        end do
        vh_correction(ip) = der%mesh%volume_element*vv
      end do
 
      call dderivatives_lapl(der, vh_correction, rho_corrected, set_bc = .false.)
 
      do ip = 1, der%mesh%np
        rho_corrected(ip) = rho(ip) + m_one/(m_four*m_pi)*rho_corrected(ip)
      end do
 
    end select
 
    call profiling_out("POISSON_CORRECT")
    pop_sub(correct_rho)
  end subroutine correct_rho
 
 
  ! ---------------------------------------------------------
  subroutine get_multipoles(this, mesh, rho, ml, mult)
    type(poisson_corr_t), intent(in)  :: this
    type(mesh_t),         intent(in)  :: mesh
    real(real64),         intent(in)  :: rho(:)
    integer,              intent(in)  :: ml
    real(real64),         intent(out) :: mult((ml+1)**2)
 
    integer :: add_lm, ll, mm
 
    push_sub(get_multipoles)
 
    mult(:) = m_zero
    add_lm = 1
    do ll = 0, ml
      do mm = -ll, ll
        mult(add_lm) = dmf_dotp(mesh, rho, this%aux(:, add_lm))
        add_lm = add_lm + 1
      end do
    end do
 
    pop_sub(get_multipoles)
  end subroutine get_multipoles
 
 
  ! ---------------------------------------------------------
  real(real64) function internal_dotp(xx, yy) result(res)
    real(real64), intent(in)    :: xx(:)
    real(real64), intent(in)    :: yy(:)
 
    push_sub(internal_dotp)
 
    res = dmf_dotp(mesh_pointer, xx, yy)
    pop_sub(internal_dotp)
  end function internal_dotp
 
 
  ! ---------------------------------------------------------
  subroutine poisson_boundary_conditions(this, mesh, rho, pot)
    type(poisson_corr_t), intent(in)    :: this
    type(mesh_t),         intent(in)    :: mesh
    real(real64),         intent(in)    :: rho(:)
    real(real64),         intent(inout) :: pot(:)
 
    integer :: ip, add_lm, ll, mm, bp_lower
    real(real64)   :: xx(mesh%box%dim), rr, s1, sa
    real(real64), allocatable :: mult(:)
 
    push_sub(poisson_boundary_conditions)
 
    assert(mesh%box%dim == 3)
 
    safe_allocate(mult(1:(this%maxl+1)**2))
 
    call get_multipoles(this, mesh, rho, this%maxl, mult)
 
    bp_lower = mesh%np + 1
    if (mesh%parallel_in_domains) then
      bp_lower = bp_lower + mesh%pv%np_ghost
    end if
 
    pot(bp_lower:mesh%np_part) = m_zero
    do ip = bp_lower, mesh%np_part ! boundary conditions
      call mesh_r(mesh, ip, rr, coords = xx)
      add_lm = 1
      do ll = 0, this%maxl
        s1 = m_four*m_pi/((m_two*ll + m_one)*rr**(ll + 1))
        do mm = -ll, ll
          call ylmr_real(xx, ll, mm, sa)
          pot(ip) = pot(ip) + sa * mult(add_lm) * s1
          add_lm = add_lm+1
        end do
      end do
    end do
 
    safe_deallocate_a(mult)
    pop_sub(poisson_boundary_conditions)
  end subroutine poisson_boundary_conditions
 
end module poisson_corrections_oct_m
 
!! Local Variables:
!! mode: f90
!! coding: utf-8
!! End: