/[MITgcm]/MITgcm_contrib/dgoldberg/streamice/streamice_cg_solve.F
ViewVC logotype

Diff of /MITgcm_contrib/dgoldberg/streamice/streamice_cg_solve.F

Parent Directory Parent Directory | Revision Log Revision Log | View Revision Graph Revision Graph | View Patch Patch

revision 1.7 by dgoldberg, Sat Apr 6 17:43:41 2013 UTC revision 1.10 by dgoldberg, Wed Aug 27 19:29:13 2014 UTC
# Line 17  CBOP Line 17  CBOP
17       I                               A_vv,       ! section of matrix that multiplies v and projects on v       I                               A_vv,       ! section of matrix that multiplies v and projects on v
18       I                               tolerance,       I                               tolerance,
19       O                               iters,       O                               iters,
20         I                               maxIter,
21       I                               myThid )       I                               myThid )
22  C     /============================================================\  C     /============================================================\
23  C     | SUBROUTINE                                                 |    C     | SUBROUTINE                                                 |  
# Line 31  C     \================================= Line 32  C     \=================================
32  #include "PARAMS.h"  #include "PARAMS.h"
33  #include "STREAMICE.h"  #include "STREAMICE.h"
34  #include "STREAMICE_CG.h"  #include "STREAMICE_CG.h"
35  #ifdef ALLOW_PETSC  
36  #include "finclude/petsc.h"  
37    
38    !#ifdef ALLOW_PETSC
39    !#include "finclude/petsc.h"
40    ! UNCOMMENT IF V3.0
41  !#include "finclude/petscvec.h"  !#include "finclude/petscvec.h"
42  !#include "finclude/petscmat.h"  !#include "finclude/petscmat.h"
43  !#include "finclude/petscksp.h"  !#include "finclude/petscksp.h"
44  !#include "finclude/petscpc.h"  !#include "finclude/petscpc.h"
45  #endif  !#endif
46  C     === Global variables ===  C     === Global variables ===
47    
48                
# Line 46  C     cg_Uin, cg_Vin - input and output Line 51  C     cg_Uin, cg_Vin - input and output
51  C     cg_Bu, cg_Bv - driving stress  C     cg_Bu, cg_Bv - driving stress
52        INTEGER myThid        INTEGER myThid
53        INTEGER iters        INTEGER iters
54          INTEGER maxIter
55        _RL tolerance        _RL tolerance
56        _RL cg_Uin (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)        _RL cg_Uin (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
57        _RL cg_Vin (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)        _RL cg_Vin (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
# Line 66  C     LOCAL VARIABLES Line 72  C     LOCAL VARIABLES
72        CHARACTER*(MAX_LEN_MBUF) msgBuf        CHARACTER*(MAX_LEN_MBUF) msgBuf
73    
74    
75  #ifdef ALLOW_PETSC  !#ifdef ALLOW_PETSC
76        INTEGER indices(2*(snx*nsx*sny*nsy))  !      INTEGER indices(2*(snx*nsx*sny*nsy))
77        INTEGER n_dofs_cum_sum (0:nPx*nPy-1), idx(1)  !      INTEGER n_dofs_cum_sum (0:nPx*nPy-1), idx(1)
78        _RL rhs_values(2*(snx*nsx*sny*nsy))  !      _RL rhs_values(2*(snx*nsx*sny*nsy))
79        _RL solution_values(2*(snx*nsx*sny*nsy))  !      _RL solution_values(2*(snx*nsx*sny*nsy))
80  !      _RL mat_values (2*Nx*Ny,2*(snx*nsx*sny*nsy))  !      _RL mat_values (2*Nx*Ny,2*(snx*nsx*sny*nsy))
81        _RL mat_values (18,1), mat_val_return(1)  !      _RL mat_values (18,1), mat_val_return(1)
82        INTEGER indices_col(18)  !      INTEGER indices_col(18)
83        INTEGER local_dofs, global_dofs, dof_index, dof_index_col  !      INTEGER local_dofs, global_dofs, dof_index, dof_index_col
84        INTEGER local_offset  !      INTEGER local_offset
85        Mat matrix  !      Mat matrix
86        KSP ksp  !      KSP ksp
87        PC  pc  !      PC  pc
88        Vec rhs  !      Vec rhs
89        Vec solution  !      Vec solution
90        PetscErrorCode ierr  !      PetscErrorCode ierr
91  #ifdef ALLOW_USE_MPI  !#ifdef ALLOW_USE_MPI
92        integer mpiRC, mpiMyWid  !      integer mpiRC, mpiMyWid
93  #endif  !#endif
94  #endif  !#endif
95    
96    
97  #ifdef ALLOW_STREAMICE  #ifdef ALLOW_STREAMICE
98    
       CALL TIMER_START ('STREAMICE_CG_SOLVE',myThid)  
   
 #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 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)  
   
   
       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)  
   
       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  
   
       call KSPGetPC(ksp, pc, ierr)  
       call KSPSetTolerances(ksp,tolerance,  
      &     PETSC_DEFAULT_DOUBLE_PRECISION,      
      &     PETSC_DEFAULT_DOUBLE_PRECISION,  
      &     streamice_max_cg_iter,ierr)  
   
       SELECT CASE (PETSC_PRECOND_TYPE)  
        CASE ('BLOCKJACOBI')  
        PRINT *, "PETSC PRECOND: SELECTED BJACOBI"  
        call PCSetType(pc, PCBJACOBI, 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 DEFAULT  
        PRINT *, "PETSC PRECOND: SELECTED DEFAULT"  
        call PCSetType(pc, PCBJACOBI, ierr)  
       END SELECT  
   
99    
       call KSPSolve(ksp, rhs, solution, ierr)  
       call KSPGetIterationNumber(ksp,iters,ierr)  
100    
101        call VecGetValues(solution,local_dofs,indices,        CALL TIMER_START ('STREAMICE_CG_SOLVE',myThid)
102       &      solution_values,ierr)  #ifndef STREAMICE_SERIAL_TRISOLVE
   
       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  
103    
104           ENDDO  #ifdef ALLOW_PETSC
         ENDDO  
        ENDDO  
       ENDDO  
105    
106          CALL STREAMICE_CG_SOLVE_PETSC(
107         U         cg_Uin,     ! x-velocities
108         U         cg_Vin,     ! y-velocities
109         I         cg_Bu,      ! force in x dir
110         I         cg_Bv,      ! force in y dir
111         I         A_uu,       ! section of matrix that multiplies u and projects on u
112         I         A_uv,       ! section of matrix that multiplies v and projects on u
113         I         A_vu,       ! section of matrix that multiplies u and projects on v
114         I         A_vv,       ! section of matrix that multiplies v and projects on v
115         I         tolerance,
116         I         maxIter,
117         O         iters,
118         I         myThid )
119    
120    
121  #else  #else  /* ALLOW_PETSC */
122    
123    
124        iters = streamice_max_cg_iter        iters = maxIter
125        conv_flag = 0        conv_flag = 0
126    
127    
# Line 442  C     FIND INITIAL RESIDUAL, and initial Line 193  C     FIND INITIAL RESIDUAL, and initial
193         ENDDO         ENDDO
194        ENDDO        ENDDO
195    
         
   
196        DO bj = myByLo(myThid), myByHi(myThid)        DO bj = myByLo(myThid), myByHi(myThid)
197         DO bi = myBxLo(myThid), myBxHi(myThid)         DO bi = myBxLo(myThid), myBxHi(myThid)
198          DO j=1,sNy          DO j=1,sNy
# Line 457  C     FIND INITIAL RESIDUAL, and initial Line 206  C     FIND INITIAL RESIDUAL, and initial
206         ENDDO         ENDDO
207        ENDDO        ENDDO
208    
209    
210        CALL GLOBAL_SUM_TILE_RL( dot_p1_tile, dot_p1, myThid )        CALL GLOBAL_SUM_TILE_RL( dot_p1_tile, dot_p1, myThid )
211        resid_0 = sqrt(dot_p1)        resid_0 = sqrt(dot_p1)
212    
213          DO bj = myByLo(myThid), myByHi(myThid)
214           DO bi = myBxLo(myThid), myBxHi(myThid)
215    
216           WRITE(msgBuf,'(A,I1,I1,E14.7)') 'CONJ GRAD INIT RESID LOCAL, ',
217         &         bi,bj, dot_p1_tile(bi,bj)
218           CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
219         &                     SQUEEZE_RIGHT , 1)
220    
221           enddo
222          enddo
223    
224        WRITE(msgBuf,'(A,E14.7)') 'CONJ GRAD INIT RESID, ',        WRITE(msgBuf,'(A,E14.7)') 'CONJ GRAD INIT RESID, ',
225       &         resid_0       &         resid_0
226         CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,         CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
# Line 515  c  ! initially, b-grid data is valid up Line 276  c  ! initially, b-grid data is valid up
276  !       IF(STREAMICE_construct_matrix) CALL STREAMICE_CG_MAKE_A(myThid)  !       IF(STREAMICE_construct_matrix) CALL STREAMICE_CG_MAKE_A(myThid)
277    
278    
279        do iter = 1, streamice_max_cg_iter        do iter = 1, maxIter
280         if (resid .gt. tolerance*resid_0) then         if (resid .gt. tolerance*resid_0) then
281    
282  c      to avoid using "exit"  c      to avoid using "exit"
# Line 739  c      to avoid using "exit" Line 500  c      to avoid using "exit"
500  c     to avoid using "exit"  c     to avoid using "exit"
501  c     if iters has reached max_iters there is no convergence  c     if iters has reached max_iters there is no convergence
502                
503        IF (iters .lt. streamice_max_cg_iter) THEN        IF (iters .lt. maxIter) THEN
504         conv_flag = 1         conv_flag = 1
505        ENDIF        ENDIF
506    
# Line 759  c     if iters has reached max_iters the Line 520  c     if iters has reached max_iters the
520  !       _EXCH_XY_RL( cg_Uin, myThid )  !       _EXCH_XY_RL( cg_Uin, myThid )
521  !       _EXCH_XY_RL( cg_Vin, myThid )          !       _EXCH_XY_RL( cg_Vin, myThid )        
522    
523    #endif /* ifndef ALLOW_PETSC */
524    
525    #else /* STREAMICE_SERIAL_TRISOLVE */
526    
527          iters = 0
528    
529          CALL STREAMICE_TRIDIAG_SOLVE(
530         U                               cg_Uin,     ! x-velocities
531         U                               cg_Vin,
532         U                               cg_Bu,      ! force in x dir
533         I                               A_uu,       ! section of matrix that multiplies u and projects on u
534         I                               STREAMICE_umask,
535         I                               myThid )
536    
537  #endif  #endif
538    
539        CALL TIMER_STOP ('STREAMICE_CG_SOLVE',myThid)        CALL TIMER_STOP ('STREAMICE_CG_SOLVE',myThid)
Legend:
Removed from v.1.7  
changed lines
  Added in v.1.10
  ViewVC Help
Powered by ViewVC 1.1.22