submodule(perp_bnd_flux_m) perp_bnd_flux_build_s !! Routines related with building perpendicular boundary flux module implicit none integer, parameter :: nstencil = 4 integer, parameter, dimension(nstencil) :: stencil_i = [-1, 1, 0, 0] integer, parameter, dimension(nstencil) :: stencil_j = [0, 0, -1, 1] ! These parameters reflect the stencil of the elliptic operator, ! i.e. horizontal and vertical, but not diagonal neighbors are respected. contains module subroutine init_perp_bnd_flux(self, equi, mesh, penalisation, dbgout) class(perp_bnd_flux_t), intent(inout) :: self class(equilibrium_t), intent(in) :: equi type(mesh_cart_t), intent(in) :: mesh type(penalisation_t), intent(in) :: penalisation integer, intent(in), optional :: dbgout integer :: outi integer :: it, ki, l, sn, lngb, np_core, np_outer, np_target integer, allocatable, dimension(:) :: inds_core, inds_outer, inds_target real(GP), dimension(mesh%get_n_points()) :: aux_charfun ! Set output level outi = 0 if (present(dbgout)) then if (is_rank_info_writer) then outi = dbgout endif if (dbgout >= 3) then outi = dbgout ! every rank writes endif endif if (outi > 0) then write(get_stdout(),*)'' write(get_stdout(),*)'Building perpendicular flux utilities' endif ! Auxiliary array that ease to find penalisation functions ! w.r.t. global mesh index l aux_charfun = 0.0_GP !$omp parallel default(none) private(ki, l) & !$omp shared(mesh, penalisation, aux_charfun) !$omp do do ki = 1, mesh%get_n_points_inner() l = mesh%inner_indices(ki) aux_charfun(l) = penalisation%charfun(ki) enddo !$omp end do !$omp end parallel do it = 1, 2 np_core = 0 np_outer = 0 np_target = 0 ! Loop over inner point and check connection of stencil to boundary points do ki = 1, mesh%get_n_points_inner() l = mesh%inner_indices(ki) do sn = 1, nstencil ! Loop over stencil neighbours lngb = mesh%get_index_neighbor(stencil_i(sn), stencil_j(sn), l) if (is_core_boundary(ki, lngb)) then np_core = np_core + 1 if (it == 2) then inds_core(np_core) = lngb endif endif if (is_walldome_boundary(ki, lngb)) then np_outer = np_outer + 1 if (it == 2) then inds_outer(np_outer) = lngb endif endif if (is_target_boundary(ki, lngb)) then np_target = np_target + 1 np_outer = np_outer + 1 if (it == 2) then inds_target(np_target) = lngb inds_outer(np_outer) = lngb endif endif enddo enddo if (it == 1) then allocate(inds_core(np_core)) allocate(inds_outer(np_outer)) allocate(inds_target(np_target)) endif enddo ! Build polygons, core polygon with clockwise direction, others counter-clockwise call self%core%init(equi, mesh, inds_core, dir_counterclockwise=.false., dbgout=outi) deallocate(inds_core) call self%outer%init(equi, mesh, inds_outer, dbgout=outi) deallocate(inds_outer) call self%target%init(equi, mesh, inds_target, dbgout=outi) deallocate(inds_target) ! Build connectivity matrices call build_connectivity_polygon_matrix(mesh, penalisation, & self%core, self%conn_core) call build_connectivity_polygon_matrix(mesh, penalisation, & self%outer, self%conn_outer) call build_connectivity_polygon_matrix(mesh, penalisation, & self%target, self%conn_target) ! Build full_to_inner map allocate(self%full_to_inner(mesh%get_n_points())) !$omp parallel do default(none) private(l, ki), shared(self, mesh) do l = 1, mesh%get_n_points() ! Initialize with invalid value self%full_to_inner(l) = 0 if (mesh%is_inner_point(l)) then do ki = 1, mesh%get_n_points_inner() if (mesh%inner_indices(ki) == l) then self%full_to_inner(l) = ki exit endif enddo endif enddo !$omp end parallel do if (outi > 0) then write(get_stdout(),*)'done' write(get_stdout(),*)'' endif contains logical function is_core_boundary(ki, lcand) !! Local function that returns true if mesh point is part of core boundary integer, intent(in) :: ki !! Index of inner point that is connected to lcand integer, intent(in) :: lcand !! Mesh point to be checked integer :: district is_core_boundary = .false. district = mesh%get_district(lcand) if ( (mesh%is_boundary_point(lcand)) .and. & (district == DISTRICT_CORE) ) then is_core_boundary = .true. endif end function logical function is_walldome_boundary(ki, lcand) !! Local function that returns true if mesh point is part of wall or dome boundary !! excluding penalisation region integer, intent(in) :: ki !! Index of inner point that is connected to lcand integer, intent(in) :: lcand !! Mesh point to be checked integer :: district is_walldome_boundary = .false. district = mesh%get_district(lcand) if ( (penalisation%charfun(ki) < GP_EPS) .and. & (mesh%is_boundary_point(lcand)) .and. & ( (district == DISTRICT_WALL) .or. (district == DISTRICT_DOME) ) ) then is_walldome_boundary = .true. endif end function logical function is_target_boundary(ki, lcand) !! Local function that returns true if mesh point is part of target boundary integer, intent(in) :: ki !! Index of inner point that is connected to lcand integer, intent(in) :: lcand !! Mesh point to be checked is_target_boundary = .false. if ( (penalisation%charfun(ki) < GP_EPS) .and. & (mesh%is_inner_point(lcand)) .and. & (aux_charfun(lcand) > GP_EPS) ) then is_target_boundary = .true. endif end function end subroutine subroutine build_connectivity_polygon_matrix(mesh, penalisation, polygon, conn) !! Builds connectivity matrix for boundary poygon with inner mesh points type(mesh_cart_t), intent(in) :: mesh !! Mesh type(penalisation_t), intent(in) :: penalisation !! Penalisation type(polygon_in_mesh_t), intent(in) :: polygon !! Polygon type(csrmat_t), intent(out) :: conn !! Connectivity matrix integer :: it, ki, nz, is, ks, ks_glob, l ,sn, lcand real(GP), dimension(mesh%get_n_points()) :: aux_charfun ! Auxiliary array that eases to find penalisation function ! w.r.t. global mesh index l aux_charfun = 0.0_GP do ki = 1, mesh%get_n_points_inner() l = mesh%inner_indices(ki) aux_charfun(l) = penalisation%charfun(ki) enddo conn%ndim = polygon%get_np_total() conn%ncol = mesh%get_n_points() allocate(conn%i(conn%ndim+1)) do it = 1, 2 ! it = 1: Determine number of non-zeros and i-array of matrix conn ! it = 2: Create j and val-array of matrix conn nz = 0 do is = 1, polygon%get_nsegs() do ks = 1, polygon%get_nps(is) ks_glob = polygon%ki_seg(is) - 1 + ks if (it == 1) then conn%i(ks_glob) = nz + 1 endif l = polygon%get_ind_on_mesh(ks, is) do sn = 1, nstencil lcand = mesh%get_index_neighbor(stencil_i(sn), stencil_j(sn), l) if ( (mesh%is_inner_point(lcand)) .and. & (aux_charfun(lcand) < GP_EPS) ) then nz = nz + 1 if (it == 2) then conn%j(nz) = lcand conn%val(nz) = 1.0_GP endif endif enddo enddo enddo if (it == 1) then conn%i(conn%ndim+1) = nz + 1 conn%nnz = nz allocate(conn%j(conn%nnz), conn%val(conn%nnz)) endif enddo end subroutine module subroutine destructor(self) !! Destructor type(perp_bnd_flux_t), intent(inout) :: self !! Instance of type deallocate(self%full_to_inner) end subroutine end submodule