module zonal_averages_factory_m !! Creates matrices to compute polar averages !! - flux-surface average (which is actually a volume average) !! - surface average (an actual surface average, e.g. for computation of fluxes) !! TODO: Not tested on 3D equilibria - phi variable only used for !! equilibrium and polar grid interfaces 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, jacobian_polar use polar_grid_m, only : polar_grid_t use list_operations_m, only : sort_and_sum implicit none public :: create_flux_surface_average_csr public :: create_surface_average_csr public :: create_flux_surface_average_csr_viaaddvols private :: create_zonal_average_csr contains subroutine create_flux_surface_average_csr(flux_surface_average_csr, equi, mesh, polar_grid, polar_map_csr, dbgout) !! Creates a matrix to compute flux surface average from Cartesian field ! wrapper around create_zonal_average_csr implicit none type(csrmat_t), allocatable, intent(out) :: flux_surface_average_csr !! Zonal average matrix 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 type(csrmat_t), intent(in) :: polar_map_csr !! Polar map matrix in CSR format integer, intent(in), optional :: dbgout !! Specifies the number of information printed on screen call create_zonal_average_csr(flux_surface_average_csr, equi, mesh, polar_grid, polar_map_csr, 2, dbgout) end subroutine subroutine create_surface_average_csr(surface_average_csr, equi, mesh, polar_grid, polar_map_csr, dbgout) !! Creates a matrix to surface average from Cartesian field ! wrapper around create_zonal_average_csr implicit none type(csrmat_t), allocatable, intent(out) :: surface_average_csr !! Zonal average matrix 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 type(csrmat_t), intent(in) :: polar_map_csr !! Polar map matrix in CSR format integer, intent(in), optional :: dbgout !! Specifies the number of information printed on screen call create_zonal_average_csr(surface_average_csr, equi, mesh, polar_grid, polar_map_csr, 1, dbgout) end subroutine subroutine create_zonal_average_csr(zonal_average_csr, equi, mesh, polar_grid, polar_map_csr, average_mode, dbgout) !! Creates a matrix to compute zonal averages from Cartesian field !! Builds zonal averages based on polar_map_csr implicit none type(csrmat_t), allocatable, intent(out) :: zonal_average_csr !! Zonal average matrix 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 type(csrmat_t), intent(in) :: polar_map_csr !! Polar map matrix in CSR format integer, intent(in) :: average_mode !! Indicates type of average, !! average_mode = 1: matrix for polar average is established (actually surface average) !! average_mode = 2: matrix for flux surface average is established (actually volume average) integer, intent(in), optional :: dbgout !! Specifies the number of information printed on screen integer, parameter :: nz_alfac = 16 integer :: outi integer :: nrho, ntheta, nz_alloc, nz, ip, jp, k, nzwrk, kwrk real(FP) :: phi, rho, theta, xp, yp, drho, dtheta real(FP) :: xl1, xl2, yl1, yl2, lpol, ltor, fac, normsrf integer, allocatable, dimension(:) :: jwrk real(FP), allocatable, dimension(:) :: valwrk ! 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 zonal_average matrix' write(get_stdout(), 100) average_mode 100 FORMAT( & X,'average_mode =',I4, / & ) endif phi = polar_grid%get_phi() nrho = polar_grid%get_nrho() ntheta = polar_grid%get_ntheta() drho = polar_grid%get_drho() dtheta = polar_grid%get_dtheta() ! Allocate matrix allocate(zonal_average_csr) zonal_average_csr%ndim = nrho zonal_average_csr%ncol = polar_map_csr%ncol zonal_average_csr%nnz = 0 nz_alloc = polar_map_csr%nnz allocate(zonal_average_csr%i(zonal_average_csr%ndim+1)) allocate(zonal_average_csr%j(nz_alloc)) allocate(zonal_average_csr%val(nz_alloc)) ! Allocate work space allocate(jwrk(nz_alfac*ntheta)) allocate(valwrk(nz_alfac*ntheta)) ! Build matrix nz = 0 do ip = 1, nrho zonal_average_csr%i(ip) = nz+1 nzwrk = 0 normsrf = 0.0_FP do jp=1,ntheta rho = polar_grid%get_rho(ip) theta = polar_grid%get_theta(jp) call polar_to_cart(equi, rho, theta, phi, xp, yp) ! Compute averging factor k = (ip-1)*ntheta+jp if (average_mode == 1) then ! polar average matrix (surface) call polar_to_cart(equi, rho, theta-dtheta/2.0_FP, phi, xl1, yl1) call polar_to_cart(equi, rho, theta+dtheta/2.0_FP, phi, xl2, yl2) lpol = sqrt( (xl1-xl2)**2 + (yl1-yl2)**2 ) ltor = equi%jacobian(xp, yp, phi) fac = lpol*ltor elseif (average_mode == 2) then ! flux surface average matrix (volume) fac = drho*dtheta*jacobian_polar(equi, xp, yp, phi) endif normsrf = normsrf + fac ! Store result in workspace do kwrk = polar_map_csr%i(k), polar_map_csr%i(k+1)-1 nzwrk = nzwrk+1 jwrk(nzwrk) = polar_map_csr%j(kwrk) valwrk(nzwrk) = fac*polar_map_csr%val(kwrk) enddo enddo ! Set and sort matrix call sort_and_sum(nzwrk, jwrk, valwrk) do kwrk = 1, nzwrk zonal_average_csr%j(nz+kwrk) = jwrk(kwrk) zonal_average_csr%val(nz+kwrk)= valwrk(kwrk) / normsrf enddo nz = nz + nzwrk enddo ! Finalise matrix zonal_average_csr%i(zonal_average_csr%ndim+1) = nz + 1 zonal_average_csr%nnz = nz if (is_master() .and. outi >= 1) then write(get_stdout(), *) & 'zonal matrix created, nz=', zonal_average_csr%nnz, & 'nz_al=', nz_alloc write(get_stdout(), *) '...................................' endif end subroutine subroutine create_flux_surface_average_csr_viaaddvols(flux_surface_average_csr, equi, mesh, polar_grid, dbgout) !! Constructs flux surface average matrix via adding volumes of adjacent surfaces implicit none type(csrmat_t), allocatable, intent(out) :: flux_surface_average_csr !! Zonal average matrix 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), optional :: dbgout !! Specifies the number of information printed on screen integer :: outi integer :: nrho, nz_alloc, nz, ip, l, k real(FP) :: phi, drho, vrho, rho_s_min, rho_s_max, x, y, rho ! 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 flux surface average matrix & &via adding finite volumes' endif phi = polar_grid%get_phi() nrho = polar_grid%get_nrho() drho = polar_grid%get_drho() ! Allocate matrix allocate(flux_surface_average_csr) flux_surface_average_csr%ndim = nrho flux_surface_average_csr%ncol = mesh%get_n_points() flux_surface_average_csr%nnz = 0 nz_alloc = mesh%get_n_points() allocate(flux_surface_average_csr%i(flux_surface_average_csr%ndim+1)) allocate(flux_surface_average_csr%j(nz_alloc)) allocate(flux_surface_average_csr%val(nz_alloc)) ! Build matrix nz = 0 do ip = 1, nrho flux_surface_average_csr%i(ip) = nz+1 ! boundaries and volume of flux surface vrho = 0.0_FP rho_s_min = polar_grid%get_rho(ip) - drho/2.0_FP rho_s_max = polar_grid%get_rho(ip) + drho/2.0_FP ! If mesh point within boundaries add entry into matrix do l = 1, mesh%get_n_points() x = mesh%get_x(l) y = mesh%get_y(l) rho = equi%rho(x, y, phi) if ((rho > rho_s_min) .and. (rho <= rho_s_max)) then nz = nz+1 flux_surface_average_csr%j(nz) = l flux_surface_average_csr%val(nz) = equi%jacobian(x, y, phi) vrho = vrho + equi%jacobian(x, y, phi) endif enddo ! Divide by volume of flux surface do k = flux_surface_average_csr%i(ip), nz flux_surface_average_csr%val(k) = flux_surface_average_csr%val(k)/vrho enddo enddo ! Finalise matrix flux_surface_average_csr%i(flux_surface_average_csr%ndim+1) = nz + 1 flux_surface_average_csr%nnz = nz if (is_master() .and. outi >= 1) then write(get_stdout(), *) & 'flux surface matrix created, nz=', & flux_surface_average_csr%nnz, 'nz_al=', nz_alloc write(get_stdout(), *) '...................................' endif end subroutine end module