Changeset 787 for trunk/NEMO
- Timestamp:
- 2008-01-11T14:33:42+01:00 (16 years ago)
- File:
-
- 1 edited
-
trunk/NEMO/OPA_SRC/SOL/solpcg.F90 (modified) (5 diffs)
Legend:
- Unmodified
- Added
- Removed
-
trunk/NEMO/OPA_SRC/SOL/solpcg.F90
r719 r787 52 52 !! where q is the preconditioning matrix = diagonal matrix of the 53 53 !! diagonal elements of a 54 !! Initialization :54 !! Initialization : 55 55 !! x(o) = gcx 56 56 !! r(o) = d(o) = pgcb - pa.x(o) 57 57 !! rr(o)= < r(o) , r(o) >_q 58 !! Iteration n : 59 !! z(n) = pa.d(n) 60 !! alp(n) = rr(n) / < z(n) , d(n) >_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) 61 69 !! x(n+1) = x(n) + alp(n) d(n) 62 70 !! r(n+1) = r(n) - alp(n) z(n) 63 !! rr(n+1)= < r(n+1) , r(n+1) >_q64 !! bet(n) = rr(n+1) / rr(n)65 !! r(n+1) = r(n+1) + bet(n+1) d(n)66 71 !! Convergence test : 67 72 !! rr(n+1) / < gcb , gcb >_q =< epsr … … 73 78 !! References : 74 79 !! Madec et al. 1988, Ocean Modelling, issue 78, 1-6. 80 !! D Azevedo et al. 1993, Computer Science Technical Report, Tennessee U. 75 81 !! 76 82 !! History : … … 82 88 !! ! 96-11 (A. Weaver) correction to preconditioning 83 89 !! 8.5 ! 02-08 (G. Madec) F90: Free form 90 !! ! 08-01 (R. Benshila) mpp optimization 84 91 !!---------------------------------------------------------------------- 85 92 !! * Arguments … … 91 98 !! * Local declarations 92 99 INTEGER :: ji, jj, jn ! dummy loop indices 93 REAL(wp) :: zgcad ! temporary scalars 100 REAL(wp) :: zgcad ! temporary scalars 101 REAL(wp), DIMENSION(2) :: zsum 102 REAL(wp), DIMENSION(jpi,jpj) :: zgcr 94 103 !!---------------------------------------------------------------------- 95 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 96 157 ! !================ 97 158 DO jn = 1, nmax ! Iterative loop 98 159 ! !================ 99 160 100 IF( jn == 1 ) THEN ! Initialization of the algorithm 101 102 CALL lbc_lnk( gcx, c_solver_pt, 1. ) ! lateral boundary condition 103 104 ! gcr = gcb-a.gcx 105 ! gcdes = gsr 106 DO jj = 2, jpjm1 107 DO ji = fs_2, fs_jpim1 ! vector opt. 108 zgcad = bmask(ji,jj) * ( gcb(ji,jj ) - gcx(ji ,jj ) & 109 & - gcp(ji,jj,1) * gcx(ji ,jj-1) & 110 & - gcp(ji,jj,2) * gcx(ji-1,jj ) & 111 & - gcp(ji,jj,3) * gcx(ji+1,jj ) & 112 & - gcp(ji,jj,4) * gcx(ji ,jj+1) ) 113 gcr (ji,jj) = zgcad 114 gcdes(ji,jj) = zgcad 115 END DO 116 END DO 117 118 rnorme = SUM( gcr(:,:) * gcdmat(:,:) * gcr(:,:) ) 119 IF( lk_mpp ) CALL mpp_sum( rnorme ) ! sum over the global domain 120 rr = rnorme 121 122 ENDIF 123 124 ! ! Algorithm 125 126 CALL lbc_lnk( gcdes, c_solver_pt, 1. ) ! lateral boundary condition 127 128 ! ... gccd = matrix . gcdes 161 CALL lbc_lnk( gcr, c_solver_pt, 1. ) ! lateral boundary condition 162 163 ! zgcr = matrix . gcr 129 164 DO jj = 2, jpjm1 130 165 DO ji = fs_2, fs_jpim1 ! vector opt. 131 gccd(ji,jj) = bmask(ji,jj)*( gcdes(ji,jj) &132 & +gcp(ji,jj,1)*gc des(ji,jj-1)+gcp(ji,jj,2)*gcdes(ji-1,jj) &133 & +gcp(ji,jj,4)*gc des(ji,jj+1)+gcp(ji,jj,3)*gcdes(ji+1,jj) )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) ) 134 169 END DO 135 170 END DO 136 171 137 ! alph = (gcr,gcr)/(gcdes,gccd)138 radd = SUM( gcdes(:,:) * gcdmat(:,:) * gccd(:,:) )139 IF( lk_mpp ) CALL mpp_sum( radd ) ! sum over the global domain140 alph = rr / radd141 142 ! gcx = gcx + alph * gcdes143 ! gcr = gcr - alph * gccd144 DO jj = 2, jpjm1145 DO ji = fs_2, fs_jpim1 ! vector opt.146 gcx(ji,jj) = bmask(ji,jj) * ( gcx(ji,jj) + alph * gcdes(ji,jj) )147 gcr(ji,jj) = bmask(ji,jj) * ( gcr(ji,jj) - alph * gccd (ji,jj) )148 END DO149 END DO150 151 172 ! rnorme = (gcr,gcr) 152 rnorme = SUM( gcr(:,:) * gcdmat(:,:) * gcr(:,:) ) 153 IF( lk_mpp ) CALL mpp_sum( rnorme ) ! sum over the global domain 154 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 155 183 ! test of convergence 156 184 IF( rnorme < epsr .OR. jn == nmax ) THEN … … 159 187 ncut = 999 160 188 ENDIF 161 189 162 190 ! beta = (rk+1,rk+1)/(rk,rk) 163 191 beta = rnorme / rr 164 rr = rnorme 165 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 166 206 ! indicator of non-convergence or explosion 167 207 IF( jn == nmax .OR. SQRT(epsr)/eps > 1.e+20 ) kindic = -2 168 208 IF( ncut == 999 ) GOTO 999 169 209 170 ! gcdes = gcr + beta * gcdes171 DO jj = 2, jpjm1172 DO ji = fs_2, fs_jpim1 ! vector opt.173 gcdes(ji,jj) = bmask(ji,jj)*( gcr(ji,jj) + beta * gcdes(ji,jj) )174 END DO175 END DO176 177 210 ! !================ 178 211 END DO ! End Loop
Note: See TracChangeset
for help on using the changeset viewer.
