submodule(mesh_cart_m) mesh_cart_s !! Ordinary mesh routines use list_operations_m, only: findindex_2d use screen_io_m, only: get_stderr, get_stdout implicit none contains pure real(FP) module function get_phi(self) class(mesh_cart_t), intent(in) :: self get_phi = self%phi end function module subroutine set_phi(self, new_phi) class(mesh_cart_t), intent(inout) :: self real(FP), intent(in) :: new_phi self%phi = new_phi end subroutine pure integer module function get_lvl(self) class(mesh_cart_t), intent(in) :: self get_lvl = self%lvl end function pure integer module function get_lvst(self) class(mesh_cart_t), intent(in) :: self get_lvst = self%lvst end function pure integer module function get_n_points(self) class(mesh_cart_t), intent(in) :: self get_n_points = self%n_points end function pure real(FP) module function get_spacing_f(self) class(mesh_cart_t), intent(in) :: self get_spacing_f = self%spacing_f end function pure real(FP) module function get_spacing_c(self) class(mesh_cart_t), intent(in) :: self get_spacing_c = self%spacing_c end function pure integer module function get_size_neighbor(self) class(mesh_cart_t), intent(in) :: self get_size_neighbor = self%size_neighbor end function pure integer module function get_cart_i(self, ind) class(mesh_cart_t), intent(in) :: self integer, intent(in) :: ind get_cart_i = self%cart_i(ind) end function pure integer module function get_cart_j(self, ind) class(mesh_cart_t), intent(in) :: self integer, intent(in) :: ind get_cart_j = self%cart_j(ind) end function pure real(FP) module function get_x(self, ind) class(mesh_cart_t), intent(in) :: self integer, intent(in) :: ind get_x = self%xmin + (self%cart_i(ind)-1) * self%spacing_f end function pure real(FP) module function get_y(self, ind) class(mesh_cart_t), intent(in) :: self integer, intent(in) :: ind get_y = self%ymin + (self%cart_j(ind)-1) * self%spacing_f end function pure integer module function get_n_points_inner(self) class(mesh_cart_t), intent(in) :: self get_n_points_inner = self%n_points_inner end function pure integer module function get_n_points_boundary(self) class(mesh_cart_t), intent(in) :: self get_n_points_boundary = self%n_points_boundary end function pure integer module function get_n_points_ghost(self) class(mesh_cart_t), intent(in) :: self get_n_points_ghost = self%n_points_ghost end function pure real(FP) module function get_xmin(self) class(mesh_cart_t), intent(in) :: self get_xmin = self%xmin end function pure real(FP) module function get_ymin(self) class(mesh_cart_t), intent(in) :: self get_ymin = self%ymin end function pure integer module function get_n_points_red(self) class(mesh_cart_t), intent(in) :: self get_n_points_red = self%n_points_red end function pure integer module function get_n_points_black(self) class(mesh_cart_t), intent(in) :: self get_n_points_black = self%n_points_black end function pure integer module function get_index_neighbor(self, ishift, jshift, ind) class(mesh_cart_t), intent(in) :: self integer, intent(in) :: ishift integer, intent(in) :: jshift integer, intent(in) :: ind get_index_neighbor = self%index_neighbor(ishift, jshift, ind) end function integer module function cart_to_mesh_index(self, i, j) result(res) class(mesh_cart_t), intent(in) :: self integer, intent(in) :: i integer, intent(in) :: j integer :: j_s j_s = j ! Correction for slab geometry if periodic in y-direction if (self%yperiodic) then if (j_s < 1) then j_s = j_s + self%ny_f elseif (j_s > self%ny_f) then j_s = j_s - self%ny_f endif endif if (.not.(allocated(self%cart_shaped_l))) then ! N.b. findindex_2d has a saved attribute and therefore cannot be ! pure. This routine is O(N^2), will be very slow for large grids call findindex_2d(1, self%n_points, self%cart_i, self%cart_j, & i, j_s, res) return endif if ((i < self%min_i) .or. (i > self%max_i) .or. & (j_s < self%min_j) .or. (j_s > self%max_j)) then res = 0 else ! If available, use the cart_shaped_l stored array -- extremely fast res = self%cart_shaped_l(i, j_s) endif end function module subroutine build_shaped_cart_arrays(self) class(mesh_cart_t), intent(inout) :: self integer :: i, j, l self%min_i = minval(self%cart_i) self%max_i = maxval(self%cart_i) self%min_j = minval(self%cart_j) self%max_j = maxval(self%cart_j) if (allocated(self%cart_shaped_l)) then deallocate(self%cart_shaped_l) endif allocate(self%cart_shaped_l(self%min_i:self%max_i, & self%min_j:self%max_j), source=0) !$omp parallel default(none) private(i, j, l) shared(self) !$omp do do l = 1, self%n_points i = self%cart_i(l) j = self%cart_j(l) self%cart_shaped_l(i, j) = l enddo !$omp enddo !$omp end parallel end subroutine module subroutine deallocate_shaped_cart_arrays(self) class(mesh_cart_t), intent(inout) :: self if (allocated(self%cart_shaped_l)) then deallocate(self%cart_shaped_l) endif end subroutine pure logical module function is_inner_point(self, ind) class(mesh_cart_t), intent(in) :: self integer, intent(in) :: ind if (self%pinfo(ind) == PINFO_INNER) then is_inner_point = .true. else is_inner_point = .false. endif end function pure logical module function is_boundary_point(self, ind) class(mesh_cart_t), intent(in) :: self integer, intent(in) :: ind if (self%pinfo(ind) == PINFO_BOUNDARY) then is_boundary_point = .true. else is_boundary_point = .false. endif end function pure logical module function is_ghost_point(self, ind) class(mesh_cart_t), intent(in) :: self integer, intent(in) :: ind if (self%pinfo(ind) == PINFO_GHOST) then is_ghost_point = .true. else is_ghost_point = .false. endif end function pure integer module function get_district(self, ind) class(mesh_cart_t), intent(in) :: self integer, intent(in) :: ind get_district = self%district(ind) end function module subroutine determine_interpolation_stencil(self, xp, yp, intorder, & ns_max_nearest, ns_stencil, intorder_actual, ind_interpol_stencil) class(mesh_cart_t), intent(in) :: self real(FP), intent(in) :: xp real(FP), intent(in) :: yp integer, intent(in) :: intorder integer, intent(in) :: ns_max_nearest integer, intent(out) :: ns_stencil integer, intent(out) :: intorder_actual integer, allocatable, dimension(:,:), intent(out) :: ind_interpol_stencil integer :: ind_nearest ! Find size of largest possible interpolation stencil ! -> ns_stencil call self%get_surrounding_indices(xp=xp, yp=yp, ns=intorder+1, & ns_complete=ns_stencil ) ! Determine final interpolation order and interpolation stencil if (ns_stencil > 0) then intorder_actual = ns_stencil - 1 allocate( ind_interpol_stencil(ns_stencil, ns_stencil) ) call self%get_surrounding_indices(xp=xp, yp=yp, & ns=ns_stencil, ind_sur=ind_interpol_stencil) elseif (ns_stencil == 0) then ! Find nearest point in (ns_max_nearest x ns_max_nearest) surr. grid call self%get_surrounding_indices(xp=xp, yp=yp, & ns=ns_max_nearest, ind_nearest=ind_nearest) if (ind_nearest /= 0) then ! Nearest neighbour intorder_actual = 0 ns_stencil = 1 allocate( ind_interpol_stencil(ns_stencil, ns_stencil) ) ind_interpol_stencil(1, 1) = ind_nearest else ! No interpolation at all intorder_actual = -1 ns_stencil = 0 allocate( ind_interpol_stencil(ns_stencil, ns_stencil) ) endif endif end subroutine module subroutine get_surrounding_indices(self, xp, yp, ns, ind_sur, ns_complete, ind_nearest) class(mesh_cart_t), intent(in) :: self real(FP), intent(in) :: xp real(FP), intent(in) :: yp integer, intent(in) :: ns integer, intent(out), optional, dimension(ns, ns) :: ind_sur integer, intent(out), optional :: ns_complete integer, intent(out), optional :: ind_nearest integer :: isw, jsw, in, jn, ia, ja, ln, k real(FP) :: xn, yn, dst, dst_probe integer, dimension(ns, ns) :: ind_sur_wrk ! TODO: routine has some limitations, check for consistent input -------------------------- if (self%lvl /= 1) then write(get_stderr(), *) & 'error(mesh%get_surrounding_indices) only working for lvl=1', & self%lvl error stop endif if ( (mod(ns,2) /= 0) .or. (ns < 2) ) then write(get_stderr(), *) & 'error(mesh%get_surrounding_indices) & &ns must be even and larger equal 2 but is',ns error stop endif ! Point soutwest of grid point isw = floor( (xp-self%xmin) / self%spacing_f ) + 1 jsw = floor( (yp-self%ymin) / self%spacing_f ) + 1 ! Southwest point of surrounding grid isw = isw - ns/2 + 1 jsw = jsw - ns/2 + 1 do jn = 1, ns do in = 1, ns ia = isw + in-1 ja = jsw + jn-1 ind_sur_wrk(in, jn) = self%cart_to_mesh_index(ia, ja) enddo enddo if (present(ns_complete)) then ! check up to which dimension surrounding grid is complete ns_complete = 0 kloop: do k = 1, ns/2 do jn = ns/2-k+1, ns/2+k do in = ns/2-k+1, ns/2+k ln = ind_sur_wrk(in, jn) if (ln == 0) then exit kloop endif enddo enddo ns_complete = ns_complete + 2 enddo kloop endif if (present(ind_nearest)) then ind_nearest = 0 dst = FP_LARGEST do jn = 1, ns do in = 1, ns ln = ind_sur_wrk(in, jn) if (ln /= 0) then xn = self%get_x(ln) yn = self%get_y(ln) dst_probe = sqrt((xp-xn)**2 + (yp-yn)**2) if (dst_probe < dst) then ind_nearest = ln dst = dst_probe endif endif enddo enddo endif if (present(ind_sur)) then ind_sur = ind_sur_wrk endif end subroutine get_surrounding_indices module subroutine normal_to_boundary(self, ind, nbx, nby) class(mesh_cart_t), intent(in) :: self integer, intent(in) :: ind integer, intent(out) :: nbx integer, intent(out) :: nby integer :: in, jn, lngb integer, dimension(-1:1, -1:1) :: disc_stencil if (.not. self%is_boundary_point(ind)) then ! This should only occur in the case of using an inner boundary with ! a 3D equilibrium write(get_stdout(), *) "Calculating boundary normal for & &non-boundary point, mesh index = ", ind endif ! Compute normal gradient to boundary via finite difference on suitable ! discrete level set function do jn = -1, 1 do in = -1, 1 lngb = self%get_index_neighbor(in, jn, ind) if (lngb == 0) then disc_stencil(in, jn) = 0 elseif (self%is_inner_point(lngb)) then disc_stencil(in, jn) = 2 elseif (self%is_boundary_point(lngb)) then disc_stencil(in, jn) = 1 elseif (self%is_ghost_point(lngb)) then disc_stencil(in, jn) = 0 endif enddo enddo nbx = disc_stencil(1, 0) - disc_stencil(-1, 0) nby = disc_stencil(0, 1) - disc_stencil( 0, -1) if(nbx == 0 .and. nby == 0) then call handle_error('Boundary normal could not be determined', & PARALLAX_ERR_MESH, __LINE__, __FILE__, & additional_info=& error_info_t('Info on point (unstructured & &index, level, x, y, phi):', & [ind, self%lvl], & [self%get_x(ind), & self%get_y(ind), & self%get_phi()])) endif end subroutine module subroutine display_mesh_cart(self) class(mesh_cart_t), intent(in) :: self if (.not. is_master()) then return endif write(get_stdout(),*) & 'Information on mesh -----------------------------------' write(get_stdout(), 101) & self%phi, & self%extend_beyond_wall, & self%lvl, & self%lvst, & self%spacing_f, & self%spacing_c, & self%n_points, & self%xmin, & self%ymin, & self%nx_f, & self%ny_f, & self%size_neighbor, & self%size_ghost_layer, & self%n_points_inner, & self%n_points_boundary, & self%n_points_ghost, & self%n_points_red, & self%n_points_black 101 FORMAT( 3X,'phi = ',ES14.6E3 /, & 3X,'extend_beyond_wall = ',L1 /, & 3X,'lvl = ',I4 /, & 3X,'lvst = ',I4 /, & 3X,'spacing_f = ',ES14.6E3 /, & 3X,'spacing_c = ',ES14.6E3 /, & 3X,'n_points = ',I9 /, & 3X,'xmin = ',ES14.6E3 /, & 3X,'ymin = ',ES14.6E3 /, & 3X,'nx_f = ',I9 /, & 3X,'nx_f = ',I9 /, & 3X,'size_neighbor = ',I4 /, & 3X,'size_ghost_layer = ',I4 /, & 3X,'n_points_inner = ',I9 /, & 3X,'n_points_boundary = ',I9 /, & 3X,'n_points_ghost = ',I9 /, & 3X,'n_points_red = ',I9 /, & 3X,'n_points_black = ',I9 ) write(get_stdout(), *) & '-------------------------------------------------------' end subroutine module subroutine expose_data(self, data_object) class(mesh_cart_t), intent(in), target :: self type(mesh_cart_data_t), intent(out) :: data_object ! Copy values data_object%extend_beyond_wall= self%extend_beyond_wall data_object%lvl = self%lvl data_object%lvst = self%lvst data_object%spacing_f = self%spacing_f data_object%spacing_c = self%spacing_c data_object%n_points = self%n_points data_object%phi = self%phi data_object%xmin = self%xmin data_object%ymin = self%ymin data_object%nx_f = self%nx_f data_object%ny_f = self%ny_f data_object%size_neighbor = self%size_neighbor data_object%size_ghost_layer = self%size_ghost_layer data_object%n_points_inner = self%n_points_inner data_object%n_points_boundary = self%n_points_boundary data_object%n_points_ghost = self%n_points_ghost data_object%n_points_red = self%n_points_red data_object%n_points_black = self%n_points_black ! Now assign pointers data_object%cart_i = c_loc(self%cart_i(1)) data_object%cart_j = c_loc(self%cart_j(1)) data_object%index_neighbor = c_loc( & self%index_neighbor(-self%size_neighbor, -self%size_neighbor, 1)) data_object%inner_indices = c_loc( & self%inner_indices(1)) data_object%boundary_indices = c_loc( & self%boundary_indices(1)) if (self%n_points_ghost >= 1) then data_object%ghost_indices = c_loc( & self%ghost_indices(1)) else data_object%ghost_indices = c_null_ptr end if data_object%pinfo = c_loc( & self%pinfo(1)) data_object%district = c_loc( & self%district(1)) data_object%redblack_indices = c_loc( & self%redblack_indices(1)) end subroutine end submodule