submodule(helmholtz_solver_m) helmholtz_solver_direct_s !! Helmholtz solver based on direct solver use screen_io_m, only : get_stderr use helmholtz_operator_m, only : build_helmholtz_csr use boundaries_perp_m, only : bnd_types_to_bnd_descrs implicit none contains module subroutine create_direct(self, multigrid, & bnd_type_core, bnd_type_wall, & bnd_type_dome, bnd_type_out, & co, lambda, xi) class(helmholtz_solver_direct_t), intent(inout) :: self type(multigrid_t), intent(inout) :: multigrid integer, intent(in) :: bnd_type_core integer, intent(in) :: bnd_type_wall integer, intent(in) :: bnd_type_dome integer, intent(in) :: bnd_type_out real(FP), dimension(multigrid%get_np(1)), & intent(in) :: co real(FP), dimension(multigrid%get_np_inner(1)), & intent(in) :: lambda real(FP), dimension(multigrid%get_np_inner(1)), & intent(in) :: xi integer, dimension(multigrid%get_np_boundary(1)) :: bnd_descrs self%mesh => multigrid%get_mesh_pointer(lvl=1) self%ndim = self%mesh%get_n_points() self%np_inner = self%mesh%get_n_points_inner() call bnd_types_to_bnd_descrs(self%mesh, bnd_type_core, bnd_type_wall, & bnd_type_dome, bnd_type_out, bnd_descrs) call build_helmholtz_csr(self%mesh, self%hcsr, & co=co, xi=xi, lambda=lambda, & bnd_descrs=bnd_descrs) end subroutine module subroutine init_direct(self, dirsolver) class(helmholtz_solver_direct_t), intent(inout) :: self class(direct_solver_t), allocatable, intent(inout) :: dirsolver ! Create solver and factorise if (.not.allocated(dirsolver)) then write(get_stderr(), *) & 'error(helmholtz_solver%init_direct): & &direct solver provided not allocated' error stop endif call move_alloc(dirsolver, self%dirslv) call self%dirslv%create(self%hcsr) call self%dirslv%factorise(self%hcsr,'ANALYSE_AND_NUMFAC') end subroutine module subroutine update_direct(self, co, lambda, xi, & bnd_type_core, bnd_type_wall, & bnd_type_dome, bnd_type_out) class(helmholtz_solver_direct_t), intent(inout) :: self real(FP), dimension(self%ndim), intent(in), target :: co real(FP), dimension(self%np_inner), intent(in), target :: lambda real(FP), dimension(self%np_inner), intent(in), target :: xi integer, intent(in), optional :: bnd_type_core integer, intent(in), optional :: bnd_type_wall integer, intent(in), optional :: bnd_type_dome integer, intent(in), optional :: bnd_type_out logical :: update_with_bnds integer, dimension(self%mesh%get_n_points_boundary()) :: bnd_descrs update_with_bnds = present(bnd_type_core) & .and. present(bnd_type_wall) & .and. present(bnd_type_dome) & .and. present(bnd_type_out) if (update_with_bnds) then call bnd_types_to_bnd_descrs(self%mesh, & bnd_type_core, bnd_type_wall, bnd_type_dome, bnd_type_out, & bnd_descrs) call build_helmholtz_csr(self%mesh, self%hcsr, & co=co, xi=xi, lambda=lambda, & bnd_descrs=bnd_descrs) else call build_helmholtz_csr(self%mesh, self%hcsr, & co=co, xi=xi, lambda=lambda) endif call self%dirslv%factorise(self%hcsr,'NUMFAC') end subroutine module subroutine solve_direct(self, rhs, sol, res, info) class(helmholtz_solver_direct_t), intent(inout) :: self real(FP), dimension(self%ndim), intent(inout), target :: rhs real(FP), dimension(self%ndim), intent(inout), target :: sol real(FP), intent(out) :: res integer, intent(out) :: info ! Solve call self%dirslv%backsub(self%hcsr, rhs, sol) ! Compute residuum call csr_residuum(self%hcsr, sol, rhs, res) if (res < 1E-8_FP) then info = 0 else info = -1 endif end subroutine module subroutine destructor_direct(self) type(helmholtz_solver_direct_t), intent(inout) :: self if (allocated(self%dirslv)) deallocate(self%dirslv) self%mesh => null() end subroutine end submodule