submodule(mesh_cart_m) mesh_cart_build_s !! Submodule related with building of mesh use list_operations_m, only: unique_tuples, sort_tuples use constants_m, only: TWO_PI use screen_io_m, only: get_stdout implicit none contains module subroutine init_mesh_cart(self, equi, phi, lvl, spacing_f, & size_neighbor, size_ghost_layer, & extend_beyond_wall, dbgout) class(mesh_cart_t), intent(out) :: self class(equilibrium_t), intent(inout) :: equi real(FP), intent(in) :: phi integer, intent(in) :: lvl real(FP), intent(inout) :: spacing_f integer, intent(in) :: size_neighbor integer, intent(in) :: size_ghost_layer logical, intent(in):: extend_beyond_wall integer, intent(in), optional ::dbgout integer :: outi integer :: nxf, npwr ! 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 self%extend_beyond_wall = extend_beyond_wall if (outi >= 1) then write(get_stdout(), *) '' write(get_stdout(), *) & '------------------------------------------------' write(get_stdout(), 201) & phi, lvl, spacing_f, size_neighbor, size_ghost_layer, & extend_beyond_wall 201 FORMAT('Creating a full mesh' & 3X, 'phi = ', ES14.6E3, / & 3X, 'lvl = ', I3, / & 3X, 'spacing_f = ', ES14.6E3, / & 3X, 'size_neighbor = ', I3, / & 3X, 'size_ghost_layer = ', I3, / & 3X, 'extend_beyond_wall = ', L1 ) endif ! Set toroidal angle self%phi = modulo(phi, TWO_PI) ! Build firstly Cartesian grid select type(equi) type is (slab_t) if (outi >= 1) then write(get_stdout(), *) '...Creating mesh for slab geometry...' endif self%yperiodic = equi%is_yperiodic() ! Compute power exponent for number of mesh points in x-direction nxf = anint( (equi%xmax - equi%xmin) / spacing_f + 1) npwr = anint(log(1.0_FP * (nxf - 1)) / log(2.0_FP)) call self%build_cart_slab(equi, lvl, npwr) spacing_f = self%spacing_f if (outi >= 1) then write(get_stdout(), *) & 'Warning, spacing_f might have been adapted to', spacing_f endif class default if (outi >= 1) then write(get_stdout(), *) '...Creating mesh...' endif self%yperiodic = .false. call self%build_cart(equi, lvl, spacing_f) end select if (outi >= 1) then write(get_stdout(), *) '...building shaped cartesian arrays...' endif call self%build_shaped_cart_arrays() ! Build connectivity if (outi >= 1) then write(get_stdout(), *) '...building connectivity...' endif call self%build_connectivity(size_neighbor) ! Build boundary if (outi >= 1) then write(get_stdout(), *) '...building boundary...' endif call self%build_boundary() ! Build ghost layer if (outi >= 1) then write(get_stdout(), *) '...building ghost layer...' endif call self%build_ghost_layer(size_ghost_layer) if (outi >= 1) then write(get_stdout(), *) '...rebuilding shaped cartesian arrays...' endif call self%build_shaped_cart_arrays() call self%build_connectivity(size_neighbor) ! Build information for individual mesh points if (outi >= 1) then write(get_stdout(), *) '...building additional info...' endif call self%build_pinfo() call self%build_district(equi) ! Patch mesh if (outi >= 1) then write(get_stdout(), *) '...patching mesh...' endif call self%build_patch() ! Find red-black indices if (outi >= 1) then write(get_stdout(), *) '...identifying red-black indices...' endif call self%build_redblack() if (outi >= 1) then write(get_stdout(), *) 'Mesh succesfully built:' call self%display() write(get_stdout(), *) '' endif end subroutine module subroutine build_cart(self, equi, lvl, spacing_f) class(mesh_cart_t), intent(inout) :: self class(equilibrium_t), intent(inout) :: equi integer, intent(in) :: lvl real(FP), intent(in) :: spacing_f integer :: i, j, i_s, j_s, district, l real(FP) :: x, y, x_s, y_s integer :: length_x, length_y, i_step, j_step integer, dimension(:, :), allocatable :: work_cart_i, work_cart_j ! Set meta-information self%lvl = lvl self%spacing_f = spacing_f self%xmin = equi%xmin self%ymin = equi%ymin self%nx_f = anint((equi%xmax - equi%xmin) / spacing_f) self%ny_f = anint((equi%ymax - equi%ymin) / spacing_f) self%lvst = 2**(lvl - 1) self%spacing_c = self%spacing_f * self%lvst ! If not using a unit step size, the number of points in a unit-spaced ! array will be (end - start) / step + start ! N.b. the integer division is intentional length_x = (self%nx_f - 1) / self%lvst + 1 length_y = (self%ny_f - 1) / self%lvst + 1 ! Allocate arrays to store the i and j indices of all possible ! grid points. Initialise these to zero. allocate(work_cart_i(length_x, length_y), source=0) allocate(work_cart_j(length_x, length_y), source=0) ! Search for points !$omp parallel default(none) private(i_step, j_step, i, j, x, y, & !$omp& i_s, j_s, x_s, y_s, district) & !$omp& shared(length_x, length_y, self, work_cart_i, work_cart_j, equi) !$omp do search_loop: do j = 1, self%ny_f, self%lvst do i = 1, self%nx_f, self%lvst ! If not using a unit step size, the i^th point will ! correspond to the ((i - start) / step + start)^th point ! in a unit-spaced array i_step = ((i - 1) / self%lvst) + 1 j_step = ((j - 1) / self%lvst) + 1 x = self%xmin + (i - 1) * self%spacing_f y = self%ymin + (j - 1) * self%spacing_f ! Add point if interior domain is located within at least ! half a grid distance neighbor_loop: do j_s = -1, 1 do i_s = -1, 1 x_s = x + 0.5_FP * i_s * self%lvst * self%spacing_f y_s = y + 0.5_FP * j_s * self%lvst * self%spacing_f district = equi%district(x_s, y_s, self%get_phi()) if (self%point_inside(district)) then ! If the district corresponds to an interior point, ! record the i and j position of the point work_cart_i(i_step, j_step) = i work_cart_j(i_step, j_step) = j exit neighbor_loop endif enddo enddo neighbor_loop enddo enddo search_loop !$omp end do !$omp end parallel ! Find the number of non-zero points in the work arrays ! and allocate arrays to store the non-zero values self%n_points = count(work_cart_i /= 0) allocate(self%cart_i(self%n_points)) allocate(self%cart_j(self%n_points)) ! Fill the cart_i and cart_j arrays with the non-zero points self%cart_i = pack(work_cart_i, work_cart_i /= 0) self%cart_j = pack(work_cart_j, work_cart_j /= 0) end subroutine module subroutine build_connectivity(self, size_neighbor) class(mesh_cart_t), intent(inout) :: self integer, intent(in) :: size_neighbor integer :: n_points, l, ic, jc, i_s, j_s, ingb, jngb n_points = self%n_points self%size_neighbor = size_neighbor if (allocated(self%index_neighbor)) then deallocate(self%index_neighbor) endif allocate(self%index_neighbor( & -size_neighbor:size_neighbor, & -size_neighbor:size_neighbor, & n_points )) !$omp parallel default(none) private(l, ic, jc, i_s, j_s, ingb, jngb) & !$omp& shared(n_points, self, size_neighbor) !$omp do do l = 1, n_points ic = self%cart_i(l) jc = self%cart_j(l) do i_s = -size_neighbor, size_neighbor do j_s = -size_neighbor, size_neighbor ingb = ic + self%lvst * i_s jngb = jc + self%lvst * j_s self%index_neighbor(i_s, j_s, l) = & self%cart_to_mesh_index(ingb, jngb) enddo enddo enddo !$omp end do !$omp end parallel end subroutine module subroutine build_boundary(self) class(mesh_cart_t), intent(inout) :: self logical :: is_on_boundary integer :: l integer, dimension(self%get_n_points()) :: work_boundary_indices, & work_inner_indices ! Work arrays used to calculate boundary points in parallel. Note: ! these do not have the same shape as the corresponding ! self%boundary_indices and self%inner_indices arrays! !$omp parallel default(none) private(is_on_boundary, l)& !$omp& shared(self, work_inner_indices, work_boundary_indices) !$omp do do l = 1, self%get_n_points() is_on_boundary = .false. ! Check if a horizontal or vertical neighbor is not existent if (self%get_index_neighbor(0, -1, l) == 0) then is_on_boundary = .true. else if (self%get_index_neighbor(-1, 0, l) == 0) then is_on_boundary = .true. else if (self%get_index_neighbor( 1, 0, l) == 0) then is_on_boundary = .true. else if (self%get_index_neighbor( 0, 1, l) == 0) then is_on_boundary = .true. endif if (is_on_boundary) then work_inner_indices(l) = 0 work_boundary_indices(l) = l else work_inner_indices(l) = l work_boundary_indices(l) = 0 endif end do !$omp end do !$omp end parallel self%n_points_boundary = count(work_boundary_indices /= 0) self%n_points_inner = count(work_inner_indices /= 0) allocate(self%boundary_indices(self%n_points_boundary)) allocate(self%inner_indices(self%n_points_inner)) self%boundary_indices = pack(work_boundary_indices, & work_boundary_indices /= 0) self%inner_indices = pack(work_inner_indices, & work_inner_indices /= 0) end subroutine module subroutine build_ghost_layer(self, size_ghost_layer) class(mesh_cart_t), intent(inout) :: self integer, intent(in) :: size_ghost_layer integer :: ki, l, p, ic, jc, i_s, j_s, ig, jg, lg integer :: n_points_ghost, n_points_new integer, allocatable, dimension(:, :) :: work_ghost_cart integer, dimension(2) :: precedence integer, dimension(:), allocatable :: original_cart_i, original_cart_j integer, dimension(:,:,:), allocatable :: work_ghost_i, work_ghost_j allocate(original_cart_i(self%get_n_points())) allocate(original_cart_j(self%get_n_points())) allocate(work_ghost_i(1:self%n_points_inner, & & -size_ghost_layer:size_ghost_layer, & & -size_ghost_layer:size_ghost_layer)) allocate(work_ghost_j(1:self%n_points_inner, & & -size_ghost_layer:size_ghost_layer, & & -size_ghost_layer:size_ghost_layer)) ! Work arrays to allow parallel calculation of ghost points. ! These arrays have a shape such that they can hold ! the maximum number of possible non-zero elements (usually, they ! are mostly unfilled -- so this method has a large memory footprint) self%size_ghost_layer = size_ghost_layer !$omp parallel default(none) private(ki, l, ic, jc, i_s, j_s, & !$omp& ig, jg, lg) & !$omp& shared(self, work_ghost_i, work_ghost_j, size_ghost_layer) !$omp do do ki = 1, self%n_points_inner l = self%inner_indices(ki) ic = self%cart_i(l) jc = self%cart_j(l) do j_s = -size_ghost_layer, size_ghost_layer do i_s = -size_ghost_layer, size_ghost_layer ig = ic + i_s * self%lvst jg = jc + j_s * self%lvst ! Correction for slab geometry if periodic in y-direction if (self%yperiodic) then if (jg < 1) then jg = jg + self%ny_f elseif (jg > self%ny_f) then jg = jg - self%ny_f endif endif lg = self%cart_to_mesh_index(ig, jg) if (lg == 0) then work_ghost_i(ki, i_s, j_s) = ig work_ghost_j(ki, i_s, j_s) = jg else work_ghost_i(ki, i_s, j_s) = 0 work_ghost_j(ki, i_s, j_s) = 0 endif enddo enddo enddo !$omp end do !$omp end parallel n_points_ghost = count(work_ghost_i /= 0 .or. work_ghost_j /= 0) allocate(work_ghost_cart(2, n_points_ghost)) work_ghost_cart(1,:) = pack(work_ghost_i, & work_ghost_i /= 0 .or. work_ghost_j /= 0) work_ghost_cart(2,:) = pack(work_ghost_j, & work_ghost_i /= 0 .or. work_ghost_j /= 0) ! From here, work_ghost_i and work_ghost_j are not used. deallocate(work_ghost_i, work_ghost_j) ! Sort and unique ghost points call unique_tuples(2, n_points_ghost, work_ghost_cart) precedence = [2, 1] call sort_tuples(n_points_ghost, 2, precedence, work_ghost_cart) ! Add ghost points to mesh n_points_new = self%n_points + n_points_ghost original_cart_i = self%cart_i original_cart_j = self%cart_j deallocate(self%cart_i) deallocate(self%cart_j) allocate(self%cart_i(n_points_new)) allocate(self%cart_j(n_points_new)) allocate(self%ghost_indices(n_points_ghost)) !$omp parallel default(none) private(p,l) & !$omp& shared(self, work_ghost_cart, original_cart_i, original_cart_j, & !$omp& n_points_ghost) !$omp do do l = 1, self%n_points self%cart_i(l) = original_cart_i(l) self%cart_j(l) = original_cart_j(l) enddo !$omp end do nowait !$omp do do p = 1, n_points_ghost l = self%n_points + p self%cart_i(l) = work_ghost_cart(1, p) self%cart_j(l) = work_ghost_cart(2, p) self%ghost_indices(p) = l enddo !$omp end do !$omp end parallel self%n_points = n_points_new self%n_points_ghost = n_points_ghost deallocate(original_cart_i, original_cart_j) end subroutine module subroutine build_pinfo(self) class(mesh_cart_t), intent(inout) :: self integer :: ki, kb, kg, l, ch allocate(self%pinfo(self%n_points), source = 0) !$omp parallel default(none) private(ki, kb, kg, l) & !$omp& shared(self) !$omp do do ki = 1, self%n_points_inner l = self%inner_indices(ki) self%pinfo(l) = PINFO_INNER enddo !$omp end do nowait !$omp do do kb = 1, self%n_points_boundary l = self%boundary_indices(kb) self%pinfo(l) = PINFO_BOUNDARY enddo !$omp end do nowait !$omp do do kg = 1, self%n_points_ghost l = self%ghost_indices(kg) self%pinfo(l) = PINFO_GHOST enddo !$omp end do !$omp end parallel ! Check if all points have some info ch = minval(self%pinfo) if (ch <= 0) then call handle_error("Pinfo not set for some point", & PARALLAX_ERR_MESH, __LINE__, __FILE__) endif end subroutine module subroutine build_district(self, equi) class(mesh_cart_t), intent(inout) :: self class(equilibrium_t), intent(in) :: equi real(FP), parameter :: inv_sqrt2 = sqrt(2.0_FP) / 2.0_FP integer :: ki, kb, kg, l, district, nbx, nby real(FP) :: x, y, phi, nbn, x_n, y_n allocate(self%district(self%n_points), source = 0) phi = self%get_phi() !$omp parallel private(ki, kb, kg, l, x, y, district, & !$omp& nbx, nby, nbn, x_n, y_n) & !$omp& default(none) shared(self, equi, phi) !$omp do do ki = 1, self%n_points_inner l = self%inner_indices(ki) x = self%get_x(l) y = self%get_y(l) district = equi%district(x, y, phi) ! Correct district for inner points that are outside equi boundary ! by shifting them by a distance 1 / sqrt(2) * spacing_c towards ! the compute region (i.e. in the direction of the boundary normal) if (.not. self%point_inside(district)) then call self%normal_to_boundary(l, nbx, nby) nbn = sqrt(1.0_FP * (nbx**2 + nby**2)) x_n = x + (1.0_FP * nbx) / nbn * self%spacing_c * inv_sqrt2 y_n = y + (1.0_FP * nby) / nbn * self%spacing_c * inv_sqrt2 district = equi%district(x_n, y_n, phi) endif if (.not. self%point_inside(district)) then !$omp critical call handle_error('Inner point is not close enough to the & &inner region!', & PARALLAX_ERR_MESH, __LINE__, __FILE__, & additional_info=& error_info_t('Info on point (unstructured & &index, district, level, & &x, y, phi):', & [l, district, self%lvl], & [x, y, phi])) !$omp end critical endif self%district(l) = district enddo !$omp end do nowait !$omp do do kb = 1, self%n_points_boundary l = self%boundary_indices(kb) x = self%get_x(l) y = self%get_y(l) district = equi%district(x, y, phi) ! Correct district for boundary points that are inside equi boundary ! by shifting them by a distance 1 / sqrt(2) * spacing_c away from ! the compute region (i.e. opposite to the boundary normal) if (self%point_inside(district)) then call self%normal_to_boundary(l, nbx, nby) nbn = sqrt(1.0_FP * (nbx**2 + nby**2)) x_n = x - (1.0_FP * nbx) / nbn * self%spacing_c * inv_sqrt2 y_n = y - (1.0_FP * nby) / nbn * self%spacing_c * inv_sqrt2 district = equi%district(x_n, y_n, phi) endif if (self%point_inside(district)) then !$omp critical call handle_error("Boundary point is not close enough to the & &outer region!", & PARALLAX_ERR_MESH, __LINE__, __FILE__, & additional_info=& error_info_t("Info on point (unstructured & &index, district, level, & &x, y, phi):", & [l, district, self%lvl], & [x, y, phi])) !$omp end critical endif self%district(l) = district enddo !$omp end do nowait !$omp do do kg = 1, self%n_points_ghost l = self%ghost_indices(kg) x = self%get_x(l) y = self%get_y(l) district = equi%district(x, y, phi) if (self%point_inside(district)) then !$omp critical call handle_error("Ghost point is inside compute region!", & PARALLAX_ERR_MESH, __LINE__, __FILE__, & additional_info=& error_info_t("Info on point (unstructured & &index, district, level, & &x, y, phi):", & [l, district, self%lvl], & [x, y, phi])) !$omp end critical endif self%district(l) = district enddo !$omp end do !$omp end parallel end subroutine module subroutine build_patch(self) class(mesh_cart_t), intent(inout) :: self logical :: patch integer :: kb, l, i_s, j_s, ln, kp, lp, kg integer, parameter :: npatch_max = 100 ! Maximum number of points to be patched, ! very unlikely that this number is exceeded (if so review your mesh ! and possibly increase this number and recompile) integer :: npatch ! number of points to patch integer, dimension(npatch_max) :: patch_list ! List of points to be patched ! Find points to be patched, i.e. boundary points that have no ! connection to inner point npatch = 0 do kb = 1, self%n_points_boundary l = self%boundary_indices(kb) patch = .true. do j_s = -1, 1 do i_s = -1, 1 ln = self%index_neighbor(i_s, j_s, l) if (ln > 0) then if (self%is_inner_point(ln)) then patch = .false. endif endif enddo enddo if (patch) then npatch = npatch + 1 if (npatch > npatch_max) then call handle_error("Maximum number of patches reached", & PARALLAX_ERR_MESH, __LINE__, __FILE__, & additional_info=error_info_t( & "npatch_max = ", [npatch_max]) ) endif patch_list(npatch) = l endif enddo ! Move boundary to ghost points if (npatch > 0) then ! write(stdout, *)' warning: patch is applied for', npatch, & ! &'points with indices' ! write(stdout, *)' ',patch_list(1:npatch) do kp = 1, npatch lp = patch_list(kp) self%pinfo(lp) = PINFO_GHOST enddo deallocate(self%boundary_indices) deallocate(self%ghost_indices) self%n_points_boundary = self%n_points_boundary - npatch self%n_points_ghost = self%n_points_ghost + npatch allocate(self%boundary_indices(self%n_points_boundary)) allocate(self%ghost_indices(self%n_points_ghost)) kb = 0 kg = 0 do l = 1, self%n_points if (self%pinfo(l) == PINFO_BOUNDARY) then kb = kb + 1 self%boundary_indices(kb) = l endif if (self%pinfo(l) == PINFO_GHOST) then kg = kg + 1 self%ghost_indices(kg) = l endif enddo endif end subroutine module subroutine build_redblack(self) class(mesh_cart_t), intent(inout) :: self integer :: ki, l, ic, jc integer :: i_start, i_end integer, dimension(self%get_n_points_inner()) :: work_red_points, & work_black_points ! Work arrays to allow parallel computation of redblack_indices allocate(self%redblack_indices(self%get_n_points())) !$omp parallel default(none) shared(self, work_red_points, & !$omp& work_black_points) private(ki,l,ic,jc) !$omp do do ki = 1, self%get_n_points_inner() work_black_points(ki) = 0 work_red_points(ki) = 0 l = self%inner_indices(ki) ic = (self%get_cart_i(l) - 1) / self%get_lvst() + 1 jc = (self%get_cart_j(l) - 1) / self%get_lvst() + 1 if (mod(ic + jc, 2) == 0) then work_red_points(ki) = l else if (mod(ic + jc, 2) == 1) then work_black_points(ki) = l endif end do !$omp end do !$omp end parallel self%n_points_red = count(work_red_points /= 0) self%n_points_black = count(work_black_points /= 0) i_start = 0 i_end = self%get_n_points_red() self%redblack_indices(i_start+1:i_end) = pack(work_red_points, & work_red_points /= 0) i_start = i_end i_end = i_end + self%get_n_points_black() self%redblack_indices(i_start+1:i_end) = pack(work_black_points, & work_black_points /= 0) i_start = i_end i_end = i_end + self%get_n_points_boundary() self%redblack_indices(i_start + 1:i_end) = self%boundary_indices i_start = i_end i_end = i_end + self%get_n_points_ghost() self%redblack_indices(i_start + 1:i_end) = self%ghost_indices end subroutine module subroutine build_cart_slab(self, equi, lvl, npwr) class(mesh_cart_t), intent(inout) :: self class(equilibrium_t), intent(inout) :: equi integer, intent(in) :: lvl integer, intent(in) :: npwr integer :: npwrc, nx_c, ny_c, l, ic, jc ! Set meta-information self%lvl = lvl self%lvst = 2**(lvl - 1) self%xmin = equi%xmin self%ymin = equi%ymin self%nx_f = 2**npwr + 1 self%spacing_f = (equi%xmax - equi%xmin) / (self%nx_f - 1) self%spacing_c = self%spacing_f * self%lvst npwrc = npwr - (lvl - 1) nx_c = 2**npwrc + 1 if (self%yperiodic) then self%ny_f = 2**npwr ny_c = 2**npwrc else self%ny_f = 2**npwr + 1 ny_c = 2**npwrc + 1 endif self%n_points = nx_c * ny_c allocate(self%cart_i(self%n_points)) allocate(self%cart_j(self%n_points)) l = 0 do jc = 1, ny_c do ic = 1, nx_c l = l + 1 self%cart_i(l) = (ic - 1) * self%lvst + 1 self%cart_j(l) = (jc - 1) * self%lvst + 1 enddo enddo end subroutine logical module function point_inside(self, district) class(mesh_cart_t), intent(in) :: self integer, intent(in) :: district if ( (district == DISTRICT_CLOSED) & .or. (district == DISTRICT_SOL) & .or. (district == DISTRICT_PRIVFLUX) ) then point_inside = .true. elseif ( (district == DISTRICT_CORE) & .or. (district == DISTRICT_OUT) ) then point_inside = .false. elseif ( (district == DISTRICT_WALL) & .or. (district == DISTRICT_DOME)) then if (self%extend_beyond_wall) then point_inside = .true. else point_inside = .false. endif else call handle_error('district not valid', & PARALLAX_ERR_MESH, __LINE__, __FILE__, & additional_info=& error_info_t('district=',[district])) endif end function end submodule