/[MITgcm]/MITgcm_contrib/heimbach/OpenAD/code_shallow_openad/cg2d.F
ViewVC logotype

Contents of /MITgcm_contrib/heimbach/OpenAD/code_shallow_openad/cg2d.F

Parent Directory Parent Directory | Revision Log Revision Log | View Revision Graph Revision Graph


Revision 1.4 - (show annotations) (download)
Tue May 8 14:04:56 2007 UTC (19 years, 6 months ago) by utke
Branch: MAIN
Changes since 1.3: +22 -19 lines
removing extra copies and synching form main

1 C $Header: /u/gcmpack/MITgcm_contrib/heimbach/OpenAD/code_shallow_openad/cg2d.F,v 1.3 2007/01/11 15:19:59 utke Exp $
2 C $Name: $
3
4 #include "CPP_OPTIONS.h"
5
6 CBOP
7 C !ROUTINE: CG2D
8 C !INTERFACE:
9 SUBROUTINE CG2D(
10 I cg2d_b,
11 U cg2d_x,
12 O firstResidual,
13 O lastResidual,
14 U numIters,
15 I myThid )
16 C !DESCRIPTION: \bv
17 C *==========================================================*
18 C | SUBROUTINE CG2D
19 C | o Two-dimensional grid problem conjugate-gradient
20 C | inverter (with preconditioner).
21 C *==========================================================*
22 C | Con. grad is an iterative procedure for solving Ax = b.
23 C | It requires the A be symmetric.
24 C | This implementation assumes A is a five-diagonal
25 C | matrix of the form that arises in the discrete
26 C | representation of the del^2 operator in a
27 C | two-dimensional space.
28 C | Notes:
29 C | ======
30 C | This implementation can support shared-memory
31 C | multi-threaded execution. In order to do this COMMON
32 C | blocks are used for many of the arrays - even ones that
33 C | are only used for intermedaite results. This design is
34 C | OK if you want to all the threads to collaborate on
35 C | solving the same problem. On the other hand if you want
36 C | the threads to solve several different problems
37 C | concurrently this implementation will not work.
38 C *==========================================================*
39 C \ev
40
41 C !USES:
42 IMPLICIT NONE
43 C === Global data ===
44 #include "SIZE.h"
45 #include "EEPARAMS.h"
46 #include "PARAMS.h"
47 #include "CG2D.h"
48 c#include "GRID.h"
49 c#include "SURFACE.h"
50
51 C !INPUT/OUTPUT PARAMETERS:
52 C === Routine arguments ===
53 C myThid :: Thread on which I am working.
54 C cg2d_b :: The source term or "right hand side"
55 C cg2d_x :: The solution
56 C firstResidual :: the initial residual before any iterations
57 C lastResidual :: the actual residual reached
58 C numIters :: Entry: the maximum number of iterations allowed
59 C Exit: the actual number of iterations used
60 _RL cg2d_b(1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
61 _RL cg2d_x(1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
62 _RL firstResidual
63 _RL lastResidual
64 INTEGER numIters
65 INTEGER myThid
66
67 C !LOCAL VARIABLES:
68 C === Local variables ====
69 C actualIts :: Number of iterations taken
70 C actualResidual :: residual
71 C bi, bj :: Block index in X and Y.
72 C eta_qrN :: Used in computing search directions
73 C eta_qrNM1 suffix N and NM1 denote current and
74 C cgBeta previous iterations respectively.
75 C alpha
76 C sumRHS :: Sum of right-hand-side. Sometimes this is a
77 C useful debuggin/trouble shooting diagnostic.
78 C For neumann problems sumRHS needs to be ~0.
79 C or they converge at a non-zero residual.
80 C err :: Measure of residual of Ax - b, usually the norm.
81 C I, J, it2d :: Loop counters ( it2d counts CG iterations )
82 INTEGER actualIts
83 _RL actualResidual
84 INTEGER bi, bj
85 INTEGER I, J, it2d
86 c INTEGER ks
87 _RL err,errTile
88 _RL eta_qrN,eta_qrNtile
89 _RL eta_qrNM1
90 _RL cgBeta
91 _RL alpha,alphaTile
92 _RL sumRHS,sumRHStile
93 _RL rhsMax
94 _RL rhsNorm
95 CEOP
96
97
98 CcnhDebugStarts
99 C CHARACTER*(MAX_LEN_FNAM) suff
100 CcnhDebugEnds
101
102
103 C-- Initialise inverter
104 eta_qrNM1 = 1. _d 0
105
106 CcnhDebugStarts
107 C _EXCH_XY_R8( cg2d_b, myThid )
108 C CALL PLOT_FIELD_XYRL( cg2d_b, 'CG2D.0 CG2D_B' , 1, myThid )
109 C suff = 'unnormalised'
110 C CALL WRITE_FLD_XY_RL ( 'cg2d_b.',suff, cg2d_b, 1, myThid)
111 C STOP
112 CcnhDebugEnds
113
114 C-- Normalise RHS
115 rhsMax = 0. _d 0
116 DO bj=myByLo(myThid),myByHi(myThid)
117 DO bi=myBxLo(myThid),myBxHi(myThid)
118 DO J=1,sNy
119 DO I=1,sNx
120 cg2d_b(I,J,bi,bj) = cg2d_b(I,J,bi,bj)*cg2dNorm
121 rhsMax = MAX(ABS(cg2d_b(I,J,bi,bj)),rhsMax)
122 ENDDO
123 ENDDO
124 ENDDO
125 ENDDO
126
127 IF (cg2dNormaliseRHS) THEN
128 C- Normalise RHS :
129 #ifdef LETS_MAKE_JAM
130 C _GLOBAL_MAX_R8( rhsMax, myThid )
131 rhsMax=1.
132 #else
133 _GLOBAL_MAX_R8( rhsMax, myThid )
134 #endif
135 rhsNorm = 1. _d 0
136 IF ( rhsMax .NE. 0. ) rhsNorm = 1. _d 0 / rhsMax
137 DO bj=myByLo(myThid),myByHi(myThid)
138 DO bi=myBxLo(myThid),myBxHi(myThid)
139 DO J=1,sNy
140 DO I=1,sNx
141 cg2d_b(I,J,bi,bj) = cg2d_b(I,J,bi,bj)*rhsNorm
142 cg2d_x(I,J,bi,bj) = cg2d_x(I,J,bi,bj)*rhsNorm
143 ENDDO
144 ENDDO
145 ENDDO
146 ENDDO
147 C- end Normalise RHS
148 ENDIF
149
150 C-- Update overlaps
151 _EXCH_XY_R8( cg2d_b, myThid )
152 _EXCH_XY_R8( cg2d_x, myThid )
153 CcnhDebugStarts
154 C CALL PLOT_FIELD_XYRL( cg2d_b, 'CG2D.1 CG2D_B' , 1, myThid )
155 C suff = 'normalised'
156 C CALL WRITE_FLD_XY_RL ( 'cg2d_b.',suff, cg2d_b, 1, myThid)
157 CcnhDebugEnds
158
159 C-- Initial residual calculation
160 err = 0. _d 0
161 sumRHS = 0. _d 0
162 DO bj=myByLo(myThid),myByHi(myThid)
163 DO bi=myBxLo(myThid),myBxHi(myThid)
164 sumRHStile = 0. _d 0
165 errTile = 0. _d 0
166 DO J=1,sNy
167 DO I=1,sNx
168 c ks = ksurfC(I,J,bi,bj)
169 cg2d_s(I,J,bi,bj) = 0.
170 cg2d_r(I,J,bi,bj) = cg2d_b(I,J,bi,bj) -
171 & (aW2d(I ,J ,bi,bj)*cg2d_x(I-1,J ,bi,bj)
172 & +aW2d(I+1,J ,bi,bj)*cg2d_x(I+1,J ,bi,bj)
173 & +aS2d(I ,J ,bi,bj)*cg2d_x(I ,J-1,bi,bj)
174 & +aS2d(I ,J+1,bi,bj)*cg2d_x(I ,J+1,bi,bj)
175 & +aC2d(I ,J ,bi,bj)*cg2d_x(I ,J ,bi,bj)
176 & )
177 c & -aW2d(I ,J ,bi,bj)*cg2d_x(I ,J ,bi,bj)
178 c & -aW2d(I+1,J ,bi,bj)*cg2d_x(I ,J ,bi,bj)
179 c & -aS2d(I ,J ,bi,bj)*cg2d_x(I ,J ,bi,bj)
180 c & -aS2d(I ,J+1,bi,bj)*cg2d_x(I ,J ,bi,bj)
181 c & -freeSurfFac*_rA(i,j,bi,bj)*deepFac2F(ks)*recip_Bo(i,j,bi,bj)
182 c & *cg2d_x(I ,J ,bi,bj)
183 c & /deltaTMom/deltaTfreesurf*cg2dNorm
184 c & )
185 errTile = errTile + cg2d_r(I,J,bi,bj)*cg2d_r(I,J,bi,bj)
186 sumRHStile = sumRHStile + cg2d_b(I,J,bi,bj)
187 ENDDO
188 ENDDO
189 sumRHS = sumRHS + sumRHStile
190 err = err + errTile
191 ENDDO
192 ENDDO
193 C _EXCH_XY_R8( cg2d_r, myThid )
194 #ifdef LETS_MAKE_JAM
195 CALL EXCH_XY_O1_R8_JAM( cg2d_r )
196 #else
197 CALL EXCH_XY_RL( cg2d_r, myThid )
198 #endif
199 C _EXCH_XY_R8( cg2d_s, myThid )
200 #ifdef LETS_MAKE_JAM
201 CALL EXCH_XY_O1_R8_JAM( cg2d_s )
202 #else
203 CALL EXCH_XY_RL( cg2d_s, myThid )
204 #endif
205 _GLOBAL_SUM_R8( sumRHS, myThid )
206 _GLOBAL_SUM_R8( err , myThid )
207 err = SQRT(err)
208 actualIts = 0
209 actualResidual = err
210
211 IF ( debugLevel .GE. debLevZero ) THEN
212 _BEGIN_MASTER( myThid )
213 WRITE(standardmessageunit,'(A,1P2E22.14)')
214 & ' cg2d: Sum(rhs),rhsMax = ', sumRHS,rhsMax
215 _END_MASTER( myThid )
216 ENDIF
217 C _BARRIER
218 c _BEGIN_MASTER( myThid )
219 c WRITE(standardmessageunit,'(A,I6,1PE30.14)')
220 c & ' CG2D iters, err = ',
221 c & actualIts, actualResidual
222 c _END_MASTER( myThid )
223 firstResidual=actualResidual
224
225 IF ( err .ge. cg2dTolerance ) then
226
227 C >>>>>>>>>>>>>>> BEGIN SOLVER <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
228 it2d=1
229
230 DO while ((it2d .le. numIters .and. err .ge. cg2dTolerance))
231
232 CcnhDebugStarts
233 C WRITE(*,*) ' CG2D: Iteration ',it2d-1,' residual = ',
234 C & actualResidual
235 CcnhDebugEnds
236 C-- Solve preconditioning equation and update
237 C-- conjugate direction vector "s".
238 eta_qrN = 0. _d 0
239 DO bj=myByLo(myThid),myByHi(myThid)
240 DO bi=myBxLo(myThid),myBxHi(myThid)
241 eta_qrNtile = 0. _d 0
242 DO J=1,sNy
243 DO I=1,sNx
244 cg2d_q(I,J,bi,bj) =
245 & pC(I ,J ,bi,bj)*cg2d_r(I ,J ,bi,bj)
246 & +pW(I ,J ,bi,bj)*cg2d_r(I-1,J ,bi,bj)
247 & +pW(I+1,J ,bi,bj)*cg2d_r(I+1,J ,bi,bj)
248 & +pS(I ,J ,bi,bj)*cg2d_r(I ,J-1,bi,bj)
249 & +pS(I ,J+1,bi,bj)*cg2d_r(I ,J+1,bi,bj)
250 CcnhDebugStarts
251 C cg2d_q(I,J,bi,bj) = cg2d_r(I ,J ,bi,bj)
252 CcnhDebugEnds
253 eta_qrNtile = eta_qrNtile
254 & +cg2d_q(I,J,bi,bj)*cg2d_r(I,J,bi,bj)
255 ENDDO
256 ENDDO
257 eta_qrN = eta_qrN + eta_qrNtile
258 ENDDO
259 ENDDO
260
261 _GLOBAL_SUM_R8(eta_qrN, myThid)
262 CcnhDebugStarts
263 C WRITE(*,*) ' CG2D: Iteration ',it2d-1,' eta_qrN = ',eta_qrN
264 CcnhDebugEnds
265 cgBeta = eta_qrN/eta_qrNM1
266 CcnhDebugStarts
267 C WRITE(*,*) ' CG2D: Iteration ',it2d-1,' beta = ',cgBeta
268 CcnhDebugEnds
269 eta_qrNM1 = eta_qrN
270
271 DO bj=myByLo(myThid),myByHi(myThid)
272 DO bi=myBxLo(myThid),myBxHi(myThid)
273 DO J=1,sNy
274 DO I=1,sNx
275 cg2d_s(I,J,bi,bj) = cg2d_q(I,J,bi,bj)
276 & + cgBeta*cg2d_s(I,J,bi,bj)
277 ENDDO
278 ENDDO
279 ENDDO
280 ENDDO
281
282 C-- Do exchanges that require messages i.e. between
283 C-- processes.
284 C _EXCH_XY_R8( cg2d_s, myThid )
285 #ifdef LETS_MAKE_JAM
286 CALL EXCH_XY_O1_R8_JAM( cg2d_s )
287 #else
288 CALL EXCH_XY_RL( cg2d_s, myThid )
289 #endif
290
291 C== Evaluate laplace operator on conjugate gradient vector
292 C== q = A.s
293 alpha = 0. _d 0
294 DO bj=myByLo(myThid),myByHi(myThid)
295 DO bi=myBxLo(myThid),myBxHi(myThid)
296 alphaTile = 0. _d 0
297 DO J=1,sNy
298 DO I=1,sNx
299 c ks = ksurfC(I,J,bi,bj)
300 cg2d_q(I,J,bi,bj) =
301 & aW2d(I ,J ,bi,bj)*cg2d_s(I-1,J ,bi,bj)
302 & +aW2d(I+1,J ,bi,bj)*cg2d_s(I+1,J ,bi,bj)
303 & +aS2d(I ,J ,bi,bj)*cg2d_s(I ,J-1,bi,bj)
304 & +aS2d(I ,J+1,bi,bj)*cg2d_s(I ,J+1,bi,bj)
305 & +aC2d(I ,J ,bi,bj)*cg2d_s(I ,J ,bi,bj)
306 c & -aW2d(I ,J ,bi,bj)*cg2d_s(I ,J ,bi,bj)
307 c & -aW2d(I+1,J ,bi,bj)*cg2d_s(I ,J ,bi,bj)
308 c & -aS2d(I ,J ,bi,bj)*cg2d_s(I ,J ,bi,bj)
309 c & -aS2d(I ,J+1,bi,bj)*cg2d_s(I ,J ,bi,bj)
310 c & -freeSurfFac*_rA(i,j,bi,bj)*deepFac2F(ks)*recip_Bo(i,j,bi,bj)
311 c & *cg2d_s(I ,J ,bi,bj)
312 c & /deltaTMom/deltaTfreesurf*cg2dNorm
313 alphaTile = alphaTile+cg2d_s(I,J,bi,bj)*cg2d_q(I,J,bi,bj)
314 ENDDO
315 ENDDO
316 alpha = alpha + alphaTile
317 ENDDO
318 ENDDO
319 _GLOBAL_SUM_R8(alpha,myThid)
320 CcnhDebugStarts
321 C WRITE(*,*) ' CG2D: Iteration ',it2d-1,' SUM(s*q)= ',alpha
322 CcnhDebugEnds
323 alpha = eta_qrN/alpha
324 CcnhDebugStarts
325 C WRITE(*,*) ' CG2D: Iteration ',it2d-1,' alpha= ',alpha
326 CcnhDebugEnds
327
328 C== Update solution and residual vectors
329 C Now compute "interior" points.
330 err = 0. _d 0
331 DO bj=myByLo(myThid),myByHi(myThid)
332 DO bi=myBxLo(myThid),myBxHi(myThid)
333 errTile = 0. _d 0
334 DO J=1,sNy
335 DO I=1,sNx
336 cg2d_x(I,J,bi,bj)=cg2d_x(I,J,bi,bj)+alpha*cg2d_s(I,J,bi,bj)
337 cg2d_r(I,J,bi,bj)=cg2d_r(I,J,bi,bj)-alpha*cg2d_q(I,J,bi,bj)
338 errTile = errTile+cg2d_r(I,J,bi,bj)*cg2d_r(I,J,bi,bj)
339 ENDDO
340 ENDDO
341 err = err + errTile
342 ENDDO
343 ENDDO
344
345 _GLOBAL_SUM_R8( err , myThid )
346 err = SQRT(err)
347 actualIts = it2d
348 actualResidual = err
349 C _EXCH_XY_R8(cg2d_r, myThid )
350 #ifdef LETS_MAKE_JAM
351 CALL EXCH_XY_O1_R8_JAM( cg2d_r )
352 #else
353 CALL EXCH_XY_RL( cg2d_r, myThid )
354 #endif
355 it2d=it2d+1
356
357 end do
358 end if
359
360 IF (cg2dNormaliseRHS) THEN
361 C-- Un-normalise the answer
362 DO bj=myByLo(myThid),myByHi(myThid)
363 DO bi=myBxLo(myThid),myBxHi(myThid)
364 DO J=1,sNy
365 DO I=1,sNx
366 cg2d_x(I ,J ,bi,bj) = cg2d_x(I ,J ,bi,bj)/rhsNorm
367 ENDDO
368 ENDDO
369 ENDDO
370 ENDDO
371 ENDIF
372
373 C The following exchange was moved up to solve_for_pressure
374 C for compatibility with TAMC.
375 C _EXCH_XY_R8(cg2d_x, myThid )
376 c _BEGIN_MASTER( myThid )
377 c WRITE(*,'(A,I6,1PE30.14)') ' CG2D iters, err = ',
378 c & actualIts, actualResidual
379 c _END_MASTER( myThid )
380
381 C-- Return parameters to caller
382 lastResidual=actualResidual
383 numIters=actualIts
384
385 CcnhDebugStarts
386 C CALL PLOT_FIELD_XYRL( cg2d_x, 'CALC_MOM_RHS CG2D_X' , 1, myThid )
387 C err = 0. _d 0
388 C DO bj=myByLo(myThid),myByHi(myThid)
389 C DO bi=myBxLo(myThid),myBxHi(myThid)
390 C DO J=1,sNy
391 C DO I=1,sNx
392 C cg2d_r(I,J,bi,bj) = cg2d_b(I,J,bi,bj) -
393 C & (aW2d(I ,J ,bi,bj)*cg2d_x(I-1,J ,bi,bj)
394 C & +aW2d(I+1,J ,bi,bj)*cg2d_x(I+1,J ,bi,bj)
395 C & +aS2d(I ,J ,bi,bj)*cg2d_x(I ,J-1,bi,bj)
396 C & +aS2d(I ,J+1,bi,bj)*cg2d_x(I ,J+1,bi,bj)
397 C & -aW2d(I ,J ,bi,bj)*cg2d_x(I ,J ,bi,bj)
398 C & -aW2d(I+1,J ,bi,bj)*cg2d_x(I ,J ,bi,bj)
399 C & -aS2d(I ,J ,bi,bj)*cg2d_x(I ,J ,bi,bj)
400 C & -aS2d(I ,J+1,bi,bj)*cg2d_x(I ,J ,bi,bj))
401 C err = err +
402 C & cg2d_r(I,J,bi,bj)*cg2d_r(I,J,bi,bj)
403 C ENDDO
404 C ENDDO
405 C ENDDO
406 C ENDDO
407 C _GLOBAL_SUM_R8( err , myThid )
408 C write(*,*) 'cg2d: Ax - b = ',SQRT(err)
409 CcnhDebugEnds
410
411 RETURN
412 END

  ViewVC Help
Powered by ViewVC 1.1.22