X_vars.f90

Go to the documentation of this file.

!------------------------------------------------------------------------------!
!> Variables pertaining to the perturbation quantities.
!------------------------------------------------------------------------------!
module x_vars
#include <PB3D_macros.h>
    use str_utilities
    use messages
    use num_vars, only: dp, max_name_ln, iu, weight_dp
    use grid_vars, only: grid_type
    
    implicit none
    
    private
    public set_nm_x, set_nn_mod, &
        &n_mod_x, prim_x, min_sec_x, max_sec_x, min_nm_x, min_n_x, max_n_x, &
        &min_m_x, max_m_x, min_r_sol, max_r_sol, mds_x, mds_sol
#if ldebug
    public n_alloc_x_1s, n_alloc_x_2s
#endif
    
    !> mode number type
    !!
    !! Type containing information about mode numbers at every flux surface:
    !!  - mode numbers: \c n, \c m,
    !!  - information of secondary modes: \c sec,
    !!      * 1: mode number
    !!      * 2: lower limit of normal range
    !!      * 3: upper limit of normal range
    !!      * 4: index in tables
    !!      where the number of modes can be bigger than the range if the safety
    !!      factor or  rotational transform is  non-monotonous so that  the same
    !!      mode number can have multiple ranges.
    !!  - Flux normal coordinate: \c r_F.
    !!
    !! These are set up in setup_nm_x().
    type, public :: modes_type
        integer, allocatable :: n(:,:)                                          !< \f$n\f$ for all modes, in total grid
        integer, allocatable :: m(:,:)                                          !< \f$m\f$ for all modes, in total grid
        integer, allocatable :: sec(:,:)                                        !< \c m or \c n for all possible modes, index and limits, in total grid
    contains
        !> deallocate
        procedure :: dealloc => dealloc_mds
    end type
    
    !> vectorial perturbation type
    !!
    !! The arrays here are of the form:
    !!  - \c U_x_i and DU_X_i: <tt>(1:angle_1,1:angle_2,1;n_mod)</tt>
    !!
    !! \see See grid_vars.grid_type for a discussion on \c ang_1 and \c ang_2.
    type, public :: x_1_type
        integer :: n_mod                                                        !< size of \f$n\f$ and \f$m\f$ (nr. of modes)
        integer :: lim_sec_x(2)                                                 !< limits of \c m_X (pol. flux) or \c n_X (tor. flux)
        integer, allocatable :: n(:,:)                                          !< vector of poloidal mode numbers
        integer, allocatable :: m(:,:)                                          !< vector of poloidal mode numbers
        complex(dp), allocatable :: u_0(:,:,:,:)                                !< \f$U_m^0\f$
        complex(dp), allocatable :: u_1(:,:,:,:)                                !< \f$U_m^1\f$
        complex(dp), allocatable :: du_0(:,:,:,:)                               !< \f$\mathcal{J}\vec{B}\cdot\nabla U_m^0\f$
        complex(dp), allocatable :: du_1(:,:,:,:)                               !< \f$\mathcal{J}\vec{B}\cdot\nabla U_m^1\f$
#if ldebug
        real(dp) :: estim_mem_usage                                             !< estimated memory usage \ldebug
