chebyshev_coefficients.F90

Go to the documentation of this file.

!! Copyright (C) 2024 S. Ohlmann
!!
!! 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 chebyshev_coefficients_oct_m
  use debug_oct_m
  use fftw_oct_m
  use fftw_params_oct_m
  use global_oct_m
  use, intrinsic :: iso_c_binding
  use, intrinsic :: iso_fortran_env
  use loct_math_oct_m
  use math_oct_m
  use messages_oct_m
  use profiling_oct_m
 
  implicit none
 
  private
  public ::                            &
    chebyshev_function_t,              &
    chebyshev_exp_t,                   &
    chebyshev_exp_imagtime_t,          &
    chebyshev_numerical_t
 
  type, abstract :: chebyshev_function_t
    private
    real(real64) :: half_span
    real(real64) :: middle_point
    real(real64) :: deltat
  contains
    procedure(chebyshev_get_coefficients), deferred :: get_coefficients
    procedure(chebyshev_get_error), deferred :: get_error
    procedure :: set_parameters => chebyshev_set_parameters
  end type chebyshev_function_t
 
  abstract interface
    subroutine chebyshev_get_coefficients(this, order, coefficients)
      import chebyshev_function_t
      import real64
      implicit none
      class(chebyshev_function_t),  intent(in)  :: this
      integer,                      intent(in)  :: order
      complex(real64), allocatable, intent(out) :: coefficients(:)
    end subroutine chebyshev_get_coefficients
 
    real(real64) function chebyshev_get_error(this, order)
      import chebyshev_function_t
      import real64
      implicit none
      class(chebyshev_function_t),  intent(in)  :: this
      integer,      intent(in)  :: order
    end function chebyshev_get_error
  end interface
 
  
  type, extends(chebyshev_function_t) :: chebyshev_exp_t
  contains
    procedure :: get_coefficients => chebyshev_exp_coefficients
    procedure :: get_error => chebyshev_exp_error
  end type chebyshev_exp_t
 
  interface chebyshev_exp_t
    procedure chebyshev_exp_constructor
  end interface chebyshev_exp_t
 
  
  type, extends(chebyshev_function_t) :: chebyshev_exp_imagtime_t
  contains
    procedure :: get_coefficients => chebyshev_exp_imagtime_coefficients
    procedure :: get_error => chebyshev_exp_imagtime_error
  end type chebyshev_exp_imagtime_t
 
  interface chebyshev_exp_imagtime_t
    procedure chebyshev_exp_imagtime_constructor
  end interface chebyshev_exp_imagtime_t
 
  ! interface needed to save function pointer in type
  abstract interface
    complex(real64) function complex_function_i(z)
      import real64
      implicit none
      complex(real64), intent(in) :: z
    end function complex_function_i
  end interface
 
  
  type, extends(chebyshev_function_t) :: chebyshev_numerical_t
    procedure(complex_function_i), pointer, nopass :: complex_function
  contains
    procedure :: get_coefficients => chebyshev_numerical_coefficients
    procedure :: get_error => chebyshev_numerical_error
  end type chebyshev_numerical_t
 
  interface chebyshev_numerical_t
    procedure chebyshev_numerical_constructor
  end interface chebyshev_numerical_t
 
