/[MITgcm]/MITgcm/model/src/cg2d.F
ViewVC logotype

Contents of /MITgcm/model/src/cg2d.F

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


Revision 1.31 - (show annotations) (download)
Tue Apr 10 22:35:25 2001 UTC (23 years, 1 month ago) by heimbach
Branch: MAIN
CVS Tags: checkpoint38
Changes since 1.30: +11 -11 lines
See doc/tag-index and doc/notes_c37_adj.txt
Preparation for stand-alone autodifferentiability.

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

  ViewVC Help
Powered by ViewVC 1.1.22