parallel_map_m.f90 Source File


Contents

Source Code


Source Code

module parallel_map_m
    !! Parallel operators and mapping
    use MPI
    use netcdf
    use perf_m, only : perf_start, perf_stop
    use comm_handler_m, only : comm_handler_t
    use parallelisation_setup_m, only : is_rank_info_writer
    use precision_grillix_m, only : GP
    use type_conversion_m, only : logical_to_integer, integer_to_logical
    use error_handling_grillix_m, only: handle_error_netcdf, handle_error, error_info_t
    use status_codes_grillix_m, only : GRILLIX_ERR_OTHER, &
                                       GRILLIX_ERR_NAMELIST, GRILLIX_ERR_ALLOC
    
    ! From PARALLAX
    use constants_m, only : two_pi
    use screen_io_m, only :  get_stdout
    use equilibrium_m, only : equilibrium_t
    use csrmat_m, only : csrmat_t, csr_transpose, csr_times_vec_single
    use mesh_cart_m, only : mesh_cart_t
    use map_factory_m, only : create_map_matrix
    use fieldline_tracer_m, only : trace
    use mpi_mapping_auxiliaries_m, only : getdata_fwdbwdplane
#if ENABLE_CUDA
    use iso_c_binding, only : c_ptr, c_null_ptr
    use allocate_data_cuda, only : allocate_memory_int, allocate_memory_double
#endif
    implicit none

#if ENABLE_CUDA
    type, private :: cuda_map_t
        !! Defines structure that contains pointers, exposing data to CUDA.
        !! Needed to evaluate parallel operators via CUDA
        integer :: nrows
        !! Number of rows of matrix
        integer :: nnz
        !! Number of columns of matrix
        integer, pointer, dimension(:) :: rows
        !! Pointer to row array of CSR matrix
        integer, pointer, dimension(:) :: cols
        !! Pointer to column indices  of CSR matrix
        double precision, pointer, dimension(:) :: vals
        !! Pointer to values of CSR matrix
        double precision, pointer, dimension(:) :: x
        !! Pointer to input array of parallel operator routine 
        double precision, pointer, dimension(:) :: y
        !! Pointer to output array of parallel operator routine 
        type(c_ptr) :: cuda_struct = c_null_ptr
        !! Pointer for exposure to CUDA ???
    contains
        final :: destructor_cuda_map_t
    end type

    interface
        module subroutine destructor_cuda_map_t(self)
            !! Frees memory associated with cuda_map_t
            type(cuda_map_t), intent(inout) :: self
            !! Instance of the type
        end subroutine
    end interface

