leq_bicgs.f

Go to the documentation of this file.

!vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvC
!                                                                      C
!  Subroutine LEQ_BICGS                                                C
!  Purpose: Solve system of linear system using BICGS method           C
!           Biconjugate gradients stabilized                           C
!                                                                      C
!  Author: Ed D'Azevedo                               Date: 21-JAN-99  C
!  Reviewer:                                          Date:            C
!                                                                      C
!  Literature/Document References:                                     C
!  Variables referenced:                                               C
!  Variables modified:                                                 C
!  Local variables:                                                    C
!                                                                      C
!^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^C
! Handan Liu wrote below:               !Jan 22 2013
! The modification is as below:
!       Adding a loop of 2D RSRS sweep and parallelizing for OpenMP.
!       Splitting the existing 3D RSRS loop into two loops for OpenMP
!               due to data dependency
!       Adding openmp directives in all loops in leq_bicgs0.
!
      SUBROUTINE leq_bicgs(VNAME, VNO, VAR, A_M, B_m, cmethod, &
                           tol, pc, itmax, ier)

!-----------------------------------------------
! Modules
!-----------------------------------------------
      USE param
      USE param1
      USE geometry
      USE compar
      USE indices
      USE leqsol
      USE funits
      IMPLICIT NONE
!-----------------------------------------------
! Dummy arguments
!-----------------------------------------------
! variable name
      CHARACTER(LEN=*), INTENT(IN) :: Vname
! variable number (not really used here; see calling subroutine)
      INTEGER, INTENT(IN) :: VNO
! variable
!     e.g., pp_g, epp, rop_g, rop_s, u_g, u_s, v_g, v_s, w_g,
!     w_s, T_g, T_s, x_g, x_s, Theta_m, scalar, K_Turb_G,
!     e_Turb_G
!      DOUBLE PRECISION, DIMENSION(ijkstart3:ijkend3), INTENT(INOUT) :: Var
      DOUBLE PRECISION, DIMENSION(DIMENSION_3), INTENT(INOUT) :: Var
! Septadiagonal matrix A_m
!      DOUBLE PRECISION, DIMENSION(ijkstart3:ijkend3,-3:3), INTENT(INOUT) :: A_m
     DOUBLE PRECISION, DIMENSION(DIMENSION_3,-3:3), INTENT(INOUT) :: A_m

! Vector b_m
!      DOUBLE PRECISION, DIMENSION(ijkstart3:ijkend3), INTENT(INOUT) :: B_m
      DOUBLE PRECISION, DIMENSION(DIMENSION_3), INTENT(INOUT) :: B_m
! Sweep direction of leq solver (leq_sweep)
!     e.g., options = 'isis', 'rsrs' (default), 'asas'
! Note: this setting only seems to matter when leq_pc='line'
      CHARACTER(LEN=*), INTENT(IN) :: CMETHOD
! convergence tolerance (generally leq_tol)
      DOUBLE PRECISION, INTENT(IN) :: TOL
! preconditioner (leq_pc)
!     options = 'line' (default), 'diag', 'none'
      CHARACTER(LEN=4), INTENT(IN) ::  PC
! maximum number of iterations (generally leq_it)
      INTEGER, INTENT(IN) :: ITMAX
! error indicator
      INTEGER, INTENT(INOUT) :: IER
!-------------------------------------------------
! Local Variables
!-------------------------------------------------

      if(pc.eq.'LINE') then   ! default
         call leq_bicgs0( vname, vno, var, a_m, b_m,  &
            cmethod, tol, itmax, leq_matvec, leq_msolve, .true., ier )
      elseif(pc.eq.'DIAG') then
         call leq_bicgs0( vname, vno, var, a_m, b_m,   &
            cmethod, tol, itmax, leq_matvec, leq_msolve1, .true., ier )
      elseif(pc.eq.'NONE') then
         call leq_bicgs0( vname, vno, var, a_m, b_m,   &
            cmethod, tol, itmax, leq_matvec, leq_msolve0, .false., ier )
      else
         IF(dmp_log)WRITE (unit_log,*) &
           'preconditioner option not found - check mfix.dat and readme'
         call mfix_exit(mype)
      endif

      return
      END SUBROUTINE leq_bicgs

!vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvC
!                                                                      C
!  Subroutine: LEQ_BICGS0                                              C
!  Purpose: Compute residual of linear system                          C
!                                                                      C
!  Author: Ed D'Azevedo                               Date: 21-JAN-99  C
!  Reviewer:                                          Date:            C
!                                                                      C
!  Literature/Document References:                                     C
!  Variables referenced:                                               C
!  Variables modified:                                                 C
!  Local variables:                                                    C
!                                                                      C
!^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^C

      SUBROUTINE leq_bicgs0(VNAME, VNO, VAR, A_M, B_m, cmethod, &
                            tol, itmax, matvec, msolve, use_pc, ier )

