module diagnostics_helpers_m !! Module containing useful diagnostic computation routines use MPI use precision_grillix_m, only : GP, MPI_GP, GP_EPS use comm_handler_m, only : comm_handler_t use error_handling_grillix_m, only: handle_error, error_info_t use status_codes_grillix_m, only : GRILLIX_ERR_OTHER use mesh_cart_m, only: mesh_cart_t use parallel_map_m, only : parallel_map_t use penalisation_m, only : penalisation_t use interpolation_m, only : interpol_coeffs use equilibrium_m, only: equilibrium_t implicit none contains real(GP) function inner_product(comm_handler, mesh, map, penalisation, & is_stag, u_in, v_in, & use_abs, use_vol_avg, use_full_domain, sum_over_planes) !! Compute inner product of u*v over all volume elements !! Only accepts input arrays of lengths n_points or n_points_inner type(comm_handler_t), intent(in) :: comm_handler !! Communicator type(mesh_cart_t), intent(in) :: mesh !! Mesh within poloidal plane type(parallel_map_t), intent(in) :: map !! Parallel map type(penalisation_t), intent(in) :: penalisation !! Penalisation logical, intent(in) :: is_stag !! Wheter computation is performed on staggered grid real(GP), dimension(:), intent(in) :: u_in !! First input quantity real(GP), dimension(:), intent(in), optional :: v_in !! Second input quantity logical, intent(in), optional :: use_abs !! Whether to use absolue u,v values in product (default: false) logical, intent(in), optional :: use_vol_avg !! Whether to divide result by total fluxbox volume (default: false) logical, intent(in), optional :: use_full_domain !! Whether to include outer mesh points and penalisation region !! in the computation (default: false) logical, intent(in), optional :: sum_over_planes !! Whether to sum over all planes (default: true) integer :: l, ki, kb, kg, ierr real(GP), dimension(:), allocatable :: u real(GP), dimension(:), allocatable :: v real(GP) :: inner_product_plane logical :: apply_abs = .false. logical :: apply_avg = .false. logical :: full_domain = .false. logical :: apply_sum_planes = .true. ! Set optional switches if (present(use_abs)) then apply_abs = use_abs endif if (present(use_vol_avg)) then apply_avg = use_vol_avg endif if (present(use_full_domain)) then full_domain = use_full_domain endif if (present(sum_over_planes)) then apply_sum_planes = sum_over_planes endif ! Assert that input array length is either n_points or n_points_inner if (size(u_in) /= mesh%get_n_points() .and. & size(u_in) /= mesh%get_n_points_inner()) then call handle_error('Input size must be n_points or n_points_inner', & GRILLIX_ERR_OTHER, __LINE__, __FILE__, & error_info_t( & 'size(u_in), n_points, n_points_inner: ', & [size(u_in), mesh%get_n_points(), & mesh%get_n_points_inner()])) endif if (.not. allocated(u)) allocate(u, source=u_in) if (present(v_in)) then ! Assert equal sized inputs if (size(u_in) /= size(v_in)) then call handle_error('Dimension mismatch', & GRILLIX_ERR_OTHER, __LINE__, __FILE__, & error_info_t('size(u_in), size(v_in): ', & [size(u_in), size(v_in)])) endif if (.not. allocated(v)) allocate(v, source=v_in) else ! If no v provided, set it to 1 everywhere if (.not. allocated(v)) allocate(v(size(u)), source=1.0_GP) endif if (apply_abs) then !$omp parallel default(none) private(l) shared(u, v) !$omp do do l = 1, size(u) u(l) = abs(u(l)) v(l) = abs(v(l)) enddo !$omp end do !$omp end parallel endif inner_product_plane = 0.0_GP ! Main computation loop if (size(u) == mesh%get_n_points_inner()) then ! Input vector contains only inner points !$omp parallel default(none) private(l, ki) & !$omp shared(mesh, full_domain, penalisation, is_stag, u, v, map) & !$omp reduction(+:inner_product_plane) !$omp do do ki = 1, mesh%get_n_points_inner() if ( (.not. full_domain) & .and. (penalisation%get_charfun_val(ki) > GP_EPS) ) then ! Skip points in penalisation region cycle else l = mesh%inner_indices(ki) if (is_stag) then inner_product_plane = inner_product_plane & + u(ki) * v(ki) * map%fluxbox_vol_stag(l) else inner_product_plane = inner_product_plane & + u(ki) * v(ki) * map%fluxbox_vol_cano(l) endif endif enddo !$omp end do !$omp end parallel else ! Input vector contains all points !$omp parallel default(none) private(l, ki) & !$omp shared(mesh, full_domain, penalisation, is_stag, u, v, map) & !$omp reduction(+:inner_product_plane) !$omp do do ki = 1, mesh%get_n_points_inner() if ( (.not. full_domain) & .and. (penalisation%get_charfun_val(ki) > GP_EPS) ) then ! Skip points in penalisation region cycle else l = mesh%inner_indices(ki) if (is_stag) then inner_product_plane = inner_product_plane & + u(l) * v(l) * map%fluxbox_vol_stag(l) else inner_product_plane = inner_product_plane & + u(l) * v(l) * map%fluxbox_vol_cano(l) endif endif enddo !$omp end do !$omp end parallel if (full_domain) then ! Treatment of outer points !$omp parallel default(none) private(l, kb, kg) & !$omp shared(mesh, is_stag, u, v, map) & !$omp reduction(+:inner_product_plane) !$omp do do kb = 1, mesh%get_n_points_boundary() l = mesh%boundary_indices(kb) if (is_stag) then inner_product_plane = inner_product_plane & + u(l) * v(l) * map%fluxbox_vol_stag(l) else inner_product_plane = inner_product_plane & + u(l) * v(l) * map%fluxbox_vol_cano(l) endif enddo !$omp end do !$omp do do kg = 1, mesh%get_n_points_ghost() l = mesh%ghost_indices(kg) if (is_stag) then inner_product_plane = inner_product_plane & + u(l) * v(l) * map%fluxbox_vol_stag(l) else inner_product_plane = inner_product_plane & + u(l) * v(l) * map%fluxbox_vol_cano(l) endif enddo !$omp end do !$omp end parallel endif endif if (apply_sum_planes) then call MPI_allreduce(inner_product_plane, inner_product, 1, MPI_GP, & MPI_SUM, comm_handler%get_comm_planes(), ierr) else inner_product = inner_product_plane endif if (apply_avg) then inner_product = inner_product / total_volume(comm_handler, mesh, & map, penalisation, is_stag, & full_domain, apply_sum_planes) endif if (allocated(u)) deallocate(u) if (allocated(v)) deallocate(v) end function real(GP) function total_volume(comm_handler, mesh, map, penalisation, is_stag, use_full_domain, sum_over_planes) !! Compute total volume over all mesh points type(comm_handler_t), intent(in) :: comm_handler !! Communicator type(mesh_cart_t), intent(in) :: mesh !! Mesh within poloidal plane type(parallel_map_t), intent(in) :: map !! Parallel map type(penalisation_t), intent(in) :: penalisation !! Penalisation logical, intent(in) :: is_stag !! Wheter computation is performed on staggered grid logical, intent(in), optional :: use_full_domain !! Whether to include outer mesh points and penalisation region !! in the computation (default: false) logical, intent(in), optional :: sum_over_planes !! Whether to sum over all planes (default: true) integer :: l, ki, kb, kg, ierr real(GP) :: total_volume_plane logical :: full_domain = .false. logical :: apply_sum_planes = .true. ! Set optional switch if (present(use_full_domain)) then full_domain = use_full_domain endif if (present(sum_over_planes)) then apply_sum_planes = sum_over_planes endif total_volume_plane = 0.0_GP !$omp parallel default(none) private(l, ki) & !$omp shared(mesh, full_domain, penalisation, is_stag, map) & !$omp reduction(+:total_volume_plane) !$omp do do ki = 1, mesh%get_n_points_inner() if ( (.not. full_domain) & .and. (penalisation%get_charfun_val(ki) > GP_EPS) ) then ! Skip points in penalisation region cycle else l = mesh%inner_indices(ki) if (is_stag) then total_volume_plane = total_volume_plane + map%fluxbox_vol_stag(l) else total_volume_plane = total_volume_plane + map%fluxbox_vol_cano(l) endif endif enddo !$omp end do !$omp end parallel if (full_domain) then !$omp parallel default(none) private(l, kb, kg) & !$omp shared(mesh, is_stag, map) & !$omp reduction(+:total_volume_plane) !$omp do do kb = 1, mesh%get_n_points_boundary() l = mesh%boundary_indices(kb) if (is_stag) then total_volume_plane = total_volume_plane + map%fluxbox_vol_stag(l) else total_volume_plane = total_volume_plane + map%fluxbox_vol_cano(l) endif enddo !$omp end do !$omp do do kg = 1, mesh%get_n_points_ghost() l = mesh%ghost_indices(kg) if (is_stag) then total_volume_plane = total_volume_plane + map%fluxbox_vol_stag(l) else total_volume_plane = total_volume_plane + map%fluxbox_vol_cano(l) endif enddo !$omp end do !$omp end parallel endif if (apply_sum_planes) then call MPI_allreduce(total_volume_plane, total_volume, 1, MPI_GP, & MPI_SUM, comm_handler%get_comm_planes(), ierr) else total_volume = total_volume_plane endif end function total_volume subroutine interpolate_lineout(mesh, n_points, int_order, x_coords, y_coords, & u_in, u_out) !! Compute interpolation along a lineout within poloidal plane type(mesh_cart_t), intent(in) :: mesh !! Mesh within poloidal plane integer, intent(in) :: n_points !! Number of points in lineout to interpolate integer, intent(in) :: int_order !! Interpolation order real(GP), dimension(n_points), intent(in) :: x_coords !! x-coordinates of lineout points to interpolate real(GP), dimension(n_points), intent(in) :: y_coords !! y-coordinates of lineout points to interpolate real(GP), dimension(mesh%get_n_points()), intent(in) :: u_in !! Input data given on entire plane real(GP), dimension(n_points), intent(out) :: u_out !! Interpolated values on the lineout points integer :: i, j, l, lsw, ns_stencil integer :: intorder_actual integer, parameter :: ns_max_nearest = 6 real(GP) :: xp, yp, xsw, ysw, xrel, yrel, u_p integer, dimension(:,:), allocatable :: ind_stencil real(GP), dimension(:,:), allocatable :: coeffs_interpol ! Iterate over points in lineout do l = 1, n_points xp = x_coords(l) yp = y_coords(l) call mesh%determine_interpolation_stencil(xp, yp, & int_order, ns_max_nearest, ns_stencil, intorder_actual, & ind_stencil) allocate(coeffs_interpol(ns_stencil, ns_stencil) ) ! Compute interpolation coefficents if (intorder_actual >= 1) then ! Polynomial interpolation of order intorder_actual lsw = ind_stencil(1, 1) xsw = mesh%get_x(lsw) ysw = mesh%get_y(lsw) 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_GP endif u_p = 0.0_GP do j = 1, ns_stencil do i = 1, ns_stencil u_p = u_p + coeffs_interpol(i,j) * u_in(ind_stencil(i,j)) enddo enddo u_out(l) = u_p deallocate(ind_stencil) deallocate(coeffs_interpol) enddo end subroutine end module