#endif
  
    type, public :: parallel_map_t
        !! Datatype for parallel_map
    
        ! >>>>>>>>>>>
        ! Metadata
        integer, private :: nplanes
        !! Number of poloidal planes of map
        real(GP), private :: iplane
        !! Index of poloidal plane
        real(GP), private :: phi_cano
        !! phi coordinate of canonical mesh
        real(GP), private :: phi_stag
        !! phi coordinate of staggered mesh
        real(GP), private :: dphi       
        !! Toroidal grid distance of map
        integer, private :: intorder
        !! Interpolation order of map
        integer, private :: xorder
        !! Order of subcell splitting for q-matrices
        logical, private :: use_fixed_stencil
        !! Indicates interpolation method for the mapped quadrature points:
        !!  ... if false (default), the interpolation stencil is chosen
        !!      individually for each mapped quadrature point
        !!  ... if true, the interpolation stencil is fixed as the mapped
        !!      midpoint's stencil
        logical, private :: use_gauss_quadrature
        !! Indicates quadrature formula for mapped integral
        !! ... if false (default), uniform mapped midpoint quadrature
        !! ... if true, mapped tensor Gauss-Legendre quadrature
        
        ! >>>>>>>>>>>
        ! Arrays
        real(GP), dimension(:), allocatable, public :: dpar_cano_stag_fwd
        !! Parallel distance from phi_cano to phi_cano+dphi/2
        
        real(GP), dimension(:), allocatable, public :: dpar_cano_stag_bwd
        !! Parallel distance from phi_cano to phi_cano-dphi/2
        
        real(GP), dimension(:), allocatable, public :: dpar_cano_cano_fwd
        !! Parallel distance from phi_cano to phi_cano+dphi
        
        real(GP), dimension(:), allocatable, public :: dpar_cano_cano_bwd
        !! Parallel distance from phi_cano to phi_cano-dphi
        
        real(GP), dimension(:), allocatable, public  :: fluxbox_vol_cano
        !! Fluxbox volume from phi_cano-dphi/2 to phi_cano+dphi/2
        
        real(GP), pointer, dimension(:), public  :: dpar_stag_cano_fwd => null()
        !! Parallel distance from phi_stag to phi_stag+dphi/2
        !! Points to dpar_stag_cano_fwd_mem for non-axisymmetric case
        !! Points to dpar_cano_stag_fwd for axisymmetric case
        real(GP), allocatable, dimension(:), private  :: dpar_stag_cano_fwd_mem
        !! Parallel distance from phi_stag to phi_stag+dphi/2
        !! Only built for non-axisymmetric case, 
        !! to be targeted by dpar_stag_cano_fwd
        
        real(GP), pointer, dimension(:), public  :: dpar_stag_cano_bwd => null()
        !! Parallel distance from phi_stag to phi_stag-dphi/2
        !! Points to dpar_stag_cano_bwd_mem for non-axisymmetric case
        !! Points to dpar_cano_stag_bwd for axisymmetric case
        real(GP), allocatable, dimension(:), private  :: dpar_stag_cano_bwd_mem
        !! Parallel distance from phi_stag to phi_stag-dphi/2
        !! Only built for non-axisymmetric case, 
        !! to be targeted by dpar_stag_cano_bwd
        
        real(GP), pointer, dimension(:), public  :: dpar_stag_stag_fwd => null()
        !! Parallel distance from phi_stag to phi_stag+dphi
        !! Points to dpar_stag_stag_fwd_mem for non-axisymmetric case
        !! Points to dpar_cano_cano_fwd for axisymmetric case
        real(GP), allocatable, dimension(:), private  :: dpar_stag_stag_fwd_mem
        !! Parallel distance from phi_stag to phi_stag+dphi
        !! Only built for non-axisymmetric case, 
        !! to be targeted by dpar_stag_stag_fwd
        
        real(GP), pointer, dimension(:), public  :: dpar_stag_stag_bwd => null()
        !! Parallel distance from phi_stag to phi_stag-dphi
        !! Points to dpar_stag_stag_bwd_mem for non-axisymmetric case
        !! Points to dpar_cano_cano_bwd for axisymmetric case
        real(GP), allocatable, dimension(:), private  :: dpar_stag_stag_bwd_mem
        !! Parallel distance from phi_stag to phi_stag-dphi
        !! Only built for non-axisymmetric case, 
        !! to be targeted by dpar_stag_stag_bwd
        
        real(GP), pointer, dimension(:), public :: fluxbox_vol_stag => null()
        !! Fluxbox volume for staggered grid from phi_stag-dphi/2 to phi_stag+dphi/2
        !! Points to fluxbox_vol_stag_mem for non-axisymmetric case
        !! Points to fluxbox_vol_cano for axisymmetric case
        real(GP), allocatable, dimension(:), private  :: fluxbox_vol_stag_mem
        !! Fluxbox volume for staggered grid from phi_stag-dphi/2 to phi_stag+dphi/2
        !! Only built for non-axisymmetric case, 
        !! to be targeted by fluxbox_vol_stag
        
        ! >>>>>>>>>>>
        ! Map matrices
        ! Note: Matrices are public, in order to be tested in unit tests,
        ! In applications, rather do not access them directly but use the corresponding 
        ! routines (par_grad, par_divb, lift_, flat_, upstream_)
        type(csrmat_t), public, allocatable :: map_stag_cano_fwd
        !! Map matrix to staggered mesh (phi_stag) 
        !! from forward canonical mesh (phi_stag+dphi/2)
        
        type(csrmat_t), public, allocatable :: map_stag_cano_bwd
        !! Map matrix to staggered mesh (phi_stag) 
        !! from backward canonical mesh (phi_stag-dphi/2)
        
        type(csrmat_t), public, allocatable :: map_cano_cano_fwd
        !! Map matrix to canonical mesh (phi_cano)
        !! from forward canonical mesh (phi_cano+dphi)

        type(csrmat_t), public, allocatable :: map_cano_cano_fwd2
        !! Map matrix to canonical mesh (phi_cano)
        !! from next but one forward canonical mesh (phi_cano+2dphi)

        type(csrmat_t), public, allocatable :: map_cano_cano_bwd
        !! Map matrix to canonical mesh (phi_cano)
        !! from backward canonical mesh (phi_cano-dphi)

        type(csrmat_t), public, allocatable :: map_cano_cano_bwd2
        !! Map matrix to canonical mesh (phi_cano)
        !! from next but one backward canonical mesh (phi_cano-2dphi)

        type(csrmat_t), public, pointer :: map_cano_stag_fwd => null()
        !! Map matrix to canonical mesh (phi_cano) 
        !! from forward staggered mesh (phi_cano+dphi/2)
        !! Points to map_cano_stag_fwd_mem for non-axisymmetric case
        !! Points to map_stag_cano_fwd for axisymmetric case
        type(csrmat_t), private, allocatable :: map_cano_stag_fwd_mem
        !! Map matrix to canonical mesh (phi_cano) 
        !! from forward staggered mesh (phi_cano+dphi/2) 
        !! Only built for non-axisymmetric case, 
        !! to be targeted by map_cano_stag_fwd

        type(csrmat_t), public, pointer :: map_cano_stag_bwd => null()
        !! Map matrix to canonical mesh (phi_cano) 
        !! from backward staggered mesh (phi_cano-dphi/2) 
        !! Points tomap_cano_stag_bwd_mem for non-axisymmetric case
        !! Points to map_stag_cano_bwd for axisymmetric case
        type(csrmat_t), private, allocatable :: map_cano_stag_bwd_mem
        !! Map matrix to canonical mesh (phi_cano) 
        !! from backward staggered mesh (phi_cano-dphi/2) 
        !! Only built for non-axisymmetric case, 
        !! to be targeted by map_cano_stag_bwd
        
        type(csrmat_t), public, pointer :: map_stag_stag_fwd => null()
        !! Map matrix to staggered mesh (phi_stag)
        !! from forward staggered mesh (phi_stag+dphi)
        !! Points to map_stag_stag_fwd_mem for non-axisymmetric case
        !! Points to map_cano_cano_fwd for axisymmetric case
        type(csrmat_t), private, allocatable :: map_stag_stag_fwd_mem
        !! Map matrix to staggered mesh (phi_stag)
        !! from forward staggered mesh (phi_stag+dphi)
        !! Only built for non-axisymmetric case, 
        !! to be targeted by map_stag_stag_fwd
        
        type(csrmat_t), public, pointer :: map_stag_stag_bwd => null()
        !! Map matrix to staggered mesh (phi_stag)
        !! from backward staggered mesh (phi_stag-dphi)
        !! Points to map_cano_cano_bwd for axisymmetric case
        type(csrmat_t), private, allocatable :: map_stag_stag_bwd_mem
        !! Map matrix to staggered mesh (phi_stag)
        !! from backward staggered mesh (phi_stag-dphi)
        !! Only built for non-axisymmetric case, 
        !! to be targeted by map_stag_stag_bwd
        
        ! >>>>>>>>>>>
        ! Operator matrices
        type(csrmat_t), public, allocatable :: grad_stag_cano_fwd
        !! Parallel gradient matrix to staggered mesh (phi_stag)
        !! from forward canonical mesh (phi_stag+dphi/2)
        
        type(csrmat_t), public, allocatable :: grad_stag_cano_bwd
        !! Parallel gradient matrix to staggered mesh (phi_stag)
        !! from backward canonical mesh (phi_stag-dphi/2)
        
        type(csrmat_t), public, allocatable :: pdiv_cano_stag_fwd
        !! Parallel divergence matrix to canonical mesh (phi_cano)
        !! from forward_staggered_mesh (phi_cano+dphi/2)
        !! ~ -fluxbox_vol_cano^(-1) * grad_stag_cano_bwd^T * fluxbox_vol_stag
        
        type(csrmat_t), public, allocatable :: pdiv_cano_stag_bwd
        !! Parallel divergence matrix to canonical mesh (phi_cano)
        !! from backward_staggered_mesh (phi_cano-dphi/2)
        !! ~ -fluxbox_vol_cano^(-1) * grad_stag_cano_fwd^T * fluxbox_vol_stag
                