contains
 
  subroutine chebyshev_set_parameters(this, half_span, middle_point, deltat)
    class(chebyshev_function_t), intent(inout) :: this
    real(real64),                intent(in)    :: half_span
    real(real64),                intent(in)    :: middle_point
    real(real64),                intent(in)    :: deltat
 
    push_sub(chebyshev_set_parameters)
 
    this%half_span = half_span
    this%middle_point = middle_point
    this%deltat = deltat
 
    pop_sub(chebyshev_set_parameters)
  end subroutine chebyshev_set_parameters
 
  function chebyshev_exp_constructor(half_span, middle_point, deltat) result(chebyshev_function)
    class(chebyshev_exp_t), pointer    :: chebyshev_function
    real(real64),           intent(in) :: half_span
    real(real64),           intent(in) :: middle_point
    real(real64),           intent(in) :: deltat
 
    push_sub(chebyshev_exp_constructor)
 
    allocate(chebyshev_function)
    call chebyshev_function%set_parameters(half_span, middle_point, deltat)
 
    pop_sub(chebyshev_exp_constructor)
  end function chebyshev_exp_constructor
 
  subroutine chebyshev_exp_coefficients(this, order, coefficients)
    class(chebyshev_exp_t),       intent(in)  :: this
    integer,                      intent(in)  :: order
    complex(real64), allocatable, intent(out) :: coefficients(:)
 
    integer :: i
 
    push_sub(chebyshev_exp_coefficients)
    safe_allocate(coefficients(0:order))
    do i = 0, order
      coefficients(i) = exp(-m_zi*this%middle_point*this%deltat) * &
        (-m_zi)**i * m_two * loct_bessel(i, this%half_span*this%deltat)
    end do
    coefficients(0) = coefficients(0) / m_two
    pop_sub(chebyshev_exp_coefficients)
  end subroutine chebyshev_exp_coefficients
 
  real(real64) function chebyshev_exp_error(this, order)
    class(chebyshev_exp_t), intent(in)  :: this
    integer,                intent(in)  :: order
 
    real(real64) :: r, inv_r
 
    push_sub(chebyshev_exp_error)
    if (order >= this%half_span*this%deltat) then
      r = m_two*(order+1)/abs(this%half_span*this%deltat)
      inv_r = m_one/r
      chebyshev_exp_error = m_two*inv_r**(order+1)/(m_one - inv_r) * exp(abs(this%half_span*this%deltat)*(r-inv_r)/m_two)
    else
      chebyshev_exp_error = -m_one
    end if
    pop_sub(chebyshev_exp_error)
  end function chebyshev_exp_error
 
  function chebyshev_exp_imagtime_constructor(half_span, middle_point, deltat) result(chebyshev_function)
    class(chebyshev_exp_imagtime_t), pointer    :: chebyshev_function
    real(real64),                    intent(in) :: half_span
    real(real64),                    intent(in) :: middle_point
    real(real64),                    intent(in) :: deltat
 
    push_sub(chebyshev_exp_imagtime_constructor)
 
    allocate(chebyshev_function)
    call chebyshev_function%set_parameters(half_span, middle_point, deltat)
 
    pop_sub(chebyshev_exp_imagtime_constructor)
  end function chebyshev_exp_imagtime_constructor
 
  subroutine chebyshev_exp_imagtime_coefficients(this, order, coefficients)
    class(chebyshev_exp_imagtime_t), intent(in)  :: this
    integer,                         intent(in)  :: order
    complex(real64), allocatable,    intent(out) :: coefficients(:)
 
    integer :: i
 
    push_sub(chebyshev_exp_imagtime_coefficients)
    safe_allocate(coefficients(0:order))
    do i = 0, order
      coefficients(i) = exp(-this%middle_point*this%deltat) * &
        m_two * loct_bessel_in(i, -this%half_span*this%deltat)
    end do
    coefficients(0) = coefficients(0) / m_two
    pop_sub(chebyshev_exp_imagtime_coefficients)
  end subroutine chebyshev_exp_imagtime_coefficients
 
  real(real64) function chebyshev_exp_imagtime_error(this, order)
    class(chebyshev_exp_imagtime_t), intent(in)  :: this
    integer,                         intent(in)  :: order
 
    real(real64) :: upper_bound, lower_bound
    real(real64), parameter :: b = 0.618, d = 0.438
 
    push_sub(chebyshev_exp_imagtime_error)
    upper_bound = this%half_span+this%middle_point
    lower_bound = -this%half_span+this%middle_point
    if (order <= upper_bound*this%deltat) then
      chebyshev_exp_imagtime_error = m_two*exp(-lower_bound*this%deltat)*exp(-b*(order+1)**2/(upper_bound*this%deltat))*&
        (m_one+sqrt(upper_bound*this%deltat*m_pi/(m_four*b))) + m_two*d**(upper_bound*this%deltat)/(m_one-d)
    else
      chebyshev_exp_imagtime_error = m_two*exp(-lower_bound*this%deltat)*d**order/(m_one-d)
    end if
    pop_sub(chebyshev_exp_imagtime_error)
  end function chebyshev_exp_imagtime_error
 
  function chebyshev_numerical_constructor(half_span, middle_point, deltat, complex_function) result(chebyshev_function)
    class(chebyshev_numerical_t),  pointer    :: chebyshev_function
    real(real64),                  intent(in) :: half_span
    real(real64),                  intent(in) :: middle_point
    real(real64),                  intent(in) :: deltat
    procedure(complex_function_i)             :: complex_function
 
    push_sub(chebyshev_numerical_constructor)
 
    allocate(chebyshev_function)
    call chebyshev_function%set_parameters(half_span, middle_point, deltat)
    chebyshev_function%complex_function => complex_function
 
    pop_sub(chebyshev_numerical_constructor)
  end function chebyshev_numerical_constructor
 
  subroutine chebyshev_numerical_coefficients(this, order, coefficients)
    class(chebyshev_numerical_t), intent(in)  :: this
    integer,                      intent(in)  :: order
    complex(real64), allocatable, intent(out) :: coefficients(:)
 
    integer :: i
    complex(real64) :: values(0:order)
    real(real64) :: re_values(0:order), im_values(0:order)
    real(real64) :: re_coefficients(0:order), im_coefficients(0:order)
    real(real64) :: theta
    type(c_ptr) :: plan
 
    push_sub(chebyshev_numerical_coefficients)
    safe_allocate(coefficients(0:order))
    ! compute values for DCT
    do i = 0, order
      theta = m_pi*i/order
      values(i) = this%complex_function(-m_zi*this%deltat*(this%half_span*cos(theta)+this%middle_point))
      re_values(i) = real(values(i), real64)
      im_values(i) = aimag(values(i))
    end do
    plan = fftw_plan_r2r_1d(order+1, re_values, re_coefficients, fftw_redft00, fftw_estimate)
    call fftw_execute_r2r(plan, re_values, re_coefficients)
    call fftw_execute_r2r(plan, im_values, im_coefficients)
    call fftw_destroy_plan(plan)
    do i = 0, order
      ! normalize correctly, also need factor 2 for expansion coefficients
      coefficients(i) = cmplx(re_coefficients(i), im_coefficients(i), real64)/order
    end do
    ! need to divide first and last coefficient by 2
    coefficients(0) = coefficients(0) / m_two
    coefficients(order) = coefficients(order) / m_two
 
    pop_sub(chebyshev_numerical_coefficients)
  end subroutine chebyshev_numerical_coefficients
 
  real(real64) function chebyshev_numerical_error(this, order)
    class(chebyshev_numerical_t), intent(in)  :: this
    integer,                      intent(in)  :: order
 
    real(real64) :: r
    real(real64) :: inv_r
 
    push_sub(chebyshev_numerical_error)
    if (order >= this%half_span*this%deltat) then
      r = m_two*order/abs(this%half_span*this%deltat)
      inv_r = m_one/r
      chebyshev_numerical_error = m_two*inv_r**order/(m_one - inv_r) * &
        real(this%complex_function(cmplx(abs(this%half_span*this%deltat)*(r-inv_r)/M_TWO, M_ZERO, real64)), real64)
    else
      chebyshev_numerical_error = -m_one
    end if
    pop_sub(chebyshev_numerical_error)
  end function chebyshev_numerical_error
end module chebyshev_coefficients_oct_m
 
!! Local Variables:
!! mode: f90
!! coding: utf-8
!! End: