source: trunk/NEMO/OPA_SRC/SOL/solpcg.F90 @ 787

Last change on this file since 787 was 787, checked in by rblod, 13 years ago

Improve PCG algortithm in MPI case, see ticket #47

  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 9.3 KB
Line 
1MODULE solpcg
2   !!======================================================================
3   !!                     ***  MODULE  solfet
4   !! Ocean solver :  preconditionned conjugate gradient solver
5   !!=====================================================================
6
7   !!----------------------------------------------------------------------
8   !!   sol_pcg    : preconditionned conjugate gradient solver
9   !!----------------------------------------------------------------------
10   !! * Modules used
11   USE oce             ! ocean dynamics and tracers variables
12   USE dom_oce         ! ocean space and time domain variables
13   USE sol_oce         ! ocean solver variables
14   USE lib_mpp         ! distributed memory computing
15   USE lbclnk          ! ocean lateral boundary conditions (or mpp link)
16   USE in_out_manager  ! I/O manager
17
18   IMPLICIT NONE
19   PRIVATE
20
21   !! * Routine accessibility
22   PUBLIC sol_pcg              ! ???
23
24   !! * Substitutions
25#  include "vectopt_loop_substitute.h90"
26   !!----------------------------------------------------------------------
27   !!----------------------------------------------------------------------
28   !!  OPA 9.0 , LOCEAN-IPSL (2005)
29   !! $Header$
30   !! This software is governed by the CeCILL licence see modipsl/doc/NEMO_CeCILL.txt
31   !!----------------------------------------------------------------------
32CONTAINS
33
34   SUBROUTINE sol_pcg( kindic )
35      !!----------------------------------------------------------------------
36      !!                  ***  ROUTINE sol_pcg  ***
37      !!                   
38      !! ** Purpose :   Solve the ellipic equation for the barotropic stream
39      !!      function system (lk_dynspg_rl=T) or the transport divergence
40      !!      system (lk_dynspg_flt=T) using a diagonal preconditionned
41      !!      conjugate gradient method.
42      !!      In the former case, the barotropic stream function trend has a
43      !!      zero boundary condition along all coastlines (i.e. continent
44      !!      as well as islands) while in the latter the boundary condition
45      !!      specification is not required.
46      !!
47      !! ** Method  :   Diagonal preconditionned conjugate gradient method.
48      !!      the algorithm is multitasked. (case of 5 points matrix)
49      !!      define              pa  = q^-1 * a
50      !!                        pgcb  = q^-1 * gcb
51      !!                 < . ; . >_q  = ( . )^t q ( . )
52      !!      where q is the preconditioning matrix = diagonal matrix of the
53      !!                                              diagonal elements of a
54      !!      Initialization  :
55      !!         x(o) = gcx
56      !!         r(o) = d(o) = pgcb - pa.x(o)
57      !!         rr(o)= < r(o) , r(o) >_q
58      !!      Iteration 1     :
59      !!         standard PCG algorithm
60      !!      Iteration n > 1 :
61      !!         s(n)   = pa.r(n)
62      !!         gam(n) = < r(n) , r(n) >_q
63      !!         del(n) = < r(n) , s(n) >_q
64      !!         bet(n) = gam(n) / gam(n-1)
65      !!         d(n)   = r(n) + bet(n) d(n-1)
66      !!         z(n)   = s(n) + bet(n) z(n-1)
67      !!         sig(n) = del(n) - bet(n)*bet(n)*sig(n-1)
68      !!         alp(n) = gam(n) / sig(n)
69      !!         x(n+1) = x(n) + alp(n) d(n)
70      !!         r(n+1) = r(n) - alp(n) z(n)
71      !!      Convergence test :
72      !!         rr(n+1) / < gcb , gcb >_q   =< epsr
73      !!
74      !! ** Action : - niter  : solver number of iteration done
75      !!             - res    : solver residu reached
76      !!             - gcx()  : solution of the elliptic system
77      !!
78      !! References :
79      !!      Madec et al. 1988, Ocean Modelling, issue 78, 1-6.
80      !!      D Azevedo et al. 1993, Computer Science Technical Report, Tennessee U.
81      !!
82      !! History :
83      !!        !  90-10  (G. Madec)  Original code
84      !!        !  91-11  (G. Madec)
85      !!        !  93-04  (M. Guyon)  loops and suppress pointers
86      !!        !  95-09  (M. Imbard, J. Escobar)  mpp exchange
87      !!        !  96-05  (G. Madec)  merge sor and pcg formulations
88      !!        !  96-11  (A. Weaver)  correction to preconditioning
89      !!   8.5  !  02-08  (G. Madec)  F90: Free form
90      !!        !  08-01  (R. Benshila) mpp optimization
91      !!----------------------------------------------------------------------
92      !! * Arguments
93      INTEGER, INTENT( inout ) ::   kindic   ! solver indicator, < 0 if the conver-
94      !                                      ! gence is not reached: the model is
95      !                                      ! stopped in step
96      !                                      ! set to zero before the call of solpcg
97
98      !! * Local declarations
99      INTEGER ::   ji, jj, jn                ! dummy loop indices
100      REAL(wp) ::  zgcad                     ! temporary scalars
101      REAL(wp), DIMENSION(2) :: zsum
102      REAL(wp), DIMENSION(jpi,jpj) :: zgcr
103      !!----------------------------------------------------------------------
104
105      ! Initialization of the algorithm with standard PCG
106      ! -------------------------------------------------
107
108      CALL lbc_lnk( gcx, c_solver_pt, 1. )   ! lateral boundary condition
109
110      ! gcr   = gcb-a.gcx
111      ! gcdes = gcr
112
113      DO jj = 2, jpjm1
114         DO ji = fs_2, fs_jpim1   ! vector opt.
115            zgcad = bmask(ji,jj) * ( gcb(ji,jj  ) -                gcx(ji  ,jj  )   &
116               &                                  - gcp(ji,jj,1) * gcx(ji  ,jj-1)   &
117               &                                  - gcp(ji,jj,2) * gcx(ji-1,jj  )   &
118               &                                  - gcp(ji,jj,3) * gcx(ji+1,jj  )   &
119               &                                  - gcp(ji,jj,4) * gcx(ji  ,jj+1)   )
120            gcr  (ji,jj) = zgcad
121            gcdes(ji,jj) = zgcad
122         END DO
123      END DO
124
125      ! rnorme = (gcr,gcr)
126      rnorme = SUM(  gcr(:,:) * gcdmat(:,:) * gcr(:,:)  )
127      IF( lk_mpp )   CALL mpp_sum( rnorme )   ! sum over the global domain
128
129      CALL lbc_lnk( gcdes, c_solver_pt, 1. )   ! lateral boundary condition
130
131      ! gccd = matrix . gcdes
132      DO jj = 2, jpjm1
133         DO ji = fs_2, fs_jpim1   ! vector opt.
134            gccd(ji,jj) = bmask(ji,jj)*( gcdes(ji,jj)   &
135               &        +gcp(ji,jj,1)*gcdes(ji,jj-1)+gcp(ji,jj,2)*gcdes(ji-1,jj)   &
136               &        +gcp(ji,jj,4)*gcdes(ji,jj+1)+gcp(ji,jj,3)*gcdes(ji+1,jj)   )
137         END DO
138      END DO 
139
140      ! alph = (gcr,gcr)/(gcdes,gccd)
141      radd = SUM(  gcdes(:,:) * gcdmat(:,:) * gccd(:,:)  )
142      IF( lk_mpp )   CALL mpp_sum( radd )   ! sum over the global domain
143      alph = rnorme /radd
144
145      ! gcx = gcx + alph * gcdes
146      ! gcr = gcr - alph * gccd
147      DO jj = 2, jpjm1
148         DO ji = fs_2, fs_jpim1   ! vector opt.
149            gcx(ji,jj) = bmask(ji,jj) * ( gcx(ji,jj) + alph * gcdes(ji,jj) )
150            gcr(ji,jj) = bmask(ji,jj) * ( gcr(ji,jj) - alph * gccd (ji,jj) )
151         END DO
152      END DO
153
154      ! Algorithm wtih Eijkhout rearrangement
155      ! -------------------------------------
156       
157      !                                                !================
158      DO jn = 1, nmax                                  ! Iterative loop
159         !                                             !================
160
161         CALL lbc_lnk( gcr, c_solver_pt, 1. )   ! lateral boundary condition
162
163         ! zgcr = matrix . gcr
164         DO jj = 2, jpjm1
165            DO ji = fs_2, fs_jpim1   ! vector opt.
166               zgcr(ji,jj) = bmask(ji,jj)*( gcr(ji,jj)   &
167                  &        +gcp(ji,jj,1)*gcr(ji,jj-1)+gcp(ji,jj,2)*gcr(ji-1,jj)   &
168                  &        +gcp(ji,jj,4)*gcr(ji,jj+1)+gcp(ji,jj,3)*gcr(ji+1,jj)   )
169            END DO
170         END DO
171 
172         ! rnorme = (gcr,gcr)
173         rr = rnorme
174         zsum(1) = SUM(  gcr(:,:) * gcdmat(:,:) * gcr(:,:)  )
175
176         ! zgcad = (zgcr,gcr)
177         zsum(2) = SUM( gcr(2:jpim1,2:jpjm1) * gcdmat(2:jpim1,2:jpjm1) * zgcr(2:jpim1,2:jpjm1) )
178
179         IF( lk_mpp )   CALL mpp_sum( zsum, 2 )   ! sum over the global domain
180         rnorme = zsum(1) 
181         zgcad  = zsum(2)
182
183         ! test of convergence
184         IF( rnorme < epsr .OR. jn == nmax ) THEN
185            res = SQRT( rnorme )
186            niter = jn
187            ncut = 999
188         ENDIF
189
190         ! beta = (rk+1,rk+1)/(rk,rk)
191         beta = rnorme / rr
192         radd = zgcad - beta*beta*radd
193         alph = rnorme / radd
194
195         ! gcx = gcx + alph * gcdes
196         ! gcr = gcr - alph * gccd
197         DO jj = 2, jpjm1
198            DO ji = fs_2, fs_jpim1   ! vector opt.
199               gcdes(ji,jj) = gcr (ji,jj) + beta * gcdes(ji,jj) 
200               gccd (ji,jj) = zgcr(ji,jj) + beta * gccd (ji,jj) 
201               gcx  (ji,jj) = gcx (ji,jj) + alph * gcdes(ji,jj) 
202               gcr  (ji,jj) = gcr (ji,jj) - alph * gccd (ji,jj) 
203            END DO
204         END DO
205       
206         ! indicator of non-convergence or explosion
207         IF( jn == nmax .OR. SQRT(epsr)/eps > 1.e+20 ) kindic = -2
208         IF( ncut == 999 ) GOTO 999
209
210         !                                             !================
211      END DO                                           !    End Loop
212      !                                                !================
213     
214999   CONTINUE
215     
216     
217      ! Output in gcx with lateral b.c. applied
218      ! ---------------------------------------
219     
220      CALL lbc_lnk( gcx, c_solver_pt, 1. )
221     
222   END SUBROUTINE sol_pcg
223
224   !!=====================================================================
225END MODULE solpcg
Note: See TracBrowser for help on using the repository browser.