#if ENABLE_CUDA
        type(cuda_map_t) :: grad_stag_cano_fwd_cuda, grad_stag_cano_bwd_cuda
        !! Parallel gradient matrix exposed to CUDA
        type(cuda_map_t) :: pdiv_cano_stag_bwd_cuda, pdiv_cano_stag_fwd_cuda
        !! Parallel divergence matrix exposed to CUDA
        type(cuda_map_t) :: map_cano_stag_bwd_cuda, map_cano_stag_fwd_cuda
        !! Map matrix exposed to CUDA
#endif

    contains   
        ! >>>>>>>>>>>
        ! Getter/Setter routines
        procedure, public :: get_nplanes
        procedure, public :: get_iplane
        procedure, public :: get_phi_cano
        procedure, public :: get_phi_stag   
        procedure, public :: get_dphi
        procedure, public :: get_intorder
        procedure, public :: set_intorder
        procedure, public :: get_xorder
        procedure, public :: set_xorder
        procedure, public :: is_fixed_stencil
        procedure, public :: is_gauss_quadrature
        procedure, public :: set_parameters => set_parameters_map
        procedure, public :: set_toroidal_location
       
        ! >>>>>>>>>>>
        ! Operator routines
        ! The single routines perform operations only for a single mesh point,
        ! the other routines on the whole mesh. These shall be called from within
        ! omp parallel regions
        procedure, public :: par_grad
        procedure, public :: par_grad_single
        procedure, public :: par_divb
        procedure, public :: par_divb_single
        procedure, public :: lift_cano_to_stag_single
        procedure, public :: lift_cano_to_stag
        procedure, public :: flat_stag_to_cano_single
        procedure, public :: flat_stag_to_cano
        procedure, public :: upstream_cano_from_cano_fwd_single
        procedure, public :: upstream_cano_from_cano_fwd
        procedure, public :: upstream_cano_from_cano_bwd_single
        procedure, public :: upstream_cano_from_cano_bwd
        procedure, public :: upstream_cano_from_cano_fwd2_single
        procedure, public :: upstream_cano_from_cano_fwd2
        procedure, public :: upstream_cano_from_cano_bwd2_single
        procedure, public :: upstream_cano_from_cano_bwd2
        procedure, public :: upstream_stag_from_stag_fwd_single
        procedure, public :: upstream_stag_from_stag_fwd
        procedure, public :: upstream_stag_from_stag_bwd_single
        procedure, public :: upstream_stag_from_stag_bwd
