| 1 | c-----------------------------------------------------------------------
|
|---|
| 2 | c
|
|---|
| 3 | c To do:
|
|---|
| 4 | c
|
|---|
| 5 | c Differing BC's imposed for ophinv, incomprn, etc.
|
|---|
| 6 | c
|
|---|
| 7 | c 1-shot Fast solver for Helmholtz and pressure
|
|---|
| 8 | c
|
|---|
| 9 | c
|
|---|
| 10 | c-----------------------------------------------------------------------
|
|---|
| 11 | subroutine induct (igeom)
|
|---|
| 12 | c
|
|---|
| 13 | c Solve the convection-diffusion equation for the B-field, with
|
|---|
| 14 | c projection onto a div-free space.
|
|---|
| 15 | c
|
|---|
| 16 | c
|
|---|
| 17 | include 'SIZE'
|
|---|
| 18 | include 'INPUT'
|
|---|
| 19 | include 'EIGEN'
|
|---|
| 20 | include 'SOLN'
|
|---|
| 21 | include 'TSTEP'
|
|---|
| 22 | include 'MASS'
|
|---|
| 23 |
|
|---|
| 24 | common /scrns/ resv1 (lx1,ly1,lz1,lelv)
|
|---|
| 25 | $ , resv2 (lx1,ly1,lz1,lelv)
|
|---|
| 26 | $ , resv3 (lx1,ly1,lz1,lelv)
|
|---|
| 27 | $ , dv1 (lx1,ly1,lz1,lelv)
|
|---|
| 28 | $ , dv2 (lx1,ly1,lz1,lelv)
|
|---|
| 29 | $ , dv3 (lx1,ly1,lz1,lelv)
|
|---|
| 30 | common /scrvh/ h1 (lx1,ly1,lz1,lelv)
|
|---|
| 31 | $ , h2 (lx1,ly1,lz1,lelv)
|
|---|
| 32 |
|
|---|
| 33 | ifield = ifldmhd
|
|---|
| 34 |
|
|---|
| 35 | if (igeom.eq.1) then ! old geometry, old velocity
|
|---|
| 36 |
|
|---|
| 37 | call makebsource_mhd
|
|---|
| 38 |
|
|---|
| 39 | else
|
|---|
| 40 |
|
|---|
| 41 | call lagbfield
|
|---|
| 42 | call lagvel
|
|---|
| 43 |
|
|---|
| 44 | call elsasserh(igeom)
|
|---|
| 45 |
|
|---|
| 46 | call vol_flow ! check for fixed flow rate
|
|---|
| 47 |
|
|---|
| 48 |
|
|---|
| 49 | endif
|
|---|
| 50 | c
|
|---|
| 51 | return
|
|---|
| 52 | end
|
|---|
| 53 | c--------------------------------------------------------------------
|
|---|
| 54 | subroutine lagbfield
|
|---|
| 55 | c
|
|---|
| 56 | c Keep old B-field(s)
|
|---|
| 57 | c
|
|---|
| 58 | include 'SIZE'
|
|---|
| 59 | include 'INPUT'
|
|---|
| 60 | include 'SOLN'
|
|---|
| 61 | include 'TSTEP'
|
|---|
| 62 | c
|
|---|
| 63 | do ilag=nbdinp-1,2,-1
|
|---|
| 64 | call opcopy
|
|---|
| 65 | $ ( bxlag(1,ilag ),bylag(1,ilag ),bzlag(1,ilag )
|
|---|
| 66 | $ , bxlag(1,ilag-1),bylag(1,ilag-1),bzlag(1,ilag-1) )
|
|---|
| 67 | enddo
|
|---|
| 68 | call opcopy (bxlag,bylag,bzlag,bx,by,bz)
|
|---|
| 69 | c
|
|---|
| 70 | return
|
|---|
| 71 | end
|
|---|
| 72 | c--------------------------------------------------------------------
|
|---|
| 73 | subroutine makebsource_mhd
|
|---|
| 74 | c
|
|---|
| 75 | c Make rhs for induction equation
|
|---|
| 76 | c
|
|---|
| 77 | include 'SIZE'
|
|---|
| 78 | include 'SOLN'
|
|---|
| 79 | include 'MASS'
|
|---|
| 80 | include 'INPUT'
|
|---|
| 81 | include 'TSTEP'
|
|---|
| 82 | include 'CTIMER'
|
|---|
| 83 |
|
|---|
| 84 | if (icalld.eq.0) tbmhd=0.0
|
|---|
| 85 | icalld = icalld+1
|
|---|
| 86 | nbmhd = icalld
|
|---|
| 87 | etime1 = dnekclock()
|
|---|
| 88 | c
|
|---|
| 89 | ifield = 1
|
|---|
| 90 | call makeuf
|
|---|
| 91 | c
|
|---|
| 92 | ifield = ifldmhd
|
|---|
| 93 | call makeufb
|
|---|
| 94 | if (ifaxis) then
|
|---|
| 95 | c do ifield = 2,3
|
|---|
| 96 | do ifield = 2,npscal+1
|
|---|
| 97 | call makeuq ! nonlinear terms
|
|---|
| 98 | enddo
|
|---|
| 99 | ifield = ifldmhd
|
|---|
| 100 | endif
|
|---|
| 101 |
|
|---|
| 102 | if (ifnav.and.(.not.ifchar)) then
|
|---|
| 103 | call advab_elsasser_fast
|
|---|
| 104 | endif
|
|---|
| 105 | if (ifchar) then
|
|---|
| 106 | write(6,*) 'No IFCHAR for MHD, yet.'
|
|---|
| 107 | call exitt
|
|---|
| 108 | endif
|
|---|
| 109 |
|
|---|
| 110 | ifield = 1
|
|---|
| 111 | if (iftran) call makeabf
|
|---|
| 112 | call makebdf
|
|---|
| 113 | c
|
|---|
| 114 | ifield = ifldmhd
|
|---|
| 115 | if (iftran) call makextb
|
|---|
| 116 | call makebdfb
|
|---|
| 117 | c
|
|---|
| 118 | tbmhd=tbmhd+(dnekclock()-etime1)
|
|---|
| 119 | return
|
|---|
| 120 | end
|
|---|
| 121 | c--------------------------------------------------------------------
|
|---|
| 122 | subroutine makeufb
|
|---|
| 123 | c
|
|---|
| 124 | c Compute and add: (1) user specified forcing function (FX,FY,FZ)
|
|---|
| 125 | c
|
|---|
| 126 | include 'SIZE'
|
|---|
| 127 | include 'SOLN'
|
|---|
| 128 | include 'MASS'
|
|---|
| 129 | include 'TSTEP'
|
|---|
| 130 | C
|
|---|
| 131 | time = time-dt
|
|---|
| 132 | call nekuf (bmx,bmy,bmz)
|
|---|
| 133 | call opcolv2 (bmx,bmy,bmz,vtrans(1,1,1,1,ifield),bm1)
|
|---|
| 134 | time = time+dt
|
|---|
| 135 | c
|
|---|
| 136 | return
|
|---|
| 137 | end
|
|---|
| 138 | c--------------------------------------------------------------------
|
|---|
| 139 | subroutine makextb
|
|---|
| 140 | c
|
|---|
| 141 | c Add extrapolation terms to magnetic source terms
|
|---|
| 142 | c
|
|---|
| 143 | c (nek5 equivalent for velocity is "makeabf")
|
|---|
| 144 | c
|
|---|
| 145 | include 'SIZE'
|
|---|
| 146 | include 'INPUT'
|
|---|
| 147 | include 'SOLN'
|
|---|
| 148 | include 'MASS'
|
|---|
| 149 | include 'TSTEP'
|
|---|
| 150 | C
|
|---|
| 151 | common /scrns/ ta1 (lx1,ly1,lz1,lelv)
|
|---|
| 152 | $ , ta2 (lx1,ly1,lz1,lelv)
|
|---|
| 153 | $ , ta3 (lx1,ly1,lz1,lelv)
|
|---|
| 154 | c
|
|---|
| 155 | ntot1 = lx1*ly1*lz1*nelv
|
|---|
| 156 | c
|
|---|
| 157 | ab0 = ab(1)
|
|---|
| 158 | ab1 = ab(2)
|
|---|
| 159 | ab2 = ab(3)
|
|---|
| 160 | call add3s2 (ta1,bbx1,bbx2,ab1,ab2,ntot1)
|
|---|
| 161 | call add3s2 (ta2,bby1,bby2,ab1,ab2,ntot1)
|
|---|
| 162 | call copy (bbx2,bbx1,ntot1)
|
|---|
| 163 | call copy (bby2,bby1,ntot1)
|
|---|
| 164 | call copy (bbx1,bmx,ntot1)
|
|---|
| 165 | call copy (bby1,bmy,ntot1)
|
|---|
| 166 | call add2s1 (bmx,ta1,ab0,ntot1)
|
|---|
| 167 | call add2s1 (bmy,ta2,ab0,ntot1)
|
|---|
| 168 | if (ldim.eq.3) then
|
|---|
| 169 | call add3s2 (ta3,bbz1,bbz2,ab1,ab2,ntot1)
|
|---|
| 170 | call copy (bbz2,bbz1,ntot1)
|
|---|
| 171 | call copy (bbz1,bmz,ntot1)
|
|---|
| 172 | call add2s1 (bmz,ta3,ab0,ntot1)
|
|---|
| 173 | endif
|
|---|
| 174 | c
|
|---|
| 175 | return
|
|---|
| 176 | end
|
|---|
| 177 | c--------------------------------------------------------------------
|
|---|
| 178 | subroutine makebdfb
|
|---|
| 179 | C
|
|---|
| 180 | C Add contributions to magnetic source from lagged BD terms.
|
|---|
| 181 | C
|
|---|
| 182 | include 'SIZE'
|
|---|
| 183 | include 'SOLN'
|
|---|
| 184 | include 'MASS'
|
|---|
| 185 | include 'GEOM'
|
|---|
| 186 | include 'INPUT'
|
|---|
| 187 | include 'TSTEP'
|
|---|
| 188 | C
|
|---|
| 189 | COMMON /SCRNS/ TA1(LX1,LY1,LZ1,LELV)
|
|---|
| 190 | $ , TA2(LX1,LY1,LZ1,LELV)
|
|---|
| 191 | $ , TA3(LX1,LY1,LZ1,LELV)
|
|---|
| 192 | $ , TB1(LX1,LY1,LZ1,LELV)
|
|---|
| 193 | $ , TB2(LX1,LY1,LZ1,LELV)
|
|---|
| 194 | $ , TB3(LX1,LY1,LZ1,LELV)
|
|---|
| 195 | $ , H2 (LX1,LY1,LZ1,LELV)
|
|---|
| 196 | C
|
|---|
| 197 | NTOT1 = lx1*ly1*lz1*NELV
|
|---|
| 198 | CONST = 1./DT
|
|---|
| 199 | CALL CMULT2(H2,vtrans(1,1,1,1,ifield),CONST,NTOT1)
|
|---|
| 200 | CALL OPCOLV3c (TB1,TB2,TB3,BX,BY,BZ,BM1,bd(2))
|
|---|
| 201 | C
|
|---|
| 202 | DO 100 ILAG=2,NBD
|
|---|
| 203 | IF (IFGEOM) THEN
|
|---|
| 204 | CALL OPCOLV3c(TA1,TA2,TA3,BXLAG (1,ILAG-1),
|
|---|
| 205 | $ BYLAG (1,ILAG-1),
|
|---|
| 206 | $ BZLAG (1,ILAG-1),
|
|---|
| 207 | $ BM1LAG(1,1,1,1,ILAG-1),bd(ilag+1))
|
|---|
| 208 | ELSE
|
|---|
| 209 | CALL OPCOLV3c(TA1,TA2,TA3,BXLAG (1,ILAG-1),
|
|---|
| 210 | $ BYLAG (1,ILAG-1),
|
|---|
| 211 | $ BZLAG (1,ILAG-1),
|
|---|
| 212 | $ BM1 ,bd(ilag+1))
|
|---|
| 213 | ENDIF
|
|---|
| 214 | CALL OPADD2 (TB1,TB2,TB3,TA1,TA2,TA3)
|
|---|
| 215 | 100 CONTINUE
|
|---|
| 216 | CALL OPADD2col (BMX,BMY,BMZ,TB1,TB2,TB3,h2)
|
|---|
| 217 | c
|
|---|
| 218 | return
|
|---|
| 219 | end
|
|---|
| 220 | c--------------------------------------------------------------------
|
|---|
| 221 | subroutine cresvib(resv1,resv2,resv3,h1,h2)
|
|---|
| 222 | c
|
|---|
| 223 | c Account for inhomogeneous Dirichlet boundary contributions
|
|---|
| 224 | c in rhs of induction eqn.
|
|---|
| 225 | c n
|
|---|
| 226 | c Also, subtract off best estimate of grad p
|
|---|
| 227 | c
|
|---|
| 228 | include 'SIZE'
|
|---|
| 229 | include 'TOTAL'
|
|---|
| 230 | real resv1 (lx1,ly1,lz1,1)
|
|---|
| 231 | real resv2 (lx1,ly1,lz1,1)
|
|---|
| 232 | real resv3 (lx1,ly1,lz1,1)
|
|---|
| 233 | real h1 (lx1,ly1,lz1,1)
|
|---|
| 234 | real h2 (lx1,ly1,lz1,1)
|
|---|
| 235 | common /scruz/ w1 (lx1,ly1,lz1,lelv)
|
|---|
| 236 | $ , w2 (lx1,ly1,lz1,lelv)
|
|---|
| 237 | $ , w3 (lx1,ly1,lz1,lelv)
|
|---|
| 238 | c
|
|---|
| 239 | c
|
|---|
| 240 | call bcdirvc (bx,by,bz,b1mask,b2mask,b3mask)
|
|---|
| 241 | call extrapp (pm,pmlag)
|
|---|
| 242 | call opgradt (resv1,resv2,resv3,pm)
|
|---|
| 243 | call opadd2 (resv1,resv2,resv3,bmx,bmy,bmz)
|
|---|
| 244 | c call opcopy (resv1,resv2,resv3,bmx,bmy,bmz)
|
|---|
| 245 | call ophx (w1,w2,w3,bx,by,bz,h1,h2)
|
|---|
| 246 | call opsub2 (resv1,resv2,resv3,w1,w2,w3)
|
|---|
| 247 | c
|
|---|
| 248 | return
|
|---|
| 249 | end
|
|---|
| 250 | c--------------------------------------------------------------------
|
|---|
| 251 | subroutine cresvibp(resv1,resv2,resv3,h1,h2)
|
|---|
| 252 | c
|
|---|
| 253 | c Account for inhomogeneous Dirichlet boundary contributions
|
|---|
| 254 | c in rhs of momentum eqn.
|
|---|
| 255 | c n
|
|---|
| 256 | c Also, subtract off best estimate of grad p
|
|---|
| 257 | c
|
|---|
| 258 | include 'SIZE'
|
|---|
| 259 | include 'TOTAL'
|
|---|
| 260 | real resv1 (lx1,ly1,lz1,1)
|
|---|
| 261 | real resv2 (lx1,ly1,lz1,1)
|
|---|
| 262 | real resv3 (lx1,ly1,lz1,1)
|
|---|
| 263 | real h1 (lx1,ly1,lz1,1)
|
|---|
| 264 | real h2 (lx1,ly1,lz1,1)
|
|---|
| 265 | common /scruz/ w1 (lx1,ly1,lz1,lelv)
|
|---|
| 266 | $ , w2 (lx1,ly1,lz1,lelv)
|
|---|
| 267 | $ , w3 (lx1,ly1,lz1,lelv)
|
|---|
| 268 | c
|
|---|
| 269 | call bcdirvc (vx,vy,vz,v1mask,v2mask,v3mask)
|
|---|
| 270 | if (ifstrs) call bcneutr
|
|---|
| 271 | call extrapp (pr,prlag)
|
|---|
| 272 | call opgradt (resv1,resv2,resv3,pr)
|
|---|
| 273 | call opadd2 (resv1,resv2,resv3,bfx,bfy,bfz)
|
|---|
| 274 | c call opcopy (resv1,resv2,resv3,bfx,bfy,bfz)
|
|---|
| 275 | call ophx (w1,w2,w3,vx,vy,vz,h1,h2)
|
|---|
| 276 | call opsub2 (resv1,resv2,resv3,w1,w2,w3)
|
|---|
| 277 | c
|
|---|
| 278 | return
|
|---|
| 279 | end
|
|---|
| 280 | c--------------------------------------------------------------------
|
|---|
| 281 | subroutine incomprn (ux,uy,uz,up)
|
|---|
| 282 | c
|
|---|
| 283 | c Project U onto the closest incompressible field
|
|---|
| 284 | c
|
|---|
| 285 | c Input: U := (ux,uy,uz)
|
|---|
| 286 | c
|
|---|
| 287 | c Output: updated values of U, iproj, proj; and
|
|---|
| 288 | c up := pressure currection req'd to impose div U = 0
|
|---|
| 289 | c
|
|---|
| 290 | c
|
|---|
| 291 | c Dependencies: ifield ==> which "density" (vtrans) is used.
|
|---|
| 292 | c
|
|---|
| 293 | c Notes 1. up is _not_ scaled by bd(1)/dt. This should be done
|
|---|
| 294 | c external to incompr().
|
|---|
| 295 | c
|
|---|
| 296 | c 2. up accounts _only_ for the perturbation pressure,
|
|---|
| 297 | c not the current pressure derived from extrapolation.
|
|---|
| 298 | c
|
|---|
| 299 | c
|
|---|
| 300 | include 'SIZE'
|
|---|
| 301 | include 'TOTAL'
|
|---|
| 302 | include 'CTIMER'
|
|---|
| 303 | c
|
|---|
| 304 | common /scrns/ w1 (lx1,ly1,lz1,lelv)
|
|---|
| 305 | $ , w2 (lx1,ly1,lz1,lelv)
|
|---|
| 306 | $ , w3 (lx1,ly1,lz1,lelv)
|
|---|
| 307 | $ , dv1 (lx1,ly1,lz1,lelv)
|
|---|
| 308 | $ , dv2 (lx1,ly1,lz1,lelv)
|
|---|
| 309 | $ , dv3 (lx1,ly1,lz1,lelv)
|
|---|
| 310 | $ , dp (lx2,ly2,lz2,lelv)
|
|---|
| 311 | common /scrvh/ h1 (lx1,ly1,lz1,lelv)
|
|---|
| 312 | $ , h2 (lx1,ly1,lz1,lelv)
|
|---|
| 313 | common /scrhi/ h2inv (lx1,ly1,lz1,lelv)
|
|---|
| 314 |
|
|---|
| 315 | parameter(nset = 1 + lbelv/lelv)
|
|---|
| 316 | common /orthov/ pset(lx2*ly2*lz2*lelv*mxprev,nset)
|
|---|
| 317 | common /orthbi/ nprv(2)
|
|---|
| 318 | logical ifprjp
|
|---|
| 319 |
|
|---|
| 320 | ifprjp=.false. ! Project out previous pressure solutions?
|
|---|
| 321 | istart=param(95)
|
|---|
| 322 | if (istep.ge.istart.and.istart.ne.0) ifprjp=.true.
|
|---|
| 323 |
|
|---|
| 324 | if (icalld.eq.0) tpres=0.0
|
|---|
| 325 | icalld = icalld+1
|
|---|
| 326 | npres = icalld
|
|---|
| 327 | etime1 = dnekclock()
|
|---|
| 328 |
|
|---|
| 329 | ntot1 = lx1*ly1*lz1*nelv
|
|---|
| 330 | ntot2 = lx2*ly2*lz2*nelv
|
|---|
| 331 | intype = 1
|
|---|
| 332 |
|
|---|
| 333 | call rzero (h1,ntot1)
|
|---|
| 334 | call copy (h2,vtrans(1,1,1,1,ifield),ntot1)
|
|---|
| 335 | call invers2 (h2inv,h2,ntot1)
|
|---|
| 336 |
|
|---|
| 337 | call opdiv (dp,ux,uy,uz)
|
|---|
| 338 |
|
|---|
| 339 | bdti = -bd(1)/dt
|
|---|
| 340 | call cmult (dp,bdti,ntot2)
|
|---|
| 341 |
|
|---|
| 342 | call add2col2(dp,bm2,usrdiv,ntot2) ! User-defined divergence.
|
|---|
| 343 |
|
|---|
| 344 | call ortho (dp)
|
|---|
| 345 |
|
|---|
| 346 | i = 1 + ifield/ifldmhd
|
|---|
| 347 | if (ifprjp) call setrhsp (dp,h1,h2,h2inv,pset(1,i),nprv(i))
|
|---|
| 348 | scaledt = dt/bd(1)
|
|---|
| 349 | scaledi = 1./scaledt
|
|---|
| 350 | call cmult(dp,scaledt,ntot2) ! scale for tol
|
|---|
| 351 | call esolver (dp,h1,h2,h2inv,intype)
|
|---|
| 352 | call cmult(dp,scaledi,ntot2)
|
|---|
| 353 | if (ifprjp) call gensolnp (dp,h1,h2,h2inv,pset(1,i),nprv(i))
|
|---|
| 354 |
|
|---|
| 355 | call add2(up,dp,ntot2)
|
|---|
| 356 |
|
|---|
| 357 | call opgradt (w1 ,w2 ,w3 ,dp)
|
|---|
| 358 | call opbinv (dv1,dv2,dv3,w1 ,w2 ,w3 ,h2inv)
|
|---|
| 359 | dtb = dt/bd(1)
|
|---|
| 360 | call opadd2cm (ux ,uy ,uz ,dv1,dv2,dv3, dtb )
|
|---|
| 361 |
|
|---|
| 362 | if (ifmhd) call chkptol ! to avoid repetition
|
|---|
| 363 |
|
|---|
| 364 | tpres=tpres+(dnekclock()-etime1)
|
|---|
| 365 |
|
|---|
| 366 | return
|
|---|
| 367 | end
|
|---|
| 368 | c-----------------------------------------------------------------------
|
|---|
| 369 | subroutine extrapp(p,plag)
|
|---|
| 370 | C
|
|---|
| 371 | C Pressure extrapolation
|
|---|
| 372 | C
|
|---|
| 373 | INCLUDE 'SIZE'
|
|---|
| 374 | INCLUDE 'SOLN'
|
|---|
| 375 | INCLUDE 'TSTEP'
|
|---|
| 376 |
|
|---|
| 377 | real p (lx2,ly2,lz2,1)
|
|---|
| 378 | $ ,plag (lx2,ly2,lz2,1)
|
|---|
| 379 |
|
|---|
| 380 | common /cgeom/ igeom
|
|---|
| 381 |
|
|---|
| 382 | ntot2 = lx2*ly2*lz2*nelv
|
|---|
| 383 |
|
|---|
| 384 | if (nbd.eq.2.and.nbdinp.gt.2.and.igeom.le.2) then
|
|---|
| 385 | call copy(plag,p,ntot2)
|
|---|
| 386 | elseif (nbd.eq.3.and.igeom.le.2) then
|
|---|
| 387 |
|
|---|
| 388 | const = dtlag(1)/dtlag(2)
|
|---|
| 389 |
|
|---|
| 390 | do i=1,ntot2
|
|---|
| 391 | pnm1 = p (i,1,1,1)
|
|---|
| 392 | pnm2 = plag(i,1,1,1)
|
|---|
| 393 | p (i,1,1,1) = pnm1 + const*(pnm1-pnm2)
|
|---|
| 394 | plag(i,1,1,1) = pnm1
|
|---|
| 395 | enddo
|
|---|
| 396 |
|
|---|
| 397 | elseif (nbd.gt.3) then
|
|---|
| 398 | WRITE (6,*) 'Pressure extrapolation cannot be completed'
|
|---|
| 399 | WRITE (6,*) 'Try a lower-order temporal scheme'
|
|---|
| 400 | call exitt
|
|---|
| 401 | endif
|
|---|
| 402 | return
|
|---|
| 403 | end
|
|---|
| 404 | c-----------------------------------------------------------------------
|
|---|
| 405 | subroutine opzero(ux,uy,uz)
|
|---|
| 406 | include 'SIZE'
|
|---|
| 407 | include 'TOTAL'
|
|---|
| 408 | real ux(1),uy(1),uz(1)
|
|---|
| 409 | c
|
|---|
| 410 | n = lx1*ly1*lz1*nelfld(ifield)
|
|---|
| 411 | call rzero(ux,n)
|
|---|
| 412 | call rzero(uy,n)
|
|---|
| 413 | if (if3d) call rzero(uz,n)
|
|---|
| 414 | c
|
|---|
| 415 | return
|
|---|
| 416 | end
|
|---|
| 417 | c-----------------------------------------------------------------------
|
|---|
| 418 | subroutine opnorm(unorm,ux,uy,uz,type3)
|
|---|
| 419 | include 'SIZE'
|
|---|
| 420 | include 'TOTAL'
|
|---|
| 421 | character*3 type3
|
|---|
| 422 | real ux(1),uy(1),uz(1)
|
|---|
| 423 | real un(3),wn(3)
|
|---|
| 424 | c
|
|---|
| 425 | n = lx1*ly1*lz1*nelfld(ifield)
|
|---|
| 426 | if (type3.eq.'L2 ') then
|
|---|
| 427 | if (if3d) then
|
|---|
| 428 | un(1) = vlsc3(ux,ux,bm1,n)
|
|---|
| 429 | un(2) = vlsc3(uy,uy,bm1,n)
|
|---|
| 430 | un(3) = vlsc3(uz,uz,bm1,n)
|
|---|
| 431 | un(1) = un(1) + un(2) + un(3)
|
|---|
| 432 | unorm = glsum(un(1),1)
|
|---|
| 433 | if (unorm.gt.0) unorm = sqrt(unorm/volvm1)
|
|---|
| 434 | else
|
|---|
| 435 | un(1) = vlsc3(ux,ux,bm1,n)
|
|---|
| 436 | un(2) = vlsc3(uy,uy,bm1,n)
|
|---|
| 437 | un(1) = un(1) + un(2)
|
|---|
| 438 | unorm = glsum(un(1),1)
|
|---|
| 439 | if (unorm.gt.0) unorm = sqrt(unorm/volvm1)
|
|---|
| 440 | endif
|
|---|
| 441 | endif
|
|---|
| 442 | c
|
|---|
| 443 | return
|
|---|
| 444 | end
|
|---|
| 445 | c-----------------------------------------------------------------------
|
|---|
| 446 | subroutine lorentz_force (lf,b1,b2,b3,w1,w2)
|
|---|
| 447 | c
|
|---|
| 448 | c Compute Lorentz force
|
|---|
| 449 | c
|
|---|
| 450 | c Input: B := (b1,b2,b3)
|
|---|
| 451 | c
|
|---|
| 452 | c Output: lf(1,ldim)
|
|---|
| 453 | c
|
|---|
| 454 | c Work arrays: w1(ltot) and w2(ltot)
|
|---|
| 455 | c
|
|---|
| 456 | c The output will not be continuous. However, it will be in
|
|---|
| 457 | c the form appropriate for incorporation as a body force term
|
|---|
| 458 | c in the variational formulation of the Navier-Stokes equations.
|
|---|
| 459 | c
|
|---|
| 460 | c (i.e., rhs(NS) = rhs(NS) + B*lf, where B is the mass matrix)
|
|---|
| 461 | c
|
|---|
| 462 | include 'SIZE'
|
|---|
| 463 | c
|
|---|
| 464 | real lf(lx1*ly1*lz1*lelv,ldim)
|
|---|
| 465 | real b1(lx1*ly1*lz1*lelv)
|
|---|
| 466 | real b2(lx1*ly1*lz1*lelv)
|
|---|
| 467 | real b3(lx1*ly1*lz1*lelv)
|
|---|
| 468 | c
|
|---|
| 469 | call curl(lf,b1,b2,b3,.false.,w1,w2)
|
|---|
| 470 | c
|
|---|
| 471 | ntot = lx1*ly1*lz1*nelv
|
|---|
| 472 | c
|
|---|
| 473 | do i=1,ntot
|
|---|
| 474 | c1 = lf(i,2)*b3(i) - lf(i,3)*b2(i)
|
|---|
| 475 | c2 = lf(i,3)*b1(i) - lf(i,1)*b3(i)
|
|---|
| 476 | c3 = lf(i,1)*b2(i) - lf(i,2)*b1(i)
|
|---|
| 477 | lf(i,1) = c1
|
|---|
| 478 | lf(i,2) = c2
|
|---|
| 479 | lf(i,3) = c3
|
|---|
| 480 | enddo
|
|---|
| 481 | c
|
|---|
| 482 | return
|
|---|
| 483 | end
|
|---|
| 484 | c-----------------------------------------------------------------------
|
|---|
| 485 | subroutine curl(vort,u,v,w,ifavg,work1,work2)
|
|---|
| 486 | c
|
|---|
| 487 | include 'SIZE'
|
|---|
| 488 | include 'TOTAL'
|
|---|
| 489 | c
|
|---|
| 490 | logical ifavg
|
|---|
| 491 | c
|
|---|
| 492 | parameter(lt=lx1*ly1*lz1*lelv)
|
|---|
| 493 | real vort(lt,3),work1(1),work2(1),u(1),v(1),w(1)
|
|---|
| 494 | c
|
|---|
| 495 | ntot = lx1*ly1*lz1*nelv
|
|---|
| 496 | if (if3d) then
|
|---|
| 497 | c work1=dw/dy ; work2=dv/dz
|
|---|
| 498 | call dudxyz(work1,w,rym1,sym1,tym1,jacm1,1,2)
|
|---|
| 499 | call dudxyz(work2,v,rzm1,szm1,tzm1,jacm1,1,3)
|
|---|
| 500 | call sub3(vort(1,1),work1,work2,ntot)
|
|---|
| 501 | c work1=du/dz ; work2=dw/dx
|
|---|
| 502 | call dudxyz(work1,u,rzm1,szm1,tzm1,jacm1,1,3)
|
|---|
| 503 | call dudxyz(work2,w,rxm1,sxm1,txm1,jacm1,1,1)
|
|---|
| 504 | call sub3(vort(1,2),work1,work2,ntot)
|
|---|
| 505 | c work1=dv/dx ; work2=du/dy
|
|---|
| 506 | call dudxyz(work1,v,rxm1,sxm1,txm1,jacm1,1,1)
|
|---|
| 507 | call dudxyz(work2,u,rym1,sym1,tym1,jacm1,1,2)
|
|---|
| 508 | call sub3(vort(1,3),work1,work2,ntot)
|
|---|
| 509 | else
|
|---|
| 510 | c work1=dv/dx ; work2=du/dy
|
|---|
| 511 | call dudxyz(work1,v,rxm1,sxm1,txm1,jacm1,1,1)
|
|---|
| 512 | call dudxyz(work2,u,rym1,sym1,tym1,jacm1,1,2)
|
|---|
| 513 | call sub3(vort(1,3),work1,work2,ntot)
|
|---|
| 514 | endif
|
|---|
| 515 | c
|
|---|
| 516 | c Avg at bndry
|
|---|
| 517 | c
|
|---|
| 518 | if (ifavg) then
|
|---|
| 519 | ifielt = ifield
|
|---|
| 520 | ifield = 1
|
|---|
| 521 | if (if3d) then
|
|---|
| 522 | do idim=1,ldim
|
|---|
| 523 | call col2 (vort(1,idim),bm1,ntot)
|
|---|
| 524 | call dssum (vort(1,idim),lx1,ly1,lz1)
|
|---|
| 525 | call col2 (vort(1,idim),binvm1,ntot)
|
|---|
| 526 | enddo
|
|---|
| 527 | else
|
|---|
| 528 | call col2 (vort(1,3),bm1,ntot) ! NOTE: This differs from
|
|---|
| 529 | call dssum (vort(1,3),lx1,ly1,lz1) ! "comp_vort", which returns
|
|---|
| 530 | call col2 (vort(1,3),binvm1,ntot) ! vorticity as 1st entry in vort
|
|---|
| 531 | endif
|
|---|
| 532 | ifield = ifielt
|
|---|
| 533 | endif
|
|---|
| 534 | c
|
|---|
| 535 | return
|
|---|
| 536 | end
|
|---|
| 537 | c-----------------------------------------------------------------------
|
|---|
| 538 | subroutine lorentz_force2(lf,b1,b2,b3)
|
|---|
| 539 | c
|
|---|
| 540 | c Compute Lorentz force
|
|---|
| 541 | c
|
|---|
| 542 | c Input: B := (b1,b2,b3)
|
|---|
| 543 | c
|
|---|
| 544 | c Output: lf(1,ldim)
|
|---|
| 545 | c
|
|---|
| 546 | c Work arrays: w1(ltot) and w2(ltot)
|
|---|
| 547 | c
|
|---|
| 548 | c The output will not be continuous. However, it will be in
|
|---|
| 549 | c the form appropriate for incorporation as a body force term
|
|---|
| 550 | c in the variational formulation of the Navier-Stokes equations.
|
|---|
| 551 | c
|
|---|
| 552 | c (i.e., rhs(NS) = rhs(NS) + B*lf, where B is the mass matrix)
|
|---|
| 553 | c
|
|---|
| 554 | include 'SIZE'
|
|---|
| 555 | c
|
|---|
| 556 | real lf(lx1*ly1*lz1*ldim,lelt)
|
|---|
| 557 | real b1(lx1*ly1*lz1,lelt)
|
|---|
| 558 | real b2(lx1*ly1*lz1,lelt)
|
|---|
| 559 | real b3(lx1*ly1*lz1,lelt)
|
|---|
| 560 | c
|
|---|
| 561 | integer e
|
|---|
| 562 | c
|
|---|
| 563 | do e = 1,nelt ! NOTE: the order is different from v. 1
|
|---|
| 564 | call lorentz_force_e(lf(1,e),b1(1,e),b2(1,e),b3(1,e),e)
|
|---|
| 565 | enddo
|
|---|
| 566 | c
|
|---|
| 567 | return
|
|---|
| 568 | end
|
|---|
| 569 | c-----------------------------------------------------------------------
|
|---|
| 570 | subroutine lorentz_force_e(lf,b1,b2,b3,e)
|
|---|
| 571 | c
|
|---|
| 572 | c Compute Lorentz force field for a single element
|
|---|
| 573 | c
|
|---|
| 574 | c Input: B := (b1,b2,b3)
|
|---|
| 575 | c
|
|---|
| 576 | c Output: lf(1,ldim)
|
|---|
| 577 | c
|
|---|
| 578 | c Work arrays: cb(lxyzd) and w2(lxyzd)
|
|---|
| 579 | c
|
|---|
| 580 | c The output will not be continuous. However, it will be in
|
|---|
| 581 | c the form appropriate for incorporation as a body force term
|
|---|
| 582 | c in the variational formulation of the Navier-Stokes equations.
|
|---|
| 583 | c (i.e., rhs(NS) = rhs(NS) + B*lf, where B is the mass matrix)
|
|---|
| 584 | c
|
|---|
| 585 | c Dealiasing strategy:
|
|---|
| 586 | c
|
|---|
| 587 | c e e -1 ~ T ~e ~ ~
|
|---|
| 588 | c lf = ( B ) I B ( I j x I b )
|
|---|
| 589 | c i ~ i
|
|---|
| 590 | c
|
|---|
| 591 | c ~
|
|---|
| 592 | c Here, I is the interpolant from N to M, where M = 3/2 N.
|
|---|
| 593 | c
|
|---|
| 594 | c ~ -1
|
|---|
| 595 | c j is a special curl (sans Jacobian )
|
|---|
| 596 | c
|
|---|
| 597 | c b is the B-field on the N-points
|
|---|
| 598 | c
|
|---|
| 599 | c ~e
|
|---|
| 600 | c B is the local mass matrix on the M-points (sans Jacabian)
|
|---|
| 601 | c
|
|---|
| 602 | c e
|
|---|
| 603 | c B is the local mass matrix on the N-points (with Jacabian)
|
|---|
| 604 | c
|
|---|
| 605 | c The last B is req'd to compensate for the subsequent multiplication
|
|---|
| 606 | c by B that takes place once the rhs is formed.
|
|---|
| 607 | c
|
|---|
| 608 | c
|
|---|
| 609 | c
|
|---|
| 610 | include 'SIZE'
|
|---|
| 611 | include 'DEALIAS'
|
|---|
| 612 | include 'GEOM'
|
|---|
| 613 | include 'WZ'
|
|---|
| 614 | c
|
|---|
| 615 | real lf(lx1*ly1*lz1,3)
|
|---|
| 616 | real b1(lx1*ly1*lz1)
|
|---|
| 617 | real b2(lx1*ly1*lz1)
|
|---|
| 618 | real b3(lx1*ly1*lz1)
|
|---|
| 619 | integer d,e
|
|---|
| 620 | c
|
|---|
| 621 | integer icalld
|
|---|
| 622 | save icalld
|
|---|
| 623 | data icalld /0/
|
|---|
| 624 | c
|
|---|
| 625 | common /ctmp1x/ lfd(lxd*lyd*lzd,3)
|
|---|
| 626 | $ , bd (lxd*lyd*lzd,3)
|
|---|
| 627 | $ , cb (lx1*ly1*lz1,3)
|
|---|
| 628 | $ , cbd(lxd*lyd*lzd,3)
|
|---|
| 629 | real lfd
|
|---|
| 630 |
|
|---|
| 631 | if (icalld .eq. 0) then
|
|---|
| 632 | write(6,*) 'CALL SET PROJ',lx1,lxd
|
|---|
| 633 | call setmap (lx1,lxd) ! Set up interpolation operators
|
|---|
| 634 | call setproj(lx1,lxd) ! Set up interpolation operators
|
|---|
| 635 | icalld = icalld + 1
|
|---|
| 636 | endif
|
|---|
| 637 | c
|
|---|
| 638 | call spec_curl_e (cb,b1,b2,b3 !local curl, w/o Jacobian
|
|---|
| 639 | $ , rxm1(1,1,1,e),rym1(1,1,1,e),rzm1(1,1,1,e)
|
|---|
| 640 | $ , sxm1(1,1,1,e),sym1(1,1,1,e),szm1(1,1,1,e)
|
|---|
| 641 | $ , txm1(1,1,1,e),tym1(1,1,1,e),tzm1(1,1,1,e) )
|
|---|
| 642 | c
|
|---|
| 643 | do d=1,ldim ! interpolate to M points
|
|---|
| 644 | call specmp(cbd(1,d),lxd,cb(1,d),lx1,im1d,im1dt,bd)
|
|---|
| 645 | enddo
|
|---|
| 646 | call specmp(bd(1,1),lxd,b1,lx1,im1d,im1dt,lfd)
|
|---|
| 647 | call specmp(bd(1,2),lxd,b2,lx1,im1d,im1dt,lfd)
|
|---|
| 648 | call specmp(bd(1,3),lxd,b3,lx1,im1d,im1dt,lfd)
|
|---|
| 649 | c
|
|---|
| 650 | nxyzd = lxd*lyd*lzd
|
|---|
| 651 | do i=1,nxyzd
|
|---|
| 652 | lfd(i,1) = cbd(i,2)*bd(i,3) - cbd(i,3)*bd(i,2) ! Curl B x B
|
|---|
| 653 | lfd(i,2) = cbd(i,3)*bd(i,1) - cbd(i,1)*bd(i,3)
|
|---|
| 654 | lfd(i,3) = cbd(i,1)*bd(i,2) - cbd(i,2)*bd(i,1)
|
|---|
| 655 | enddo
|
|---|
| 656 | c
|
|---|
| 657 | c Project back and simultaneous collocate with local quadrature weights
|
|---|
| 658 | c
|
|---|
| 659 | c ~ ~ ~
|
|---|
| 660 | c P := Pz x Py x Px = Iz*B x Iy*B x Ix*B
|
|---|
| 661 | c
|
|---|
| 662 | c
|
|---|
| 663 | c ~ M
|
|---|
| 664 | c where B = diag (rho )
|
|---|
| 665 | c i
|
|---|
| 666 | c
|
|---|
| 667 | call specmp(lf(1,1),lx1,lfd(1,1),lxd,pmd1,pmd1t,cbd)
|
|---|
| 668 | call specmp(lf(1,2),lx1,lfd(1,2),lxd,pmd1,pmd1t,cbd)
|
|---|
| 669 | call specmp(lf(1,3),lx1,lfd(1,3),lxd,pmd1,pmd1t,cbd)
|
|---|
| 670 | c
|
|---|
| 671 | c Finally, divide by local mass matrix in anticipation of subsequent
|
|---|
| 672 | c multiply by BM1.
|
|---|
| 673 | c
|
|---|
| 674 | nxyz = lx1*ly1*lz1
|
|---|
| 675 | do i=1,nxyz
|
|---|
| 676 | c scale = 1./(w3m1(i,1,1)*jacm1(i,1,1,e))
|
|---|
| 677 | scale = 1./jacm1(i,1,1,e)
|
|---|
| 678 | lf(i,1) = scale*lf(i,1)
|
|---|
| 679 | lf(i,2) = scale*lf(i,2)
|
|---|
| 680 | lf(i,3) = scale*lf(i,3)
|
|---|
| 681 | enddo
|
|---|
| 682 | c
|
|---|
| 683 | return
|
|---|
| 684 | end
|
|---|
| 685 | c-----------------------------------------------------------------------
|
|---|
| 686 | subroutine spec_curl_e (cb,b1,b2,b3,rx,ry,rz,sx,sy,sz,tx,ty,tz)
|
|---|
| 687 | c
|
|---|
| 688 | c local curl, multiplied by Jacobian
|
|---|
| 689 | c
|
|---|
| 690 | include 'SIZE'
|
|---|
| 691 | include 'DXYZ'
|
|---|
| 692 | c
|
|---|
| 693 | real cb(lx1*ly1*lz1,3) ! Output J*curl B (J:=Jacobian)
|
|---|
| 694 | real b1(1),b2(1),b3(1) ! Input B-field
|
|---|
| 695 | c
|
|---|
| 696 | real rx(1),ry(1),rz(1) ! Metrics
|
|---|
| 697 | real sx(1),sy(1),sz(1)
|
|---|
| 698 | real tx(1),ty(1),tz(1)
|
|---|
| 699 | c
|
|---|
| 700 | common /ctmp0x/ br(lx1*ly1*lz1),bs(lx1*ly1*lz1),bt(lx1*ly1*lz1)
|
|---|
| 701 | c
|
|---|
| 702 | c / db3 db2 \
|
|---|
| 703 | c cb1 = J | --- - --- | ! Keep J ( J:=Jacobian)
|
|---|
| 704 | c \ dy dz /
|
|---|
| 705 | c
|
|---|
| 706 | c
|
|---|
| 707 | c / db1 db3 \
|
|---|
| 708 | c cb2 = J | --- - --- | ! Keep J ( J:=Jacobian)
|
|---|
| 709 | c \ dz dx /
|
|---|
| 710 | c
|
|---|
| 711 | c
|
|---|
| 712 | c / db2 db1 \
|
|---|
| 713 | c cb3 = J | --- - --- | ! Keep J ( J:=Jacobian)
|
|---|
| 714 | c \ dx dy /
|
|---|
| 715 | c
|
|---|
| 716 | c
|
|---|
| 717 | c Note:
|
|---|
| 718 | c
|
|---|
| 719 | c db2 db2 dr db2 ds db2 dt
|
|---|
| 720 | c J --- = --- J -- + --- J -- + --- J --
|
|---|
| 721 | c dz dr dz ds dz dt dz
|
|---|
| 722 | c
|
|---|
| 723 | c etc.
|
|---|
| 724 | c
|
|---|
| 725 | c
|
|---|
| 726 | c
|
|---|
| 727 | nxyz = lx1*ly1*lz1
|
|---|
| 728 | c
|
|---|
| 729 | N=lx1-1
|
|---|
| 730 | call local_grad3(br,bs,bt,b1,N,1,dxm1,dxtm1)
|
|---|
| 731 | do i=1,nxyz
|
|---|
| 732 | cb(i,2) = (br(i)*rz(i)+bs(i)*sz(i)+bt(i)*tz(i))
|
|---|
| 733 | cb(i,3) = -(br(i)*ry(i)+bs(i)*sy(i)+bt(i)*ty(i))
|
|---|
| 734 | enddo
|
|---|
| 735 | c
|
|---|
| 736 | call local_grad3(br,bs,bt,b2,N,1,dxm1,dxtm1)
|
|---|
| 737 | do i=1,nxyz
|
|---|
| 738 | cb(i,1) = -(br(i)*rz(i)+bs(i)*sz(i)+bt(i)*tz(i))
|
|---|
| 739 | cb(i,3) = (br(i)*rx(i)+bs(i)*sx(i)+bt(i)*tx(i))
|
|---|
| 740 | $ + cb(i,3)
|
|---|
| 741 | enddo
|
|---|
| 742 | c
|
|---|
| 743 | call local_grad3(br,bs,bt,b3,N,1,dxm1,dxtm1)
|
|---|
| 744 | do i=1,nxyz
|
|---|
| 745 | cb(i,1) = (br(i)*ry(i)+bs(i)*sy(i)+bt(i)*ty(i))
|
|---|
| 746 | $ + cb(i,1)
|
|---|
| 747 | cb(i,2) = -(br(i)*rx(i)+bs(i)*sx(i)+bt(i)*tx(i))
|
|---|
| 748 | $ + cb(i,2)
|
|---|
| 749 | enddo
|
|---|
| 750 | c
|
|---|
| 751 | return
|
|---|
| 752 | end
|
|---|
| 753 | c-----------------------------------------------------------------------
|
|---|
| 754 | subroutine specx(b,nb,a,na,ba,ab,w)
|
|---|
| 755 | c
|
|---|
| 756 | include 'SIZE'
|
|---|
| 757 | include 'INPUT'
|
|---|
| 758 | real b(1),a(1)
|
|---|
| 759 | real w(1)
|
|---|
| 760 | c
|
|---|
| 761 | n=na*na*na
|
|---|
| 762 | do i=1,n
|
|---|
| 763 | b(i) = a(i)
|
|---|
| 764 | enddo
|
|---|
| 765 | c
|
|---|
| 766 | return
|
|---|
| 767 | end
|
|---|
| 768 | c-----------------------------------------------------------------------
|
|---|
| 769 | subroutine phys_to_elsasser(u1,u2,u3,b1,b2,b3,n)
|
|---|
| 770 | c
|
|---|
| 771 | real u1(1),u2(1),u3(1),b1(1),b2(1),b3(1)
|
|---|
| 772 | c
|
|---|
| 773 | do i=1,n
|
|---|
| 774 | c
|
|---|
| 775 | zpx = u1(i) + b1(i)
|
|---|
| 776 | zpy = u2(i) + b2(i)
|
|---|
| 777 | zpz = u3(i) + b3(i)
|
|---|
| 778 | c
|
|---|
| 779 | zmx = u1(i) - b1(i)
|
|---|
| 780 | zmy = u2(i) - b2(i)
|
|---|
| 781 | zmz = u3(i) - b3(i)
|
|---|
| 782 | c
|
|---|
| 783 | u1(i) = zpx
|
|---|
| 784 | u2(i) = zpy
|
|---|
| 785 | u3(i) = zpz
|
|---|
| 786 | c
|
|---|
| 787 | b1(i) = zmx
|
|---|
| 788 | b2(i) = zmy
|
|---|
| 789 | b3(i) = zmz
|
|---|
| 790 | c
|
|---|
| 791 | enddo
|
|---|
| 792 | c
|
|---|
| 793 | return
|
|---|
| 794 | end
|
|---|
| 795 | c-----------------------------------------------------------------------
|
|---|
| 796 | subroutine elsasser_to_phys(u1,u2,u3,b1,b2,b3,n)
|
|---|
| 797 | c
|
|---|
| 798 | real u1(1),u2(1),u3(1),b1(1),b2(1),b3(1)
|
|---|
| 799 | c
|
|---|
| 800 | do i=1,n
|
|---|
| 801 | c
|
|---|
| 802 | zpx = 0.5*( u1(i) + b1(i) )
|
|---|
| 803 | zpy = 0.5*( u2(i) + b2(i) )
|
|---|
| 804 | zpz = 0.5*( u3(i) + b3(i) )
|
|---|
| 805 | c
|
|---|
| 806 | zmx = 0.5*( u1(i) - b1(i) )
|
|---|
| 807 | zmy = 0.5*( u2(i) - b2(i) )
|
|---|
| 808 | zmz = 0.5*( u3(i) - b3(i) )
|
|---|
| 809 | c
|
|---|
| 810 | u1(i) = zpx
|
|---|
| 811 | u2(i) = zpy
|
|---|
| 812 | u3(i) = zpz
|
|---|
| 813 | c
|
|---|
| 814 | b1(i) = zmx
|
|---|
| 815 | b2(i) = zmy
|
|---|
| 816 | b3(i) = zmz
|
|---|
| 817 | c
|
|---|
| 818 | enddo
|
|---|
| 819 | c
|
|---|
| 820 | return
|
|---|
| 821 | end
|
|---|
| 822 | c-----------------------------------------------------------------------
|
|---|
| 823 | subroutine phys_to_elsasser2(u1,b1,n)
|
|---|
| 824 | c
|
|---|
| 825 | real u1(1),b1(1)
|
|---|
| 826 | c
|
|---|
| 827 | do i=1,n
|
|---|
| 828 | zpx = u1(i) + b1(i)
|
|---|
| 829 | zmx = u1(i) - b1(i)
|
|---|
| 830 | u1(i) = zpx
|
|---|
| 831 | b1(i) = zmx
|
|---|
| 832 | enddo
|
|---|
| 833 | c
|
|---|
| 834 | return
|
|---|
| 835 | end
|
|---|
| 836 | c-----------------------------------------------------------------------
|
|---|
| 837 | subroutine elsasser_to_phys2(u1,b1,n)
|
|---|
| 838 | c
|
|---|
| 839 | real u1(1),b1(1)
|
|---|
| 840 | c
|
|---|
| 841 | do i=1,n
|
|---|
| 842 | zpx = 0.5*( u1(i) + b1(i) )
|
|---|
| 843 | zmx = 0.5*( u1(i) - b1(i) )
|
|---|
| 844 | u1(i) = zpx
|
|---|
| 845 | b1(i) = zmx
|
|---|
| 846 | enddo
|
|---|
| 847 | c
|
|---|
| 848 | return
|
|---|
| 849 | end
|
|---|
| 850 | c-----------------------------------------------------------------------
|
|---|
| 851 | subroutine elsasserh(igeom)
|
|---|
| 852 | c
|
|---|
| 853 | c
|
|---|
| 854 | c Solve MHD in Elsasser variables
|
|---|
| 855 | c
|
|---|
| 856 | c
|
|---|
| 857 | include 'SIZE'
|
|---|
| 858 | include 'INPUT'
|
|---|
| 859 | include 'EIGEN'
|
|---|
| 860 | include 'SOLN'
|
|---|
| 861 | include 'TSTEP'
|
|---|
| 862 | include 'MASS'
|
|---|
| 863 | include 'GEOM'
|
|---|
| 864 | C
|
|---|
| 865 | common /scrnt/ besv1 (lbx1,lby1,lbz1,lbelv)
|
|---|
| 866 | $ , besv2 (lbx1,lby1,lbz1,lbelv)
|
|---|
| 867 | $ , besv3 (lbx1,lby1,lbz1,lbelv)
|
|---|
| 868 | COMMON /SCRNS/ RESV1 (LX1,LY1,LZ1,LELV)
|
|---|
| 869 | $ , RESV2 (LX1,LY1,LZ1,LELV)
|
|---|
| 870 | $ , RESV3 (LX1,LY1,LZ1,LELV)
|
|---|
| 871 | $ , DV1 (LX1,LY1,LZ1,LELV)
|
|---|
| 872 | $ , DV2 (LX1,LY1,LZ1,LELV)
|
|---|
| 873 | $ , DV3 (LX1,LY1,LZ1,LELV)
|
|---|
| 874 | COMMON /SCRVH/ H1 (LX1,LY1,LZ1,LELV)
|
|---|
| 875 | $ , H2 (LX1,LY1,LZ1,LELV)
|
|---|
| 876 | c
|
|---|
| 877 | n = lx1*ly1*lz1*nelv
|
|---|
| 878 | c
|
|---|
| 879 | c New geometry, new velocity
|
|---|
| 880 | c
|
|---|
| 881 | intype = -1
|
|---|
| 882 | c
|
|---|
| 883 | ifield = 1
|
|---|
| 884 | call sethlm (h1,h2,intype)
|
|---|
| 885 | call cresvibp (resv1,resv2,resv3,h1,h2)
|
|---|
| 886 | c
|
|---|
| 887 | ifield = ifldmhd
|
|---|
| 888 | call sethlm (h1,h2,intype)
|
|---|
| 889 | call cresvib (besv1,besv2,besv3,h1,h2)
|
|---|
| 890 | c
|
|---|
| 891 | c
|
|---|
| 892 | ifield = 1
|
|---|
| 893 | call sethlm (h1,h2,intype)
|
|---|
| 894 |
|
|---|
| 895 | call ophinv (dv1,dv2,dv3,resv1,resv2,resv3,h1,h2,tolhv,nmxv)
|
|---|
| 896 |
|
|---|
| 897 | call opadd2 (vx,vy,vz,dv1,dv2,dv3)
|
|---|
| 898 |
|
|---|
| 899 | if (filterType.eq.1) then
|
|---|
| 900 | alpha_filt=param(103) ! Optional Filtering
|
|---|
| 901 | call q_filter(alpha_filt)
|
|---|
| 902 | endif
|
|---|
| 903 |
|
|---|
| 904 | call incomprn (vx,vy,vz,pr) ! project U onto div-free space
|
|---|
| 905 |
|
|---|
| 906 | ifield = ifldmhd
|
|---|
| 907 | call sethlm (h1,h2,intype)
|
|---|
| 908 |
|
|---|
| 909 | call ophinv (dv1,dv2,dv3,besv1,besv2,besv3,h1,h2,tolhv,nmxv)
|
|---|
| 910 | call opadd2 (bx,by,bz,dv1,dv2,dv3)
|
|---|
| 911 |
|
|---|
| 912 | call incomprn (bx,by,bz,pm) ! project B onto div-free space
|
|---|
| 913 |
|
|---|
| 914 | return
|
|---|
| 915 | end
|
|---|
| 916 | c--------------------------------------------------------------------
|
|---|
| 917 | subroutine compute_cfl(cfl,u,v,w,dt)
|
|---|
| 918 | c
|
|---|
| 919 | c Given velocity field (u,v,w) and dt, compute current CFL number.
|
|---|
| 920 | c
|
|---|
| 921 | include 'SIZE'
|
|---|
| 922 | include 'GEOM'
|
|---|
| 923 | include 'INPUT'
|
|---|
| 924 | include 'WZ'
|
|---|
| 925 | include 'SOLN'
|
|---|
| 926 | c
|
|---|
| 927 | real u(lx1,ly1,lz1,nelv),v(lx1,ly1,lz1,nelv),w(lx1,ly1,lz1,nelv)
|
|---|
| 928 | c
|
|---|
| 929 | c Store the inverse jacobian to speed up this operation
|
|---|
| 930 | c
|
|---|
| 931 | common /cfldx/ dri(lx1),dsi(ly1),dti(lz1)
|
|---|
| 932 | c
|
|---|
| 933 | integer e
|
|---|
| 934 | c
|
|---|
| 935 | integer icalld
|
|---|
| 936 | save icalld
|
|---|
| 937 | data icalld /0/
|
|---|
| 938 | c
|
|---|
| 939 | if (icalld.eq.0) then
|
|---|
| 940 | icalld=1
|
|---|
| 941 | call getdr(dri,zgm1(1,1),lx1)
|
|---|
| 942 | call getdr(dsi,zgm1(1,2),ly1)
|
|---|
| 943 | if (if3d) call getdr(dti,zgm1(1,3),lz1)
|
|---|
| 944 | endif
|
|---|
| 945 |
|
|---|
| 946 | cfl = 0.
|
|---|
| 947 | l = 0
|
|---|
| 948 |
|
|---|
| 949 | if (if3d) then
|
|---|
| 950 | nxyz = lx1*ly1*lz1
|
|---|
| 951 | do e=1,nelv
|
|---|
| 952 | do k=1,lz1
|
|---|
| 953 | do j=1,ly1
|
|---|
| 954 | do i=1,lx1
|
|---|
| 955 | l = l+1
|
|---|
| 956 | ur = ( u(i,j,k,e)*rxm1(i,j,k,e)
|
|---|
| 957 | $ + v(i,j,k,e)*rym1(i,j,k,e)
|
|---|
| 958 | $ + w(i,j,k,e)*rzm1(i,j,k,e) ) * jacmi(l,1)
|
|---|
| 959 | us = ( u(i,j,k,e)*sxm1(i,j,k,e)
|
|---|
| 960 | $ + v(i,j,k,e)*sym1(i,j,k,e)
|
|---|
| 961 | $ + w(i,j,k,e)*szm1(i,j,k,e) ) * jacmi(l,1)
|
|---|
| 962 | ut = ( u(i,j,k,e)*txm1(i,j,k,e)
|
|---|
| 963 | $ + v(i,j,k,e)*tym1(i,j,k,e)
|
|---|
| 964 | $ + w(i,j,k,e)*tzm1(i,j,k,e) ) * jacmi(l,1)
|
|---|
| 965 |
|
|---|
| 966 | cflr = abs(dt*ur*dri(i))
|
|---|
| 967 | cfls = abs(dt*us*dsi(j))
|
|---|
| 968 | cflt = abs(dt*ut*dti(k))
|
|---|
| 969 |
|
|---|
| 970 | cflm = cflr + cfls + cflt
|
|---|
| 971 | cfl = max(cfl,cflm)
|
|---|
| 972 |
|
|---|
| 973 | cflf(i,j,k,e) = cflm
|
|---|
| 974 |
|
|---|
| 975 | enddo
|
|---|
| 976 | enddo
|
|---|
| 977 | enddo
|
|---|
| 978 | enddo
|
|---|
| 979 | else
|
|---|
| 980 | nxyz = lx1*ly1
|
|---|
| 981 | do e=1,nelv
|
|---|
| 982 | do j=1,ly1
|
|---|
| 983 | do i=1,lx1
|
|---|
| 984 | l = l+1
|
|---|
| 985 | ur = ( u(i,j,1,e)*rxm1(i,j,1,e)
|
|---|
| 986 | $ + v(i,j,1,e)*rym1(i,j,1,e) ) * jacmi(l,1)
|
|---|
| 987 | us = ( u(i,j,1,e)*sxm1(i,j,1,e)
|
|---|
| 988 | $ + v(i,j,1,e)*sym1(i,j,1,e) ) * jacmi(l,1)
|
|---|
| 989 |
|
|---|
| 990 | cflr = abs(dt*ur*dri(i))
|
|---|
| 991 | cfls = abs(dt*us*dsi(j))
|
|---|
| 992 |
|
|---|
| 993 | cflm = cflr + cfls
|
|---|
| 994 | cfl = max(cfl,cflm)
|
|---|
| 995 |
|
|---|
| 996 | cflf(i,j,1,e) = cflm
|
|---|
| 997 |
|
|---|
| 998 | enddo
|
|---|
| 999 | enddo
|
|---|
| 1000 | enddo
|
|---|
| 1001 | endif
|
|---|
| 1002 | c
|
|---|
| 1003 | cfl = glmax(cfl,1)
|
|---|
| 1004 | c
|
|---|
| 1005 | return
|
|---|
| 1006 | end
|
|---|
| 1007 | c-----------------------------------------------------------------------
|
|---|
| 1008 | subroutine getdr(dri,zgm1,lx1)
|
|---|
| 1009 | real dri(lx1),zgm1(lx1)
|
|---|
| 1010 | c
|
|---|
| 1011 | dri(1) = zgm1(2) - zgm1(1) ! Compute 1/Dx
|
|---|
| 1012 | do i=2,lx1-1
|
|---|
| 1013 | dri(i) = 0.5*( zgm1(i+1) - zgm1(i-1) )
|
|---|
| 1014 | enddo
|
|---|
| 1015 | dri(lx1) = zgm1(lx1) - zgm1(lx1-1)
|
|---|
| 1016 | c
|
|---|
| 1017 | call invcol1(dri,lx1)
|
|---|
| 1018 | c
|
|---|
| 1019 | return
|
|---|
| 1020 | end
|
|---|
| 1021 | c-----------------------------------------------------------------------
|
|---|
| 1022 | subroutine ophinv(o1,o2,o3,i1,i2,i3,h1,h2,tolh,nmxhi)
|
|---|
| 1023 | C
|
|---|
| 1024 | C Ok = (H1*A+H2*B)-1 * Ik (implicit)
|
|---|
| 1025 | C
|
|---|
| 1026 | include 'SIZE'
|
|---|
| 1027 | include 'INPUT'
|
|---|
| 1028 | include 'MASS'
|
|---|
| 1029 | include 'SOLN'
|
|---|
| 1030 | include 'VPROJ'
|
|---|
| 1031 | include 'TSTEP'
|
|---|
| 1032 | c logical ifproj
|
|---|
| 1033 |
|
|---|
| 1034 | real o1 (lx1,ly1,lz1,1) , o2 (lx1,ly1,lz1,1) , o3 (lx1,ly1,lz1,1)
|
|---|
| 1035 | real i1 (lx1,ly1,lz1,1) , i2 (lx1,ly1,lz1,1) , i3 (lx1,ly1,lz1,1)
|
|---|
| 1036 | real h1 (lx1,ly1,lz1,1) , h2 (lx1,ly1,lz1,1)
|
|---|
| 1037 |
|
|---|
| 1038 | c ifproj = .false.
|
|---|
| 1039 | c if (param(94).gt.0) ifproj = .true.
|
|---|
| 1040 | c if (ifprojfld(ifield)) ifproj = .true.
|
|---|
| 1041 | c
|
|---|
| 1042 | c if (.not.ifproj) then
|
|---|
| 1043 | c if (ifield.eq.1) call ophinv
|
|---|
| 1044 | c $ (o1,o2,o3,i1,i2,i3,h1,h2,tolh,nmxhi)
|
|---|
| 1045 | c if (ifield.eq.ifldmhd) call ophinv
|
|---|
| 1046 | c $ (o1,o2,o3,i1,i2,i3,h1,h2,tolh,nmxhi)
|
|---|
| 1047 | c return
|
|---|
| 1048 | c endif
|
|---|
| 1049 |
|
|---|
| 1050 | mtmp = param(93)
|
|---|
| 1051 | do i=1,2*ldim
|
|---|
| 1052 | ivproj(1,i) = min(mxprev,mtmp) - 1
|
|---|
| 1053 | enddo
|
|---|
| 1054 |
|
|---|
| 1055 | imesh = 1
|
|---|
| 1056 |
|
|---|
| 1057 | if (ifstrs) then
|
|---|
| 1058 | matmod = 0
|
|---|
| 1059 | call hmhzsf ('NOMG',o1,o2,o3,i1,i2,i3,h1,h2,
|
|---|
| 1060 | $ v1mask,v2mask,v3mask,vmult,
|
|---|
| 1061 | $ tolh,nmxhi,matmod)
|
|---|
| 1062 | else
|
|---|
| 1063 | if (ifield.eq.1) then
|
|---|
| 1064 | call hsolve ('VELX',o1,i1,h1,h2,v1mask,vmult
|
|---|
| 1065 | $ ,imesh,tolh,nmxhi,1
|
|---|
| 1066 | $ ,vproj(1,1),ivproj(1,1),binvm1)
|
|---|
| 1067 | call hsolve ('VELY',o2,i2,h1,h2,v2mask,vmult
|
|---|
| 1068 | $ ,imesh,tolh,nmxhi,2
|
|---|
| 1069 | $ ,vproj(1,2),ivproj(1,2),binvm1)
|
|---|
| 1070 | if (if3d)
|
|---|
| 1071 | $ call hsolve ('VELZ',o3,i3,h1,h2,v3mask,vmult
|
|---|
| 1072 | $ ,imesh,tolh,nmxhi,3
|
|---|
| 1073 | $ ,vproj(1,3),ivproj(1,3),binvm1)
|
|---|
| 1074 | elseif (ifield.eq.ifldmhd) then ! B-field
|
|---|
| 1075 | call hsolve (' BX ',o1,i1,h1,h2,b1mask,vmult
|
|---|
| 1076 | $ ,imesh,tolh,nmxhi,1
|
|---|
| 1077 | $ ,vproj(1,4),ivproj(1,4),binvm1)
|
|---|
| 1078 | call hsolve (' BY ',o2,i2,h1,h2,b2mask,vmult
|
|---|
| 1079 | $ ,imesh,tolh,nmxhi,2
|
|---|
| 1080 | $ ,vproj(1,5),ivproj(1,5),binvm1)
|
|---|
| 1081 | if (if3d)
|
|---|
| 1082 | $ call hsolve (' BZ ',o3,i3,h1,h2,b3mask,vmult
|
|---|
| 1083 | $ ,imesh,tolh,nmxhi,3
|
|---|
| 1084 | $ ,vproj(1,6),ivproj(1,6),binvm1)
|
|---|
| 1085 | endif
|
|---|
| 1086 | endif
|
|---|
| 1087 | C
|
|---|
| 1088 | return
|
|---|
| 1089 | end
|
|---|
| 1090 | c--------------------------------------------------------------------
|
|---|
| 1091 | subroutine set_ifbcor
|
|---|
| 1092 | include 'SIZE'
|
|---|
| 1093 | include 'GEOM'
|
|---|
| 1094 | include 'INPUT'
|
|---|
| 1095 | c include 'TSTEP' ! ifield?
|
|---|
| 1096 |
|
|---|
| 1097 | common /nekcb/ cb
|
|---|
| 1098 | character cb*3
|
|---|
| 1099 |
|
|---|
| 1100 | ifbcor = .true.
|
|---|
| 1101 | ifield = ifldmhd
|
|---|
| 1102 |
|
|---|
| 1103 | nface = 2*ldim
|
|---|
| 1104 | do iel=1,nelv
|
|---|
| 1105 | do ifc=1,nface
|
|---|
| 1106 | cb = cbc(ifc,iel,ifield)
|
|---|
| 1107 | if (cb.eq.'ndd' .or. cb.eq.'dnd' .or. cb.eq.'ddn')
|
|---|
| 1108 | $ ifbcor = .false.
|
|---|
| 1109 | enddo
|
|---|
| 1110 | enddo
|
|---|
| 1111 |
|
|---|
| 1112 | call gllog(ifbcor , .false.)
|
|---|
| 1113 |
|
|---|
| 1114 | if (nio.eq.0) write (6,*) 'IFBCOR =',ifbcor
|
|---|
| 1115 |
|
|---|
| 1116 | return
|
|---|
| 1117 | end
|
|---|
| 1118 | c--------------------------------------------------------------------
|
|---|
| 1119 | c subroutine set_ifbcor_old(ifbcor)
|
|---|
| 1120 | cc
|
|---|
| 1121 | cc This is a quick hack for the rings problem - it is not general,
|
|---|
| 1122 | cc but will also work fine for the periodic in Z problem
|
|---|
| 1123 | cc
|
|---|
| 1124 | c include 'SIZE'
|
|---|
| 1125 | c include 'TOTAL'
|
|---|
| 1126 | c
|
|---|
| 1127 | c logical ifbcor
|
|---|
| 1128 | c
|
|---|
| 1129 | c integer e,f
|
|---|
| 1130 | c character*1 cb1(3)
|
|---|
| 1131 | c
|
|---|
| 1132 | c itest = 0
|
|---|
| 1133 | c
|
|---|
| 1134 | c do e=1,nelv
|
|---|
| 1135 | c do f=1,2*ldim
|
|---|
| 1136 | c call chcopy(cb1,cbc(f,e,ifldmhd),3)
|
|---|
| 1137 | c if (cb1(3).eq.'n'.or.cb1(3).eq.'N') itest = 1
|
|---|
| 1138 | c enddo
|
|---|
| 1139 | c enddo
|
|---|
| 1140 | c
|
|---|
| 1141 | c itest = iglmax(itest,1)
|
|---|
| 1142 | c
|
|---|
| 1143 | c ifbcor = .true. ! adjust mean pressure to rmv hydrostatic mode
|
|---|
| 1144 | c
|
|---|
| 1145 | c
|
|---|
| 1146 | c if (itest.eq.1) ifbcor = .false. ! do not adjust mean pressure
|
|---|
| 1147 | c
|
|---|
| 1148 | c return
|
|---|
| 1149 | c end
|
|---|
| 1150 | c--------------------------------------------------------------------
|
|---|
| 1151 | subroutine setrhsp(p,h1,h2,h2inv,pset,niprev)
|
|---|
| 1152 | C
|
|---|
| 1153 | C Project soln onto best fit in the "E" norm.
|
|---|
| 1154 | C
|
|---|
| 1155 | include 'SIZE'
|
|---|
| 1156 | include 'INPUT'
|
|---|
| 1157 | include 'MASS'
|
|---|
| 1158 | include 'SOLN'
|
|---|
| 1159 | include 'TSTEP'
|
|---|
| 1160 |
|
|---|
| 1161 | real p (lx2,ly2,lz2,lelv)
|
|---|
| 1162 | real h1 (lx1,ly1,lz1,lelv)
|
|---|
| 1163 | real h2 (lx1,ly1,lz1,lelv)
|
|---|
| 1164 | real h2inv(lx1,ly1,lz1,lelv)
|
|---|
| 1165 | real pset (lx2*ly2*lz2*lelv,mxprev)
|
|---|
| 1166 |
|
|---|
| 1167 | parameter (ltot2=lx2*ly2*lz2*lelv)
|
|---|
| 1168 | common /orthox/ pbar(ltot2),pnew(ltot2)
|
|---|
| 1169 | common /orthos/ alpha(mxprev),work(mxprev)
|
|---|
| 1170 | common /orthoi/ nprev,mprev
|
|---|
| 1171 |
|
|---|
| 1172 | nprev = niprev
|
|---|
| 1173 | if (nprev.eq.0) return
|
|---|
| 1174 |
|
|---|
| 1175 | c Diag to see how much reduction in the residual is attained.
|
|---|
| 1176 | ntot2 = lx2*ly2*lz2*nelv
|
|---|
| 1177 | alpha1 = glsc3(p,p,bm2inv,ntot2)
|
|---|
| 1178 | if (alpha1.gt.0) then
|
|---|
| 1179 | alpha1 = sqrt(alpha1/volvm2)
|
|---|
| 1180 | else
|
|---|
| 1181 | return
|
|---|
| 1182 | endif
|
|---|
| 1183 |
|
|---|
| 1184 | call updrhse(p,h1,h2,h2inv,ierr) ! Update rhs's if E-matrix has changed
|
|---|
| 1185 | c if (ierr.eq.1) Nprev=0 ! Doesn't happen w/ new formulation
|
|---|
| 1186 |
|
|---|
| 1187 | do i=1,nprev ! Perform Gram-Schmidt for previous soln's.
|
|---|
| 1188 | alpha(i) = vlsc2(p,pset(1,i),ntot2)
|
|---|
| 1189 | enddo
|
|---|
| 1190 | call gop(alpha,work,'+ ',nprev)
|
|---|
| 1191 |
|
|---|
| 1192 | call rzero(pbar,ntot2)
|
|---|
| 1193 | do i=1,nprev
|
|---|
| 1194 | call add2s2(pbar,pset(1,i),alpha(i),ntot2)
|
|---|
| 1195 | enddo
|
|---|
| 1196 | C
|
|---|
| 1197 | intetype = 1
|
|---|
| 1198 | call cdabdtp(pnew,pbar,h1,h2,h2inv,intetype)
|
|---|
| 1199 | call sub2 (p,pnew,ntot2)
|
|---|
| 1200 |
|
|---|
| 1201 | c ................................................................
|
|---|
| 1202 | alpha2 = glsc3(p,p,bm2inv,ntot2) ! Diagnostics
|
|---|
| 1203 | if (alpha2.gt.0) then
|
|---|
| 1204 | alpha2 = sqrt(alpha2/volvm2)
|
|---|
| 1205 | ratio = alpha1/alpha2
|
|---|
| 1206 | n10=min(10,nprev)
|
|---|
| 1207 | c if (nio.eq.0) write(6,11) istep,nprev,(alpha(i),i=1,n10)
|
|---|
| 1208 | c if (nio.eq.0) write(6,12) istep,nprev,alpha1,alpha2,ratio
|
|---|
| 1209 | 11 format(2i5,' alpha:',1p10e12.4)
|
|---|
| 1210 | 12 format(i6,i4,1p3e12.4,' alph12')
|
|---|
| 1211 |
|
|---|
| 1212 | if (nio.eq.0) write(6,13) istep,' Project PRES ',
|
|---|
| 1213 | & alpha2,alpha1,ratio,nprev,mxprev
|
|---|
| 1214 | 13 format(i11,a,6x,1p3e13.4,i4,i4)
|
|---|
| 1215 | endif
|
|---|
| 1216 | c ................................................................
|
|---|
| 1217 |
|
|---|
| 1218 | return
|
|---|
| 1219 | end
|
|---|
| 1220 | c-----------------------------------------------------------------------
|
|---|
| 1221 | subroutine gensolnp(p,h1,h2,h2inv,pset,nprev)
|
|---|
| 1222 | C
|
|---|
| 1223 | C Reconstruct the solution to the original problem by adding back
|
|---|
| 1224 | C the previous solutions
|
|---|
| 1225 | C
|
|---|
| 1226 | include 'SIZE'
|
|---|
| 1227 | include 'INPUT'
|
|---|
| 1228 |
|
|---|
| 1229 | real p (lx2,ly2,lz2,lelv)
|
|---|
| 1230 | real h1 (lx1,ly1,lz1,lelv)
|
|---|
| 1231 | real h2 (lx1,ly1,lz1,lelv)
|
|---|
| 1232 | real h2inv(lx1,ly1,lz1,lelv)
|
|---|
| 1233 | real pset (lx2*ly2*lz2*lelv,mxprev)
|
|---|
| 1234 |
|
|---|
| 1235 | parameter (ltot2=lx2*ly2*lz2*lelv)
|
|---|
| 1236 | common /orthox/ pbar(ltot2),pnew(ltot2)
|
|---|
| 1237 | common /orthos/ alpha(mxprev),work(mxprev)
|
|---|
| 1238 |
|
|---|
| 1239 | mprev=param(93)
|
|---|
| 1240 | mprev=min(mprev,mxprev)
|
|---|
| 1241 |
|
|---|
| 1242 | ntot2=lx2*ly2*lz2*nelv
|
|---|
| 1243 |
|
|---|
| 1244 | if (nprev.lt.mprev) then
|
|---|
| 1245 | nprev = nprev+1
|
|---|
| 1246 | call copy (pset(1,nprev),p,ntot2) ! Save current solution
|
|---|
| 1247 | call add2 (p,pbar,ntot2) ! Reconstruct solution.
|
|---|
| 1248 | call econjp(pset,nprev,h1,h2,h2inv,ierr) ! Orthonormalize set
|
|---|
| 1249 |
|
|---|
| 1250 | if (ierr.eq.1) then
|
|---|
| 1251 | nprev = 1
|
|---|
| 1252 | call copy (pset(1,nprev),p,ntot2) ! Save current solution
|
|---|
| 1253 | call econjp(pset,nprev,h1,h2,h2inv,ierr) ! and orthonormalize.
|
|---|
| 1254 | endif
|
|---|
| 1255 |
|
|---|
| 1256 | else ! (uses pnew).
|
|---|
| 1257 | nprev = 1
|
|---|
| 1258 | call add2 (p,pbar,ntot2) ! Reconstruct solution.
|
|---|
| 1259 | call copy (pset(1,nprev),p,ntot2) ! Save current solution
|
|---|
| 1260 | call econjp(pset,nprev,h1,h2,h2inv,ierr) ! and orthonormalize.
|
|---|
| 1261 | endif
|
|---|
| 1262 |
|
|---|
| 1263 | return
|
|---|
| 1264 | end
|
|---|
| 1265 | c-----------------------------------------------------------------------
|
|---|
| 1266 | subroutine econjp(pset,nprev,h1,h2,h2inv,ierr)
|
|---|
| 1267 |
|
|---|
| 1268 | c Orthogonalize the soln wrt previous soln's for which we already
|
|---|
| 1269 | c know the soln.
|
|---|
| 1270 |
|
|---|
| 1271 | include 'SIZE'
|
|---|
| 1272 | include 'INPUT'
|
|---|
| 1273 |
|
|---|
| 1274 | real p (lx2,ly2,lz2,lelv)
|
|---|
| 1275 | real h1 (lx1,ly1,lz1,lelv)
|
|---|
| 1276 | real h2 (lx1,ly1,lz1,lelv)
|
|---|
| 1277 | real h2inv(lx1,ly1,lz1,lelv)
|
|---|
| 1278 | real pset (lx2*ly2*lz2*lelv,mxprev)
|
|---|
| 1279 |
|
|---|
| 1280 | parameter (ltot2=lx2*ly2*lz2*lelv)
|
|---|
| 1281 | common /orthox/ pbar(ltot2),pnew(ltot2)
|
|---|
| 1282 | common /orthos/ alpha(mxprev),work(mxprev)
|
|---|
| 1283 |
|
|---|
| 1284 | ierr = 0
|
|---|
| 1285 |
|
|---|
| 1286 | ntot2=lx2*ly2*lz2*nelv
|
|---|
| 1287 |
|
|---|
| 1288 | C
|
|---|
| 1289 | C Gram Schmidt, w re-orthogonalization
|
|---|
| 1290 | C
|
|---|
| 1291 | npass=1
|
|---|
| 1292 | c if (abs(param(105)).eq.2) npass=2
|
|---|
| 1293 | do ipass=1,npass
|
|---|
| 1294 |
|
|---|
| 1295 | intetype=1
|
|---|
| 1296 | call cdabdtp(pnew,pset(1,nprev),h1,h2,h2inv,intetype)
|
|---|
| 1297 | alphad = glsc2(pnew,pset(1,nprev),ntot2) ! compute part of the norm
|
|---|
| 1298 |
|
|---|
| 1299 | nprev1 = nprev-1
|
|---|
| 1300 | do i=1,nprev1 ! Gram-Schmidt
|
|---|
| 1301 | alpha(i) = vlsc2(pnew,pset(1,i),ntot2)
|
|---|
| 1302 | enddo
|
|---|
| 1303 | if (nprev1.gt.0) call gop(alpha,work,'+ ',nprev1)
|
|---|
| 1304 |
|
|---|
| 1305 | do i=1,nprev1
|
|---|
| 1306 | alpham = -alpha(i)
|
|---|
| 1307 | call add2s2(pset(1,nprev),pset(1,i),alpham,ntot2)
|
|---|
| 1308 | alphad = alphad - alpha(i)**2
|
|---|
| 1309 | enddo
|
|---|
| 1310 | enddo
|
|---|
| 1311 | C
|
|---|
| 1312 | C .Normalize new element in P~
|
|---|
| 1313 | C
|
|---|
| 1314 | if (alphad.le.0) then
|
|---|
| 1315 | write(6,*) 'ERROR: alphad .le. 0 in econjp',alphad,nprev
|
|---|
| 1316 | ierr = 1
|
|---|
| 1317 | return
|
|---|
| 1318 | endif
|
|---|
| 1319 | alphad = 1./sqrt(alphad)
|
|---|
| 1320 | call cmult(pset(1,nprev),alphad,ntot2)
|
|---|
| 1321 |
|
|---|
| 1322 | return
|
|---|
| 1323 | end
|
|---|
| 1324 | c-----------------------------------------------------------------------
|
|---|
| 1325 | subroutine advab_elsasser_fast
|
|---|
| 1326 | C
|
|---|
| 1327 | C Eulerian scheme, add convection term to forcing function
|
|---|
| 1328 | C at current time step.
|
|---|
| 1329 | C
|
|---|
| 1330 | include 'SIZE'
|
|---|
| 1331 | include 'INPUT'
|
|---|
| 1332 | include 'GEOM'
|
|---|
| 1333 | include 'SOLN'
|
|---|
| 1334 | include 'MASS'
|
|---|
| 1335 | include 'TSTEP'
|
|---|
| 1336 | c
|
|---|
| 1337 | parameter (lxy=lx1*ly1*lz1,ltd=lxd*lyd*lzd)
|
|---|
| 1338 | common /scrns/ wk(2*ltd)
|
|---|
| 1339 | $ , fx(lxy),fy(lxy),fz(lxy)
|
|---|
| 1340 | $ , gx(lxy),gy(lxy),gz(lxy)
|
|---|
| 1341 | $ , zr(ltd),zs(ltd),zt(ltd)
|
|---|
| 1342 | $ , tr(ltd,3),zp(ltd,3),zm(ltd,3)
|
|---|
| 1343 |
|
|---|
| 1344 | integer e
|
|---|
| 1345 | integer icalld
|
|---|
| 1346 | save icalld
|
|---|
| 1347 | data icalld /0/
|
|---|
| 1348 |
|
|---|
| 1349 | if (icalld.eq.0) call set_dealias_rx
|
|---|
| 1350 | icalld = icalld + 1
|
|---|
| 1351 |
|
|---|
| 1352 | nxyz1 = lx1*ly1*lz1
|
|---|
| 1353 | nxyzd = lxd*lyd*lzd
|
|---|
| 1354 |
|
|---|
| 1355 |
|
|---|
| 1356 | if (if3d) then
|
|---|
| 1357 | c
|
|---|
| 1358 | do e=1,nelv
|
|---|
| 1359 |
|
|---|
| 1360 | c Interpolate z+ and z- into fine mesh, translate to r-s-t coords
|
|---|
| 1361 |
|
|---|
| 1362 | c write(6,*) nid,' inside fast',e,nxyz1
|
|---|
| 1363 |
|
|---|
| 1364 | call add3(wk,vx(1,1,1,e),bx(1,1,1,e),nxyz1)
|
|---|
| 1365 | call intp_rstd(zp(1,1),wk,lx1,lxd,if3d,0) ! 0 --> forward
|
|---|
| 1366 |
|
|---|
| 1367 | call add3(wk,vy(1,1,1,e),by(1,1,1,e),nxyz1)
|
|---|
| 1368 | call intp_rstd(zp(1,2),wk,lx1,lxd,if3d,0)
|
|---|
| 1369 |
|
|---|
| 1370 | call add3(wk,vz(1,1,1,e),bz(1,1,1,e),nxyz1)
|
|---|
| 1371 | call intp_rstd(zp(1,3),wk,lx1,lxd,if3d,0)
|
|---|
| 1372 |
|
|---|
| 1373 | call sub3(wk,vx(1,1,1,e),bx(1,1,1,e),nxyz1)
|
|---|
| 1374 | call intp_rstd(zm(1,1),wk,lx1,lxd,if3d,0)
|
|---|
| 1375 |
|
|---|
| 1376 | call sub3(wk,vy(1,1,1,e),by(1,1,1,e),nxyz1)
|
|---|
| 1377 | call intp_rstd(zm(1,2),wk,lx1,lxd,if3d,0)
|
|---|
| 1378 |
|
|---|
| 1379 | call sub3(wk,vz(1,1,1,e),bz(1,1,1,e),nxyz1)
|
|---|
| 1380 | call intp_rstd(zm(1,3),wk,lx1,lxd,if3d,0)
|
|---|
| 1381 |
|
|---|
| 1382 | do i=1,nxyzd ! Convert convector (zm) to r-s-t coordinates
|
|---|
| 1383 | tr(i,1)=
|
|---|
| 1384 | $ rx(i,1,e)*zm(i,1)+rx(i,2,e)*zm(i,2)+rx(i,3,e)*zm(i,3)
|
|---|
| 1385 | tr(i,2)=
|
|---|
| 1386 | $ rx(i,4,e)*zm(i,1)+rx(i,5,e)*zm(i,2)+rx(i,6,e)*zm(i,3)
|
|---|
| 1387 | tr(i,3)=
|
|---|
| 1388 | $ rx(i,7,e)*zm(i,1)+rx(i,8,e)*zm(i,2)+rx(i,9,e)*zm(i,3)
|
|---|
| 1389 | enddo
|
|---|
| 1390 |
|
|---|
| 1391 |
|
|---|
| 1392 | call grad_rst(zr,zs,zt,zp(1,1),lxd,if3d)
|
|---|
| 1393 | do i=1,nxyzd ! mass matrix included, per DFM (4.8.5)
|
|---|
| 1394 | wk(i) = tr(i,1)*zr(i)+tr(i,2)*zs(i)+tr(i,3)*zt(i)
|
|---|
| 1395 | enddo
|
|---|
| 1396 | call intp_rstd(fx,wk,lx1,lxd,if3d,1) ! Project back to coarse
|
|---|
| 1397 |
|
|---|
| 1398 | call grad_rst(zr,zs,zt,zp(1,2),lxd,if3d)
|
|---|
| 1399 | do i=1,nxyzd ! mass matrix included, per DFM (4.8.5)
|
|---|
| 1400 | wk(i) = tr(i,1)*zr(i)+tr(i,2)*zs(i)+tr(i,3)*zt(i)
|
|---|
| 1401 | enddo
|
|---|
| 1402 | call intp_rstd(fy,wk,lx1,lxd,if3d,1) ! Project back to coarse
|
|---|
| 1403 |
|
|---|
| 1404 | call grad_rst(zr,zs,zt,zp(1,3),lxd,if3d)
|
|---|
| 1405 | do i=1,nxyzd ! mass matrix included, per DFM (4.8.5)
|
|---|
| 1406 | wk(i) = tr(i,1)*zr(i)+tr(i,2)*zs(i)+tr(i,3)*zt(i)
|
|---|
| 1407 | enddo
|
|---|
| 1408 | call intp_rstd(fz,wk,lx1,lxd,if3d,1) ! Project back to coarse
|
|---|
| 1409 |
|
|---|
| 1410 |
|
|---|
| 1411 | do i=1,nxyzd ! Convert convector (zp) to r-s-t coordinates
|
|---|
| 1412 | tr(i,1)=
|
|---|
| 1413 | $ rx(i,1,e)*zp(i,1)+rx(i,2,e)*zp(i,2)+rx(i,3,e)*zp(i,3)
|
|---|
| 1414 | tr(i,2)=
|
|---|
| 1415 | $ rx(i,4,e)*zp(i,1)+rx(i,5,e)*zp(i,2)+rx(i,6,e)*zp(i,3)
|
|---|
| 1416 | tr(i,3)=
|
|---|
| 1417 | $ rx(i,7,e)*zp(i,1)+rx(i,8,e)*zp(i,2)+rx(i,9,e)*zp(i,3)
|
|---|
| 1418 | enddo
|
|---|
| 1419 |
|
|---|
| 1420 |
|
|---|
| 1421 | call grad_rst(zr,zs,zt,zm(1,1),lxd,if3d)
|
|---|
| 1422 | do i=1,nxyzd ! mass matrix included, per DFM (4.8.5)
|
|---|
| 1423 | wk(i) = tr(i,1)*zr(i)+tr(i,2)*zs(i)+tr(i,3)*zt(i)
|
|---|
| 1424 | enddo
|
|---|
| 1425 | call intp_rstd(gx,wk,lx1,lxd,if3d,1) ! Project back to coarse
|
|---|
| 1426 |
|
|---|
| 1427 | call grad_rst(zr,zs,zt,zm(1,2),lxd,if3d)
|
|---|
| 1428 | do i=1,nxyzd ! mass matrix included, per DFM (4.8.5)
|
|---|
| 1429 | wk(i) = tr(i,1)*zr(i)+tr(i,2)*zs(i)+tr(i,3)*zt(i)
|
|---|
| 1430 | enddo
|
|---|
| 1431 | call intp_rstd(gy,wk,lx1,lxd,if3d,1) ! Project back to coarse
|
|---|
| 1432 |
|
|---|
| 1433 | call grad_rst(zr,zs,zt,zm(1,3),lxd,if3d)
|
|---|
| 1434 | do i=1,nxyzd ! mass matrix included, per DFM (4.8.5)
|
|---|
| 1435 | wk(i) = tr(i,1)*zr(i)+tr(i,2)*zs(i)+tr(i,3)*zt(i)
|
|---|
| 1436 | enddo
|
|---|
| 1437 | call intp_rstd(gz,wk,lx1,lxd,if3d,1) ! Project back to coarse
|
|---|
| 1438 |
|
|---|
| 1439 | tmp = -0.5 ! vtrans() assumed to be 1.0 !
|
|---|
| 1440 | do i=1,nxyz1
|
|---|
| 1441 | bfx(i,1,1,e) = bfx(i,1,1,e) + tmp*( fx(i) + gx(i) )
|
|---|
| 1442 | bmx(i,1,1,e) = bmx(i,1,1,e) + tmp*( fx(i) - gx(i) )
|
|---|
| 1443 |
|
|---|
| 1444 | bfy(i,1,1,e) = bfy(i,1,1,e) + tmp*( fy(i) + gy(i) )
|
|---|
| 1445 | bmy(i,1,1,e) = bmy(i,1,1,e) + tmp*( fy(i) - gy(i) )
|
|---|
| 1446 |
|
|---|
| 1447 | bfz(i,1,1,e) = bfz(i,1,1,e) + tmp*( fz(i) + gz(i) )
|
|---|
| 1448 | bmz(i,1,1,e) = bmz(i,1,1,e) + tmp*( fz(i) - gz(i) )
|
|---|
| 1449 | enddo
|
|---|
| 1450 |
|
|---|
| 1451 | enddo
|
|---|
| 1452 |
|
|---|
| 1453 | else ! 2D
|
|---|
| 1454 | c
|
|---|
| 1455 | do e=1,nelv
|
|---|
| 1456 |
|
|---|
| 1457 | c Interpolate z+ and z- into fine mesh, translate to r-s-t coords
|
|---|
| 1458 |
|
|---|
| 1459 | call add3(wk,vx(1,1,1,e),bx(1,1,1,e),nxyz1)
|
|---|
| 1460 | call intp_rstd(zp(1,1),wk,lx1,lxd,if3d,0) ! 0 --> forward
|
|---|
| 1461 |
|
|---|
| 1462 | call add3(wk,vy(1,1,1,e),by(1,1,1,e),nxyz1)
|
|---|
| 1463 | call intp_rstd(zp(1,2),wk,lx1,lxd,if3d,0)
|
|---|
| 1464 |
|
|---|
| 1465 | call sub3(wk,vx(1,1,1,e),bx(1,1,1,e),nxyz1)
|
|---|
| 1466 | call intp_rstd(zm(1,1),wk,lx1,lxd,if3d,0)
|
|---|
| 1467 |
|
|---|
| 1468 | call sub3(wk,vy(1,1,1,e),by(1,1,1,e),nxyz1)
|
|---|
| 1469 | call intp_rstd(zm(1,2),wk,lx1,lxd,if3d,0)
|
|---|
| 1470 |
|
|---|
| 1471 | do i=1,nxyzd ! Convert convector (zm) to r-s-t coordinates
|
|---|
| 1472 | tr(i,1)=
|
|---|
| 1473 | $ rx(i,1,e)*zm(i,1)+rx(i,2,e)*zm(i,2)
|
|---|
| 1474 | tr(i,2)=
|
|---|
| 1475 | $ rx(i,3,e)*zm(i,1)+rx(i,4,e)*zm(i,2)
|
|---|
| 1476 | enddo
|
|---|
| 1477 |
|
|---|
| 1478 | call grad_rst(zr,zs,zt,zp(1,1),lxd,if3d)
|
|---|
| 1479 | do i=1,nxyzd ! mass matrix included, per DFM (4.8.5)
|
|---|
| 1480 | wk(i) = tr(i,1)*zr(i)+tr(i,2)*zs(i)
|
|---|
| 1481 | enddo
|
|---|
| 1482 | call intp_rstd(fx,wk,lx1,lxd,if3d,1) ! Project back to coarse
|
|---|
| 1483 |
|
|---|
| 1484 | call grad_rst(zr,zs,zt,zp(1,2),lxd,if3d)
|
|---|
| 1485 | do i=1,nxyzd ! mass matrix included, per DFM (4.8.5)
|
|---|
| 1486 | wk(i) = tr(i,1)*zr(i)+tr(i,2)*zs(i)
|
|---|
| 1487 | enddo
|
|---|
| 1488 | call intp_rstd(fy,wk,lx1,lxd,if3d,1) ! Project back to coarse
|
|---|
| 1489 |
|
|---|
| 1490 | do i=1,nxyzd ! Convert convector (zp) to r-s-t coordinates
|
|---|
| 1491 | tr(i,1)=
|
|---|
| 1492 | $ rx(i,1,e)*zp(i,1)+rx(i,2,e)*zp(i,2)
|
|---|
| 1493 | tr(i,2)=
|
|---|
| 1494 | $ rx(i,3,e)*zp(i,1)+rx(i,4,e)*zp(i,2)
|
|---|
| 1495 | enddo
|
|---|
| 1496 |
|
|---|
| 1497 | call grad_rst(zr,zs,zt,zm(1,1),lxd,if3d)
|
|---|
| 1498 | do i=1,nxyzd ! mass matrix included, per DFM (4.8.5)
|
|---|
| 1499 | wk(i) = tr(i,1)*zr(i)+tr(i,2)*zs(i)
|
|---|
| 1500 | enddo
|
|---|
| 1501 | call intp_rstd(gx,wk,lx1,lxd,if3d,1) ! Project back to coarse
|
|---|
| 1502 |
|
|---|
| 1503 | call grad_rst(zr,zs,zt,zm(1,2),lxd,if3d)
|
|---|
| 1504 | do i=1,nxyzd ! mass matrix included, per DFM (4.8.5)
|
|---|
| 1505 | wk(i) = tr(i,1)*zr(i)+tr(i,2)*zs(i)
|
|---|
| 1506 | enddo
|
|---|
| 1507 | call intp_rstd(gy,wk,lx1,lxd,if3d,1) ! Project back to coarse
|
|---|
| 1508 |
|
|---|
| 1509 | tmp = -0.5 ! vtrans() assumed to be 1.0 !
|
|---|
| 1510 | do i=1,nxyz1
|
|---|
| 1511 | bfx(i,1,1,e) = bfx(i,1,1,e) + tmp*( fx(i) + gx(i) )
|
|---|
| 1512 | bmx(i,1,1,e) = bmx(i,1,1,e) + tmp*( fx(i) - gx(i) )
|
|---|
| 1513 |
|
|---|
| 1514 | bfy(i,1,1,e) = bfy(i,1,1,e) + tmp*( fy(i) + gy(i) )
|
|---|
| 1515 | bmy(i,1,1,e) = bmy(i,1,1,e) + tmp*( fy(i) - gy(i) )
|
|---|
| 1516 | enddo
|
|---|
| 1517 | enddo
|
|---|
| 1518 | endif
|
|---|
| 1519 |
|
|---|
| 1520 | return
|
|---|
| 1521 | end
|
|---|
| 1522 | c-----------------------------------------------------------------------
|
|---|
| 1523 | subroutine cfl_check
|
|---|
| 1524 | c
|
|---|
| 1525 | include 'SIZE'
|
|---|
| 1526 | include 'INPUT'
|
|---|
| 1527 | include 'SOLN'
|
|---|
| 1528 | include 'MASS'
|
|---|
| 1529 | include 'TSTEP'
|
|---|
| 1530 | c
|
|---|
| 1531 | common /scrns/ ta1 (lx1,ly1,lz1,lelv)
|
|---|
| 1532 | $ , ta2 (lx1,ly1,lz1,lelv)
|
|---|
| 1533 | $ , ta3 (lx1,ly1,lz1,lelv)
|
|---|
| 1534 | $ , tb1 (lx1,ly1,lz1,lelv)
|
|---|
| 1535 | $ , tb2 (lx1,ly1,lz1,lelv)
|
|---|
| 1536 | $ , tb3 (lx1,ly1,lz1,lelv)
|
|---|
| 1537 |
|
|---|
| 1538 |
|
|---|
| 1539 | call opcopy(ta1,ta2,ta3,vx,vy,vz)
|
|---|
| 1540 | call opcopy(tb1,tb2,tb3,bx,by,bz)
|
|---|
| 1541 |
|
|---|
| 1542 | ntot1 = lx1*ly1*lz1*nelv
|
|---|
| 1543 | call phys_to_elsasser(ta1,ta2,ta3,tb1,tb2,tb3,ntot1) ! crude, but effective
|
|---|
| 1544 |
|
|---|
| 1545 | call compute_cfl(cflp,ta1,ta2,ta3,dt) ! vx = U+B
|
|---|
| 1546 | call compute_cfl(cflm,tb1,tb2,tb3,dt) ! bx = U-B
|
|---|
| 1547 |
|
|---|
| 1548 | courno = max(cflp,cflm)
|
|---|
| 1549 |
|
|---|
| 1550 | if (nio.eq.0) write(6,1) istep,time,dt,cflp,cflm
|
|---|
| 1551 | 1 format(i9,1p4e15.7,' CFL')
|
|---|
| 1552 |
|
|---|
| 1553 | return
|
|---|
| 1554 | end
|
|---|
| 1555 | c--------------------------------------------------------------------
|
|---|
