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 @ 3837

Last change on this file since 3837 was 3837, checked in by trackstand2, 11 years ago

Merge of finiss
  • Property svn:keywords set to Id
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   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!      USE arpdebugging, ONLY: dump_array
92      !!
93      INTEGER, INTENT(inout) ::   kindic   ! solver indicator, < 0 if the conver-
94      !                                    ! gence is not reached: the model is stopped in step
95      !                                    ! set to zero before the call of solpcg
96      !!
97      INTEGER  ::   ji, jj, jn   ! dummy loop indices
98      REAL(wp) ::   zgcad        ! temporary scalars
99      REAL(wp), DIMENSION(2) ::   zsum
100      INTEGER, SAVE :: istep = 0 ! ARPDBG
101      !!----------------------------------------------------------------------
102     
103      IF( wrk_in_use(2, 1) )THEN
104         CALL ctl_stop('sol_pcg: requested workspace array is unavailable')   ;   RETURN
105      ENDIF
106
107      ! Initialization of the algorithm with standard PCG
108      ! -------------------------------------------------
109      zgcr = 0._wp
110      gcr  = 0._wp
111!      CALL dump_array(istep, 'gcx_pre_lbc', gcx, withHalos=.TRUE.)
112
113      CALL lbc_lnk( gcx, c_solver_pt, 1. )   ! lateral boundary condition
114
115      istep = istep + 1
116!      CALL dump_array(istep, 'gcx', gcx, withHalos=.TRUE.)
117!      CALL dump_array(istep, 'gcp', gcp(:,:,1), withHalos=.TRUE.)
118!      CALL dump_array(istep, 'gcb', gcb, withHalos=.TRUE.)
119!      CALL dump_array(istep, 'ua', ua(:,:,1), withHalos=.TRUE.)
120
121      ! gcr   = gcb-a.gcx
122      ! gcdes = gcr
123      DO jj = 2, jpjm1
124         DO ji = fs_2, fs_jpim1   ! vector opt.
125            zgcad = bmask(ji,jj) * ( gcb(ji,jj  ) -                gcx(ji  ,jj  )   &
126               &                                  - gcp(ji,jj,1) * gcx(ji  ,jj-1)   &
127               &                                  - gcp(ji,jj,2) * gcx(ji-1,jj  )   &
128               &                                  - gcp(ji,jj,3) * gcx(ji+1,jj  )   &
129               &                                  - gcp(ji,jj,4) * gcx(ji  ,jj+1)   )
130            gcr  (ji,jj) = zgcad
131            gcdes(ji,jj) = zgcad
132         END DO
133      END DO
134
135      ! rnorme = (gcr,gcr)
136      rnorme = glob_sum(  gcr(:,:) * gcdmat(:,:) * gcr(:,:)  )
137
138      CALL lbc_lnk( gcdes, c_solver_pt, 1. )   ! lateral boundary condition
139
140      ! gccd = matrix . gcdes
141      DO jj = 2, jpjm1
142         DO ji = fs_2, fs_jpim1   ! vector opt.
143            gccd(ji,jj) = bmask(ji,jj)*( gcdes(ji,jj)   &
144               &        +gcp(ji,jj,1)*gcdes(ji,jj-1)+gcp(ji,jj,2)*gcdes(ji-1,jj)   &
145               &        +gcp(ji,jj,4)*gcdes(ji,jj+1)+gcp(ji,jj,3)*gcdes(ji+1,jj)   )
146         END DO
147      END DO 
148
149      ! alph = (gcr,gcr)/(gcdes,gccd)
150      radd = glob_sum(  gcdes(:,:) * gcdmat(:,:) * gccd(:,:)  )
151      alph = rnorme /radd
152
153      ! gcx = gcx + alph * gcdes
154      ! gcr = gcr - alph * gccd
155      DO jj = 2, jpjm1
156         DO ji = fs_2, fs_jpim1   ! vector opt.
157            gcx(ji,jj) = bmask(ji,jj) * ( gcx(ji,jj) + alph * gcdes(ji,jj) )
158            gcr(ji,jj) = bmask(ji,jj) * ( gcr(ji,jj) - alph * gccd (ji,jj) )
159         END DO
160      END DO
161
162      ! Algorithm wtih Eijkhout rearrangement
163      ! -------------------------------------
164       
165      !                                                !================
166      DO jn = 1, nn_nmax                               ! Iterative loop
167         !                                             !================
168
169         CALL lbc_lnk( gcr, c_solver_pt, 1. )   ! lateral boundary condition
170
171         ! zgcr = matrix . gcr
172         DO jj = 2, jpjm1
173            DO ji = fs_2, fs_jpim1   ! vector opt.
174               zgcr(ji,jj) = bmask(ji,jj)*( gcr(ji,jj)   &
175                  &        +gcp(ji,jj,1)*gcr(ji,jj-1)+gcp(ji,jj,2)*gcr(ji-1,jj)   &
176                  &        +gcp(ji,jj,4)*gcr(ji,jj+1)+gcp(ji,jj,3)*gcr(ji+1,jj)   )
177            END DO
178         END DO
179 
180         ! rnorme = (gcr,gcr)
181         rr = rnorme
182
183         ! zgcad = (zgcr,gcr)
184         zsum(1) = glob_sum(gcr(:,:) * gcdmat(:,:) * gcr(:,:))
185         zsum(2) = glob_sum(gcr(:,:) * gcdmat(:,:) * zgcr(:,:) * bmask(:,:))
186
187         !!RB we should gather the 2 glob_sum
188         rnorme = zsum(1) 
189         zgcad  = zsum(2)
190         ! test of convergence
191         IF( rnorme < epsr .OR. jn == nn_nmax ) THEN
192            res = SQRT( rnorme )
193            niter = jn
194            ncut = 999
195         ENDIF
196
197         ! beta = (rk+1,rk+1)/(rk,rk)
198         beta = rnorme / rr
199         radd = zgcad - beta*beta*radd
200         alph = rnorme / radd
201
202         ! gcx = gcx + alph * gcdes
203         ! gcr = gcr - alph * gccd
204         DO jj = 2, jpjm1
205            DO ji = fs_2, fs_jpim1   ! vector opt.
206               gcdes(ji,jj) = gcr (ji,jj) + beta * gcdes(ji,jj) 
207               gccd (ji,jj) = zgcr(ji,jj) + beta * gccd (ji,jj) 
208               gcx  (ji,jj) = gcx (ji,jj) + alph * gcdes(ji,jj) 
209               gcr  (ji,jj) = gcr (ji,jj) - alph * gccd (ji,jj) 
210            END DO
211         END DO
212       
213         ! indicator of non-convergence or explosion
214         IF( jn == nn_nmax .OR. SQRT(epsr)/eps > 1.e+20 ) kindic = -2
215         IF( ncut == 999 ) GOTO 999
216
217         !                                             !================
218      END DO                                           !    End Loop
219      !                                                !================
220999   CONTINUE
221         
222      CALL lbc_lnk( gcx, c_solver_pt, 1. )      ! Output in gcx with lateral b.c. applied
223      !
224      IF( wrk_not_released(2, 1) )   CALL ctl_stop('sol_pcg: failed to release workspace array')
225      !
226   END SUBROUTINE sol_pcg
227
228   !!=====================================================================
229END MODULE solpcg
Note: See TracBrowser for help on using the repository browser.