module polar_map_factory_m !! Creates matrix mapping from Cartesian to Polar grid use precision_m, only : FP use screen_io_m, only : get_stdout use comm_handling_m, only : is_master use csrmat_m, only : csrmat_t use equilibrium_m, only : equilibrium_t use mesh_cart_m, only : mesh_cart_t use coords_polar_m, only : polar_to_cart, cart_to_polar use polar_grid_m, only : polar_grid_t use interpolation_m, only : interpol_coeffs use descriptors_m, only : DISTRICT_CLOSED, DISTRICT_CORE implicit none public :: create_polar_map_matrix public :: create_cart_map_matrix contains subroutine create_polar_map_matrix(polar_map_csr, equi, mesh, polar_grid, intorder, dbgout) !! Creates matrix that maps a field from Cartesian to polar mesh type(csrmat_t), allocatable, intent(out) :: polar_map_csr !! Created polar map matrix in CSR format class(equilibrium_t), intent(inout) :: equi !! Equilibrium type(mesh_cart_t), intent(in) :: mesh !! Cartesian mesh type(polar_grid_t), intent(in) :: polar_grid !! Polar grid integer, intent(in) :: intorder !! Interpolation order of the mapping integer, intent(in), optional :: dbgout !! Specifies the number of information printed on screen integer, parameter :: ns_max_nearest = 6 !! Size of area that is searched for nearest point integer :: outi integer :: nrho, ntheta, nz_alloc, nz, ip, jp, k, i, j integer :: ns_stencil, intorder_actual, ind_nearest, lsw real(FP) :: rho, theta, phi, xp, yp real(FP) :: xsw, ysw, xrel, yrel integer, allocatable, dimension(:,:) :: ind_interpol_stencil real(FP), allocatable, dimension(:,:) :: coeffs_interpol ! set output level outi = 0 if (present(dbgout)) then if (is_master()) then outi = dbgout endif if (dbgout >= 3) then outi = dbgout ! every rank writes endif endif ! Print to console if (is_master() .and. outi >= 1) then write(get_stdout(), *) '...................................' write(get_stdout(), *) 'Building a cart->polar map matrix' write(get_stdout(), 100) intorder 100 FORMAT( & X,'intorder =',I4, / & ) endif phi = polar_grid%get_phi() ! Allocate matrix nrho = polar_grid%get_nrho() ntheta = polar_grid%get_ntheta() allocate(polar_map_csr) polar_map_csr%ndim = nrho*ntheta polar_map_csr%ncol = mesh%get_n_points() polar_map_csr%nnz = 0 nz_alloc = nrho*ntheta * (intorder+1)**2 allocate(polar_map_csr%i(polar_map_csr%ndim+1)) allocate(polar_map_csr%j(nz_alloc)) allocate(polar_map_csr%val(nz_alloc)) nz = 0 do ip = 1, nrho do jp = 1, ntheta k = (ip-1)*ntheta+jp polar_map_csr%i(k) = nz+1 rho = polar_grid%get_rho(ip) theta = polar_grid%get_theta(jp) call polar_to_cart(equi, rho, theta, phi, xp, yp) ! Find interpolation stencil ---------------------------------- ! find size of largest possible interpolation stencil -> ns_stencil call mesh%get_surrounding_indices(xp = xp, & yp = yp, & ns = intorder + 1, & ns_complete = ns_stencil ) ! determine final interpolation order and interpolation stencil if (ns_stencil > 0) then intorder_actual = ns_stencil - 1 allocate( ind_interpol_stencil(ns_stencil, ns_stencil) ) call mesh%get_surrounding_indices(xp = xp, & yp = yp, & ns = ns_stencil, & ind_sur = ind_interpol_stencil) elseif (ns_stencil == 0) then ! find nearest point in 6x6 surrounding grid call mesh%get_surrounding_indices(xp = xp, & yp = yp, & ns = ns_max_nearest, & ind_nearest = ind_nearest) if (ind_nearest /= 0) then ! nearest neighbour intorder_actual = 0 ns_stencil = 1 allocate( ind_interpol_stencil(ns_stencil, ns_stencil) ) ind_interpol_stencil(1,1) = ind_nearest else ! no interpolation at all intorder_actual = -1 ns_stencil = 0 allocate( ind_interpol_stencil(ns_stencil, ns_stencil) ) endif endif ! Compute interpolation coefficients -------------------------- allocate( coeffs_interpol(ns_stencil, ns_stencil) ) if (intorder_actual >= 1) then ! polynomial interpolation of order cintorer ! get south-west point of interpolation stamp lsw = ind_interpol_stencil(1,1) xsw = mesh%get_x(lsw) ysw = mesh%get_y(lsw) ! get interpolation coefficients xrel = (xp-xsw)/mesh%get_spacing_f() yrel = (yp-ysw)/mesh%get_spacing_f() call interpol_coeffs(intorder_actual, xrel, yrel, coeffs_interpol) elseif (intorder_actual == 0) then ! nearest neighbour interpolation coeffs_interpol(1, 1) = 1.0_FP endif ! fill matrix ------------------------------------------------- do j = 1, ns_stencil do i = 1, ns_stencil if (ind_interpol_stencil(i, j) /= 0) then nz = nz+1 polar_map_csr%j(nz) = ind_interpol_stencil(i, j) polar_map_csr%val(nz) = coeffs_interpol(i,j) endif enddo enddo deallocate(coeffs_interpol) deallocate(ind_interpol_stencil) enddo enddo ! Finalize matrix polar_map_csr%i(nrho*ntheta+1) = nz+1 polar_map_csr%nnz = polar_map_csr%i(nrho*ntheta+1)-1 call polar_map_csr%sortsum() if (is_master() .and. outi >= 1) then write(get_stdout(), *) & 'communication of cart->polar map matrix done, nz=', & polar_map_csr%nnz, 'nz_al=', nz_alloc write(get_stdout(), *) '...................................' endif end subroutine subroutine create_cart_map_matrix(cart_map_csr, equi, mesh, polar_grid, & intorder, dbgout) !! Creates matrix that maps a field from polar mesh to Cartesian mesh !! (only for closed field line region) type(csrmat_t), allocatable, intent(out) :: cart_map_csr !! Created Cartesian map matrix in CSR format class(equilibrium_t), intent(inout) :: equi !! Equilibrium type(mesh_cart_t), intent(in) :: mesh !! Cartesian mesh type(polar_grid_t), intent(in) :: polar_grid !! Polar grid integer, intent(in) :: intorder !! Interpolation order of the mapping integer, intent(in), optional :: dbgout !! Specifies the number of information printed on screen integer :: l, nrho, ntheta, nz, outi, & district, in, jn, kpol, ks, st_hsize real(FP) :: x, y, rho, theta, phi, frho, ftheta, xn, yn integer, dimension(intorder+1) :: st_ip, st_jp real(FP), dimension(intorder+1, intorder+1) :: coeffs ! set output level outi = 0 if (present(dbgout)) then if (is_master()) then outi = dbgout endif if (dbgout >= 3) then outi = dbgout ! every rank writes endif endif ! Print to console if (is_master() .and. outi >= 1) then write(get_stdout(), *) '...................................' write(get_stdout(), *) 'Building a polar->cart map matrix' write(get_stdout(), 100) intorder 100 FORMAT( & X,'intorder =',I4, / & ) endif phi = mesh%get_phi() ! Allocate matrix nrho = polar_grid%get_nrho() ntheta = polar_grid%get_ntheta() allocate(cart_map_csr) cart_map_csr%ndim = mesh%get_n_points() cart_map_csr%ncol = nrho*ntheta cart_map_csr%nnz = 0 ! Compute numbers of nonzero entries nz = 0 do l = 1, mesh%get_n_points() district = mesh%get_district(l) if ((district == DISTRICT_CLOSED) .or. & (district == DISTRICT_CORE)) then nz = nz + (intorder + 1)**2 endif enddo allocate(cart_map_csr%i(cart_map_csr%ndim+1)) allocate(cart_map_csr%j(nz)) allocate(cart_map_csr%val(nz), source=0.0_FP) ! Half size of interpolation stencil st_hsize = (intorder + 1) / 2 nz = 0 do l = 1, mesh%get_n_points() cart_map_csr%i(l) = nz + 1 district = mesh%get_district(l) ! Skip points in SOL and private flux region if ((district == DISTRICT_CLOSED) .or. & (district == DISTRICT_CORE)) then x = mesh%get_x(l) y = mesh%get_y(l) call cart_to_polar(equi, x, y, phi, rho, theta) ! Find indices of interpolation stencil in rho direction >>>>>> frho = (rho - polar_grid%get_rhopol_min()) & / polar_grid%get_drho() ! Coordinate of interpolation point in ! normalised stencil units xn = frho - floor(frho) + 1.0_FP * (st_hsize - 1) do ks = 1, intorder + 1 st_ip(ks) = floor(frho) - st_hsize + 1 + ks enddo ! Move interpolation stencil into valid region for left edge do while (st_ip(1) <=0 ) st_ip(:) = st_ip(:) + 1 xn = xn - 1.0_FP enddo ! Move interpolation stencil into valid region for right edge do while (st_ip(intorder+1) > polar_grid%get_nrho()) st_ip(:) = st_ip(:) - 1 xn = xn + 1.0_FP enddo ! Find indices of interpolation stencil in theta direction ftheta = theta / polar_grid%get_dtheta() yn = ftheta - floor(ftheta) + 1.0_FP * (st_hsize - 1) do ks = 1, intorder + 1 st_jp(ks) = floor(ftheta) - st_hsize + 1 + ks ! Deal with periodicity if (st_jp(ks) <= 0) then st_jp(ks) = st_jp(ks) + polar_grid%get_ntheta() endif if (st_jp(ks) > polar_grid%get_ntheta()) then st_jp(ks) = st_jp(ks) - polar_grid%get_ntheta() endif enddo call interpol_coeffs(intorder, xn, yn, coeffs) do in = 1, intorder + 1 do jn = 1, intorder + 1 nz = nz + 1 kpol = (st_ip(in)-1) * polar_grid%get_ntheta() + & st_jp(jn) cart_map_csr%j(nz) = kpol cart_map_csr%val(nz) = coeffs(in, jn) enddo enddo endif enddo ! Finalize matrix cart_map_csr%nnz = nz cart_map_csr%i(mesh%get_n_points()+1) = nz + 1 call cart_map_csr%sortsum() if (is_master() .and. outi >= 1) then write(get_stdout(), *) & 'build of polar->cart map matrix done, nz=', cart_map_csr%nnz write(get_stdout(), *) '...................................' endif end subroutine end module