pipecg_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 pipecg_sx
  use krylov, only : ksp_t, ksp_monitor_t, ksp_max_iter
  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 bc, only : bc_list_t, bc_list_apply
  use math, only : glsc3, abscmp
  use comm
  implicit none
  private
 
  type, public, extends(ksp_t) :: sx_pipecg_t
     real(kind=rp), allocatable :: p(:)
     real(kind=rp), allocatable :: q(:)
     real(kind=rp), allocatable :: r(:)
     real(kind=rp), allocatable :: s(:)
     real(kind=rp), allocatable :: u(:)
     real(kind=rp), allocatable :: w(:)
     real(kind=rp), allocatable :: z(:)
     real(kind=rp), allocatable :: mi(:)
     real(kind=rp), allocatable :: ni(:)
   contains
     procedure, pass(this) :: init => sx_pipecg_init
     procedure, pass(this) :: free => sx_pipecg_free
     procedure, pass(this) :: solve => sx_pipecg_solve
     procedure, pass(this) :: solve_coupled => sx_pipecg_solve_coupled
  end type sx_pipecg_t
 
contains
 
  subroutine sx_pipecg_init(this, n, max_iter, M, rel_tol, abs_tol, monitor)
    class(sx_pipecg_t), intent(inout) :: this
    class(pc_t), optional, intent(in), target :: M
    integer, intent(in) :: n
    integer, intent(in) :: max_iter
    real(kind=rp), optional, intent(in) :: rel_tol
    real(kind=rp), optional, intent(in) :: abs_tol
    logical, optional, intent(in) :: monitor
 
    call this%free()
 
    allocate(this%p(n))
    allocate(this%q(n))
    allocate(this%r(n))
    allocate(this%s(n))
    allocate(this%u(n))
    allocate(this%w(n))
    allocate(this%z(n))
    allocate(this%mi(n))
    allocate(this%ni(n))
    if (present(m)) then
       this%M => m
    end if
 
    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_pipecg_init
 
  subroutine sx_pipecg_free(this)
    class(sx_pipecg_t), intent(inout) :: this
 
    call this%ksp_free()
 
    if (allocated(this%p)) then
       deallocate(this%p)
    end if
    if (allocated(this%q)) then
       deallocate(this%q)
    end if
    if (allocated(this%r)) then
       deallocate(this%r)
    end if
    if (allocated(this%s)) then
       deallocate(this%s)
    end if
    if (allocated(this%u)) then
       deallocate(this%u)
    end if
    if (allocated(this%w)) then
       deallocate(this%w)
    end if
    if (allocated(this%z)) then
       deallocate(this%z)
    end if
    if (allocated(this%mi)) then
       deallocate(this%mi)
    end if
    if (allocated(this%ni)) then
       deallocate(this%ni)
    end if
 
    nullify(this%M)
 
 
  end subroutine sx_pipecg_free
 
  function sx_pipecg_solve(this, Ax, x, f, n, coef, blst, gs_h, niter) result(ksp_results)
    class(sx_pipecg_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
    type(bc_list_t), intent(in) :: blst
    type(gs_t), intent(inout) :: gs_h
    type(ksp_monitor_t) :: ksp_results
    integer, optional, intent(in) :: niter
    integer :: iter, max_iter, i, ierr
    real(kind=rp) :: rnorm, rtr, reduction(3), norm_fac
    real(kind=rp) :: alpha, beta, gamma1, gamma2, delta
    real(kind=rp) :: tmp1, tmp2, tmp3
    type(mpi_request) :: request
    type(mpi_status) :: status
 
    if (present(niter)) then
       max_iter = niter
    else
       max_iter = this%max_iter
    end if
    norm_fac = 1.0_rp / sqrt(coef%volume)
 
    do i = 1, n
       x%x(i,1,1,1) = 0.0_rp
       this%z(i) = 0.0_rp
       this%q(i) = 0.0_rp
       this%p(i) = 0.0_rp
       this%s(i) = 0.0_rp
       this%r(i) = f(i)
    end do
 
    call this%M%solve(this%u, this%r, n)
    call ax%compute(this%w, this%u, coef, x%msh, x%Xh)
    call gs_h%op(this%w, n, gs_op_add)
    call bc_list_apply(blst, this%w, n)
 
    rtr = glsc3(this%r, coef%mult, this%r, n)
    rnorm = sqrt(rtr)*norm_fac
    ksp_results%res_start = rnorm
    ksp_results%res_final = rnorm
    ksp_results%iter = 0
    if(abscmp(rnorm, 0.0_rp)) return
 
    gamma1 = 0.0_rp
 
    call this%monitor_start('PipeCG')
    do iter = 1, max_iter
 
       tmp1 = 0.0_rp
       tmp2 = 0.0_rp
       tmp3 = 0.0_rp
       do i = 1, n
          tmp1 = tmp1 + this%r(i) * coef%mult(i,1,1,1) * this%u(i)
          tmp2 = tmp2 + this%w(i) * coef%mult(i,1,1,1) * this%u(i)
          tmp3 = tmp3 + this%r(i) * coef%mult(i,1,1,1) * this%r(i)
       end do
       reduction(1) = tmp1
       reduction(2) = tmp2
       reduction(3) = tmp3
 
       call mpi_iallreduce(mpi_in_place, reduction, 3, &
            mpi_real_precision, mpi_sum, neko_comm, request, ierr)
 
       call this%M%solve(this%mi, this%w, n)
       call ax%compute(this%ni, this%mi, coef, x%msh, x%Xh)
       call gs_h%op(this%ni, n, gs_op_add)
       call bc_list_apply(blst, this%ni, n)
 
       call mpi_wait(request, status, ierr)
       gamma2 = gamma1
       gamma1 = reduction(1)
       delta = reduction(2)
       rtr = reduction(3)
 
       rnorm = sqrt(rtr)*norm_fac
       call this%monitor_iter(iter, rnorm)
       if (rnorm .lt. this%abs_tol) then
          exit
       end if
 
       if (iter .gt. 1) then
          beta = gamma1 / gamma2
          alpha = gamma1 / (delta - (beta * gamma1/alpha))
       else
          beta = 0.0_rp
          alpha = gamma1/delta
       end if
 
       do i = 1, n
          this%z(i) = beta * this%z(i) + this%ni(i)
          this%q(i) = beta * this%q(i) + this%mi(i)
          this%s(i) = beta * this%s(i) + this%w(i)
          this%p(i) = beta * this%p(i) + this%u(i)
       end do
 
       do i = 1, n
          x%x(i,1,1,1) = x%x(i,1,1,1) + alpha * this%p(i)
          this%r(i) = this%r(i) - alpha * this%s(i)
          this%u(i) = this%u(i) - alpha * this%q(i)
          this%w(i) = this%w(i) - alpha * this%z(i)
       end do
 
    end do
    call this%monitor_stop()
    ksp_results%res_final = rnorm
    ksp_results%iter = iter
 
  end function sx_pipecg_solve
 
  function sx_pipecg_solve_coupled(this, Ax, x, y, z, fx, fy, fz, &
       n, coef, blstx, blsty, blstz, gs_h, niter) result(ksp_results)
    class(sx_pipecg_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
    type(bc_list_t), intent(in) :: blstx
    type(bc_list_t), intent(in) :: blsty
    type(bc_list_t), intent(in) :: blstz
    type(gs_t), intent(inout) :: gs_h
    type(ksp_monitor_t), dimension(3) :: ksp_results
    integer, optional, intent(in) :: niter
 
    ksp_results(1) =  this%solve(ax, x, fx, n, coef, blstx, gs_h, niter)
    ksp_results(2) =  this%solve(ax, y, fy, n, coef, blsty, gs_h, niter)
    ksp_results(3) =  this%solve(ax, z, fz, n, coef, blstz, gs_h, niter)
 
  end function sx_pipecg_solve_coupled
 
end module pipecg_sx