!-----------------------------------------------
! Modules
!-----------------------------------------------
      USE param
      USE param1
      USE geometry
      USE compar
      USE mpi_utility
      USE sendrecv
      USE indices
      USE leqsol
      USE cutcell
      USE functions

      IMPLICIT NONE
!-----------------------------------------------
! Dummy arguments/procedure
!-----------------------------------------------
! variable name
      CHARACTER(LEN=*), INTENT(IN) :: Vname
! variable number (not really used here-see calling subroutine)
      INTEGER, INTENT(IN) :: VNO
! variable
!     e.g., pp_g, epp, rop_g, rop_s, u_g, u_s, v_g, v_s, w_g,
!     w_s, T_g, T_s, x_g, x_s, Theta_m, scalar, K_Turb_G,
!     e_Turb_G
!      DOUBLE PRECISION, INTENT(INOUT) :: Var(ijkstart3:ijkend3)
      DOUBLE PRECISION, DIMENSION(DIMENSION_3), INTENT(INOUT) :: Var
! Septadiagonal matrix A_m
!      DOUBLE PRECISION, INTENT(INOUT) :: A_m(ijkstart3:ijkend3,-3:3)
      DOUBLE PRECISION, DIMENSION(DIMENSION_3,-3:3), INTENT(INOUT) :: A_m
! Vector b_m
!      DOUBLE PRECISION, INTENT(INOUT) :: B_m(ijkstart3:ijkend3)
      DOUBLE PRECISION, DIMENSION(DIMENSION_3), INTENT(INOUT) :: B_m
! Sweep direction of leq solver (leq_sweep)
!     e.g., options = 'isis', 'rsrs' (default), 'asas'
      CHARACTER(LEN=*), INTENT(IN) :: CMETHOD
! convergence tolerance (generally leq_tol)
      DOUBLE PRECISION, INTENT(IN) :: TOL
! maximum number of iterations (generally leq_it)
      INTEGER, INTENT(IN) :: ITMAX
! indicate whether to use preconditioner
      LOGICAL, INTENT(IN) :: USE_PC
! error indicator
      INTEGER, INTENT(INOUT) :: IER
! dummy arguments/procedures set as indicated
!     matvec->leq_matvec
! for preconditioner (leq_pc)
!    'line' msolve->leq_msolve  (default)
!    'diag' msolve->leq_msolve1
!    'none' msolve->leq_msolve0
!-----------------------------------------------
! Local parameters
!-----------------------------------------------
      INTEGER, PARAMETER :: idebugl = 0
      DOUBLE PRECISION, PARAMETER :: ratiotol = 0.2
      LOGICAL, PARAMETER :: do_unit_scaling = .true.
!-----------------------------------------------
! Local variables
!-----------------------------------------------

      DOUBLE PRECISION, DIMENSION(:), allocatable :: R,Rtilde, Tvec,V
      DOUBLE PRECISION, DIMENSION(:), allocatable, target :: P, P_preconditioned
      DOUBLE PRECISION, DIMENSION(:), allocatable, target :: Svec, Svec_preconditioned

      ! Phat points to either preconditioned value of P, or P itself (to avoid copying for efficiency)
      DOUBLE PRECISION, POINTER :: Phat(:), Shat(:)

      DOUBLE PRECISION, DIMENSION(0:ITMAX+1) :: &
                        alpha, beta, omega, rho
      DOUBLE PRECISION :: TxS, TxT, RtildexV, &
                          aijmax, oam
      DOUBLE PRECISION :: Rnorm, Rnorm0, Snorm, TOLMIN, pnorm
      LOGICAL :: isconverged
      INTEGER :: i, j, k, ijk
      INTEGER :: iter
      DOUBLE PRECISION, DIMENSION(2) :: TxS_TxT
!-----------------------------------------------

! Initialize the error flag.
      ier=0

! Scale matrix to have unit diagonal
! ---------------------------------------------------------------->>>
      if (do_unit_scaling) then

         IF(re_indexing) THEN  ! Loop only over active cells
!$omp parallel do default(shared) private(ijk,oam,aijmax)
            DO ijk = ijkstart3,ijkend3
               aijmax = maxval(abs(a_m(ijk,:)) )
               if (aijmax > 0.0)then
                  oam = one/aijmax
                  a_m(ijk,:) = a_m(ijk,:)*oam
                  b_m(ijk) = b_m(ijk)*oam
               else
                  ier = -2
               endif
            ENDDO

         ELSE

