C $Header: /u/gcmpack/MITgcm/pkg/streamice/streamice_cg_solve_petsc.F,v 1.5 2015/06/15 14:34:38 dgoldberg Exp $
C $Name:  $

#include "STREAMICE_OPTIONS.h"

C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|

CBOP
      SUBROUTINE STREAMICE_CG_SOLVE_PETSC( 
     U                               cg_Uin,     ! x-velocities
     U                               cg_Vin,     ! y-velocities
     I                               cg_Bu,      ! force in x dir
     I                               cg_Bv,      ! force in y dir
     I                               A_uu,       ! section of matrix that multiplies u and projects on u
     I                               A_uv,       ! section of matrix that multiplies v and projects on u
     I                               A_vu,       ! section of matrix that multiplies u and projects on v
     I                               A_vv,       ! section of matrix that multiplies v and projects on v
     I                               tolerance, 
     O                               iters,
     I                               maxIter,
     I                               myThid )

C     /============================================================\
C     | SUBROUTINE                                                 |   
C     | o                                                          |
C     |============================================================|
C     |                                                            |
C     \============================================================/
      IMPLICIT NONE

#include "SIZE.h"
#include "EEPARAMS.h"
#include "PARAMS.h"
#include "STREAMICE.h"
#include "STREAMICE_CG.h"



