preconditioners.F90

Go to the documentation of this file.

!! Copyright (C) 2006 X. Andrade, 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 preconditioners_oct_m
  use batch_oct_m
  use batch_ops_oct_m
  use boundaries_oct_m
  use debug_oct_m
  use derivatives_oct_m
  use global_oct_m
  use grid_oct_m
  use hamiltonian_elec_oct_m
  use hamiltonian_elec_base_oct_m
  use, intrinsic :: iso_fortran_env
  use lalg_basic_oct_m
  use mesh_oct_m
  use messages_oct_m
  use multicomm_oct_m
  use multigrid_oct_m
  use namespace_oct_m
  use nl_operator_oct_m
  use parser_oct_m
  use poisson_oct_m
  use profiling_oct_m
  use stencil_star_oct_m
  use space_oct_m
  use varinfo_oct_m
  use wfs_elec_oct_m
 
  implicit none
  private
 
  integer, public, parameter ::     &
    PRE_NONE      = 0,              &
    pre_filter    = 1,              &
    pre_jacobi    = 2,              &
    pre_poisson   = 3,              &
    pre_multigrid = 7
 
  public ::                            &
    preconditioner_t,                  &
    preconditioner_init,               &
    preconditioner_end,                &
    dpreconditioner_apply,             &
    zpreconditioner_apply,             &
    dpreconditioner_apply_batch,       &
    zpreconditioner_apply_batch,       &
    preconditioner_obsolete_variables
 
  type preconditioner_t
    private
    integer :: which = pre_none
 
    type(nl_operator_t), allocatable :: op_array(:)
    type(nl_operator_t), pointer, public :: op
 
    real(real64), allocatable  :: diag_lapl(:)
    integer             :: npre, npost, nmiddle
 
    type(multigrid_t) :: mgrid  ! multigrid object
 
    type(derivatives_t), pointer :: der => null()
  end type preconditioner_t
 
