/[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.4 by dgoldberg, Thu Jul 19 18:46:56 2012 UTC revision 1.8 by dgoldberg, Tue May 28 22:32:39 2013 UTC
# Line 26  C     | Line 26  C     |
26  C     \============================================================/  C     \============================================================/
27        IMPLICIT NONE        IMPLICIT NONE
28    
 C     === Global variables ===  
29  #include "SIZE.h"  #include "SIZE.h"
30  #include "EEPARAMS.h"  #include "EEPARAMS.h"
31  #include "PARAMS.h"  #include "PARAMS.h"
32  #include "STREAMICE.h"  #include "STREAMICE.h"
33  #include "STREAMICE_CG.h"  #include "STREAMICE_CG.h"
34    
35    
36    
37    #ifdef ALLOW_PETSC
38    #include "finclude/petsc.h"
39    #include "finclude/petscvec.h"
40    #include "finclude/petscmat.h"
41    #include "finclude/petscksp.h"
42    #include "finclude/petscpc.h"
43    #endif
44    C     === Global variables ===
45    
46                
47  C     !INPUT/OUTPUT ARGUMENTS  C     !INPUT/OUTPUT ARGUMENTS
48  C     cg_Uin, cg_Vin - input and output velocities  C     cg_Uin, cg_Vin - input and output velocities
# Line 56  C     LOCAL VARIABLES Line 66  C     LOCAL VARIABLES
66        _RL dot_p1, dot_p2, alpha_k, beta_k, resid, resid_0        _RL dot_p1, dot_p2, alpha_k, beta_k, resid, resid_0
67        _RL dot_p1_tile (nSx,nSy)        _RL dot_p1_tile (nSx,nSy)
68        _RL dot_p2_tile (nSx,nSy)        _RL dot_p2_tile (nSx,nSy)
69          CHARACTER*(MAX_LEN_MBUF) msgBuf
70    
71    
72    #ifdef ALLOW_PETSC
73          INTEGER indices(2*(snx*nsx*sny*nsy))
74          INTEGER n_dofs_cum_sum (0:nPx*nPy-1), idx(1)
75          _RL rhs_values(2*(snx*nsx*sny*nsy))
76          _RL solution_values(2*(snx*nsx*sny*nsy))
77    !      _RL mat_values (2*Nx*Ny,2*(snx*nsx*sny*nsy))
78          _RL mat_values (18,1), mat_val_return(1)
79          INTEGER indices_col(18)
80          INTEGER local_dofs, global_dofs, dof_index, dof_index_col
81          INTEGER local_offset
82          Mat matrix
83          KSP ksp
84          PC  pc
85          Vec rhs
86          Vec solution
87          PetscErrorCode ierr
88    #ifdef ALLOW_USE_MPI
89          integer mpiRC, mpiMyWid
90    #endif
91    #endif
92    
       iters = streamice_max_cg_iter  
93    
94  #ifdef ALLOW_STREAMICE  #ifdef ALLOW_STREAMICE
95    
96    
97    
98          CALL TIMER_START ('STREAMICE_CG_SOLVE',myThid)
99    #ifndef STREAMICE_SERIAL_TRISOLVE
100    
101    #ifdef ALLOW_PETSC
102    
103    #ifdef ALLOW_USE_MPI
104    
105    
106          CALL MPI_COMM_RANK( MPI_COMM_WORLD, mpiMyWId, mpiRC )
107          local_dofs = n_dofs_process (mpiMyWid)
108          global_dofs = 0
109          
110          n_dofs_cum_sum(0) = 0
111          DO i=0,nPx*nPy-1
112           global_dofs = global_dofs + n_dofs_process (i)
113           if (i.ge.1) THEN
114             n_dofs_cum_sum(i) = n_dofs_cum_sum(i-1)+
115         &                     n_dofs_process(i-1)
116           endif
117          ENDDO
118          local_offset = n_dofs_cum_sum(mpimywid)
119    
120    #else
121    
122          local_dofs = n_dofs_process (0)
123          global_dofs = local_dofs
124          local_offset = 0
125    
126    #endif
127    
128    !      call petscInitialize(PETSC_NULL_CHARACTER,ierr)
129    
130    !----------------------
131    
132          call VecCreate(PETSC_COMM_WORLD, rhs, ierr)
133          call VecSetSizes(rhs, local_dofs, global_dofs, ierr)
134          call VecSetType(rhs, VECMPI, ierr)
135    
136          call VecCreate(PETSC_COMM_WORLD, solution, ierr)
137          call VecSetSizes(solution, local_dofs, global_dofs, ierr)
138          call VecSetType(solution, VECMPI, ierr)
139    
140          do i=1,local_dofs
141            indices(i) = i-1 + local_offset
142          end do
143          do i=1,2*nSx*nSy*sNx*sNy
144            rhs_values (i) = 0. _d 0
145            solution_values (i) = 0. _d 0
146          enddo
147    
148    ! gather rhs and initial guess values to populate petsc vectors
149    
150          DO bj = myByLo(myThid), myByHi(myThid)
151           DO bi = myBxLo(myThid), myBxHi(myThid)
152            DO j=1,sNy
153             DO i=1,sNx
154    
155              dof_index = INT(streamice_petsc_dofs_u(i,j,bi,bj))
156         &                - local_offset
157    
158              if (dof_index.ge.0) THEN
159              
160               rhs_values(dof_index+1) = cg_Bu(i,j,bi,bj)
161               solution_values(dof_index+1) = cg_Uin(i,j,bi,bj)
162    
163              endif
164    
165    !---------------
166    
167              dof_index = INT(streamice_petsc_dofs_v(i,j,bi,bj))
168         &                - local_offset
169    
170              if (dof_index.ge.0) THEN
171    
172               rhs_values(dof_index+1) = cg_Bv(i,j,bi,bj)
173               solution_values(dof_index+1) = cg_Vin(i,j,bi,bj)
174    
175              endif
176    
177             ENDDO
178            ENDDO
179           ENDDO
180          ENDDO
181    
182    
183          call VecSetValues(rhs, local_dofs, indices, rhs_values,
184         &                  INSERT_VALUES, ierr)
185          call VecAssemblyBegin(rhs, ierr)
186          call VecAssemblyEnd(rhs, ierr)
187    
188    
189          call VecSetValues(solution, local_dofs, indices,
190         &                  solution_values, INSERT_VALUES, ierr)
191          call VecAssemblyBegin(solution, ierr)
192          call VecAssemblyEnd(solution, ierr)
193    
194    
195          call MatCreateMPIAIJ (PETSC_COMM_WORLD,
196         &                      local_dofs, local_dofs,
197         &                      global_dofs, global_dofs,
198         &                      18, PETSC_NULL_INTEGER,
199         &                      18, PETSC_NULL_INTEGER,
200         &                      matrix, ierr)
201    
202    
203    ! populate petsc matrix
204    
205          DO bj = myByLo(myThid), myByHi(myThid)
206           DO bi = myBxLo(myThid), myBxHi(myThid)
207            DO j=1,sNy
208             DO i=1,sNx
209    
210              dof_index = INT(streamice_petsc_dofs_u(i,j,bi,bj))
211    !     &                - local_offset
212    
213              IF (dof_index .ge. 0) THEN
214    
215               DO k=1,18
216                indices_col(k) = 0
217                mat_values(k,1) = 0. _d 0
218               ENDDO
219               k=0
220    
221               DO coly=-1,1
222                DO colx=-1,1
223    
224                 dof_index_col = streamice_petsc_dofs_u(i+colx,j+coly,bi,bj)
225    
226                 if (dof_index_col.ge.0) THEN
227    !               pscal = A_uu(i,j,bi,bj,colx,coly)
228    !               CALL MatSetValue (matrix,dof_index, dof_index_col,
229    !     &              pscal,INSERT_VALUES,ierr)
230                    k=k+1
231                    mat_values (k,1) = A_uu(i,j,bi,bj,colx,coly)
232                    indices_col (k) = dof_index_col
233                 endif
234                  
235                 dof_index_col = streamice_petsc_dofs_v(i+colx,j+coly,bi,bj)
236    
237                 if (dof_index_col.ge.0) THEN
238    !               CALL MatSetValue (matrix,dof_index, dof_index_col,
239    !     &              A_uv(i,j,bi,bj,colx,coly),INSERT_VALUES,ierr)
240                    k=k+1
241                    mat_values (k,1) = A_uv(i,j,bi,bj,colx,coly)
242                    indices_col (k) = dof_index_col
243                 endif
244                
245                ENDDO
246               ENDDO
247    
248               call matSetValues (matrix, 1, dof_index, k, indices_col,
249         &                        mat_values,INSERT_VALUES,ierr)
250    
251    
252              ENDIF
253    
254    ! ----------------------------------------------
255    
256              dof_index = INT(streamice_petsc_dofs_v(i,j,bi,bj))
257    !     &                - local_offset
258    
259              IF (dof_index .ge. 0) THEN
260    
261               DO k=1,18
262                indices_col(k) = 0
263                mat_values(k,1) = 0. _d 0
264               ENDDO
265               k=0
266    
267               DO coly=-1,1
268                DO colx=-1,1
269    
270                 dof_index_col = streamice_petsc_dofs_u(i+colx,j+coly,bi,bj)
271    
272                 if (dof_index_col.ge.0) THEN
273    !               CALL MatSetValue (matrix,dof_index, dof_index_col,
274    !     &              A_vu(i,j,bi,bj,colx,coly),INSERT_VALUES,ierr)
275                    k=k+1
276                    mat_values (k,1) = A_vu(i,j,bi,bj,colx,coly)
277                    indices_col (k) = dof_index_col
278                 endif
279    
280                 dof_index_col = streamice_petsc_dofs_v(i+colx,j+coly,bi,bj)
281    
282                 if (dof_index_col.ge.0) THEN
283    !               CALL MatSetValue (matrix,dof_index, dof_index_col,
284    !     &              A_vv(i,j,bi,bj,colx,coly),INSERT_VALUES,ierr)
285                    k=k+1
286                    mat_values (k,1) = A_vv(i,j,bi,bj,colx,coly)
287                    indices_col (k) = dof_index_col
288                 endif
289    
290                ENDDO
291               ENDDO
292    
293               call matSetValues (matrix, 1, dof_index, k, indices_col,
294         &                        mat_values,INSERT_VALUES,ierr)
295              ENDIF
296    
297             ENDDO
298            ENDDO
299           ENDDO
300          ENDDO
301    
302          call MatAssemblyBegin(matrix,MAT_FINAL_ASSEMBLY,ierr)
303          call MatAssemblyEnd(matrix,MAT_FINAL_ASSEMBLY,ierr)
304    
305    
306          call KSPCreate(PETSC_COMM_WORLD, ksp, ierr)
307          call KSPSetOperators(ksp, matrix, matrix,
308         &                     DIFFERENT_NONZERO_PATTERN, ierr)
309    
310          SELECT CASE (PETSC_SOLVER_TYPE)
311           CASE ('CG')
312           PRINT *, "PETSC SOLVER: SELECTED CG"
313           call KSPSetType(ksp, KSPCG, ierr)
314           CASE ('GMRES')
315           PRINT *, "PETSC SOLVER: SELECTED GMRES"
316           call KSPSetType(ksp, KSPGMRES, ierr)
317           CASE ('BICG')
318           PRINT *, "PETSC SOLVER: SELECTED BICG"
319           call KSPSetType(ksp, KSPBICG, ierr)
320           CASE DEFAULT
321           PRINT *, "PETSC SOLVER: SELECTED DEFAULT"
322           call KSPSetType(ksp, KSPCG, ierr)
323          END SELECT
324    
325          call KSPGetPC(ksp, pc, ierr)
326          call KSPSetTolerances(ksp,tolerance,
327         &     PETSC_DEFAULT_DOUBLE_PRECISION,    
328         &     PETSC_DEFAULT_DOUBLE_PRECISION,
329         &     streamice_max_cg_iter,ierr)
330    
331          SELECT CASE (PETSC_PRECOND_TYPE)
332           CASE ('BLOCKJACOBI')
333           PRINT *, "PETSC PRECOND: SELECTED BJACOBI"
334           call PCSetType(pc, PCBJACOBI, ierr)
335           CASE ('JACOBI')
336           PRINT *, "PETSC PRECOND: SELECTED JACOBI"
337           call PCSetType(pc, PCJACOBI, ierr)
338           CASE ('ILU')
339           PRINT *, "PETSC PRECOND: SELECTED ILU"
340           call PCSetType(pc, PCILU, ierr)
341           CASE DEFAULT
342           PRINT *, "PETSC PRECOND: SELECTED DEFAULT"
343           call PCSetType(pc, PCBJACOBI, ierr)
344          END SELECT
345    
346    
347          call KSPSolve(ksp, rhs, solution, ierr)
348          call KSPGetIterationNumber(ksp,iters,ierr)
349    
350          call VecGetValues(solution,local_dofs,indices,
351         &      solution_values,ierr)
352    
353          DO bj = myByLo(myThid), myByHi(myThid)
354           DO bi = myBxLo(myThid), myBxHi(myThid)
355            DO j=1,sNy
356             DO i=1,sNx
357    
358              dof_index = INT(streamice_petsc_dofs_u(i,j,bi,bj))
359         &                - local_offset
360              if (dof_index.ge.0) THEN
361               cg_Uin(i,j,bi,bj) = solution_values(dof_index+1)
362              endif
363              
364              dof_index = INT(streamice_petsc_dofs_v(i,j,bi,bj))
365         &                - local_offset
366              if (dof_index.ge.0) THEN
367               cg_Vin(i,j,bi,bj) = solution_values(dof_index+1)
368              endif
369    
370             ENDDO
371            ENDDO
372           ENDDO
373          ENDDO
374    
375    
376    
377    #else  /* ALLOW_PETSC */
378    
379    
380          iters = streamice_max_cg_iter
381        conv_flag = 0        conv_flag = 0
382    
383    
# Line 151  C     FIND INITIAL RESIDUAL, and initial Line 467  C     FIND INITIAL RESIDUAL, and initial
467        CALL GLOBAL_SUM_TILE_RL( dot_p1_tile, dot_p1, myThid )        CALL GLOBAL_SUM_TILE_RL( dot_p1_tile, dot_p1, myThid )
468        resid_0 = sqrt(dot_p1)        resid_0 = sqrt(dot_p1)
469    
470          WRITE(msgBuf,'(A,E14.7)') 'CONJ GRAD INIT RESID, ',
471         &         resid_0
472           CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
473         &                     SQUEEZE_RIGHT , 1)
474    
475  C    CCCCCCCCCCCCCCCCCCCC  C    CCCCCCCCCCCCCCCCCCCC
476    
477        DO bj = myByLo(myThid), myByHi(myThid)        DO bj = myByLo(myThid), myByHi(myThid)
# Line 194  c  !!!!!!!!!!!!!!!!!! Line 515  c  !!!!!!!!!!!!!!!!!!
515        
516  c  ! initially, b-grid data is valid up to 3 halo nodes out -- right? (check for MITgcm!!)  c  ! initially, b-grid data is valid up to 3 halo nodes out -- right? (check for MITgcm!!)
517    
518        print *, "BEGINNING MAIN CG LOOP"        WRITE(msgBuf,'(A)') 'BEGINNING MAIN CG LOOP'
519           CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
520         &                     SQUEEZE_RIGHT , 1)
521    
522  !       IF(STREAMICE_construct_matrix) CALL STREAMICE_CG_MAKE_A(myThid)  !       IF(STREAMICE_construct_matrix) CALL STREAMICE_CG_MAKE_A(myThid)
523    
# Line 443  c     if iters has reached max_iters the Line 766  c     if iters has reached max_iters the
766  !       _EXCH_XY_RL( cg_Uin, myThid )  !       _EXCH_XY_RL( cg_Uin, myThid )
767  !       _EXCH_XY_RL( cg_Vin, myThid )          !       _EXCH_XY_RL( cg_Vin, myThid )        
768    
769  #endif  #endif /* ifndef ALLOW_PETSC */
       RETURN  
       END  
770    
771    #else /* STREAMICE_SERIAL_TRISOLVE */
772    
773                CALL STREAMICE_TRIDIAG_SOLVE(
774         U                               cg_Uin,     ! x-velocities
775         U                               cg_Vin,
776         U                               cg_Bu,      ! force in x dir
777         I                               A_uu,       ! section of matrix that multiplies u and projects on u
778         I                               STREAMICE_umask,
779         I                               myThid )
780    
781    #endif
782    
783          CALL TIMER_STOP ('STREAMICE_CG_SOLVE',myThid)
784    
785    
786    #endif
787          RETURN
788          END
Legend:
Removed from v.1.4  
changed lines
  Added in v.1.8
  ViewVC Help
Powered by ViewVC 1.1.22