#ifdef ALLOW_PETSC
#include "finclude/petsc.h"
#include "STREAMICE_PETSC.h"
! UNCOMMENT IF V3.0
!#include "finclude/petscvec.h"
!#include "finclude/petscmat.h"
!#include "finclude/petscksp.h"
!#include "finclude/petscpc.h"
#endif
C     === Global variables ===

      
C     !INPUT/OUTPUT ARGUMENTS
C     cg_Uin, cg_Vin - input and output velocities
C     cg_Bu, cg_Bv - driving stress
      INTEGER myThid
      INTEGER iters, maxiter
      _RL tolerance
      _RL cg_Uin (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
      _RL cg_Vin (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
      _RL cg_Bu (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
      _RL cg_Bv (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
      _RL 
     & A_uu (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy,-1:1,-1:1),
     & A_vu (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy,-1:1,-1:1),
     & A_uv (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy,-1:1,-1:1),
     & A_vv (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy,-1:1,-1:1)

#if (defined (ALLOW_OPENAD)  defined (ALLOW_STREAMICE_OAD_FP))
       LOGICAL create_mat, destroy_mat
#endif

C     LOCAL VARIABLES
      INTEGER i, j, bi, bj, cg_halo, conv_flag
      INTEGER iter, is, js, ie, je, colx, coly, k
      _RL dot_p1, dot_p2, alpha_k, beta_k, resid, resid_0
      _RL dot_p1_tile (nSx,nSy)
      _RL dot_p2_tile (nSx,nSy)
      CHARACTER*(MAX_LEN_MBUF) msgBuf


#ifdef ALLOW_PETSC
      INTEGER indices(2*(snx*nsx*sny*nsy))
      INTEGER n_dofs_cum_sum (0:nPx*nPy-1), idx(1)
      _RL rhs_values(2*(snx*nsx*sny*nsy))
      _RL solution_values(2*(snx*nsx*sny*nsy))
!      _RL mat_values (2*Nx*Ny,2*(snx*nsx*sny*nsy))
      _RL mat_values (18,1), mat_val_return(1)
      INTEGER indices_col(18)
      INTEGER local_dofs, global_dofs, dof_index, dof_index_col
      INTEGER local_offset
!      Mat matrix
!      KSP ksp
!      PC  pc
      PC  subpc
      KSP subksp(1)
      Vec rhs
      Vec solution
      PetscErrorCode ierr
#ifdef ALLOW_USE_MPI
      integer mpiRC, mpiMyWid
#endif 
#endif


#ifdef ALLOW_STREAMICE




#ifdef ALLOW_PETSC

#ifdef ALLOW_USE_MPI


      CALL MPI_COMM_RANK( MPI_COMM_WORLD, mpiMyWId, mpiRC )
      local_dofs = n_dofs_process (mpiMyWid)
      global_dofs = 0
      
      n_dofs_cum_sum(0) = 0
      DO i=0,nPx*nPy-1
       global_dofs = global_dofs + n_dofs_process (i)
       if (i.ge.1) THEN
         n_dofs_cum_sum(i) = n_dofs_cum_sum(i-1)+
     &                     n_dofs_process(i-1)
       endif
      ENDDO
      local_offset = n_dofs_cum_sum(mpimywid)

#else

      local_dofs = n_dofs_process (0)
      global_dofs = local_dofs
      local_offset = 0

#endif 

!      call petscInitialize(PETSC_NULL_CHARACTER,ierr)

!----------------------

      CALL TIMER_START ('STREAMICE_PETSC_SETUP',myThid)

      call VECCREATE(PETSC_COMM_WORLD, rhs, ierr)
      call VECSETSIZES(rhs, local_dofs, global_dofs, ierr)
      call VECSETTYPE(rhs, VECMPI, ierr)

      call VECCREATE(PETSC_COMM_WORLD, solution, ierr)
      call VECSETSIZES(solution, local_dofs, global_dofs, ierr) 
      call VECSETTYPE(solution, VECMPI, ierr)

      do i=1,local_dofs
        indices(i) = i-1 + local_offset
      end


do do i=1,2*nSx*nSy*sNx*sNy rhs_values (i) = 0. _d 0 solution_values (i) = 0. _d 0 enddo ! gather rhs and initial guess values to populate petsc vectors DO bj = myByLo(myThid), myByHi(myThid) DO bi = myBxLo(myThid), myBxHi(myThid) DO j=1,sNy DO i=1,sNx dof_index = INT(streamice_petsc_dofs_u(i,j,bi,bj)) & - local_offset if (dof_index.ge.0) THEN rhs_values(dof_index+1) = cg_Bu(i,j,bi,bj) solution_values(dof_index+1) = cg_Uin(i,j,bi,bj) endif !--------------- dof_index = INT(streamice_petsc_dofs_v(i,j,bi,bj)) & - local_offset if (dof_index.ge.0) THEN rhs_values(dof_index+1) = cg_Bv(i,j,bi,bj) solution_values(dof_index+1) = cg_Vin(i,j,bi,bj) endif ENDDO ENDDO ENDDO ENDDO call VECSETVALUES(rhs, local_dofs, indices, rhs_values, & INSERT_VALUES, ierr) call VECASSEMBLYBEGIN(rhs, ierr) call VECASSEMBLYEND(rhs, ierr) call VECSETVALUES(solution, local_dofs, indices, & solution_values, INSERT_VALUES, ierr) call VECASSEMBLYBEGIN(solution, ierr) call VECASSEMBLYEND(solution, ierr) #if (defined (ALLOW_OPENAD) defined (ALLOW_STREAMICE_OAD_FP)) #ifdef ALLOW_PETSC if (STREAMICE_need2createmat) then #endif #endif ! IF USING v3.0 THEN ! call MatCreateMPIAIJ (PETSC_COMM_WORLD, call MATCREATEAIJ (PETSC_COMM_WORLD, & local_dofs, local_dofs, & global_dofs, global_dofs, & 18, PETSC_NULL_INTEGER, & 18, PETSC_NULL_INTEGER, & matrix, ierr) ! populate petsc matrix DO bj = myByLo(myThid), myByHi(myThid) DO bi = myBxLo(myThid), myBxHi(myThid) DO j=1,sNy DO i=1,sNx dof_index = INT(streamice_petsc_dofs_u(i,j,bi,bj)) ! & - local_offset IF (dof_index .ge. 0) THEN DO k=1,18 indices_col(k) = 0 mat_values(k,1) = 0. _d 0 ENDDO k=0 DO coly=-1,1 DO colx=-1,1 dof_index_col = streamice_petsc_dofs_u(i+colx,j+coly,bi,bj) if (dof_index_col.ge.0) THEN ! pscal = A_uu(i,j,bi,bj,colx,coly) ! CALL MatSetValue (matrix,dof_index, dof_index_col, ! & pscal,INSERT_VALUES,ierr) k=k+1 mat_values (k,1) = A_uu(i,j,bi,bj,colx,coly) indices_col (k) = dof_index_col endif dof_index_col = streamice_petsc_dofs_v(i+colx,j+coly,bi,bj) if (dof_index_col.ge.0) THEN ! CALL MatSetValue (matrix,dof_index, dof_index_col, ! & A_uv(i,j,bi,bj,colx,coly),INSERT_VALUES,ierr) k=k+1 mat_values (k,1) = A_uv(i,j,bi,bj,colx,coly) indices_col (k) = dof_index_col endif ENDDO ENDDO call MATSETVALUES (matrix, 1, dof_index, k, indices_col, & mat_values,INSERT_VALUES,ierr) ENDIF ! ---------------------------------------------- dof_index = INT(streamice_petsc_dofs_v(i,j,bi,bj)) ! & - local_offset IF (dof_index .ge. 0) THEN DO k=1,18 indices_col(k) = 0 mat_values(k,1) = 0. _d 0 ENDDO k=0 DO coly=-1,1 DO colx=-1,1 dof_index_col = streamice_petsc_dofs_u(i+colx,j+coly,bi,bj) if (dof_index_col.ge.0) THEN ! CALL MatSetValue (matrix,dof_index, dof_index_col, ! & A_vu(i,j,bi,bj,colx,coly),INSERT_VALUES,ierr) k=k+1 mat_values (k,1) = A_vu(i,j,bi,bj,colx,coly) indices_col (k) = dof_index_col endif dof_index_col = streamice_petsc_dofs_v(i+colx,j+coly,bi,bj) if (dof_index_col.ge.0) THEN ! CALL MatSetValue (matrix,dof_index, dof_index_col, ! & A_vv(i,j,bi,bj,colx,coly),INSERT_VALUES,ierr) k=k+1 mat_values (k,1) = A_vv(i,j,bi,bj,colx,coly) indices_col (k) = dof_index_col endif ENDDO ENDDO call MATSETVALUES (matrix, 1, dof_index, k, indices_col, & mat_values,INSERT_VALUES,ierr) ENDIF ENDDO ENDDO ENDDO ENDDO call MATASSEMBLYBEGIN(matrix,MAT_FINAL_ASSEMBLY,ierr) call MATASSEMBLYEND(matrix,MAT_FINAL_ASSEMBLY,ierr) call KSPCREATE(PETSC_COMM_WORLD, ksp, ierr) call KSPSETOPERATORS(ksp, matrix, matrix, & DIFFERENT_NONZERO_PATTERN, ierr) IF (PETSC_PRECOND_TYPE.eq.'MUMPS') then call KSPSETTYPE(ksp,KSPPREONLY,ierr) ELSE SELECT CASE (PETSC_SOLVER_TYPE) CASE ('CG') PRINT *, "PETSC SOLVER: SELECTED CG" call KSPSETTYPE(ksp, KSPCG, ierr) CASE ('GMRES') PRINT *, "PETSC SOLVER: SELECTED GMRES" call KSPSETTYPE(ksp, KSPGMRES, ierr) CASE ('BICG') PRINT *, "PETSC SOLVER: SELECTED BICG" call KSPSETTYPE(ksp, KSPBICG, ierr) CASE DEFAULT PRINT *, "PETSC SOLVER: SELECTED DEFAULT" call KSPSETTYPE(ksp, KSPCG, ierr) END


SELECT ENDIF call KSPGETPC(ksp, pc, ierr) call KSPSETTOLERANCES(ksp,tolerance, & PETSC_DEFAULT_DOUBLE_PRECISION, & PETSC_DEFAULT_DOUBLE_PRECISION, & maxiter,ierr) SELECT CASE (PETSC_PRECOND_TYPE) CASE ('BLOCKJACOBI') PRINT *, "PETSC PRECOND: SELECTED BJACOBI" call PCSETTYPE(pc, PCBJACOBI, ierr) call KSPSETUP (ksp, ierr) call PCBJACOBIGETSUBKSP (pc,PETSC_NULL_INTEGER, & PETSC_NULL_INTEGER, & subksp,ierr) call KSPGETPC (subksp, subpc, ierr) call PCSETTYPE (subpc, PCILU, ierr) ! call PCFactorSetLevels(subpc,0,ierr) CASE ('JACOBI') PRINT *, "PETSC PRECOND: SELECTED JACOBI" call PCSETTYPE(pc, PCJACOBI, ierr) CASE ('ILU') PRINT *, "PETSC PRECOND: SELECTED ILU" call PCSETTYPE(pc, PCILU, ierr) CASE ('GAMG') PRINT *, "PETSC PRECOND: SELECTED GAMG" call PCSETTYPE(pc, PCGAMG, ierr) call PCGAMGSETCOARSEEQLIM(pc,10000,ierr) call PCGAMGSETSYMGRAPH(pc, PETSC_TRUE,ierr) call PCGAMGSETNSMOOTHS(pc, 0,ierr) call PCGAMGSETTHRESHOLD(pc, .001,ierr) ! call PCGAMGSetReuseProl(pc,PETSC_FALSE,ierr) call KSPSETUP (ksp, ierr) CASE ('MUMPS') PRINT *, "PETSC PRECOND: SELECTED MUMPS" call PCSETTYPE(pc,PCLU,ierr) call PCFACTORSETMATSOLVERPACKAGE(pc,MATSOLVERMUMPS,ierr) call PCFACTORSETUPMATSOLVERPACKAGE(pc,ierr) call PCFACTORGETMATRIX(pc,mumpsFac,ierr) call MATMUMPSSETICNTL(mumpsfac,7,2,ierr) call MATMUMPSSETICNTL(mumpsfac,24,1,ierr) ! call MatMumpsSetCntl(mumpsfac,3,2.e-6,ierr) call KSPSETUP (ksp, ierr) CASE DEFAULT PRINT *, "PETSC PRECOND: SELECTED DEFAULT" call PCSETTYPE(pc, PCBJACOBI, ierr) END


SELECT CALL TIMER_STOP ('STREAMICE_PETSC_SETUP',myThid) #if (defined (ALLOW_OPENAD) defined (ALLOW_STREAMICE_OAD_FP)) #ifdef ALLOW_PETSC endif #endif #endif CALL TIMER_START ('STREAMICE_PETSC_SOLVE',myThid) call KSPSOLVE(ksp, rhs, solution, ierr) CALL TIMER_STOP ('STREAMICE_PETSC_SOLVE',myThid) call KSPGETITERATIONNUMBER(ksp,iters,ierr) call VECGETVALUES(solution,local_dofs,indices, & solution_values,ierr) DO bj = myByLo(myThid), myByHi(myThid) DO bi = myBxLo(myThid), myBxHi(myThid) DO j=1,sNy DO i=1,sNx dof_index = INT(streamice_petsc_dofs_u(i,j,bi,bj)) & - local_offset if (dof_index.ge.0) THEN cg_Uin(i,j,bi,bj) = solution_values(dof_index+1) endif dof_index = INT(streamice_petsc_dofs_v(i,j,bi,bj)) & - local_offset if (dof_index.ge.0) THEN cg_Vin(i,j,bi,bj) = solution_values(dof_index+1) endif ENDDO ENDDO ENDDO ENDDO #if (defined (ALLOW_OPENAD) defined (ALLOW_STREAMICE_OAD_FP)) if (streamice_need2destroymat) then #endif call KSPDESTROY (ksp, ierr) call MATDESTROY (matrix, ierr) #if (defined (ALLOW_OPENAD) defined (ALLOW_STREAMICE_OAD_FP)) endif #endif call VECDESTROY (rhs, ierr) call VECDESTROY (solution, ierr) ! call PetscFinalize(ierr) #endif #endif RETURN END