parallel_map_build_s.f90 Source File


Contents


Source Code

submodule (parallel_map_m) parallel_map_build_s
    !! Rotines related with build of map
    implicit none

contains

    module subroutine build_parallel_map(self, comm_handler, equi, mesh_cano, mesh_stag, dbgout)
        class(parallel_map_t), intent(inout), target :: self
        type(comm_handler_t),  intent(in) :: comm_handler
        class(equilibrium_t) :: equi
        type(mesh_cart_t), intent(in) :: mesh_cano
        type(mesh_cart_t), intent(in) :: mesh_stag
        integer, intent(in), optional :: dbgout    
    
        integer ::  l, iz, n_points, n_points_stag, outi
        real(GP) :: phi_cano, phi_stag, dphi, dphi_half, x, y
        real(GP) :: fluxexp_fwd, fluxexp_bwd, vb
        
        type(mesh_cart_t), allocatable :: mesh_target
        
        self%nplanes = comm_handler%get_nplanes()
        self%iplane = comm_handler%get_plane()
        self%phi_cano = mesh_cano%get_phi()
        self%phi_stag = mesh_stag%get_phi()
        self%dphi = 2.0_GP * (self%phi_stag - self%phi_cano)
        
        dphi = self%dphi
        dphi_half = dphi/2.0_GP
        phi_cano = self%phi_cano
        phi_stag = self%phi_stag 

        ! set output level
        outi = 0
        if (present(dbgout)) then
            if (is_rank_info_writer) then
                outi = dbgout
            endif
            if (dbgout >= 3) then
                outi = dbgout ! every rank writes
            endif
        endif
             
        
        ! Build Map matrices to stag from cano
        call perf_start('../map_build_matrices')
        
        allocate(mesh_target)
        call mesh_cano%exchange_mesh_mpi(comm_handler%get_comm_planes(), &
                                         step=1, mesh_received=mesh_target)
        
        allocate(self%map_stag_cano_fwd)                  
        call create_map_matrix(map_matrix = self%map_stag_cano_fwd,     &
                                comm = comm_handler%get_comm_planes(),  &
                                equi = equi,                            &
                                mesh_base = mesh_stag,                  &
                                mesh_target = mesh_target,              &
                                phi_direction = +1,                     &
                                intorder =  self%intorder,              &
                                xorder   =  0,                          &
                                flux_box_surface_integral = .false.,    &
                                use_fixed_stencil = self%use_fixed_stencil, &
                                use_gauss_quadrature = self%use_gauss_quadrature, &
                                dbgout = outi                           )
        deallocate(mesh_target)
               
        allocate(mesh_target)
        call mesh_cano%exchange_mesh_mpi(comm_handler%get_comm_planes(), &
                                         step=0, mesh_received=mesh_target)
        allocate(self%map_stag_cano_bwd)
        call create_map_matrix(map_matrix = self%map_stag_cano_bwd,     &
                                comm = comm_handler%get_comm_planes(),  &
                                equi = equi,                            &
                                mesh_base = mesh_stag,                  &
                                mesh_target = mesh_target,              &
                                phi_direction = -1,                     &
                                intorder =  self%intorder,              &
                                xorder   =  0,                          &
                                flux_box_surface_integral = .false.,    &
                                use_fixed_stencil = self%use_fixed_stencil, &
                                use_gauss_quadrature = self%use_gauss_quadrature, &
                                dbgout = outi                           )
        deallocate(mesh_target)
        
        ! Build map matrices to cano from stag
        if (equi%is_axisymmetric()) then
            self%map_cano_stag_fwd => self%map_stag_cano_fwd
            self%map_cano_stag_bwd => self%map_stag_cano_bwd
        else
            allocate(mesh_target)
            call mesh_stag%exchange_mesh_mpi(comm_handler%get_comm_planes(), &
                                             step=0, mesh_received=mesh_target)
            allocate(self%map_cano_stag_fwd_mem)
            self%map_cano_stag_fwd => self%map_cano_stag_fwd_mem  
            call create_map_matrix(map_matrix = self%map_cano_stag_fwd,     &
                                    comm = comm_handler%get_comm_planes(),  &
                                    equi = equi,                            &
                                    mesh_base = mesh_cano,                  &
                                    mesh_target = mesh_target,              &
                                    phi_direction = +1,                     &
                                    intorder =  self%intorder,              &
                                    xorder   =  0,                          &
                                    flux_box_surface_integral = .false.,    &
                                    use_fixed_stencil = self%use_fixed_stencil, &
                                    use_gauss_quadrature = self%use_gauss_quadrature, &
                                    dbgout = outi                           )
            deallocate(mesh_target)
            
            allocate(mesh_target)
            call mesh_stag%exchange_mesh_mpi(comm_handler%get_comm_planes(), &
                                             step=-1, mesh_received=mesh_target)
            allocate(self%map_cano_stag_bwd_mem)
            self%map_cano_stag_bwd => self%map_cano_stag_bwd_mem
            call create_map_matrix(map_matrix = self%map_cano_stag_bwd,     &
                                    comm = comm_handler%get_comm_planes(),  &
                                    equi = equi,                            &
                                    mesh_base = mesh_cano,                  &
                                    mesh_target = mesh_target,              &
                                    phi_direction = -1,                     &
                                    intorder =  self%intorder,              &
                                    xorder   =  0,                          &
                                    flux_box_surface_integral = .false.,    &
                                    use_fixed_stencil = self%use_fixed_stencil, &
                                    use_gauss_quadrature = self%use_gauss_quadrature, &
                                    dbgout = outi                           )
            deallocate(mesh_target)       
        endif                         
        
        ! Build parallel gradient matrices to stag from cano
        allocate(mesh_target)
        call mesh_cano%exchange_mesh_mpi(comm_handler%get_comm_planes(), &
                                         step=+1, mesh_received=mesh_target)
        allocate(self%grad_stag_cano_fwd)
        call create_map_matrix(map_matrix = self%grad_stag_cano_fwd,    &
                                comm = comm_handler%get_comm_planes(),  &
                                equi = equi,                            &
                                mesh_base = mesh_stag,                  &
                                mesh_target = mesh_target,              &
                                phi_direction = +1,                     &
                                intorder =  self%intorder,              &
                                xorder   =  self%xorder,                &
                                flux_box_surface_integral = .true.,     &
                                use_fixed_stencil = self%use_fixed_stencil, &
                                use_gauss_quadrature = self%use_gauss_quadrature, &
                                dbgout = outi                           )
        deallocate(mesh_target)
        
        allocate(mesh_target)
        call mesh_cano%exchange_mesh_mpi(comm_handler%get_comm_planes(), &
                                         step=0, mesh_received=mesh_target)
        allocate(self%grad_stag_cano_bwd)
        call create_map_matrix(map_matrix = self%grad_stag_cano_bwd,    &
                                comm = comm_handler%get_comm_planes(),  &
                                equi = equi,                            &
                                mesh_base = mesh_stag,                  &
                                mesh_target = mesh_target,              &
                                phi_direction = -1,                     &
                                intorder =  self%intorder,              &
                                xorder   =  self%xorder,                &
                                flux_box_surface_integral = .true.,     &
                                use_fixed_stencil = self%use_fixed_stencil, &
                                use_gauss_quadrature = self%use_gauss_quadrature, &
                                dbgout = outi                           )
        deallocate(mesh_target)

        ! Build map matrices to cano from cano
        allocate(mesh_target)
        call mesh_cano%exchange_mesh_mpi(comm_handler%get_comm_planes(), &
                                    step=+1, mesh_received=mesh_target)
        allocate(self%map_cano_cano_fwd)
        call create_map_matrix(map_matrix = self%map_cano_cano_fwd,     &
                                comm = comm_handler%get_comm_planes(),  &
                                equi = equi,                            &
                                mesh_base = mesh_cano,                  &
                                mesh_target = mesh_target,              &
                                phi_direction = +1,                     &
                                intorder =  self%intorder,              &
                                xorder   =  0,                          &
                                flux_box_surface_integral = .false.,    &
                                use_fixed_stencil = self%use_fixed_stencil, &
                                use_gauss_quadrature = self%use_gauss_quadrature, &
                                dbgout = outi                           )
        deallocate(mesh_target)
        
        allocate(mesh_target)
        call mesh_cano%exchange_mesh_mpi(comm_handler%get_comm_planes(), &
                                    step=-1, mesh_received=mesh_target)
        allocate(self%map_cano_cano_bwd)
        call create_map_matrix(map_matrix = self%map_cano_cano_bwd,     &
                                comm = comm_handler%get_comm_planes(),  &
                                equi = equi,                            &
                                mesh_base = mesh_cano,                  &
                                mesh_target = mesh_target,              &
                                phi_direction = -1,                     &
                                intorder =  self%intorder,              &
                                xorder   =  0,                          &
                                flux_box_surface_integral = .false.,    &
                                use_fixed_stencil = self%use_fixed_stencil, &
                                use_gauss_quadrature = self%use_gauss_quadrature, &
                                dbgout = outi                           )
        deallocate(mesh_target)

        ! Build map matrices to cano from next but one cano
        allocate(mesh_target)
        call mesh_cano%exchange_mesh_mpi(comm_handler%get_comm_planes(), &
                                    step=+2, mesh_received=mesh_target)
        allocate(self%map_cano_cano_fwd2)
        call create_map_matrix(map_matrix = self%map_cano_cano_fwd2,    &
                                comm = comm_handler%get_comm_planes(),  &
                                equi = equi,                            &
                                mesh_base = mesh_cano,                  &
                                mesh_target = mesh_target,              &
                                phi_direction = +1,                     &
                                intorder =  self%intorder,              &
                                xorder   =  0,                          &
                                flux_box_surface_integral = .false.,    &
                                use_fixed_stencil = self%use_fixed_stencil, &
                                use_gauss_quadrature = self%use_gauss_quadrature, &
                                dbgout = outi                           )
        deallocate(mesh_target)

        allocate(mesh_target)
        call mesh_cano%exchange_mesh_mpi(comm_handler%get_comm_planes(), &
                                    step=-2, mesh_received=mesh_target)
        allocate(self%map_cano_cano_bwd2)
        call create_map_matrix(map_matrix = self%map_cano_cano_bwd2,    &
                                comm = comm_handler%get_comm_planes(),  &
                                equi = equi,                            &
                                mesh_base = mesh_cano,                  &
                                mesh_target = mesh_target,              &
                                phi_direction = -1,                     &
                                intorder =  self%intorder,              &
                                xorder   =  0,                          &
                                flux_box_surface_integral = .false.,    &
                                use_fixed_stencil = self%use_fixed_stencil, &
                                use_gauss_quadrature = self%use_gauss_quadrature, &
                                dbgout = outi                           )
        deallocate(mesh_target)

        ! Build map matrices to stag from stag
        if (equi%is_axisymmetric()) then
            self%map_stag_stag_fwd => self%map_cano_cano_fwd
            self%map_stag_stag_bwd => self%map_cano_cano_bwd
        else
            allocate(mesh_target)
            call mesh_stag%exchange_mesh_mpi(comm_handler%get_comm_planes(), &
                                        step=+1, mesh_received=mesh_target)
            allocate(self%map_stag_stag_fwd_mem)
            self%map_stag_stag_fwd => self%map_stag_stag_fwd_mem
            call create_map_matrix(map_matrix = self%map_stag_stag_fwd,     &
                                    comm = comm_handler%get_comm_planes(),  &
                                    equi = equi,                            &
                                    mesh_base = mesh_stag,                  &
                                    mesh_target = mesh_target,              &
                                    phi_direction = +1,                     &
                                    intorder =  self%intorder,              &
                                    xorder   =  0,                          &
                                    flux_box_surface_integral = .false.,    &
                                    use_fixed_stencil = self%use_fixed_stencil, &
                                    use_gauss_quadrature = self%use_gauss_quadrature, &
                                    dbgout = outi                           )
            deallocate(mesh_target)
            
            allocate(mesh_target)
            call mesh_stag%exchange_mesh_mpi(comm_handler%get_comm_planes(), &
                                        step=-1, mesh_received=mesh_target)
            allocate(self%map_stag_stag_bwd_mem)
            self%map_stag_stag_bwd => self%map_stag_stag_bwd_mem
            call create_map_matrix(map_matrix = self%map_stag_stag_bwd,     &
                                    comm = comm_handler%get_comm_planes(),  &
                                    equi = equi,                            &
                                    mesh_base = mesh_stag,                  &
                                    mesh_target = mesh_target,              &
                                    phi_direction = -1,                     &
                                    intorder =  self%intorder,              &
                                    xorder   =  0,                          &
                                    flux_box_surface_integral = .false.,    &
                                    use_fixed_stencil = self%use_fixed_stencil, &
                                    use_gauss_quadrature = self%use_gauss_quadrature, &
                                    dbgout = outi                           )
            deallocate(mesh_target)
        endif
        
        call perf_stop('../map_build_matrices')

        call perf_start('../map_build_meta')
        ! Compute parallel distances and flux box volumes for canonical mesh
        n_points = mesh_cano%get_n_points()

        allocate(self%dpar_cano_stag_fwd(n_points))  
        allocate(self%dpar_cano_stag_bwd(n_points))
        allocate(self%dpar_cano_cano_fwd(n_points))  
        allocate(self%dpar_cano_cano_bwd(n_points))
        allocate(self%fluxbox_vol_cano(n_points)) 

        !$omp parallel &
        !$omp default(none) &
        !$omp shared(n_points, self, phi_cano, dphi_half, dphi, equi, mesh_cano) &
        !$omp private(l, x, y, fluxexp_fwd, fluxexp_bwd)
        !$omp do
        do l = 1, n_points
            x = mesh_cano%get_x(l)
            y = mesh_cano%get_y(l)
            
            call trace(x_init = x, &
                y_init = y, &
                phi_init = phi_cano, &
                dphi = dphi_half, &
                equi = equi, &
                arclen = self%dpar_cano_stag_fwd(l), &
                fluxexp = fluxexp_fwd)

            call trace(x_init = x, &
                y_init = y, &
                phi_init = phi_cano, &
                dphi = -dphi_half, &
                equi = equi, &
                arclen = self%dpar_cano_stag_bwd(l), &
                fluxexp = fluxexp_bwd)
            
            self%fluxbox_vol_cano(l) = abs( equi%btor(x, y, phi_cano) * &
                        mesh_cano%get_spacing_f()**2 * (fluxexp_fwd + fluxexp_bwd) )
            
            call trace(x_init = x, &
                y_init = y, &
                phi_init = phi_cano, &
                dphi = dphi, &
                equi = equi, &
                arclen = self%dpar_cano_cano_fwd(l))

            call trace(x_init = x, &
                y_init = y, &
                phi_init = phi_cano, &
                dphi = -dphi, &
                equi = equi, &
                arclen = self%dpar_cano_cano_bwd(l))

        enddo
        !$omp end do
        !$omp end parallel
        
        ! Parallel distances and flux box volumes for staggerd mesh
        n_points_stag = mesh_stag%get_n_points()
        if (equi%is_axisymmetric()) then
            ! For axisymmetric case quantities for staggered mesh 
            ! are the same asfor canonical mesh
            self%dpar_stag_cano_fwd => self%dpar_cano_stag_fwd
            self%dpar_stag_cano_bwd => self%dpar_cano_stag_bwd
            self%dpar_stag_stag_fwd => self%dpar_cano_cano_fwd
            self%dpar_stag_stag_bwd => self%dpar_cano_cano_bwd
            self%fluxbox_vol_stag => self%fluxbox_vol_cano
        else
            ! For non axisymetric case compute quantities for staggered mesh
            allocate(self%dpar_stag_cano_fwd_mem(n_points_stag))  
            self%dpar_stag_cano_fwd => self%dpar_stag_cano_fwd_mem
            
            allocate(self%dpar_stag_cano_bwd_mem(n_points_stag))
            self%dpar_stag_cano_bwd => self%dpar_stag_cano_bwd_mem
            
            allocate(self%dpar_stag_stag_fwd_mem(n_points_stag)) 
            self%dpar_stag_stag_fwd => self%dpar_stag_stag_fwd_mem
             
            allocate(self%dpar_stag_stag_bwd_mem(n_points_stag))
            self%dpar_stag_stag_bwd => self%dpar_stag_stag_bwd_mem
            
            allocate(self%fluxbox_vol_stag_mem(n_points_stag))
            self%fluxbox_vol_stag => self%fluxbox_vol_stag_mem
            !$omp parallel &
            !$omp default(none) &
            !$omp shared(n_points_stag, self, phi_stag, dphi_half, dphi, equi, mesh_stag) &
            !$omp private(l, x, y, fluxexp_fwd, fluxexp_bwd)
            !$omp do        
            do l = 1, n_points_stag           
                x = mesh_stag%get_x(l)
                y = mesh_stag%get_y(l)
               
                call trace(x_init = x, &
                    y_init = y, &
                    phi_init = phi_stag, &
                    dphi = dphi_half, &
                    equi = equi, &
                    arclen = self%dpar_stag_cano_fwd(l), &
                    fluxexp = fluxexp_fwd)
                    
                call trace(x_init = x, &
                    y_init = y, &
                    phi_init = phi_stag, &
                    dphi = -dphi_half, &
                    equi = equi, &
                    arclen = self%dpar_stag_cano_bwd(l), &
                    fluxexp = fluxexp_bwd)
                
                self%fluxbox_vol_stag(l) = abs( equi%btor(x, y, phi_stag) * &
                                mesh_stag%get_spacing_f()**2 * (fluxexp_fwd + fluxexp_bwd) )
                
                call trace(x_init = x, &
                    y_init = y, &
                    phi_init = phi_stag, &
                    dphi = dphi, &
                    equi = equi, &
                    arclen = self%dpar_stag_stag_fwd(l))
                
                call trace(x_init = x, &
                    y_init = y, &
                    phi_init = phi_stag, &
                    dphi = -dphi, &
                    equi = equi, &
                    arclen = self%dpar_stag_stag_bwd(l))
            enddo
            !$omp end do
            !$omp end parallel
        endif
        
        ! Divide matrices Q by vb
        !$omp parallel &
        !$omp default(none) &
        !$omp shared(n_points_stag, self, phi_stag, equi, mesh_stag) &
        !$omp private(l, x, y, vb, iz)
        !$omp do        
        do l = 1, n_points_stag  
            x = mesh_stag%get_x(l)
            y = mesh_stag%get_y(l)
            vb = (self%dpar_stag_cano_fwd(l) + self%dpar_stag_cano_bwd(l)) &
                 * equi%btor(x, y, phi_stag) ! Length along field line           
            do iz = self%grad_stag_cano_fwd%i(l), self%grad_stag_cano_fwd%i(l+1)-1
                self%grad_stag_cano_fwd%val(iz) = self%grad_stag_cano_fwd%val(iz) / vb  
            enddo 
            do iz = self%grad_stag_cano_bwd%i(l), self%grad_stag_cano_bwd%i(l+1)-1
                self%grad_stag_cano_bwd%val(iz) = self%grad_stag_cano_bwd%val(iz) / vb  
            enddo 
        enddo
        !$omp end do
        !$omp end parallel
        
        ! Build derived operators, i.e. parallel divergence matrices pdiv_cano_stag_fwd and pdiv_cano_stag_bwd
        call self%build_derived_operators(comm_handler)
        
