gmres_sx.f90

Go to the documentation of this file.

! Copyright (c) 2021-2024, The Neko Authors
! All rights reserved.
!
! Redistribution and use in source and binary forms, with or without
! modification, are permitted provided that the following conditions
! are met:
!
!   * Redistributions of source code must retain the above copyright
!     notice, this list of conditions and the following disclaimer.
!
!   * Redistributions in binary form must reproduce the above
!     copyright notice, this list of conditions and the following
!     disclaimer in the documentation and/or other materials provided
!     with the distribution.
!
!   * Neither the name of the authors nor the names of its
!     contributors may be used to endorse or promote products derived
!     from this software without specific prior written permission.
!
! THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
! "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
! LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS
! FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE
! COPYRIGHT OWNER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT,
! INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING,
! BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES;
! LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER
! CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT
! LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN
! ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
! POSSIBILITY OF SUCH DAMAGE.
!
module gmres_sx
  use krylov, only : ksp_t, ksp_monitor_t
  use precon, only : pc_t
  use ax_product, only : ax_t
  use num_types, only : rp
  use field, only : field_t
  use coefs, only : coef_t
  use gather_scatter, only : gs_t, gs_op_add
  use scalar_bc_projector, only : scalar_bc_projector_t
  use vector_bc_projector, only : vector_bc_projector_t, &
       vector_bc_projector_components
  use math, only : glsc3, rzero, rone, copy, cmult2, col2, col3, add2s2, abscmp
  use comm, only : neko_comm, mpi_real_precision
  use mpi_f08
  implicit none
  private
 
  type, public, extends(ksp_t) :: sx_gmres_t
     integer :: lgmres = 30
     real(kind=rp), allocatable :: w(:)
     real(kind=rp), allocatable :: c(:)
     real(kind=rp), allocatable :: r(:)
     real(kind=rp), allocatable :: z(:,:)
     real(kind=rp), allocatable :: h(:,:)
     real(kind=rp), allocatable :: ml(:)
     real(kind=rp), allocatable :: v(:,:)
     real(kind=rp), allocatable :: s(:)
     real(kind=rp), allocatable :: mu(:)
     real(kind=rp), allocatable :: gam(:)
     real(kind=rp), allocatable :: wk1(:)
     real(kind=rp) :: rnorm
   contains
     procedure, pass(this) :: init => sx_gmres_init
     procedure, pass(this) :: free => sx_gmres_free
     procedure, pass(this) :: solve => sx_gmres_solve
     procedure, pass(this) :: solve_coupled => sx_gmres_solve_coupled
  end type sx_gmres_t
 
