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