| 1 | c-----------------------------------------------------------------------
|
|---|
| 2 | subroutine fluidp (igeom)
|
|---|
| 3 | c
|
|---|
| 4 | c Driver for perturbation velocity
|
|---|
| 5 | c
|
|---|
| 6 | include 'SIZE'
|
|---|
| 7 | include 'INPUT'
|
|---|
| 8 | include 'TSTEP'
|
|---|
| 9 | include 'SOLN'
|
|---|
| 10 |
|
|---|
| 11 | do jp=1,npert
|
|---|
| 12 |
|
|---|
| 13 | if (nio.eq.0.and.igeom.eq.2) write(6,1) istep,time,jp
|
|---|
| 14 | 1 format(i9,1pe14.7,' Perturbation Solve:',i5)
|
|---|
| 15 |
|
|---|
| 16 | call perturbv (igeom)
|
|---|
| 17 |
|
|---|
| 18 | enddo
|
|---|
| 19 |
|
|---|
| 20 | jp=0 ! set jp to zero, for baseline flow
|
|---|
| 21 |
|
|---|
| 22 | return
|
|---|
| 23 | end
|
|---|
| 24 | c-----------------------------------------------------------------------
|
|---|
| 25 | subroutine perturbv (igeom)
|
|---|
| 26 | c
|
|---|
| 27 | c
|
|---|
| 28 | c Solve the convection-diffusion equation for the perturbation field,
|
|---|
| 29 | c with projection onto a div-free space.
|
|---|
| 30 | c
|
|---|
| 31 | c
|
|---|
| 32 | include 'SIZE'
|
|---|
| 33 | include 'INPUT'
|
|---|
| 34 | include 'EIGEN'
|
|---|
| 35 | include 'SOLN'
|
|---|
| 36 | include 'TSTEP'
|
|---|
| 37 | include 'MASS'
|
|---|
| 38 | C
|
|---|
| 39 | COMMON /SCRNS/ RESV1 (LX1,LY1,LZ1,LELV)
|
|---|
| 40 | $ , RESV2 (LX1,LY1,LZ1,LELV)
|
|---|
| 41 | $ , RESV3 (LX1,LY1,LZ1,LELV)
|
|---|
| 42 | $ , DV1 (LX1,LY1,LZ1,LELV)
|
|---|
| 43 | $ , DV2 (LX1,LY1,LZ1,LELV)
|
|---|
| 44 | $ , DV3 (LX1,LY1,LZ1,LELV)
|
|---|
| 45 | COMMON /SCRVH/ H1 (LX1,LY1,LZ1,LELV)
|
|---|
| 46 | $ , H2 (LX1,LY1,LZ1,LELV)
|
|---|
| 47 | c
|
|---|
| 48 | ifield = 1
|
|---|
| 49 | c
|
|---|
| 50 | if (igeom.eq.1) then
|
|---|
| 51 | c
|
|---|
| 52 | c Old geometry, old velocity
|
|---|
| 53 | c
|
|---|
| 54 | call makefp
|
|---|
| 55 | call lagfieldp
|
|---|
| 56 | c
|
|---|
| 57 | else
|
|---|
| 58 | c
|
|---|
| 59 | c New geometry, new velocity
|
|---|
| 60 | c
|
|---|
| 61 | intype = -1
|
|---|
| 62 | call sethlm (h1,h2,intype)
|
|---|
| 63 | call cresvipp (resv1,resv2,resv3,h1,h2)
|
|---|
| 64 | call ophinv (dv1,dv2,dv3,resv1,resv2,resv3,h1,h2,tolhv,nmxv)
|
|---|
| 65 | call opadd2 (vxp(1,jp),vyp(1,jp),vzp(1,jp),dv1,dv2,dv3)
|
|---|
| 66 | call incomprp (vxp(1,jp),vyp(1,jp),vzp(1,jp),prp(1,jp))
|
|---|
| 67 | c
|
|---|
| 68 | endif
|
|---|
| 69 | c
|
|---|
| 70 | return
|
|---|
| 71 | end
|
|---|
| 72 | c-----------------------------------------------------------------------
|
|---|
| 73 | subroutine lagfieldp
|
|---|
| 74 | c
|
|---|
| 75 | c Keep old Vp-field(s)
|
|---|
| 76 | c
|
|---|
| 77 | include 'SIZE'
|
|---|
| 78 | include 'INPUT'
|
|---|
| 79 | include 'SOLN'
|
|---|
| 80 | include 'TSTEP'
|
|---|
| 81 | c
|
|---|
| 82 | do ilag=nbdinp-1,2,-1
|
|---|
| 83 | call opcopy
|
|---|
| 84 | $ (vxlagp(1,ilag ,jp),vylagp(1,ilag ,jp),vzlagp(1,ilag ,jp)
|
|---|
| 85 | $ ,vxlagp(1,ilag-1,jp),vylagp(1,ilag-1,jp),vzlagp(1,ilag-1,jp))
|
|---|
| 86 | enddo
|
|---|
| 87 | call opcopy(vxlagp(1,1,jp),vylagp(1,1,jp),vzlagp(1,1,jp)
|
|---|
| 88 | $ ,vxp (1,jp) ,vyp (1,jp) ,vzp (1,jp) )
|
|---|
| 89 | c
|
|---|
| 90 | return
|
|---|
| 91 | end
|
|---|
| 92 | c-----------------------------------------------------------------------
|
|---|
| 93 | subroutine makefp
|
|---|
| 94 | c
|
|---|
| 95 | c Make rhs for velocity perturbation equation
|
|---|
| 96 | c
|
|---|
| 97 | include 'SIZE'
|
|---|
| 98 | include 'SOLN'
|
|---|
| 99 | include 'MASS'
|
|---|
| 100 | include 'INPUT'
|
|---|
| 101 | include 'TSTEP'
|
|---|
| 102 | include 'ADJOINT'
|
|---|
| 103 |
|
|---|
| 104 | call makeufp
|
|---|
| 105 | if (filterType.eq.2) call make_hpf
|
|---|
| 106 | if (ifnav.and.(.not.ifchar).and.(.not.ifadj))call advabp
|
|---|
| 107 | if (ifnav.and.(.not.ifchar).and.( ifadj))call advabp_adjoint
|
|---|
| 108 | if (iftran) call makextp
|
|---|
| 109 | call makebdfp
|
|---|
| 110 | c
|
|---|
| 111 | return
|
|---|
| 112 | end
|
|---|
| 113 | c--------------------------------------------------------------------
|
|---|
| 114 | subroutine makeufp
|
|---|
| 115 | c
|
|---|
| 116 | c Compute and add: (1) user specified forcing function (FX,FY,FZ)
|
|---|
| 117 | c
|
|---|
| 118 | include 'SIZE'
|
|---|
| 119 | include 'SOLN'
|
|---|
| 120 | include 'MASS'
|
|---|
| 121 | include 'TSTEP'
|
|---|
| 122 | C
|
|---|
| 123 | time = time-dt
|
|---|
| 124 | call nekuf (bfxp(1,jp),bfyp(1,jp),bfzp(1,jp))
|
|---|
| 125 | call opcolv2 (bfxp(1,jp),bfyp(1,jp),bfzp(1,jp)
|
|---|
| 126 | $ ,vtrans(1,1,1,1,ifield),bm1)
|
|---|
| 127 | time = time+dt
|
|---|
| 128 | c
|
|---|
| 129 | return
|
|---|
| 130 | end
|
|---|
| 131 | c--------------------------------------------------------------------
|
|---|
| 132 | subroutine advabp
|
|---|
| 133 | C
|
|---|
| 134 | C Eulerian scheme, add convection term to forcing function
|
|---|
| 135 | C at current time step.
|
|---|
| 136 | C
|
|---|
| 137 | include 'SIZE'
|
|---|
| 138 | include 'INPUT'
|
|---|
| 139 | include 'SOLN'
|
|---|
| 140 | include 'MASS'
|
|---|
| 141 | include 'TSTEP'
|
|---|
| 142 | C
|
|---|
| 143 | COMMON /SCRNS/ TA1 (LX1*LY1*LZ1*LELV)
|
|---|
| 144 | $ , TA2 (LX1*LY1*LZ1*LELV)
|
|---|
| 145 | $ , TA3 (LX1*LY1*LZ1*LELV)
|
|---|
| 146 | $ , TB1 (LX1*LY1*LZ1*LELV)
|
|---|
| 147 | $ , TB2 (LX1*LY1*LZ1*LELV)
|
|---|
| 148 | $ , TB3 (LX1*LY1*LZ1*LELV)
|
|---|
| 149 | C
|
|---|
| 150 | ntot1 = lx1*ly1*lz1*nelv
|
|---|
| 151 | ntot2 = lx2*ly2*lz2*nelv
|
|---|
| 152 | c
|
|---|
| 153 | if (if3d) then
|
|---|
| 154 | call opcopy (tb1,tb2,tb3,vx,vy,vz) ! Save velocity
|
|---|
| 155 | call opcopy (vx,vy,vz,vxp(1,jp),vyp(1,jp),vzp(1,jp)) ! U <-- dU
|
|---|
| 156 | call convop (ta1,tb1) ! du.grad U
|
|---|
| 157 | call convop (ta2,tb2)
|
|---|
| 158 | call convop (ta3,tb3)
|
|---|
| 159 | call opcopy (vx,vy,vz,tb1,tb2,tb3) ! Restore velocity
|
|---|
| 160 | c
|
|---|
| 161 | do i=1,ntot1
|
|---|
| 162 | tmp = bm1(i,1,1,1)*vtrans(i,1,1,1,ifield)
|
|---|
| 163 | bfxp(i,jp) = bfxp(i,jp)-tmp*ta1(i)
|
|---|
| 164 | bfyp(i,jp) = bfyp(i,jp)-tmp*ta2(i)
|
|---|
| 165 | bfzp(i,jp) = bfzp(i,jp)-tmp*ta3(i)
|
|---|
| 166 | enddo
|
|---|
| 167 | c
|
|---|
| 168 | call convop (ta1,vxp(1,jp)) ! U.grad dU
|
|---|
| 169 | call convop (ta2,vyp(1,jp))
|
|---|
| 170 | call convop (ta3,vzp(1,jp))
|
|---|
| 171 | c
|
|---|
| 172 | do i=1,ntot1
|
|---|
| 173 | tmp = bm1(i,1,1,1)*vtrans(i,1,1,1,ifield)
|
|---|
| 174 | bfxp(i,jp) = bfxp(i,jp)-tmp*ta1(i)
|
|---|
| 175 | bfyp(i,jp) = bfyp(i,jp)-tmp*ta2(i)
|
|---|
| 176 | bfzp(i,jp) = bfzp(i,jp)-tmp*ta3(i)
|
|---|
| 177 | enddo
|
|---|
| 178 | c
|
|---|
| 179 | else ! 2D
|
|---|
| 180 | c
|
|---|
| 181 | call opcopy (tb1,tb2,tb3,vx,vy,vz) ! Save velocity
|
|---|
| 182 | call opcopy (vx,vy,vz,vxp(1,jp),vyp(1,jp),vzp(1,jp)) ! U <-- dU
|
|---|
| 183 | call convop (ta1,tb1) ! du.grad U
|
|---|
| 184 | call convop (ta2,tb2)
|
|---|
| 185 | call opcopy (vx,vy,vz,tb1,tb2,tb3) ! Restore velocity
|
|---|
| 186 | c
|
|---|
| 187 | do i=1,ntot1
|
|---|
| 188 | tmp = bm1(i,1,1,1)*vtrans(i,1,1,1,ifield)
|
|---|
| 189 | bfxp(i,jp) = bfxp(i,jp)-tmp*ta1(i)
|
|---|
| 190 | bfyp(i,jp) = bfyp(i,jp)-tmp*ta2(i)
|
|---|
| 191 | enddo
|
|---|
| 192 | c
|
|---|
| 193 | call convop (ta1,vxp(1,jp)) ! U.grad dU
|
|---|
| 194 | call convop (ta2,vyp(1,jp))
|
|---|
| 195 | c
|
|---|
| 196 | do i=1,ntot1
|
|---|
| 197 | tmp = bm1(i,1,1,1)*vtrans(i,1,1,1,ifield)
|
|---|
| 198 | bfxp(i,jp) = bfxp(i,jp)-tmp*ta1(i)
|
|---|
| 199 | bfyp(i,jp) = bfyp(i,jp)-tmp*ta2(i)
|
|---|
| 200 | enddo
|
|---|
| 201 | c
|
|---|
| 202 | endif
|
|---|
| 203 | c
|
|---|
| 204 | return
|
|---|
| 205 | end
|
|---|
| 206 | c--------------------------------------------------------------------
|
|---|
| 207 | subroutine advabp_adjoint
|
|---|
| 208 | C
|
|---|
| 209 | C Eulerian scheme, add convection term to forcing function
|
|---|
| 210 | C at current time step for backward part of adjoint:
|
|---|
| 211 | C Convective term is now (U.Grad)u - (Grad U)^T .u
|
|---|
| 212 | C instead of (u.Grad)U + (U.Grad)u in above subroutine advabp
|
|---|
| 213 | C
|
|---|
| 214 | include 'SIZE'
|
|---|
| 215 | include 'INPUT'
|
|---|
| 216 | include 'SOLN'
|
|---|
| 217 | include 'MASS'
|
|---|
| 218 | include 'TSTEP'
|
|---|
| 219 | include 'GEOM'
|
|---|
| 220 | include 'ADJOINT'
|
|---|
| 221 | C
|
|---|
| 222 | COMMON /SCRNS/ TA1 (LX1*LY1*LZ1*LELV)
|
|---|
| 223 | $ , TA2 (LX1*LY1*LZ1*LELV)
|
|---|
| 224 | $ , TA3 (LX1*LY1*LZ1*LELV)
|
|---|
| 225 | $ , TB1 (LX1*LY1*LZ1*LELV)
|
|---|
| 226 | $ , TB2 (LX1*LY1*LZ1*LELV)
|
|---|
| 227 | $ , TB3 (LX1*LY1*LZ1*LELV)
|
|---|
| 228 |
|
|---|
| 229 |
|
|---|
| 230 | real fact,factx,facty
|
|---|
| 231 | C
|
|---|
| 232 | ntot1 = lx1*ly1*lz1*nelv
|
|---|
| 233 | ntot2 = lx2*ly2*lz2*nelv !dimensionn arrays
|
|---|
| 234 | NTOT = lx1*ly1*lz1*NELT
|
|---|
| 235 |
|
|---|
| 236 | c
|
|---|
| 237 | if (if3d) then
|
|---|
| 238 | call opcopy (tb1,tb2,tb3,vx,vy,vz) ! Save velocity
|
|---|
| 239 | call opcopy (vx,vy,vz,vxp(1,jp),vyp(1,jp),vzp(1,jp)) ! U <-- u
|
|---|
| 240 | c
|
|---|
| 241 | call convop_adj (ta1,ta2,ta3,tb1,tb2,tb3,vx,vy,vz) ! u.grad U^T
|
|---|
| 242 |
|
|---|
| 243 | call opcopy (vx,vy,vz,tb1,tb2,tb3) ! Restore velocity
|
|---|
| 244 | c
|
|---|
| 245 | do i=1,ntot1
|
|---|
| 246 | tmp = bm1(i,1,1,1)*vtrans(i,1,1,1,ifield)
|
|---|
| 247 | bfxp(i,jp) = bfxp(i,jp)-tmp*ta1(i)
|
|---|
| 248 | bfyp(i,jp) = bfyp(i,jp)-tmp*ta2(i)
|
|---|
| 249 | bfzp(i,jp) = bfzp(i,jp)-tmp*ta3(i)
|
|---|
| 250 | enddo
|
|---|
| 251 | c
|
|---|
| 252 | c
|
|---|
| 253 | c
|
|---|
| 254 | call convop (ta1,vxp(1,jp)) ! U.grad u
|
|---|
| 255 | call convop (ta2,vyp(1,jp))
|
|---|
| 256 | call convop (ta3,vzp(1,jp))
|
|---|
| 257 | c
|
|---|
| 258 | do i=1,ntot1
|
|---|
| 259 | tmp = bm1(i,1,1,1)*vtrans(i,1,1,1,ifield)
|
|---|
| 260 | bfxp(i,jp) = bfxp(i,jp)+tmp*ta1(i)
|
|---|
| 261 | bfyp(i,jp) = bfyp(i,jp)+tmp*ta2(i)
|
|---|
| 262 | bfzp(i,jp) = bfzp(i,jp)+tmp*ta3(i)
|
|---|
| 263 | enddo
|
|---|
| 264 | c
|
|---|
| 265 | if (ifheat) then ! dt.(grad T)
|
|---|
| 266 | c ! term coming from the temperature convection
|
|---|
| 267 | call opcolv3 (ta1,ta2,ta3,dTdx,dTdy,dTdz,tp)
|
|---|
| 268 | c
|
|---|
| 269 | do i=1,ntot1
|
|---|
| 270 | tmp = bm1(i,1,1,1)*vtrans(i,1,1,1,2)
|
|---|
| 271 | bfxp(i,jp) = bfxp(i,jp)-tmp*ta1(i)
|
|---|
| 272 | bfyp(i,jp) = bfyp(i,jp)-tmp*ta2(i)
|
|---|
| 273 | bfzp(i,jp) = bfzp(i,jp)-tmp*ta3(i)
|
|---|
| 274 | enddo
|
|---|
| 275 | endif
|
|---|
| 276 | c
|
|---|
| 277 | else ! 2D
|
|---|
| 278 |
|
|---|
| 279 | call opcopy (tb1,tb2,tb3,vx,vy,vz) ! Save velocity
|
|---|
| 280 | call opcopy (vx,vy,vz,vxp(1,jp),vyp(1,jp),vzp(1,jp))
|
|---|
| 281 |
|
|---|
| 282 | call convop_adj (ta1,ta2,ta3,tb1,tb2,tb3,vx,vy,vz) ! u.((grad U)^T)
|
|---|
| 283 |
|
|---|
| 284 | call opcopy (vx,vy,vz,tb1,tb2,tb3) ! Restore velocity
|
|---|
| 285 |
|
|---|
| 286 | do i=1,ntot1
|
|---|
| 287 | tmp = bm1(i,1,1,1)*vtrans(i,1,1,1,ifield)
|
|---|
| 288 | bfxp(i,jp) = bfxp(i,jp)-tmp*ta1(i)
|
|---|
| 289 | bfyp(i,jp) = bfyp(i,jp)-tmp*ta2(i)
|
|---|
| 290 | enddo
|
|---|
| 291 |
|
|---|
| 292 | call convop (ta1,vxp(1,jp)) ! U.grad u
|
|---|
| 293 | call convop (ta2,vyp(1,jp))
|
|---|
| 294 | c
|
|---|
| 295 | do i=1,ntot1
|
|---|
| 296 | tmp = bm1(i,1,1,1)*vtrans(i,1,1,1,ifield)
|
|---|
| 297 | bfxp(i,jp) = bfxp(i,jp)+tmp*ta1(i)
|
|---|
| 298 | bfyp(i,jp) = bfyp(i,jp)+tmp*ta2(i)
|
|---|
| 299 | enddo
|
|---|
| 300 |
|
|---|
| 301 | if (ifheat) then ! dt.(grad T)^T
|
|---|
| 302 | ! term coming from the temperature convection
|
|---|
| 303 | call opcolv3 (ta1,ta2,ta3,dTdx,dTdy,dTdz,tp)
|
|---|
| 304 | c
|
|---|
| 305 | do i=1,ntot1
|
|---|
| 306 | tmp = bm1(i,1,1,1)*vtrans(i,1,1,1,2)
|
|---|
| 307 | bfxp(i,jp) = bfxp(i,jp)-tmp*ta1(i)
|
|---|
| 308 | bfyp(i,jp) = bfyp(i,jp)-tmp*ta2(i)
|
|---|
| 309 | enddo
|
|---|
| 310 | endif
|
|---|
| 311 |
|
|---|
| 312 | endif
|
|---|
| 313 | c
|
|---|
| 314 | return
|
|---|
| 315 | end
|
|---|
| 316 | c--------------------------------------------------------------------
|
|---|
| 317 | subroutine convop_adj (bdux,bduy,bduz,udx,udy,udz,cx,cy,cz)
|
|---|
| 318 |
|
|---|
| 319 | include 'SIZE'
|
|---|
| 320 | include 'TOTAL'
|
|---|
| 321 |
|
|---|
| 322 | parameter (lxy=lx1*ly1*lz1,ltd=lxd*lyd*lzd)
|
|---|
| 323 | common /scrcv/ fx(ltd),fy(ltd),fz(ltd)
|
|---|
| 324 | $ , uf1(ltd),uf2(ltd),uf3(ltd),uf4(ltd),uf5(ltd),uf6(ltd)
|
|---|
| 325 | real urx(ltd),usx(ltd),utx(ltd)
|
|---|
| 326 | real ury(ltd),usy(ltd),uty(ltd)
|
|---|
| 327 | real urz(ltd),usz(ltd),utz(ltd)
|
|---|
| 328 | real bdux(1),bduy(1),bduz(1),u(1),cx(1),cy(1),cz(1)
|
|---|
| 329 | real udx(1),udy(1),udz(1)
|
|---|
| 330 | logical ifuf,ifcf ! u and/or c already on fine mesh?
|
|---|
| 331 | integer e
|
|---|
| 332 | real bdrx(1), bdry(1),bdrz (1)
|
|---|
| 333 |
|
|---|
| 334 | call set_dealias_rx
|
|---|
| 335 | nxyz1 = lx1*ly1*lz1
|
|---|
| 336 | c AM DING DING
|
|---|
| 337 | nxyzd = lxd*lyd*lzd
|
|---|
| 338 | nxyzu = nxyz1
|
|---|
| 339 | nxyzc = nxyz1
|
|---|
| 340 | ntot1=lx1*ly1*lz1*nelv
|
|---|
| 341 | ic = 1 ! pointer to vector field C
|
|---|
| 342 | iu = 1 ! pointer to scalar field u
|
|---|
| 343 | ib = 1 ! pointer to scalar field Bdu
|
|---|
| 344 | if(if3d)then
|
|---|
| 345 | do e=1,nelv
|
|---|
| 346 | ! map coarse velocity to fine mesh (C-->F)
|
|---|
| 347 | call intp_rstd(fx,cx(ic),lx1,lxd,if3d,0) ! 0 --> forward
|
|---|
| 348 | call intp_rstd(fy,cy(ic),lx1,lxd,if3d,0)
|
|---|
| 349 | call intp_rstd(fz,cz(ic),lx1,lxd,if3d,0)
|
|---|
| 350 |
|
|---|
| 351 | call intp_rstd(uf1,udx(iu),lx1,lxd,if3d,0) ! 0 --> forward
|
|---|
| 352 | call grad_rst(urx,usx,utx,uf1,lxd,if3d)
|
|---|
| 353 |
|
|---|
| 354 | call intp_rstd(uf2,udy(iu),lx1,lxd,if3d,0)
|
|---|
| 355 | call grad_rst(ury,usy,uty,uf2,lxd,if3d)
|
|---|
| 356 |
|
|---|
| 357 | call intp_rstd(uf3,udz(iu),lx1,lxd,if3d,0)
|
|---|
| 358 | call grad_rst(urz,usz,utz,uf3,lxd,if3d)
|
|---|
| 359 |
|
|---|
| 360 | do i=1,nxyzd ! mass matrix included, per DFM (4.8.5)
|
|---|
| 361 | uf4(i)=fx(i)*(rx(i,1,e)*urx(i)+rx(i,4,e)*usx(i)
|
|---|
| 362 | $ +rx(i,7,e)*utx(i))+
|
|---|
| 363 | $ fy(i)*(rx(i,1,e)*ury(i)+rx(i,4,e)*usy(i)
|
|---|
| 364 | $ +rx(i,7,e)*uty(i))+
|
|---|
| 365 | $ fz(i)*(rx(i,1,e)*urz(i)+rx(i,4,e)*usz(i)
|
|---|
| 366 | $ +rx(i,7,e)*utz(i))
|
|---|
| 367 | uf5(i)=fx(i)*(rx(i,2,e)*urx(i)+rx(i,5,e)*usx(i)
|
|---|
| 368 | $ +rx(i,8,e)*utx(i))+
|
|---|
| 369 | $ fy(i)*(rx(i,2,e)*ury(i)+rx(i,5,e)*usy(i)
|
|---|
| 370 | $ +rx(i,8,e)*uty(i))+
|
|---|
| 371 | $ fz(i)*(rx(i,2,e)*urz(i)+rx(i,5,e)*usz(i)
|
|---|
| 372 | $ +rx(i,8,e)*utz(i))
|
|---|
| 373 | uf6(i)=fx(i)*(rx(i,3,e)*urx(i)+rx(i,6,e)*usx(i)
|
|---|
| 374 | $ +rx(i,9,e)*utx(i))+
|
|---|
| 375 | $ fy(i)*(rx(i,3,e)*ury(i)+rx(i,6,e)*usy(i)
|
|---|
| 376 | $ +rx(i,9,e)*uty(i))+
|
|---|
| 377 | $ fz(i)*(rx(i,3,e)*urz(i)+rx(i,6,e)*usz(i)
|
|---|
| 378 | $ +rx(i,9,e)*utz(i))
|
|---|
| 379 | enddo
|
|---|
| 380 |
|
|---|
| 381 | call intp_rstd(bdux(ib),uf4,lx1,lxd,if3d,1) ! Project back to coarse
|
|---|
| 382 | call intp_rstd(bduy(ib),uf5,lx1,lxd,if3d,1)
|
|---|
| 383 | call intp_rstd(bduz(ib),uf6,lx1,lxd,if3d,1)
|
|---|
| 384 |
|
|---|
| 385 | ic = ic + nxyzc
|
|---|
| 386 | iu = iu + nxyzu
|
|---|
| 387 | ib = ib + nxyz1
|
|---|
| 388 | enddo
|
|---|
| 389 | call invcol2 (bdux,bm1,ntot1) ! local mass inverse
|
|---|
| 390 | call invcol2 (bduy,bm1,ntot1) ! local mass inverse
|
|---|
| 391 | call invcol2 (bduz,bm1,ntot1) ! local mass inverse
|
|---|
| 392 | else
|
|---|
| 393 | do e=1,nelv
|
|---|
| 394 |
|
|---|
| 395 | ! map coarse velocity to fine mesh (C-->F)
|
|---|
| 396 | call intp_rstd(fx,cx(ic),lx1,lxd,if3d,0) ! 0 --> forward
|
|---|
| 397 | call intp_rstd(fy,cy(ic),lx1,lxd,if3d,0)
|
|---|
| 398 |
|
|---|
| 399 | call intp_rstd(uf1,udx(iu),lx1,lxd,if3d,0)
|
|---|
| 400 | call grad_rst(urx,usx,utx,uf1,lxd,if3d)
|
|---|
| 401 |
|
|---|
| 402 | call intp_rstd(uf2,udy(iu),lx1,lxd,if3d,0)
|
|---|
| 403 | call grad_rst(ury,usy,uty,uf2,lxd,if3d)
|
|---|
| 404 |
|
|---|
| 405 | do i=1,nxyzd
|
|---|
| 406 | uf4(i) = fx(i)*(rx(i,1,e)*urx(i)+rx(i,3,e)*usx(i))+
|
|---|
| 407 | $ fy(i)*(rx(i,1,e)*ury(i)+rx(i,3,e)*usy(i))
|
|---|
| 408 | uf5(i) = fx(i)*(rx(i,2,e)*urx(i)+rx(i,4,e)*usx(i))+
|
|---|
| 409 | $ fy(i)*(rx(i,2,e)*ury(i)+rx(i,4,e)*usy(i))
|
|---|
| 410 | enddo
|
|---|
| 411 |
|
|---|
| 412 | call intp_rstd(bdux(ib),uf4,lx1,lxd,if3d,1)
|
|---|
| 413 | call intp_rstd(bduy(ib),uf5,lx1,lxd,if3d,1)
|
|---|
| 414 |
|
|---|
| 415 | ic = ic + nxyzc
|
|---|
| 416 | iu = iu + nxyzu
|
|---|
| 417 | ib = ib + nxyz1
|
|---|
| 418 | end do
|
|---|
| 419 | call invcol2 (bdux,bm1,ntot1) ! local mass inverse
|
|---|
| 420 | call invcol2 (bduy,bm1,ntot1) ! local mass inverse
|
|---|
| 421 | end if
|
|---|
| 422 |
|
|---|
| 423 |
|
|---|
| 424 | return
|
|---|
| 425 | end
|
|---|
| 426 | c--------------------------------------------------------------------
|
|---|
| 427 | subroutine makextp
|
|---|
| 428 | c
|
|---|
| 429 | c Add extrapolation terms to perturbation source terms
|
|---|
| 430 | c
|
|---|
| 431 | c (nek5 equivalent for velocity is "makeabf")
|
|---|
| 432 | c
|
|---|
| 433 | include 'SIZE'
|
|---|
| 434 | include 'INPUT'
|
|---|
| 435 | include 'SOLN'
|
|---|
| 436 | include 'MASS'
|
|---|
| 437 | include 'TSTEP'
|
|---|
| 438 | C
|
|---|
| 439 | common /scrns/ ta1 (lx1,ly1,lz1,lelv)
|
|---|
| 440 | $ , ta2 (lx1,ly1,lz1,lelv)
|
|---|
| 441 | $ , ta3 (lx1,ly1,lz1,lelv)
|
|---|
| 442 | c
|
|---|
| 443 | ntot1 = lx1*ly1*lz1*nelv
|
|---|
| 444 | c
|
|---|
| 445 | ab0 = ab(1)
|
|---|
| 446 | ab1 = ab(2)
|
|---|
| 447 | ab2 = ab(3)
|
|---|
| 448 | call add3s2 (ta1,exx1p(1,jp),exx2p(1,jp),ab1,ab2,ntot1)
|
|---|
| 449 | call add3s2 (ta2,exy1p(1,jp),exy2p(1,jp),ab1,ab2,ntot1)
|
|---|
| 450 | call copy (exx2p(1,jp),exx1p(1,jp),ntot1)
|
|---|
| 451 | call copy (exy2p(1,jp),exy1p(1,jp),ntot1)
|
|---|
| 452 | call copy (exx1p(1,jp),bfxp (1,jp),ntot1)
|
|---|
| 453 | call copy (exy1p(1,jp),bfyp (1,jp),ntot1)
|
|---|
| 454 | call add2s1 (bfxp(1,jp),ta1,ab0,ntot1)
|
|---|
| 455 | call add2s1 (bfyp(1,jp),ta2,ab0,ntot1)
|
|---|
| 456 | if (if3d) then
|
|---|
| 457 | call add3s2 (ta3,exz1p(1,jp),exz2p(1,jp),ab1,ab2,ntot1)
|
|---|
| 458 | call copy (exz2p(1,jp),exz1p(1,jp),ntot1)
|
|---|
| 459 | call copy (exz1p(1,jp),bfzp (1,jp),ntot1)
|
|---|
| 460 | call add2s1 (bfzp(1,jp),ta3,ab0,ntot1)
|
|---|
| 461 | endif
|
|---|
| 462 | c
|
|---|
| 463 | return
|
|---|
| 464 | end
|
|---|
| 465 | c--------------------------------------------------------------------
|
|---|
| 466 | subroutine makebdfp
|
|---|
| 467 | C
|
|---|
| 468 | C Add contributions to perturbation source from lagged BD terms.
|
|---|
| 469 | C
|
|---|
| 470 | include 'SIZE'
|
|---|
| 471 | include 'SOLN'
|
|---|
| 472 | include 'MASS'
|
|---|
| 473 | include 'GEOM'
|
|---|
| 474 | include 'INPUT'
|
|---|
| 475 | include 'TSTEP'
|
|---|
| 476 | C
|
|---|
| 477 | COMMON /SCRNS/ TA1(LX1,LY1,LZ1,LELV)
|
|---|
| 478 | $ , TA2(LX1,LY1,LZ1,LELV)
|
|---|
| 479 | $ , TA3(LX1,LY1,LZ1,LELV)
|
|---|
| 480 | $ , TB1(LX1,LY1,LZ1,LELV)
|
|---|
| 481 | $ , TB2(LX1,LY1,LZ1,LELV)
|
|---|
| 482 | $ , TB3(LX1,LY1,LZ1,LELV)
|
|---|
| 483 | $ , H2 (LX1,LY1,LZ1,LELV)
|
|---|
| 484 | C
|
|---|
| 485 | ntot1 = lx1*ly1*lz1*nelv
|
|---|
| 486 | const = 1./dt
|
|---|
| 487 | call cmult2(h2,vtrans(1,1,1,1,ifield),const,ntot1)
|
|---|
| 488 | call opcolv3c (tb1,tb2,tb3
|
|---|
| 489 | $ ,vxp(1,jp),vyp(1,jp),vzp(1,jp),bm1,bd(2))
|
|---|
| 490 | C
|
|---|
| 491 | do ilag=2,nbd
|
|---|
| 492 | if (ifgeom) then
|
|---|
| 493 | call opcolv3c(ta1,ta2,ta3,vxlagp(1,ilag-1,jp),
|
|---|
| 494 | $ vylagp(1,ilag-1,jp),
|
|---|
| 495 | $ vzlagp(1,ilag-1,jp),
|
|---|
| 496 | $ bm1lag(1,1,1,1,ilag-1),bd(ilag+1))
|
|---|
| 497 | else
|
|---|
| 498 | call opcolv3c(ta1,ta2,ta3,vxlagp(1,ilag-1,jp),
|
|---|
| 499 | $ vylagp(1,ilag-1,jp),
|
|---|
| 500 | $ vzlagp(1,ilag-1,jp),
|
|---|
| 501 | $ bm1 ,bd(ilag+1))
|
|---|
| 502 | endif
|
|---|
| 503 | call opadd2 (tb1,tb2,tb3,ta1,ta2,ta3)
|
|---|
| 504 | enddo
|
|---|
| 505 | call opadd2col (bfxp(1,jp),bfyp(1,jp),bfzp(1,jp),tb1,tb2,tb3,h2)
|
|---|
| 506 | c
|
|---|
| 507 | return
|
|---|
| 508 | end
|
|---|
| 509 | c-----------------------------------------------------------------------
|
|---|
| 510 | subroutine cresvipp(resv1,resv2,resv3,h1,h2)
|
|---|
| 511 | c
|
|---|
| 512 | c Account for inhomogeneous Dirichlet boundary contributions
|
|---|
| 513 | c in rhs of perturbation eqn.
|
|---|
| 514 | c n
|
|---|
| 515 | c Also, subtract off best estimate of grad p
|
|---|
| 516 | c
|
|---|
| 517 | include 'SIZE'
|
|---|
| 518 | include 'TOTAL'
|
|---|
| 519 | real resv1 (lx1,ly1,lz1,1)
|
|---|
| 520 | real resv2 (lx1,ly1,lz1,1)
|
|---|
| 521 | real resv3 (lx1,ly1,lz1,1)
|
|---|
| 522 | real h1 (lx1,ly1,lz1,1)
|
|---|
| 523 | real h2 (lx1,ly1,lz1,1)
|
|---|
| 524 | common /scruz/ w1 (lx1,ly1,lz1,lelv)
|
|---|
| 525 | $ , w2 (lx1,ly1,lz1,lelv)
|
|---|
| 526 | $ , w3 (lx1,ly1,lz1,lelv)
|
|---|
| 527 | $ , prextr(lx2,ly2,lz2,lelv)
|
|---|
| 528 | c
|
|---|
| 529 | ntot1 = lx1*ly1*lz1*nelv
|
|---|
| 530 | ntot2 = lx2*ly2*lz2*nelv
|
|---|
| 531 | c
|
|---|
| 532 | call bcdirvc (vxp(1,jp),vyp(1,jp),vzp(1,jp)
|
|---|
| 533 | $ ,v1mask,v2mask,v3mask)
|
|---|
| 534 | call extrapprp(prextr)
|
|---|
| 535 | call opgradt(resv1,resv2,resv3,prextr)
|
|---|
| 536 | call opadd2(resv1,resv2,resv3,bfxp(1,jp),bfyp(1,jp),bfzp(1,jp))
|
|---|
| 537 | call ophx (w1,w2,w3,vxp(1,jp),vyp(1,jp),vzp(1,jp),h1,h2)
|
|---|
| 538 | call opsub2(resv1,resv2,resv3,w1,w2,w3)
|
|---|
| 539 | c
|
|---|
| 540 | return
|
|---|
| 541 | end
|
|---|
| 542 | c--------------------------------------------------------------------
|
|---|
| 543 | subroutine heatp (igeom)
|
|---|
| 544 | C
|
|---|
| 545 | C Driver for temperature or passive scalar.
|
|---|
| 546 | C
|
|---|
| 547 | C Current version:
|
|---|
| 548 | C (1) Varaiable properties.
|
|---|
| 549 | C (2) Implicit time stepping.
|
|---|
| 550 | C (3) User specified tolerance for the Helmholtz solver
|
|---|
| 551 | C (not based on eigenvalues).
|
|---|
| 552 | C (4) A passive scalar can be defined on either the
|
|---|
| 553 | C temperatur or the velocity mesh.
|
|---|
| 554 | C (5) A passive scalar has its own multiplicity (B.C.).
|
|---|
| 555 | C
|
|---|
| 556 | INCLUDE 'SIZE'
|
|---|
| 557 | INCLUDE 'INPUT'
|
|---|
| 558 | INCLUDE 'TSTEP'
|
|---|
| 559 | INCLUDE 'SOLN'
|
|---|
| 560 | C
|
|---|
| 561 | do jp=1,npert
|
|---|
| 562 | do ifield=2,nfield
|
|---|
| 563 | INTYPE = -1
|
|---|
| 564 | IF (.NOT.IFTMSH(IFIELD)) IMESH = 1
|
|---|
| 565 | IF ( IFTMSH(IFIELD)) IMESH = 2
|
|---|
| 566 | CALL UNORM
|
|---|
| 567 | CALL SETTOLT
|
|---|
| 568 | CALL CDSCALP (IGEOM)
|
|---|
| 569 | enddo
|
|---|
| 570 | enddo
|
|---|
| 571 | c
|
|---|
| 572 | jp=0 ! set jp to zero, for baseline flow
|
|---|
| 573 | c
|
|---|
| 574 | return
|
|---|
| 575 | end
|
|---|
| 576 | C
|
|---|
| 577 | C-----------------------------------------------------------------------
|
|---|
| 578 | subroutine cdscalp (igeom)
|
|---|
| 579 | C-----------------------------------------------------------------------
|
|---|
| 580 | C
|
|---|
| 581 | C Solve the convection-diffusion equation for passive scalar IPSCAL
|
|---|
| 582 | C
|
|---|
| 583 | C-----------------------------------------------------------------------
|
|---|
| 584 | INCLUDE 'SIZE'
|
|---|
| 585 | INCLUDE 'INPUT'
|
|---|
| 586 | INCLUDE 'GEOM'
|
|---|
| 587 | INCLUDE 'MVGEOM'
|
|---|
| 588 | INCLUDE 'SOLN'
|
|---|
| 589 | INCLUDE 'MASS'
|
|---|
| 590 | INCLUDE 'TSTEP'
|
|---|
| 591 | COMMON /CPRINT/ IFPRINT
|
|---|
| 592 | LOGICAL IFPRINT
|
|---|
| 593 | LOGICAL IFCONV
|
|---|
| 594 | C
|
|---|
| 595 | COMMON /SCRNS/ TA(LX1,LY1,LZ1,LELT)
|
|---|
| 596 | $ ,TB(LX1,LY1,LZ1,LELT)
|
|---|
| 597 | COMMON /SCRVH/ H1(LX1,LY1,LZ1,LELT)
|
|---|
| 598 | $ ,H2(LX1,LY1,LZ1,LELT)
|
|---|
| 599 | c
|
|---|
| 600 | include 'ORTHOT'
|
|---|
| 601 | ifld1 = ifield-1
|
|---|
| 602 | napproxt(1,ifld1) = laxtt
|
|---|
| 603 | C
|
|---|
| 604 | IF (IGEOM.EQ.1) THEN
|
|---|
| 605 | C
|
|---|
| 606 | C Old geometry
|
|---|
| 607 | C
|
|---|
| 608 | CALL MAKEQP
|
|---|
| 609 | CALL LAGSCALP
|
|---|
| 610 | C
|
|---|
| 611 | ELSE
|
|---|
| 612 | C
|
|---|
| 613 | IF (IFPRINT) THEN
|
|---|
| 614 | IF (IFIELD.EQ.2.AND.NIO.EQ.0) THEN
|
|---|
| 615 | WRITE (6,*) ' Temperature/Passive scalar solution'
|
|---|
| 616 | ENDIF
|
|---|
| 617 | ENDIF
|
|---|
| 618 | if1=ifield-1
|
|---|
| 619 | write(name4t,1) if1
|
|---|
| 620 | 1 format('TEM',i1)
|
|---|
| 621 | C
|
|---|
| 622 | C New geometry
|
|---|
| 623 | C
|
|---|
| 624 | NEL = NELFLD(IFIELD)
|
|---|
| 625 | NTOT = lx1*ly1*lz1*NEL
|
|---|
| 626 | C
|
|---|
| 627 | INTYPE = 0
|
|---|
| 628 | IF (IFTRAN) INTYPE = -1
|
|---|
| 629 | CALL SETHLM (H1,H2,INTYPE)
|
|---|
| 630 | CALL BCNEUSC (TA,-1)
|
|---|
| 631 | CALL ADD2 (H2,TA,NTOT)
|
|---|
| 632 | CALL BCDIRSC (TP(1,IFIELD-1,jp))
|
|---|
| 633 | CALL AXHELM (TA,TP (1,IFIELD-1,jp),H1,H2,IMESH,1)
|
|---|
| 634 | CALL SUB3 (TB,BQP(1,IFIELD-1,jp),TA,NTOT)
|
|---|
| 635 | CALL BCNEUSC (TA,1)
|
|---|
| 636 | CALL ADD2 (TB,TA,NTOT)
|
|---|
| 637 | c
|
|---|
| 638 | CALL HMHOLTZ (name4t,TA,TB,H1,H2
|
|---|
| 639 | $ ,TMASK(1,1,1,1,IFIELD-1)
|
|---|
| 640 | $ ,TMULT(1,1,1,1,IFIELD-1)
|
|---|
| 641 | $ ,IMESH,TOLHT(IFIELD),NMXT(ifield-1),1)
|
|---|
| 642 | c
|
|---|
| 643 | CALL ADD2 (TP(1,IFIELD-1,jp),TA,NTOT)
|
|---|
| 644 | C
|
|---|
| 645 | CALL BCNEUSC (TA,1)
|
|---|
| 646 | CALL ADD2 (BQP(1,IFIELD-1,jp),TA,NTOT)
|
|---|
| 647 | C
|
|---|
| 648 | ENDIF
|
|---|
| 649 | C
|
|---|
| 650 | return
|
|---|
| 651 | end
|
|---|
| 652 | C
|
|---|
| 653 | subroutine makeqp
|
|---|
| 654 | C----------------------------------------------------------------------
|
|---|
| 655 | C
|
|---|
| 656 | C Generate forcing function for the solution of a passive scalar.
|
|---|
| 657 | C !! NOTE: Do not change the content of the array BQ until the current
|
|---|
| 658 | C time step is completed
|
|---|
| 659 | C
|
|---|
| 660 | C----------------------------------------------------------------------
|
|---|
| 661 | INCLUDE 'SIZE'
|
|---|
| 662 | INCLUDE 'GEOM'
|
|---|
| 663 | INCLUDE 'INPUT'
|
|---|
| 664 | INCLUDE 'TSTEP'
|
|---|
| 665 | include 'SOLN'
|
|---|
| 666 | CALL MAKEUQP
|
|---|
| 667 | IF (FILTERTYPE.EQ.2) CALL MAKE_HPF
|
|---|
| 668 | IF (IFADVC(IFIELD).AND.(.NOT.IFCHAR)) CALL CONVABP
|
|---|
| 669 | IF (IFTRAN) CALL MAKEABQP
|
|---|
| 670 | IF ((IFTRAN.AND..NOT.IFCHAR).OR.
|
|---|
| 671 | $ (IFTRAN.AND..NOT.IFADVC(IFIELD).AND.IFCHAR)) CALL MAKEBDQP
|
|---|
| 672 | c IF (IFADVC(IFIELD).AND.IFCHAR) CALL CONVCHP
|
|---|
| 673 | IF (IFADVC(IFIELD).AND.IFCHAR) write(6,*) 'no convchp'
|
|---|
| 674 | IF (IFADVC(IFIELD).AND.IFCHAR) call exitt
|
|---|
| 675 | c
|
|---|
| 676 | return
|
|---|
| 677 | end
|
|---|
| 678 | C
|
|---|
| 679 | subroutine makeuqp
|
|---|
| 680 | C---------------------------------------------------------------------
|
|---|
| 681 | C
|
|---|
| 682 | C Fill up user defined forcing function and collocate will the
|
|---|
| 683 | C mass matrix on the Gauss-Lobatto mesh.
|
|---|
| 684 | C
|
|---|
| 685 | C---------------------------------------------------------------------
|
|---|
| 686 | INCLUDE 'SIZE'
|
|---|
| 687 | INCLUDE 'MASS'
|
|---|
| 688 | INCLUDE 'SOLN'
|
|---|
| 689 | INCLUDE 'TSTEP'
|
|---|
| 690 | C
|
|---|
| 691 | NTOT = lx1*ly1*lz1*NELFLD(IFIELD)
|
|---|
| 692 | C
|
|---|
| 693 | time = time-dt ! time is tn
|
|---|
| 694 | c
|
|---|
| 695 | call rzero ( bqp(1,ifield-1,jp) , ntot)
|
|---|
| 696 | call setqvol ( bqp(1,ifield-1,jp) )
|
|---|
| 697 | call col2 ( bqp(1,ifield-1,jp) ,bm1,ntot)
|
|---|
| 698 | c
|
|---|
| 699 | time = time+dt ! restore time
|
|---|
| 700 | C
|
|---|
| 701 | return
|
|---|
| 702 | end
|
|---|
| 703 | C
|
|---|
| 704 | subroutine convabp
|
|---|
| 705 | C---------------------------------------------------------------
|
|---|
| 706 | C
|
|---|
| 707 | C Eulerian scheme, add convection term to forcing function
|
|---|
| 708 | C at current time step.
|
|---|
| 709 | C
|
|---|
| 710 | C---------------------------------------------------------------
|
|---|
| 711 | include 'SIZE'
|
|---|
| 712 | include 'ADJOINT'
|
|---|
| 713 | include 'SOLN'
|
|---|
| 714 | include 'MASS'
|
|---|
| 715 | include 'TSTEP'
|
|---|
| 716 | c
|
|---|
| 717 | common /scruz/ ta (lx1,ly1,lz1,lelt)
|
|---|
| 718 | $ , ua (lx1,ly1,lz1,lelt)
|
|---|
| 719 | $ , ub (lx1,ly1,lz1,lelt)
|
|---|
| 720 | $ , uc (lx1,ly1,lz1,lelt)
|
|---|
| 721 | real coeff
|
|---|
| 722 | c
|
|---|
| 723 | nel = nelfld(ifield)
|
|---|
| 724 | ntot1 = lx1*ly1*lz1*nel
|
|---|
| 725 | c
|
|---|
| 726 | if (.not.ifadj) then
|
|---|
| 727 | call opcopy(ua,ub,uc,vx,vy,vz)
|
|---|
| 728 | call opcopy(vx,vy,vz,vxp(1,jp),vyp(1,jp),vzp(1,jp))
|
|---|
| 729 | call convop(ta,t(1,1,1,1,ifield-1)) ! dU.grad T
|
|---|
| 730 | call opcopy(vx,vy,vz,ua,ub,uc)
|
|---|
| 731 | endif
|
|---|
| 732 | c
|
|---|
| 733 | call convop (ua,tp(1,ifield-1,jp)) ! U.grad dT
|
|---|
| 734 | c
|
|---|
| 735 | if (.not.ifadj) then
|
|---|
| 736 | call add2 (ta,ua,ntot1) ! U.grad dT + dU.grad T
|
|---|
| 737 | else
|
|---|
| 738 | call copy (ta,ua,ntot1)
|
|---|
| 739 | coeff=-1.0
|
|---|
| 740 | call cmult(ta,coeff,ntot1) ! -U.grad dT
|
|---|
| 741 | ! the second term depends on the buoyancy
|
|---|
| 742 | endif
|
|---|
| 743 | c
|
|---|
| 744 | call col2 (ta,vtrans(1,1,1,1,ifield),ntot1) !vtrans (U.grad dT + dU.grad T)
|
|---|
| 745 | call subcol3 (bqp(1,ifield-1,jp),bm1,ta,ntot1)
|
|---|
| 746 | c
|
|---|
| 747 | return
|
|---|
| 748 | end
|
|---|
| 749 | C
|
|---|
| 750 | subroutine makeabqp
|
|---|
| 751 | C-----------------------------------------------------------------------
|
|---|
| 752 | C
|
|---|
| 753 | C Sum up contributions to 3rd order Adams-Bashforth scheme.
|
|---|
| 754 | C
|
|---|
| 755 | C-----------------------------------------------------------------------
|
|---|
| 756 | INCLUDE 'SIZE'
|
|---|
| 757 | INCLUDE 'SOLN'
|
|---|
| 758 | INCLUDE 'TSTEP'
|
|---|
| 759 | C
|
|---|
| 760 | COMMON /SCRUZ/ TA (LX1,LY1,LZ1,LELT)
|
|---|
| 761 | C
|
|---|
| 762 | AB0 = AB(1)
|
|---|
| 763 | AB1 = AB(2)
|
|---|
| 764 | AB2 = AB(3)
|
|---|
| 765 | NEL = NELFLD(IFIELD)
|
|---|
| 766 | NTOT1 = lx1*ly1*lz1*NEL
|
|---|
| 767 | C
|
|---|
| 768 | CALL ADD3S2 (TA,VGRADT1P(1,IFIELD-1,jp),
|
|---|
| 769 | $ VGRADT2P(1,IFIELD-1,jp),AB1,AB2,NTOT1)
|
|---|
| 770 | CALL COPY ( VGRADT2P(1,IFIELD-1,jp),
|
|---|
| 771 | $ VGRADT1P(1,IFIELD-1,jp),NTOT1)
|
|---|
| 772 | CALL COPY ( VGRADT1P(1,IFIELD-1,jp),
|
|---|
| 773 | $ BQP(1,IFIELD-1,jp),NTOT1)
|
|---|
| 774 | CALL ADD2S1 (BQP(1,IFIELD-1,jp),TA,AB0,NTOT1)
|
|---|
| 775 | C
|
|---|
| 776 | return
|
|---|
| 777 | end
|
|---|
| 778 | C
|
|---|
| 779 | subroutine makebdqp
|
|---|
| 780 | C-----------------------------------------------------------------------
|
|---|
| 781 | C
|
|---|
| 782 | C Add contributions to Q from lagged BD terms.
|
|---|
| 783 | C
|
|---|
| 784 | C-----------------------------------------------------------------------
|
|---|
| 785 | INCLUDE 'SIZE'
|
|---|
| 786 | INCLUDE 'SOLN'
|
|---|
| 787 | INCLUDE 'MASS'
|
|---|
| 788 | INCLUDE 'GEOM'
|
|---|
| 789 | INCLUDE 'INPUT'
|
|---|
| 790 | INCLUDE 'TSTEP'
|
|---|
| 791 | C
|
|---|
| 792 | COMMON /SCRNS/ TA (LX1,LY1,LZ1,LELT)
|
|---|
| 793 | $ , TB (LX1,LY1,LZ1,LELT)
|
|---|
| 794 | $ , H2 (LX1,LY1,LZ1,LELT)
|
|---|
| 795 | C
|
|---|
| 796 | NEL = NELFLD(IFIELD)
|
|---|
| 797 | NTOT1 = lx1*ly1*lz1*NEL
|
|---|
| 798 | CONST = 1./DT
|
|---|
| 799 | CALL COPY (H2,VTRANS(1,1,1,1,IFIELD),NTOT1)
|
|---|
| 800 | CALL CMULT (H2,CONST,NTOT1)
|
|---|
| 801 | C
|
|---|
| 802 | CALL COL3 (TB,BM1,TP(1,IFIELD-1,jp),NTOT1)
|
|---|
| 803 | CALL CMULT (TB,BD(2),NTOT1)
|
|---|
| 804 | C
|
|---|
| 805 | DO 100 ILAG=2,NBD
|
|---|
| 806 | IF (IFGEOM) THEN
|
|---|
| 807 | CALL COL3 (TA,BM1LAG(1,1,1,1,ILAG-1),
|
|---|
| 808 | $ TLAGP (1,ILAG-1,IFIELD-1,jp),NTOT1)
|
|---|
| 809 | ELSE
|
|---|
| 810 | CALL COL3 (TA,BM1,
|
|---|
| 811 | $ TLAGP (1,ILAG-1,IFIELD-1,jp),NTOT1)
|
|---|
| 812 | ENDIF
|
|---|
| 813 | CALL ADD2S2(TB,TA,BD(ilag+1),NTOT1)
|
|---|
| 814 | 100 CONTINUE
|
|---|
| 815 | C
|
|---|
| 816 | CALL COL2 (TB,H2,NTOT1)
|
|---|
| 817 | CALL ADD2 (BQP(1,IFIELD-1,jp),TB,NTOT1)
|
|---|
| 818 | C
|
|---|
| 819 | return
|
|---|
| 820 | end
|
|---|
| 821 | C
|
|---|
| 822 | subroutine lagscalp
|
|---|
| 823 | C-----------------------------------------------------------------------
|
|---|
| 824 | C
|
|---|
| 825 | C Keep old passive scalar field(s)
|
|---|
| 826 | C
|
|---|
| 827 | C-----------------------------------------------------------------------
|
|---|
| 828 | INCLUDE 'SIZE'
|
|---|
| 829 | INCLUDE 'INPUT'
|
|---|
| 830 | INCLUDE 'SOLN'
|
|---|
| 831 | INCLUDE 'TSTEP'
|
|---|
| 832 | C
|
|---|
| 833 | NTOT1 = lx1*ly1*lz1*NELFLD(IFIELD)
|
|---|
| 834 | C
|
|---|
| 835 | DO 100 ILAG=NBDINP-1,2,-1
|
|---|
| 836 | CALL COPY (TLAGP(1,ILAG ,IFIELD-1,jp),
|
|---|
| 837 | $ TLAGP(1,ILAG-1,IFIELD-1,jp),NTOT1)
|
|---|
| 838 | 100 CONTINUE
|
|---|
| 839 | C
|
|---|
| 840 | CALL COPY (TLAGP(1,1,IFIELD-1,jp),TP(1,IFIELD-1,jp),NTOT1)
|
|---|
| 841 | C
|
|---|
| 842 | return
|
|---|
| 843 | end
|
|---|
| 844 | c-----------------------------------------------------------------------
|
|---|
| 845 | subroutine incomprp (ux,uy,uz,up)
|
|---|
| 846 | c
|
|---|
| 847 | c Project U onto the closest incompressible field
|
|---|
| 848 | c
|
|---|
| 849 | c Input: U := (ux,uy,uz)
|
|---|
| 850 | c
|
|---|
| 851 | c Output: updated values of U, iproj, proj; and
|
|---|
| 852 | c up := pressure correction req'd to impose div U = 0
|
|---|
| 853 | c
|
|---|
| 854 | c
|
|---|
| 855 | c Dependencies: ifield ==> which "density" (vtrans) is used.
|
|---|
| 856 | c
|
|---|
| 857 | c
|
|---|
| 858 | include 'SIZE'
|
|---|
| 859 | include 'TOTAL'
|
|---|
| 860 | include 'CTIMER'
|
|---|
| 861 | c
|
|---|
| 862 | common /scrns/ w1 (lx1,ly1,lz1,lelv)
|
|---|
| 863 | $ , w2 (lx1,ly1,lz1,lelv)
|
|---|
| 864 | $ , w3 (lx1,ly1,lz1,lelv)
|
|---|
| 865 | $ , dv1 (lx1,ly1,lz1,lelv)
|
|---|
| 866 | $ , dv2 (lx1,ly1,lz1,lelv)
|
|---|
| 867 | $ , dv3 (lx1,ly1,lz1,lelv)
|
|---|
| 868 | $ , dp (lx2,ly2,lz2,lelv)
|
|---|
| 869 | common /scrvh/ h1 (lx1,ly1,lz1,lelv)
|
|---|
| 870 | $ , h2 (lx1,ly1,lz1,lelv)
|
|---|
| 871 | common /scrhi/ h2inv (lx1,ly1,lz1,lelv)
|
|---|
| 872 | COMMON /SCRCH/ PREXTR(LX2,LY2,LZ2,LELV)
|
|---|
| 873 | logical ifprjp
|
|---|
| 874 |
|
|---|
| 875 | c
|
|---|
| 876 | if (icalld.eq.0) tpres=0.0
|
|---|
| 877 | icalld=icalld+1
|
|---|
| 878 | npres=icalld
|
|---|
| 879 | c
|
|---|
| 880 | ntot1 = lx1*ly1*lz1*nelv
|
|---|
| 881 | ntot2 = lx2*ly2*lz2*nelv
|
|---|
| 882 | intype = 1
|
|---|
| 883 | dtbd = bd(1)/dt
|
|---|
| 884 |
|
|---|
| 885 | call rzero (h1,ntot1)
|
|---|
| 886 | c call copy (h2,vtrans(1,1,1,1,ifield),ntot1)
|
|---|
| 887 | call cmult2 (h2,vtrans(1,1,1,1,ifield),dtbd,ntot1)
|
|---|
| 888 | call invers2 (h2inv,h2,ntot1)
|
|---|
| 889 |
|
|---|
| 890 | call opdiv (dp,ux,uy,uz)
|
|---|
| 891 | call chsign (dp,ntot2)
|
|---|
| 892 | call ortho (dp)
|
|---|
| 893 |
|
|---|
| 894 |
|
|---|
| 895 | C******************************************************************
|
|---|
| 896 |
|
|---|
| 897 |
|
|---|
| 898 | ifprjp=.false. ! project out previous pressure solutions?
|
|---|
| 899 | istart=param(95)
|
|---|
| 900 | if (istep.ge.istart.and.istart.ne.0) ifprjp=.true.
|
|---|
| 901 |
|
|---|
| 902 | ! Most likely, the following can be commented out. (pff, 1/6/2010)
|
|---|
| 903 | c if (npert.gt.1.or.ifbase) ifprjp=.false.
|
|---|
| 904 | cpff if (ifprjp) call setrhs (dp,h1,h2,h2inv)
|
|---|
| 905 |
|
|---|
| 906 | call esolver (dp,h1,h2,h2inv,intype)
|
|---|
| 907 |
|
|---|
| 908 | cpff if (ifprjp) call gensoln (dp,h1,h2,h2inv)
|
|---|
| 909 |
|
|---|
| 910 | cNOTE: The "cpff" comments added 11/24/17 to avoid old-style projection,
|
|---|
| 911 | cNOTE: which should be replaced with something more updated.
|
|---|
| 912 |
|
|---|
| 913 | C******************************************************************
|
|---|
| 914 |
|
|---|
| 915 | call opgradt (w1 ,w2 ,w3 ,dp)
|
|---|
| 916 | call opbinv (dv1,dv2,dv3,w1 ,w2 ,w3 ,h2inv)
|
|---|
| 917 | call opadd2 (ux ,uy ,uz ,dv1,dv2,dv3)
|
|---|
| 918 |
|
|---|
| 919 | call extrapprp(prextr)
|
|---|
| 920 | call lagpresp
|
|---|
| 921 | call add3(up,prextr,dp,ntot2)
|
|---|
| 922 |
|
|---|
| 923 | return
|
|---|
| 924 | end
|
|---|
| 925 | c------------------------------------------------------------------------
|
|---|
| 926 | subroutine extrapprp (prextr)
|
|---|
| 927 | C
|
|---|
| 928 | C Pressure extrapolation
|
|---|
| 929 | C
|
|---|
| 930 | INCLUDE 'SIZE'
|
|---|
| 931 | INCLUDE 'SOLN'
|
|---|
| 932 | INCLUDE 'TSTEP'
|
|---|
| 933 | COMMON /CTMP0/ DPR (LX2,LY2,LZ2,LELV)
|
|---|
| 934 | REAL PREXTR (LX2,LY2,LZ2,LELV)
|
|---|
| 935 | C
|
|---|
| 936 | ntot2 = lx2*ly2*lz2*nelv
|
|---|
| 937 | if (nbd.eq.1.or.nbd.eq.2) then
|
|---|
| 938 | call copy (prextr,prp(1,JP),ntot2)
|
|---|
| 939 | elseif (nbd.eq.3) then
|
|---|
| 940 | const = dtlag(1)/dtlag(2)
|
|---|
| 941 | call sub3 (dpr,prp(1,JP),prlagp(1,1,JP),ntot2)
|
|---|
| 942 | call cmult(dpr,const,ntot2)
|
|---|
| 943 | call add3 (prextr,prp(1,JP),dpr,ntot2)
|
|---|
| 944 | elseif (nbd.gt.3) then
|
|---|
| 945 | write (6,*) 'Pressure extrapolation cannot be completed'
|
|---|
| 946 | write (6,*) 'Try a lower-order temporal scheme'
|
|---|
| 947 | call exitt
|
|---|
| 948 | endif
|
|---|
| 949 | return
|
|---|
| 950 | end
|
|---|
| 951 | C-------------------------------------------------------------------
|
|---|
| 952 | subroutine lagpresp
|
|---|
| 953 | C
|
|---|
| 954 | C save old pressure values
|
|---|
| 955 | C
|
|---|
| 956 | INCLUDE 'SIZE'
|
|---|
| 957 | INCLUDE 'SOLN'
|
|---|
| 958 | INCLUDE 'TSTEP'
|
|---|
| 959 | C
|
|---|
| 960 | if (nbdinp.eq.3) then
|
|---|
| 961 | ntot2 = lx2*ly2*lz2*nelv
|
|---|
| 962 | call copy (prlagp(1,1,JP),prp(1,JP),ntot2)
|
|---|
| 963 | endif
|
|---|
| 964 | return
|
|---|
| 965 | end
|
|---|
| 966 | C-------------------------------------------------------------------
|
|---|
| 967 | subroutine lyap_scale ! Rescale / orthogonalize perturbation fields
|
|---|
| 968 | c
|
|---|
| 969 | c
|
|---|
| 970 | include 'SIZE'
|
|---|
| 971 | include 'TOTAL'
|
|---|
| 972 | c
|
|---|
| 973 | real sigma(0:lpert)
|
|---|
| 974 | c
|
|---|
| 975 | ntotv = lx1*ly1*lz1*nelv
|
|---|
| 976 | ntotp = lx2*ly2*lz2*nelv
|
|---|
| 977 | ntott = lx1*ly1*lz1*nelt
|
|---|
| 978 | c
|
|---|
| 979 | do j=1,npert
|
|---|
| 980 | call normvc(h1,semi,pl2,plinf,vxp(1,j),vyp(1,j),vzp(1,j))
|
|---|
| 981 | sigma(j) = pl2
|
|---|
| 982 | if (pl2.gt.0) then
|
|---|
| 983 | write(6,*) 'this is pl2:',pl2
|
|---|
| 984 | scale = 1./pl2
|
|---|
| 985 | call opcmult(vxp(1,j),vyp(1,j),vzp(1,j),scale)
|
|---|
| 986 | call cmult(tp(1,1,j),scale,ntott)
|
|---|
| 987 | call cmult(prp(1,j) ,scale,ntotp)
|
|---|
| 988 | endif
|
|---|
| 989 | c
|
|---|
| 990 | c Have to do lag terms as well, etc
|
|---|
| 991 | c
|
|---|
| 992 | c Also, must orthogonalize
|
|---|
| 993 | enddo
|
|---|
| 994 | c
|
|---|
| 995 | if (nio.eq.0) then
|
|---|
| 996 | if (npert.eq.1) write(6,1) istep,time,(sigma(j),j=1,npert)
|
|---|
| 997 | if (npert.eq.2) write(6,2) istep,time,(sigma(j),j=1,npert)
|
|---|
| 998 | if (npert.eq.3) write(6,3) istep,time,(sigma(j),j=1,npert)
|
|---|
| 999 | if (npert.eq.4) write(6,4) istep,time,(sigma(j),j=1,npert)
|
|---|
| 1000 | if (npert.eq.5) write(6,5) istep,time,(sigma(j),j=1,npert)
|
|---|
| 1001 | if (npert.eq.6) write(6,6) istep,time,(sigma(j),j=1,npert)
|
|---|
| 1002 | if (npert.eq.7) write(6,7) istep,time,(sigma(j),j=1,npert)
|
|---|
| 1003 | if (npert.eq.8) write(6,8) istep,time,(sigma(j),j=1,npert)
|
|---|
| 1004 | if (npert.eq.9) write(6,9) istep,time,(sigma(j),j=1,npert)
|
|---|
| 1005 | endif
|
|---|
| 1006 | 1 format(i8,1p2e12.4,' pgrow')
|
|---|
| 1007 | 2 format(i8,1p3e12.4,' pgrow')
|
|---|
| 1008 | 3 format(i8,1p4e12.4,' pgrow')
|
|---|
| 1009 | 4 format(i8,1p5e12.4,' pgrow')
|
|---|
| 1010 | 5 format(i8,1p6e12.4,' pgrow')
|
|---|
| 1011 | 6 format(i8,1p7e12.4,' pgrow')
|
|---|
| 1012 | 7 format(i8,1p8e12.4,' pgrow')
|
|---|
| 1013 | 8 format(i8,1p9e12.4,' pgrow')
|
|---|
| 1014 | 9 format(i8,1p10e12.4,' pgrow')
|
|---|
| 1015 | c
|
|---|
| 1016 | return
|
|---|
| 1017 | end
|
|---|
| 1018 | c-----------------------------------------------------------------------
|
|---|
| 1019 | subroutine out_pert ! dump perturbation .fld files
|
|---|
| 1020 | c
|
|---|
| 1021 | include 'SIZE'
|
|---|
| 1022 | include 'TOTAL'
|
|---|
| 1023 | include 'NEKUSE'
|
|---|
| 1024 | c
|
|---|
| 1025 | character*3 s3
|
|---|
| 1026 | c
|
|---|
| 1027 | do jpp=1,npert
|
|---|
| 1028 | write(s3,3) jpp
|
|---|
| 1029 | 3 format('p',i2.2)
|
|---|
| 1030 | call outpost2
|
|---|
| 1031 | $ (vxp(1,jpp),vyp(1,jpp),vzp(1,jpp),prp(1,jpp),tp(1,1,jpp),1,s3)
|
|---|
| 1032 | c call writehist
|
|---|
| 1033 | c $ (vxp(1,jpp),vyp(1,jpp),vzp(1,jpp),tp(1,1,jpp),jpp)
|
|---|
| 1034 | enddo
|
|---|
| 1035 | c
|
|---|
| 1036 | return
|
|---|
| 1037 | end
|
|---|
| 1038 | c-----------------------------------------------------------------------
|
|---|
| 1039 | subroutine pert_add2s2(i,j,scale) ! xi = xi + scale * xj
|
|---|
| 1040 | include 'SIZE'
|
|---|
| 1041 | include 'TOTAL'
|
|---|
| 1042 |
|
|---|
| 1043 | ntotp = lx2*ly2*lz2*nelv
|
|---|
| 1044 | ntotv = lx1*ly1*lz1*nelv
|
|---|
| 1045 | ntott = lx1*ly1*lz1*nelt
|
|---|
| 1046 |
|
|---|
| 1047 | call add2s2(vxp(1,i),vxp(1,j),scale,ntotv)
|
|---|
| 1048 | call add2s2(vyp(1,i),vyp(1,j),scale,ntotv)
|
|---|
| 1049 | if (if3d) call add2s2(vzp(1,i),vzp(1,j),scale,ntotv)
|
|---|
| 1050 | call add2s2(prp(1,i),prp(1,j),scale,ntotp)
|
|---|
| 1051 |
|
|---|
| 1052 | do l=1,lorder-1
|
|---|
| 1053 | call add2s2(vxlagp(1,l,i),vxlagp(1,l,j),scale,ntotv)
|
|---|
| 1054 | call add2s2(vylagp(1,l,i),vylagp(1,l,j),scale,ntotv)
|
|---|
| 1055 | if (if3d) call add2s2(vzlagp(1,l,i),vzlagp(1,l,j),scale,ntotv)
|
|---|
| 1056 | if (l.le.lorder-2)
|
|---|
| 1057 | $ call add2s2(prlagp(1,l,i),prlagp(1,l,j),scale,ntotp)
|
|---|
| 1058 | enddo
|
|---|
| 1059 |
|
|---|
| 1060 | call add2s2(exx1p(1,i),exx1p(1,j),scale,ntotv)
|
|---|
| 1061 | call add2s2(exy1p(1,i),exy1p(1,j),scale,ntotv)
|
|---|
| 1062 | if (if3d) call add2s2(exz1p(1,i),exz1p(1,j),scale,ntotv)
|
|---|
| 1063 |
|
|---|
| 1064 | call add2s2(exx2p(1,i),exx2p(1,j),scale,ntotv)
|
|---|
| 1065 | call add2s2(exy2p(1,i),exy2p(1,j),scale,ntotv)
|
|---|
| 1066 | if (if3d) call add2s2(exz2p(1,i),exz2p(1,j),scale,ntotv)
|
|---|
| 1067 |
|
|---|
| 1068 | if (ifheat) then
|
|---|
| 1069 | do k=0,npscal
|
|---|
| 1070 | k1=k+1
|
|---|
| 1071 | ntotk = lx1*ly1*lz1*nelfld(k+2)
|
|---|
| 1072 | call add2s2(tp(1,k1,i),tp(1,k1,j),scale,ntotk)
|
|---|
| 1073 | do l=1,lorder-1
|
|---|
| 1074 | call add2s2(tlagp(1,l,k1,i),tlagp(1,l,k1,j),scale,ntotk)
|
|---|
| 1075 | enddo
|
|---|
| 1076 | call add2s2(vgradt1p(1,k1,i),vgradt1p(1,k1,j),scale,ntotk)
|
|---|
| 1077 | call add2s2(vgradt2p(1,k1,i),vgradt2p(1,k1,j),scale,ntotk)
|
|---|
| 1078 | enddo
|
|---|
| 1079 | endif
|
|---|
| 1080 |
|
|---|
| 1081 | return
|
|---|
| 1082 | end
|
|---|
| 1083 | c-----------------------------------------------------------------------
|
|---|
| 1084 | function pert_inner_prod(i,j) ! weighted inner product vi^T vj
|
|---|
| 1085 | include 'SIZE'
|
|---|
| 1086 | include 'TOTAL'
|
|---|
| 1087 |
|
|---|
| 1088 | common/normset/pran, ray, rayc
|
|---|
| 1089 |
|
|---|
| 1090 | ntotv=lx1*ly1*lz1*nelv
|
|---|
| 1091 | ntott=lx1*ly1*lz1*nelt
|
|---|
| 1092 |
|
|---|
| 1093 | s1 = rayc*glsc3(vxp(1,i),bm1,vxp(1,j),ntotv)
|
|---|
| 1094 | s2 = rayc*glsc3(vyp(1,i),bm1,vyp(1,j),ntotv)
|
|---|
| 1095 | s3 = 0
|
|---|
| 1096 | if (if3d) s3 = rayc*glsc3(vzp(1,i),bm1,vzp(1,j),ntotv)
|
|---|
| 1097 |
|
|---|
| 1098 | t1 = 0
|
|---|
| 1099 | if (ifheat) t1=pran*ray*ray*glsc3(tp(1,1,i),bm1,tp(1,1,j),ntott)
|
|---|
| 1100 |
|
|---|
| 1101 | pert_inner_prod = (s1+s2+s3+t1)/volvm1
|
|---|
| 1102 |
|
|---|
| 1103 | return
|
|---|
| 1104 | end
|
|---|
| 1105 | c-----------------------------------------------------------------------
|
|---|
| 1106 | subroutine pert_ortho_norm ! orthogonalize and rescale pert. arrays
|
|---|
| 1107 | include 'SIZE'
|
|---|
| 1108 | include 'TOTAL'
|
|---|
| 1109 |
|
|---|
| 1110 | do k=1,npert
|
|---|
| 1111 | call pert_ortho_norm1(k)
|
|---|
| 1112 | enddo
|
|---|
| 1113 |
|
|---|
| 1114 | return
|
|---|
| 1115 | end
|
|---|
| 1116 | c-----------------------------------------------------------------------
|
|---|
| 1117 | subroutine pert_ortho_norm1 (k) ! orthogonalize k against 1...k-1
|
|---|
| 1118 | include 'SIZE'
|
|---|
| 1119 | include 'TOTAL'
|
|---|
| 1120 |
|
|---|
| 1121 | do j=1,k-1
|
|---|
| 1122 | scale = -pert_inner_prod(j,k)
|
|---|
| 1123 | call pert_add2s2(k,j,scale) ! xi = xi + scale * xj
|
|---|
| 1124 | enddo
|
|---|
| 1125 | scale = pert_inner_prod(k,k)
|
|---|
| 1126 | if (scale.gt.0) scale = 1./scale
|
|---|
| 1127 | call rescalepert(pertnorm,scale,k)
|
|---|
| 1128 |
|
|---|
| 1129 | return
|
|---|
| 1130 | end
|
|---|
| 1131 | c-----------------------------------------------------------------------
|
|---|
| 1132 | function opnorm2(v1,v2,v3)
|
|---|
| 1133 | include 'SIZE'
|
|---|
| 1134 | include 'TOTAL'
|
|---|
| 1135 | c
|
|---|
| 1136 | real v1(1) , v2(1), v3(1)
|
|---|
| 1137 | real normsq1,normsq2,normsq3,opnorm
|
|---|
| 1138 | c
|
|---|
| 1139 | ntotv=lx1*ly1*lz1*nelv
|
|---|
| 1140 | normsq1=glsc3(v1,bm1,v1,ntotv)
|
|---|
| 1141 | normsq2=glsc3(v2,bm1,v2,ntotv)
|
|---|
| 1142 | if(if3d) then
|
|---|
| 1143 | normsq3=glsc3(v3,bm1,v3,ntotv)
|
|---|
| 1144 | else
|
|---|
| 1145 | normsq3=0
|
|---|
| 1146 | endif
|
|---|
| 1147 |
|
|---|
| 1148 | opnorm2=normsq1+normsq2+normsq3
|
|---|
| 1149 | if (opnorm2.gt.0) opnorm2=sqrt(opnorm2/volvm1)
|
|---|
| 1150 | return
|
|---|
| 1151 | end
|
|---|
| 1152 | c-----------------------------------------------------------------------
|
|---|
| 1153 |
|
|---|
| 1154 | function Tnorm(temp)
|
|---|
| 1155 | include 'SIZE'
|
|---|
| 1156 | include 'TOTAL'
|
|---|
| 1157 |
|
|---|
| 1158 | real temp(*)
|
|---|
| 1159 | c
|
|---|
| 1160 | ntotv = lx1*ly1*lz1*nelv
|
|---|
| 1161 | Tnorm = sqrt( glsc3(temp,BM1,temp,ntotv) /voltm1)
|
|---|
| 1162 | c
|
|---|
| 1163 | return
|
|---|
| 1164 | end
|
|---|
| 1165 | c--------------------------------------------
|
|---|
| 1166 | function dmnorm1(v1,v2,v3,temp)
|
|---|
| 1167 | c Norm weighted by mass matrix
|
|---|
| 1168 | include 'SIZE'
|
|---|
| 1169 | include 'TOTAL'
|
|---|
| 1170 |
|
|---|
| 1171 | real v1(1),v2(1),v3(1),temp(1)
|
|---|
| 1172 | real normsq1,normsq2,normsq3,normsqT,dMnorm
|
|---|
| 1173 | common/normset/pran, ray, rayc
|
|---|
| 1174 |
|
|---|
| 1175 | ntotv=lx1*ly1*lz1*nelv
|
|---|
| 1176 | normsq1=(rayc)*glsc3(v1,BM1,v1,ntotv)
|
|---|
| 1177 | normsq2=(rayc)*glsc3(v2,BM1,v2,ntotv)
|
|---|
| 1178 | if(if3d) then
|
|---|
| 1179 | normsq3=(rayc)*glsc3(v3,BM1,v3,ntotv)
|
|---|
| 1180 | else
|
|---|
| 1181 | normsq3=0
|
|---|
| 1182 | endif
|
|---|
| 1183 |
|
|---|
| 1184 | if(ifheat) then
|
|---|
| 1185 | normsqT = (pran*ray*ray)*glsc3(temp,BM1,temp,ntotv)
|
|---|
| 1186 | else
|
|---|
| 1187 | normsqT = 0
|
|---|
| 1188 | endif
|
|---|
| 1189 |
|
|---|
| 1190 | dmnorm1=sqrt((normsq1+normsq2+normsq3+normsqT)/volvm1)
|
|---|
| 1191 |
|
|---|
| 1192 | return
|
|---|
| 1193 | end
|
|---|
| 1194 |
|
|---|
| 1195 | c---------------------------------------------------------------
|
|---|
| 1196 | subroutine opscale(v1,v2,v3,temp,alpha)
|
|---|
| 1197 | c v = alpha*v
|
|---|
| 1198 | include 'SIZE'
|
|---|
| 1199 | include 'INPUT'
|
|---|
| 1200 |
|
|---|
| 1201 | real alpha
|
|---|
| 1202 | real v1(1),v2(1),v3(1),temp(1)
|
|---|
| 1203 |
|
|---|
| 1204 | ltotv=lx1*ly1*lz1*lelv
|
|---|
| 1205 | ltott=lx1*ly1*lz1*lelt
|
|---|
| 1206 |
|
|---|
| 1207 | call cmult(v1,alpha,ltotv)
|
|---|
| 1208 | call cmult(v2,alpha,ltotv)
|
|---|
| 1209 | if (if3d) call cmult(v3,alpha,ltotv)
|
|---|
| 1210 | if (ifheat) call cmult(temp,alpha,ltott*ldimt)
|
|---|
| 1211 |
|
|---|
| 1212 | return
|
|---|
| 1213 | end
|
|---|
| 1214 |
|
|---|
| 1215 | c---------------------------------------------------------------
|
|---|
| 1216 | subroutine opscaleV(v1,v2,v3,alpha)
|
|---|
| 1217 | c v = alpha*v
|
|---|
| 1218 | include 'SIZE'
|
|---|
| 1219 | include 'INPUT'
|
|---|
| 1220 | real alpha
|
|---|
| 1221 | real v1(*),v2(*),v3(*)
|
|---|
| 1222 |
|
|---|
| 1223 | ntotv=lx1*ly1*lz1*nelv
|
|---|
| 1224 |
|
|---|
| 1225 | call cmult(v1,alpha,ntotv)
|
|---|
| 1226 | call cmult(v2,alpha,ntotv)
|
|---|
| 1227 |
|
|---|
| 1228 | if (if3d) call cmult(v3,alpha,ntotv)
|
|---|
| 1229 | c
|
|---|
| 1230 | return
|
|---|
| 1231 | end
|
|---|
| 1232 |
|
|---|
| 1233 | c-----------------------------------------------------------------------
|
|---|
| 1234 | subroutine computelyap
|
|---|
| 1235 | include 'SIZE'
|
|---|
| 1236 | include 'TOTAL'
|
|---|
| 1237 |
|
|---|
| 1238 | do jpp=1,npert ! Loop through each Lyapunov eigenvalue
|
|---|
| 1239 | call computelyap1
|
|---|
| 1240 | $ (vxp(1,jpp),vyp(1,jpp),vzp(1,jpp),tp(1,1,jpp),jpp)
|
|---|
| 1241 | enddo
|
|---|
| 1242 |
|
|---|
| 1243 | return
|
|---|
| 1244 | end
|
|---|
| 1245 | c-----------------------------------------------------------------------
|
|---|
| 1246 | subroutine computelyap1(vxq,vyq,vzq,tq,jpp)
|
|---|
| 1247 | include 'SIZE'
|
|---|
| 1248 | include 'TOTAL'
|
|---|
| 1249 |
|
|---|
| 1250 | real vxq(1),vyq(1),vzq(1),tq(1)
|
|---|
| 1251 | real lyapsum,twt
|
|---|
| 1252 | common /pertsave/ timestart,tinitial
|
|---|
| 1253 |
|
|---|
| 1254 | integer icount
|
|---|
| 1255 | save icount
|
|---|
| 1256 | data icount /0/
|
|---|
| 1257 |
|
|---|
| 1258 | logical if_restart,if_ortho_lyap
|
|---|
| 1259 | common/restar/if_restart,if_ortho_lyap
|
|---|
| 1260 |
|
|---|
| 1261 | character*132 lyprestart
|
|---|
| 1262 | common/restflename/lyprestart !file for restart data
|
|---|
| 1263 |
|
|---|
| 1264 | twt = param(126) !time to wait to start computing exponents
|
|---|
| 1265 |
|
|---|
| 1266 | if (nio.eq.0)
|
|---|
| 1267 | $ write(6,8) istep,icount,time,twt,(lyap(k,jpp),k=1,3),jpp
|
|---|
| 1268 | 8 format(i9,i4,2f8.4,1p3e12.4,i3,'clyap')
|
|---|
| 1269 |
|
|---|
| 1270 | if(time.lt.twt) then
|
|---|
| 1271 | c
|
|---|
| 1272 | c For a fresh simulation, then we wait 5 vertical diffusion
|
|---|
| 1273 | c times before we start measuring, so in this case just rescale
|
|---|
| 1274 | c
|
|---|
| 1275 | pertnorm = dmnorm1(vxq,vyq,vzq,tq)
|
|---|
| 1276 | pertinvnorm = 1.0/pertnorm
|
|---|
| 1277 | call rescalepert(pertnorm,pertinvnorm,jpp)
|
|---|
| 1278 | lyap(3,jpp) = pertnorm !store latest norm
|
|---|
| 1279 | timestart = time
|
|---|
| 1280 | tinitial = time
|
|---|
| 1281 | icount = 0
|
|---|
| 1282 | return
|
|---|
| 1283 | else
|
|---|
| 1284 | if (jpp.eq.1) icount = icount + 1
|
|---|
| 1285 | endif
|
|---|
| 1286 |
|
|---|
| 1287 | irescale = param(128)
|
|---|
| 1288 | if (icount.eq.irescale) then
|
|---|
| 1289 |
|
|---|
| 1290 | lyapsum = lyap(2,jpp)
|
|---|
| 1291 | oldpertnorm = lyap(3,jpp)
|
|---|
| 1292 | pertnorm=dmnorm1(vxq,vyq,vzq,tq)
|
|---|
| 1293 | c
|
|---|
| 1294 | lyap(1,jpp) = log(pertnorm/oldpertnorm)/(time-timestart)
|
|---|
| 1295 | lyapsum = lyapSum + log(pertnorm/oldpertnorm)
|
|---|
| 1296 | lyap(2,jpp) = lyapSum
|
|---|
| 1297 |
|
|---|
| 1298 | if(nid.eq.0) then ! write out results to the .lyp file
|
|---|
| 1299 |
|
|---|
| 1300 | write(6 ,1) istep,time,lyap(1,jpp),lyapsum,pertnorm,jpp
|
|---|
| 1301 | write(79,2) time,lyap(1,jpp),lyapsum,pertnorm,oldpertnorm,jpp
|
|---|
| 1302 | 1 format(i9,1p4e17.8,i4,'lyap')
|
|---|
| 1303 | 2 format(1p5e17.8,i4,'lyap')
|
|---|
| 1304 | c call flushbuffer(79)
|
|---|
| 1305 |
|
|---|
| 1306 | if (jpp.eq.1) open(unit=96,file=lyprestart)
|
|---|
| 1307 | write(96,9) lyapsum,timestart,jpp
|
|---|
| 1308 | 9 format(1p2e19.11,i9)
|
|---|
| 1309 | if (jpp.eq.npert) close(96)
|
|---|
| 1310 | endif
|
|---|
| 1311 |
|
|---|
| 1312 | pertinvnorm = 1.0/pertnorm
|
|---|
| 1313 | call rescalepert(pertnorm,pertinvnorm,jpp)
|
|---|
| 1314 | lyap(3,jpp) = pertnorm !save current pertnorm as old pertnorm
|
|---|
| 1315 |
|
|---|
| 1316 | if (jpp.eq.npert) then
|
|---|
| 1317 | icount = 0
|
|---|
| 1318 | timestart = time
|
|---|
| 1319 | endif
|
|---|
| 1320 |
|
|---|
| 1321 | if_ortho_lyap = .false.
|
|---|
| 1322 | if (param(125).gt.0) if_ortho_lyap = .true.
|
|---|
| 1323 | if (jpp.eq.npert .and. if_ortho_lyap) call pert_ortho_norm
|
|---|
| 1324 |
|
|---|
| 1325 | endif
|
|---|
| 1326 |
|
|---|
| 1327 | return
|
|---|
| 1328 | end
|
|---|
| 1329 | c-----------------------------------------------------------------------
|
|---|
| 1330 |
|
|---|
| 1331 | subroutine rescalepert(pertnorm,pertinvnorm,jpp)
|
|---|
| 1332 | include 'SIZE'
|
|---|
| 1333 | include 'TOTAL'
|
|---|
| 1334 |
|
|---|
| 1335 | ntotp = lx2*ly2*lz2*nelv
|
|---|
| 1336 |
|
|---|
| 1337 | call opscale !normalize vectors to unit norm
|
|---|
| 1338 | $ (vxp(1,jpp),vyp(1,jpp),vzp(1,jpp),tp(1,1,jpp),pertinvnorm)
|
|---|
| 1339 | call cmult(prp(1,jpp),pertinvnorm,ntotp)
|
|---|
| 1340 |
|
|---|
| 1341 | call opscale(exx1p(1,jpp),exy1p(1,jpp),exz1p(1,jpp)
|
|---|
| 1342 | $ ,vgradt1p(1,1,jpp),pertinvnorm)
|
|---|
| 1343 | call opscale(exx2p(1,jpp),exy2p(1,jpp),exz2p(1,jpp)
|
|---|
| 1344 | $ ,vgradt2p(1,1,jpp),pertinvnorm)
|
|---|
| 1345 |
|
|---|
| 1346 | ltotv = lx1*ly1*lz1*lelv
|
|---|
| 1347 | ltotp = lx2*ly2*lz2*lelv
|
|---|
| 1348 |
|
|---|
| 1349 | call cmult( tlagp(1,1,1,jpp),pertinvnorm,ltotv*(lorder-1)*ldimt)
|
|---|
| 1350 | call cmult(vxlagp(1,1,jpp),pertinvnorm,ltotv*(lorder-1))
|
|---|
| 1351 | call cmult(vylagp(1,1,jpp),pertinvnorm,ltotv*(lorder-1))
|
|---|
| 1352 | call cmult(vzlagp(1,1,jpp),pertinvnorm,ltotv*(lorder-1))
|
|---|
| 1353 | call cmult(prlagp(1,1,jpp),pertinvnorm,ltotp*(Lorder-2))
|
|---|
| 1354 |
|
|---|
| 1355 | if (nio.eq.0) write(6,1) istep,pertnorm,pertinvnorm,jpp,'PNORM'
|
|---|
| 1356 | 1 format(i4,1p2e12.4,i4,a5)
|
|---|
| 1357 | pertnorm = pertnorm*pertinvnorm
|
|---|
| 1358 |
|
|---|
| 1359 | return
|
|---|
| 1360 | end
|
|---|
| 1361 | c-----------------------------------------------------------------------
|
|---|
