curv_modine.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 curv_modine_oct_m
  use coordinate_system_oct_m
  use debug_oct_m
  use global_oct_m
  use, intrinsic :: iso_fortran_env
  use lalg_adv_oct_m
  use messages_oct_m
  use namespace_oct_m
  use parser_oct_m
  use profiling_oct_m
  use root_solver_oct_m
  use unit_oct_m
  use unit_system_oct_m
 
  implicit none
 
  private
  public ::                     &
    curv_modine_t,              &
    curv_modine_copy
 
  type, extends(coordinate_system_t) :: curv_modine_t
    private
    real(real64), allocatable :: lsize(:)
    real(real64)       :: xbar
    real(real64)       :: Jbar
 
    real(real64), allocatable :: Jlocal(:)
    real(real64), allocatable :: Jrange(:)
 
    real(real64), allocatable :: chi_atoms(:,:)
    real(real64), allocatable :: csi(:,:)
 
    integer :: npos
    type(root_solver_t) :: rs
  contains
    procedure :: to_cartesian => curv_modine_to_cartesian
    procedure :: from_cartesian => curv_modine_from_cartesian
    procedure :: jacobian => curv_modine_jacobian
    procedure :: jacobian_inverse => curv_modine_jacobian_inverse
    procedure :: write_info => curv_modine_write_info
    procedure :: surface_element => curv_modine_surface_element
    final :: curv_modine_finalize
  end type curv_modine_t
 
  interface curv_modine_t
    procedure curv_modine_constructor
  end interface curv_modine_t
 
  ! Auxiliary variables for the root solver.
  class(curv_modine_t), pointer :: modine_p
  real(real64), allocatable :: x_p(:)
 
