/[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.6 by dgoldberg, Thu Jul 26 16:22:58 2012 UTC revision 1.7 by dgoldberg, Sat Apr 6 17:43:41 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    #ifdef ALLOW_PETSC
35    #include "finclude/petsc.h"
36    !#include "finclude/petscvec.h"
37    !#include "finclude/petscmat.h"
38    !#include "finclude/petscksp.h"
39    !#include "finclude/petscpc.h"
40    #endif
41    C     === Global variables ===
42    
43                
44  C     !INPUT/OUTPUT ARGUMENTS  C     !INPUT/OUTPUT ARGUMENTS
# Line 58  C     LOCAL VARIABLES Line 65  C     LOCAL VARIABLES
65        _RL dot_p2_tile (nSx,nSy)        _RL dot_p2_tile (nSx,nSy)
66        CHARACTER*(MAX_LEN_MBUF) msgBuf        CHARACTER*(MAX_LEN_MBUF) msgBuf
67    
68        iters = streamice_max_cg_iter  
69    #ifdef ALLOW_PETSC
70          INTEGER indices(2*(snx*nsx*sny*nsy))
71          INTEGER n_dofs_cum_sum (0:nPx*nPy-1), idx(1)
72          _RL rhs_values(2*(snx*nsx*sny*nsy))
73          _RL solution_values(2*(snx*nsx*sny*nsy))
74    !      _RL mat_values (2*Nx*Ny,2*(snx*nsx*sny*nsy))
75          _RL mat_values (18,1), mat_val_return(1)
76          INTEGER indices_col(18)
77          INTEGER local_dofs, global_dofs, dof_index, dof_index_col
78          INTEGER local_offset
79          Mat matrix
80          KSP ksp
81          PC  pc
82          Vec rhs
83          Vec solution
84          PetscErrorCode ierr
85    #ifdef ALLOW_USE_MPI
86          integer mpiRC, mpiMyWid
87    #endif
88    #endif
89    
90    
91  #ifdef ALLOW_STREAMICE  #ifdef ALLOW_STREAMICE
92    
93          CALL TIMER_START ('STREAMICE_CG_SOLVE',myThid)
94    
95    #ifdef ALLOW_PETSC
96    
97    #ifdef ALLOW_USE_MPI
98    
99    
100          CALL MPI_COMM_RANK( MPI_COMM_WORLD, mpiMyWId, mpiRC )
101          local_dofs = n_dofs_process (mpiMyWid)
102          global_dofs = 0
103          
104          n_dofs_cum_sum(0) = 0
105          DO i=0,nPx*nPy-1
106           global_dofs = global_dofs + n_dofs_process (i)
107           if (i.ge.1) THEN
108             n_dofs_cum_sum(i) = n_dofs_cum_sum(i-1)+
109         &                     n_dofs_process(i-1)
110           endif
111          ENDDO
112          local_offset = n_dofs_cum_sum(mpimywid)
113    
114    #else
115    
116          local_dofs = n_dofs_process (0)
117          global_dofs = local_dofs
118          local_offset = 0
119    
120    #endif
121    
122    !      call petscInitialize(PETSC_NULL_CHARACTER,ierr)
123    
124    !----------------------
125    
126          call VecCreate(PETSC_COMM_WORLD, rhs, ierr)
127          call VecSetSizes(rhs, local_dofs, global_dofs, ierr)
128          call VecSetType(rhs, VECMPI, ierr)
129    
130          call VecCreate(PETSC_COMM_WORLD, solution, ierr)
131          call VecSetSizes(solution, local_dofs, global_dofs, ierr)
132          call VecSetType(solution, VECMPI, ierr)
133    
134          do i=1,local_dofs
135            indices(i) = i-1 + local_offset
136          end do
137          do i=1,2*nSx*nSy*sNx*sNy
138            rhs_values (i) = 0. _d 0
139            solution_values (i) = 0. _d 0
140          enddo
141    
142    ! gather rhs and initial guess values to populate petsc vectors
143    
144          DO bj = myByLo(myThid), myByHi(myThid)
145           DO bi = myBxLo(myThid), myBxHi(myThid)
146            DO j=1,sNy
147             DO i=1,sNx
148    
149              dof_index = INT(streamice_petsc_dofs_u(i,j,bi,bj))
150         &                - local_offset
151    
152              if (dof_index.ge.0) THEN
153              
154               rhs_values(dof_index+1) = cg_Bu(i,j,bi,bj)
155               solution_values(dof_index+1) = cg_Uin(i,j,bi,bj)
156    
157              endif
158    
159    !---------------
160    
161              dof_index = INT(streamice_petsc_dofs_v(i,j,bi,bj))
162         &                - local_offset
163    
164              if (dof_index.ge.0) THEN
165    
166               rhs_values(dof_index+1) = cg_Bv(i,j,bi,bj)
167               solution_values(dof_index+1) = cg_Vin(i,j,bi,bj)
168    
169              endif
170    
171             ENDDO
172            ENDDO
173           ENDDO
174          ENDDO
175    
176    
177          call VecSetValues(rhs, local_dofs, indices, rhs_values,
178         &                  INSERT_VALUES, ierr)
179          call VecAssemblyBegin(rhs, ierr)
180          call VecAssemblyEnd(rhs, ierr)
181    
182    
183          call VecSetValues(solution, local_dofs, indices,
184         &                  solution_values, INSERT_VALUES, ierr)
185          call VecAssemblyBegin(solution, ierr)
186          call VecAssemblyEnd(solution, ierr)
187    
188    
189          call MatCreateAIJ (PETSC_COMM_WORLD,
190         &                      local_dofs, local_dofs,
191         &                      global_dofs, global_dofs,
192         &                      18, PETSC_NULL_INTEGER,
193         &                      18, PETSC_NULL_INTEGER,
194         &                      matrix, ierr)
195    
196    ! populate petsc matrix
197    
198          DO bj = myByLo(myThid), myByHi(myThid)
199           DO bi = myBxLo(myThid), myBxHi(myThid)
200            DO j=1,sNy
201             DO i=1,sNx
202    
203              dof_index = INT(streamice_petsc_dofs_u(i,j,bi,bj))
204    !     &                - local_offset
205    
206              IF (dof_index .ge. 0) THEN
207    
208               DO k=1,18
209                indices_col(k) = 0
210                mat_values(k,1) = 0. _d 0
211               ENDDO
212               k=0
213    
214               DO coly=-1,1
215                DO colx=-1,1
216    
217                 dof_index_col = streamice_petsc_dofs_u(i+colx,j+coly,bi,bj)
218    
219                 if (dof_index_col.ge.0) THEN
220    !               pscal = A_uu(i,j,bi,bj,colx,coly)
221    !               CALL MatSetValue (matrix,dof_index, dof_index_col,
222    !     &              pscal,INSERT_VALUES,ierr)
223                    k=k+1
224                    mat_values (k,1) = A_uu(i,j,bi,bj,colx,coly)
225                    indices_col (k) = dof_index_col
226                 endif
227                  
228                 dof_index_col = streamice_petsc_dofs_v(i+colx,j+coly,bi,bj)
229    
230                 if (dof_index_col.ge.0) THEN
231    !               CALL MatSetValue (matrix,dof_index, dof_index_col,
232    !     &              A_uv(i,j,bi,bj,colx,coly),INSERT_VALUES,ierr)
233                    k=k+1
234                    mat_values (k,1) = A_uv(i,j,bi,bj,colx,coly)
235                    indices_col (k) = dof_index_col
236                 endif
237                
238                ENDDO
239               ENDDO
240    
241               call matSetValues (matrix, 1, dof_index, k, indices_col,
242         &                        mat_values,INSERT_VALUES,ierr)
243    
244    
245              ENDIF
246    
247    ! ----------------------------------------------
248    
249              dof_index = INT(streamice_petsc_dofs_v(i,j,bi,bj))
250    !     &                - local_offset
251    
252              IF (dof_index .ge. 0) THEN
253    
254               DO k=1,18
255                indices_col(k) = 0
256                mat_values(k,1) = 0. _d 0
257               ENDDO
258               k=0
259    
260               DO coly=-1,1
261                DO colx=-1,1
262    
263                 dof_index_col = streamice_petsc_dofs_u(i+colx,j+coly,bi,bj)
264    
265                 if (dof_index_col.ge.0) THEN
266    !               CALL MatSetValue (matrix,dof_index, dof_index_col,
267    !     &              A_vu(i,j,bi,bj,colx,coly),INSERT_VALUES,ierr)
268                    k=k+1
269                    mat_values (k,1) = A_vu(i,j,bi,bj,colx,coly)
270                    indices_col (k) = dof_index_col
271                 endif
272    
273                 dof_index_col = streamice_petsc_dofs_v(i+colx,j+coly,bi,bj)
274    
275                 if (dof_index_col.ge.0) THEN
276    !               CALL MatSetValue (matrix,dof_index, dof_index_col,
277    !     &              A_vv(i,j,bi,bj,colx,coly),INSERT_VALUES,ierr)
278                    k=k+1
279                    mat_values (k,1) = A_vv(i,j,bi,bj,colx,coly)
280                    indices_col (k) = dof_index_col
281                 endif
282    
283                ENDDO
284               ENDDO
285    
286               call matSetValues (matrix, 1, dof_index, k, indices_col,
287         &                        mat_values,INSERT_VALUES,ierr)
288              ENDIF
289    
290             ENDDO
291            ENDDO
292           ENDDO
293          ENDDO
294    
295          call MatAssemblyBegin(matrix,MAT_FINAL_ASSEMBLY,ierr)
296          call MatAssemblyEnd(matrix,MAT_FINAL_ASSEMBLY,ierr)
297    
298    
299          call KSPCreate(PETSC_COMM_WORLD, ksp, ierr)
300          call KSPSetOperators(ksp, matrix, matrix,
301         &                     DIFFERENT_NONZERO_PATTERN, ierr)
302    
303          SELECT CASE (PETSC_SOLVER_TYPE)
304           CASE ('CG')
305           PRINT *, "PETSC SOLVER: SELECTED CG"
306           call KSPSetType(ksp, KSPCG, ierr)
307           CASE ('GMRES')
308           PRINT *, "PETSC SOLVER: SELECTED GMRES"
309           call KSPSetType(ksp, KSPGMRES, ierr)
310           CASE ('BICG')
311           PRINT *, "PETSC SOLVER: SELECTED BICG"
312           call KSPSetType(ksp, KSPBICG, ierr)
313           CASE DEFAULT
314           PRINT *, "PETSC SOLVER: SELECTED DEFAULT"
315           call KSPSetType(ksp, KSPCG, ierr)
316          END SELECT
317    
318          call KSPGetPC(ksp, pc, ierr)
319          call KSPSetTolerances(ksp,tolerance,
320         &     PETSC_DEFAULT_DOUBLE_PRECISION,    
321         &     PETSC_DEFAULT_DOUBLE_PRECISION,
322         &     streamice_max_cg_iter,ierr)
323    
324          SELECT CASE (PETSC_PRECOND_TYPE)
325           CASE ('BLOCKJACOBI')
326           PRINT *, "PETSC PRECOND: SELECTED BJACOBI"
327           call PCSetType(pc, PCBJACOBI, ierr)
328           CASE ('JACOBI')
329           PRINT *, "PETSC PRECOND: SELECTED JACOBI"
330           call PCSetType(pc, PCJACOBI, ierr)
331           CASE ('ILU')
332           PRINT *, "PETSC PRECOND: SELECTED ILU"
333           call PCSetType(pc, PCILU, ierr)
334           CASE DEFAULT
335           PRINT *, "PETSC PRECOND: SELECTED DEFAULT"
336           call PCSetType(pc, PCBJACOBI, ierr)
337          END SELECT
338    
339    
340          call KSPSolve(ksp, rhs, solution, ierr)
341          call KSPGetIterationNumber(ksp,iters,ierr)
342    
343          call VecGetValues(solution,local_dofs,indices,
344         &      solution_values,ierr)
345    
346          DO bj = myByLo(myThid), myByHi(myThid)
347           DO bi = myBxLo(myThid), myBxHi(myThid)
348            DO j=1,sNy
349             DO i=1,sNx
350    
351              dof_index = INT(streamice_petsc_dofs_u(i,j,bi,bj))
352         &                - local_offset
353              if (dof_index.ge.0) THEN
354               cg_Uin(i,j,bi,bj) = solution_values(dof_index+1)
355              endif
356              
357              dof_index = INT(streamice_petsc_dofs_v(i,j,bi,bj))
358         &                - local_offset
359              if (dof_index.ge.0) THEN
360               cg_Vin(i,j,bi,bj) = solution_values(dof_index+1)
361              endif
362    
363             ENDDO
364            ENDDO
365           ENDDO
366          ENDDO
367    
368    
369    
370    #else
371    
372    
373          iters = streamice_max_cg_iter
374        conv_flag = 0        conv_flag = 0
375    
376    
# Line 452  c     if iters has reached max_iters the Line 760  c     if iters has reached max_iters the
760  !       _EXCH_XY_RL( cg_Vin, myThid )          !       _EXCH_XY_RL( cg_Vin, myThid )        
761    
762  #endif  #endif
       RETURN  
       END  
763    
764          CALL TIMER_STOP ('STREAMICE_CG_SOLVE',myThid)
         
   
   
   
765    
766    
767    #endif
768          RETURN
769          END
Legend:
Removed from v.1.6  
changed lines
  Added in v.1.7
  ViewVC Help
Powered by ViewVC 1.1.22