!$omp parallel do default(shared) private(ijk,i,j,k,oam,aijmax)
            do k = kstart2,kend2
               do i = istart2,iend2
                  do j = jstart2,jend2
                     ijk = funijk(i,j,k)
                     aijmax = maxval(abs(a_m(ijk,:)) )
                     if (aijmax > 0.0) then
                        oam = one/aijmax
                        a_m(ijk,:) = a_m(ijk,:)*oam
                        b_m(ijk) = b_m(ijk)*oam
                     else
                        ier = -2
                     endif
                  enddo
               enddo
            enddo

         ENDIF

      endif

! A singlular matrix was detected.
      if(ier /= 0) RETURN
! ----------------------------------------------------------------<<<



      allocate(r(dimension_3))
      allocate(rtilde(dimension_3))
      allocate(p(dimension_3))
      allocate(p_preconditioned(dimension_3))
      allocate(svec(dimension_3))
      allocate(svec_preconditioned(dimension_3))
      allocate(tvec(dimension_3))
      allocate(v(dimension_3))

! these scalars should not be necessary to initialize but done as failsafe
      rnorm = zero
      rnorm0 = zero
      snorm = zero
      pnorm = zero

! initializing
      alpha(:)  = zero
      beta(:)   = zero
      omega(:)  = zero
      rho(:)    = zero

!$omp parallel sections
      r(:) = zero
!$omp section
      rtilde(:) = zero
!$omp section
      p(:) = zero
!$omp section
      p_preconditioned(:) = zero
!$omp section
      svec(:) = zero
!$omp section
      svec_preconditioned(:) = zero
!$omp section
      tvec(:) = zero
!$omp section
      v(:) = zero
!$omp end parallel sections

      tolmin = epsilon( one )

! Compute initial residual (R = b-A*x) for Ax=b
!    assume initial guess in Var
!    rtilde = r
! ---------------------------------------------------------------->>>
      call matvec(vname, var, a_m, r)   ! returns R=A*Var

!$omp parallel workshare
      r(:) = b_m(:) - r(:)
!$omp end parallel workshare

      call send_recv(r,nlayers_bicgs)

      rnorm0 = sqrt( dot_product_par( r, r ) )

! determine an initial guess for the residual = residual + small random
! number (so it could be set to anything). note that since random_number
! is used to supply the guess, this line could potentially be the source
! of small differences between runs.  the random number is shifted below
! between -1 and 1 and then scaled by factor 1.0D-6*Rnorm0
      call random_number(rtilde(:))

! Shift random number array to be consistent with case when RE_INDEXING is .FALSE.
       IF(re_indexing) CALL shift_dp_array(rtilde)

!$omp parallel workshare
       rtilde(:) = r(:) + (2.0d0*rtilde(:)-1.0d0)*1.0d-6*rnorm0
!$omp end parallel workshare

      if (idebugl >= 1) then
         if(mype.eq.0) print*,'leq_bicgs, initial: ', vname,' resid ', rnorm0
      endif
! ----------------------------------------------------------------<<<


! Main loop
! ---------------------------------------------------------------->>>
      iter = 1
      do i=1,itmax

         rho(i-1) = dot_product_par( rtilde, r )

         if (rho(i-1) .eq. zero) then
            if(i /= 1)then
! Method fails
! --------------------------------
               ier = -2
            else
! Method converged.  residual is already zero
! --------------------------------
               ier = 0
            endif
            call send_recv(var,2)
            return
         endif ! rho(i-1).eq.0

         if (i .eq. 1) then
!$omp parallel workshare
            p(:) = r(:)
!$omp end parallel workshare
         else
            beta(i-1) = ( rho(i-1)/rho(i-2) )*( alpha(i-1) / omega(i-1) )
!$omp parallel workshare
            p(:) = r(:) + beta(i-1)*( p(:) - omega(i-1)*v(:) )
!$omp end parallel workshare
         endif ! i.eq.1

! Solve A*Phat(:) = P(:)
! V(:) = A*Phat(:)
! --------------------------------
         if (use_pc) then
            call msolve(vname, p, a_m, p_preconditioned, cmethod) ! returns P_preconditioned
            phat => p_preconditioned
         else
            phat => p
         endif

         call matvec(vname, phat, a_m, v)   ! returns V=A*Phat

         rtildexv = dot_product_par( rtilde, v )

! compute alpha
! --------------------------------
         alpha(i) = rho(i-1) / rtildexv

! compute Svec
! --------------------------------
!$omp parallel workshare
         svec(:) = r(:) - alpha(i) * v(:)
!$omp end parallel workshare

! Check norm of Svec(:); if small enough:
! set X(:) = X(:) + alpha(i)*Phat(:) and stop
! --------------------------------
         if(.not.minimize_dotproducts) then
            snorm = sqrt( dot_product_par( svec, svec ) )

            if (snorm <= tolmin) then
!$omp parallel workshare
               var(:) = var(:) + alpha(i)*phat(:)
!$omp end parallel workshare

! Recompute residual norm
! --------------------------------
               if (idebugl >= 1) then
                  call matvec(vname, var, a_m, r)   ! returns R=A*Var
!                  Rnorm = sqrt( dot_product_par( Var, Var ) )
!                  print*,'leq_bicgs, initial: ', Vname,' Vnorm ', Rnorm

!$omp parallel workshare
                  r(:) = b_m(:) - r(:)
!$omp end parallel workshare

                  rnorm = sqrt( dot_product_par( r, r ) )
               endif            ! idebugl >= 1
               isconverged = .true.
               EXIT
            endif               ! end if (Snorm <= TOLMIN)
         endif                  ! end if (.not.minimize_dotproducts)

! Solve A*Shat(:) = Svec(:)
! Tvec(:) = A*Shat(:)
! --------------------------------

         if (use_pc) then
            call msolve(vname, svec, a_m, svec_preconditioned, cmethod) ! returns S_preconditioned
            shat => svec_preconditioned
         else
            shat => svec
         endif

         call matvec( vname, shat, a_m, tvec )   ! returns Tvec=A*Shat

         if(.not.minimize_dotproducts) then
!!     $omp parallel sections
            txs = dot_product_par( tvec, svec )
!!     $omp section
            txt = dot_product_par( tvec, tvec )
!!     $omp end parallel sections
         else
            txs_txt = dot_product_par2(tvec, svec, tvec, tvec )
            txs = txs_txt(1)
            txt = txs_txt(2)
         endif

         IF(txt.eq.zero) txt = small_number

! compute omega
! --------------------------------
         omega(i) = txs / txt

! compute new guess for Var
! --------------------------------

!$omp parallel sections
            var(:) = var(:) + alpha(i)*phat(:) + omega(i)*shat(:)
!$omp section
            r(:) = svec(:) - omega(i)*tvec(:)
!$omp end parallel sections

! --------------------------------
         if(.not.minimize_dotproducts.or.(mod(iter,icheck_bicgs).eq.0)) then
            rnorm = sqrt( dot_product_par(r, r) )

            if (idebugl.ge.1) then
               if (mype.eq.pe_io) then
                  print*,'iter, Rnorm ', iter, rnorm, snorm
                  print*,'alpha(i), omega(i) ', alpha(i), omega(i)
                  print*,'TxS, TxT ', txs, txt
                  print*,'RtildexV, rho(i-1) ', rtildexv, rho(i-1)
               endif
            endif

!           call mfix_exit(myPE)

! Check convergence; continue if necessary
! for continuation, it is necessary that omega(i) .ne. 0
            isconverged = (rnorm <= tol*rnorm0)

            if (isconverged) then
               iter_tot(vno) = iter_tot(vno) + iter + 1
               EXIT
            endif
         endif                  ! end if(.not.minimize_dotproducts)

! Advance the iteration count
         iter = iter + 1

      enddo   ! end do i=1,itmax
! end of linear solver loop
! ----------------------------------------------------------------<<<

      if (idebugl >= 1) then
         call matvec(vname, var, a_m, r)   ! returns R=A*Var

!$omp parallel workshare
         r(:) = r(:) - b_m(:)
!$omp end parallel workshare

         rnorm = sqrt( dot_product_par( r,r) )

         if(mype.eq.0) print*,'leq_bicgs: final Rnorm ', rnorm

         if(mype.eq.0)  print*,'leq_bicgs ratio : ', vname,' ',iter,  &
         ' L-2', rnorm/rnorm0
      endif   ! end if(idebugl >=1)

!      isconverged = (real(Rnorm) <= TOL*Rnorm0);
      if(.NOT.isconverged) isconverged = (real(Rnorm) <= tol*rnorm0);
!     write(*,*) '***',iter, isconverged, Rnorm, TOL, Rnorm0, myPE
      ier = 0
      if (.not.isconverged) then
         ier = -1
         iter_tot(vno) = iter_tot(vno) + iter
         if (real(Rnorm) >= ratiotol*real(rnorm0)) then
            ier = -2
         endif
      endif

      call send_recv(var,2)

      deallocate(r)
      deallocate(rtilde)
      deallocate(p)
      deallocate(p_preconditioned)
      deallocate(svec)
      deallocate(svec_preconditioned)
      deallocate(tvec)
      deallocate(v)

      return
      end subroutine leq_bicgs0