submodule (map_factory_m) map_factory_s use MPI use screen_io_m, only: progress_bar, get_stderr, get_stdout use comm_handling_m, only: is_master use precision_m, only: MPI_FP use fieldline_tracer_m, only: trace use interpolation_m, only: interpol_coeffs use list_operations_m, only: sort_and_sum use mpi_mapping_auxiliaries_m, only: get_mpipattern, communicate_matrixmpi use gauss_quadrature_m, only: gauss_legendre use dommaschk_equilibrium_m, only: dommaschk_t implicit none contains module subroutine build_map_row(equi, & mesh_base, mesh_target, dphi, x_base, y_base, & intorder, xorder, nz, jrow, valrow, arclength, & use_fixed_stencil, use_gauss_quadrature) class(equilibrium_t), intent(inout) :: equi type(mesh_cart_t), intent(in) :: mesh_base type(mesh_cart_t), intent(in) :: mesh_target real(FP), intent(in) :: dphi real(FP), intent(in) :: x_base real(FP), intent(in) :: y_base integer, intent(in) :: intorder integer, intent(in) :: xorder integer, intent(inout) :: nz integer, intent(out), dimension(nz)::jrow real(FP), intent(out), dimension(nz)::valrow real(FP), intent(out) :: arclength logical, intent(in), optional :: use_fixed_stencil logical, intent(in), optional :: use_gauss_quadrature integer, parameter :: ns_max_nearest = 6 !! Size of area that is searched for nearest point real(FP) :: phi_base, phi_target real(FP) :: x_target, y_target integer :: nzwrk, nsubx, nsuby, p1, p2, i, j, k real(FP) :: dh, xsub_base, ysub_base, xsub_target, ysub_target real(FP) :: area_target, xsw, ysw, xrel, yrel integer :: ns_stencil, intorder_actual, lsw integer, allocatable, dimension(:,:) :: ind_interpol_stencil real(FP), allocatable, dimension(:,:) :: coeffs_interpol integer, dimension(2**(xorder*2)*(intorder+1)**2) :: jwrk real(FP),dimension(2**(xorder*2)*(intorder+1)**2) :: valwrk real(FP), allocatable, dimension(:) :: gauss_x, gauss_w logical :: fixst_actual logical :: gq_actual real(FP) :: PHI_TO_DO ! Check for consistent input if ( (mod(intorder,2) /= 1) .or. (intorder < 1) ) then write(get_stderr(), *) 'error(build_map_row), intorder must be & &odd and > 1, but is',intorder error stop endif phi_base = mesh_base%get_phi() phi_target = mesh_target%get_phi() ! Adjust optional parameters, if necessary fixst_actual = .false. if (present(use_fixed_stencil)) fixst_actual = use_fixed_stencil gq_actual = .false. if (present(use_gauss_quadrature)) gq_actual = use_gauss_quadrature ! Subpoints if (.not. gq_actual) then nsubx = 2**xorder nsuby = nsubx dh = mesh_base%get_spacing_f() / nsubx else nsubx = 2**xorder nsuby = nsubx allocate( gauss_x(nsubx) ) allocate( gauss_w(nsubx) ) call gauss_legendre(gauss_x, gauss_w) endif nzwrk = 0 jwrk = 0 valwrk = 0.0_FP ! Field line tracing ! TODO: ! For axisymmetric equilibria tracing shoud be independent ! of phi_init, but only depend on dphi. ! However the trace routine seems to be a little sensitive to ! phi_init. So passing for axisymmetric equilibria ! phi_init=phi_base instead of phi_init=0 gives different ! results according to trace precision, which spoils CI/CD. ! So as workaround phi_init=0 is passed to trace for ! axisymmetric equilibria. ! --> Trace should be checked for consistency in ! axisymmetric equilibria. if (equi%is_axisymmetric()) then PHI_TO_DO = 0.0_FP else PHI_TO_DO = phi_base endif if (fixst_actual) then ! Midpoint stencil is reused for every mapped quadrature node call trace(x_init = x_base, & y_init = y_base, & phi_init = PHI_TO_DO, & ! TODO !!! dphi = dphi, & equi = equi, & x_end = x_target, & y_end = y_target) call mesh_target%determine_interpolation_stencil( & xp=x_target, yp=y_target, intorder=intorder, & ns_max_nearest=ns_max_nearest, & ns_stencil=ns_stencil, & intorder_actual=intorder_actual, & ind_interpol_stencil=ind_interpol_stencil) allocate( coeffs_interpol(ns_stencil, ns_stencil) ) endif do p1 = 1, nsubx do p2 = 1, nsuby if (.not. gq_actual) then xsub_base = (x_base - mesh_base%get_spacing_f() / 2.0_FP & + dh / 2.0_FP) + dh * (p1 - 1) ysub_base = (y_base - mesh_base%get_spacing_f() / 2.0_FP & + dh / 2.0_FP) + dh * (p2 - 1) else xsub_base = x_base + gauss_x(p1) & * mesh_base%get_spacing_f() / 2.0_FP ysub_base = y_base + gauss_x(p2) & * mesh_base%get_spacing_f() / 2.0_FP endif call trace(x_init = xsub_base, & y_init = ysub_base, & phi_init = PHI_TO_DO, & ! TODO !!! dphi = dphi, & equi = equi, & x_end = xsub_target, & y_end = ysub_target, & arclen = arclength) if (.not. fixst_actual) then ! Stencil is recomputed for every mapped quadrature node call mesh_target%determine_interpolation_stencil( & xp=xsub_target, yp=ysub_target, intorder=intorder, & ns_max_nearest=ns_max_nearest, & ns_stencil=ns_stencil, & intorder_actual=intorder_actual, & ind_interpol_stencil=ind_interpol_stencil) allocate( coeffs_interpol(ns_stencil, ns_stencil) ) endif if (intorder_actual >= 1) then ! Polynomial interpolation of order intorder_actual ! Get south-west point of interpolation stencil lsw = ind_interpol_stencil(1, 1) xsw = mesh_target%get_x(lsw) ysw = mesh_target%get_y(lsw) ! Get interpolation coefficients xrel = (xsub_target - xsw) / mesh_target%get_spacing_f() yrel = (ysub_target - ysw) / mesh_target%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 ! Mapped area area_target = abs(equi%btor(xsub_base, ysub_base, phi_base)) & / abs(equi%btor(xsub_target, ysub_target, phi_target)) if (.not. gq_actual) then area_target = area_target / (nsubx * nsuby) else area_target = area_target & * gauss_w(p1) * gauss_w(p2) / 4.0_FP endif ! Fill matrix do j = 1, ns_stencil do i = 1, ns_stencil if (ind_interpol_stencil(i, j) /= 0) then nzwrk = nzwrk+1 jwrk(nzwrk) = ind_interpol_stencil(i, j) valwrk(nzwrk) = coeffs_interpol(i, j) & * equi%btor(xsub_target, ysub_target, & phi_target) * area_target else write(get_stderr(), *) & 'error(build_fieldlinemap_matrixrow): & &ind_interpol_stencil(i, j) = 0' error stop endif enddo enddo if (.not. fixst_actual) then ! Stencil is recomputed for every mapped quadrature node deallocate(coeffs_interpol) deallocate(ind_interpol_stencil) endif enddo enddo if (allocated(coeffs_interpol)) deallocate(coeffs_interpol) if (allocated(ind_interpol_stencil)) deallocate(ind_interpol_stencil) if (allocated(gauss_x)) deallocate(gauss_x) if (allocated(gauss_w)) deallocate(gauss_w) ! Sort and sum entries call sort_and_sum(nzwrk, jwrk, valwrk) if (nzwrk > nz) then write(get_stderr(), *) & 'error(build_fieldlinemap_matrixrow): nz > allocated' error stop endif nz = nzwrk do k = 1, nzwrk jrow(k) = jwrk(k) valrow(k) = valwrk(k) enddo end subroutine module subroutine create_map_matrix(map_matrix, comm, equi, & mesh_base, mesh_target, phi_direction, & intorder, xorder, flux_box_surface_integral, & use_fixed_stencil, use_gauss_quadrature, only_first_field_period, & n_turns, dbgout, arclength_array, only_inner_points) type(csrmat_t), intent(out) :: map_matrix integer, intent(in) :: comm class(equilibrium_t), intent(inout) :: equi type(mesh_cart_t), intent(in) :: mesh_base type(mesh_cart_t), intent(in) :: mesh_target integer, intent(in) :: phi_direction integer, intent(in) :: intorder integer, intent(in), optional :: xorder logical, intent(in), optional :: flux_box_surface_integral logical, intent(in), optional :: use_fixed_stencil logical, intent(in), optional :: use_gauss_quadrature logical, intent(in), optional :: only_first_field_period integer, intent(in), optional :: n_turns integer, intent(in), optional :: dbgout real(FP), intent(out), dimension(mesh_base%get_n_points()), & optional :: arclength_array logical, intent(in), optional :: only_inner_points integer, parameter :: nz_alfac = 16*4 !! Maximum number of non-zero entries per row of map matrices integer, dimension(mesh_base%get_n_points()) :: nz_wrk integer, dimension(nz_alfac, mesh_base%get_n_points()) :: j_wrk real(FP), dimension(nz_alfac, mesh_base%get_n_points()) :: val_wrk real(FP) :: phi_base, phi_target, dphi logical :: fbsi_actual integer :: np, np_loc, lloc_start, lloc_end, nnz_local, nnz_total, l, ipnt, k real(FP) :: x_base, y_base integer :: outi, xorder_actual logical :: fixst_actual logical :: gq_actual real(FP), dimension(mesh_base%get_n_points()) :: arclength_local logical :: only_first_field_period_actual integer :: n_turns_actual logical :: only_inner_actual real(FP) :: dphi_max integer :: ierr ! Adjust optional parameters, if necessary xorder_actual = 0 if (present(xorder)) xorder_actual = xorder fbsi_actual = .false. if (present(flux_box_surface_integral)) fbsi_actual & = flux_box_surface_integral fixst_actual = .false. if (present(use_fixed_stencil)) fixst_actual = use_fixed_stencil gq_actual = .false. if (present(use_gauss_quadrature)) gq_actual = use_gauss_quadrature n_turns_actual = 0 if (present(n_turns)) n_turns_actual = n_turns only_first_field_period_actual = .false. if (present(only_first_field_period)) only_first_field_period_actual = & only_first_field_period only_inner_actual = .false. if (present(only_inner_points)) only_inner_actual = only_inner_points ! Check for consistency of input if (mesh_base%get_lvl() /= 1) then write(get_stderr(), *) 'error(create_map_matrix): mesh level & &must be 1, but is', mesh_base%get_lvl() error stop endif if ((xorder_actual > 0) .and. (.not. fbsi_actual) ) then ! for xorder > 0 only integral formulation allowed write(get_stderr(), *) & 'error(create_map_matrix):: For xorder > 0, & &only integral formulation is allowed' error stop endif if (present(arclength_array) .and. xorder_actual /= 0) then ! Field line lengths are only returned if xorder is 0, i.e. if the ! trace starts from the mesh points themselves and not from subcells write(get_stderr(), *) & 'error(create_map_matrix):: To return field line lengths, & &xorder must equal 0!' error stop endif ! Calculate dphi phi_base = mesh_base%get_phi() phi_target = mesh_target%get_phi() if(only_first_field_period_actual) then select type(equi) type is(dommaschk_t) dphi_max = TWO_PI / equi%get_num_field_periods() class default write(get_stderr(), *) & 'error(create_map_matrix): The only_first_field_period & &option can only be specified for Dommaschk equilibria!' error stop end select else dphi_max = TWO_PI endif dphi = calc_dphi(phi_base, phi_target, phi_direction, & dphi_max=dphi_max, n_turns=n_turns_actual) ! 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 map matrix' write(get_stdout(), 100) & phi_base, dphi, intorder, xorder_actual, fbsi_actual, & fixst_actual, gq_actual, present(arclength_array) 100 FORMAT( & X,'phi_base =',ES14.6E3, / & X,'dphi =',ES14.6E3, / & X,'intorder =',I4, / & X,'xorder =',I4, / & X,'integral_method =',L1 / & X,'use_fixed_stencil =',L1 / & X,'use_gauss_quadrature =',L1 / & X,'arclength_array =',L1 & ) endif if (equi%is_axisymmetric()) then call get_mpipattern(comm, 1, mesh_base%get_n_points(), & np, np_loc, lloc_start, lloc_end) else call get_mpipattern(MPI_COMM_SELF, 1, mesh_base%get_n_points(), & np, np_loc, lloc_start, lloc_end) endif ! Initialize field line length buffer to zero for MPI communication arclength_local = 0.0_FP ! Build rows of matrix !$omp parallel default(none) & !$omp shared(lloc_start, lloc_end, nz_wrk, j_wrk, val_wrk, & !$omp arclength_local, equi, mesh_base, mesh_target, dphi, & !$omp intorder, xorder_actual, fixst_actual, gq_actual, & !$omp np_loc, outi, only_inner_actual) & !$omp private(l, x_base, y_base) !$omp do do l = lloc_start, lloc_end if (is_master() .and. outi >= 1) then call progress_bar(l-lloc_start+1, np_loc, 20) endif if (only_inner_actual .and. .not. mesh_base%is_inner_point(l)) then nz_wrk(l) = 0 else x_base = mesh_base%get_x(l) y_base = mesh_base%get_y(l) ! Build row of matrix nz_wrk(l) = nz_alfac call build_map_row(equi, mesh_base, mesh_target, dphi, & x_base, y_base, intorder, xorder_actual, & nz_wrk(l), j_wrk(:, l), val_wrk(:, l), & arclength_local(l), fixst_actual, gq_actual) endif enddo !$omp end do !$omp end parallel ! Get total number of nonzero elements of map matrix nnz_local = sum(nz_wrk(lloc_start:lloc_end)) if (equi%is_axisymmetric()) then call MPI_allreduce(nnz_local, nnz_total, 1, MPI_INTEGER, & MPI_sum, comm, ierr) else nnz_total = nnz_local endif ! Allocate map_matrix (CSR) map_matrix%ndim = np map_matrix%ncol = mesh_target%get_n_points() map_matrix%nnz = nnz_local allocate(map_matrix%i(mesh_base%get_n_points()+1)) allocate(map_matrix%j(nnz_total)) allocate(map_matrix%val(nnz_total)) ! Establish matrix in CSR format !$omp parallel default(none) & !$omp shared(lloc_start, lloc_end, nz_wrk, j_wrk, val_wrk, & !$omp equi, mesh_base, fbsi_actual, phi_base, map_matrix) & !$omp private(l, x_base, y_base, ipnt, k) !$omp do do l = lloc_start, lloc_end x_base = mesh_base%get_x(l) y_base = mesh_base%get_y(l) map_matrix%i(l) = sum(nz_wrk(lloc_start:l-1)) + 1 do k = 1, nz_wrk(l) ipnt = map_matrix%i(l) map_matrix%j(ipnt+k-1) = j_wrk(k, l) if (.not. fbsi_actual) then map_matrix%val(ipnt+k-1) = val_wrk(k, l) & / equi%btor(x_base, y_base, phi_base) else map_matrix%val(ipnt+k-1) = val_wrk(k, l) endif enddo enddo !$omp end do !$omp end parallel map_matrix%i(mesh_base%get_n_points()+1) = map_matrix%nnz + 1 ! Gather matrix and optional fieldline lengths if (equi%is_axisymmetric()) then call communicate_matrixmpi(comm, map_matrix%ndim, & np_loc, nnz_total, map_matrix%i, & map_matrix%j, map_matrix%val) map_matrix%nnz = map_matrix%i(mesh_base%get_n_points()+1) - 1 if (present(arclength_array)) then call MPI_allreduce(arclength_local, arclength_array, & mesh_base%get_n_points(), MPI_FP, & MPI_sum, comm, ierr) endif else if (present(arclength_array)) arclength_array = arclength_local endif if (is_master() .and. outi >= 1) then write(get_stdout(), *) & 'communication of map matrix done, nz = ', map_matrix%nnz write(get_stdout(), *) '...................................' endif end subroutine end submodule