contains
 
  ! ---------------------------------------------------------
  function curv_modine_constructor(namespace, dim, npos, pos, lsize, spacing) result(modine)
    type(namespace_t), intent(in)  :: namespace
    integer,           intent(in)  :: dim
    integer,           intent(in)  :: npos
    real(real64),      intent(in)  :: pos(1:dim,1:npos)
    real(real64),      intent(in)  :: lsize(1:dim)
    real(real64),      intent(in)  :: spacing(1:dim)
    class(curv_modine_t), pointer :: modine
 
    push_sub(curv_modine_constructor)
 
    safe_allocate(modine)
 
    modine%dim = dim
    modine%local_basis = .true.
    modine%orthogonal = .true. ! This needs to be checked.
 
    modine%npos = npos
 
    !%Variable CurvModineXBar
    !%Type float
    !%Default 1/3
    !%Section Mesh::Curvilinear::Modine
    !%Description
    !% Size of central flat region (in units of <tt>Lsize</tt>). Must be between 0 and 1.
    !% See N. A. Modine, G. Zumbach, and E. Kaxiras, <i>Phys. Rev. B</i> <b>55</b>, 10289-10301 (1997).
    !%End
    call parse_variable(namespace, 'CurvModineXBar', m_one/m_three, modine%xbar)
 
    !%Variable CurvModineJBar
    !%Type float
    !%Default 1/2
    !%Section Mesh::Curvilinear::Modine
    !%Description
    !% Increase in density of points is inverse of this parameter.
    !% See N. A. Modine, G. Zumbach, and E. Kaxiras, <i>Phys. Rev. B</i> <b>55</b>, 10289-10301 (1997).
    !%End
    call parse_variable(namespace, 'CurvModineJBar', m_half, modine%Jbar)
 
    safe_allocate(modine%lsize(1:dim))
    modine%lsize = lsize / modine%Jbar
 
    if (modine%xbar < m_zero .or. modine%xbar > m_one) then
      message(1) = 'The parameter "CurvModineXBar" must lie between 0 and 1.'
      call messages_fatal(1, namespace=namespace)
    end if
 
    safe_allocate(modine%Jlocal(1:npos))
    safe_allocate(modine%Jrange(1:npos))
 
    ! \warning: the reading has to be done for each atom kind
 
    !%Variable CurvModineJlocal
    !%Type float
    !%Default 0.25
    !%Section Mesh::Curvilinear::Modine
    !%Description
    !% Local refinement around the atoms. Must be between 0 and 1.
    !% See N. A. Modine, G. Zumbach, and E. Kaxiras, <i>Phys. Rev. B</i> <b>55</b>, 10289-10301 (1997).
    !%End
    call parse_variable(namespace, 'CurvModineJlocal', 0.25_real64, modine%Jlocal(1))
 
    !%Variable CurvModineJrange
    !%Type float
    !%Default 2 b
    !%Section Mesh::Curvilinear::Modine
    !%Description
    !% Local refinement range (a length).
    !% See N. A. Modine, G. Zumbach, and E. Kaxiras, <i>Phys. Rev. B</i> <b>55</b>, 10289-10301 (1997).
    !%End
    call parse_variable(namespace, 'CurvModineJrange', m_two, modine%Jrange(1), units_inp%length)
 
    if (modine%Jlocal(1)<m_zero.or.modine%Jlocal(1)>m_one) then
      message(1) = 'The parameter "CurvModineJlocal" must lie between 0 and 1.'
      call messages_fatal(1, namespace=namespace)
    end if
 
    modine%Jlocal(:) = modine%Jlocal(1)
    modine%Jrange(:) = modine%Jrange(1)
 
    ! initialize root solver for the optimization
    call root_solver_init(modine%rs, namespace, dim, &
      solver_type = root_newton, maxiter = 500, abs_tolerance = 1.0e-10_real64)
 
    call find_atom_points()
    call optimize()
 
    modine%min_mesh_scaling_product = modine%Jbar**modine%dim
 
    pop_sub(curv_modine_constructor)
  contains
 
    subroutine find_atom_points()
      integer :: iatom, jj
 
      push_sub(curv_modine_constructor.find_atom_points)
 
      ! Initialize csi
      safe_allocate(modine%csi(1:modine%dim, 1:modine%npos))
      modine%csi = npos
 
      ! get first estimate for chi_atoms
      safe_allocate(modine%chi_atoms(1:modine%dim, 1:modine%npos))
      do jj = 1, 10  ! \warning: make something better
        do iatom = 1, modine%npos
          modine%chi_atoms(:, iatom) = modine%from_cartesian(pos(:, iatom))
        end do
        modine%csi(:,:) = modine%chi_atoms(:,:)
      end do
 
      do iatom = 1, modine%npos
        ! These are the chi positions where we want the atoms.
        modine%chi_atoms(:, iatom) = nint(modine%chi_atoms(:, iatom) / spacing(:)) * spacing(:)
      end do
 
      pop_sub(curv_modine_constructor.find_atom_points)
    end subroutine find_atom_points
 
    subroutine optimize()
      logical :: conv
      integer :: iatom, idim, index
      real(real64), allocatable :: my_csi(:), start_csi(:)
 
      push_sub(curv_modine_constructor.optimize)
 
      modine_p  => modine
 
      safe_allocate(x_p(1:modine%dim*modine%npos))
      safe_allocate(my_csi(1:modine%dim*modine%npos))
      safe_allocate(start_csi(1:modine%dim*modine%npos))
 
      do iatom = 1, modine%npos
        do idim = 1, modine%dim
          index = (iatom-1)*dim + idim
          x_p(index)       = pos(idim, iatom)
          start_csi(index) = modine%chi_atoms(idim, iatom)
        end do
      end do
 
      call droot_solver_run(modine%rs, getf2, my_csi, conv, startval=start_csi)
 
      if (.not. conv) then
        message(1) = "During the construction of the adaptive grid, the Newton-Raphson"
        message(2) = "method did not converge."
        call messages_fatal(2, namespace=namespace)
      end if
 
      ! Now set csi to the new values
      do iatom = 1, modine%npos
        do idim = 1, modine%dim
          index = (iatom-1)*modine%dim + idim
          modine_p%csi(idim, iatom) = my_csi(index)
        end do
      end do
 
      safe_deallocate_a(x_p)
      safe_deallocate_a(my_csi)
      safe_deallocate_a(start_csi)
 
      nullify(modine_p)
 
      pop_sub(curv_modine_constructor.optimize)
    end subroutine optimize
 
  end function curv_modine_constructor
 
  ! ---------------------------------------------------------
  subroutine curv_modine_copy(this_out, this_in)
    type(curv_modine_t), intent(inout) :: this_out
    type(curv_modine_t), intent(in)    :: this_in
 
    push_sub(curv_modine_copy)
 
    this_out%dim = this_in%dim
    this_out%local_basis = this_in%local_basis
    safe_allocate_source_a(this_out%lsize, this_in%lsize)
    this_out%xbar = this_in%xbar
    this_out%Jbar = this_in%Jbar
    safe_allocate_source_a(this_out%Jlocal, this_in%Jlocal)
    safe_allocate_source_a(this_out%Jrange, this_in%Jrange)
    safe_allocate_source_a(this_out%chi_atoms, this_in%chi_atoms)
    safe_allocate_source_a(this_out%csi, this_in%csi)
    this_out%npos = this_in%npos
 
    pop_sub(curv_modine_copy)
  end subroutine curv_modine_copy
 
  ! ---------------------------------------------------------
  subroutine curv_modine_finalize(this)
    type(curv_modine_t), intent(inout) :: this
 
    push_sub(curv_modine_finalize)
 
    safe_deallocate_a(this%lsize)
    safe_deallocate_a(this%Jlocal)
    safe_deallocate_a(this%Jrange)
    safe_deallocate_a(this%chi_atoms)
    safe_deallocate_a(this%csi)
 
    pop_sub(curv_modine_finalize)
  end subroutine curv_modine_finalize
 
  ! ---------------------------------------------------------
  function curv_modine_to_cartesian(this, chi) result(xx)
    class(curv_modine_t), target, intent(in)  :: this
    real(real64),                 intent(in)  :: chi(:)
    real(real64) :: xx(1:this%dim)
 
    real(real64) :: chi2(this%dim), rr2, dd
    integer :: iatom
 
    ! no PUSH_SUB, called too often
 
    call curv_modine_chi2chi2(this, chi, chi2)
 
    xx(:) = chi2(:)
    do iatom = 1, this%npos
      rr2 = sum((chi2(:) - this%csi(:, iatom))**2)
      dd = exp(-rr2/(m_two*this%Jrange(iatom)**2))
 
      xx(:) = xx(:) - this%Jlocal(iatom)*(chi2(:) - this%csi(:, iatom)) * dd
    end do
 
  end function curv_modine_to_cartesian
 
  ! ---------------------------------------------------------
  function curv_modine_from_cartesian(this, xx) result(chi)
    class(curv_modine_t), target, intent(in)  :: this
    real(real64),                 intent(in)  :: xx(:)
    real(real64) :: chi(1:this%dim)
 
    logical :: conv
 
    ! no PUSH_SUB, called too often
 
    modine_p  => this
    safe_allocate(x_p(1:this%dim))
    x_p(:) = xx(:)
 
    call droot_solver_run(this%rs, getf, chi, conv, startval = xx)
 
    safe_deallocate_a(x_p)
    nullify(modine_p)
 
    if (.not. conv) then
      message(1) = "During the construction of the adaptive grid, the Newton-Raphson"
      message(2) = "method did not converge for point:"
      write(message(3),'(3f14.6)') xx(1:this%dim)
      call messages_fatal(3)
    end if
 
  end function curv_modine_from_cartesian
 
  ! ---------------------------------------------------------
  subroutine curv_modine_write_info(this, iunit, namespace)
    class(curv_modine_t),         intent(in) :: this
    integer,            optional, intent(in) :: iunit
    type(namespace_t),  optional, intent(in) :: namespace
 
    push_sub(curv_modine_write_info)
 
    write(message(1), '(a)') '  Curvilinear Method = modine'
    call messages_info(1, iunit=iunit, namespace=namespace)
 
    pop_sub(curv_modine_write_info)
  end subroutine curv_modine_write_info
 
  ! ---------------------------------------------------------
  real(real64) function curv_modine_surface_element(this, idir) result(ds)
    class(curv_modine_t), intent(in) :: this
    integer,              intent(in) :: idir
 
    ds = m_zero
    message(1) = 'Surface element with modine curvilinear coordinates not implemented'
    call messages_fatal(1)
 
  end function curv_modine_surface_element
 
  ! ---------------------------------------------------------
  pure subroutine curv_modine_chi2chi2(this, chi_, chi2, Jac)
    class(curv_modine_t), intent(in)  :: this
    real(real64),         intent(in)  :: chi_(:)
    real(real64),         intent(out) :: chi2(:)
    real(real64),      optional, intent(out) :: Jac(:)
 
    integer, parameter :: qq = 3
    real(real64) :: chibar(this%dim), rr, chi
    logical :: neg
    integer :: i
 
    ! no PUSH_SUB, called too often
 
    chibar = this%xbar*this%lsize
 
    do i = 1, this%dim
      neg = (chi_(i) < 0)
      chi = abs(chi_(i))
 
      chi2(i)  = this%Jbar * chi
      if (present(jac)) jac(i) = this%Jbar
 
      if (chi > chibar(i)) then
        rr = (chi - chibar(i))/(this%lsize(i) - chibar(i))
 
        chi2(i)  = chi2(i) + this%lsize(i)/m_two*(1 - this%Jbar) * rr**qq * &
          (qq + m_one - (qq - m_one)*rr)
 
        if (present(jac)) then
          jac(i) = jac(i) + this%lsize(i)/m_two*(m_one - this%Jbar) * rr**(qq - 1)/(this%lsize(i) - chibar(i)) * &
            (qq*(qq + m_one) - (qq**2 - m_one)*rr)
        end if
      end if
 
      if (neg) chi2(i) = -chi2(i)
      ! CHECK if Jacobian does not have to be negated!
    end do
 
  end subroutine curv_modine_chi2chi2
 
  ! ---------------------------------------------------------
  pure subroutine curv_modine_jacobian_inv(this, chi, xx, Jac)
    type(curv_modine_t), intent(in)  :: this
    real(real64),        intent(in)  :: chi(:)
    real(real64),        intent(out) :: xx(:)
    real(real64),        intent(out) :: jac(:, :)
 
    real(real64) :: chi2(this%dim), rr2, dd, j2(this%dim)
    integer :: iatom, idim, idim2
 
    ! no PUSH_SUB, called too often
 
    call curv_modine_chi2chi2(this, chi, chi2, j2)
 
    ! initialize both xx and the Jacobian
    xx(:) = chi2(:)
    jac(:,:) = m_zero
    do idim = 1, this%dim
      jac(idim, idim) = m_one
    end do
 
    do iatom = 1, this%npos
      rr2 = sum((chi2(:) - this%csi(:, iatom))**2)
      dd = exp(-rr2/(m_two*this%Jrange(iatom)**2))
 
      xx(:) = xx(:) -  this%Jlocal(iatom)*(chi2(:) - this%csi(:, iatom)) * dd
 
      do idim = 1, this%dim
        jac(idim, idim) = jac(idim, idim) - this%Jlocal(iatom) * dd
        do idim2 = 1, this%dim
          jac(idim, idim2) = jac(idim, idim2) + &
            this%Jlocal(iatom)*(chi2(idim) - this%csi(idim, iatom))*(chi2(idim2) - this%csi(idim2, iatom)) * &
            m_two/(m_two*this%Jrange(iatom)**2) * dd
        end do
      end do
    end do
 
    do idim = 1, this%dim
      jac(idim, :) = jac(idim, :) * j2(:)
    end do
 
  end subroutine curv_modine_jacobian_inv
 
  ! ---------------------------------------------------------
  pure subroutine getf(yy, ff, jf)
    real(real64), intent(in)  :: yy(:)
    real(real64), intent(out) :: ff(:), jf(:, :)
 
    ! no PUSH_SUB, called too often
 
    call curv_modine_jacobian_inv(modine_p, yy, ff, jf)
    ff(:) = ff(:) - x_p(:)
 
  end subroutine getf
 
  ! ---------------------------------------------------------
  subroutine getf2(csi, ff, jf)
    real(real64), intent(in)  :: csi(:)
    real(real64), intent(out) :: ff(:), jf(:, :)
 
    integer :: i1, j1, i2, j2, index1, index2
    real(real64) :: rr2, dd, dd2
    real(real64), allocatable :: xx(:), chi2(:)
 
    ! no PUSH_SUB, called too often
 
    safe_allocate(xx(1:modine_p%dim))
    safe_allocate(chi2(1:modine_p%dim))
 
    ! first we fill in coord_system%csi with the values we have
    index1 = 1
    do i1 = 1, modine_p%npos
      do j1 = 1, modine_p%dim
        modine_p%csi(j1, i1) = csi(index1)
        index1 = index1 + 1
      end do
    end do
 
    ! get ff and jf
    jf(:,:) = m_zero
    do i1 = 1, modine_p%npos
      call curv_modine_chi2chi2(modine_p, modine_p%chi_atoms(:,i1), chi2)
      xx(:) = chi2(:)
 
      do i2 = 1, modine_p%npos
        rr2 = sum((chi2 - modine_p%csi(:,i2))**2)
        dd = exp(-rr2/(m_two*modine_p%Jrange(i2)**2))
 
        xx = xx - modine_p%Jlocal(i2)*(chi2 - modine_p%csi(:,i2)) * dd
      end do
 
      do j1 = 1, modine_p%dim
        index1 = (i1 - 1)*modine_p%dim + j1
        ff(index1) = xx(j1) - x_p(index1)
 
        do i2 = 1, modine_p%npos
          rr2  = sum((chi2 - modine_p%csi(:,i2))**2)
          dd  = exp(-rr2/(m_two*modine_p%Jrange(i2)**2))
          dd2 = -m_two/(m_two*modine_p%Jrange(i2)**2)*dd
 
          index2 = (i2 - 1)*modine_p%dim + j1
          jf(index1, index2) = modine_p%Jlocal(i2) * dd
 
          do j2 = 1, modine_p%dim
            index2 = (i2 - 1)*modine_p%dim + j2
 
            jf(index1, index2) =  jf(index1, index2) + modine_p%Jlocal(i2) * dd2 * &
              (chi2(j1) - modine_p%csi(j1,i2))*(chi2(j2) - modine_p%csi(j2,i2))
          end do
        end do
      end do
    end do
 
    safe_deallocate_a(xx)
    safe_deallocate_a(chi2)
 
  end subroutine getf2
 
  ! ---------------------------------------------------------
  function curv_modine_jacobian(this, chi) result(jacobian)
    class(curv_modine_t), intent(in)  :: this
    real(real64),         intent(in)  :: chi(:)
    real(real64) :: jacobian(1:this%dim, 1:this%dim)
 
    real(real64) :: dummy(this%dim)
 
    call curv_modine_jacobian_inv(this, chi, dummy, jacobian)
 
  end function curv_modine_jacobian
 
  ! ---------------------------------------------------------
  function curv_modine_jacobian_inverse(this, chi) result(jacobian_inverse)
    class(curv_modine_t), intent(in)  :: this
    real(real64),         intent(in)  :: chi(:)
    real(real64) :: jacobian_inverse(1:this%dim, 1:this%dim)
 
    real(real64) :: dummy(this%dim)
 
    call curv_modine_jacobian_inv(this, chi, dummy, jacobian_inverse)
    call lalg_inverse(this%dim, jacobian_inverse, "dir")
 
  end function curv_modine_jacobian_inverse
 
end module curv_modine_oct_m
 
!! Local Variables:
!! mode: f90
!! coding: utf-8
!! End: