| 1 | c-----------------------------------------------------------------------
|
|---|
| 2 | c
|
|---|
| 3 | c TO DO: Need to monitor dssum when using DG. (11/29/15, pff)
|
|---|
| 4 | c
|
|---|
| 5 | c Specifically - line 837 in navier4.f
|
|---|
| 6 | c - line 537 in hmholtz.f
|
|---|
| 7 | c
|
|---|
| 8 | c-----------------------------------------------------------------------
|
|---|
| 9 | subroutine setrhs(p,h1,h2,h2inv)
|
|---|
| 10 | C
|
|---|
| 11 | C Project rhs onto best fit in the "E" norm.
|
|---|
| 12 | C
|
|---|
| 13 | include 'SIZE'
|
|---|
| 14 | include 'INPUT'
|
|---|
| 15 | include 'MASS'
|
|---|
| 16 | include 'SOLN'
|
|---|
| 17 | include 'TSTEP'
|
|---|
| 18 | C
|
|---|
| 19 | REAL P (LX2,LY2,LZ2,LELV)
|
|---|
| 20 | REAL H1 (LX1,LY1,LZ1,LELV)
|
|---|
| 21 | REAL H2 (LX1,LY1,LZ1,LELV)
|
|---|
| 22 | REAL H2INV(LX1,LY1,LZ1,LELV)
|
|---|
| 23 | C
|
|---|
| 24 | logical ifdump
|
|---|
| 25 | save ifdump
|
|---|
| 26 | data ifdump /.false./
|
|---|
| 27 | C
|
|---|
| 28 | PARAMETER (LTOT2=LX2*LY2*LZ2*LELV)
|
|---|
| 29 | COMMON /ORTHOV/ RHS(LTOT2,MXPREV)
|
|---|
| 30 | COMMON /ORTHOX/ Pbar(LTOT2),Pnew(LTOT2)
|
|---|
| 31 | COMMON /ORTHOS/ ALPHA(Mxprev), WORK(Mxprev), ALPHAN, DTLAST
|
|---|
| 32 | COMMON /ORTHOI/ Nprev,Mprev
|
|---|
| 33 | REAL ALPHA,WORK
|
|---|
| 34 | C
|
|---|
| 35 | C
|
|---|
| 36 | integer icalld
|
|---|
| 37 | save icalld
|
|---|
| 38 | data icalld/0/
|
|---|
| 39 | C
|
|---|
| 40 | C First call, we have no vectors to orthogonalize against.
|
|---|
| 41 | IF (ICALLD.EQ.0) THEN
|
|---|
| 42 | icalld=icalld+1
|
|---|
| 43 | Nprev=0
|
|---|
| 44 | Mprev=param(93)
|
|---|
| 45 | Mprev=min(Mprev,Mxprev)
|
|---|
| 46 | endif
|
|---|
| 47 | C
|
|---|
| 48 | C Diag to see how much reduction in the residual is attained.
|
|---|
| 49 | C
|
|---|
| 50 | NTOT2 = lx2*ly2*lz2*NELV
|
|---|
| 51 | ALPHA1 = GLSC3(p,p,bm2inv,NTOT2)
|
|---|
| 52 | if (alpha1.gt.0) ALPHA1 = sqrt(alpha1/volvm2)
|
|---|
| 53 | C
|
|---|
| 54 | C Update rhs's if E-matrix has changed
|
|---|
| 55 | C
|
|---|
| 56 | CALL UPDRHSE(P,H1,H2,H2INV,ierr)
|
|---|
| 57 | if (ierr.eq.1) Nprev=0
|
|---|
| 58 | C
|
|---|
| 59 | C Perform Gram-Schmidt for previous rhs's.
|
|---|
| 60 | C
|
|---|
| 61 | DO 10 I=1,Nprev
|
|---|
| 62 | ALPHA(i) = VLSC2(P,RHS(1,i),NTOT2)
|
|---|
| 63 | 10 CONTINUE
|
|---|
| 64 | C
|
|---|
| 65 | IF (Nprev.GT.0) CALL gop(alpha,WORK,'+ ',Nprev)
|
|---|
| 66 | C
|
|---|
| 67 | CALL RZERO(Pbar,NTOT2)
|
|---|
| 68 | DO 20 I=1,Nprev
|
|---|
| 69 | alphas = alpha(i)
|
|---|
| 70 | CALL ADD2S2(Pbar,RHS(1,i),alphas,NTOT2)
|
|---|
| 71 | 20 CONTINUE
|
|---|
| 72 | C
|
|---|
| 73 | if (Nprev.gt.0) then
|
|---|
| 74 | INTETYPE = 1
|
|---|
| 75 | CALL CDABDTP(Pnew,Pbar,H1,H2,H2INV,INTETYPE)
|
|---|
| 76 | CALL SUB2 (P,Pnew,NTOT2)
|
|---|
| 77 | C ................................................................
|
|---|
| 78 | C Diag.
|
|---|
| 79 | ALPHA2 = GLSC3(p,p,bm2inv,NTOT2)
|
|---|
| 80 | if (alpha2.gt.0) ALPHA2 = sqrt(alpha2/volvm2)
|
|---|
| 81 | ratio = alpha1/alpha2
|
|---|
| 82 | n10=min(10,nprev)
|
|---|
| 83 | c IF (NIO.EQ.0) WRITE(6,11)ISTEP,Nprev,(ALPHA(I),I=1,N10)
|
|---|
| 84 | c IF (NIO.EQ.0) WRITE(6,12) ISTEP,nprev,ALPHA1,ALPHA2,ratio
|
|---|
| 85 | 11 FORMAT(2I5,' alpha:',1p10e12.4)
|
|---|
| 86 | 12 FORMAT(I6,i4,1p3e12.4,' alph12')
|
|---|
| 87 | C
|
|---|
| 88 | c if (alpha1.gt.0.001 .and. .not.ifdump) then
|
|---|
| 89 | c IF (NID.EQ.0) WRITE(6,*) 'alph1 large ... '
|
|---|
| 90 | c if (istep.gt.10) then
|
|---|
| 91 | c IF (NID.EQ.0) WRITE(6,*) ' ... dumping'
|
|---|
| 92 | c call prepost (.true.,' ')
|
|---|
| 93 | c ifdump = .true.
|
|---|
| 94 | c else
|
|---|
| 95 | c IF (NID.EQ.0) WRITE(6,*) ' ... doing nothing'
|
|---|
| 96 | c endif
|
|---|
| 97 | c endif
|
|---|
| 98 | c if (alpha1.gt.0.1.and.istep.gt.10) then
|
|---|
| 99 | c IF (NID.EQ.0) WRITE(6,*) 'alph1 too large ... aborting'
|
|---|
| 100 | c call prepost (.true.,' ')
|
|---|
| 101 | c call exitt
|
|---|
| 102 | c endif
|
|---|
| 103 | C ................................................................
|
|---|
| 104 | endif
|
|---|
| 105 | C
|
|---|
| 106 | C
|
|---|
| 107 | return
|
|---|
| 108 | end
|
|---|
| 109 | c-----------------------------------------------------------------------
|
|---|
| 110 | subroutine gensoln(p,h1,h2,h2inv)
|
|---|
| 111 | C
|
|---|
| 112 | C Reconstruct the solution to the original problem by adding back
|
|---|
| 113 | C the previous solutions
|
|---|
| 114 | C know the soln.
|
|---|
| 115 | C
|
|---|
| 116 | include 'SIZE'
|
|---|
| 117 | PARAMETER (LTOT2=LX2*LY2*LZ2*LELV)
|
|---|
| 118 | COMMON /ORTHOV/ RHS(LTOT2,MXPREV)
|
|---|
| 119 | COMMON /ORTHOX/ Pbar(LTOT2),Pnew(LTOT2)
|
|---|
| 120 | COMMON /ORTHOS/ ALPHA(Mxprev), WORK(Mxprev), ALPHAN, DTLAST
|
|---|
| 121 | COMMON /ORTHOI/ Nprev,Mprev
|
|---|
| 122 | REAL ALPHA,WORK
|
|---|
| 123 |
|
|---|
| 124 | REAL P (LX2,LY2,LZ2,LELV)
|
|---|
| 125 | REAL H1 (LX1,LY1,LZ1,LELV)
|
|---|
| 126 | REAL H2 (LX1,LY1,LZ1,LELV)
|
|---|
| 127 | REAL H2INV(LX1,LY1,LZ1,LELV)
|
|---|
| 128 | C
|
|---|
| 129 | NTOT2=lx2*ly2*lz2*NELV
|
|---|
| 130 | C
|
|---|
| 131 | C First, save current solution
|
|---|
| 132 | C
|
|---|
| 133 | CALL COPY (Pnew,P,NTOT2)
|
|---|
| 134 | C
|
|---|
| 135 | C Reconstruct solution
|
|---|
| 136 | C
|
|---|
| 137 | CALL ADD2(P,Pbar,NTOT2)
|
|---|
| 138 | C
|
|---|
| 139 | C Update the set of <p,rhs>
|
|---|
| 140 | C
|
|---|
| 141 | CALL UPDTSET(P,H1,H2,H2INV,ierr)
|
|---|
| 142 | if (ierr.eq.1) Nprev = 0
|
|---|
| 143 | c
|
|---|
| 144 | return
|
|---|
| 145 | end
|
|---|
| 146 | c-----------------------------------------------------------------------
|
|---|
| 147 | subroutine updtset(p,h1,h2,h2inv,IERR)
|
|---|
| 148 | C
|
|---|
| 149 | C Update the set of rhs's and the corresponding p-set:
|
|---|
| 150 | C
|
|---|
| 151 | C . Standard case is to add P_new, and RHS_new = E*P_new
|
|---|
| 152 | C
|
|---|
| 153 | C . However, when Nprev=Mprev (max. allowed), we throw out
|
|---|
| 154 | C the old set, and set P_1 = P, RHS_1=E*P_1
|
|---|
| 155 | C
|
|---|
| 156 | C . Other schemes are possible, e.g., let's save a bunch of
|
|---|
| 157 | C old vectors, perhaps chosen wisely via P.O.D.
|
|---|
| 158 | C
|
|---|
| 159 | C
|
|---|
| 160 | include 'SIZE'
|
|---|
| 161 | include 'INPUT'
|
|---|
| 162 | include 'MASS'
|
|---|
| 163 | PARAMETER (LTOT2=LX2*LY2*LZ2*LELV)
|
|---|
| 164 | COMMON /ORTHOV/ RHS(LTOT2,MXPREV)
|
|---|
| 165 | COMMON /ORTHOX/ Pbar(LTOT2),Pnew(LTOT2)
|
|---|
| 166 | COMMON /ORTHOS/ ALPHA(Mxprev), WORK(Mxprev), ALPHAN, DTLAST
|
|---|
| 167 | COMMON /ORTHOI/ Nprev,Mprev
|
|---|
| 168 |
|
|---|
| 169 | REAL ALPHA,WORK
|
|---|
| 170 |
|
|---|
| 171 | REAL P (LX2,LY2,LZ2,LELV)
|
|---|
| 172 | REAL H1 (LX1,LY1,LZ1,LELV)
|
|---|
| 173 | REAL H2 (LX1,LY1,LZ1,LELV)
|
|---|
| 174 | REAL H2INV(LX1,LY1,LZ1,LELV)
|
|---|
| 175 | C
|
|---|
| 176 | NTOT2=lx2*ly2*lz2*NELV
|
|---|
| 177 | C
|
|---|
| 178 | IF (Nprev.EQ.Mprev) THEN
|
|---|
| 179 | CALL COPY(Pnew,P,NTOT2)
|
|---|
| 180 | Nprev=0
|
|---|
| 181 | endif
|
|---|
| 182 | C
|
|---|
| 183 | C Increment solution set
|
|---|
| 184 | Nprev = Nprev+1
|
|---|
| 185 | C
|
|---|
| 186 | CALL COPY (RHS(1,Nprev),Pnew,NTOT2)
|
|---|
| 187 |
|
|---|
| 188 | C Orthogonalize rhs against previous rhs and normalize
|
|---|
| 189 |
|
|---|
| 190 | write(6,*) istep,nprev,' call econj2'
|
|---|
| 191 | call econj (nprev,h1,h2,h2inv,ierr)
|
|---|
| 192 | C call echeck(nprev,h1,h2,h2inv,intetype)
|
|---|
| 193 | C
|
|---|
| 194 | c Save last sol'n
|
|---|
| 195 | CALL COPY(Pnew,P,NTOT2)
|
|---|
| 196 | C
|
|---|
| 197 | return
|
|---|
| 198 | end
|
|---|
| 199 | c-----------------------------------------------------------------------
|
|---|
| 200 | subroutine econj(kprev,h1,h2,h2inv,ierr)
|
|---|
| 201 | C
|
|---|
| 202 | C Orthogonalize the rhs wrt previous rhs's for which we already
|
|---|
| 203 | C know the soln.
|
|---|
| 204 | C
|
|---|
| 205 | include 'SIZE'
|
|---|
| 206 | include 'INPUT'
|
|---|
| 207 | include 'MASS'
|
|---|
| 208 | include 'SOLN'
|
|---|
| 209 | include 'TSTEP'
|
|---|
| 210 | C
|
|---|
| 211 | REAL H1 (LX1,LY1,LZ1,LELV)
|
|---|
| 212 | REAL H2 (LX1,LY1,LZ1,LELV)
|
|---|
| 213 | REAL H2INV(LX1,LY1,LZ1,LELV)
|
|---|
| 214 | C
|
|---|
| 215 | PARAMETER (LTOT2=LX2*LY2*LZ2*LELV)
|
|---|
| 216 | COMMON /ORTHOV/ RHS(LTOT2,MXPREV)
|
|---|
| 217 | COMMON /ORTHOX/ Pbar(LTOT2),Pnew(LTOT2),Pbrr(ltot2)
|
|---|
| 218 | COMMON /ORTHOS/ ALPHA(Mxprev), WORK(Mxprev), ALPHAN, DTLAST
|
|---|
| 219 | COMMON /ORTHOI/ Nprev,Mprev
|
|---|
| 220 | REAL ALPHA,WORK
|
|---|
| 221 | real ALPHAd
|
|---|
| 222 | C
|
|---|
| 223 | C
|
|---|
| 224 | ierr = 0
|
|---|
| 225 | NTOT2 = lx2*ly2*lz2*NELV
|
|---|
| 226 | INTETYPE=1
|
|---|
| 227 | C
|
|---|
| 228 | C Gram Schmidt, w re-orthogonalization
|
|---|
| 229 | C
|
|---|
| 230 | npass=1
|
|---|
| 231 | if (abs(param(105)).eq.2) npass=2
|
|---|
| 232 | do ipass=1,npass
|
|---|
| 233 | c
|
|---|
| 234 | CALL CDABDTP(Pbrr,RHS(1,Kprev),H1,H2,H2INV,INTETYPE)
|
|---|
| 235 | C
|
|---|
| 236 | C Compute part of the norm
|
|---|
| 237 | Alphad = GLSC2(RHS(1,Kprev),Pbrr,NTOT2)
|
|---|
| 238 | C
|
|---|
| 239 | C Gram-Schmidt
|
|---|
| 240 | Kprev1=Kprev-1
|
|---|
| 241 | DO 10 I=1,Kprev1
|
|---|
| 242 | ALPHA(I) = VLSC2(Pbrr,RHS(1,i),NTOT2)
|
|---|
| 243 | 10 CONTINUE
|
|---|
| 244 | IF (Kprev1.GT.0) CALL gop(alpha,WORK,'+ ',Kprev1)
|
|---|
| 245 | C
|
|---|
| 246 | DO 20 I=1,Kprev1
|
|---|
| 247 | alpham = -alpha(i)
|
|---|
| 248 | CALL ADD2S2(RHS(1,Kprev),RHS (1,i),alpham,NTOT2)
|
|---|
| 249 | Alphad = Alphad - alpha(i)**2
|
|---|
| 250 | 20 CONTINUE
|
|---|
| 251 | enddo
|
|---|
| 252 | C
|
|---|
| 253 | C .Normalize new element in P~
|
|---|
| 254 | C
|
|---|
| 255 | if (ALPHAd.le.0.0) then
|
|---|
| 256 | if (nio.eq.0)
|
|---|
| 257 | $ write(6,*) 'ERROR: alphad .le. 0 in econj',alphad,Kprev
|
|---|
| 258 | ierr = 1
|
|---|
| 259 | return
|
|---|
| 260 | endif
|
|---|
| 261 | ALPHAd = 1.0/SQRT(ALPHAd)
|
|---|
| 262 | ALPHAN = Alphad
|
|---|
| 263 | CALL CMULT(RHS (1,Kprev),alphan,NTOT2)
|
|---|
| 264 | C
|
|---|
| 265 | return
|
|---|
| 266 | end
|
|---|
| 267 | c-----------------------------------------------------------------------
|
|---|
| 268 | subroutine chkptol
|
|---|
| 269 | C--------------------------------------------------------------------
|
|---|
| 270 | C
|
|---|
| 271 | C Check pressure tolerance for transient case.
|
|---|
| 272 | C
|
|---|
| 273 | C pff 6/20/92
|
|---|
| 274 | C This routine has been modified for diagnostic purposes only.
|
|---|
| 275 | C It can be replaced with the standard nekton version.
|
|---|
| 276 | C
|
|---|
| 277 | C--------------------------------------------------------------------
|
|---|
| 278 | include 'SIZE'
|
|---|
| 279 | include 'SOLN'
|
|---|
| 280 | include 'MASS'
|
|---|
| 281 | include 'INPUT'
|
|---|
| 282 | include 'TSTEP'
|
|---|
| 283 | COMMON /CTOLPR/ DIVEX
|
|---|
| 284 | COMMON /CPRINT/ IFPRINT
|
|---|
| 285 | LOGICAL IFPRINT
|
|---|
| 286 | C
|
|---|
| 287 | COMMON /SCRUZ/ DIVV (LX2,LY2,LZ2,LELV)
|
|---|
| 288 | $ , BDIVV(LX2,LY2,LZ2,LELV)
|
|---|
| 289 | C
|
|---|
| 290 | if (ifsplit) return
|
|---|
| 291 | IF (param(102).eq.0.and.(TOLPDF.NE.0. .OR. ISTEP.LE.5)) return
|
|---|
| 292 | five = 5.0
|
|---|
| 293 | if (param(102).ne.0.0) five=param(102)
|
|---|
| 294 | C
|
|---|
| 295 | NTOT2 = lx2*ly2*lz2*NELV
|
|---|
| 296 | if (ifield.eq.1) then ! avo: sub arguments?
|
|---|
| 297 | CALL OPDIV (BDIVV,VX,VY,VZ)
|
|---|
| 298 | else
|
|---|
| 299 | CALL OPDIV (BDIVV,BX,BY,BZ)
|
|---|
| 300 | endif
|
|---|
| 301 | CALL COL3 (DIVV,BDIVV,BM2INV,NTOT2)
|
|---|
| 302 | DNORM = SQRT(GLSC2(DIVV,BDIVV,NTOT2)/VOLVM2)
|
|---|
| 303 | C
|
|---|
| 304 | if (nio.eq.0) WRITE (6,*) istep,' DNORM, DIVEX',DNORM,DIVEX
|
|---|
| 305 | C
|
|---|
| 306 | c IF (istep.gt.10.and.DNORM.GT.(1.01*DIVEX).AND.DIVEX.GT.0.) then
|
|---|
| 307 | c if (DNORM.gt.1e-8) then
|
|---|
| 308 | c if (nid.eq.0) WRITE(6,*) 'DNORM-DIVEX div. ... aborting'
|
|---|
| 309 | c call prepost (.true.,' ')
|
|---|
| 310 | c call exitt
|
|---|
| 311 | c else
|
|---|
| 312 | c if (nid.eq.0) WRITE(6,*) 'DNORM-DIVEX div. ... small'
|
|---|
| 313 | c endif
|
|---|
| 314 | c endif
|
|---|
| 315 | C
|
|---|
| 316 | c IF (DNORM.GT.(1.2*DIVEX).AND.DIVEX.GT.0.) TOLPDF = 5.*DNORM
|
|---|
| 317 | IF (istep.gt.5.and.tolpdf.eq.0.0.and.
|
|---|
| 318 | $ DNORM.GT.(1.2*DIVEX).AND.DIVEX.GT.0.)
|
|---|
| 319 | $ TOLPDF = FIVE*DNORM
|
|---|
| 320 | C
|
|---|
| 321 | return
|
|---|
| 322 | end
|
|---|
| 323 | FUNCTION VLSC3(X,Y,B,N)
|
|---|
| 324 | C
|
|---|
| 325 | C local inner product, with weight
|
|---|
| 326 | C
|
|---|
| 327 | DIMENSION X(1),Y(1),B(1)
|
|---|
| 328 | REAL DT
|
|---|
| 329 | C
|
|---|
| 330 | include 'OPCTR'
|
|---|
| 331 | C
|
|---|
| 332 | #ifdef TIMER
|
|---|
| 333 | C
|
|---|
| 334 | if (isclld.eq.0) then
|
|---|
| 335 | isclld=1
|
|---|
| 336 | nrout=nrout+1
|
|---|
| 337 | myrout=nrout
|
|---|
| 338 | rname(myrout) = 'VLSC3 '
|
|---|
| 339 | endif
|
|---|
| 340 | isbcnt = 3*n
|
|---|
| 341 | dct(myrout) = dct(myrout) + dfloat(isbcnt)
|
|---|
| 342 | ncall(myrout) = ncall(myrout) + 1
|
|---|
| 343 | dcount = dcount + dfloat(isbcnt)
|
|---|
| 344 | #endif
|
|---|
| 345 | C
|
|---|
| 346 | DT = 0.0
|
|---|
| 347 | DO 10 I=1,N
|
|---|
| 348 | T = X(I)*Y(I)*B(I)
|
|---|
| 349 | DT = DT+T
|
|---|
| 350 | 10 CONTINUE
|
|---|
| 351 | T=DT
|
|---|
| 352 | VLSC3 = T
|
|---|
| 353 | return
|
|---|
| 354 | end
|
|---|
| 355 | c-----------------------------------------------------------------------
|
|---|
| 356 | subroutine updrhse(p,h1,h2,h2inv,ierr)
|
|---|
| 357 | C
|
|---|
| 358 | C Update rhs's if E-matrix has changed
|
|---|
| 359 | C
|
|---|
| 360 | C
|
|---|
| 361 | include 'SIZE'
|
|---|
| 362 | include 'INPUT'
|
|---|
| 363 | include 'MASS'
|
|---|
| 364 | include 'TSTEP'
|
|---|
| 365 | C
|
|---|
| 366 | PARAMETER (LTOT2=LX2*LY2*LZ2*LELV)
|
|---|
| 367 | COMMON /ORTHOV/ RHS(LTOT2,MXPREV)
|
|---|
| 368 | COMMON /ORTHOX/ Pbar(LTOT2),Pnew(LTOT2)
|
|---|
| 369 | COMMON /ORTHOS/ ALPHA(Mxprev), WORK(Mxprev), ALPHAN, DTLAST
|
|---|
| 370 | COMMON /ORTHOI/ Nprev,Mprev
|
|---|
| 371 | COMMON /ORTHOL/ IFNEWE
|
|---|
| 372 | REAL ALPHA,WORK
|
|---|
| 373 | LOGICAL IFNEWE
|
|---|
| 374 | C
|
|---|
| 375 | C
|
|---|
| 376 | REAL P (LX2,LY2,LZ2,LELV)
|
|---|
| 377 | REAL H1 (LX1,LY1,LZ1,LELV)
|
|---|
| 378 | REAL H2 (LX1,LY1,LZ1,LELV)
|
|---|
| 379 | REAL H2INV(LX1,LY1,LZ1,LELV)
|
|---|
| 380 | C
|
|---|
| 381 | integer icalld
|
|---|
| 382 | save icalld
|
|---|
| 383 | data icalld/0/
|
|---|
| 384 |
|
|---|
| 385 | ntot2=lx2*ly2*lz2*nelv
|
|---|
| 386 |
|
|---|
| 387 |
|
|---|
| 388 | C First, we have to decide if the E matrix has changed.
|
|---|
| 389 |
|
|---|
| 390 |
|
|---|
| 391 | if (icalld.eq.0) then
|
|---|
| 392 | icalld=1
|
|---|
| 393 | dtlast=dt
|
|---|
| 394 | endif
|
|---|
| 395 |
|
|---|
| 396 | ifnewe=.false.
|
|---|
| 397 | if (ifmvbd) then
|
|---|
| 398 | ifnewe=.true.
|
|---|
| 399 | call invers2(bm2inv,bm2,ntot2)
|
|---|
| 400 | endif
|
|---|
| 401 | if (dtlast.ne.dt) then
|
|---|
| 402 | ifnewe=.true.
|
|---|
| 403 | dtlast=dt
|
|---|
| 404 | endif
|
|---|
| 405 | c if (ifnewe.and.nio.eq.0) write(6,*) istep,'reorthogo:',nprev
|
|---|
| 406 |
|
|---|
| 407 |
|
|---|
| 408 | C
|
|---|
| 409 | C Next, we reconstruct a new rhs set.
|
|---|
| 410 | C
|
|---|
| 411 | if (ifnewe) then
|
|---|
| 412 | c
|
|---|
| 413 | c new idea...
|
|---|
| 414 | c if (nprev.gt.0) nprev=1
|
|---|
| 415 | c call copy(rhs,pnew,ntot2)
|
|---|
| 416 | c
|
|---|
| 417 | Nprevt = Nprev
|
|---|
| 418 | DO 100 Iprev=1,Nprevt
|
|---|
| 419 | C Orthogonalize this rhs w.r.t. previous rhs's
|
|---|
| 420 | call econj (iprev,h1,h2,h2inv,ierr)
|
|---|
| 421 | if (ierr.eq.1) then
|
|---|
| 422 | if (nio.eq.0) write(6,*) istep,ierr,' econj error'
|
|---|
| 423 | nprev = 0
|
|---|
| 424 | return
|
|---|
| 425 | endif
|
|---|
| 426 | 100 CONTINUE
|
|---|
| 427 | C
|
|---|
| 428 | endif
|
|---|
| 429 | C
|
|---|
| 430 | return
|
|---|
| 431 | end
|
|---|
| 432 | c-----------------------------------------------------------------------
|
|---|
| 433 | subroutine hconj(approx,k,h1,h2,vml,vmk,ws,name4,ierr)
|
|---|
| 434 | c
|
|---|
| 435 | c Orthonormalize the kth vector against vector set
|
|---|
| 436 | c
|
|---|
| 437 | include 'SIZE'
|
|---|
| 438 | include 'INPUT'
|
|---|
| 439 | include 'MASS'
|
|---|
| 440 | include 'SOLN'
|
|---|
| 441 | include 'PARALLEL'
|
|---|
| 442 | include 'TSTEP'
|
|---|
| 443 | c
|
|---|
| 444 | parameter (lt=lx1*ly1*lz1*lelt)
|
|---|
| 445 | real approx(lt,0:1),h1(1),h2(1),vml(1),vmk(1),ws(1)
|
|---|
| 446 | character*4 name4
|
|---|
| 447 | c
|
|---|
| 448 | ierr=0
|
|---|
| 449 | ntot=lx1*ly1*lz1*nelfld(ifield)
|
|---|
| 450 | c
|
|---|
| 451 | call axhelm (approx(1,0),approx(1,k),h1,h2,1,1)
|
|---|
| 452 | call col2 (approx(1,0),vmk,ntot)
|
|---|
| 453 | c
|
|---|
| 454 | c Compute part of the norm (Note: a(0) already scaled by vml)
|
|---|
| 455 | c
|
|---|
| 456 | alpha = glsc2(approx(1,0),approx(1,k),ntot)
|
|---|
| 457 | alph1 = alpha
|
|---|
| 458 | c
|
|---|
| 459 | c Gram-Schmidt
|
|---|
| 460 | c
|
|---|
| 461 | km1=k-1
|
|---|
| 462 | do i=1,km1
|
|---|
| 463 | ws(i) = vlsc2(approx(1,0),approx(1,i),ntot)
|
|---|
| 464 | enddo
|
|---|
| 465 | if (km1.gt.0) call gop(ws,ws(k),'+ ',km1)
|
|---|
| 466 | c
|
|---|
| 467 | do i=1,km1
|
|---|
| 468 | alpham = -ws(i)
|
|---|
| 469 | call add2s2(approx(1,k),approx(1,i),alpham,ntot)
|
|---|
| 470 | alpha = alpha - ws(i)**2
|
|---|
| 471 | enddo
|
|---|
| 472 | c
|
|---|
| 473 | c .Normalize new element in approximation space
|
|---|
| 474 | c
|
|---|
| 475 | eps = 1.e-7
|
|---|
| 476 | if (wdsize.eq.8) eps = 1.e-15
|
|---|
| 477 | ratio = alpha/alph1
|
|---|
| 478 | c
|
|---|
| 479 | if (ratio.le.0) then
|
|---|
| 480 | ierr=1
|
|---|
| 481 | if (nio.eq.0) write(6,12) istep,name4,k,alpha,alph1
|
|---|
| 482 | 12 format(I6,1x,a4,' alpha b4 sqrt:',i4,1p2e12.4)
|
|---|
| 483 | elseif (ratio.le.eps) then
|
|---|
| 484 | ierr=2
|
|---|
| 485 | if (nio.eq.0) write(6,12) istep,name4,k,alpha,alph1
|
|---|
| 486 | else
|
|---|
| 487 | ierr=0
|
|---|
| 488 | alpha = 1.0/sqrt(alpha)
|
|---|
| 489 | call cmult(approx(1,k),alpha,ntot)
|
|---|
| 490 | endif
|
|---|
| 491 | c
|
|---|
| 492 | if (ierr.ne.0) then
|
|---|
| 493 | call axhelm (approx(1,0),approx(1,k),h1,h2,1,1)
|
|---|
| 494 | call col2 (approx(1,0),vmk,ntot)
|
|---|
| 495 | c
|
|---|
| 496 | c Compute part of the norm (Note: a(0) already scaled by vml)
|
|---|
| 497 | c
|
|---|
| 498 | alpha = glsc2(approx(1,0),approx(1,k),ntot)
|
|---|
| 499 | if (nio.eq.0) write(6,12) istep,name4,k,alpha,alph1
|
|---|
| 500 | if (alpha.le.0) then
|
|---|
| 501 | ierr=3
|
|---|
| 502 | if (nio.eq.0) write(6,12) istep,name4,k,alpha,alph1
|
|---|
| 503 | return
|
|---|
| 504 | endif
|
|---|
| 505 | alpha = 1.0/sqrt(alpha)
|
|---|
| 506 | call cmult(approx(1,k),alpha,ntot)
|
|---|
| 507 | ierr = 0
|
|---|
| 508 | endif
|
|---|
| 509 | c
|
|---|
| 510 | return
|
|---|
| 511 | end
|
|---|
| 512 | c-----------------------------------------------------------------------
|
|---|
| 513 | subroutine hmhzpf(name,u,r,h1,h2,mask,mult,imesh,tli,maxit,isd,bi)
|
|---|
| 514 | include 'SIZE'
|
|---|
| 515 | include 'INPUT'
|
|---|
| 516 | include 'MASS'
|
|---|
| 517 | include 'FDMH1'
|
|---|
| 518 | include 'CTIMER'
|
|---|
| 519 | c
|
|---|
| 520 | CHARACTER*4 NAME
|
|---|
| 521 | REAL U (LX1,LY1,LZ1,1)
|
|---|
| 522 | REAL R (LX1,LY1,LZ1,1)
|
|---|
| 523 | REAL H1 (LX1,LY1,LZ1,1)
|
|---|
| 524 | REAL H2 (LX1,LY1,LZ1,1)
|
|---|
| 525 | REAL MASK (LX1,LY1,LZ1,1)
|
|---|
| 526 | REAL MULT (LX1,LY1,LZ1,1)
|
|---|
| 527 | REAL bi (LX1,LY1,LZ1,1)
|
|---|
| 528 | COMMON /CTMP0/ W1 (LX1,LY1,LZ1,LELT)
|
|---|
| 529 | $ , W2 (LX1,LY1,LZ1,LELT)
|
|---|
| 530 | c
|
|---|
| 531 | etime1=dnekclock()
|
|---|
| 532 | c
|
|---|
| 533 | IF (IMESH.EQ.1) NTOT = lx1*ly1*lz1*NELV
|
|---|
| 534 | IF (IMESH.EQ.2) NTOT = lx1*ly1*lz1*NELT
|
|---|
| 535 | c
|
|---|
| 536 | tol = tli
|
|---|
| 537 | if (param(22).ne.0) tol = abs(param(22))
|
|---|
| 538 | CALL CHKTCG1 (TOL,R,H1,H2,MASK,MULT,IMESH,ISD)
|
|---|
| 539 | c
|
|---|
| 540 | c
|
|---|
| 541 | c Set flags for overlapping Schwarz preconditioner (pff 11/12/98)
|
|---|
| 542 | c
|
|---|
| 543 | kfldfdm = -1
|
|---|
| 544 | c if (name.eq.'TEMP') kfldfdm = 0
|
|---|
| 545 | c if (name.eq.'VELX') kfldfdm = 1
|
|---|
| 546 | c if (name.eq.'VELY') kfldfdm = 2
|
|---|
| 547 | c if (name.eq.'VELZ') kfldfdm = 3
|
|---|
| 548 | if (name.eq.'PRES') kfldfdm = ldim+1
|
|---|
| 549 |
|
|---|
| 550 | if (ifdg) then
|
|---|
| 551 | call cggo_dg (u,r,h1,h2,bi,mask,name,tol,maxit)
|
|---|
| 552 | else
|
|---|
| 553 | call cggo
|
|---|
| 554 | $ (u,r,h1,h2,mask,mult,imesh,tol,maxit,isd,bi,name)
|
|---|
| 555 | endif
|
|---|
| 556 | thmhz=thmhz+(dnekclock()-etime1)
|
|---|
| 557 | c
|
|---|
| 558 | c
|
|---|
| 559 | return
|
|---|
| 560 | end
|
|---|
| 561 | c-----------------------------------------------------------------------
|
|---|
| 562 | subroutine hsolve(name,u,r,h1,h2,vmk,vml,imsh,tol,maxit,isd
|
|---|
| 563 | $ ,approx,napprox,bi)
|
|---|
| 564 | c
|
|---|
| 565 | c Either std. Helmholtz solve, or a projection + Helmholtz solve
|
|---|
| 566 | c
|
|---|
| 567 | include 'SIZE'
|
|---|
| 568 | include 'TOTAL'
|
|---|
| 569 | include 'CTIMER'
|
|---|
| 570 | c
|
|---|
| 571 | CHARACTER*4 NAME
|
|---|
| 572 | REAL U (LX1,LY1,LZ1,1)
|
|---|
| 573 | REAL R (LX1,LY1,LZ1,1)
|
|---|
| 574 | REAL H1 (LX1,LY1,LZ1,1)
|
|---|
| 575 | REAL H2 (LX1,LY1,LZ1,1)
|
|---|
| 576 | REAL vmk (LX1,LY1,LZ1,1)
|
|---|
| 577 | REAL vml (LX1,LY1,LZ1,1)
|
|---|
| 578 | REAL bi (LX1,LY1,LZ1,1)
|
|---|
| 579 | REAL approx (1)
|
|---|
| 580 | integer napprox(1)
|
|---|
| 581 | common /ctmp2/ w1 (lx1,ly1,lz1,lelt)
|
|---|
| 582 | common /ctmp3/ w2 (2+2*mxprev)
|
|---|
| 583 |
|
|---|
| 584 | logical ifstdh
|
|---|
| 585 | character*4 cname
|
|---|
| 586 | character*6 name6
|
|---|
| 587 |
|
|---|
| 588 | logical ifwt,ifvec
|
|---|
| 589 |
|
|---|
| 590 | call chcopy(cname,name,4)
|
|---|
| 591 | call capit (cname,4)
|
|---|
| 592 |
|
|---|
| 593 | ifstdh = .true. ! default is no projection
|
|---|
| 594 | if (ifprojfld(ifield)) then
|
|---|
| 595 | ifstdh = .false.
|
|---|
| 596 | endif
|
|---|
| 597 |
|
|---|
| 598 | p945 = param(94)
|
|---|
| 599 | if (cname.eq.'PRES') then
|
|---|
| 600 | ifstdh = .false.
|
|---|
| 601 | p945 = param(95)
|
|---|
| 602 | endif
|
|---|
| 603 |
|
|---|
| 604 | if (ifield.gt.ldimt_proj+1) ifstdh = .true.
|
|---|
| 605 | if (param(93).eq.0) ifstdh = .true.
|
|---|
| 606 | if (p945.eq.0) ifstdh = .true.
|
|---|
| 607 | if (istep.lt.p945) ifstdh = .true.
|
|---|
| 608 |
|
|---|
| 609 | if (ifstdh) then
|
|---|
| 610 | call hmholtz(name,u,r,h1,h2,vmk,vml,imsh,tol,maxit,isd)
|
|---|
| 611 | else
|
|---|
| 612 |
|
|---|
| 613 | n = lx1*ly1*lz1*nelfld(ifield)
|
|---|
| 614 |
|
|---|
| 615 | call col2 (r,vmk,n)
|
|---|
| 616 | call dssum (r,lx1,ly1,lz1)
|
|---|
| 617 |
|
|---|
| 618 | call blank (name6,6)
|
|---|
| 619 | call chcopy(name6,name,4)
|
|---|
| 620 | ifwt = .true.
|
|---|
| 621 | ifvec = .false.
|
|---|
| 622 |
|
|---|
| 623 | call project1
|
|---|
| 624 | $ (r,n,approx,napprox,h1,h2,vmk,vml,ifwt,ifvec,name6)
|
|---|
| 625 |
|
|---|
| 626 | call hmhzpf (name,u,r,h1,h2,vmk,vml,imsh,tol,maxit,isd,bi)
|
|---|
| 627 |
|
|---|
| 628 | call project2
|
|---|
| 629 | $ (u,n,approx,napprox,h1,h2,vmk,vml,ifwt,ifvec,name6)
|
|---|
| 630 |
|
|---|
| 631 | endif
|
|---|
| 632 |
|
|---|
| 633 | return
|
|---|
| 634 | end
|
|---|
| 635 | c-----------------------------------------------------------------------
|
|---|
| 636 | subroutine project1(b,n,rvar,ivar,h1,h2,msk,w,ifwt,ifvec,name6)
|
|---|
| 637 |
|
|---|
| 638 | c 1. Compute the projection of x onto X
|
|---|
| 639 |
|
|---|
| 640 | c 2. Re-orthogonalize the X basis set and corresponding B=A*X
|
|---|
| 641 | c vectors if A has changed.
|
|---|
| 642 |
|
|---|
| 643 | c Output: b = b - projection of b onto B
|
|---|
| 644 |
|
|---|
| 645 | c Input: n = length of field (or multifields, when ifvec=true)
|
|---|
| 646 | c rvar = real array of field values, including old h1,h2, etc.
|
|---|
| 647 | c ivar = integer array of pointers, etc.
|
|---|
| 648 | c h1 = current h1, for Axhelm(.,.,h1,h2,...)
|
|---|
| 649 | c h2 = current h2
|
|---|
| 650 | c msk = mask for Dirichlet BCs
|
|---|
| 651 | c w = weight for inner products (typ. w=vmult, tmult, etc.)
|
|---|
| 652 | c ifwt = use weighted inner products when ifwt=.true.
|
|---|
| 653 | c ifvec = are x and b vectors, or scalar fields?
|
|---|
| 654 | c name6 = discriminator for action of A*x
|
|---|
| 655 |
|
|---|
| 656 | c The idea here is to have one pair of projection routines for
|
|---|
| 657 | c constructing the new rhs (project1) and reconstructing the new
|
|---|
| 658 | c solution (x = xbar + dx) plus updating the approximation space.
|
|---|
| 659 | c The latter functions are done in project2.
|
|---|
| 660 | c
|
|---|
| 661 | c The approximation space X and corresponding right-hand sides,
|
|---|
| 662 | c B := A*X are stored in rvar, as well as h1old and h2old and a
|
|---|
| 663 | c couple of other auxiliary arrays.
|
|---|
| 664 |
|
|---|
| 665 | c In this new code, we retain both X and B=A*X and we re-orthogonalize
|
|---|
| 666 | c at each timestep (with no extra matrix-vector products, but O(nm)
|
|---|
| 667 | c work. The idea is to retain fresh vectors by injecting the most
|
|---|
| 668 | c recent solution and pushing the oldest off the stack, hopefully
|
|---|
| 669 | c keeping the number of vectors, m, small.
|
|---|
| 670 |
|
|---|
| 671 |
|
|---|
| 672 | include 'SIZE' ! For nid/nio
|
|---|
| 673 | include 'TSTEP' ! For istep
|
|---|
| 674 | include 'CTIMER'
|
|---|
| 675 |
|
|---|
| 676 | real b(n),rvar(n,1),h1(n),h2(n),w(n),msk(n)
|
|---|
| 677 | integer ivar(1)
|
|---|
| 678 | character*6 name6
|
|---|
| 679 | logical ifwt,ifvec
|
|---|
| 680 |
|
|---|
| 681 | etime0 = dnekclock()
|
|---|
| 682 |
|
|---|
| 683 | nn = n
|
|---|
| 684 | if (ifvec) nn = n*ldim
|
|---|
| 685 |
|
|---|
| 686 | call proj_get_ivar
|
|---|
| 687 | $ (m,mmx,ixb,ibb,ix,ib,ih1,ih2,ivar,n,ifvec,name6)
|
|---|
| 688 |
|
|---|
| 689 | if (m.le.0) return
|
|---|
| 690 |
|
|---|
| 691 | ireset=iproj_chk(rvar(ih1,1),rvar(ih2,1),h1,h2,n) ! Updated matrix?
|
|---|
| 692 |
|
|---|
| 693 | bb4 = glsc3(b,w,b,n)
|
|---|
| 694 | bb4 = sqrt(bb4)
|
|---|
| 695 |
|
|---|
| 696 |
|
|---|
| 697 | c Re-orthogonalize basis set w.r.t. new vectors if space has changed.
|
|---|
| 698 |
|
|---|
| 699 | if (ireset.eq.1) then
|
|---|
| 700 |
|
|---|
| 701 | do j=0,m-1 ! First, set B := A*X
|
|---|
| 702 | jb = ib+j*nn !Iterate through xs and bs
|
|---|
| 703 | jx = ix+j*nn
|
|---|
| 704 | call proj_matvec (rvar(jb,1),rvar(jx,1),n,h1,h2,msk,name6)
|
|---|
| 705 | enddo
|
|---|
| 706 |
|
|---|
| 707 | if (nio.eq.0 .and. loglevel.gt.2)
|
|---|
| 708 | $ write(6,'(13x,A)') 'Reorthogonalize Basis'
|
|---|
| 709 |
|
|---|
| 710 | call proj_ortho_full
|
|---|
| 711 | $ (rvar(ix,1),rvar(ib,1),n,m,w,ifwt,ifvec,name6)
|
|---|
| 712 |
|
|---|
| 713 | ivar(2) = m ! Update number of saved vectors
|
|---|
| 714 |
|
|---|
| 715 | endif
|
|---|
| 716 |
|
|---|
| 717 | c ixb is pointer to xbar, ibb is pointer to bbar := A*xbar
|
|---|
| 718 |
|
|---|
| 719 | call project1_a(rvar(ixb,1),rvar(ibb,1),b,rvar(ix,1),rvar(ib,1)
|
|---|
| 720 | $ ,n,m,w,ifwt,ifvec)
|
|---|
| 721 |
|
|---|
| 722 | baf = glsc3(b,w,b,n)
|
|---|
| 723 | baf = sqrt(baf)
|
|---|
| 724 | ratio = bb4/baf
|
|---|
| 725 |
|
|---|
| 726 | tproj = tproj + dnekclock() - etime0
|
|---|
| 727 |
|
|---|
| 728 | if (nio.eq.0) write(6,1) istep,' Project ' // name6,
|
|---|
| 729 | & baf,bb4,ratio,m,mmx
|
|---|
| 730 | 1 format(i11,a,6x,1p3e13.4,i4,i4)
|
|---|
| 731 |
|
|---|
| 732 | return
|
|---|
| 733 | end
|
|---|
| 734 | c-----------------------------------------------------------------------
|
|---|
| 735 | subroutine project1_a(xbar,bbar,b,xx,bb,n,m,w,ifwt,ifvec)
|
|---|
| 736 |
|
|---|
| 737 | c xbar is best fit in xx, bbar = A*xbar
|
|---|
| 738 | c b <-- b - bbar
|
|---|
| 739 |
|
|---|
| 740 | include 'SIZE'
|
|---|
| 741 | include 'TSTEP'
|
|---|
| 742 | include 'INPUT'
|
|---|
| 743 | include 'PARALLEL'
|
|---|
| 744 | parameter(lt=lx1*ly1*lz1*lelv)
|
|---|
| 745 | real xbar(n),bbar(n),b(n),xx(n,m),bb(n,m),w(n)
|
|---|
| 746 | logical ifwt,ifvec
|
|---|
| 747 | integer e,eg
|
|---|
| 748 | parameter(lxyz=lx1*ly1*lz1)
|
|---|
| 749 | real tt(lxyz,lelt)
|
|---|
| 750 |
|
|---|
| 751 | real work(mxprev),alpha(mxprev)
|
|---|
| 752 |
|
|---|
| 753 | if (m.le.0) return
|
|---|
| 754 |
|
|---|
| 755 | !First round of CGS
|
|---|
| 756 | do k = 1, m
|
|---|
| 757 | if(ifwt) then
|
|---|
| 758 | alpha(k) = vlsc3(xx(1,k),w,b,n)
|
|---|
| 759 | else
|
|---|
| 760 | alpha(k) = vlsc2(xx(1,k),b,n)
|
|---|
| 761 | endif
|
|---|
| 762 | enddo
|
|---|
| 763 | !First one outside loop to avoid zeroing xbar and bbar
|
|---|
| 764 | call gop(alpha,work,'+ ',m)
|
|---|
| 765 | call cmult2(xbar,xx(1,1),alpha(1),n)
|
|---|
| 766 | call cmult2(bbar,bb(1,1),alpha(1),n)
|
|---|
| 767 | call add2s2(b,bb(1,1),-alpha(1),n)
|
|---|
| 768 | do k = 2,m
|
|---|
| 769 | call add2s2(xbar,xx(1,k),alpha(k),n)
|
|---|
| 770 | call add2s2(bbar,bb(1,k),alpha(k),n)
|
|---|
| 771 | call add2s2(b,bb(1,k),-alpha(k),n)
|
|---|
| 772 | enddo
|
|---|
| 773 | !Second round of CGS
|
|---|
| 774 | do k = 1, m
|
|---|
| 775 | if(ifwt) then
|
|---|
| 776 | alpha(k) = vlsc3(xx(1,k),w,b,n)
|
|---|
| 777 | else
|
|---|
| 778 | alpha(k) = vlsc2(xx(1,k),b,n)
|
|---|
| 779 | endif
|
|---|
| 780 | enddo
|
|---|
| 781 | call gop(alpha,work,'+ ',m)
|
|---|
| 782 | do k = 1,m
|
|---|
| 783 | call add2s2(xbar,xx(1,k),alpha(k),n)
|
|---|
| 784 | call add2s2(bbar,bb(1,k),alpha(k),n)
|
|---|
| 785 | call add2s2(b,bb(1,k),-alpha(k),n)
|
|---|
| 786 | enddo
|
|---|
| 787 |
|
|---|
| 788 | return
|
|---|
| 789 | end
|
|---|
| 790 | c-----------------------------------------------------------------------
|
|---|
| 791 | function iproj_chk(h1old,h2old,h1,h2,n)
|
|---|
| 792 | include 'SIZE'
|
|---|
| 793 | include 'TOTAL'
|
|---|
| 794 |
|
|---|
| 795 | c Matrix has changed if h1/h2 differ from old values
|
|---|
| 796 |
|
|---|
| 797 | real h1(n),h2(n),h1old(n),h2old(n)
|
|---|
| 798 |
|
|---|
| 799 | iproj_chk = 0
|
|---|
| 800 |
|
|---|
| 801 | if (ifmvbd) then
|
|---|
| 802 | iproj_chk = 1
|
|---|
| 803 | return
|
|---|
| 804 | endif
|
|---|
| 805 |
|
|---|
| 806 | dh1 = 0.
|
|---|
| 807 | dh2 = 0.
|
|---|
| 808 | do i=1,n
|
|---|
| 809 | dh1 = max(dh1,abs(h1(i)-h1old(i)))
|
|---|
| 810 | dh2 = max(dh2,abs(h2(i)-h2old(i)))
|
|---|
| 811 | enddo
|
|---|
| 812 | dh = max(dh1,dh2)
|
|---|
| 813 | dh = glmax(dh,1) ! Max across all processors
|
|---|
| 814 |
|
|---|
| 815 | if (dh.gt.0) then
|
|---|
| 816 |
|
|---|
| 817 | call copy(h1old,h1,n) ! Save old h1 / h2 values
|
|---|
| 818 | call copy(h2old,h2,n)
|
|---|
| 819 |
|
|---|
| 820 | iproj_chk = 1 ! Force re-orthogonalization of basis
|
|---|
| 821 |
|
|---|
| 822 | endif
|
|---|
| 823 |
|
|---|
| 824 | return
|
|---|
| 825 | end
|
|---|
| 826 | c-----------------------------------------------------------------------
|
|---|
| 827 | subroutine proj_matvec(b,x,n,h1,h2,msk,name6)
|
|---|
| 828 | include 'SIZE'
|
|---|
| 829 | include 'TOTAL'
|
|---|
| 830 | real b(n),x(n),h1(n),h2(n),msk(n)
|
|---|
| 831 | character*6 name6
|
|---|
| 832 |
|
|---|
| 833 | c This is the default matvec for nekcem.
|
|---|
| 834 |
|
|---|
| 835 | c The code can later be updated to support different matvec
|
|---|
| 836 | c implementations, which would be discriminated by the character
|
|---|
| 837 | c string "name6"
|
|---|
| 838 |
|
|---|
| 839 | isd = 1 ! This probably won't work for axisymmetric
|
|---|
| 840 | imsh = 1
|
|---|
| 841 | if (iftmsh(ifield)) imsh=2
|
|---|
| 842 |
|
|---|
| 843 | call axhelm (b,x,h1,h2,imsh,isd) ! b = A x
|
|---|
| 844 | call dssum (b,lx1,ly1,lz1)
|
|---|
| 845 | call col2 (b,msk,n)
|
|---|
| 846 |
|
|---|
| 847 | return
|
|---|
| 848 | end
|
|---|
| 849 | c-----------------------------------------------------------------------
|
|---|
| 850 | !O(nm) method for updating projection space
|
|---|
| 851 | !See James Lotte's note or Nicholas Christensen's master's thesis
|
|---|
| 852 | subroutine proj_ortho(xx,bb,n,m,w,ifwt,ifvec,name6)
|
|---|
| 853 |
|
|---|
| 854 | include 'SIZE' ! nio
|
|---|
| 855 | include 'TSTEP' ! istep
|
|---|
| 856 | include 'PARALLEL' ! wdsize
|
|---|
| 857 |
|
|---|
| 858 | real xx(n,1), bb(n,1), w(n)
|
|---|
| 859 | character*6 name6
|
|---|
| 860 | logical ifwt, ifvec
|
|---|
| 861 | real tol, nrm, scl1, scl2, c, s
|
|---|
| 862 | real work(mxprev), alpha(mxprev), beta(mxprev)
|
|---|
| 863 | integer h
|
|---|
| 864 |
|
|---|
| 865 | if(m.le.0) return !No vectors to ortho-normalize
|
|---|
| 866 |
|
|---|
| 867 | ! AX = B
|
|---|
| 868 | ! Calculate dx, db: dx = x-XX^Tb, db=b-BX^Tb
|
|---|
| 869 |
|
|---|
| 870 | do k = 1, m !First round CGS
|
|---|
| 871 | if(ifwt) then
|
|---|
| 872 | alpha(k) = 0.5*(vlsc3(xx(1,k),w,bb(1,m),n)
|
|---|
| 873 | $ + vlsc3(bb(1,k),w,xx(1,m),n))
|
|---|
| 874 | else
|
|---|
| 875 | alpha(k) = 0.5*(vlsc2(xx(1,k),bb(1,m),n)
|
|---|
| 876 | $ + vlsc2(bb(1,k),xx(1,m),n))
|
|---|
| 877 | endif
|
|---|
| 878 | enddo
|
|---|
| 879 | call gop(alpha,work,'+ ',m)
|
|---|
| 880 | nrm = sqrt(alpha(m)) !Calculate A-norm of new vector
|
|---|
| 881 | do k = 1,m-1
|
|---|
| 882 | call add2s2(xx(1,m),xx(1,k),-alpha(k),n)
|
|---|
| 883 | call add2s2(bb(1,m),bb(1,k),-alpha(k),n)
|
|---|
| 884 | enddo
|
|---|
| 885 |
|
|---|
| 886 | do k = 1, m-1 !Second round CGS
|
|---|
| 887 | if(ifwt) then
|
|---|
| 888 | beta(k) = 0.5*(vlsc3(xx(1,k),w,bb(1,m),n)
|
|---|
| 889 | $ + vlsc3(bb(1,k),w,xx(1,m),n))
|
|---|
| 890 | else
|
|---|
| 891 | beta(k) = 0.5*(vlsc2(xx(1,k),bb(1,m),n)
|
|---|
| 892 | $ + vlsc2(bb(1,k),xx(1,m),n))
|
|---|
| 893 | endif
|
|---|
| 894 | enddo
|
|---|
| 895 | call gop(beta,work,'+ ',m-1)
|
|---|
| 896 | do k = 1,m-1
|
|---|
| 897 | call add2s2(xx(1,m),xx(1,k),-beta(k),n)
|
|---|
| 898 | call add2s2(bb(1,m),bb(1,k),-beta(k),n)
|
|---|
| 899 | !While we're at it,
|
|---|
| 900 | !Sum weights from each round to get the total alpha
|
|---|
| 901 | alpha(k) = alpha(k) + beta(k)
|
|---|
| 902 | enddo
|
|---|
| 903 |
|
|---|
| 904 | !Calculate A-norm of newest solution
|
|---|
| 905 | if(ifwt) then
|
|---|
| 906 | alpha(m) = glsc3(xx(1,m), w, bb(1,m), n)
|
|---|
| 907 | else
|
|---|
| 908 | alpha(m) = glsc2(xx(1,m), bb(1,m), n)
|
|---|
| 909 | endif
|
|---|
| 910 | alpha(m) = sqrt(alpha(m))
|
|---|
| 911 | !dx and db now stored in last column of xx and bb
|
|---|
| 912 |
|
|---|
| 913 | c Set tolerance for linear independence
|
|---|
| 914 | tol = 1.e-7
|
|---|
| 915 | if (wdsize.eq.4) tol=1.e-3
|
|---|
| 916 |
|
|---|
| 917 | c Check for linear independence.
|
|---|
| 918 | if(alpha(m).gt.tol*nrm) then !New vector is linearly independent
|
|---|
| 919 |
|
|---|
| 920 | !Normalize dx and db
|
|---|
| 921 | scl1 = 1.0/alpha(m)
|
|---|
| 922 | call cmult(xx(1,m), scl1, n)
|
|---|
| 923 | call cmult(bb(1,m), scl1, n)
|
|---|
| 924 |
|
|---|
| 925 | !We want to throw away the oldest information
|
|---|
| 926 | !The below propagates newest information to first vector.
|
|---|
| 927 | !This will make the first vector a scalar
|
|---|
| 928 | !multiple of x.
|
|---|
| 929 | do k = m, 2, -1
|
|---|
| 930 | h = k - 1
|
|---|
| 931 | call givens_rotation(alpha(h),alpha(k),c,s,nrm)
|
|---|
| 932 | alpha(h) = nrm
|
|---|
| 933 | do i = 1, n !Apply rotation to xx and bb
|
|---|
| 934 | scl1 = c*xx(i,h) + s*xx(i,k)
|
|---|
| 935 | xx(i,k) = -s*xx(i,h) + c*xx(i,k)
|
|---|
| 936 | xx(i,h) = scl1
|
|---|
| 937 | scl2 = c*bb(i,h) + s*bb(i,k)
|
|---|
| 938 | bb(i,k) = -s*bb(i,h) + c*bb(i,k)
|
|---|
| 939 | bb(i,h) = scl2
|
|---|
| 940 | enddo
|
|---|
| 941 | enddo
|
|---|
| 942 |
|
|---|
| 943 | else !New vector is not linearly independent, forget about it
|
|---|
| 944 | k = m !location of rank deficient column
|
|---|
| 945 |
|
|---|
| 946 | if (nio.eq.0 .and. loglevel.gt.2)
|
|---|
| 947 | $ write(6,1) istep,k,m,name6,alpha(m),tol
|
|---|
| 948 |
|
|---|
| 949 | 1 format(i9,'proj_ortho: ',2i4,1x,a6,
|
|---|
| 950 | $ ' Detect rank deficiency:',1p2e12.4)
|
|---|
| 951 |
|
|---|
| 952 | m = m - 1 !Remove column
|
|---|
| 953 | endif
|
|---|
| 954 |
|
|---|
| 955 | return
|
|---|
| 956 | end
|
|---|
| 957 | c-----------------------------------------------------------------------
|
|---|
| 958 | !Function to switch between mgs and cgs2 full reorthogonalization
|
|---|
| 959 | subroutine proj_ortho_full(xx,bb,n,m,w,ifwt,ifvec,name6)
|
|---|
| 960 |
|
|---|
| 961 | include 'SIZE'
|
|---|
| 962 |
|
|---|
| 963 | real xx(n,1),bb(n,1),w(n)
|
|---|
| 964 | character*6 name6
|
|---|
| 965 | logical ifwt,ifvec
|
|---|
| 966 | integer flag(mxprev)
|
|---|
| 967 |
|
|---|
| 968 | !call proj_ortho_full_mgs(xx,bb,n,m,w,ifwt,ifvec,name6)
|
|---|
| 969 | call proj_ortho_full_cgs2(xx,bb,n,m,w,ifwt,ifvec,name6)
|
|---|
| 970 |
|
|---|
| 971 | return
|
|---|
| 972 | end
|
|---|
| 973 | c-----------------------------------------------------------------------
|
|---|
| 974 | c Full MGS reorthogonalization
|
|---|
| 975 | subroutine proj_ortho_full_mgs(xx,bb,n,m,w,ifwt,ifvec,name6)
|
|---|
| 976 |
|
|---|
| 977 | include 'SIZE' ! nio
|
|---|
| 978 | include 'TSTEP' ! istep
|
|---|
| 979 | include 'PARALLEL' ! wdsize
|
|---|
| 980 |
|
|---|
| 981 | real xx(n,1),bb(n,1),w(n)
|
|---|
| 982 | character*6 name6
|
|---|
| 983 | logical ifwt,ifvec
|
|---|
| 984 | integer flag(mxprev)
|
|---|
| 985 | real normk,normp,alpha
|
|---|
| 986 |
|
|---|
| 987 | if (m.le.0) return
|
|---|
| 988 |
|
|---|
| 989 | if ( ifwt) alpha = glsc3(xx(1,m),w,bb(1,m),n)
|
|---|
| 990 | if (.not. ifwt) alpha = glsc2(xx(1,m),bb(1,m),n)
|
|---|
| 991 | if (alpha.eq.0) return
|
|---|
| 992 |
|
|---|
| 993 | scale = 1./sqrt(alpha)
|
|---|
| 994 | call cmult(xx(1,m),scale,n)
|
|---|
| 995 | call cmult(bb(1,m),scale,n)
|
|---|
| 996 | flag(m) = 1
|
|---|
| 997 |
|
|---|
| 998 | do k=m-1,1,-1 ! Reorthogonalize, starting with latest solution
|
|---|
| 999 |
|
|---|
| 1000 | if ( ifwt) normk = glsc3(xx(1,k),w,bb(1,k),n)
|
|---|
| 1001 | if (.not. ifwt) normk = glsc2(xx(1,k),bb(1,k),n)
|
|---|
| 1002 | normk=sqrt(normk)
|
|---|
| 1003 |
|
|---|
| 1004 | do j=m,k+1,-1 ! Modified GS
|
|---|
| 1005 | if(flag(j).eq.1) then
|
|---|
| 1006 | alpha = 0.
|
|---|
| 1007 | if (ifwt) then
|
|---|
| 1008 | alpha = alpha + .5*(vlsc3(xx(1,j),w,bb(1,k),n)
|
|---|
| 1009 | $ + vlsc3(bb(1,j),w,xx(1,k),n))
|
|---|
| 1010 | else
|
|---|
| 1011 | alpha = alpha + .5*(vlsc2(xx(1,j),bb(1,k),n)
|
|---|
| 1012 | $ + vlsc2(bb(1,j),xx(1,k),n))
|
|---|
| 1013 | endif
|
|---|
| 1014 | scale = -glsum(alpha,1)
|
|---|
| 1015 | call add2s2(xx(1,k),xx(1,j),scale,n)
|
|---|
| 1016 | call add2s2(bb(1,k),bb(1,j),scale,n)
|
|---|
| 1017 | endif
|
|---|
| 1018 | enddo
|
|---|
| 1019 | if ( ifwt) normp = glsc3(xx(1,k),w,bb(1,k),n)
|
|---|
| 1020 | if (.not. ifwt) normp = glsc2(xx(1,k),bb(1,k),n)
|
|---|
| 1021 | if (normp.gt.0.0) normp=sqrt(normp)
|
|---|
| 1022 |
|
|---|
| 1023 | tol = 1.e-7
|
|---|
| 1024 | if (wdsize.eq.4) tol=1.e-3
|
|---|
| 1025 |
|
|---|
| 1026 | if (normp.gt.tol*normk) then ! linearly independent vectors
|
|---|
| 1027 | scale = 1./normp
|
|---|
| 1028 | call cmult(xx(1,k),scale,n)
|
|---|
| 1029 | call cmult(bb(1,k),scale,n)
|
|---|
| 1030 | flag(k) = 1
|
|---|
| 1031 | c if (nio.eq.0) write(6,2) istep,k,m,name6,normp,normk
|
|---|
| 1032 | c 2 format(i9,'proj_ortho: ',2i4,1x,a6,' project ok.'
|
|---|
| 1033 | c $ ,1p2e12.4)
|
|---|
| 1034 | else
|
|---|
| 1035 | flag(k) = 0
|
|---|
| 1036 | c if (nio.eq.0) write(6,1) istep,k,m,name6,normp,normk
|
|---|
| 1037 | c 1 format(i9,'proj_ortho: ',2i4,1x,a6,' Detect rank deficiency:'
|
|---|
| 1038 | c $ ,1p2e12.4)
|
|---|
| 1039 | endif
|
|---|
| 1040 |
|
|---|
| 1041 | enddo
|
|---|
| 1042 |
|
|---|
| 1043 | k=0
|
|---|
| 1044 | do j=1,m
|
|---|
| 1045 | if (flag(j).eq.1) then
|
|---|
| 1046 | k=k+1
|
|---|
| 1047 | if (k.lt.j) then
|
|---|
| 1048 | call copy(xx(1,k),xx(1,j),n)
|
|---|
| 1049 | call copy(bb(1,k),bb(1,j),n)
|
|---|
| 1050 | endif
|
|---|
| 1051 | endif
|
|---|
| 1052 | enddo
|
|---|
| 1053 | m = k
|
|---|
| 1054 |
|
|---|
| 1055 | return
|
|---|
| 1056 | end
|
|---|
| 1057 | c-----------------------------------------------------------------------
|
|---|
| 1058 | !CGS2 version of full reorthogonalization, possibly more stable in
|
|---|
| 1059 | !certain instances. Much faster for large m.
|
|---|
| 1060 | subroutine proj_ortho_full_cgs2(xx,bb,n,m,w,ifwt,ifvec,name6)
|
|---|
| 1061 |
|
|---|
| 1062 | include 'SIZE' ! nio
|
|---|
| 1063 | include 'TSTEP' ! istep
|
|---|
| 1064 | include 'PARALLEL' ! wdsize
|
|---|
| 1065 |
|
|---|
| 1066 | real xx(n,1),bb(n,1),w(n)
|
|---|
| 1067 | character*6 name6
|
|---|
| 1068 | logical ifwt,ifvec
|
|---|
| 1069 | integer flag(mxprev)
|
|---|
| 1070 | real normk,normp,alpha(mxprev),work(mxprev),scl1,tol
|
|---|
| 1071 |
|
|---|
| 1072 | if (m.le.0) return
|
|---|
| 1073 |
|
|---|
| 1074 | tol = 1.e-7
|
|---|
| 1075 | if (wdsize.eq.4) tol=1.e-3
|
|---|
| 1076 |
|
|---|
| 1077 | do i = 1, 2 !Do this twice for CGS2
|
|---|
| 1078 |
|
|---|
| 1079 | do k = m, 1, -1
|
|---|
| 1080 | do j = m, k, -1
|
|---|
| 1081 | alpha(j) = 0.0
|
|---|
| 1082 | if(ifwt) then
|
|---|
| 1083 | alpha(j) = .5*(vlsc3(xx(1,j),w,bb(1,k),n)
|
|---|
| 1084 | $ + vlsc3(bb(1,j),w,xx(1,k),n))
|
|---|
| 1085 | else
|
|---|
| 1086 | alpha(j) = .5*(vlsc2(xx(1,j),bb(1,k),n)
|
|---|
| 1087 | $ + vlsc2(bb(1,j),xx(1,k),n))
|
|---|
| 1088 | endif
|
|---|
| 1089 | enddo
|
|---|
| 1090 | call gop(alpha(k), work,'+ ',(m - k) + 1)
|
|---|
| 1091 | do j = m, k+1, -1
|
|---|
| 1092 | call add2s2(xx(1,k),xx(1,j),-alpha(j),n)
|
|---|
| 1093 | call add2s2(bb(1,k),bb(1,j),-alpha(j),n)
|
|---|
| 1094 | enddo
|
|---|
| 1095 | normp = sqrt(alpha(k))
|
|---|
| 1096 | if(ifwt) then
|
|---|
| 1097 | normk = glsc3(xx(1,k),w,bb(1,k),n)
|
|---|
| 1098 | else
|
|---|
| 1099 | normk = glsc2(xx(1,k),bb(1,k),n)
|
|---|
| 1100 | endif
|
|---|
| 1101 | normk = sqrt(normk)
|
|---|
| 1102 | if(normk.gt.tol*normp) then
|
|---|
| 1103 | scl1 = 1.0/normk
|
|---|
| 1104 | call cmult(xx(1,k), scl1, n)
|
|---|
| 1105 | call cmult(bb(1,k), scl1, n)
|
|---|
| 1106 | flag(k) = 1
|
|---|
| 1107 | else
|
|---|
| 1108 | flag(k) = 0
|
|---|
| 1109 | endif
|
|---|
| 1110 | enddo
|
|---|
| 1111 |
|
|---|
| 1112 | enddo
|
|---|
| 1113 |
|
|---|
| 1114 | k=0
|
|---|
| 1115 | do j=1,m
|
|---|
| 1116 | if (flag(j).eq.1) then
|
|---|
| 1117 | k=k+1
|
|---|
| 1118 | if (k.lt.j) then
|
|---|
| 1119 | call copy(xx(1,k),xx(1,j),n)
|
|---|
| 1120 | call copy(bb(1,k),bb(1,j),n)
|
|---|
| 1121 | endif
|
|---|
| 1122 | endif
|
|---|
| 1123 | enddo
|
|---|
| 1124 | m = k
|
|---|
| 1125 |
|
|---|
| 1126 | return
|
|---|
| 1127 | end
|
|---|
| 1128 | c-----------------------------------------------------------------------
|
|---|
| 1129 | subroutine project2(x,n,rvar,ivar,h1,h2,msk,w,ifwt,ifvec,name6)
|
|---|
| 1130 |
|
|---|
| 1131 | include 'CTIMER'
|
|---|
| 1132 |
|
|---|
| 1133 | real x(n),b(n),rvar(n,1),h1(n),h2(n),w(n),msk(n)
|
|---|
| 1134 | integer ivar(1)
|
|---|
| 1135 | character*6 name6
|
|---|
| 1136 | logical ifwt,ifvec
|
|---|
| 1137 |
|
|---|
| 1138 | etime0 = dnekclock()
|
|---|
| 1139 |
|
|---|
| 1140 | call proj_get_ivar(m,mmx,ixb,ibb,ix,ib,ih1,ih2,
|
|---|
| 1141 | $ ivar,n,ifvec,name6)
|
|---|
| 1142 |
|
|---|
| 1143 | c ix is pointer to X, ib is pointer to B
|
|---|
| 1144 | c ixb is pointer to xbar, ibb is pointer to bbar := A*xbar
|
|---|
| 1145 |
|
|---|
| 1146 | call project2_a(x,rvar(ixb,1),rvar(ix,1),rvar(ib,1)
|
|---|
| 1147 | $ ,n,m,mmx,h1,h2,msk,w,ifwt,ifvec,name6)
|
|---|
| 1148 |
|
|---|
| 1149 | ivar(2) = m ! Update number of saved vectors
|
|---|
| 1150 |
|
|---|
| 1151 | tproj = tproj + dnekclock() - etime0
|
|---|
| 1152 |
|
|---|
| 1153 | return
|
|---|
| 1154 | end
|
|---|
| 1155 | c-----------------------------------------------------------------------
|
|---|
| 1156 | subroutine project2_a(x,xbar,xx,bb,n,m,mmx,h1,h2,msk,w,ifwt,ifvec,
|
|---|
| 1157 | $ name6)
|
|---|
| 1158 |
|
|---|
| 1159 | include 'SIZE'
|
|---|
| 1160 |
|
|---|
| 1161 | real x(n),xbar(n),xx(n,1),bb(n,1),h1(n),h2(n),w(n),msk(n)
|
|---|
| 1162 | character*6 name6
|
|---|
| 1163 | logical ifwt,ifvec
|
|---|
| 1164 |
|
|---|
| 1165 | nn = n
|
|---|
| 1166 | if (ifvec) nn=ldim*n
|
|---|
| 1167 |
|
|---|
| 1168 | if (m.gt.0) call add2(x,xbar,n) ! Restore desired solution
|
|---|
| 1169 |
|
|---|
| 1170 | !Uncomment this if using full reorthogonalization
|
|---|
| 1171 | c if (m.eq.mmx) then ! Push old vector off the stack
|
|---|
| 1172 | c do k=2,mmx
|
|---|
| 1173 | c call copy (xx(1,k-1),xx(1,k),nn)
|
|---|
| 1174 | c call copy (bb(1,k-1),bb(1,k),nn)
|
|---|
| 1175 | c enddo
|
|---|
| 1176 | c endif
|
|---|
| 1177 |
|
|---|
| 1178 | m = min(m+1,mmx)
|
|---|
| 1179 | !print *, "m", m
|
|---|
| 1180 | call copy (xx(1,m),x,nn) ! Update (X,B)
|
|---|
| 1181 | call proj_matvec (bb(1,m),xx(1,m),n,h1,h2,msk,name6)
|
|---|
| 1182 | call proj_ortho (xx,bb,n,m,w,ifwt,ifvec,name6) !Update orthogonalization
|
|---|
| 1183 | !Uncomment the if block above if using full reorthogonalization
|
|---|
| 1184 | c call proj_ortho_full (xx,bb,n,m,w,ifwt,ifvec,name6) !Fully reorthogonalize
|
|---|
| 1185 |
|
|---|
| 1186 | return
|
|---|
| 1187 | end
|
|---|
| 1188 | c-----------------------------------------------------------------------
|
|---|
| 1189 | subroutine proj_get_ivar
|
|---|
| 1190 | $ (m,mmx,ixb,ibb,ix,ib,ih1,ih2,ivar,n,ifvec,name6)
|
|---|
| 1191 |
|
|---|
| 1192 | include 'SIZE'
|
|---|
| 1193 | include 'TSTEP'
|
|---|
| 1194 |
|
|---|
| 1195 | logical ifvec
|
|---|
| 1196 | character*6 name6
|
|---|
| 1197 |
|
|---|
| 1198 | integer ivar(10)
|
|---|
| 1199 |
|
|---|
| 1200 | integer icalld
|
|---|
| 1201 | save icalld
|
|---|
| 1202 | data icalld/0/
|
|---|
| 1203 |
|
|---|
| 1204 | m = ivar(2)
|
|---|
| 1205 | mmx = (mxprev-4)/2 ! ivar=0 --> mxprev array
|
|---|
| 1206 | ivar(1) = mmx
|
|---|
| 1207 |
|
|---|
| 1208 | nn = n
|
|---|
| 1209 | if (ifvec) nn = n*ldim ! Number of entries in a vector
|
|---|
| 1210 |
|
|---|
| 1211 |
|
|---|
| 1212 | ih1 = 1
|
|---|
| 1213 | ih2 = ih1 + n
|
|---|
| 1214 | ixb = ih2 + n ! pointer to xbar
|
|---|
| 1215 | ibb = ixb + nn ! " to bbar
|
|---|
| 1216 | ix = ibb + nn ! " to X
|
|---|
| 1217 | ib = ix + nn*mmx ! " to B
|
|---|
| 1218 |
|
|---|
| 1219 | return
|
|---|
| 1220 | end
|
|---|
| 1221 | c-----------------------------------------------------------------------
|
|---|
| 1222 | subroutine laplacep(name,u,mask,mult,ifld,tol,maxi,approx,napprox)
|
|---|
| 1223 | c
|
|---|
| 1224 | c Solve Laplace's equation, with projection onto previous solutions.
|
|---|
| 1225 | c
|
|---|
| 1226 | c Boundary condition strategy:
|
|---|
| 1227 | c
|
|---|
| 1228 | c u = u0 + ub
|
|---|
| 1229 | c
|
|---|
| 1230 | c u0 = 0 on Dirichlet boundaries
|
|---|
| 1231 | c ub = u on Dirichlet boundaries
|
|---|
| 1232 | c
|
|---|
| 1233 | c _
|
|---|
| 1234 | c A ( u0 + ub ) = 0
|
|---|
| 1235 | c
|
|---|
| 1236 | c _ _
|
|---|
| 1237 | c A u0 = - A ub
|
|---|
| 1238 | c
|
|---|
| 1239 | c _ _
|
|---|
| 1240 | c MAM u0 = -M A ub, M is the mask
|
|---|
| 1241 | c
|
|---|
| 1242 | c _
|
|---|
| 1243 | c A u0 = -M A ub , Helmholtz solve with SPD matrix A
|
|---|
| 1244 | c
|
|---|
| 1245 | c u = u0+ub
|
|---|
| 1246 | c
|
|---|
| 1247 | include 'SIZE'
|
|---|
| 1248 | include 'TOTAL'
|
|---|
| 1249 | include 'CTIMER'
|
|---|
| 1250 | c
|
|---|
| 1251 | character*4 name
|
|---|
| 1252 | real u(1),mask(1),mult(1),approx (1)
|
|---|
| 1253 | integer napprox(1)
|
|---|
| 1254 |
|
|---|
| 1255 | parameter (lt=lx1*ly1*lz1*lelt)
|
|---|
| 1256 | common /scrvh/ h1(lt),h2(lt)
|
|---|
| 1257 | common /scruz/ r (lt),ub(lt)
|
|---|
| 1258 |
|
|---|
| 1259 | logical ifstdh
|
|---|
| 1260 | character*4 cname
|
|---|
| 1261 | character*6 name6
|
|---|
| 1262 |
|
|---|
| 1263 | logical ifwt,ifvec
|
|---|
| 1264 |
|
|---|
| 1265 | call chcopy(cname,name,4)
|
|---|
| 1266 | call capit (cname,4)
|
|---|
| 1267 |
|
|---|
| 1268 | call blank (name6,6)
|
|---|
| 1269 | call chcopy(name6,name,4)
|
|---|
| 1270 | ifwt = .true.
|
|---|
| 1271 | ifvec = .false.
|
|---|
| 1272 | isd = 1
|
|---|
| 1273 | imsh = 1
|
|---|
| 1274 | nel = nelfld(ifld)
|
|---|
| 1275 |
|
|---|
| 1276 | n = lx1*ly1*lz1*nel
|
|---|
| 1277 |
|
|---|
| 1278 | call copy (ub,u,n) ! ub = u on boundary
|
|---|
| 1279 | call dsavg(ub) ! Make certain ub is in H1
|
|---|
| 1280 | call rone (h1,n)
|
|---|
| 1281 | call rzero(h2,n)
|
|---|
| 1282 | ! _
|
|---|
| 1283 | call axhelm (r,ub,h1,h2,1,1) ! r = A*ub
|
|---|
| 1284 |
|
|---|
| 1285 | do i=1,n ! _
|
|---|
| 1286 | r(i)=-r(i)*mask(i) ! r = -M*A*ub
|
|---|
| 1287 | enddo
|
|---|
| 1288 |
|
|---|
| 1289 | call dssum (r,lx1,ly1,lz1) ! dssum rhs
|
|---|
| 1290 |
|
|---|
| 1291 | call project1
|
|---|
| 1292 | $ (r,n,approx,napprox,h1,h2,mask,mult,ifwt,ifvec,name6)
|
|---|
| 1293 |
|
|---|
| 1294 | if (nel.eq.nelv) then
|
|---|
| 1295 | call hmhzpf (name,u,r,h1,h2,mask,mult,imsh,tol,maxi,isd,binvm1)
|
|---|
| 1296 | else
|
|---|
| 1297 | call hmhzpf (name,u,r,h1,h2,mask,mult,imsh,tol,maxi,isd,bintm1)
|
|---|
| 1298 | endif
|
|---|
| 1299 |
|
|---|
| 1300 | call project2
|
|---|
| 1301 | $ (u,n,approx,napprox,h1,h2,mask,mult,ifwt,ifvec,name6)
|
|---|
| 1302 |
|
|---|
| 1303 | call add2(u,ub,n)
|
|---|
| 1304 |
|
|---|
| 1305 | return
|
|---|
| 1306 | end
|
|---|
| 1307 | c-----------------------------------------------------------------------
|
|---|
| 1308 | subroutine givens_rotation(a, b, c, s, r)
|
|---|
| 1309 |
|
|---|
| 1310 | real a, b, c, s, r, h, d
|
|---|
| 1311 |
|
|---|
| 1312 | if(b.ne.0.0) then
|
|---|
| 1313 | h = hypot(a,b) !Can use library version or below implementation
|
|---|
| 1314 | d = 1.0/h
|
|---|
| 1315 | c = abs(a)*d
|
|---|
| 1316 | s = sign(d,a)*b
|
|---|
| 1317 | r = sign(1.0,a)*h
|
|---|
| 1318 | else
|
|---|
| 1319 | c = 1.0
|
|---|
| 1320 | s = 0.0
|
|---|
| 1321 | r = a
|
|---|
| 1322 | endif
|
|---|
| 1323 |
|
|---|
| 1324 | return
|
|---|
| 1325 | end
|
|---|
| 1326 | c-----------------------------------------------------------------------
|
|---|
| 1327 | ! Compilers with Fortran 2008 support should have a
|
|---|
| 1328 | ! library implementation of this (gfortran and ifort do).
|
|---|
| 1329 |
|
|---|
| 1330 | real function hypot(a, b) !Does not handle a = b = 0 case
|
|---|
| 1331 |
|
|---|
| 1332 | real a, b, t, x, c, d, ix
|
|---|
| 1333 |
|
|---|
| 1334 | c = abs(a)
|
|---|
| 1335 | d = abs(b)
|
|---|
| 1336 | x = max(c,d)
|
|---|
| 1337 | ix = 1.0/x
|
|---|
| 1338 | t = min(c,d)
|
|---|
| 1339 | t = ix*t
|
|---|
| 1340 |
|
|---|
| 1341 | hypot = x*sqrt(1.0+t*t)
|
|---|
| 1342 |
|
|---|
| 1343 | return
|
|---|
| 1344 | end
|
|---|
| 1345 | c-----------------------------------------------------------------------
|
|---|