#if ENABLE_CUDA
        call self%build_cuda_struct()
#endif
        
        call perf_stop('../map_build_meta')
        
        if (outi >0) then
            call self%display()
        endif

    end subroutine

    module subroutine build_derived_operators(self, comm_handler)
        class(parallel_map_t), intent(inout) :: self
        type(comm_handler_t),  intent(in) :: comm_handler
        
        integer :: n_points, n_points_stag, n_points_stag_bwd, l, iz, j
        real(GP), allocatable, dimension(:) :: fluxbox_vol_stag_bwd
        type(csrmat_t), allocatable :: wrk_csr

        ! Compute parallel divergence matrices 
        ! pdiv_cano_stag_fwd = -V^(-1)*grad_stag_cano_bwd^T*V_stag
        ! pdiv_cano_stag_bwd = V^(-1)*grad_stag_cano_fwd^T*V_stag_bwd 

        ! Transposition
        call getdata_fwdbwdplane(comm_handler%get_comm_planes(), step=-1, &
                                 acsr_send=self%grad_stag_cano_fwd, &
                                 acsr_recv=wrk_csr)
        call csr_transpose(wrk_csr, self%pdiv_cano_stag_bwd)
        deallocate(wrk_csr)
        call csr_transpose(self%grad_stag_cano_bwd, self%pdiv_cano_stag_fwd)
        
        n_points = self%pdiv_cano_stag_fwd%ndim
        n_points_stag = self%pdiv_cano_stag_fwd%ncol
        

        ! Get flux box volumes from backward staggered mesh
        call getdata_fwdbwdplane(comm_handler%get_comm_planes(), step=-1, &
                                 nsend=n_points_stag, &
                                 usend=self%fluxbox_vol_stag, &
                                 nrecv=n_points_stag_bwd, &
                                 urecv=fluxbox_vol_stag_bwd)
        
        ! Correct with flux box volume factors
        !$omp parallel &
        !$omp default(none) &
        !$omp shared(n_points, self, fluxbox_vol_stag_bwd) &
        !$omp private(l, iz, j)
        !$omp do        
        do l = 1, n_points
            do iz = self%pdiv_cano_stag_fwd%i(l), self%pdiv_cano_stag_fwd%i(l+1)-1
                j = self%pdiv_cano_stag_fwd%j(iz)
                self%pdiv_cano_stag_fwd%val(iz) = -self%pdiv_cano_stag_fwd%val(iz) * self%fluxbox_vol_stag(j) / self%fluxbox_vol_cano(l)                
            enddo 
            do iz = self%pdiv_cano_stag_bwd%i(l), self%pdiv_cano_stag_bwd%i(l+1)-1
                j = self%pdiv_cano_stag_bwd%j(iz)
                self%pdiv_cano_stag_bwd%val(iz) = -self%pdiv_cano_stag_bwd%val(iz) * fluxbox_vol_stag_bwd(j) / self%fluxbox_vol_cano(l)                
            enddo         
        enddo
        !$omp end do
        !$omp end parallel
        
        deallocate(fluxbox_vol_stag_bwd)

    end subroutine

    module subroutine destructor_parallel_map_t(self)
        type(parallel_map_t), intent(inout) :: self
        
        if (associated(self%map_stag_stag_bwd)) nullify(self%map_stag_stag_bwd)
        if (allocated(self%map_stag_stag_bwd_mem)) deallocate(self%map_stag_stag_bwd_mem)
        
        if (associated(self%map_stag_stag_fwd)) nullify(self%map_stag_stag_fwd)
        if (allocated(self%map_stag_stag_fwd_mem)) deallocate(self%map_stag_stag_fwd_mem)
        
        if (allocated(self%map_cano_cano_bwd)) deallocate(self%map_cano_cano_bwd)
        
        if (allocated(self%map_cano_cano_fwd)) deallocate(self%map_cano_cano_fwd)
        
        if (allocated(self%pdiv_cano_stag_bwd)) deallocate(self%pdiv_cano_stag_bwd)
        
        if (allocated(self%pdiv_cano_stag_fwd)) deallocate(self%pdiv_cano_stag_fwd)
        
        if (allocated(self%grad_stag_cano_bwd)) deallocate(self%grad_stag_cano_bwd)
        
        if (allocated(self%grad_stag_cano_fwd)) deallocate(self%grad_stag_cano_fwd)
        
        if (associated(self%map_cano_stag_bwd)) nullify(self%map_cano_stag_bwd)
        if (allocated(self%map_cano_stag_bwd_mem)) deallocate(self%map_cano_stag_bwd_mem)
        
        if (associated(self%map_cano_stag_fwd)) nullify(self%map_cano_stag_fwd)
        if (allocated(self%map_cano_stag_fwd_mem)) deallocate(self%map_cano_stag_fwd_mem)
        
        if (allocated(self%map_stag_cano_bwd)) deallocate(self%map_stag_cano_bwd)
        
        if (allocated(self%map_stag_cano_fwd)) deallocate(self%map_stag_cano_fwd)

        if (associated(self%fluxbox_vol_stag)) nullify(self%fluxbox_vol_stag)
        if (allocated(self%fluxbox_vol_stag_mem)) deallocate(self%fluxbox_vol_stag_mem)
        
        if (associated(self%dpar_stag_stag_bwd)) nullify(self%dpar_stag_stag_bwd)
        if (allocated(self%dpar_stag_stag_bwd_mem)) deallocate(self%dpar_stag_stag_bwd_mem)
        
        if (associated(self%dpar_stag_stag_fwd)) nullify(self%dpar_stag_stag_fwd)  
        if (allocated(self%dpar_stag_stag_fwd_mem)) deallocate(self%dpar_stag_stag_fwd_mem)
        
        if (associated(self%dpar_stag_cano_bwd)) nullify(self%dpar_stag_cano_bwd)
        if (allocated(self%dpar_stag_cano_bwd_mem)) deallocate(self%dpar_stag_cano_bwd_mem)
        
        if (associated(self%dpar_stag_cano_fwd)) nullify(self%dpar_stag_cano_fwd)
        if (allocated(self%dpar_stag_cano_fwd_mem)) deallocate(self%dpar_stag_cano_fwd_mem)
          
        if (allocated(self%fluxbox_vol_cano)) deallocate(self%fluxbox_vol_cano)
        
        if (allocated(self%dpar_cano_cano_bwd)) deallocate(self%dpar_cano_cano_bwd)
        
        if (allocated(self%dpar_cano_cano_fwd)) deallocate(self%dpar_cano_cano_fwd)
        
        if (allocated(self%dpar_cano_stag_bwd)) deallocate(self%dpar_cano_stag_bwd)
        
        if (allocated(self%dpar_cano_stag_fwd)) deallocate(self%dpar_cano_stag_fwd)  

    end subroutine destructor_parallel_map_t