#if ENABLE_CUDA
        procedure, public :: par_grad_cuda
        procedure, public :: par_divb_cuda
#endif
        
        ! >>>>>>>>>>>
        ! I/O routines
        procedure, public :: write_netcdf => write_netcdf_parallel_map
        procedure, public :: read_netcdf => read_netcdf_parallel_map
        procedure, public :: display => display_parallel_map
            
        ! >>>>>>>>>>>
        ! Build routines
        procedure, public :: build => build_parallel_map
        procedure, private :: build_derived_operators
#if ENABLE_CUDA
        procedure, public :: build_cuda_struct
#endif
        final :: destructor_parallel_map_t
        
    end type parallel_map_t

    interface

        module subroutine set_parameters_map(self, comm_handler, filename, &
                                             intorder_in, xorder_in, &
                                             use_fixed_stencil_in, use_gauss_quadrature_in)
            !! Sets parameters for map, either via namelist from file, or setting explicitly
            implicit none
            class(parallel_map_t), intent(inout) :: self
            !! Instance of the type
            type(comm_handler_t), intent(in) :: comm_handler
            !! Communicators
            character(*), intent(in), optional :: filename
            !! Filename where parameter are read from
            integer, intent(in), optional :: intorder_in
            !! Interpolation order
            integer, intent(in), optional :: xorder_in
            !! Order for subcell splitting 
            logical, intent(in), optional :: use_fixed_stencil_in
            !! Indicates if fixed stencil is used for interpolation
            logical, intent(in), optional :: use_gauss_quadrature_in
            !! Indicates if gauss-quadrature or midpoint rule is used for integration
        end subroutine
        
        module subroutine set_toroidal_location(self, new_iplane, new_phi_cano, new_phi_stag)
            !! Sets toroidal location
            class(parallel_map_t), intent(inout) :: self
            !! Instance of the type
            integer, intent(in) :: new_iplane
            !! New toroidal index
            real(GP), intent(in) :: new_phi_cano
            !! New phi_cano
            real(GP), intent(in) :: new_phi_stag
            !! New phi_stag
        end subroutine
        
        module subroutine build_parallel_map(self, comm_handler, equi, mesh_cano, mesh_stag, dbgout)
            !! Builds parallel map matrices, parallel gradient  and map metadata
            class(parallel_map_t), intent(inout), target :: self
            !! Instance of the type
            type(comm_handler_t), intent(in) :: comm_handler
            !! Communicators
            class(equilibrium_t) :: equi
            !! Equilibrium (not changed)
            type(mesh_cart_t), intent(in) :: mesh_cano
            !! Mesh (canonical)
            type(mesh_cart_t), intent(in) :: mesh_stag
            !! Mesh (staggered)
            integer, intent(in), optional :: dbgout
            !! Debug output level      
        end subroutine

        module subroutine build_derived_operators(self, comm_handler)
            !! Builds derived operator, i.e. parallel divergence matrices
            class(parallel_map_t), intent(inout) :: self
            !! Instance of the type
            type(comm_handler_t), intent(in) :: comm_handler
            !! Communicators
        end subroutine

#if ENABLE_CUDA
        module subroutine build_cuda_struct(self)
            !! Exposes matrices to CUDA
            class(parallel_map_t), intent(inout), target :: self
            !! Instance of the type
        end subroutine
