submodule(helmholtz_solver_m) helmholtz_solver_mgmres_cpu_s !! Helmholtz solver based on GMRES with multigrid preconditioner !! Accelerates and stabilises pure multigrid solver !! Employs MKL RCI_ISS solver based on reverse communication interface use screen_io_m, only : get_stderr, get_stdout #ifdef ENABLE_MKL use mkl_rci #endif implicit none contains module subroutine create_mgmres_cpu(self, multigrid, & bnd_type_core, bnd_type_wall, & bnd_type_dome, bnd_type_out, & co, lambda, xi) class(helmholtz_solver_mgmres_cpu_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 self%ndim = multigrid%get_np(1) self%np_inner = multigrid%get_np_inner(1) ! Create multigrid solver allocate(self%multigrid_solver) call self%multigrid_solver%create(multigrid, bnd_type_core, & bnd_type_wall, bnd_type_dome, bnd_type_out, & co, lambda, xi) end subroutine module subroutine init_mgmres_cpu(self, rtol, restol_zero, & nrestart, niter_max, & mgr_ncycle, & mgr_smoother, mgr_npresmooth, & mgr_npostsmooth, & mgr_dirsolver, & dbgout) class(helmholtz_solver_mgmres_cpu_t), intent(inout) :: self real(FP), intent(in) :: rtol real(FP), intent(in) :: restol_zero integer, intent(in) :: nrestart integer, intent(in) :: niter_max integer, intent(in) :: mgr_ncycle class(splitting_t), intent(inout), allocatable :: mgr_smoother integer, intent(in) :: mgr_npresmooth integer, intent(in) :: mgr_npostsmooth class(direct_solver_t), allocatable, intent(inout) :: mgr_dirsolver integer, intent(in) :: dbgout ! Set parameters self%rtol = rtol self%restol_zero = restol_zero self%nrestart = nrestart self%niter_max = niter_max self%dbgout = dbgout ! Initialise multigrid solver call self%multigrid_solver%init(mgr_ncycle, & mgr_smoother, mgr_npresmooth, mgr_npostsmooth, & mgr_dirsolver) end subroutine module subroutine update_mgmres_cpu(self, co, lambda, xi, & bnd_type_core, bnd_type_wall, & bnd_type_dome, bnd_type_out) class(helmholtz_solver_mgmres_cpu_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 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 self%multigrid_solver%update(co, lambda, xi, & bnd_type_core, bnd_type_wall, & bnd_type_dome, bnd_type_out) else call self%multigrid_solver%update(co, lambda, xi) endif end subroutine module subroutine solve_mgmres_cpu(self, rhs, sol, res, info) class(helmholtz_solver_mgmres_cpu_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 type(csrmat_t), pointer :: hcsr !! Matrix on finest level integer :: rci_request integer, dimension(128) :: ipar real(FP), dimension(128) :: dpar real(FP), dimension((2*self%nrestart+1)*self%ndim& +self%nrestart*(self%nrestart+9)/2 + 1) :: tmp integer :: i real(FP), dimension(self%ndim) :: wrk ! Get matrix on finest level and dimension of problem hcsr => self%multigrid_solver%get_hcsr_finest_pointer() ! Set parameters for RCI_ISS solver #if defined ENABLE_MKL && defined DOUBLE_PREC call dfgmres_init(self%ndim, sol, rhs, rci_request, ipar, dpar, tmp) #elif !defined ENABLE_MKL write(get_stderr(), *) & 'error(helmholtz%solve_mgmres): illegal call in MKL-free build' error stop #else write(get_stderr(), *) & 'error(helmholtz%solve_mgmres): not available for single precision' error stop #endif ipar(1) = self%ndim ! Specifies the size of the problem ipar(2) = 6 ! Output unit for warning messages ipar(5) = self%niter_max ! Maximum number of iterations ipar(8) = 1 ! dfgmres checks number of maximum iterations ipar(9) = 0 ! dfgmres does not compute stop criterion (residual) ipar(10) = 1 ! User computes stop criterion (residual) -> rci_request = 2 ipar(11) = 1 ! Use preconditioner -> rci_request = 3 ipar(12) = 0 ! User performs the test for zero norm -> rci_request = 4 ipar(15) = self%nrestart ! Number of iterations before restart dpar(1) = self%rtol ! Specifies the relative tolerance ! Check parameters #if defined ENABLE_MKL && defined DOUBLE_PREC ipar(7) = 0 ! Deactivate warning message for check as warning is print out due to manual setting of nrestart call dfgmres_check(self%ndim, sol, rhs, rci_request, ipar, dpar, tmp) ipar(7) = 1 ! Activate warning messages in the following again #endif if (rci_request /= 0) then if (rci_request == -1001) then ! This code is only a warning message caused by setting ! nrestart manually. Code may continue normally ! See https://community.intel.com/t5/Intel-oneAPI-Math-Kernel-Library/RCI-ISS-solver-FGMRES-working-well-in-IPS-XE-2020-Update-4-does/m-p/1271247/highlight/true else write(get_stderr(),*) & 'error(helmholtz%solve_mgmres): & ¶meter check not succesfull', rci_request error stop endif endif !$OMP PARALLEL !$OMP DO do i = 1, self%ndim wrk(i) = rhs(i) enddo !$OMP END DO !$OMP END PARALLEL gmresloop: do #if defined ENABLE_MKL && defined DOUBLE_PREC call dfgmres(self%ndim, sol, wrk, rci_request, ipar, dpar, tmp) #endif if (rci_request == 0) then ! The routine completes task normally, the solution is found. ! This occurs only if the stopping tests are fully automatic. exit gmresloop elseif (rci_request == 1) then ! Multiply the matrix by tmp(ipar(22)) and put the result ! in tmp(ipar(23)) !$OMP PARALLEL call csr_times_vec(hcsr, tmp(ipar(22):ipar(22)-1+self%ndim), & tmp(ipar(23):ipar(23)-1+self%ndim)) !$OMP END PARALLEL elseif (rci_request == 2) then ! Get solution into wrk and check residual ! (perform the stopping tests, if fail continue). ipar(13) = 1 #if defined ENABLE_MKL && defined DOUBLE_PREC call dfgmres_get(self%ndim, sol, wrk, rci_request, ipar, & dpar, tmp, info) #endif call csr_residuum(hcsr, wrk, rhs, res, epsb=self%restol_zero) if (self%dbgout > 0) then write(get_stdout(), *) & '#iteration =', info, 'Residual = ', res endif if (res < dpar(1)) then exit gmresloop endif elseif (rci_request == 3) then ! Precondition tmp(ipar(23) = MGR(tmp(ipar(22))) !$OMP PARALLEL !$OMP DO do i = 1, self%ndim tmp(ipar(23)+i-1) = 0.0_FP enddo !$OMP END DO !$OMP END PARALLEL call self%multigrid_solver%cycle(tmp(ipar(23):ipar(23)-1+self%ndim), & tmp(ipar(22):ipar(22)-1+self%ndim), 1) elseif (rci_request == 4) then ! Check the norm of the next orthogonal vector, ! which is contained in dpar(7). ! dpar(8) - contains the tolerance for the zero norm of ! the currently generated vector. The default value is 1.0D-12. if (dpar(7) <= dpar(8)) then ! Compute residuum ipar(13) = 1 #if defined ENABLE_MKL && defined DOUBLE_PREC call dfgmres_get(self%ndim, sol, wrk, rci_request, & ipar, dpar, tmp, info) #endif call csr_residuum(hcsr, wrk, rhs, res, & epsb=self%restol_zero) exit gmresloop endif else write(get_stderr(),*) & 'error(helmholtz%solve_mgmres): gmres failed', rci_request error stop endif enddo gmresloop ipar(13) = 0 ! Get solution into u #if defined ENABLE_MKL && defined DOUBLE_PREC call dfgmres_get(self%ndim, sol, wrk, rci_request, ipar, dpar, & tmp, info) #endif if (res > dpar(1)) then write(get_stdout(), *) 'GMRES not converged returning with info=-1' info = -1 endif end subroutine module subroutine destructor_mgmres_cpu(self) type(helmholtz_solver_mgmres_cpu_t), intent(inout) :: self if (allocated(self%multigrid_solver)) deallocate(self%multigrid_solver) self%ndim = 0 self%rtol = 0.0_FP self%restol_zero = 0.0_FP self%nrestart = 0 self%niter_max = 0 self%dbgout = 0 end subroutine end submodule