#if ENABLE_CUDA

module subroutine build_cuda_struct(self)
        class(parallel_map_t), intent(inout), target :: self

        call allocate_cuda_struct(self%grad_stag_cano_fwd_cuda, self%grad_stag_cano_fwd%ndim, &
          self%grad_stag_cano_fwd%nnz, self%grad_stag_cano_fwd%i, self%grad_stag_cano_fwd%j, &
          self%grad_stag_cano_fwd%val)
        call allocate_cuda_struct(self%grad_stag_cano_bwd_cuda, self%grad_stag_cano_bwd%ndim, &
          self%grad_stag_cano_bwd%nnz, self%grad_stag_cano_bwd%i, self%grad_stag_cano_bwd%j, &
          self%grad_stag_cano_bwd%val)
        call allocate_cuda_struct(self%pdiv_cano_stag_bwd_cuda, self%pdiv_cano_stag_bwd%ndim, &
          self%pdiv_cano_stag_bwd%nnz, self%pdiv_cano_stag_bwd%i, self%pdiv_cano_stag_bwd%j, &
          self%pdiv_cano_stag_bwd%val)
        call allocate_cuda_struct(self%pdiv_cano_stag_fwd_cuda, self%pdiv_cano_stag_fwd%ndim, &
          self%pdiv_cano_stag_fwd%nnz, self%pdiv_cano_stag_fwd%i, self%pdiv_cano_stag_fwd%j, &
          self%pdiv_cano_stag_fwd%val)
        call allocate_cuda_struct(self%map_cano_stag_bwd_cuda, self%map_cano_stag_bwd%ndim, &
          self%map_cano_stag_bwd%nnz, self%map_cano_stag_bwd%i, self%map_cano_stag_bwd%j, &
          self%map_cano_stag_bwd%val)
        call allocate_cuda_struct(self%map_cano_stag_fwd_cuda, self%map_cano_stag_fwd%ndim, &
          self%map_cano_stag_fwd%nnz, self%map_cano_stag_fwd%i, self%map_cano_stag_fwd%j, &
          self%map_cano_stag_fwd%val)

    end subroutine

    subroutine allocate_cuda_struct(cuda_map, nrows, nnz, rows, cols, vals)
        !! Exposes CSR matrix to CUDA
        type(cuda_map_t), intent(inout) :: cuda_map
        !! Contains exposure information of CSR matrix
        integer, intent(in) :: nrows
        !! Number of rows of matrix
        integer, intent(in) :: nnz
        !! Number of non-zeros of matrix
        integer, dimension(nrows+1), intent(in) :: rows
        !! Row array of CSR matrix
        integer, dimension(nnz), intent(in) :: cols
        !! Column indices of CSR matrix
        real(GP), dimension(nnz), intent(in) :: vals
        !! Values of CSR matrix

        cuda_map%nrows = nrows
        cuda_map%nnz = nnz
        call allocate_memory_int(cuda_map%cols, 1, cuda_map%nnz)
        call allocate_memory_int(cuda_map%rows, 1, cuda_map%nrows+1)
        call allocate_memory_double(cuda_map%vals, 1, cuda_map%nnz)
        call allocate_memory_double(cuda_map%x, 1, cuda_map%nrows)
        call allocate_memory_double(cuda_map%y, 1, cuda_map%nrows)
        
        cuda_map%rows = rows - 1 
        cuda_map%cols = cols - 1 
        cuda_map%vals = vals

        call allocate_cusparse_struct(cuda_map%cuda_struct, cuda_map%nrows, &
                cuda_map%nrows, cuda_map%nnz, cuda_map%rows, &
                cuda_map%cols, cuda_map%vals, cuda_map%x, cuda_map%y)
  
    end subroutine

    module subroutine destructor_cuda_map_t(self)
        type(cuda_map_t), intent(inout) :: self

        call freeint(self%rows)
        call freeint(self%cols)
        call freedouble(self%vals)
        call freedouble(self%x)
        call freedouble(self%y)
        
        self%cuda_struct = c_null_ptr

    end subroutine

#endif

end submodule