#endif

        module subroutine display_parallel_map(self)
            !! Displays information about parallel_map
            class(parallel_map_t), intent(in) :: self
            !! Instance of the type
        end subroutine
        
        pure real(GP) module function par_grad_single(self, ucano, ucano_fwd, l_stag)
            !! Computes parallel gradient on mesh_stag index l_stag
            class(parallel_map_t), intent(in) :: self
            !! Instance of the type
            real(GP), dimension(self%grad_stag_cano_bwd%ncol), intent(in) :: ucano
            !! Field on mesh_cano at phi_cano
            real(GP), dimension(self%grad_stag_cano_fwd%ncol), intent(in) :: ucano_fwd
            !! Field on mesh_cano_fwd at phi_cano+dphi
            integer, intent(in) :: l_stag
            !! Index of mesh_stag at phi_stag, where parallel gradient 
            !! is requested
        end function

        module subroutine par_grad(self, ucano, ucano_fwd, par_grad_u_stag)
            !! Computes parallel gradient on mesh_stag
            class(parallel_map_t), intent(in) :: self
            !! Instance of the type
            real(GP), dimension(self%grad_stag_cano_bwd%ncol), intent(in) :: ucano
            !! Field on mesh_cano at phi_cano
            real(GP), dimension(self%grad_stag_cano_fwd%ncol), intent(in) :: ucano_fwd
            !! Field on mesh_cano_fwd at phi_cano+dphi
            real(GP), dimension(self%grad_stag_cano_fwd%ndim), intent(out) :: par_grad_u_stag
            !! Parallel gradient on mesh_stag at phi_stag
        end subroutine
        
        pure real(GP) module function par_divb_single(self, ustag, ustag_bwd, l_cano)
            !! Computes parallel divergence on mesh_cano at index l_cano
            class(parallel_map_t), intent(in) :: self
            !! Instance of the type
            real(GP), dimension(self%pdiv_cano_stag_fwd%ncol), intent(in) :: ustag
            !! Field on mesh_stag at phi_stag
            real(GP), dimension(self%pdiv_cano_stag_bwd%ncol), intent(in) :: ustag_bwd
            !! Field on mesh_stag_bwd at phi_cano-dphi
            integer, intent(in) :: l_cano
            !! Index of mesh_cano at phi_cano, where parallel divergence 
            !! is requested
        end function
        
        module subroutine par_divb(self, ustag, ustag_bwd, par_pdiv_u_cano)
            !! Computes parallel divergence on mesh_cano
            class(parallel_map_t), intent(in) :: self
            !! Instance of the type
            real(GP), dimension(self%pdiv_cano_stag_fwd%ncol), intent(in) :: ustag
            !! Field on mesh_stag at phi_stag
            real(GP), dimension(self%pdiv_cano_stag_bwd%ncol), intent(in) :: ustag_bwd
            !! Field on mesh_stag_bwd at phi_stag-dphi
            real(GP), dimension(self%pdiv_cano_stag_fwd%ndim), intent(out) :: par_pdiv_u_cano
            !! Parallel divergence on mesh_cano at phi_cano
        end subroutine
        
        pure real(GP) module function lift_cano_to_stag_single(self, ucano, ucano_fwd, l_stag)
            !! Lifts (maps) quantity from mesh_cano to mesh_stag at index l_stag
            class(parallel_map_t), intent(in) :: self
            !! Instance of the type
            real(GP), dimension(self%map_stag_cano_bwd%ncol), intent(in) :: ucano
            !! Field on mesh_cano at phi_cano
            real(GP), dimension(self%map_stag_cano_fwd%ncol), intent(in) :: ucano_fwd
            !! Field on mesh_cano_fwd at phi_cano+dphi
            integer, intent(in) :: l_stag
            !! Index of mesh_stag at phi_stag, where lift is requested
        end function

        module subroutine lift_cano_to_stag(self, ucano, ucano_fwd, lift_u_stag)
            !! Lifts (maps) quantity from mesh_cano to mesh_stag
            class(parallel_map_t), intent(in) :: self
            !! Instance of the type
            real(GP), dimension(self%map_stag_cano_bwd%ncol), intent(in) :: ucano
            !! Field on mesh_cano at phi_cano
            real(GP), dimension(self%map_stag_cano_fwd%ncol), intent(in) :: ucano_fwd
            !! Field on mesh_cano_fwd at phi_cano+dphi
            real(GP), dimension(self%map_stag_cano_fwd%ndim), intent(out) :: lift_u_stag
            !! Quantity lifted to mesh_stag at phi_stag
        end subroutine
        
        pure real(GP) module function flat_stag_to_cano_single(self, ustag, ustag_bwd, l_cano)
            !! Flats (maps) quantity from mesh_stag to mesh_cano at index l_cano
            class(parallel_map_t), intent(in) :: self
            !! Instance of the type
            real(GP), dimension(self%pdiv_cano_stag_fwd%ncol), intent(in) :: ustag
            !! Field on mesh_stag at phi_stag
            real(GP), dimension(self%pdiv_cano_stag_bwd%ncol), intent(in) :: ustag_bwd
            !! Field on mesh_stag_bwd at phi_cano-dphi
            integer, intent(in) :: l_cano
            !! Index of mesh_cano at phi_cano, where flat is requested
        end function
        
        module subroutine flat_stag_to_cano(self, ustag, ustag_bwd, flat_u_cano)
            !! Flats (maps) quantity from mesh_stag to mesh_cano
            class(parallel_map_t), intent(in) :: self
            !! Instance of the type
            real(GP), dimension(self%map_cano_stag_fwd%ncol), intent(in) :: ustag
            !! Field on mesh_stag at phi_stag
            real(GP), dimension(self%map_cano_stag_bwd%ncol), intent(in) :: ustag_bwd
            !! Field on mesh_stag_bwd at phi_stag-dphi
            real(GP), dimension(self%map_cano_stag_fwd%ndim), intent(out) :: flat_u_cano
            !! Quantity flatted to mesh_cano at phi_cano
        end subroutine
        
        pure real(GP) module function upstream_cano_from_cano_fwd_single(self, ucano_fwd, l_cano)
            !! Maps upstream values to mesh_cano from mesh_cano_fwd at phi_cano+dphi at index l_cano
            class(parallel_map_t), intent(in) :: self
            !! Instance of the type
            real(GP), dimension(self%map_cano_cano_fwd%ncol), intent(in) :: ucano_fwd
            !! Field on mesh_cano_fwd at phi_cano+dphi
            integer, intent(in) :: l_cano
            !! Index of mesh_cano, where value is requested
        end function
        
        module subroutine upstream_cano_from_cano_fwd(self, ucano_fwd, upstream_cano_fwd)
            !! Maps upstream values to mesh_cano from mesh_cano_fwd at phi_cano+dphi
            class(parallel_map_t), intent(in) :: self
            !! Instance of the type
            real(GP), dimension(self%map_cano_cano_fwd%ncol), intent(in) :: ucano_fwd
            !! Field on mesh_cano_fwd at phi_cano+dphi
            real(GP), dimension(self%map_cano_cano_fwd%ndim), intent(out) :: upstream_cano_fwd
            !! Upstream qauntity mapped to mesh_cano
        end subroutine

        pure real(GP) module function upstream_cano_from_cano_fwd2_single(self, ucano_fwd2, l_cano)
            !! Maps upstream values to mesh_cano from mesh_cano_fwd2 at phi_cano+2dphi at index l_cano
            class(parallel_map_t), intent(in) :: self
            !! Instance of the type
            real(GP), dimension(self%map_cano_cano_fwd2%ncol), intent(in) :: ucano_fwd2
            !! Field on mesh_cano_fwd2 at phi_cano+2dphi
            integer, intent(in) :: l_cano
            !! Index of mesh_cano, where value is requested
        end function

        module subroutine upstream_cano_from_cano_fwd2(self, ucano_fwd2, upstream_cano_fwd2)
            !! Maps upstream values to mesh_cano from mesh_cano_fwd2 at phi_cano+2dphi
            class(parallel_map_t), intent(in) :: self
            !! Instance of the type
            real(GP), dimension(self%map_cano_cano_fwd2%ncol), intent(in) :: ucano_fwd2
            !! Field on mesh_cano_fwd2 at phi_cano+2dphi
            real(GP), dimension(self%map_cano_cano_fwd2%ndim), intent(out) :: upstream_cano_fwd2
            !! Upstream qauntity mapped to mesh_cano
        end subroutine

        pure real(GP) module function upstream_cano_from_cano_bwd_single(self, ucano_bwd, l_cano)
            !! Maps upstream values to mesh_cano from mesh_cano_bwd at phi_cano-dphi at index l_cano
            class(parallel_map_t), intent(in) :: self
            !! Instance of the type
            real(GP), dimension(self%map_cano_cano_bwd%ncol), intent(in) :: ucano_bwd
            !! Field on mesh_cano_bwd at phi_cano-dphi
            integer, intent(in) :: l_cano
            !! Index of mesh_cano, where value is requested
        end function

        module subroutine upstream_cano_from_cano_bwd(self, ucano_bwd, upstream_cano_bwd)
            !! Maps upstream values to mesh_cano from mesh_cano_bwd at phi_cano-dphi
            class(parallel_map_t), intent(in) :: self
            !! Instance of the type
            real(GP), dimension(self%map_cano_cano_bwd%ncol), intent(in) :: ucano_bwd
            !! Field on mesh_cano_bwd at phi_cano-dphi
            real(GP), dimension(self%map_cano_cano_bwd%ndim), intent(out) :: upstream_cano_bwd
            !! Upstream qauntity mapped to mesh_cano
        end subroutine

        pure real(GP) module function upstream_cano_from_cano_bwd2_single(self, ucano_bwd2, l_cano)
            !! Maps upstream values to mesh_cano from mesh_cano_bwd2 at phi_cano-2dphi at index l_cano
            class(parallel_map_t), intent(in) :: self
            !! Instance of the type
            real(GP), dimension(self%map_cano_cano_bwd2%ncol), intent(in) :: ucano_bwd2
            !! Field on mesh_cano_bwd2 at phi_cano-2dphi
            integer, intent(in) :: l_cano
            !! Index of mesh_cano, where value is requested
        end function

        module subroutine upstream_cano_from_cano_bwd2(self, ucano_bwd2, upstream_cano_bwd2)
            !! Maps upstream values to mesh_cano from mesh_cano_bwd at phi_cano-dphi
            class(parallel_map_t), intent(in) :: self
            !! Instance of the type
            real(GP), dimension(self%map_cano_cano_bwd2%ncol), intent(in) :: ucano_bwd2
            !! Field on mesh_cano_bwd2 at phi_cano-2dphi
            real(GP), dimension(self%map_cano_cano_bwd2%ndim), intent(out) :: upstream_cano_bwd2
            !! Upstream qauntity mapped to mesh_cano
        end subroutine

        pure real(GP) module function upstream_stag_from_stag_fwd_single(self, ustag_fwd, l_stag)
            !! Maps upstream values to mesh_stag from mesh_stag_fwd at phi_stag+dphi at index l_stag
            class(parallel_map_t), intent(in) :: self
            !! Instance of the type
            real(GP), dimension(self%map_stag_stag_fwd%ncol), intent(in) :: ustag_fwd
            !! Field on mesh_stag_fwd at phi_stag+dphi
            integer, intent(in) :: l_stag
            !! Index of mesh_stag, where value is requested
        end function
        
        module subroutine upstream_stag_from_stag_fwd(self, ustag_fwd, upstream_stag_fwd)
            !! Maps upstream values to mesh_stag from mesh_stag_fwd at phi_stag+dphi
            class(parallel_map_t), intent(in) :: self
            !! Instance of the type
            real(GP), dimension(self%map_stag_stag_fwd%ncol), intent(in) :: ustag_fwd
            !! Field on mesh_stag_fwd at phi_stag+dphi
            real(GP), dimension(self%map_stag_stag_fwd%ndim), intent(out) :: upstream_stag_fwd
            !! Upstream qauntity mapped to mesh_stag
        end subroutine
  
        pure real(GP) module function upstream_stag_from_stag_bwd_single(self, ustag_bwd, l_stag)
            !! Maps upstream values to mesh_stag from mesh_stag_bwd at phi_stag-dphi at index l_stag
            class(parallel_map_t), intent(in) :: self
            !! Instance of the type
            real(GP), dimension(self%map_stag_stag_bwd%ncol), intent(in) :: ustag_bwd
            !! Field on mesh_stag_bwd at phi_stag-dphi
            integer, intent(in) :: l_stag
            !! Index of mesh_stag, where value is requested
        end function
        
        module subroutine upstream_stag_from_stag_bwd(self, ustag_bwd, upstream_stag_bwd)
            !! Maps upstream values to mesh_stag from mesh_stag_bwd at phi_stag-dphi
            class(parallel_map_t), intent(in) :: self
            !! Instance of the type
            real(GP), dimension(self%map_stag_stag_bwd%ncol), intent(in) :: ustag_bwd
            !! Field on mesh_stag_bwd at phi_stag-dphi
            real(GP), dimension(self%map_stag_stag_bwd%ndim), intent(out) :: upstream_stag_bwd
            !! Upstream qauntity mapped to mesh_stag
        end subroutine

        module subroutine par_grad_cuda(self, ucano, ucano_fwd, par_grad_u_stag)
            !! Computes parallel gradient on mesh_stag,
            !! based on CUDA implementation, not to be called from OMP parallel region
            class(parallel_map_t), intent(in) :: self
            !! Instance of the type
            real(GP), dimension(self%grad_stag_cano_bwd%ncol), intent(in) :: ucano
            !! Field on mesh_cano at phi_cano
            real(GP), dimension(self%grad_stag_cano_fwd%ncol), intent(in) :: ucano_fwd
            !! Field on mesh_cano_fwd at phi_cano+dphi
            real(GP), dimension(self%grad_stag_cano_fwd%ndim), intent(out) :: par_grad_u_stag
            !! Parallel gradient on mesh_stag at phi_stag
        end subroutine
                
        module subroutine par_divb_cuda(self, ustag, ustag_bwd, par_pdiv_u_cano)
            !! Computes parallel divergence on mesh_cano,
            !! based on CUDA implementation, not to be called from OMP parallel region
            class(parallel_map_t), intent(in) :: self
            !! Instance of the type
            real(GP), dimension(self%pdiv_cano_stag_fwd%ncol), intent(in) :: ustag
            !! Field on mesh_stag at phi_stag
            real(GP), dimension(self%pdiv_cano_stag_bwd%ncol), intent(in) :: ustag_bwd
            !! Field on mesh_stag_bwd at phi_stag-dphi
            real(GP), dimension(self%pdiv_cano_stag_fwd%ndim), intent(out) :: par_pdiv_u_cano
            !! Parallel divergence on mesh_cano at phi_cano
        end subroutine
       
        module subroutine write_netcdf_parallel_map(self, axisymmetric, fgid)
            !! Writes parallel map to netcdf file
            class(parallel_map_t), intent(in) :: self
            !! Instance of the type
            logical, intent(in) :: axisymmetric
            !! Indicates if output is written based on axisymmetric equilibria,
            !! writes out only non-redundant data
            integer, intent(in) :: fgid
            !! netcdf file or group id
        end subroutine

        module subroutine read_netcdf_parallel_map(self, axisymmetric, fgid)
            !! Reads parallel map to netcdf file
            class(parallel_map_t), intent(inout), target :: self
            !! Instance of the type
            logical, intent(in) :: axisymmetric
            !! Indicates if input is read based on axisymmetric equilibria
            integer, intent(in) :: fgid
            !! netcdf file or group id
        end subroutine

        module subroutine destructor_parallel_map_t(self)
            !! Frees memory associated with parallel_map
            type(parallel_map_t), intent(inout) :: self
            !! Instance of the type
        end subroutine

    end interface
    
