cg_coupled.f90

Go to the documentation of this file.

! Copyright (c) 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 cg_cpld
  use num_types, only: rp, xp
  use krylov, only : ksp_t, ksp_monitor_t, ksp_max_iter
  use precon, only : pc_t
  use ax_product, only : ax_t
  use field, only : field_t
  use coefs, only : coef_t
  use gather_scatter, only : gs_t, gs_op_add
  use bc_list, only : bc_list_t
  use math, only : glsc2, abscmp
  use comm, only : mpi_extra_precision, neko_comm
  use mpi_f08, only : mpi_allreduce, mpi_in_place, mpi_sum
  use utils, only : neko_error
  use operators, only : rotate_cyc
  implicit none
  private
 
  type, public, extends(ksp_t) :: cg_cpld_t
     real(kind=rp), allocatable :: w1(:)
     real(kind=rp), allocatable :: w2(:)
     real(kind=rp), allocatable :: w3(:)
     real(kind=rp), allocatable :: r1(:)
     real(kind=rp), allocatable :: r2(:)
     real(kind=rp), allocatable :: r3(:)
     real(kind=rp), allocatable :: p1(:)
     real(kind=rp), allocatable :: p2(:)
     real(kind=rp), allocatable :: p3(:)
     real(kind=rp), allocatable :: z1(:)
     real(kind=rp), allocatable :: z2(:)
     real(kind=rp), allocatable :: z3(:)
     real(kind=rp), allocatable :: tmp(:)
   contains
     procedure, pass(this) :: init => cg_cpld_init
     procedure, pass(this) :: free => cg_cpld_free
     procedure, pass(this) :: solve => cg_cpld_nop
     procedure, pass(this) :: solve_coupled => cg_cpld_solve
  end type cg_cpld_t
 
