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/UKMO/dev_r5518_GO6_package/NEMOGCM/NEMO/OPA_SRC/SOL – NEMO

source: branches/UKMO/dev_r5518_GO6_package/NEMOGCM/NEMO/OPA_SRC/SOL/solpcg.F90 @ 7771

Last change on this file since 7771 was 7771, checked in by frrh, 7 years ago

Apply optimisations to various areas of code replacing the use of
allocated pointers with straightforward direct ALLOCATE and DEALLOCATE
operations.

These optimisations largely have an impact in models featuring MEDUSA,
i.e. those with significant numbers of tracers, although they are
expected to have a small impact in all configurations.

Code developed and tested in NEMO branch branches/UKMO/dev_r5518_optim_GO6_alloc
Tested in stand-alone GO6-GSI8, GO6-GSI8-MEDUSA and UKESM coupled models.
NEMO ticket #1821 documents this change further.
File size: 8.8 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     ! Fortran routines library
17   USE wrk_nemo        ! Memory allocation
18   USE timing          ! Timing
19
20   IMPLICIT NONE
21   PRIVATE
22
23   PUBLIC   sol_pcg    !
24
25   !! * Substitutions
26#  include "vectopt_loop_substitute.h90"
27   !!----------------------------------------------------------------------
28   !! NEMO/OPA 3.3 , NEMO Consortium (2010)
29   !! $Id$
30   !! Software governed by the CeCILL licence     (NEMOGCM/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 transport
39      !!      divergence system  using a diagonal preconditionned
40      !!      conjugate gradient method.
41      !!
42      !! ** Method  :   Diagonal preconditionned conjugate gradient method.
43      !!      the algorithm is multitasked. (case of 5 points matrix)
44      !!      define              pa  = q^-1 * a
45      !!                        pgcb  = q^-1 * gcb
46      !!                 < . ; . >_q  = ( . )^t q ( . )
47      !!      where q is the preconditioning matrix = diagonal matrix of the
48      !!                                              diagonal elements of a
49      !!      Initialization  :
50      !!         x(o) = gcx
51      !!         r(o) = d(o) = pgcb - pa.x(o)
52      !!         rr(o)= < r(o) , r(o) >_q
53      !!      Iteration 1     :
54      !!         standard PCG algorithm
55      !!      Iteration n > 1 :
56      !!         s(n)   = pa.r(n)
57      !!         gam(n) = < r(n) , r(n) >_q
58      !!         del(n) = < r(n) , s(n) >_q
59      !!         bet(n) = gam(n) / gam(n-1)
60      !!         d(n)   = r(n) + bet(n) d(n-1)
61      !!         z(n)   = s(n) + bet(n) z(n-1)
62      !!         sig(n) = del(n) - bet(n)*bet(n)*sig(n-1)
63      !!         alp(n) = gam(n) / sig(n)
64      !!         x(n+1) = x(n) + alp(n) d(n)
65      !!         r(n+1) = r(n) - alp(n) z(n)
66      !!      Convergence test :
67      !!         rr(n+1) / < gcb , gcb >_q   =< epsr
68      !!
69      !! ** Action : - niter  : solver number of iteration done
70      !!             - res    : solver residu reached
71      !!             - gcx()  : solution of the elliptic system
72      !!
73      !! References :
74      !!      Madec et al. 1988, Ocean Modelling, issue 78, 1-6.
75      !!      D Azevedo et al. 1993, Computer Science Technical Report, Tennessee U.
76      !!
77      !! History :
78      !!        !  90-10  (G. Madec)  Original code
79      !!        !  91-11  (G. Madec)
80      !!        !  93-04  (M. Guyon)  loops and suppress pointers
81      !!        !  95-09  (M. Imbard, J. Escobar)  mpp exchange
82      !!        !  96-05  (G. Madec)  merge sor and pcg formulations
83      !!        !  96-11  (A. Weaver)  correction to preconditioning
84      !!   8.5  !  02-08  (G. Madec)  F90: Free form
85      !!        !  08-01  (R. Benshila) mpp optimization
86      !!----------------------------------------------------------------------
87      !!
88      INTEGER, INTENT(inout) ::   kindic   ! solver indicator, < 0 if the conver-
89      !                                    ! gence is not reached: the model is stopped in step
90      !                                    ! set to zero before the call of solpcg
91      !!
92      INTEGER  ::   ji, jj, jn   ! dummy loop indices
93      REAL(wp) ::   zgcad        ! temporary scalars
94      REAL(wp), DIMENSION(2) ::   zsum
95      REAL(wp), ALLOCATABLE, DIMENSION(:,:) ::   zgcr
96      !!----------------------------------------------------------------------
97      !
98      IF( nn_timing == 1 )  CALL timing_start('sol_pcg')
99      !
100      ALLOCATE( zgcr(1:jpi,1:jpj) )
101      !
102      ! Initialization of the algorithm with standard PCG
103      ! -------------------------------------------------
104      zgcr = 0._wp
105      gcr  = 0._wp
106
107      CALL lbc_lnk( gcx, c_solver_pt, 1. )   ! lateral boundary condition
108
109      ! gcr   = gcb-a.gcx
110      ! gcdes = gcr
111      DO jj = 2, jpjm1
112         DO ji = fs_2, fs_jpim1   ! vector opt.
113            zgcad = bmask(ji,jj) * ( gcb(ji,jj  ) -                gcx(ji  ,jj  )   &
114               &                                  - gcp(ji,jj,1) * gcx(ji  ,jj-1)   &
115               &                                  - gcp(ji,jj,2) * gcx(ji-1,jj  )   &
116               &                                  - gcp(ji,jj,3) * gcx(ji+1,jj  )   &
117               &                                  - gcp(ji,jj,4) * gcx(ji  ,jj+1)   )
118            gcr  (ji,jj) = zgcad
119            gcdes(ji,jj) = zgcad
120         END DO
121      END DO
122
123      ! rnorme = (gcr,gcr)
124      rnorme = glob_sum(  gcr(:,:) * gcdmat(:,:) * gcr(:,:)  )
125
126      CALL lbc_lnk( gcdes, c_solver_pt, 1. )   ! lateral boundary condition
127
128      ! gccd = matrix . gcdes
129      DO jj = 2, jpjm1
130         DO ji = fs_2, fs_jpim1   ! vector opt.
131            gccd(ji,jj) = bmask(ji,jj)*( gcdes(ji,jj)   &
132               &        +gcp(ji,jj,1)*gcdes(ji,jj-1)+gcp(ji,jj,2)*gcdes(ji-1,jj)   &
133               &        +gcp(ji,jj,4)*gcdes(ji,jj+1)+gcp(ji,jj,3)*gcdes(ji+1,jj)   )
134         END DO
135      END DO 
136
137      ! alph = (gcr,gcr)/(gcdes,gccd)
138      radd = glob_sum(  gcdes(:,:) * gcdmat(:,:) * gccd(:,:)  )
139      alph = rnorme /radd
140
141      ! gcx = gcx + alph * gcdes
142      ! gcr = gcr - alph * gccd
143      DO jj = 2, jpjm1
144         DO ji = fs_2, fs_jpim1   ! vector opt.
145            gcx(ji,jj) = bmask(ji,jj) * ( gcx(ji,jj) + alph * gcdes(ji,jj) )
146            gcr(ji,jj) = bmask(ji,jj) * ( gcr(ji,jj) - alph * gccd (ji,jj) )
147         END DO
148      END DO
149
150      ! Algorithm wtih Eijkhout rearrangement
151      ! -------------------------------------
152       
153      !                                                !================
154      DO jn = 1, nn_nmax                               ! Iterative loop
155         !                                             !================
156
157         CALL lbc_lnk( gcr, c_solver_pt, 1. )   ! lateral boundary condition
158
159         ! zgcr = matrix . gcr
160         DO jj = 2, jpjm1
161            DO ji = fs_2, fs_jpim1   ! vector opt.
162               zgcr(ji,jj) = bmask(ji,jj)*( gcr(ji,jj)   &
163                  &        +gcp(ji,jj,1)*gcr(ji,jj-1)+gcp(ji,jj,2)*gcr(ji-1,jj)   &
164                  &        +gcp(ji,jj,4)*gcr(ji,jj+1)+gcp(ji,jj,3)*gcr(ji+1,jj)   )
165            END DO
166         END DO
167 
168         ! rnorme = (gcr,gcr)
169         rr = rnorme
170
171         ! zgcad = (zgcr,gcr)
172         zsum(1) = glob_sum(gcr(:,:) * gcdmat(:,:) * gcr(:,:))
173         zsum(2) = glob_sum(gcr(:,:) * gcdmat(:,:) * zgcr(:,:) * bmask(:,:))
174
175         !!RB we should gather the 2 glob_sum
176         rnorme = zsum(1) 
177         zgcad  = zsum(2)
178         ! test of convergence
179         IF( rnorme < epsr .OR. jn == nn_nmax ) THEN
180            res = SQRT( rnorme )
181            niter = jn
182            ncut = 999
183         ENDIF
184
185         ! beta = (rk+1,rk+1)/(rk,rk)
186         beta = rnorme / rr
187         radd = zgcad - beta*beta*radd
188         alph = rnorme / radd
189
190         ! gcx = gcx + alph * gcdes
191         ! gcr = gcr - alph * gccd
192         DO jj = 2, jpjm1
193            DO ji = fs_2, fs_jpim1   ! vector opt.
194               gcdes(ji,jj) = gcr (ji,jj) + beta * gcdes(ji,jj) 
195               gccd (ji,jj) = zgcr(ji,jj) + beta * gccd (ji,jj) 
196               gcx  (ji,jj) = gcx (ji,jj) + alph * gcdes(ji,jj) 
197               gcr  (ji,jj) = gcr (ji,jj) - alph * gccd (ji,jj) 
198            END DO
199         END DO
200       
201         ! indicator of non-convergence or explosion
202         IF( jn == nn_nmax .OR. SQRT(epsr)/eps > 1.e+20 ) kindic = -2
203         IF( ncut == 999 ) GOTO 999
204
205         !                                             !================
206      END DO                                           !    End Loop
207      !                                                !================
208999   CONTINUE
209         
210      CALL lbc_lnk( gcx, c_solver_pt, 1. )      ! Output in gcx with lateral b.c. applied
211      !
212      DEALLOCATE ( zgcr )
213      !
214      IF( nn_timing == 1 )  CALL timing_stop('sol_pcg')
215      !
216   END SUBROUTINE sol_pcg
217
218   !!=====================================================================
219END MODULE solpcg
Note: See TracBrowser for help on using the repository browser.