contains

    pure integer function get_nplanes(self)
        !! Returns number of poloidal planes used for map
        class(parallel_map_t), intent(in) :: self
        !! Instance of the type
        get_nplanes = self%nplanes
    end function
    
    pure integer function get_iplane(self)
        !! Returns index of poloidal plane (0 to nplanes-1)
        class(parallel_map_t), intent(in) :: self
        !! Instance of the type
        get_iplane = self%iplane  
    end function

    pure real(GP) function get_dphi(self)
        !! Returns toroidal grid distance used for map
        class(parallel_map_t), intent(in) :: self
        !! Instance of the type
        get_dphi = self%dphi  
    end function

    real(GP) function get_phi_cano(self)
        !! Returns phi coordinate of current plane
        class(parallel_map_t), intent(in) :: self
        !! Instance of the type
        get_phi_cano = self%phi_cano
    end function

    real(GP) function get_phi_stag(self)
        !! Returns phi coordinate of current staggered plane
        class(parallel_map_t), intent(in) :: self
        !! Instance of the type
        get_phi_stag = self%phi_stag
    end function

    pure integer function get_intorder(self)
        !! Returns interpolation order used for map
        class(parallel_map_t), intent(in) :: self
        !! Instance of the type
        get_intorder = self%intorder
    end function

    subroutine set_intorder(self, intorder)
        !! Set interpolation order
        class(parallel_map_t), intent(inout) :: self
        !! Instance of the type
        integer, intent(in) :: intorder
        self%intorder = intorder
    end subroutine

    pure integer function get_xorder(self)
        !! Returns order of subcell splitting used for map
        class(parallel_map_t), intent(in) :: self
        !! Instance of the type
        get_xorder = self%xorder
    end function  

    subroutine set_xorder(self, xorder)
        !! Set order of subcell splitting
        class(parallel_map_t), intent(inout) :: self
        !! Instance of the type
        integer, intent(in) :: xorder
        self%xorder = xorder
    end subroutine

    pure logical function is_fixed_stencil(self)
        !! Returns if fixed stencil is used for interpolation
        class(parallel_map_t), intent(in) :: self
        !! Instance of the type
        is_fixed_stencil = self%use_fixed_stencil
    end function
    
    pure logical function is_gauss_quadrature(self)
        !! Returns if gauss-quadrature or midpoint rule is used for integration
        class(parallel_map_t), intent(in) :: self
        !! Instance of the type
        is_gauss_quadrature = self%use_gauss_quadrature
    end function 
    
end module