contains
 
  subroutine cg_cpld_init(this, n, max_iter, M, rel_tol, abs_tol, monitor)
    class(cg_cpld_t), target, intent(inout) :: this
    integer, intent(in) :: max_iter
    class(pc_t), optional, intent(in), target :: M
    integer, intent(in) :: n
    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%w1(n))
    allocate(this%w2(n))
    allocate(this%w3(n))
    allocate(this%r1(n))
    allocate(this%r2(n))
    allocate(this%r3(n))
    allocate(this%p1(n))
    allocate(this%p2(n))
    allocate(this%p3(n))
    allocate(this%z1(n))
    allocate(this%z2(n))
    allocate(this%z3(n))
    allocate(this%tmp(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 cg_cpld_init
 
  subroutine cg_cpld_free(this)
    class(cg_cpld_t), intent(inout) :: this
 
    call this%ksp_free()
 
    if (allocated(this%w1)) then
       deallocate(this%w1)
    end if
 
    if (allocated(this%w2)) then
       deallocate(this%w2)
    end if
 
    if (allocated(this%w3)) then
       deallocate(this%w3)
    end if
 
    if (allocated(this%r1)) then
       deallocate(this%r1)
    end if
 
    if (allocated(this%r2)) then
       deallocate(this%r2)
    end if
 
    if (allocated(this%r3)) then
       deallocate(this%r3)
    end if
 
    if (allocated(this%p1)) then
       deallocate(this%p1)
    end if
 
    if (allocated(this%p2)) then
       deallocate(this%p2)
    end if
 
    if (allocated(this%p3)) then
       deallocate(this%p3)
    end if
 
    if (allocated(this%z1)) then
       deallocate(this%z1)
    end if
 
    if (allocated(this%z2)) then
       deallocate(this%z2)
    end if
 
    if (allocated(this%z3)) then
       deallocate(this%z3)
    end if
 
    if (allocated(this%tmp)) then
       deallocate(this%tmp)
    end if
 
    nullify(this%M)
 
  end subroutine cg_cpld_free
 
  function cg_cpld_nop(this, Ax, x, f, n, coef, blst, gs_h, niter) &
       result(ksp_results)
    class(cg_cpld_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(inout) :: blst
    type(gs_t), intent(inout) :: gs_h
    type(ksp_monitor_t) :: ksp_results
    integer, optional, intent(in) :: niter
 
    ! Throw and error
    call neko_error('The cpldcg solver is only defined for coupled solves')
 
    ksp_results%res_final = 0.0
    ksp_results%iter = 0
  end function cg_cpld_nop
 
  function cg_cpld_solve(this, Ax, x, y, z, fx, fy, fz, &
       n, coef, blstx, blsty, blstz, gs_h, niter) result(ksp_results)
    class(cg_cpld_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(inout) :: blstx
    type(bc_list_t), intent(inout) :: blsty
    type(bc_list_t), intent(inout) :: blstz
    type(gs_t), intent(inout) :: gs_h
    type(ksp_monitor_t), dimension(3) :: ksp_results
    integer, optional, intent(in) :: niter
    integer :: i, iter, max_iter, ierr
    real(kind=rp) :: rnorm, rtr, rtr0, rtz2, rtz1
    real(kind=rp) :: beta, pap, alpha, norm_fac
    real(kind=xp) :: tmp_xp, r1_xp, r2_xp, r3_xp, mult_xp
 
    if (present(niter)) then
       max_iter = niter
    else
       max_iter = this%max_iter
    end if
    norm_fac = 1.0_rp / sqrt(coef%volume)
 
    associate(p1 => this%p1, p2 => this%p2, p3 => this%p3, z1 => this%z1, &
         z2 => this%z2, z3 => this%z3, r1 => this%r1, r2 => this%r2, &
         r3 => this%r3, tmp => this%tmp, w1 => this%w1, w2 => this%w2, &
         w3 => this%w3)
 
      rtz1 = 1.0_rp
      do concurrent(i = 1:n)
         x%x(i,1,1,1) = 0.0_rp
         y%x(i,1,1,1) = 0.0_rp
         z%x(i,1,1,1) = 0.0_rp
         p1(i) = 0.0_rp
         p2(i) = 0.0_rp
         p3(i) = 0.0_rp
         z1(i) = 0.0_rp
         z2(i) = 0.0_rp
         z3(i) = 0.0_rp
         r1(i) = fx(i)
         r2(i) = fy(i)
         r3(i) = fz(i)
         tmp(i) = r1(i)**2 + r2(i)**2 + r3(i)**2
      end do
 
      rtr = glsc2(tmp, coef%mult, 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)) then
         ksp_results%converged = .true.
         return
      end if
 
      call this%monitor_start('cpldCG')
      do iter = 1, max_iter
         call this%M%solve(z1, this%r1, n)
         call this%M%solve(z2, this%r2, n)
         call this%M%solve(z3, this%r3, n)
         rtz2 = rtz1
 
         do concurrent(i = 1:n)
            this%tmp(i) = z1(i) * r1(i) &
                 + z2(i) * r2(i) &
                 + z3(i) * r3(i)
         end do
 
         rtz1 = glsc2(tmp, coef%mult, n)
 
         beta = rtz1 / rtz2
         if (iter .eq. 1) beta = 0.0_rp
         do concurrent(i = 1:n)
            p1(i) = p1(i) * beta + z1(i)
            p2(i) = p2(i) * beta + z2(i)
            p3(i) = p3(i) * beta + z3(i)
         end do
 
         call ax%compute_vector(w1, w2, w3, p1, p2, p3, coef, x%msh, x%Xh)
 
         call rotate_cyc(w1, w2, w3, 1, coef)
         call gs_h%op(w1, n, gs_op_add)
         call gs_h%op(w2, n, gs_op_add)
         call gs_h%op(w3, n, gs_op_add)
         call rotate_cyc(w1, w2, w3, 0, coef)
 
         call blstx%apply_scalar(w1, n)
         call blsty%apply_scalar(w2, n)
         call blstz%apply_scalar(w3, n)
 
         do concurrent(i = 1:n)
            tmp(i) = w1(i) * p1(i) &
                 + w2(i) * p2(i) &
                 + w3(i) * p3(i)
         end do
 
         pap = glsc2(tmp, coef%mult, n)
 
         alpha = rtz1 / pap
         do concurrent(i = 1:n)
            x%x(i,1,1,1) = x%x(i,1,1,1) + alpha * p1(i)
            y%x(i,1,1,1) = y%x(i,1,1,1) + alpha * p2(i)
            z%x(i,1,1,1) = z%x(i,1,1,1) + alpha * p3(i)
         end do
 
         tmp = 0.0_rp
         tmp_xp = 0.0_xp
         !$omp parallel do reduction(+:tmp_xp)
         do i = 1, n
            r1(i) = r1(i) - alpha * w1(i)
            r2(i) = r2(i) - alpha * w2(i)
            r3(i) = r3(i) - alpha * w3(i)
            r1_xp = real(r1(i), kind=xp)
            r2_xp = real(r2(i), kind=xp)
            r3_xp = real(r3(i), kind=xp)
            mult_xp = real(coef%mult(i,1,1,1), kind=xp)
            tmp_xp = tmp_xp + &
                 (r1_xp * r1_xp + r2_xp * r2_xp + r3_xp * r3_xp) * mult_xp
         end do
         !$omp end parallel do
 
         call mpi_allreduce(mpi_in_place, tmp_xp, 1, mpi_extra_precision, &
              mpi_sum, neko_comm, ierr)
         rtr = tmp_xp
 
         if (iter .eq. 1) rtr0 = rtr
         rnorm = sqrt(rtr) * norm_fac
         call this%monitor_iter(iter, rnorm)
         if (rnorm .lt. this%abs_tol) then
            exit
         end if
      end do
    end associate
    call this%monitor_stop()
    ksp_results%res_final = rnorm
    ksp_results%iter = iter
    ksp_results%converged = this%is_converged(iter, rnorm)
  end function cg_cpld_solve
 
end module cg_cpld