/[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.9 by dgoldberg, Sat Jun 8 22:15:33 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 34  C     \================================= Line 35  C     \=================================
35    
36    
37    
38  #ifdef ALLOW_PETSC  !#ifdef ALLOW_PETSC
39  #include "finclude/petsc.h"  !#include "finclude/petsc.h"
40  ! UNCOMMENT IF V3.0  ! 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 50  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 70  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
# Line 101  C     LOCAL VARIABLES Line 103  C     LOCAL VARIABLES
103    
104  #ifdef ALLOW_PETSC  #ifdef ALLOW_PETSC
105    
106  #ifdef ALLOW_USE_MPI        CALL STREAMICE_CG_SOLVE_PETSC(
107         U         cg_Uin,     ! x-velocities
108         U         cg_Vin,     ! y-velocities
109        CALL MPI_COMM_RANK( MPI_COMM_WORLD, mpiMyWId, mpiRC )       I         cg_Bu,      ! force in x dir
110        local_dofs = n_dofs_process (mpiMyWid)       I         cg_Bv,      ! force in y dir
111        global_dofs = 0       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        n_dofs_cum_sum(0) = 0       I         A_vu,       ! section of matrix that multiplies u and projects on v
114        DO i=0,nPx*nPy-1       I         A_vv,       ! section of matrix that multiplies v and projects on v
115         global_dofs = global_dofs + n_dofs_process (i)       I         tolerance,
116         if (i.ge.1) THEN       I         maxIter,
117           n_dofs_cum_sum(i) = n_dofs_cum_sum(i-1)+       O         iters,
118       &                     n_dofs_process(i-1)       I         myThid )
        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)  
   
 !     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)  
   
       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  
   
   
       call KSPSolve(ksp, rhs, solution, ierr)  
       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  
   
119    
120    
121  #else  /* ALLOW_PETSC */  #else  /* ALLOW_PETSC */
122    
123    
124        iters = streamice_max_cg_iter        iters = maxIter
125        conv_flag = 0        conv_flag = 0
126    
127    
# Line 451  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 466  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 524  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 748  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 772  c     if iters has reached max_iters the Line 524  c     if iters has reached max_iters the
524    
525  #else /* STREAMICE_SERIAL_TRISOLVE */  #else /* STREAMICE_SERIAL_TRISOLVE */
526    
527          iters = 0
528    
529        CALL STREAMICE_TRIDIAG_SOLVE(        CALL STREAMICE_TRIDIAG_SOLVE(
530       U                               cg_Uin,     ! x-velocities       U                               cg_Uin,     ! x-velocities
531       U                               cg_Vin,       U                               cg_Vin,
Legend:
Removed from v.1.9  
changed lines
  Added in v.1.10
  ViewVC Help
Powered by ViewVC 1.1.22