#endif
    contains
        !> initialize
        procedure :: init => init_x_1
        !> copy
        procedure :: copy => copy_x_1
        !> deallocate
        procedure :: dealloc => dealloc_x_1
    end type
    
    !> tensorial perturbation type
    !!
    !! The arrays here are of the form:
    !!  - \c PV_i and KV_i: <tt>(1:angle_1,1:angle_2,1;n_mod^2)</tt>
    !!
    !! \see See grid_vars.grid_type for a discussion on \c ang_1 and \c ang_2.
    !!
    !! \note This  type is also  used for field-averaged  tensorial perturbation
    !! variables, with \c angle_1 of size 1.
    type, public :: x_2_type
        integer :: n_mod(2)                                                     !< size of \f$n\f$ and \f$m\f$ (nr. of modes)
        integer :: lim_sec_x(2,2)                                               !< limits of \c m_X (pol. flux) or \c n_X (tor. flux)
        integer, allocatable :: n_1(:,:)                                        !< vector of toroidal mode numbers of dimension 1
        integer, allocatable :: n_2(:,:)                                        !< vector of toroidal mode numbers of dimension 2
        integer, allocatable :: m_1(:,:)                                        !< vector of poloidal mode numbers of dimension 1
        integer, allocatable :: m_2(:,:)                                        !< vector of poloidal mode numbers of dimension 2
        complex(dp), allocatable :: pv_0(:,:,:,:)                               !< \f$\widetilde{PV}^0\f$ coefficient
        complex(dp), allocatable :: pv_1(:,:,:,:)                               !< \f$\widetilde{PV}^1\f$ coefficient
        complex(dp), allocatable :: pv_2(:,:,:,:)                               !< \f$\widetilde{PV}^2\f$ coefficient
        complex(dp), allocatable :: kv_0(:,:,:,:)                               !< \f$\widetilde{KV}^0\f$ coefficient
        complex(dp), allocatable :: kv_1(:,:,:,:)                               !< \f$\widetilde{KV}^1\f$ coefficient
        complex(dp), allocatable :: kv_2(:,:,:,:)                               !< \f$\widetilde{KV}^2\f$ coefficient
#if ldebug
        real(dp) :: estim_mem_usage                                             !< estimated memory usage \ldebug
#endif
    contains
        !> initialize
        procedure :: init => init_x_2
        !> copy
        procedure :: copy => copy_x_2
        !> deallocate
        procedure :: dealloc => dealloc_x_2
    end type
    
    ! interfaces
    
    !> \public Sets \c n_X and \c m_X.
    !!
    !! By default, this is done using  by default global \c X_vars variables but
    !! optionally different  limits for the  secondary mode numbers (\c  m_X for
    !! poloidal flux or \c n_X for toroidal flux).
    !!
    !! \note \c n_X and \c m_X need to  have been set up with the same limits as
    !! the grid used here. This is done in setup_nm_X().
    interface set_nm_x
        !> \public
        module procedure set_nm_x_1
        !> \public
        module procedure set_nm_x_2
    end interface
    
    ! global variables
    type(modes_type) :: mds_x                                                   !< modes variables for perturbation grid
    type(modes_type) :: mds_sol                                                 !< modes variables for solution grid
    integer :: prim_x                                                           !< \c n_X (pol. flux) or \c m_X (tor. flux)
    integer :: min_sec_x                                                        !< \c m_X (pol. flux) or \c n_X (tor. flux) (only for \c X style 1)
    integer :: max_sec_x                                                        !< \c m_X (pol. flux) or \c n_X (tor. flux) (only for\ c X style 1)
    integer :: n_mod_x                                                          !< size of \c m_X (pol. flux) or \c n_X (tor. flux)
    integer :: min_nm_x = 5                                                     !< minimum for the high-n theory (debable)
    integer, allocatable :: min_n_x(:)                                          !< lowest poloidal mode number \c m_X, in total eq grid
    integer, allocatable :: max_n_x(:)                                          !< highest poloidal mode number \c m_X, in total eq grid
    integer, allocatable :: min_m_x(:)                                          !< lowest poloidal mode number \c m_X, in total eq grid
    integer, allocatable :: max_m_x(:)                                          !< highest poloidal mode number \c m_X, in total eq grid
    real(dp) :: min_r_sol                                                       !< min. normal range for pert.
    real(dp) :: max_r_sol                                                       !< max. normal range for pert.
#if ldebug
    integer :: n_alloc_x_1s                                                     !< nr. of allocated <tt>X_1</tt>'s \ldebug
    integer :: n_alloc_x_2s                                                     !< nr. of allocated <tt>X_2</tt>'s \ldebug
#endif
    
