New URL for NEMO forge!   http://forge.nemo-ocean.eu

Since March 2022 along with NEMO 4.2 release, the code development moved to a self-hosted GitLab.
This present forge is now archived and remained online for history.
solpcg.F90 in branches/2011/DEV_r2739_STFC_dCSE/NEMOGCM/NEMO/OPA_SRC/SOL – NEMO

source: branches/2011/DEV_r2739_STFC_dCSE/NEMOGCM/NEMO/OPA_SRC/SOL/solpcg.F90 @ 3211

Last change on this file since 3211 was 3211, checked in by spickles2, 13 years ago

Stephen Pickles, 11 Dec 2011

Commit to bring the rest of the DCSE NEMO development branch
in line with the latest development version. This includes
array index re-ordering of all OPA_SRC/.
  • Property svn:keywords set to Id
File size: 8.9 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   USE oce             ! ocean dynamics and tracers variables
11   USE dom_oce         ! ocean space and time domain variables
12   USE sol_oce         ! ocean solver variables
13   USE lib_mpp         ! distributed memory computing
14   USE lbclnk          ! ocean lateral boundary conditions (or mpp link)
15   USE in_out_manager  ! I/O manager
16   USE lib_fortran
17
18   IMPLICIT NONE
19   PRIVATE
20
21   PUBLIC   sol_pcg    !
22
23   !! * Control permutation of array indices
24#  include "oce_ftrans.h90"
25#  include "dom_oce_ftrans.h90"
26
27   !! * Substitutions
28#  include "vectopt_loop_substitute.h90"
29   !!----------------------------------------------------------------------
30   !! NEMO/OPA 3.3 , NEMO Consortium (2010)
31   !! $Id$
32   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
33   !!----------------------------------------------------------------------
34CONTAINS
35
36   SUBROUTINE sol_pcg( kindic )
37      !!----------------------------------------------------------------------
38      !!                  ***  ROUTINE sol_pcg  ***
39      !!                   
40      !! ** Purpose :   Solve the ellipic equation for the transport
41      !!      divergence system  using a diagonal preconditionned
42      !!      conjugate gradient method.
43      !!
44      !! ** Method  :   Diagonal preconditionned conjugate gradient method.
45      !!      the algorithm is multitasked. (case of 5 points matrix)
46      !!      define              pa  = q^-1 * a
47      !!                        pgcb  = q^-1 * gcb
48      !!                 < . ; . >_q  = ( . )^t q ( . )
49      !!      where q is the preconditioning matrix = diagonal matrix of the
50      !!                                              diagonal elements of a
51      !!      Initialization  :
52      !!         x(o) = gcx
53      !!         r(o) = d(o) = pgcb - pa.x(o)
54      !!         rr(o)= < r(o) , r(o) >_q
55      !!      Iteration 1     :
56      !!         standard PCG algorithm
57      !!      Iteration n > 1 :
58      !!         s(n)   = pa.r(n)
59      !!         gam(n) = < r(n) , r(n) >_q
60      !!         del(n) = < r(n) , s(n) >_q
61      !!         bet(n) = gam(n) / gam(n-1)
62      !!         d(n)   = r(n) + bet(n) d(n-1)
63      !!         z(n)   = s(n) + bet(n) z(n-1)
64      !!         sig(n) = del(n) - bet(n)*bet(n)*sig(n-1)
65      !!         alp(n) = gam(n) / sig(n)
66      !!         x(n+1) = x(n) + alp(n) d(n)
67      !!         r(n+1) = r(n) - alp(n) z(n)
68      !!      Convergence test :
69      !!         rr(n+1) / < gcb , gcb >_q   =< epsr
70      !!
71      !! ** Action : - niter  : solver number of iteration done
72      !!             - res    : solver residu reached
73      !!             - gcx()  : solution of the elliptic system
74      !!
75      !! References :
76      !!      Madec et al. 1988, Ocean Modelling, issue 78, 1-6.
77      !!      D Azevedo et al. 1993, Computer Science Technical Report, Tennessee U.
78      !!
79      !! History :
80      !!        !  90-10  (G. Madec)  Original code
81      !!        !  91-11  (G. Madec)
82      !!        !  93-04  (M. Guyon)  loops and suppress pointers
83      !!        !  95-09  (M. Imbard, J. Escobar)  mpp exchange
84      !!        !  96-05  (G. Madec)  merge sor and pcg formulations
85      !!        !  96-11  (A. Weaver)  correction to preconditioning
86      !!   8.5  !  02-08  (G. Madec)  F90: Free form
87      !!        !  08-01  (R. Benshila) mpp optimization
88      !!----------------------------------------------------------------------
89      USE wrk_nemo, ONLY:   wrk_in_use, wrk_not_released
90      USE wrk_nemo, ONLY:   zgcr => wrk_2d_1
91      !!
92      INTEGER, INTENT(inout) ::   kindic   ! solver indicator, < 0 if the conver-
93      !                                    ! gence is not reached: the model is stopped in step
94      !                                    ! set to zero before the call of solpcg
95      !!
96      INTEGER  ::   ji, jj, jn   ! dummy loop indices
97      REAL(wp) ::   zgcad        ! temporary scalars
98      REAL(wp), DIMENSION(2) ::   zsum
99      !!----------------------------------------------------------------------
100     
101      IF( wrk_in_use(2, 1) )THEN
102         CALL ctl_stop('sol_pcg: requested workspace array is unavailable')   ;   RETURN
103      ENDIF
104
105      ! Initialization of the algorithm with standard PCG
106      ! -------------------------------------------------
107      zgcr = 0._wp
108      gcr  = 0._wp
109
110      CALL lbc_lnk( gcx, c_solver_pt, 1. )   ! lateral boundary condition
111
112      ! gcr   = gcb-a.gcx
113      ! gcdes = gcr
114      DO jj = 2, jpjm1
115         DO ji = fs_2, fs_jpim1   ! vector opt.
116            zgcad = bmask(ji,jj) * ( gcb(ji,jj  ) -                gcx(ji  ,jj  )   &
117               &                                  - gcp(ji,jj,1) * gcx(ji  ,jj-1)   &
118               &                                  - gcp(ji,jj,2) * gcx(ji-1,jj  )   &
119               &                                  - gcp(ji,jj,3) * gcx(ji+1,jj  )   &
120               &                                  - gcp(ji,jj,4) * gcx(ji  ,jj+1)   )
121            gcr  (ji,jj) = zgcad
122            gcdes(ji,jj) = zgcad
123         END DO
124      END DO
125
126      ! rnorme = (gcr,gcr)
127      rnorme = glob_sum(  gcr(:,:) * gcdmat(:,:) * gcr(:,:)  )
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 = glob_sum(  gcdes(:,:) * gcdmat(:,:) * gccd(:,:)  )
142      alph = rnorme /radd
143
144      ! gcx = gcx + alph * gcdes
145      ! gcr = gcr - alph * gccd
146      DO jj = 2, jpjm1
147         DO ji = fs_2, fs_jpim1   ! vector opt.
148            gcx(ji,jj) = bmask(ji,jj) * ( gcx(ji,jj) + alph * gcdes(ji,jj) )
149            gcr(ji,jj) = bmask(ji,jj) * ( gcr(ji,jj) - alph * gccd (ji,jj) )
150         END DO
151      END DO
152
153      ! Algorithm wtih Eijkhout rearrangement
154      ! -------------------------------------
155       
156      !                                                !================
157      DO jn = 1, nn_nmax                               ! Iterative loop
158         !                                             !================
159
160         CALL lbc_lnk( gcr, c_solver_pt, 1. )   ! lateral boundary condition
161
162         ! zgcr = matrix . gcr
163         DO jj = 2, jpjm1
164            DO ji = fs_2, fs_jpim1   ! vector opt.
165               zgcr(ji,jj) = bmask(ji,jj)*( gcr(ji,jj)   &
166                  &        +gcp(ji,jj,1)*gcr(ji,jj-1)+gcp(ji,jj,2)*gcr(ji-1,jj)   &
167                  &        +gcp(ji,jj,4)*gcr(ji,jj+1)+gcp(ji,jj,3)*gcr(ji+1,jj)   )
168            END DO
169         END DO
170 
171         ! rnorme = (gcr,gcr)
172         rr = rnorme
173
174         ! zgcad = (zgcr,gcr)
175         zsum(1) = glob_sum(gcr(:,:) * gcdmat(:,:) * gcr(:,:))
176         zsum(2) = glob_sum(gcr(:,:) * gcdmat(:,:) * zgcr(:,:) * bmask(:,:))
177
178         !!RB we should gather the 2 glob_sum
179         rnorme = zsum(1) 
180         zgcad  = zsum(2)
181         ! test of convergence
182         IF( rnorme < epsr .OR. jn == nn_nmax ) THEN
183            res = SQRT( rnorme )
184            niter = jn
185            ncut = 999
186         ENDIF
187
188         ! beta = (rk+1,rk+1)/(rk,rk)
189         beta = rnorme / rr
190         radd = zgcad - beta*beta*radd
191         alph = rnorme / radd
192
193         ! gcx = gcx + alph * gcdes
194         ! gcr = gcr - alph * gccd
195         DO jj = 2, jpjm1
196            DO ji = fs_2, fs_jpim1   ! vector opt.
197               gcdes(ji,jj) = gcr (ji,jj) + beta * gcdes(ji,jj) 
198               gccd (ji,jj) = zgcr(ji,jj) + beta * gccd (ji,jj) 
199               gcx  (ji,jj) = gcx (ji,jj) + alph * gcdes(ji,jj) 
200               gcr  (ji,jj) = gcr (ji,jj) - alph * gccd (ji,jj) 
201            END DO
202         END DO
203       
204         ! indicator of non-convergence or explosion
205         IF( jn == nn_nmax .OR. SQRT(epsr)/eps > 1.e+20 ) kindic = -2
206         IF( ncut == 999 ) GOTO 999
207
208         !                                             !================
209      END DO                                           !    End Loop
210      !                                                !================
211999   CONTINUE
212         
213      CALL lbc_lnk( gcx, c_solver_pt, 1. )      ! Output in gcx with lateral b.c. applied
214      !
215      IF( wrk_not_released(2, 1) )   CALL ctl_stop('sol_pcg: failed to release workspace array')
216      !
217   END SUBROUTINE sol_pcg
218
219   !!=====================================================================
220END MODULE solpcg
Note: See TracBrowser for help on using the repository browser.