contains
 
  ! ---------------------------------------------------------
  subroutine preconditioner_init(this, namespace, gr, mc, space)
    type(preconditioner_t), target, intent(out)    :: this
    type(namespace_t),              intent(in)     :: namespace
    type(grid_t), target,           intent(in)     :: gr
    type(multicomm_t),              intent(in)     :: mc
    class(space_t),                 intent(in)     :: space
 
    real(real64) :: alpha, default_alpha, omega
    real(real64) :: vol
    integer :: default
    integer :: maxp, is, ns, ip, ip2
    character(len=32) :: name
 
    push_sub(preconditioner_init)
 
    safe_allocate(this%op_array(1))
    this%op => this%op_array(1)
 
    this%der => gr%der
 
    !%Variable Preconditioner
    !%Type integer
    !%Section SCF::Eigensolver
    !%Description
    !% Which preconditioner to use in order to solve the Kohn-Sham
    !% equations or the linear-response equations. The default is
    !% pre_filter, except for curvilinear coordinates, where no
    !% preconditioner is applied by default.
    !%Option no 0
    !% Do not apply preconditioner.
    !%Option pre_filter 1
    !% Filter preconditioner.
    !%Option pre_jacobi 2
    !% Jacobi preconditioner. Only the local part of the pseudopotential is used.
    !% Not very helpful.
    !%Option pre_poisson 3
    !% Uses the full Laplacian as preconditioner. The inverse is calculated through
    !% the solution of the Poisson equation. This is, of course, very slow.
    !%Option pre_multigrid 7
    !% Multigrid preconditioner.
    !%End
 
    if (gr%use_curvilinear) then
      default = pre_none
    else
      default = pre_filter
    end if
 
    call parse_variable(namespace, 'Preconditioner', default, this%which)
    if (.not. varinfo_valid_option('Preconditioner', this%which)) call messages_input_error(namespace, 'Preconditioner')
    call messages_print_var_option('Preconditioner', this%which, namespace=namespace)
 
    select case (this%which)
    case (pre_filter)
      ! the smoothing is performed uing the same stencil as the Laplacian
      name = "Preconditioner"
      call derivatives_get_lapl(gr%der, namespace, this%op_array, space, name, 1)
 
      !%Variable PreconditionerFilterFactor
      !%Type float
      !%Section SCF::Eigensolver
      !%Description
      !% This variable controls how much filter preconditioner is
      !% applied. A value of 1.0 means no preconditioning, 0.5 is the
      !% standard.
      !%
      !% The default is 0.5, except for periodic systems where the
      !% default is 0.6.
      !%
      !% If you observe that the first eigenvectors are not converging
      !% properly, especially for periodic systems, you should
      !% increment this value.
      !%
      !% The allowed range for this parameter is between 0.5 and 1.0.
      !% For other values, the SCF may converge to wrong results.
      !%End
      default_alpha = m_half
      if (space%is_periodic()) default_alpha = 0.6_real64
 
      call parse_variable(namespace, 'PreconditionerFilterFactor', default_alpha, alpha)
 
      call messages_print_var_value('PreconditionerFilterFactor', alpha, namespace=namespace)
 
      ! check for correct interval of alpha
      if (alpha < m_half .or. alpha > m_one) then
        call messages_input_error(namespace, 'PreconditionerFilterFactor')
      end if
 
      ns = this%op%stencil%size
 
      if (this%op%const_w) then
        maxp = 1
      else
        maxp = gr%np
      end if
 
      !We change the weights to be the one of the kinetic energy operator
      this%op%w = -m_half * this%op%w
 
      safe_allocate(this%diag_lapl(1:this%op%np))
      call dnl_operator_operate_diag(this%op, this%diag_lapl)
 
      do ip = 1, maxp
 
        if (gr%use_curvilinear) vol = sum(gr%vol_pp(ip + this%op%ri(1:ns, this%op%rimap(ip))))
 
        !The filter preconditioner is given by two iterations of the Relaxation Jacobi method
        !This leads to \tilde{\psi} =  \omega D^{-1}(2\psi - \omega D^{-1} (-0.5\Laplacian) \psi),
        ! where \omega = 2(1-\alpha), and D is the diagonal part of (-0.5\Laplacian)
        !In order to have this to work in all cases, such as different spacings, nonorthogonal cells, ...
        !We directly apply this formula to renormalize the weights of the Laplacian
        !and get the correct preconditioner
 
        omega = m_two * (m_one-alpha)
 
        do is = 1, ns
          this%op%w(is, ip) = - omega / this%diag_lapl(ip) * this%op%w(is, ip)
          if (is == this%op%stencil%center) then
            this%op%w(is, ip) = this%op%w(is, ip) + m_two
          end if
          this%op%w(is, ip) = this%op%w(is, ip) * omega / this%diag_lapl(ip)
 
          if (gr%use_curvilinear) then
            ip2 = ip + this%op%ri(is, this%op%rimap(ip))
            this%op%w(is, ip) = this%op%w(is, ip)*(ns*gr%vol_pp(ip2)/vol)
          end if
        end do
      end do
 
      safe_deallocate_a(this%diag_lapl)
 
      call nl_operator_update_gpu_buffers(this%op)
      call nl_operator_output_weights(this%op)
 
    case (pre_jacobi, pre_multigrid)
      safe_allocate(this%diag_lapl(1:gr%np))
      call derivatives_lapl_diag(gr%der, this%diag_lapl)
      call lalg_scal(gr%np, -m_half, this%diag_lapl(:))
    end select
 
    if (this%which == pre_multigrid) then
      !%Variable PreconditionerIterationsPre
      !%Type integer
      !%Section SCF::Eigensolver
      !%Description
      !% This variable is the number of pre-smoothing iterations for the multigrid
      !% preconditioner. The default is 1.
      !%End
      call parse_variable(namespace, 'PreconditionerIterationsPre', 1, this%npre)
 
      !%Variable PreconditionerIterationsMiddle
      !%Type integer
      !%Section SCF::Eigensolver
      !%Description
      !% This variable is the number of smoothing iterations on the coarsest grid for the multigrid
      !% preconditioner. The default is 1.
      !%End
      call parse_variable(namespace, 'PreconditionerIterationsMiddle', 1, this%nmiddle)
 
      !%Variable PreconditionerIterationsPost
      !%Type integer
      !%Section SCF::Eigensolver
      !%Description
      !% This variable is the number of post-smoothing iterations for the multigrid
      !% preconditioner. The default is 2.
      !%End
      call parse_variable(namespace, 'PreconditionerIterationsPost', 2, this%npost)
 
      call multigrid_init(this%mgrid, namespace, space, gr, gr%der, gr%stencil, mc, nlevels=3)
    end if
 
    pop_sub(preconditioner_init)
  end subroutine preconditioner_init
 
  ! ---------------------------------------------------------
  subroutine preconditioner_end(this)
    type(preconditioner_t), intent(inout) :: this
 
    push_sub(preconditioner_end)
 
    select case (this%which)
    case (pre_filter)
      call nl_operator_end(this%op)
 
    case (pre_jacobi)
      safe_deallocate_a(this%diag_lapl)
 
    case (pre_multigrid)
      safe_deallocate_a(this%diag_lapl)
      call multigrid_end(this%mgrid)
 
    end select
 
    nullify(this%op)
    safe_deallocate_a(this%op_array)
 
    pop_sub(preconditioner_end)
  end subroutine preconditioner_end
 
  ! ---------------------------------------------------------
  subroutine preconditioner_obsolete_variables(namespace, old_prefix, new_prefix)
    type(namespace_t),   intent(in)    :: namespace
    character(len=*),    intent(in)    :: old_prefix
    character(len=*),    intent(in)    :: new_prefix
 
    call messages_obsolete_variable(namespace, trim(old_prefix)//'Preconditioner', trim(new_prefix)//'Preconditioner')
  end subroutine preconditioner_obsolete_variables
 
#include "undef.F90"
#include "complex.F90"
#include "preconditioners_inc.F90"
 
#include "undef.F90"
#include "real.F90"
#include "preconditioners_inc.F90"
 
end module preconditioners_oct_m
 
!! Local Variables:
!! mode: f90
!! coding: utf-8
!! End: