module boundaries_perp_m !! Setting of perpendicular boundary conditions use precision_m, only : FP, FP_NAN, FP_EPS use screen_io_m, only : get_stderr use error_handling_m, only : handle_error, error_info_t use status_codes_m, only : PARALLAX_ERR_BOUNDARIES use descriptors_m, only : BND_TYPE_NONE, BND_TYPE_DIRICHLET_ZERO, & BND_TYPE_NEUMANN, DISTRICT_CORE, DISTRICT_WALL, & DISTRICT_DOME, DISTRICT_OUT use mesh_cart_m, only : mesh_cart_t implicit none public :: set_perpbnds public :: extrapolate_ghost_points public :: compute_bndnmn_matrix_row public :: bnd_types_to_bnd_descrs private :: find_points_upstream_normal contains subroutine set_perpbnds(mesh, bnd_descrs, u, bnd_vals) !! Sets non-homogeneous perpendicular boundary condition for some quantity u type(mesh_cart_t), intent(in) :: mesh !! Mesh integer, intent(in), dimension(mesh%get_n_points_boundary()) :: bnd_descrs !! Descriptors for boundary condition at individual boundary points real(FP), dimension(mesh%get_n_points()), intent(inout) :: u !! Variable for which boundary conditions are set real(FP), dimension(mesh%get_n_points_boundary()), intent(in), optional :: bnd_vals !! Values (or normal gradient) for boundary condition, default homogeneous=0 integer :: kb, p real(FP) :: bval, u_upstream, dist_upstream integer :: nbx, nby, ind_b1, ind_b2 !$OMP DO do kb = 1, mesh%get_n_points_boundary() p = mesh%boundary_indices(kb) if (present(bnd_vals)) then bval = bnd_vals(kb) else bval = 0.0_FP endif if (bnd_descrs(kb) == BND_TYPE_DIRICHLET_ZERO) then u(p) = bval elseif (bnd_descrs(kb) == BND_TYPE_NEUMANN) then ! Get discrete normal direction to boundary call mesh%normal_to_boundary(p, nbx, nby) ! Find indices of points involved in setting (Neumann) boundary conditions call find_points_upstream_normal(mesh, p, nbx, nby, ind_b1, ind_b2, dist_upstream) if ( (ind_b1 /= 0) .and. (ind_b2 == 0) ) then u_upstream = u(ind_b1) elseif ( (ind_b1 /= 0) .and. (ind_b2 /= 0) ) then u_upstream = ( u(ind_b1) + u(ind_b2) ) / 2.0_FP else write(get_stderr(),*)'error(set_perpbnds): upstream points could not be found',kb, nbx, nby error stop endif dist_upstream = dist_upstream * mesh%get_spacing_c() u(p) = u_upstream - dist_upstream * bval elseif (bnd_descrs(kb) == BND_TYPE_NONE) then ! Do not change value else write(get_stderr(),*)'error(set_perpbnds): bnd_descrs not valid', kb, bnd_descrs(kb) error stop endif enddo !$OMP END DO end subroutine subroutine compute_bndnmn_matrix_row(mesh, ind, cols, vals, nz) !! Computes single row of Neumann boundary matrix type(mesh_cart_t), intent(in) :: mesh !! Mesh integer, intent(in) :: ind !! Mesh index (=row of boundary matrix) !! MUST be a boundary index integer, intent(out), dimension(3) :: cols !! Column indices, where non-zero values may be present real(FP), intent(out), dimension(3) :: vals !! Values of non-zero elements integer, intent(out) :: nz !! Number of non-zero elements integer :: nbx, nby, ind_b1, ind_b2 real(FP) :: dist_upstream ! Get discrete normal direction to boundary call mesh%normal_to_boundary(ind, nbx, nby) ! Find indices of points involved in setting (Neumann) boundary conditions call find_points_upstream_normal(mesh, ind, nbx, nby, ind_b1, ind_b2, dist_upstream) dist_upstream = dist_upstream * mesh%get_spacing_c() cols(1) = ind vals(1) = -1.0_FP / dist_upstream if ( (ind_b1 /= 0) .and. (ind_b2 == 0) ) then nz = 2 cols(2) = ind_b1 vals(2) = 1.0_FP / dist_upstream elseif ( (ind_b1 /= 0) .and. (ind_b2 /= 0) ) then nz = 3 cols(2) = ind_b1 vals(2) = 0.5_FP / dist_upstream cols(3) = ind_b2 vals(3) = 0.5_FP / dist_upstream else write(get_stderr(),*)'error(compute_bndnmn_matrix_row): upstream points could not be found',ind, nbx, nby error stop endif end subroutine subroutine bnd_types_to_bnd_descrs(mesh, bnd_type_core, bnd_type_wall, & bnd_type_dome, bnd_type_out, bnd_descrs) ! Converts boundary types to descriptor array type(mesh_cart_t), intent(in) :: mesh !! Mesh integer, intent(in) :: bnd_type_core !! Boundary type for core boundary integer, intent(in) :: bnd_type_wall !! Boundary type for wall boundary integer, intent(in) :: bnd_type_dome !! Boundary type for dome boundary integer, intent(in) :: bnd_type_out !! Boundary type for outer(mask) boundary integer, dimension(mesh%get_n_points_boundary()), & intent(out) :: bnd_descrs !! Boundary descriptors on boundary points integer :: kb, l, district !$omp parallel default(none) private(kb, l, district) & !$omp shared(mesh, bnd_descrs, & !$omp bnd_type_core, bnd_type_wall, bnd_type_dome, bnd_type_out) !$omp do do kb = 1, mesh%get_n_points_boundary() l = mesh%boundary_indices(kb) district = mesh%get_district(l) if (district == DISTRICT_CORE) then bnd_descrs(kb) = bnd_type_core elseif (district == DISTRICT_WALL) then bnd_descrs(kb) = bnd_type_wall elseif (district == DISTRICT_DOME) then bnd_descrs(kb) = bnd_type_dome elseif (district == DISTRICT_OUT) then bnd_descrs(kb) = bnd_type_out else !$omp critical call handle_error('district not valid', & PARALLAX_ERR_BOUNDARIES, __LINE__, __FILE__, & additional_info=& error_info_t('kb, district=',[kb, district])) !$omp end critical endif enddo !$omp end do !$omp end parallel end subroutine subroutine extrapolate_ghost_points(mesh, u) !! Extrapolates ghost points via inverse distance extrapolation type(mesh_cart_t), intent(in) :: mesh !! Mesh real(FP), dimension(mesh%get_n_points()), intent(inout) :: u !! Variable for which ghost points are set integer :: sng, kg, l, in, jn, ln real(FP) :: wgh, nrm sng = mesh%get_size_neighbor() !$OMP DO do kg = 1, mesh%get_n_points_ghost() l = mesh%ghost_indices(kg) u(l) = 0.0_FP nrm = 0.0_FP do jn = -sng, sng do in = -sng, sng ln = mesh%get_index_neighbor(in,jn,l) if ((ln /= 0).and.(ln /= l)) then if (.not.mesh%is_ghost_point(ln)) then wgh = 1.0_FP/sqrt(1.0_FP*(in**2+jn**2)) nrm = nrm+wgh u(l) = u(l) + wgh*u(ln) endif endif enddo enddo if (nrm > FP_EPS) then u(l) = u(l)/nrm endif enddo !$OMP END DO end subroutine pure subroutine find_points_upstream_normal(mesh, p, nbx, nby, ind_b1, ind_b2, dist_discrete) !! Finds mesh indices of points that are upstream to normal of boundary point !! Required to set Neumann boundary conditions type(mesh_cart_t), intent(in) :: mesh !! Mesh integer, intent(in) :: p !! Mesh index of boundary point integer, intent(in) :: nbx !! x-component of discrete normal vector integer, intent(in) :: nby !! y-component of discrete normal vector integer, intent(out) :: ind_b1 !! First upstream point integer, intent(out) :: ind_b2 !! Second upstream point (=0 if not present) real(FP), intent(out) :: dist_discrete !! Distance to upstream point in units of mesh distance if ( (abs(nbx) == 2).and.(abs(nby)==0) ) then ind_b1 = mesh%get_index_neighbor(nbx/2, 0, p) ind_b2 = 0 dist_discrete = 1.0_FP elseif ( (abs(nbx) == 0).and.(abs(nby)==2) ) then ind_b1 = mesh%get_index_neighbor(0, nby/2, p) ind_b2 = 0 dist_discrete = 1.0_FP elseif ( (abs(nbx) == 2).and.(abs(nby)==2) ) then ind_b1 = mesh%get_index_neighbor(nbx/2, nby/2, p) ind_b2 = 0 dist_discrete = sqrt(2.0_FP) elseif ( (abs(nbx) == 2).and.(abs(nby)==1) ) then ind_b1 = mesh%get_index_neighbor(nbx/2, 0, p) ind_b2 = mesh%get_index_neighbor(nbx/2, nby, p) dist_discrete = sqrt(1.0_FP**2+0.5_FP**2) elseif ( (abs(nbx) == 1).and.(abs(nby)==2) ) then ind_b1 = mesh%get_index_neighbor(0, nby/2, p) ind_b2 = mesh%get_index_neighbor(nbx, nby/2, p) dist_discrete = sqrt(1.0_FP**2+0.5_FP**2) elseif ( (abs(nbx) == 1).and.(abs(nby)==1) ) then ind_b1 = mesh%get_index_neighbor(nbx, nby, p) ind_b2 = 0 dist_discrete = sqrt(2.0_FP) else ind_b1 = 0 ind_b2 = 0 dist_discrete = FP_NAN endif end subroutine end module