File indexing completed on 2024-09-20 05:10:16 UTC
view on githubraw file Latest commit e6e223b2 on 2024-09-19 22:01:15 UTC
1eadaea85b Jean*0001 #include "PACKAGES_CONFIG.h"
05b6d6742b Patr*0002 #include "CPP_OPTIONS.h"
517dbdc414 Jean*0003 #ifdef ALLOW_AUTODIFF
0004 # include "AUTODIFF_OPTIONS.h"
0005 #endif
aecc8b0f47 Mart*0006 #ifdef ALLOW_CTRL
0007 # include "CTRL_OPTIONS.h"
0008 #endif
05b6d6742b Patr*0009
0010
0011 #undef ALLOW_LOOP_DIRECTIVE
0012
0013
0014
08061169a5 Jean*0015 SUBROUTINE CG2D_NSA(
0539a6785e Jean*0016 U cg2d_b, cg2d_x,
0017 O firstResidual, minResidualSq, lastResidual,
0018 U numIters, nIterMin,
05b6d6742b Patr*0019 I myThid )
0020
aecc8b0f47 Mart*0021
08061169a5 Jean*0022
aecc8b0f47 Mart*0023
0024
05b6d6742b Patr*0025
aecc8b0f47 Mart*0026
0027
0028
08061169a5 Jean*0029
0030
aecc8b0f47 Mart*0031
0032
0033
0034
05b6d6742b Patr*0035
0036
0037
0038 IMPLICIT NONE
0039
0040 #include "SIZE.h"
0041 #include "EEPARAMS.h"
0042 #include "PARAMS.h"
0043 #include "CG2D.h"
7c50f07931 Mart*0044 #ifdef ALLOW_AUTODIFF_TAMC
05b6d6742b Patr*0045 # include "tamc.h"
0046 #endif
0047
0048
0539a6785e Jean*0049
0050
08061169a5 Jean*0051
0539a6785e Jean*0052
08061169a5 Jean*0053
0539a6785e Jean*0054
0055
0056
0057
9f5240b52a Jean*0058
05b6d6742b Patr*0059 _RL cg2d_b(1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
0060 _RL cg2d_x(1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
0061 _RL firstResidual
0539a6785e Jean*0062 _RL minResidualSq
05b6d6742b Patr*0063 _RL lastResidual
0064 INTEGER numIters
0539a6785e Jean*0065 INTEGER nIterMin
05b6d6742b Patr*0066 INTEGER myThid
0067
0068 #ifdef ALLOW_CG2D_NSA
0069
0539a6785e Jean*0070
0071
0072
0073
0074
0075
08061169a5 Jean*0076
0539a6785e Jean*0077
0078
13785d3a37 Jean*0079
0539a6785e Jean*0080
0081
0082
0083
aecc8b0f47 Mart*0084
0085
0086
0087
0088
0089
0090
08061169a5 Jean*0091 INTEGER bi, bj
0539a6785e Jean*0092 INTEGER i, j, it2d
0093 INTEGER actualIts
7c50f07931 Mart*0094 #ifdef ALLOW_AUTODIFF_TAMC
edb6656069 Mart*0095
0096 INTEGER ikey
9f5240b52a Jean*0097 #else
0098 _RL rhsNorm
7c50f07931 Mart*0099 #endif
13785d3a37 Jean*0100 _RL err_sq, errTile(nSx,nSy)
0101 _RL eta_qrN, eta_qrNtile(nSx,nSy)
0102 _RL eta_qrNM1, recip_eta_qrNM1
0103 _RL cgBeta, alpha
0104 _RL alphaSum,alphaTile(nSx,nSy)
0105 _RL sumRHS, sumRHStile(nSx,nSy)
9f5240b52a Jean*0106 _RL rhsMax
0539a6785e Jean*0107 CHARACTER*(MAX_LEN_MBUF) msgBuf
0108 LOGICAL printResidual
d68b36a485 Mart*0109 _RL cg2d_z(1:sNx,1:sNy,nSx,nSy)
0110 _RL cg2d_q(1:sNx,1:sNy,nSx,nSy)
0111
aecc8b0f47 Mart*0112
0113 _RL cg2d_r(1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
0114 _RL cg2d_s(1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
05b6d6742b Patr*0115
0116
0117 #ifdef ALLOW_AUTODIFF_TAMC
0118 IF ( numIters .GT. numItersMax ) THEN
13785d3a37 Jean*0119 WRITE(msgBuf,'(A,I10)')
0120 & 'CG2D_NSA: numIters > numItersMax =', numItersMax
0121 CALL PRINT_ERROR( msgBuf, myThid )
0122 STOP 'ABNORMAL END: S/R CG2D_NSA'
0123 ENDIF
0124 IF ( cg2dNormaliseRHS ) THEN
0125 WRITE(msgBuf,'(A)') 'CG2D_NSA: cg2dNormaliseRHS is disabled'
0126 CALL PRINT_ERROR( msgBuf, myThid )
0127 WRITE(msgBuf,'(A)')
0128 & 'set cg2dTargetResWunit (instead of cg2dTargetResidual)'
0129 CALL PRINT_ERROR( msgBuf, myThid )
0130 STOP 'ABNORMAL END: S/R CG2D_NSA'
05b6d6742b Patr*0131 ENDIF
0132 #endif /* ALLOW_AUTODIFF_TAMC */
0133
0539a6785e Jean*0134
0135 minResidualSq = -1. _d 0
0136 eta_qrNM1 = 1. _d 0
0137 recip_eta_qrNM1= 1. _d 0
05b6d6742b Patr*0138
0139
0140 rhsMax = 0. _d 0
0141 DO bj=myByLo(myThid),myByHi(myThid)
0142 DO bi=myBxLo(myThid),myBxHi(myThid)
0539a6785e Jean*0143 DO j=1,sNy
0144 DO i=1,sNx
0145 cg2d_b(i,j,bi,bj) = cg2d_b(i,j,bi,bj)*cg2dNorm
0146 rhsMax = MAX(ABS(cg2d_b(i,j,bi,bj)),rhsMax)
05b6d6742b Patr*0147 ENDDO
0148 ENDDO
0149 ENDDO
0150 ENDDO
0151
13785d3a37 Jean*0152 #ifndef ALLOW_AUTODIFF_TAMC
05b6d6742b Patr*0153 IF (cg2dNormaliseRHS) THEN
13785d3a37 Jean*0154
0155 _GLOBAL_MAX_RL( rhsMax, myThid )
05b6d6742b Patr*0156 rhsNorm = 1. _d 0
13785d3a37 Jean*0157 IF ( rhsMax .NE. 0. ) rhsNorm = 1. _d 0 / rhsMax
0158 DO bj=myByLo(myThid),myByHi(myThid)
0159 DO bi=myBxLo(myThid),myBxHi(myThid)
0160 DO j=1,sNy
0161 DO i=1,sNx
0162 cg2d_b(i,j,bi,bj) = cg2d_b(i,j,bi,bj)*rhsNorm
0163 cg2d_x(i,j,bi,bj) = cg2d_x(i,j,bi,bj)*rhsNorm
0164 ENDDO
05b6d6742b Patr*0165 ENDDO
0166 ENDDO
0167 ENDDO
13785d3a37 Jean*0168
05b6d6742b Patr*0169 ENDIF
13785d3a37 Jean*0170 #endif /* ndef ALLOW_AUTODIFF_TAMC */
05b6d6742b Patr*0171
0172
08061169a5 Jean*0173 CALL EXCH_XY_RL( cg2d_x, myThid )
05b6d6742b Patr*0174
0175
0176 #ifdef ALLOW_AUTODIFF_TAMC
20dee61641 Mart*0177
0178
05b6d6742b Patr*0179 #endif /* ALLOW_AUTODIFF_TAMC */
0180 DO bj=myByLo(myThid),myByHi(myThid)
0181 DO bi=myBxLo(myThid),myBxHi(myThid)
ff51fd05c6 Jean*0182
d68b36a485 Mart*0183 DO j=1,sNy
0184 DO i=1,sNx
ff51fd05c6 Jean*0185 cg2d_q(i,j,bi,bj) = 0. _d 0
0186 cg2d_z(i,j,bi,bj) = 0. _d 0
0187 ENDDO
0188 ENDDO
0189 DO j=1-OLy,sNy+OLy
0190 DO i=1-OLx,sNx+OLx
0191 cg2d_r(i,j,bi,bj) = 0. _d 0
0192 cg2d_s(i,j,bi,bj) = 0. _d 0
08061169a5 Jean*0193 ENDDO
0194 ENDDO
ff51fd05c6 Jean*0195 errTile(bi,bj) = 0. _d 0
0196 sumRHStile(bi,bj) = 0. _d 0
0539a6785e Jean*0197 DO j=1,sNy
0198 DO i=1,sNx
0199 cg2d_r(i,j,bi,bj) = cg2d_b(i,j,bi,bj) -
0200 & (aW2d(i ,j ,bi,bj)*cg2d_x(i-1,j ,bi,bj)
0201 & +aW2d(i+1,j ,bi,bj)*cg2d_x(i+1,j ,bi,bj)
0202 & +aS2d(i ,j ,bi,bj)*cg2d_x(i ,j-1,bi,bj)
0203 & +aS2d(i ,j+1,bi,bj)*cg2d_x(i ,j+1,bi,bj)
0204 & +aC2d(i ,j ,bi,bj)*cg2d_x(i ,j ,bi,bj)
05b6d6742b Patr*0205 & )
13785d3a37 Jean*0206 errTile(bi,bj) = errTile(bi,bj)
0207 & + cg2d_r(i,j,bi,bj)*cg2d_r(i,j,bi,bj)
0208 sumRHStile(bi,bj) = sumRHStile(bi,bj) + cg2d_b(i,j,bi,bj)
05b6d6742b Patr*0209 ENDDO
0210 ENDDO
0211 ENDDO
0212 ENDDO
0213
08061169a5 Jean*0214
0215 CALL EXCH_XY_RL ( cg2d_r, myThid )
13785d3a37 Jean*0216 CALL GLOBAL_SUM_TILE_RL( errTile, err_sq, myThid )
0217 CALL GLOBAL_SUM_TILE_RL( sumRHStile, sumRHS, myThid )
0539a6785e Jean*0218 actualIts = 0
0219 IF ( err_sq .NE. 0. ) THEN
0220 firstResidual = SQRT(err_sq)
0221 ELSE
0222 firstResidual = 0.
0223 ENDIF
08061169a5 Jean*0224
0539a6785e Jean*0225 printResidual = .FALSE.
08061169a5 Jean*0226 IF ( debugLevel .GE. debLevZero ) THEN
0227 _BEGIN_MASTER( myThid )
0539a6785e Jean*0228 printResidual = printResidualFreq.GE.1
e6e223b277 Jean*0229 WRITE(standardMessageUnit,'(A,1P2E22.14)')
aecc8b0f47 Mart*0230 & ' cg2d_nsa: Sum(rhs),rhsMax = ', sumRHS,rhsMax
08061169a5 Jean*0231 _END_MASTER( myThid )
0232 ENDIF
05b6d6742b Patr*0233
0234
0235
0236 #if ((defined ALLOW_AUTODIFF_TAMC) && (defined ALLOW_LOOP_DIRECTIVE))
0237
0238 #endif /* ALLOW_AUTODIFF_TAMC and ALLOW_LOOP_DIRECTIVE */
0239 DO it2d=1, numIters
0240 #ifdef ALLOW_LOOP_DIRECTIVE
0241
08061169a5 Jean*0242
05b6d6742b Patr*0243
0244 #endif /* ALLOW_LOOP_DIRECTIVE */
0245
0246 #ifdef ALLOW_AUTODIFF_TAMC
edb6656069 Mart*0247 ikey = it2d + (ikey_dynamics-1)*numItersMax
0248
05b6d6742b Patr*0249 #endif /* ALLOW_AUTODIFF_TAMC */
aecc8b0f47 Mart*0250 IF ( it2d .LE. cg2dMinItersNSA .OR.
0251 & err_sq .GE. cg2dTolerance_sq ) THEN
05b6d6742b Patr*0252
0253
0254
0255 #ifdef ALLOW_AUTODIFF_TAMC
edb6656069 Mart*0256
0257
05b6d6742b Patr*0258 #endif /* ALLOW_AUTODIFF_TAMC */
0259 DO bj=myByLo(myThid),myByHi(myThid)
0260 DO bi=myBxLo(myThid),myBxHi(myThid)
13785d3a37 Jean*0261 eta_qrNtile(bi,bj) = 0. _d 0
0539a6785e Jean*0262 DO j=1,sNy
0263 DO i=1,sNx
0264 cg2d_z(i,j,bi,bj) =
0265 & pC(i ,j ,bi,bj)*cg2d_r(i ,j ,bi,bj)
0266 & +pW(i ,j ,bi,bj)*cg2d_r(i-1,j ,bi,bj)
0267 & +pW(i+1,j ,bi,bj)*cg2d_r(i+1,j ,bi,bj)
0268 & +pS(i ,j ,bi,bj)*cg2d_r(i ,j-1,bi,bj)
0269 & +pS(i ,j+1,bi,bj)*cg2d_r(i ,j+1,bi,bj)
13785d3a37 Jean*0270 eta_qrNtile(bi,bj) = eta_qrNtile(bi,bj)
0539a6785e Jean*0271 & +cg2d_z(i,j,bi,bj)*cg2d_r(i,j,bi,bj)
05b6d6742b Patr*0272 ENDDO
0273 ENDDO
0274 ENDDO
0275 ENDDO
0276
13785d3a37 Jean*0277 CALL GLOBAL_SUM_TILE_RL( eta_qrNtile,eta_qrN,myThid )
05b6d6742b Patr*0278 #ifdef ALLOW_AUTODIFF_TAMC
edb6656069 Mart*0279
0280
05b6d6742b Patr*0281 #endif /* ALLOW_AUTODIFF_TAMC */
0282
0283 cgBeta = eta_qrN*recip_eta_qrNM1
13785d3a37 Jean*0284
05b6d6742b Patr*0285
0286
aecc8b0f47 Mart*0287 recip_eta_qrNM1 = 0. _d 0
0288 IF ( eta_qrN .NE. 0. _d 0 ) recip_eta_qrNM1 = 1. _d 0/eta_qrN
05b6d6742b Patr*0289
0290 DO bj=myByLo(myThid),myByHi(myThid)
0291 DO bi=myBxLo(myThid),myBxHi(myThid)
0539a6785e Jean*0292 DO j=1,sNy
0293 DO i=1,sNx
0294 cg2d_s(i,j,bi,bj) = cg2d_z(i,j,bi,bj)
0295 & + cgBeta*cg2d_s(i,j,bi,bj)
05b6d6742b Patr*0296 ENDDO
0297 ENDDO
0298 ENDDO
0299 ENDDO
0300
0301
0302
08061169a5 Jean*0303
0304 CALL EXCH_XY_RL ( cg2d_s, myThid )
05b6d6742b Patr*0305
0306
0307
0308 #ifdef ALLOW_AUTODIFF_TAMC
0309 #ifndef ALLOW_LOOP_DIRECTIVE
edb6656069 Mart*0310
05b6d6742b Patr*0311 #endif /* not ALLOW_LOOP_DIRECTIVE */
0312 #endif /* ALLOW_AUTODIFF_TAMC */
0313 DO bj=myByLo(myThid),myByHi(myThid)
0314 DO bi=myBxLo(myThid),myBxHi(myThid)
13785d3a37 Jean*0315 alphaTile(bi,bj) = 0. _d 0
0539a6785e Jean*0316 DO j=1,sNy
0317 DO i=1,sNx
0318 cg2d_q(i,j,bi,bj) =
0319 & aW2d(i ,j ,bi,bj)*cg2d_s(i-1,j ,bi,bj)
0320 & +aW2d(i+1,j ,bi,bj)*cg2d_s(i+1,j ,bi,bj)
0321 & +aS2d(i ,j ,bi,bj)*cg2d_s(i ,j-1,bi,bj)
0322 & +aS2d(i ,j+1,bi,bj)*cg2d_s(i ,j+1,bi,bj)
0323 & +aC2d(i ,j ,bi,bj)*cg2d_s(i ,j ,bi,bj)
13785d3a37 Jean*0324 alphaTile(bi,bj) = alphaTile(bi,bj)
0325 & + cg2d_s(i,j,bi,bj)*cg2d_q(i,j,bi,bj)
05b6d6742b Patr*0326 ENDDO
0327 ENDDO
0328 ENDDO
0329 ENDDO
13785d3a37 Jean*0330 CALL GLOBAL_SUM_TILE_RL( alphaTile, alphaSum, myThid )
aecc8b0f47 Mart*0331 alpha = 0. _d 0
0332 IF ( alphaSum .NE. 0 _d 0 ) alpha = eta_qrN/alphaSum
08061169a5 Jean*0333
0539a6785e Jean*0334
05b6d6742b Patr*0335
0336 #ifdef ALLOW_AUTODIFF_TAMC
0337 #ifndef ALLOW_LOOP_DIRECTIVE
edb6656069 Mart*0338
05b6d6742b Patr*0339 #endif /* ALLOW_LOOP_DIRECTIVE */
0340 #endif /* ALLOW_AUTODIFF_TAMC */
0341 DO bj=myByLo(myThid),myByHi(myThid)
0342 DO bi=myBxLo(myThid),myBxHi(myThid)
13785d3a37 Jean*0343 errTile(bi,bj) = 0. _d 0
0539a6785e Jean*0344 DO j=1,sNy
0345 DO i=1,sNx
0346 cg2d_x(i,j,bi,bj)=cg2d_x(i,j,bi,bj)+alpha*cg2d_s(i,j,bi,bj)
0347 cg2d_r(i,j,bi,bj)=cg2d_r(i,j,bi,bj)-alpha*cg2d_q(i,j,bi,bj)
13785d3a37 Jean*0348 errTile(bi,bj) = errTile(bi,bj)
0349 & + cg2d_r(i,j,bi,bj)*cg2d_r(i,j,bi,bj)
05b6d6742b Patr*0350 ENDDO
0351 ENDDO
0352 ENDDO
0353 ENDDO
0539a6785e Jean*0354 actualIts = it2d
05b6d6742b Patr*0355
13785d3a37 Jean*0356 CALL GLOBAL_SUM_TILE_RL( errTile, err_sq, myThid )
0539a6785e Jean*0357 IF ( printResidual ) THEN
0358 IF ( MOD( it2d-1, printResidualFreq ).EQ.0 ) THEN
0359 WRITE(msgBuf,'(A,I6,A,1PE21.14)')
aecc8b0f47 Mart*0360 & ' cg2d_nsa: iter=', it2d, ' ; resid.= ', SQRT(err_sq)
0539a6785e Jean*0361 CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
0362 & SQUEEZE_RIGHT, myThid )
0363 ENDIF
0364 ENDIF
05b6d6742b Patr*0365
08061169a5 Jean*0366
0367 CALL EXCH_XY_RL ( cg2d_r, myThid )
05b6d6742b Patr*0368
0539a6785e Jean*0369
05b6d6742b Patr*0370 ENDIF
0371 ENDDO
0372
13785d3a37 Jean*0373 #ifndef ALLOW_AUTODIFF_TAMC
05b6d6742b Patr*0374 IF (cg2dNormaliseRHS) THEN
0375
0376 DO bj=myByLo(myThid),myByHi(myThid)
0377 DO bi=myBxLo(myThid),myBxHi(myThid)
0539a6785e Jean*0378 DO j=1,sNy
0379 DO i=1,sNx
0380 cg2d_x(i,j,bi,bj) = cg2d_x(i,j,bi,bj)/rhsNorm
05b6d6742b Patr*0381 ENDDO
0382 ENDDO
0383 ENDDO
0384 ENDDO
0385 ENDIF
13785d3a37 Jean*0386 #endif /* ndef ALLOW_AUTODIFF_TAMC */
05b6d6742b Patr*0387
0388
0539a6785e Jean*0389 IF ( err_sq .NE. 0. ) THEN
0390 lastResidual = SQRT(err_sq)
0391 ELSE
0392 lastResidual = 0.
0393 ENDIF
0394 numIters = actualIts
05b6d6742b Patr*0395
0396 #endif /* ALLOW_CG2D_NSA */
0397 RETURN
0398 END
0399
0539a6785e Jean*0400 #if ((defined ALLOW_AUTODIFF_TAMC) && (defined ALLOW_LOOP_DIRECTIVE))
13785d3a37 Jean*0401
05b6d6742b Patr*0402
0403
13785d3a37 Jean*0404
aecc8b0f47 Mart*0405 SUBROUTINE ADSTORE(chardum,int1,idow,int2,int3,icount)
05b6d6742b Patr*0406
aecc8b0f47 Mart*0407 IMPLICIT NONE
05b6d6742b Patr*0408
0409 #include "SIZE.h"
0410 #include "tamc.h"
0411
aecc8b0f47 Mart*0412 CHARACTER*(*) chardum
0413 INTEGER int1, int2, int3, idow, icount
05b6d6742b Patr*0414
08061169a5 Jean*0415
05b6d6742b Patr*0416
aecc8b0f47 Mart*0417 INTEGER nidow
05b6d6742b Patr*0418 #ifdef ALLOW_TAMC_CHECKPOINTING
aecc8b0f47 Mart*0419 PARAMETER ( nidow = 2*nchklev_1*nchklev_2*nchklev_3 )
05b6d6742b Patr*0420 #else
aecc8b0f47 Mart*0421 PARAMETER ( nidow = 1000000 )
05b6d6742b Patr*0422 #endif /* ALLOW_TAMC_CHECKPOINTING */
aecc8b0f47 Mart*0423 INTEGER istoreidow(nidow)
0424 COMMON /istorecommon/ istoreidow
05b6d6742b Patr*0425
aecc8b0f47 Mart*0426 PRINT *, 'ADSTORE: ', chardum, int1, idow, int2, int3, icount
05b6d6742b Patr*0427
aecc8b0f47 Mart*0428 IF ( icount .GT. nidow ) THEN
0429 PRINT *, 'adstore: error: icount > nidow = ', nidow
0430 STOP 'ABNORMAL STOP in ADSTORE'
0431 ENDIF
05b6d6742b Patr*0432
0433 istoreidow(icount) = idow
0434
aecc8b0f47 Mart*0435 RETURN
0436 END
05b6d6742b Patr*0437
aecc8b0f47 Mart*0438 SUBROUTINE ADRESTO(chardum,int1,idow,int2,int3,icount)
05b6d6742b Patr*0439
aecc8b0f47 Mart*0440 IMPLICIT NONE
05b6d6742b Patr*0441
0442 #include "SIZE.h"
0443 #include "tamc.h"
0444
aecc8b0f47 Mart*0445 CHARACTER*(*) chardum
0446 INTEGER int1, int2, int3, idow, icount
05b6d6742b Patr*0447
08061169a5 Jean*0448
05b6d6742b Patr*0449
aecc8b0f47 Mart*0450 INTEGER nidow
05b6d6742b Patr*0451 #ifdef ALLOW_TAMC_CHECKPOINTING
aecc8b0f47 Mart*0452 PARAMETER ( nidow = 2*nchklev_1*nchklev_2*nchklev_3 )
05b6d6742b Patr*0453 #else
aecc8b0f47 Mart*0454 PARAMETER ( nidow = 1000000 )
05b6d6742b Patr*0455 #endif /* ALLOW_TAMC_CHECKPOINTING */
aecc8b0f47 Mart*0456 INTEGER istoreidow(nidow)
0457 COMMON /istorecommon/ istoreidow
05b6d6742b Patr*0458
aecc8b0f47 Mart*0459 PRINT *, 'ADRESTO: ', chardum, int1, idow, int2, int3, icount
05b6d6742b Patr*0460
aecc8b0f47 Mart*0461 IF ( icount .GT. nidow ) THEN
0462 PRINT *, 'ADRESTO: error: icount > nidow = ', nidow
0463 STOP 'ABNORMAL STOP in ADRESTO'
0464 ENDIF
05b6d6742b Patr*0465
0466 idow = istoreidow(icount)
0467
aecc8b0f47 Mart*0468 RETURN
0469 END
0539a6785e Jean*0470 #endif /* ALLOW_AUTODIFF_TAMC and ALLOW_LOOP_DIRECTIVE */