contains
    !> \private vectorial version
    subroutine set_nm_x_1(mds,grid_X,lim_sec_X_o,n_X_loc,m_X_loc,lim_sec_X_i)
        
        ! input / output
        type(modes_type), intent(in) :: mds                                     !< general modes variables
        type(grid_type), intent(in) :: grid_X                                   !< perturbation grid
        integer, intent(inout) :: lim_sec_X_o(2)                                !< limits on secondary mode numbers
        integer, intent(inout), allocatable :: n_X_loc(:,:)                     !< toroidal mode numbers
        integer, intent(inout), allocatable :: m_X_loc(:,:)                     !< poloidal mode numbers
        integer, intent(in), optional :: lim_sec_X_i(2)                         !< optional limits on secondary mode numbers
        
        ! set output lim_sec_X
        lim_sec_x_o = [1,n_mod_x]
        if (present(lim_sec_x_i)) lim_sec_x_o = lim_sec_x_i
        
        ! set n and m
        allocate(n_x_loc(grid_x%loc_n_r,lim_sec_x_o(2)-lim_sec_x_o(1)+1))
        allocate(m_x_loc(grid_x%loc_n_r,lim_sec_x_o(2)-lim_sec_x_o(1)+1))
        n_x_loc = mds%n(grid_x%i_min:grid_x%i_max,&
            &lim_sec_x_o(1):lim_sec_x_o(2))
        m_x_loc = mds%m(grid_x%i_min:grid_x%i_max,&
            &lim_sec_x_o(1):lim_sec_x_o(2))
    end subroutine set_nm_x_1
    !> \private tensorial version
    subroutine set_nm_x_2(mds,grid_X,lim_sec_X_o,n_X_1,m_X_1,n_X_2,m_X_2,&
        &lim_sec_X_i)
        
        ! input / output
        type(modes_type), intent(in) :: mds                                     !< general modes variables
        type(grid_type), intent(in) :: grid_X                                   !< perturbation grid
        integer, intent(inout) :: lim_sec_X_o(2,2)                              !< limits on secondary mode numbers
        integer, intent(inout), allocatable :: n_X_1(:,:)                       !< toroidal mode numbers for dimension 1
        integer, intent(inout), allocatable :: m_X_1(:,:)                       !< poloidal mode numbers for dimension 1
        integer, intent(inout), allocatable :: n_X_2(:,:)                       !< toroidal mode numbers for dimension 2
        integer, intent(inout), allocatable :: m_X_2(:,:)                       !< poloidal mode numbers for dimension 2
        integer, intent(in), optional :: lim_sec_X_i(2,2)                       !< optional limits on secondary mode numbers
        
        ! call vectorial version
        if (present(lim_sec_x_i)) then
            call set_nm_x(mds,grid_x,lim_sec_x_o(:,1),n_x_1,m_x_1,&
                &lim_sec_x_i(:,1))
            call set_nm_x(mds,grid_x,lim_sec_x_o(:,2),n_x_2,m_x_2,&
                &lim_sec_x_i(:,2))
        else
            call set_nm_x(mds,grid_x,lim_sec_x_o(:,1),n_x_1,m_x_1)
            call set_nm_x(mds,grid_x,lim_sec_x_o(:,2),n_x_2,m_x_2)
        end if
    end subroutine set_nm_x_2
    
    !> \public Initializes a vectorial perturbation.
    !!
    !! Allocates the variables, the number of modes, as well as \c n and \c m .
    !!
    !! Optionally, the secondary mode numbers can be specified (\c m if poloidal
    !! flux is used and \c n if  toroidal flux). By default, they are taken from
    !! the global \c X_vars variables.
    !!
    !! \note If the lowest limits of  the grid is not 1 (e.g. <tt>grid_sol%i_min
    !! = 1</tt> for first process), the input variable \c i_min should be set to
    !! set correctly. For a full grid, it should be set to 1.
    subroutine init_x_1(X,mds,grid_X,lim_sec_X)
#if ldebug
        use num_vars, only: print_mem_usage, rank
#endif
        
        ! input / output
        class(x_1_type), intent(inout) :: X                                     !< vectorial perturbation variables
        type(modes_type), intent(in) :: mds                                     !< general modes variables
        type(grid_type), intent(in) :: grid_X                                   !< perturbation grid
        integer, intent(in), optional :: lim_sec_X(2)                           !< limits of \c m_X (pol. flux) or \c n_X (tor. flux)
        
        ! local variables
        integer :: loc_n_r                                                      ! local nr. of normal points
        integer :: n_par, n_geo                                                 ! tot. nr. of angular points in parallel and geodesic direction
        
        ! set local variables
        loc_n_r = grid_x%loc_n_r
        n_par = grid_x%n(1)
        n_geo = grid_x%n(2)
        
