/[MITgcm]/MITgcm_contrib/dgoldberg/code_cg3d_petsc/cg3d_petsc.F
ViewVC logotype

Annotation of /MITgcm_contrib/dgoldberg/code_cg3d_petsc/cg3d_petsc.F

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


Revision 1.1 - (hide annotations) (download)
Fri Jul 1 20:40:30 2016 UTC (9 years ago) by dgoldberg
Branch: MAIN
test files for using cg3d with petsc

1 dgoldberg 1.1 C $Header: /u/gcmpack/MITgcm/pkg/streamice/streamice_cg_solve_petsc.F,v 1.5 2015/06/15 14:34:38 dgoldberg Exp $
2     C $Name: $
3    
4     #include "CPP_OPTIONS.h"
5    
6     C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
7    
8     CBOP
9     SUBROUTINE CG3D_PETSC(
10     U cg3d_x, ! solution vector
11     U cg3d_b, ! rhs
12     I tolerance,
13     U maxIter,
14     I myIter,
15     I myThid )
16    
17     C /============================================================\
18     C | SUBROUTINE |
19     C | o |
20     C |============================================================|
21     C | |
22     C \============================================================/
23     IMPLICIT NONE
24    
25     #include "SIZE.h"
26     #include "EEPARAMS.h"
27     #include "PARAMS.h"
28     #include "CG3D.h"
29    
30     #ifdef ALLOW_PETSC
31     #include "finclude/petsc.h"
32     #include "CG3D_PETSC.h"
33     ! UNCOMMENT IF V3.0
34     !#include "finclude/petscvec.h"
35     !#include "finclude/petscmat.h"
36     !#include "finclude/petscksp.h"
37     !#include "finclude/petscpc.h"
38     #endif
39     C === Global variables ===
40    
41    
42     C !INPUT/OUTPUT ARGUMENTS
43     C cg_Uin, cg_Vin - input and output velocities
44     C cg_Bu, cg_Bv - driving stress
45     INTEGER myThid
46     INTEGER maxiter
47     INTEGER myIter
48     _RL tolerance
49     _RL cg3d_x (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr,nSx,nSy)
50     _RL cg3d_b (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr,nSx,nSy)
51    
52    
53     C LOCAL VARIABLES
54     INTEGER i, j, bi, bj, cg_halo, conv_flag
55     INTEGER iter, col_iter, k
56     _RL dot_p1, dot_p2, alpha_k, beta_k, resid, resid_0
57     _RL dot_p1_tile (nSx,nSy)
58     _RL dot_p2_tile (nSx,nSy)
59     CHARACTER*(MAX_LEN_MBUF) msgBuf
60     ! LOGICAL create_mat, destroy_mat
61    
62    
63     #ifdef ALLOW_PETSC
64     INTEGER indices(Nr*(snx*nsx*sny*nsy))
65     INTEGER cg3d_dofs_cum_sum (0:nPx*nPy*MAX_CG3D_PETSC_CPUINVERT-1),
66     & idx(1)
67     _RL rhs_values(Nr*(snx*nsx*sny*nsy))
68     _RL solution_values(Nr*(snx*nsx*sny*nsy))
69     _RL mat_values (7,1), mat_val_return(1)
70     INTEGER indices_col(7)
71     INTEGER local_dofs, global_dofs, dof_index, dof_index_col
72     INTEGER local_offset
73     INTEGER iters
74     INTEGER color, rank
75     PC subpc
76     KSP subksp(1)
77     Vec rhs
78     Vec solution
79     PetscErrorCode ierr
80     #ifdef ALLOW_USE_MPI
81     integer mpiRC, mpiMyWid
82     #endif
83     #endif
84    
85    
86     ! print *, "GOT HERE CG3D", myByLo(mythid), myByHi(myThid),
87     ! & myBxLo(mythid), myBxHi(myThid)
88     ! RETURN
89    
90     #ifdef ALLOW_PETSC
91    
92     IF ((myIter.eq.1) .OR. (.not.cg3d_petsc_reuse_mat)) THEN
93     cg3d_need2create_mat = .TRUE.
94     ELSE
95     cg3d_need2create_mat = .FALSE.
96     ENDIF
97    
98     IF ((myIter.eq.nTimeSteps) .OR. (.not.cg3d_petsc_reuse_mat)) THEN
99     cg3d_need2destroy_mat = .TRUE.
100     ELSE
101     cg3d_need2destroy_mat = .FALSE.
102     ENDIF
103    
104    
105     #ifdef ALLOW_USE_MPI
106    
107    
108     CALL MPI_COMM_RANK( MPI_COMM_WORLD, mpiMyWId, mpiRC )
109     local_dofs = cg3d_dofs_process (mpiMyWid)
110     global_dofs = 0
111    
112    
113     cg3d_dofs_cum_sum(0) = 0
114    
115     DO i=1,nPx*nPy*MAX_CG3D_PETSC_CPUINVERT-1
116     IF (i.le.cg3d_petsc_cpuInVert) THEN
117     global_dofs = global_dofs + cg3d_dofs_process (i)
118     if (i.ge.1) then
119     cg3d_dofs_cum_sum(i) = cg3d_dofs_cum_sum(i-1)+
120     & cg3d_dofs_process(i-1)
121     endif
122     ENDIF
123     ENDDO
124    
125     local_offset = cg3d_dofs_cum_sum(mpimywid)
126    
127     #else
128    
129     local_dofs = cg3d_dofs_process (0)
130     global_dofs = local_dofs
131     local_offset = 0
132    
133     #endif
134    
135    
136     !----------------------
137    
138    
139     call VecCreate(PETSC_COMM_WORLD, rhs, ierr)
140     call VecSetSizes(rhs, local_dofs, global_dofs, ierr)
141     call VecSetType(rhs, VECMPI, ierr)
142    
143     call VecCreate(PETSC_COMM_WORLD, solution, ierr)
144     call VecSetSizes(solution, local_dofs, global_dofs, ierr)
145     call VecSetType(solution, VECMPI, ierr)
146    
147     do i=1,local_dofs
148     indices(i) = i-1 + local_offset
149     end do
150     do i=1,Nr*nSx*nSy*sNx*sNy
151     rhs_values (i) = 0. _d 0
152     solution_values (i) = 0. _d 0
153     enddo
154    
155     ! gather rhs and initial guess values to populate petsc vectors
156    
157     DO bj = myByLo(myThid), myByHi(myThid)
158     DO bi = myBxLo(myThid), myBxHi(myThid)
159     DO j=1,sNy
160     DO i=1,sNx
161     DO k=1,Nr
162    
163     color = INT(cg3d_petsc_color(i,j,k,bi,bj))
164     rank = cg3d_color_rank (color)
165     dof_index = INT(cg3d_petsc_dofs(i,j,k,bi,bj))
166     & - local_offset
167    
168     if (dof_index.ge.0 .and. rank.eq.mpimywid) THEN
169    
170     rhs_values(dof_index+1) = cg3d_b(i,j,k,bi,bj)
171     solution_values(dof_index+1) = cg3d_x(i,j,k,bi,bj)
172    
173     endif
174    
175     ENDDO
176     ENDDO
177     ENDDO
178     ENDDO
179     ENDDO
180    
181     call VecSetValues(rhs, local_dofs, indices, rhs_values,
182     & INSERT_VALUES, ierr)
183     call VecAssemblyBegin(rhs, ierr)
184     call VecAssemblyEnd(rhs, ierr)
185    
186    
187     call VecSetValues(solution, local_dofs, indices,
188     & solution_values, INSERT_VALUES, ierr)
189     call VecAssemblyBegin(solution, ierr)
190     call VecAssemblyEnd(solution, ierr)
191    
192    
193     if (cg3d_need2create_mat) then
194     CALL TIMER_START ('CG3D_PETSC_MATCREATE',myThid)
195    
196     ! IF USING v3.0 THEN
197     ! call MatCreateMPIAIJ (PETSC_COMM_WORLD,
198     call MatCreateAIJ (PETSC_COMM_WORLD,
199     & local_dofs, local_dofs,
200     & global_dofs, global_dofs,
201     & 7, PETSC_NULL_INTEGER,
202     & 7, PETSC_NULL_INTEGER,
203     & matrix_petsc_cg3d, ierr)
204    
205    
206     ! populate petsc matrix
207    
208     DO bj = myByLo(myThid), myByHi(myThid)
209     DO bi = myBxLo(myThid), myBxHi(myThid)
210     DO j=1,sNy
211     DO i=1,sNx
212     DO k=1,Nr
213    
214     dof_index = INT(cg3d_petsc_dofs(i,j,k,bi,bj))
215     ! & - local_offset
216     color = INT(cg3d_petsc_color(i,j,k,bi,bj))
217     rank = cg3d_color_rank (color)
218    
219     IF (dof_index .ge. 0 .and. rank.eq.mpimywid) THEN
220    
221     DO col_iter=1,7
222     indices_col(col_iter) = 0
223     mat_values(col_iter,1) = 0. _d 0
224     ENDDO
225     col_iter = 0
226    
227     ! ASSUME THE STENCIL
228     ! ac3d (i,j,k) --> 1
229     ! aw3d (i,j,k) --> 2
230     ! as3d (i,j,k) --> 3
231     ! av3d (i,j,k) --> 4
232     ! aw3d (i+1,j,k) --> 5
233     ! as3d (i,j+1,k) --> 6
234     ! av3d (i,j,k+1) --> 7
235    
236     dof_index_col = cg3d_petsc_dofs(i,j,k,bi,bj)
237     if (dof_index_col.ge.0) then
238     col_iter = col_iter + 1
239     mat_values(col_iter,1) = ac3d(i,j,k,bi,bj)
240     indices_col(col_iter) = dof_index_col
241     endif
242    
243     dof_index_col = cg3d_petsc_dofs(i-1,j,k,bi,bj)
244     if (dof_index_col.ge.0) then
245     col_iter = col_iter + 1
246     mat_values(col_iter,1) = aw3d(i,j,k,bi,bj)
247     indices_col(col_iter) = dof_index_col
248     endif
249    
250     dof_index_col = cg3d_petsc_dofs(i,j-1,k,bi,bj)
251     if (dof_index_col.ge.0) then
252     col_iter = col_iter + 1
253     mat_values(col_iter,1) = as3d(i,j,k,bi,bj)
254     indices_col(col_iter) = dof_index_col
255     endif
256    
257     if (k.gt.1) then
258     dof_index_col = cg3d_petsc_dofs(i,j,k-1,bi,bj)
259     if (dof_index_col.ge.0) then
260     col_iter = col_iter + 1
261     mat_values(col_iter,1) = av3d(i,j,k,bi,bj)
262     indices_col(col_iter) = dof_index_col
263     endif
264     endif
265    
266     dof_index_col = cg3d_petsc_dofs(i+1,j,k,bi,bj)
267     if (dof_index_col.ge.0) then
268     col_iter = col_iter + 1
269     mat_values(col_iter,1) = aw3d(i+1,j,k,bi,bj)
270     indices_col(col_iter) = dof_index_col
271     endif
272    
273     dof_index_col = cg3d_petsc_dofs(i,j+1,k,bi,bj)
274     if (dof_index_col.ge.0) then
275     col_iter = col_iter + 1
276     mat_values(col_iter,1) = as3d(i,j+1,k,bi,bj)
277     indices_col(col_iter) = dof_index_col
278     endif
279    
280     if (k.lt.Nr) then
281     dof_index_col = cg3d_petsc_dofs(i,j,k+1,bi,bj)
282     if (dof_index_col.ge.0) then
283     col_iter = col_iter + 1
284     mat_values(col_iter,1) = av3d(i,j,k+1,bi,bj)
285     indices_col(col_iter) = dof_index_col
286     endif
287     endif
288    
289     call matSetValues (matrix_petsc_cg3d, 1, dof_index, col_iter,
290     & indices_col,
291     & mat_values,INSERT_VALUES,ierr)
292    
293    
294     ENDIF
295    
296    
297     ENDDO
298     ENDDO
299     ENDDO
300     ENDDO
301     ENDDO
302    
303     call MatAssemblyBegin(matrix_petsc_cg3d,MAT_FINAL_ASSEMBLY,ierr)
304     call MatAssemblyEnd(matrix_petsc_cg3d,MAT_FINAL_ASSEMBLY,ierr)
305    
306    
307     call KSPCreate(PETSC_COMM_WORLD, ksp_cg3d, ierr)
308     call KSPSetOperators(ksp_cg3d,
309     & matrix_petsc_cg3d, matrix_petsc_cg3d,
310     & DIFFERENT_NONZERO_PATTERN, ierr)
311    
312     IF (CG3D_PETSC_PRECOND_TYPE.eq.'MUMPS') then
313     call KSPSetType(ksp_cg3d,KSPPREONLY,ierr)
314     ELSE
315     SELECT CASE (CG3D_PETSC_SOLVER_TYPE)
316     CASE ('CG')
317     PRINT *, "PETSC SOLVER: SELECTED CG"
318     call KSPSetType(ksp_cg3d, KSPCG, ierr)
319     CASE ('GMRES')
320     PRINT *, "PETSC SOLVER: SELECTED GMRES"
321     call KSPSetType(ksp_cg3d, KSPGMRES, ierr)
322     CASE ('BICG')
323     PRINT *, "PETSC SOLVER: SELECTED BICG"
324     call KSPSetType(ksp_cg3d, KSPBICG, ierr)
325     CASE DEFAULT
326     PRINT *, "PETSC SOLVER: SELECTED DEFAULT"
327     call KSPSetType(ksp_cg3d, KSPCG, ierr)
328     END SELECT
329     ENDIF
330    
331     call KSPGetPC(ksp_cg3d, pc_cg3d, ierr)
332     call KSPSetTolerances(ksp_cg3d,tolerance,
333     & PETSC_DEFAULT_DOUBLE_PRECISION,
334     & PETSC_DEFAULT_DOUBLE_PRECISION,
335     & maxiter,ierr)
336    
337     SELECT CASE (CG3D_PETSC_PRECOND_TYPE)
338     CASE ('BLOCKJACOBI')
339     PRINT *, "PETSC PRECOND: SELECTED BJACOBI"
340     call PCSetType(pc_cg3d, PCBJACOBI, ierr)
341     call kspsetup (ksp_cg3d, ierr)
342     call PCBJacobiGetSubKSP (pc_cg3d,PETSC_NULL_INTEGER,
343     & PETSC_NULL_INTEGER,
344     & subksp,ierr);
345     call KSPGetPC (subksp, subpc, ierr)
346     call PCSetType (subpc, PCILU, ierr)
347     ! call PCFactorSetLevels(subpc,0,ierr)
348     CASE ('JACOBI')
349     PRINT *, "PETSC PRECOND: SELECTED JACOBI"
350     call PCSetType(pc_cg3d, PCJACOBI, ierr)
351     CASE ('ILU')
352     PRINT *, "PETSC PRECOND: SELECTED ILU"
353     call PCSetType(pc_cg3d, PCILU, ierr)
354    
355     CASE ('GAMG')
356     PRINT *, "PETSC PRECOND: SELECTED GAMG"
357     call PCSetType(pc_cg3d, PCGAMG, ierr)
358     call PCGAMGSetCoarseEqLim(pc_cg3d,10000,ierr)
359     call PCGAMGSetSymGraph(pc_cg3d, PETSC_TRUE,ierr)
360     call PCGAMGSetNSmooths(pc_cg3d, 0,ierr)
361     call PCGAMGSetThreshold(pc_cg3d, .001,ierr)
362     ! call PCGAMGSetReuseProl(pc,PETSC_FALSE,ierr)
363     call kspsetup (ksp_cg3d, ierr)
364    
365     CASE ('MUMPS')
366     PRINT *, "PETSC PRECOND: SELECTED MUMPS"
367     call PCSetType(pc_cg3d,PCLU,ierr)
368     call PCFactorSetMatSolverPackage(pc_cg3d,MATSOLVERMUMPS,ierr)
369     call PCFactorSetUpMatSolverPackage(pc_cg3d,ierr)
370     call PCFactorGetMatrix(pc_cg3d,mumpsFac_cg3d,ierr)
371     call MatMumpsSetIcntl(mumpsfac_cg3d,7,2,ierr)
372     call MatMumpsSetIcntl(mumpsfac_cg3d,24,1,ierr)
373     ! call MatMumpsSetCntl(mumpsfac,3,2.e-6,ierr)
374     call kspsetup (ksp_cg3d, ierr)
375    
376     CASE DEFAULT
377     PRINT *, "PETSC PRECOND: SELECTED DEFAULT"
378     call PCSetType(pc_cg3d, PCBJACOBI, ierr)
379     END SELECT
380    
381    
382     CALL TIMER_STOP ('CG3D_PETSC_MATCREATE',myThid)
383     endif ! need3create_mat
384    
385    
386     CALL TIMER_START ('CG3D_PETSC_SOLVE',myThid)
387    
388     call KSPSolve(ksp_cg3d, rhs, solution, ierr)
389    
390     CALL TIMER_STOP ('CG3D_PETSC_SOLVE',myThid)
391    
392     call KSPGetIterationNumber(ksp_cg3d,iters,ierr)
393     print *, "GOT HERE CG3D PETSC SOLVE COMPLETE", iters
394    
395     maxIter = iters
396    
397     call VecGetValues(solution,local_dofs,indices,
398     & solution_values,ierr)
399    
400     DO bj = myByLo(myThid), myByHi(myThid)
401     DO bi = myBxLo(myThid), myBxHi(myThid)
402     DO j=1,sNy
403     DO i=1,sNx
404     DO k=1,Nr
405    
406     dof_index = INT(cg3d_petsc_dofs(i,j,k,bi,bj))
407     & - local_offset
408     color = INT(cg3d_petsc_color(i,j,k,bi,bj))
409     rank = cg3d_color_rank (color)
410     if (dof_index.ge.0.and.rank.eq.mpimywid) THEN
411     cg3d_x(i,j,k,bi,bj) = solution_values(dof_index+1)
412     endif
413    
414     ENDDO
415     ENDDO
416     ENDDO
417     ENDDO
418     ENDDO
419    
420     if (cg3d_need2destroy_mat) then
421     call KSPDestroy (ksp_cg3d, ierr)
422     call MatDestroy (matrix_petsc_cg3d, ierr)
423     endif
424     call VecDestroy (rhs, ierr)
425     call VecDestroy (solution, ierr)
426    
427    
428    
429    
430     #endif
431    
432    
433    
434     RETURN
435     END

  ViewVC Help
Powered by ViewVC 1.1.22