module map_factory_m !! Creates the parallel mapping matrices. These matrices are needed for the !! construction of parallel operators. use mesh_cart_m, only: mesh_cart_t use equilibrium_m, only: equilibrium_t use precision_m, only: FP use constants_m, only: TWO_PI use csrmat_m, only: csrmat_t use gauss_quadrature_m, only: gauss_legendre implicit none public :: create_map_matrix private :: build_map_row public :: calc_dphi interface 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) !! Creates map matrix for the parallel operators !! mapping from mesh_base to mesh_target. !! The resulting map_matrix can be used to compute !! umap = u(x_target, y_target, phi_target), where !! [x,y,phi]_target is map point in mesh_target. type(csrmat_t), intent(out) :: map_matrix !! Created map matrix integer, intent(in) :: comm !! MPI communicator used for parallel creation of the map matrix class(equilibrium_t), intent(inout) :: equi !! Equilibrium type(mesh_cart_t), intent(in) :: mesh_base !! Base mesh where mapping starts type(mesh_cart_t), intent(in) :: mesh_target !! Target mesh where map matrix ends integer, intent(in) :: phi_direction !! Option whose sign determines whether the mapping should procede !! from base to target in the positive or negative phi direction integer, intent(in) :: intorder !! Interpolation order of the map matrix integer, intent(in), optional :: xorder !! Subcell splitting for flux-box integration 2^(xorder) !! Default xorder = 0 !! Advanced topic, see: GRILLIX documentation, section 7.1.4 logical, intent(in), optional :: flux_box_surface_integral !! Indicates if !! ... point interpolation (default = false) !! ... or flux box surface integral is used (true) !! Advanced topic, see: GRILLIX documentation, section 7.1.4 logical, intent(in), optional :: use_fixed_stencil !! Indicates interpolation method for the mapped quadrature points: !! ... if false (default), the interpolation stencil is chosen !! individually for each mapped quadrature point !! ... if true, the interpolation stencil is fixed as the mapped !! midpoint's stencil logical, intent(in), optional :: use_gauss_quadrature !! Indicates quadrature formula for mapped integral !! ... if false (default), uniform mapped midpoint quadrature !! ... if true, mapped tensor Gauss-Legendre quadrature logical, intent(in), optional :: only_first_field_period !! Create a map matrix for simulations in only the first field !! period of a 3D equilibrium integer, optional, intent(in) :: n_turns !! Add additional turns to the difference in toroidal angle between !! mesh_base and mesh_target, default 0. Note that the sign is !! determined by phi_direction, and only the absolute value of !! n_turns is used. ! NOTE: This option is useful for simulations where the number of ! planes is comparable to the stencil size in phi. For ! example, if n_planes = 2, the dphi to the second neighbor ! with n_turns = 1 is correctly given as dphi = 2pi, rather ! than the default of dphi = 0. integer, intent(in), optional :: dbgout !! Specifies the number of information printed on screen. !! 0 represents no output, 1 minimal output, 2 debug output and 3 !! debug output by every MPI process real(FP), intent(out), dimension(mesh_base%get_n_points()), & optional :: arclength_array !! Length of the field lines traced from the points of mesh_base !! to mesh_target. Only valid if xorder = 0, i.e. if the trace !! starts from the mesh points themselves and not from subcells. logical, intent(in), optional :: only_inner_points !! Compute the mapping only on inner mesh_base points. The map !! matrix (and the arclength array, if given) for boundary and ghost !! points will return zero. end subroutine 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) !! Builds single row of field line map matrix !! !! Surface integral over quantity u can be computed !! with established information, i.e.: !! !! int_{A_target} Btor u dA (1) !! !! A_target is mapped area of area A_base, !! where A_base is area of size spacing^2 !! around point (x_base, y_base). !! !! This expression (1) can then be computed according to !! (~sparse matrix): !! !! (1) = sum_{k=1,nz} u(jrow(k)) * valrow(k). !! !! Interpolation order is reduced (down to nearest point) !! if interpolation stamp is not complete. class(equilibrium_t), intent(inout) :: equi !! Equilibrium type(mesh_cart_t), intent(in) :: mesh_base !! Base mesh where mapping starts starts type(mesh_cart_t), intent(in) :: mesh_target !! Target mesh where map matrix ends real(FP), intent(in) :: dphi !! Toroidal span between base and target mesh real(FP), intent(in) :: x_base !! x-coordinate of point to be mapped real(FP), intent(in) :: y_base !! y-coordinate of point to be mapped integer, intent(in) :: intorder !! Desired interpolation order, must be odd number > 0 !! Actual interpolation order might be reduced depending !! on availability of interpolation stencil integer, intent(in) :: xorder !! Subcell splitting for flux-box integration 2^(xorder) integer, intent(inout) :: nz !! Number of points involved in interpolation !! (on input size of workspace) integer, intent(out), dimension(nz)::jrow !! Indices of points involved in interpolation real(FP), intent(out), dimension(nz)::valrow !! Map values ~ interpolation coefficients real(FP), intent(out) :: arclength !! Field line length from the mesh point (x_base, y_base, phi_base) !! traced through the angle dphi. This value is returned for all !! parameter combinations, but is only meaningful if xorder = 0. logical, intent(in), optional :: use_fixed_stencil !! Indicates interpolation method for the mapped quadrature points: !! ... if false (default), the interpolation stencil is chosen !! individually for each mapped quadrature point !! ... if true, the interpolation stencil is fixed as the mapped !! midpoint's stencil logical, intent(in), optional :: use_gauss_quadrature !! Indicates quadrature formula for mapped integral !! ... if false (default), uniform mapped midpoint quadrature !! ... if true, mapped tensor Gauss-Legendre quadrature end subroutine end interface contains pure real(FP) function calc_dphi(phi_base, phi_target, phi_direction, & dphi_max, n_turns) result(dphi) !! Calculates the difference in toroidal angle dphi between phi_base !! and phi_target, in the direction phi_direction. The maximum dphi is !! set by dphi_max, and extra turns can be added with n_turns. ! NOTE: This is a separate function to more easily unit test it real(FP), intent(in) :: phi_base !! Initial toroidal angle real(FP), intent(in) :: phi_target !! Final toroidal angle integer, intent(in) :: phi_direction !! The sign of this argument determines whether the difference in !! phi should be taken from base to target in the positive or !! negative phi direction real(FP), optional, intent(in) :: dphi_max !! Maximum difference in toroidal angle, default 2pi. This option can !! also be thought of as setting where the periodic boundary is in phi. integer, optional, intent(in) :: n_turns !! Add additional full turns +/- dphi_max to dphi, default zero. Note !! that the sign is determined by phi_direction, and only the absolute !! value of n_turns is used. real(FP) :: dphi_max_local if(present(dphi_max)) then dphi_max_local = dphi_max else dphi_max_local = TWO_PI endif ! Calculate dphi in the range [0, +/- dphi_max), with the sign matching ! that of phi_direction dphi = modulo(phi_target - phi_base, & sign(1, phi_direction) * dphi_max_local) ! If requested, add turns dphi_max if(present(n_turns)) then dphi = dphi + dphi_max_local * sign(1, phi_direction) * abs(n_turns) endif end function end module