#if ldebug
        ! initialize memory usage
        if (print_mem_usage) x%estim_mem_usage = 0._dp
#endif
        
        ! set mode numbers
        call set_nm_x(mds,grid_x,x%lim_sec_X,x%n,x%m,lim_sec_x)
        
        ! set n_mod
        x%n_mod = size(x%n,2)
        
        ! allocate U_0
        allocate(x%U_0(n_par,n_geo,loc_n_r,x%n_mod))
        
        ! allocate U_1
        allocate(x%U_1(n_par,n_geo,loc_n_r,x%n_mod))
        
        ! allocate DU_0
        allocate(x%DU_0(n_par,n_geo,loc_n_r,x%n_mod))
        
        ! allocate DU_1
        allocate(x%DU_1(n_par,n_geo,loc_n_r,x%n_mod))
        
#if ldebug
        ! set estimated memory usage
        if (print_mem_usage) x%estim_mem_usage = x%estim_mem_usage + &
            &(n_par*n_geo*loc_n_r)*(x%n_mod*4)
        
        ! increment n_alloc_X_1s
        n_alloc_x_1s = n_alloc_x_1s + 1
        
        ! print memory usage
        if (print_mem_usage) call writo('[rank '//trim(i2str(rank))//&
            &' - Expected memory usage of X_1: '//&
            &trim(r2strt(x%estim_mem_usage*weight_dp*2))//' kB]',alert=.true.)
#endif
    end subroutine init_x_1
    !> \public Initializes a tensorial perturbation.
    !!
    !! \see See init_X_1().
    !!
    !! \note The tensorial perturbation type  can also be used for field-aligned
    !! variables, in which  case the first index is assumed  to have dimension 1
    !! only.  This can  be triggered  using  \c is_field_averaged.  There is  no
    !! difference  between a  tensorial  perturbation type  with  size of  first
    !! dimension  set to  one through  the use  of this  flag, or  through other
    !! means.
    subroutine init_x_2(X,mds,grid_X,lim_sec_X,is_field_averaged)
#if ldebug
        use num_vars, only: print_mem_usage, rank
#endif
        
        ! input / output
        class(x_2_type), intent(inout) :: X                                     !< tensorial perturbation variables
        type(modes_type), intent(in) :: mds                                     !< general modes variables
        type(grid_type), intent(in) :: grid_X                                   !< perturbation grid
        integer, intent(in), optional :: lim_sec_X(2,2)                         !< limits of \c m_X (pol. flux) or \c n_X (tor. flux) for both dimensions
        logical, intent(in), optional :: is_field_averaged                      !< if field-aligned, only one dimension for first index
        
        ! local variables
        integer :: loc_n_r                                                      ! local nr. of normal points
        integer :: n_par, n_geo                                                 ! tot. nr. of angular points in parallel and geodesic direction
        integer :: nn_mod                                                       ! either n_mod^2, n_mod*(n_mod+1)/2 or 0
        
        ! set local variables
        n_par = grid_x%n(1)
        if (present(is_field_averaged)) then
            if (is_field_averaged) n_par = 1
        end if
        n_geo = grid_x%n(2)
        loc_n_r = grid_x%loc_n_r
        
#if ldebug
        ! initialize memory usage
        if (print_mem_usage) x%estim_mem_usage = 0._dp
#endif
        
        ! set mode numbers
        call set_nm_x(mds,grid_x,x%lim_sec_X,x%n_1,x%m_1,x%n_2,x%m_2,lim_sec_x)
        
        ! set n_mod
        x%n_mod(1) = size(x%n_1,2)
        x%n_mod(2) = size(x%n_2,2)
        
        ! set nnmod for symmetric quantities
        nn_mod = set_nn_mod(.true.,lim_sec_x)
        
        ! allocate PV_i
        allocate(x%PV_0(n_par,n_geo,loc_n_r,nn_mod))                            ! symmetric
        allocate(x%PV_1(n_par,n_geo,loc_n_r,product(x%n_mod)))                  ! not symmetric
        allocate(x%PV_2(n_par,n_geo,loc_n_r,nn_mod))                            ! symmetric
        
        ! allocate KV_i
        allocate(x%KV_0(n_par,n_geo,loc_n_r,nn_mod))                            ! symmetric
        allocate(x%KV_1(n_par,n_geo,loc_n_r,product(x%n_mod)))                  ! not symmetric
        allocate(x%KV_2(n_par,n_geo,loc_n_r,nn_mod))                            ! symmetric
        
#if ldebug
        ! set estimated memory usage
        if (print_mem_usage) x%estim_mem_usage = &
            &x%estim_mem_usage + (n_par*n_geo*loc_n_r)*&
            &(nn_mod*4+product(x%n_mod)*2) + product(x%n_mod)
        
        ! increment n_alloc_X_2s
        n_alloc_x_2s = n_alloc_x_2s + 1
        
        ! print memory usage
        if (print_mem_usage) call writo('[rank '//trim(i2str(rank))//&
            &' - Expected memory usage of X_2: '//&
            &trim(r2strt(x%estim_mem_usage*weight_dp*2))//' kB]',alert=.true.)
#endif
    end subroutine init_x_2
    
    !> \public Deep copy of vectorial perturbation variables.
    subroutine copy_x_1(X_i,mds,grid_i,X_o)
        use grid_vars, only: grid_type
        
        ! input / output
        class(x_1_type), intent(in) :: X_i                                      !< X_1 to be copied
        type(modes_type), intent(in) :: mds                                     !< general modes variables
        type(grid_type), intent(in) :: grid_i                                   !< grid of eq_i
        type(x_1_type), intent(inout) :: X_o                                    !< copied X_1
        
        call x_o%init(mds,grid_i,lim_sec_x=x_i%lim_sec_X)
        
        ! U_i
        x_o%U_0 = x_i%U_0
        x_o%U_1 = x_i%U_1
        
        ! DU_i
        x_o%DU_0 = x_i%DU_0
        x_o%DU_1 = x_i%DU_1
    end subroutine copy_x_1
    
    !> \public Deep copy of tensorial perturbation variables.
    subroutine copy_x_2(X_i,mds,grid_i,X_o)
        use grid_vars, only: grid_type
        
        ! input / output
        class(x_2_type), intent(in) :: X_i                                      !< X_2 to be copied
        type(modes_type), intent(in) :: mds                                     !< general modes variables
        type(grid_type), intent(in) :: grid_i                                   !< grid of eq_i
        type(x_2_type), intent(inout) :: X_o                                    !< copied X_1
        
        call x_o%init(mds,grid_i,lim_sec_x=x_i%lim_sec_X,&
            &is_field_averaged=size(x_i%PV_0,1).eq.1)
        
        ! PV_i
        x_o%PV_0 = x_i%PV_0
        x_o%PV_1 = x_i%PV_1
        x_o%PV_2 = x_i%PV_2
        
        ! KV_i
        x_o%KV_0 = x_i%KV_0
        x_o%KV_1 = x_i%KV_1
        x_o%KV_2 = x_i%KV_2
    end subroutine copy_x_2
    
    !> Sets number of entries for tensorial perturbation variables.
    integer function set_nn_mod(sym,lim_sec_X) result(nn_mod)
        ! input / output
        logical, intent(in) :: sym                                              !< whether the variable is symmetric
        integer, intent(in), optional :: lim_sec_x(2,2)                         !< limits of \c m_X (pol flux) or \c n_X (tor flux) for both dimensions
        
        ! local variables
        integer :: id, jd                                                       ! counters
        integer :: lim_sec_x_loc(2,2)                                           ! local version of lim_sec_X
        
        ! set local lim_sec_X
        lim_sec_x_loc(:,1) = [1,n_mod_x]
        lim_sec_x_loc(:,2) = [1,n_mod_x]
        if (present(lim_sec_x)) lim_sec_x_loc = lim_sec_x
        
        if (sym) then
            ! set nnmod for symmetric quantities: discard indices above diagonal
            nn_mod = 0
            do jd = lim_sec_x_loc(1,2),lim_sec_x_loc(2,2)
                do id = lim_sec_x_loc(1,1),lim_sec_x_loc(2,1)
                    if (id.ge.jd) nn_mod = nn_mod + 1
                end do
            end do
        else
            ! set nn_mod for asymmetric quantities: don't discard anything
            nn_mod = product(lim_sec_x_loc(2,:)-lim_sec_x_loc(1,:)+1)
        end if
    end function set_nn_mod
    
    !> \public Deallocates modes variables
    subroutine dealloc_mds(mds)
        ! input / output
        class(modes_type), intent(inout) :: mds                                 !< modes variables to be deallocated
        
        ! deallocate allocatable variables
        call dealloc_mds_final(mds)
    contains
        ! Note: intent(out) automatically deallocates the variable
        !> \private
        subroutine dealloc_mds_final(mds)
            ! input / output
            type(modes_type), intent(out) :: mds                                ! modes variables to be deallocated
        end subroutine dealloc_mds_final
    end subroutine dealloc_mds
    
    !> \public Deallocates vectorial perturbation variables.
    subroutine dealloc_x_1(X)
#if ldebug
        use num_vars, only: rank, print_mem_usage
#endif
        
        ! input / output
        class(x_1_type), intent(inout) :: X                                     !< perturbation variables to be deallocated
        
#if ldebug
        ! local variables
        integer :: mem_diff                                                     ! difference in memory
        real(dp) :: estim_mem_usage                                             ! estimated memory usage
        
        ! memory usage before deallocation
        if (print_mem_usage) then
            mem_diff = get_mem_usage()
            estim_mem_usage = x%estim_mem_usage
        end if
#endif
        
        ! deallocate allocatable variables
        call dealloc_x_1_final(x)
        
#if ldebug
        ! decrement n_alloc_X_1s
        n_alloc_x_1s = n_alloc_x_1s - 1
        
        ! memory usage difference after deallocation
        if (print_mem_usage) then
            mem_diff = mem_diff - get_mem_usage()
            call writo('[Rank '//trim(i2str(rank))//' - liberated '//&
                &trim(i2str(mem_diff))//'kB deallocating X_1 ('//&
                &trim(i2str(nint(100*mem_diff/&
                &(estim_mem_usage*weight_dp*2))))//&
                &'% of estimated)]',alert=.true.)
        end if
#endif
    contains
        ! Note: intent(out) automatically deallocates the variable
        !> \private
        subroutine dealloc_x_1_final(X)
            ! input / output
            type(x_1_type), intent(out) :: X                                    ! equilibrium to be deallocated
        end subroutine dealloc_x_1_final
    end subroutine dealloc_x_1
    !> \public Deallocates tensorial perturbation variables.
    subroutine dealloc_x_2(X)
#if ldebug
        use num_vars, only: rank, print_mem_usage
#endif
        
        ! input / output
        class(x_2_type), intent(inout) :: X                                     !< perturbation variables to be deallocated
        
#if ldebug
        ! local variables
        integer :: mem_diff                                                     ! difference in memory
        real(dp) :: estim_mem_usage                                             ! estimated memory usage
        
        ! memory usage before deallocation
        if (print_mem_usage) then
            mem_diff = get_mem_usage()
            estim_mem_usage = x%estim_mem_usage
        end if
#endif
        
        ! deallocate allocatable variables
        call dealloc_x_2_final(x)
        
#if ldebug
        ! decrement n_alloc_X_2s
        n_alloc_x_2s = n_alloc_x_2s - 1
        
        ! memory usage difference after deallocation
        if (print_mem_usage) then
            mem_diff = mem_diff - get_mem_usage()
            call writo('[Rank '//trim(i2str(rank))//' - liberated '//&
                &trim(i2str(mem_diff))//'kB deallocating X_2 ('//&
                &trim(i2str(nint(100*mem_diff/&
                &(estim_mem_usage*weight_dp*2))))//&
                &'% of estimated)]',alert=.true.)
        end if
#endif
    contains
        ! Note: intent(out) automatically deallocates the variable
        !> \private
        subroutine dealloc_x_2_final(X)
            ! input / output
            type(x_2_type), intent(out) :: X                                    ! equilibrium to be deallocated
        end subroutine dealloc_x_2_final
    end subroutine dealloc_x_2
end module x_vars