6 |
SUBROUTINE CG2D( |
SUBROUTINE CG2D( |
7 |
I cg2d_b, |
I cg2d_b, |
8 |
U cg2d_x, |
U cg2d_x, |
9 |
|
U tolerance, |
10 |
|
O residual, |
11 |
|
U numIters, |
12 |
I myThid ) |
I myThid ) |
13 |
C /==========================================================\ |
C /==========================================================\ |
14 |
C | SUBROUTINE CG2D | |
C | SUBROUTINE CG2D | |
43 |
#include "SURFACE.h" |
#include "SURFACE.h" |
44 |
|
|
45 |
C === Routine arguments === |
C === Routine arguments === |
46 |
C myThid - Thread on which I am working. |
C myThid - Thread on which I am working. |
47 |
INTEGER myThid |
C cg2d_b - The source term or "right hand side" |
48 |
_RL cg2d_x(1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy) |
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) |
_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 ==== |
C === Local variables ==== |
62 |
C actualIts - Number of iterations taken |
C actualIts - Number of iterations taken |
74 |
C err - Measure of residual of Ax - b, usually the norm. |
C err - Measure of residual of Ax - b, usually the norm. |
75 |
C I, J, N - Loop counters ( N counts CG iterations ) |
C I, J, N - Loop counters ( N counts CG iterations ) |
76 |
INTEGER actualIts |
INTEGER actualIts |
77 |
_RL actualResidual |
_RL actualResidual,initialResidual |
78 |
INTEGER bi, bj |
INTEGER bi, bj |
79 |
INTEGER I, J, it2d |
INTEGER I, J, it2d |
80 |
_RL err |
_RL err |
223 |
actualIts = 0 |
actualIts = 0 |
224 |
actualResidual = SQRT(err) |
actualResidual = SQRT(err) |
225 |
C _BARRIER |
C _BARRIER |
226 |
_BEGIN_MASTER( myThid ) |
c _BEGIN_MASTER( myThid ) |
227 |
WRITE(0,'(A,I6,1PE30.14)') ' CG2D iters, err = ', |
c WRITE(0,'(A,I6,1PE30.14)') ' CG2D iters, err = ', |
228 |
& actualIts, actualResidual |
c & actualIts, actualResidual |
229 |
_END_MASTER( ) |
c _END_MASTER( ) |
230 |
|
initialResidual=actualResidual |
231 |
|
|
232 |
C >>>>>>>>>>>>>>> BEGIN SOLVER <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< |
C >>>>>>>>>>>>>>> BEGIN SOLVER <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< |
233 |
DO 10 it2d=1, cg2dMaxIters |
DO 10 it2d=1, numIters |
234 |
|
|
235 |
CcnhDebugStarts |
CcnhDebugStarts |
236 |
C WRITE(0,*) ' CG2D: Iteration ',it2d-1,' residual = ', |
C WRITE(0,*) ' CG2D: Iteration ',it2d-1,' residual = ', |
237 |
C & actualResidual |
C & actualResidual |
238 |
CcnhDebugEnds |
CcnhDebugEnds |
239 |
IF ( err .LT. cg2dTargetResidual ) GOTO 11 |
IF ( err .LT. tolerance ) GOTO 11 |
240 |
C-- Solve preconditioning equation and update |
C-- Solve preconditioning equation and update |
241 |
C-- conjugate direction vector "s". |
C-- conjugate direction vector "s". |
242 |
eta_qrN = 0. _d 0 |
eta_qrN = 0. _d 0 |
386 |
C The following exchange was moved up to solve_for_pressure |
C The following exchange was moved up to solve_for_pressure |
387 |
C for compatibility with TAMC. |
C for compatibility with TAMC. |
388 |
C _EXCH_XY_R8(cg2d_x, myThid ) |
C _EXCH_XY_R8(cg2d_x, myThid ) |
389 |
_BEGIN_MASTER( myThid ) |
c _BEGIN_MASTER( myThid ) |
390 |
WRITE(0,'(A,I6,1PE30.14)') ' CG2D iters, err = ', |
c WRITE(0,'(A,I6,1PE30.14)') ' CG2D iters, err = ', |
391 |
& actualIts, actualResidual |
c & actualIts, actualResidual |
392 |
_END_MASTER( ) |
c _END_MASTER( ) |
393 |
|
|
394 |
|
C-- Return parameters to caller |
395 |
|
tolerance=initialResidual |
396 |
|
residual=actualResidual |
397 |
|
numIters=actualIts |
398 |
|
|
399 |
CcnhDebugStarts |
CcnhDebugStarts |
400 |
C CALL PLOT_FIELD_XYRL( cg2d_x, 'CALC_MOM_RHS CG2D_X' , 1, myThid ) |
C CALL PLOT_FIELD_XYRL( cg2d_x, 'CALC_MOM_RHS CG2D_X' , 1, myThid ) |