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

source: branches/UKMO/dev_r5518_GO6_package_OMP/NEMOGCM/NEMO/OPA_SRC/SOL/solpcg.F90 @ 9176

Last change on this file since 9176 was 9176, checked in by andmirek, 6 years ago

#2001: OMP directives
File size: 9.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     ! 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), POINTER, DIMENSION(:,:) ::   zgcr
96      REAL(wp),     DIMENSION(jpi, jpj) ::   tmp1, tmp2
97      !!----------------------------------------------------------------------
98      !
99      IF( nn_timing == 1 )  CALL timing_start('sol_pcg')
100      !
101      CALL wrk_alloc( jpi, jpj, zgcr )
102      !
103      ! Initialization of the algorithm with standard PCG
104      ! -------------------------------------------------
105      zgcr = 0._wp
106      gcr  = 0._wp
107
108      CALL lbc_lnk( gcx, c_solver_pt, 1. )   ! lateral boundary condition
109
110      ! gcr   = gcb-a.gcx
111      ! gcdes = gcr
112!$OMP PARALLEL DO PRIVATE(zgcad)
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!$OMP END PARALLEL DO
125
126      ! rnorme = (gcr,gcr)
127      tmp1 = 0.
128!$OMP PARALLEL DO
129      DO jj = 2, jpjm1
130         DO ji = fs_2, fs_jpim1   ! vector opt.
131            tmp1(ji, jj) = gcr(ji, jj) * gcdmat(ji, jj) * gcr(ji, jj)
132         END DO
133      END DO     
134!$OMP END PARALLEL DO
135      rnorme = glob_sum(  tmp1(:,:)  )
136      CALL lbc_lnk( gcdes, c_solver_pt, 1. )   ! lateral boundary condition
137
138      ! gccd = matrix . gcdes
139      gccd = 0.
140!$OMP PARALLEL DO
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!$OMP END PARALLEL DO
149      ! alph = (gcr,gcr)/(gcdes,gccd)
150!$OMP PARALLEL DO
151      DO jj = 1, jpj
152         DO ji = 1, jpi
153            tmp1(ji, jj) = gcdes(ji, jj) * gcdmat(ji, jj) * gccd(ji, jj)
154         END DO
155      END DO     
156!$OMP END PARALLEL DO
157      radd = glob_sum(  tmp1  )
158      alph = rnorme /radd
159
160      ! gcx = gcx + alph * gcdes
161      ! gcr = gcr - alph * gccd
162!$OMP PARALLEL DO
163      DO jj = 2, jpjm1
164         DO ji = fs_2, fs_jpim1   ! vector opt.
165            gcx(ji,jj) = bmask(ji,jj) * ( gcx(ji,jj) + alph * gcdes(ji,jj) )
166            gcr(ji,jj) = bmask(ji,jj) * ( gcr(ji,jj) - alph * gccd (ji,jj) )
167         END DO
168      END DO
169!$OMP END PARALLEL DO
170      ! Algorithm wtih Eijkhout rearrangement
171      ! -------------------------------------
172       
173      !                                                !================
174      DO jn = 1, nn_nmax                               ! Iterative loop
175         !                                             !================
176
177         CALL lbc_lnk( gcr, c_solver_pt, 1. )   ! lateral boundary condition
178
179         ! zgcr = matrix . gcr
180!$OMP PARALLEL DO
181         DO jj = 2, jpjm1
182            DO ji = fs_2, fs_jpim1   ! vector opt.
183               zgcr(ji,jj) = bmask(ji,jj)*( gcr(ji,jj)   &
184                  &        +gcp(ji,jj,1)*gcr(ji,jj-1)+gcp(ji,jj,2)*gcr(ji-1,jj)   &
185                  &        +gcp(ji,jj,4)*gcr(ji,jj+1)+gcp(ji,jj,3)*gcr(ji+1,jj)   )
186            END DO
187         END DO
188 
189         ! rnorme = (gcr,gcr)
190         rr = rnorme
191         
192         ! zgcad = (zgcr,gcr)
193      tmp2 = 0.
194!$OMP PARALLEL
195!$OMP DO
196      DO jj = 1, jpj
197         DO ji = 1, jpi
198            tmp2(ji, jj) = gcr(ji, jj) * gcdmat(ji, jj)
199            tmp1(ji, jj) = tmp2(ji, jj) * gcr(ji, jj)
200         END DO
201      END DO     
202!$OMP END DO
203!$OMP DO
204!DIR$ IVDEP
205      DO jj = 1, jpj
206!DIR$ IVDEP
207         DO ji = 1, jpi
208            tmp2(ji, jj) = tmp2(ji, jj) * zgcr(ji, jj) * bmask(ji, jj)
209         END DO
210      END DO     
211!$OMP END DO
212!$OMP END PARALLEL
213 
214!        zsum(1) = glob_sum(gcr(:,:) * gcdmat(:,:) * gcr(:,:))
215!        zsum(2) = glob_sum(gcr(:,:) * gcdmat(:,:) * zgcr(:,:) * bmask(:,:))
216         zsum = glob_asum_2d(tmp1, tmp2)
217
218         !!RB we should gather the 2 glob_sum
219         rnorme = zsum(1) 
220         zgcad  = zsum(2)
221         ! test of convergence
222         IF( rnorme < epsr .OR. jn == nn_nmax ) THEN
223            res = SQRT( rnorme )
224            niter = jn
225            ncut = 999
226         ENDIF
227
228         ! beta = (rk+1,rk+1)/(rk,rk)
229         beta = rnorme / rr
230         radd = zgcad - beta*beta*radd
231         alph = rnorme / radd
232
233         ! gcx = gcx + alph * gcdes
234         ! gcr = gcr - alph * gccd
235!$OMP PARALLEL DO
236         DO jj = 2, jpjm1
237            DO ji = fs_2, fs_jpim1   ! vector opt.
238               gcdes(ji,jj) = gcr (ji,jj) + beta * gcdes(ji,jj) 
239               gccd (ji,jj) = zgcr(ji,jj) + beta * gccd (ji,jj) 
240               gcx  (ji,jj) = gcx (ji,jj) + alph * gcdes(ji,jj) 
241               gcr  (ji,jj) = gcr (ji,jj) - alph * gccd (ji,jj) 
242            END DO
243         END DO
244       
245         ! indicator of non-convergence or explosion
246         IF( jn == nn_nmax .OR. SQRT(epsr)/eps > 1.e+20 ) kindic = -2
247         IF( ncut == 999 ) GOTO 999
248
249         !                                             !================
250      END DO                                           !    End Loop
251      !                                                !================
252999   CONTINUE
253         
254      CALL lbc_lnk( gcx, c_solver_pt, 1. )      ! Output in gcx with lateral b.c. applied
255      !
256      CALL wrk_dealloc( jpi, jpj, zgcr )
257      !
258      IF( nn_timing == 1 )  CALL timing_stop('sol_pcg')
259      !
260   END SUBROUTINE sol_pcg
261
262   !!=====================================================================
263END MODULE solpcg
Note: See TracBrowser for help on using the repository browser.