contains
 
  subroutine sx_gmres_init(this, n, max_iter, M, rel_tol, abs_tol, monitor)
    class(sx_gmres_t), target, intent(inout) :: this
    integer, intent(in) :: n
    integer, intent(in) :: max_iter
    class(pc_t), optional, intent(in), target :: M
    real(kind=rp), optional, intent(in) :: rel_tol
    real(kind=rp), optional, intent(in) :: abs_tol
    logical, optional, intent(in) :: monitor
 
    call this%free()
 
    if (present(m)) then
       this%M => m
    end if
 
    allocate(this%w(n))
    allocate(this%r(n))
    allocate(this%ml(n))
    allocate(this%mu(n))
    allocate(this%wk1(n))
 
    allocate(this%c(this%lgmres))
    allocate(this%s(this%lgmres))
    allocate(this%gam(this%lgmres + 1))
 
    allocate(this%z(n, this%lgmres))
    allocate(this%v(n, this%lgmres))
 
    allocate(this%h(this%lgmres, this%lgmres))
 
 
    if (present(rel_tol) .and. present(abs_tol) .and. present(monitor)) then
       call this%ksp_init(max_iter, rel_tol, abs_tol, monitor = monitor)
    else if (present(rel_tol) .and. present(abs_tol)) then
       call this%ksp_init(max_iter, rel_tol, abs_tol)
    else if (present(monitor) .and. present(abs_tol)) then
       call this%ksp_init(max_iter, abs_tol = abs_tol, monitor = monitor)
    else if (present(rel_tol) .and. present(monitor)) then
       call this%ksp_init(max_iter, rel_tol, monitor = monitor)
    else if (present(rel_tol)) then
       call this%ksp_init(max_iter, rel_tol = rel_tol)
    else if (present(abs_tol)) then
       call this%ksp_init(max_iter, abs_tol = abs_tol)
    else if (present(monitor)) then
       call this%ksp_init(max_iter, monitor = monitor)
    else
       call this%ksp_init(max_iter)
    end if
 
  end subroutine sx_gmres_init
 
  subroutine sx_gmres_free(this)
    class(sx_gmres_t), intent(inout) :: this
 
    call this%ksp_free()
 
    if (allocated(this%w)) then
       deallocate(this%w)
    end if
 
    if (allocated(this%c)) then
       deallocate(this%c)
    end if
 
    if (allocated(this%r)) then
       deallocate(this%r)
    end if
 
    if (allocated(this%z)) then
       deallocate(this%z)
    end if
 
    if (allocated(this%h)) then
       deallocate(this%h)
    end if
 
    if (allocated(this%ml)) then
       deallocate(this%ml)
    end if
 
    if (allocated(this%v)) then
       deallocate(this%v)
    end if
 
    if (allocated(this%s)) then
       deallocate(this%s)
    end if
 
    if (allocated(this%mu)) then
       deallocate(this%mu)
    end if
 
    if (allocated(this%gam)) then
       deallocate(this%gam)
    end if
 
    if (allocated(this%wk1)) then
       deallocate(this%wk1)
    end if
 
    nullify(this%M)
 
  end subroutine sx_gmres_free
 
  function sx_gmres_solve(this, Ax, x, f, n, coef, bc_projector, gs_h, niter) &
       result(ksp_results)
    class(sx_gmres_t), intent(inout) :: this
    class(ax_t), intent(in) :: ax
    type(field_t), intent(inout) :: x
    integer, intent(in) :: n
    real(kind=rp), dimension(n), intent(in) :: f
    type(coef_t), intent(inout) :: coef
    class(scalar_bc_projector_t), intent(inout) :: bc_projector
    type(gs_t), intent(inout) :: gs_h
    type(ksp_monitor_t) :: ksp_results
    integer, optional, intent(in) :: niter
    integer :: iter, max_iter, glb_n
    integer :: i, j, k, ierr
    real(kind=rp), parameter :: one = 1.0
    real(kind=rp) :: rnorm
    real(kind=rp) :: alpha, temp, l
    real(kind=rp) :: ratio, div0, norm_fac
    logical :: conv
    integer outer
 
    conv = .false.
    iter = 0
    rnorm = 0.0_rp
    glb_n = n / x%msh%nelv * x%msh%glb_nelv
 
    if (present(niter)) then
       max_iter = niter
    else
       max_iter = this%max_iter
    end if
 
    call rone(this%ml, n)
    call rone(this%mu, n)
    norm_fac = one / sqrt(coef%volume)
    call rzero(x%x, n)
    call rzero(this%gam, this%lgmres + 1)
    call rone(this%s, this%lgmres)
    call rone(this%c, this%lgmres)
    call rzero(this%h, this%lgmres * this%lgmres)
    outer = 0
    call this%monitor_start('GMRES')
    do while (.not. conv .and. iter .lt. max_iter)
       outer = outer + 1
 
       if (iter .eq. 0) then
          call col3(this%r, this%ml, f, n)
       else
          !update residual
          call copy (this%r,f,n)
          call ax%compute(this%w, x%x, coef, x%msh, x%Xh)
          call gs_h%op(this%w, n, gs_op_add)
          call bc_projector%apply(this%w, n)
          call add2s2(this%r, this%w, -one, n)
          call col2(this%r, this%ml, n)
       end if
       this%gam(1) = sqrt(glsc3(this%r, this%r, coef%mult, n))
       if (iter .eq. 0) then
          div0 = this%gam(1) * norm_fac
          ksp_results%res_start = div0
       end if
 
       if (abscmp(this%gam(1), 0.0_rp)) exit
 
       rnorm = 0.0_rp
       temp = one / this%gam(1)
       call cmult2(this%v(1,1), this%r, temp, n)
       do j = 1, this%lgmres
          iter = iter+1
          call col3(this%w, this%mu, this%v(1,j), n)
 
          !Apply precond
          call this%M%solve(this%z(1,j), this%w, n)
 
          call ax%compute(this%w, this%z(1,j), coef, x%msh, x%Xh)
          call gs_h%op(this%w, n, gs_op_add)
          call bc_projector%apply(this%w, n)
          call col2(this%w, this%ml, n)
 
          do i = 1, j
             this%h(i,j) = 0.0_rp
             do k = 1, n
                this%h(i,j) = this%h(i,j) + &
                     this%w(k) * this%v(k,i) * coef%mult(k,1,1,1)
             end do
          end do
 
          !Could probably be done inplace...
          call mpi_allreduce(this%h(1,j), this%wk1, j, &
               mpi_real_precision, mpi_sum, neko_comm, ierr)
          call copy(this%h(1,j), this%wk1, j)
 
          do i = 1, j
             do k = 1, n
                this%w(k) = this%w(k) - this%h(i,j) * this%v(k,i)
             end do
          end do
 
          !apply Givens rotations to new column
          do i = 1, j-1
             temp = this%h(i,j)
             this%h(i ,j) = this%c(i)*temp + this%s(i)*this%h(i+1,j)
             this%h(i+1,j) = -this%s(i)*temp + this%c(i)*this%h(i+1,j)
          end do
 
          alpha = sqrt(glsc3(this%w, this%w, coef%mult, n))
          rnorm = 0.0_rp
          if (abscmp(alpha, 0.0_rp)) then
             conv = .true.
             exit
          end if
          l = sqrt(this%h(j,j) * this%h(j,j) + alpha**2)
          temp = one / l
          this%c(j) = this%h(j,j) * temp
          this%s(j) = alpha * temp
          this%h(j,j) = l
          this%gam(j+1) = -this%s(j) * this%gam(j)
          this%gam(j) = this%c(j) * this%gam(j)
 
          rnorm = abs(this%gam(j+1)) * norm_fac
          call this%monitor_iter(iter, rnorm)
          ratio = rnorm / div0
          if (rnorm .lt. this%abs_tol) then
             conv = .true.
             exit
          end if
 
          if (iter + 1 .gt. max_iter) exit
 
          if (j .lt. this%lgmres) then
             temp = one / alpha
             call cmult2(this%v(1,j+1), this%w, temp, n)
          end if
       end do
       j = min(j, this%lgmres)
       !back substitution
       do k = j, 1, -1
          temp = this%gam(k)
          do i = j, k+1, -1
             temp = temp - this%h(k,i) * this%c(i)
          end do
          this%c(k) = temp / this%h(k,k)
       end do
       !sum up Arnoldi vectors
       do i = 1, j
          do k = 1, n
             x%x(k,1,1,1) = x%x(k,1,1,1) + this%c(i) * this%z(k,i)
          end do
       end do
    end do
    call this%monitor_stop()
    ksp_results%res_final = rnorm
    ksp_results%iter = iter
    ksp_results%converged = this%is_converged(iter, rnorm)
  end function sx_gmres_solve
 
  function sx_gmres_solve_coupled(this, Ax, x, y, z, fx, fy, fz, &
       n, coef, bc_projector, gs_h, niter) result(ksp_results)
    class(sx_gmres_t), intent(inout) :: this
    class(ax_t), intent(in) :: ax
    type(field_t), intent(inout) :: x
    type(field_t), intent(inout) :: y
    type(field_t), intent(inout) :: z
    integer, intent(in) :: n
    real(kind=rp), dimension(n), intent(in) :: fx
    real(kind=rp), dimension(n), intent(in) :: fy
    real(kind=rp), dimension(n), intent(in) :: fz
    type(coef_t), intent(inout) :: coef
    class(vector_bc_projector_t), intent(inout) :: bc_projector
    type(gs_t), intent(inout) :: gs_h
    type(ksp_monitor_t), dimension(3) :: ksp_results
    integer, optional, intent(in) :: niter
    type(scalar_bc_projector_t), pointer :: bc_x, bc_y, bc_z
 
    call vector_bc_projector_components(bc_projector, bc_x, bc_y, bc_z)
    ksp_results(1) = this%solve(ax, x, fx, n, coef, bc_x, gs_h, niter)
    ksp_results(2) = this%solve(ax, y, fy, n, coef, bc_y, gs_h, niter)
    ksp_results(3) = this%solve(ax, z, fz, n, coef, bc_z, gs_h, niter)
 
  end function sx_gmres_solve_coupled
 
end module gmres_sx