/[MITgcm]/MITgcm/pkg/streamice/streamice_cg_solve_petsc.F
ViewVC logotype

Annotation of /MITgcm/pkg/streamice/streamice_cg_solve_petsc.F

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


Revision 1.5 - (hide annotations) (download)
Mon Jun 15 14:34:38 2015 UTC (8 years, 11 months ago) by dgoldberg
Branch: MAIN
CVS Tags: checkpoint66g, checkpoint66f, checkpoint66e, checkpoint66d, checkpoint66c, checkpoint66b, checkpoint66a, checkpoint66o, checkpoint66n, checkpoint66m, checkpoint66l, checkpoint66k, checkpoint66j, checkpoint66i, checkpoint66h, checkpoint65z, checkpoint65x, checkpoint65y, checkpoint65r, checkpoint65s, checkpoint65p, checkpoint65q, checkpoint65v, checkpoint65w, checkpoint65t, checkpoint65u, checkpoint65n, checkpoint65o, HEAD
Changes since 1.4: +49 -8 lines
functionality to persist PETSc matrix if appropriate switches are set

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

  ViewVC Help
Powered by ViewVC 1.1.22