module helmholtz_operator_m !! Routines related with Helmholtz operator use precision_m, only : FP use screen_io_m, only : get_stderr use error_handling_m, only : handle_error, error_info_t use status_codes_m, only : PARALLAX_ERR_HELMHOLTZ use descriptors_m, only : BND_TYPE_DIRICHLET_ZERO, BND_TYPE_NEUMANN use csrmat_m, only : csrmat_t use mesh_cart_m, only : mesh_cart_t use boundaries_perp_m, only : compute_bndnmn_matrix_row use list_operations_m, only : rank_list implicit none public :: helmholtz_single_inner public :: helmholtz_single_boundary public :: helmholtz_single_ghost public :: build_helmholtz_csr contains pure function helmholtz_single_inner(mesh, u, ind, co, xiv, lambdav) result(res) !! Computes helmholtz operator for single inner grid point type(mesh_cart_t), intent(in) :: mesh !! Mesh real(FP), dimension(mesh%n_points), intent(in) :: u !! Variable integer, intent(in) :: ind !! Index where helmholtz operator should be computed !! MUST be index of an inner grid pointt (not checked for performance reasons) real(FP), dimension(mesh%n_points), intent(in), optional :: co !! Coefficient real(FP), intent(in), optional :: xiv !! Coefficient at mesh point real(FP), intent(in), optional :: lambdav !! Coefficient at mesh point real(FP) :: res !! result integer :: ind_lft, ind_rgt, ind_btm, ind_top real(FP) :: co_lft, co_rgt, co_btm, co_top real(FP) :: u_slf, u_lft, u_rgt, u_btm, u_top real(FP) :: two_hsqrd ind_lft = mesh%index_neighbor(-1, 0, ind) ind_rgt = mesh%index_neighbor(+1, 0, ind) ind_btm = mesh%index_neighbor( 0, -1, ind) ind_top = mesh%index_neighbor( 0, +1, ind) u_slf = u(ind) u_lft = u(ind_lft) u_rgt = u(ind_rgt) u_btm = u(ind_btm) u_top = u(ind_top) ! Two times co on staggered points if (present(co)) then co_lft = co(ind_lft) + co(ind) co_rgt = co(ind_rgt) + co(ind) co_btm = co(ind_btm) + co(ind) co_top = co(ind_top) + co(ind) else co_lft = 2.0_FP co_rgt = 2.0_FP co_btm = 2.0_FP co_top = 2.0_FP endif two_hsqrd = 2.0_FP*mesh%spacing_c**2 res = ( (co_btm + co_lft + co_rgt + co_top) * u_slf & - co_btm * u_btm & - co_lft * u_lft & - co_rgt * u_rgt & - co_top * u_top ) & / (two_hsqrd) if (present(xiv)) then res = res * xiv endif if (present(lambdav)) then res = res + lambdav*u_slf endif end function function helmholtz_single_boundary(mesh, u, ind, bnd_descr) result(res) !! Computes helmholtz operator for single boundary point type(mesh_cart_t), intent(in) :: mesh !! Mesh real(FP), dimension(mesh%n_points), intent(in) :: u !! Variable integer, intent(in) :: ind !! Index where helmholtz operator should be computed !! MUST be index of a boundary point (not checked for performance reasons) integer, intent(in) :: bnd_descr !! Descriptor for boundary condition real(FP) :: res !! result integer :: nz_nmn, k integer, dimension(3) :: cols_nmn real(FP), dimension(3) :: vals_nmn if (bnd_descr == BND_TYPE_DIRICHLET_ZERO) then res = u(ind) elseif (bnd_descr == BND_TYPE_NEUMANN) then call compute_bndnmn_matrix_row(mesh, ind, cols_nmn, vals_nmn, nz_nmn) res = 0.0_FP do k = 1, nz_nmn res = res + vals_nmn(k) * u(cols_nmn(k)) enddo else write(get_stderr(),*)& 'error(helmholtz_single_boundary): & &boundary descriptor not valid', ind, bnd_descr error stop endif end function pure function helmholtz_single_ghost(mesh, u, ind) result(res) !! Computes helmholtz operator for single boundary point type(mesh_cart_t), intent(in) :: mesh !! Mesh real(FP), dimension(mesh%n_points), intent(in) :: u !! Variable integer, intent(in) :: ind !! Index where helmholtz operator should be computed !! MUST be index of a ghost point (not checked for performance reasons) real(FP) :: res !! result res = u(ind) end function subroutine build_helmholtz_csr(mesh, hcsr, hdiag_inv, & co, xi, lambda, bnd_descrs) !! Builds Helmholtz matrix type(mesh_cart_t), intent(in) :: mesh !! Mesh type(csrmat_t), intent(inout) :: hcsr !! Helmholtz matrix. !! If hcsr%i is already allocated on input, the matrix is assumed !! to have the provided non-zero structure with all fields allocated, !! otherwise structure is allocated and created from scratch. real(FP), dimension(mesh%get_n_points()), & intent(out), optional :: hdiag_inv !! Inverse of diagonal of Helmholtz matrix real(FP), dimension(mesh%get_n_points()), & intent(in), optional :: co !! Coefficient within Helmholtz operator real(FP), dimension(mesh%get_n_points_inner()), & intent(in), optional :: xi !! Coefficient within Helmholtz operator real(FP), dimension(mesh%get_n_points_inner()), & intent(in), optional :: lambda !! Coefficient within Helmholtz operator integer, dimension(mesh%get_n_points_boundary()), & intent(in), optional :: bnd_descrs !! Descriptors_m for type of boundary conditions integer :: ki, kb, kg, row, nz, is, k, nz_nmn, cols_nmn(3) integer, dimension(mesh%get_n_points()) :: nnz_per_row integer, dimension(5) :: inds_stencil, rank_stencil real(FP), dimension(5) :: co_stencil real(FP) :: xiv, lambdav, two_hsqrd_inv, vals_nmn(3) if (.not. allocated(hcsr%i)) then ! Create non-zero structure, i.e. ndim, ncol, nnz and i-array, if (.not.present(bnd_descrs)) then call handle_error('Boundary descriptors must be present if matrix is created from scratch', & PARALLAX_ERR_HELMHOLTZ, __LINE__, __FILE__) endif ! Allocate j and vals array hcsr%ndim = mesh%get_n_points() hcsr%ncol = mesh%get_n_points() ! Determine number of non-zeros per row !$omp parallel default(none) & !$omp private(ki, row, cols_nmn, vals_nmn, nz_nmn) & !$omp shared(mesh, nnz_per_row, bnd_descrs) !$omp do do ki = 1, mesh%get_n_points_inner() row = mesh%inner_indices(ki) nnz_per_row(row) = 5 enddo !$omp end do !$omp do do kb = 1, mesh%get_n_points_boundary() row = mesh%boundary_indices(kb) call compute_bndnmn_matrix_row(mesh, row, & cols_nmn, vals_nmn, nz_nmn) nnz_per_row(row) = nz_nmn enddo !$omp end do !$omp do do kg = 1, mesh%get_n_points_ghost() row = mesh%ghost_indices(kg) nnz_per_row(row) = 1 enddo !$omp end do !$omp end parallel ! Create i-array of Helmholtz matrix allocate(hcsr%i(1:hcsr%ndim+1)) hcsr%i(1) = 1 do row = 2, hcsr%ndim+1 hcsr%i(row) = hcsr%i(row - 1) + nnz_per_row(row - 1) enddo hcsr%nnz = hcsr%i(hcsr%ndim + 1) - 1 allocate(hcsr%j(hcsr%nnz)) allocate(hcsr%val(hcsr%nnz)) endif two_hsqrd_inv = 1.0_FP / (2.0_FP * mesh%get_spacing_c()**2) !$omp parallel default(none) & !$omp private(ki, kb, kg, row, nz, is, xiv, lambdav, & !$omp inds_stencil, rank_stencil, co_stencil, & !$omp k, nz_nmn, cols_nmn, vals_nmn) & !$omp shared(mesh, hcsr, hdiag_inv, bnd_descrs, & !$omp co, xi, lambda, two_hsqrd_inv) !$omp do do ki = 1, mesh%get_n_points_inner() ! 5 point stencil on inner mesh row = mesh%inner_indices(ki) nz = hcsr%i(row) - 1 inds_stencil(1) = mesh%get_index_neighbor(-1, 0, row) ! left inds_stencil(2) = mesh%get_index_neighbor(+1, 0, row) ! right inds_stencil(3) = mesh%get_index_neighbor( 0, -1, row) ! bottom inds_stencil(4) = mesh%get_index_neighbor( 0, +1, row) ! top inds_stencil(5) = row ! Rank stencil, in order to build matrix in ascending column order call rank_list(inds_stencil, rank_stencil) if (present(xi)) then xiv = xi(ki) else xiv = 1.0_FP endif if (present(lambda)) then lambdav = lambda(ki) else lambdav = 1.0_FP endif if (present(co)) then co_stencil(1) = co(inds_stencil(1)) + co(inds_stencil(5)) co_stencil(2) = co(inds_stencil(2)) + co(inds_stencil(5)) co_stencil(3) = co(inds_stencil(3)) + co(inds_stencil(5)) co_stencil(4) = co(inds_stencil(4)) + co(inds_stencil(5)) else co_stencil(1:4) = 2.0_FP endif co_stencil(5) = -(co_stencil(1) + co_stencil(2) + & co_stencil(3) + co_stencil(4)) do is = 1, 5 nz = nz + 1 hcsr%j(nz) = inds_stencil(rank_stencil(is)) hcsr%val(nz) = -xiv * co_stencil(rank_stencil(is)) & * two_hsqrd_inv if (rank_stencil(is) == 5) then ! Adaptions for diagonal entry hcsr%val(nz) = lambdav + hcsr%val(nz) if (present(hdiag_inv)) then hdiag_inv(row) = 1.0_FP / hcsr%val(nz) endif endif enddo enddo !$omp end do if (present(bnd_descrs)) then !$omp do do kb = 1, mesh%get_n_points_boundary() ! boundary points row = mesh%boundary_indices(kb) nz = hcsr%i(row) - 1 call compute_bndnmn_matrix_row(mesh, row, & cols_nmn, vals_nmn, nz_nmn) call rank_list(cols_nmn(1:nz_nmn), rank_stencil(1:nz_nmn)) do k = 1, nz_nmn nz = nz + 1 hcsr%j(nz) = cols_nmn(rank_stencil(k)) if (bnd_descrs(kb) == BND_TYPE_DIRICHLET_ZERO) then if (cols_nmn(rank_stencil(k)) == row) then hcsr%val(nz) = 1.0_FP else hcsr%val(nz) = 0.0_FP endif elseif (bnd_descrs(kb) == BND_TYPE_NEUMANN) then hcsr%val(nz) = vals_nmn(rank_stencil(k)) else !$omp critical call handle_error('boundary descriptor not valid', & PARALLAX_ERR_HELMHOLTZ, __LINE__, __FILE__, & additional_info=& error_info_t('kb, bnd_descr=',[kb, bnd_descrs(kb)])) !$omp end critical endif if (present(hdiag_inv)) then if (cols_nmn(rank_stencil(k)) == row) then hdiag_inv(row) = 1.0_FP / hcsr%val(nz) endif endif enddo enddo !$omp end do !$omp do do kg = 1, mesh%get_n_points_ghost() ! identitity on ghost points row = mesh%ghost_indices(kg) nz = hcsr%i(row) - 1 nz = nz + 1 hcsr%j(nz) = row hcsr%val(nz) = 1.0_FP if (present(hdiag_inv)) then hdiag_inv(row) = 1.0_FP endif enddo !$omp end do endif !$omp end parallel end subroutine end module