| 1 | c-----------------------------------------------------------------------
|
|---|
| 2 | subroutine cdscal (igeom)
|
|---|
| 3 | C
|
|---|
| 4 | C Solve the convection-diffusion equation for passive scalar IPSCAL
|
|---|
| 5 | C
|
|---|
| 6 | include 'SIZE'
|
|---|
| 7 | include 'INPUT'
|
|---|
| 8 | include 'GEOM'
|
|---|
| 9 | include 'MVGEOM'
|
|---|
| 10 | include 'SOLN'
|
|---|
| 11 | include 'MASS'
|
|---|
| 12 | include 'TSTEP'
|
|---|
| 13 | COMMON /CPRINT/ IFPRINT
|
|---|
| 14 | LOGICAL IFPRINT
|
|---|
| 15 | LOGICAL IFCONV
|
|---|
| 16 |
|
|---|
| 17 | COMMON /SCRNS/ TA(LX1,LY1,LZ1,LELT)
|
|---|
| 18 | $ ,TB(LX1,LY1,LZ1,LELT)
|
|---|
| 19 | COMMON /SCRVH/ H1(LX1,LY1,LZ1,LELT)
|
|---|
| 20 | $ ,H2(LX1,LY1,LZ1,LELT)
|
|---|
| 21 |
|
|---|
| 22 | include 'ORTHOT'
|
|---|
| 23 |
|
|---|
| 24 | if (ifdgfld(ifield)) then
|
|---|
| 25 | call cdscal_dg(igeom)
|
|---|
| 26 | return
|
|---|
| 27 | endif
|
|---|
| 28 |
|
|---|
| 29 |
|
|---|
| 30 | ifld1 = ifield-1
|
|---|
| 31 | napproxt(1,ifld1) = laxtt
|
|---|
| 32 |
|
|---|
| 33 | nel = nelfld(ifield)
|
|---|
| 34 | n = lx1*ly1*lz1*nel
|
|---|
| 35 |
|
|---|
| 36 | if (igeom.eq.1) then ! geometry at t^{n-1}
|
|---|
| 37 |
|
|---|
| 38 | call makeq
|
|---|
| 39 | call lagscal
|
|---|
| 40 |
|
|---|
| 41 | else ! geometry at t^n
|
|---|
| 42 |
|
|---|
| 43 | IF (IFPRINT) THEN
|
|---|
| 44 | IF (IFIELD.EQ.2.AND.NID.EQ.0)
|
|---|
| 45 | $ WRITE (6,*) ' Temperature/Passive scalar solution'
|
|---|
| 46 | ENDIF
|
|---|
| 47 |
|
|---|
| 48 | if1=ifield-1
|
|---|
| 49 | write(name4t,1) if1-1
|
|---|
| 50 | 1 format('PS',i2)
|
|---|
| 51 | if(ifield.eq.2) write(name4t,'(A4)') 'TEMP'
|
|---|
| 52 |
|
|---|
| 53 | C
|
|---|
| 54 | C New geometry
|
|---|
| 55 | C
|
|---|
| 56 | isd = 1
|
|---|
| 57 | if (ifaxis.and.ifaziv.and.ifield.eq.2) isd = 2
|
|---|
| 58 | c if (ifaxis.and.ifmhd) isd = 2 !This is a problem if T is to be T!
|
|---|
| 59 |
|
|---|
| 60 | do 1000 iter=1,nmxnl ! iterate for nonlin. prob. (e.g. radiation b.c.)
|
|---|
| 61 |
|
|---|
| 62 | intype = 0
|
|---|
| 63 | if (iftran) intype = -1
|
|---|
| 64 | call sethlm (h1,h2,intype)
|
|---|
| 65 | call bcneusc (ta,-1)
|
|---|
| 66 | call add2 (h2,ta,n)
|
|---|
| 67 | call bcdirsc (t(1,1,1,1,ifield-1))
|
|---|
| 68 | call axhelm (ta,t(1,1,1,1,ifield-1),h1,h2,imesh,ISD)
|
|---|
| 69 | call sub3 (tb,bq(1,1,1,1,ifield-1),ta,n)
|
|---|
| 70 | call bcneusc (ta,1)
|
|---|
| 71 | call add2 (tb,ta,n)
|
|---|
| 72 |
|
|---|
| 73 | if(iftmsh(ifield)) then
|
|---|
| 74 | call hsolve (name4t,TA,TB,H1,H2
|
|---|
| 75 | $ ,tmask(1,1,1,1,ifield-1)
|
|---|
| 76 | $ ,tmult(1,1,1,1,ifield-1)
|
|---|
| 77 | $ ,imesh,tolht(ifield),nmxt(ifield-1),1
|
|---|
| 78 | $ ,approxt(1,0,ifld1),napproxt(1,ifld1),bintm1)
|
|---|
| 79 | else
|
|---|
| 80 | call hsolve (name4t,TA,TB,H1,H2
|
|---|
| 81 | $ ,tmask(1,1,1,1,ifield-1)
|
|---|
| 82 | $ ,tmult(1,1,1,1,ifield-1)
|
|---|
| 83 | $ ,imesh,tolht(ifield),nmxt(ifield-1),1
|
|---|
| 84 | $ ,approxt(1,0,ifld1),napproxt(1,ifld1),binvm1)
|
|---|
| 85 | endif
|
|---|
| 86 |
|
|---|
| 87 | call add2 (t(1,1,1,1,ifield-1),ta,n)
|
|---|
| 88 |
|
|---|
| 89 | call cvgnlps (ifconv) ! Check convergence for nonlinear problem
|
|---|
| 90 | if (ifconv) goto 2000
|
|---|
| 91 |
|
|---|
| 92 | C Radiation case, smooth convergence, avoid flip-flop (ER).
|
|---|
| 93 | call cmult (ta,0.5,n)
|
|---|
| 94 | call sub2 (t(1,1,1,1,ifield-1),ta,n)
|
|---|
| 95 |
|
|---|
| 96 | 1000 continue
|
|---|
| 97 | 2000 continue
|
|---|
| 98 |
|
|---|
| 99 | endif
|
|---|
| 100 |
|
|---|
| 101 | return
|
|---|
| 102 | end
|
|---|
| 103 |
|
|---|
| 104 | c-----------------------------------------------------------------------
|
|---|
| 105 | subroutine makeuq
|
|---|
| 106 |
|
|---|
| 107 | c Fill up user defined forcing function and collocate will the
|
|---|
| 108 | c mass matrix on the Gauss-Lobatto mesh.
|
|---|
| 109 |
|
|---|
| 110 | include 'SIZE'
|
|---|
| 111 | include 'INPUT'
|
|---|
| 112 | include 'MASS'
|
|---|
| 113 | include 'SOLN'
|
|---|
| 114 | include 'TSTEP'
|
|---|
| 115 |
|
|---|
| 116 | n = lx1*ly1*lz1*nelfld(ifield)
|
|---|
| 117 |
|
|---|
| 118 | if (.not.ifcvfld(ifield)) time = time-dt ! Set time to t^n-1 for user function
|
|---|
| 119 |
|
|---|
| 120 | if (nio.eq.0.and.loglevel.gt.2)
|
|---|
| 121 | $ write(6,*) 'makeuq', ifield, time
|
|---|
| 122 | call setqvol(bq(1,1,1,1,ifield-1))
|
|---|
| 123 | call col2 (bq(1,1,1,1,ifield-1) ,bm1,n)
|
|---|
| 124 |
|
|---|
| 125 | if (.not.ifcvfld(ifield)) time = time+dt ! Restore time
|
|---|
| 126 |
|
|---|
| 127 | return
|
|---|
| 128 | end
|
|---|
| 129 | c-----------------------------------------------------------------------
|
|---|
| 130 | subroutine setqvol(bql)
|
|---|
| 131 |
|
|---|
| 132 | c Set user specified volumetric forcing function (e.g. heat source).
|
|---|
| 133 |
|
|---|
| 134 | include 'SIZE'
|
|---|
| 135 | include 'INPUT'
|
|---|
| 136 | include 'SOLN'
|
|---|
| 137 | include 'TSTEP'
|
|---|
| 138 | include 'CTIMER'
|
|---|
| 139 |
|
|---|
| 140 | real bql(lx1*ly1*lz1,lelt)
|
|---|
| 141 |
|
|---|
| 142 | #ifdef TIMER
|
|---|
| 143 | etime1=dnekclock()
|
|---|
| 144 | #endif
|
|---|
| 145 |
|
|---|
| 146 | nel = nelfld(ifield)
|
|---|
| 147 | nxyz1 = lx1*ly1*lz1
|
|---|
| 148 | n = nxyz1*nel
|
|---|
| 149 |
|
|---|
| 150 | do iel=1,nel
|
|---|
| 151 |
|
|---|
| 152 | call nekuq (bql,iel) ! ONLY SUPPORT USERQ - pff, 3/08/16
|
|---|
| 153 |
|
|---|
| 154 | c igrp = igroup(iel)
|
|---|
| 155 | c if (matype(igrp,ifield).eq.1) then ! constant source within a group
|
|---|
| 156 | c cqvol = cpgrp(igrp,ifield,3)
|
|---|
| 157 | c call cfill (bql(1,iel),cqvol,nxyz1)
|
|---|
| 158 | c else ! pff 2/6/96 ............ default is to look at userq
|
|---|
| 159 | c call nekuq (bql,iel)
|
|---|
| 160 | c endif
|
|---|
| 161 |
|
|---|
| 162 | enddo
|
|---|
| 163 | c
|
|---|
| 164 | c 101 FORMAT(' Wrong material type (',I3,') for group',I3,', field',I2
|
|---|
| 165 | c $ ,/,' Aborting in SETQVOL.')
|
|---|
| 166 | C
|
|---|
| 167 |
|
|---|
| 168 | #ifdef TIMER
|
|---|
| 169 | tusfq=tusfq+(dnekclock()-etime1)
|
|---|
| 170 | #endif
|
|---|
| 171 |
|
|---|
| 172 | return
|
|---|
| 173 | end
|
|---|
| 174 | C
|
|---|
| 175 | subroutine nekuq (bql,iel)
|
|---|
| 176 | C------------------------------------------------------------------
|
|---|
| 177 | C
|
|---|
| 178 | C Generate user-specified volumetric source term (temp./p.s.)
|
|---|
| 179 | C
|
|---|
| 180 | C------------------------------------------------------------------
|
|---|
| 181 | include 'SIZE'
|
|---|
| 182 | include 'SOLN'
|
|---|
| 183 | include 'MASS'
|
|---|
| 184 | include 'PARALLEL'
|
|---|
| 185 | include 'TSTEP'
|
|---|
| 186 | include 'NEKUSE'
|
|---|
| 187 | include 'INPUT'
|
|---|
| 188 | c
|
|---|
| 189 | real bql(lx1,ly1,lz1,lelt)
|
|---|
| 190 | c
|
|---|
| 191 | ielg = lglel(iel)
|
|---|
| 192 | do 10 k=1,lz1
|
|---|
| 193 | do 10 j=1,ly1
|
|---|
| 194 | do 10 i=1,lx1
|
|---|
| 195 | if (optlevel.le.2) call nekasgn (i,j,k,iel)
|
|---|
| 196 | qvol = 0.0
|
|---|
| 197 | call userq (i,j,k,ielg)
|
|---|
| 198 | bql(i,j,k,iel) = qvol
|
|---|
| 199 | 10 continue
|
|---|
| 200 |
|
|---|
| 201 | return
|
|---|
| 202 | end
|
|---|
| 203 | c-----------------------------------------------------------------------
|
|---|
| 204 | subroutine convab
|
|---|
| 205 | C---------------------------------------------------------------
|
|---|
| 206 | C
|
|---|
| 207 | C Eulerian scheme, add convection term to forcing function
|
|---|
| 208 | C at current time step.
|
|---|
| 209 | C
|
|---|
| 210 | C---------------------------------------------------------------
|
|---|
| 211 | include 'SIZE'
|
|---|
| 212 | include 'SOLN'
|
|---|
| 213 | include 'MASS'
|
|---|
| 214 | include 'TSTEP'
|
|---|
| 215 |
|
|---|
| 216 | common /scruz/ ta (lx1*ly1*lz1*lelt)
|
|---|
| 217 |
|
|---|
| 218 | nel = nelfld(ifield)
|
|---|
| 219 | n = lx1*ly1*lz1*nel
|
|---|
| 220 |
|
|---|
| 221 | call convop (ta,t(1,1,1,1,ifield-1))
|
|---|
| 222 | do i=1,n
|
|---|
| 223 | bq(i,1,1,1,ifield-1) = bq (i,1,1,1,ifield-1)
|
|---|
| 224 | $ - bm1(i,1,1,1)*ta(i)*vtrans(i,1,1,1,ifield)
|
|---|
| 225 | enddo
|
|---|
| 226 |
|
|---|
| 227 | return
|
|---|
| 228 | end
|
|---|
| 229 | c-----------------------------------------------------------------------
|
|---|
| 230 | subroutine makeabq
|
|---|
| 231 | C
|
|---|
| 232 | C Sum up contributions to 3rd order Adams-Bashforth scheme.
|
|---|
| 233 | C
|
|---|
| 234 | include 'SIZE'
|
|---|
| 235 | include 'SOLN'
|
|---|
| 236 | include 'TSTEP'
|
|---|
| 237 |
|
|---|
| 238 | ab0 = ab(1)
|
|---|
| 239 | ab1 = ab(2)
|
|---|
| 240 | ab2 = ab(3)
|
|---|
| 241 | nel = nelfld(ifield)
|
|---|
| 242 | n = lx1*ly1*lz1*nel
|
|---|
| 243 |
|
|---|
| 244 | do i=1,n
|
|---|
| 245 | ta=ab1*vgradt1(i,1,1,1,ifield-1)+ab2*vgradt2(i,1,1,1,ifield-1)
|
|---|
| 246 | vgradt2(i,1,1,1,ifield-1)=vgradt1(i,1,1,1,ifield-1)
|
|---|
| 247 | vgradt1(i,1,1,1,ifield-1)=bq (i,1,1,1,ifield-1)
|
|---|
| 248 | bq (i,1,1,1,ifield-1)=bq (i,1,1,1,ifield-1)*ab0+ta
|
|---|
| 249 | enddo
|
|---|
| 250 |
|
|---|
| 251 | return
|
|---|
| 252 | end
|
|---|
| 253 | c-----------------------------------------------------------------------
|
|---|
| 254 | subroutine makebdq
|
|---|
| 255 | C-----------------------------------------------------------------------
|
|---|
| 256 | C
|
|---|
| 257 | C Add contributions to F from lagged BD terms.
|
|---|
| 258 | C
|
|---|
| 259 | C-----------------------------------------------------------------------
|
|---|
| 260 | include 'SIZE'
|
|---|
| 261 | include 'TOTAL'
|
|---|
| 262 |
|
|---|
| 263 | parameter (lt=lx1*ly1*lz1*lelt)
|
|---|
| 264 | common /scrns/ tb(lt),h2(lt)
|
|---|
| 265 |
|
|---|
| 266 | nel = nelfld(ifield)
|
|---|
| 267 | n = lx1*ly1*lz1*nel
|
|---|
| 268 |
|
|---|
| 269 | const = 1./dt
|
|---|
| 270 | do i=1,n
|
|---|
| 271 | h2(i)=const*vtrans(i,1,1,1,ifield)
|
|---|
| 272 | tb(i)=bd(2)*bm1(i,1,1,1)*t(i,1,1,1,ifield-1)
|
|---|
| 273 | enddo
|
|---|
| 274 |
|
|---|
| 275 | do ilag=2,nbd
|
|---|
| 276 | if (ifgeom) then
|
|---|
| 277 | do i=1,n
|
|---|
| 278 | ta=bm1lag(i,1,1,1,ilag-1)*tlag(i,1,1,1,ilag-1,ifield-1)
|
|---|
| 279 | tb(i)=tb(i)+ta*bd(ilag+1)
|
|---|
| 280 | enddo
|
|---|
| 281 | else
|
|---|
| 282 | do i=1,n
|
|---|
| 283 | ta=bm1(i,1,1,1)*tlag(i,1,1,1,ilag-1,ifield-1)
|
|---|
| 284 | tb(i)=tb(i)+ta*bd(ilag+1)
|
|---|
| 285 | enddo
|
|---|
| 286 | endif
|
|---|
| 287 | enddo
|
|---|
| 288 |
|
|---|
| 289 | call addcol3 (bq(1,1,1,1,ifield-1),tb,h2,n)
|
|---|
| 290 |
|
|---|
| 291 | return
|
|---|
| 292 | end
|
|---|
| 293 | c-----------------------------------------------------------------------
|
|---|
| 294 | subroutine lagscal ! Keep old passive scalar field(s)
|
|---|
| 295 |
|
|---|
| 296 | include 'SIZE'
|
|---|
| 297 | include 'TOTAL'
|
|---|
| 298 |
|
|---|
| 299 | n = lx1*ly1*lz1*nelfld(ifield)
|
|---|
| 300 |
|
|---|
| 301 | do ilag=nbdinp-1,2,-1
|
|---|
| 302 | call copy (tlag(1,1,1,1,ilag ,ifield-1),
|
|---|
| 303 | $ tlag(1,1,1,1,ilag-1,ifield-1),n)
|
|---|
| 304 | enddo
|
|---|
| 305 |
|
|---|
| 306 | call copy (tlag(1,1,1,1,1,ifield-1),t(1,1,1,1,ifield-1),n)
|
|---|
| 307 |
|
|---|
| 308 | return
|
|---|
| 309 | end
|
|---|
| 310 | c-----------------------------------------------------------------------
|
|---|
| 311 | subroutine outfldrq (x,txt10,ichk)
|
|---|
| 312 | include 'SIZE'
|
|---|
| 313 | include 'TSTEP'
|
|---|
| 314 | real x(lx1,ly1,lz1,lelt)
|
|---|
| 315 | character*10 txt10
|
|---|
| 316 | c
|
|---|
| 317 | integer idum,e
|
|---|
| 318 | save idum
|
|---|
| 319 | data idum /-3/
|
|---|
| 320 |
|
|---|
| 321 | if (idum.lt.0) return
|
|---|
| 322 | c
|
|---|
| 323 | C
|
|---|
| 324 | mtot = lx1*ly1*lz1*nelv
|
|---|
| 325 | if (lx1.gt.8.or.nelv.gt.16) return
|
|---|
| 326 | xmin = glmin(x,mtot)
|
|---|
| 327 | xmax = glmax(x,mtot)
|
|---|
| 328 | c
|
|---|
| 329 | nell = nelt
|
|---|
| 330 | rnel = nell
|
|---|
| 331 | snel = sqrt(rnel)+.1
|
|---|
| 332 | ne = snel
|
|---|
| 333 | ne1 = nell-ne+1
|
|---|
| 334 | k = 1
|
|---|
| 335 | do ie=1,1
|
|---|
| 336 | ne = 0
|
|---|
| 337 | write(6,116) txt10,k,ie,xmin,xmax,istep,time
|
|---|
| 338 | do l=0,1
|
|---|
| 339 | write(6,117)
|
|---|
| 340 | do j=ly1,1,-1
|
|---|
| 341 | if (lx1.eq.2) write(6,102) ((x(i,j,k,e+l),i=1,lx1),e=1,1)
|
|---|
| 342 | if (lx1.eq.3) write(6,103) ((x(i,j,k,e+l),i=1,lx1),e=1,1)
|
|---|
| 343 | if (lx1.eq.4) write(6,104) ((x(i,j,k,e+l),i=1,lx1),e=1,1)
|
|---|
| 344 | if (lx1.eq.5) write(6,105) ((x(i,j,k,e+l),i=1,lx1),e=1,1)
|
|---|
| 345 | if (lx1.eq.6) write(6,106) ((x(i,j,k,e+l),i=1,lx1),e=1,1)
|
|---|
| 346 | if (lx1.eq.7) write(6,107) ((x(i,j,k,e+l),i=1,lx1),e=1,1)
|
|---|
| 347 | if (lx1.eq.8) write(6,118) ((x(i,j,k,e+l),i=1,lx1),e=1,1)
|
|---|
| 348 | enddo
|
|---|
| 349 | enddo
|
|---|
| 350 | enddo
|
|---|
| 351 | C
|
|---|
| 352 | 102 FORMAT(4(2f9.5,2x))
|
|---|
| 353 | 103 FORMAT(4(3f9.5,2x))
|
|---|
| 354 | 104 FORMAT(4(4f7.3,2x))
|
|---|
| 355 | 105 FORMAT(5f9.5,10x,5f9.5)
|
|---|
| 356 | 106 FORMAT(6f9.5,5x,6f9.5)
|
|---|
| 357 | 107 FORMAT(7f8.4,5x,7f8.4)
|
|---|
| 358 | 108 FORMAT(8f8.4,4x,8f8.4)
|
|---|
| 359 | 118 FORMAT(8f12.9)
|
|---|
| 360 | c
|
|---|
| 361 | 116 FORMAT( /,5X,' ^ ',/,
|
|---|
| 362 | $ 5X,' Y | ',/,
|
|---|
| 363 | $ 5X,' | ',A10,/,
|
|---|
| 364 | $ 5X,' +----> ','Plane = ',I2,'/',I2,2x,2e12.4,/,
|
|---|
| 365 | $ 5X,' X ','Step =',I9,f15.5)
|
|---|
| 366 | 117 FORMAT(' ')
|
|---|
| 367 | c
|
|---|
| 368 | if (ichk.eq.1.and.idum.gt.0) call checkit(idum)
|
|---|
| 369 | return
|
|---|
| 370 | end
|
|---|
| 371 | c-----------------------------------------------------------------------
|
|---|
| 372 | subroutine cdscal_expl (igeom)
|
|---|
| 373 | C
|
|---|
| 374 | C explicit convection-diffusion equation for passive scalar
|
|---|
| 375 | C
|
|---|
| 376 | include 'SIZE'
|
|---|
| 377 | include 'INPUT'
|
|---|
| 378 | include 'GEOM'
|
|---|
| 379 | include 'MVGEOM'
|
|---|
| 380 | include 'SOLN'
|
|---|
| 381 | include 'MASS'
|
|---|
| 382 | include 'TSTEP'
|
|---|
| 383 | common /cprint/ ifprint
|
|---|
| 384 | logical ifprint
|
|---|
| 385 | logical ifconv
|
|---|
| 386 |
|
|---|
| 387 | common /scrns/ ta(lx1,ly1,lz1,lelt)
|
|---|
| 388 | $ ,tb(lx1,ly1,lz1,lelt)
|
|---|
| 389 | common /scrvh/ h1(lx1,ly1,lz1,lelt)
|
|---|
| 390 | $ ,h2(lx1,ly1,lz1,lelt)
|
|---|
| 391 |
|
|---|
| 392 |
|
|---|
| 393 | c QUESTIONABLE support for Robin BC's at this point! (5/15/08)
|
|---|
| 394 |
|
|---|
| 395 | nel = nelfld(ifield)
|
|---|
| 396 | n = lx1*ly1*lz1*nel
|
|---|
| 397 |
|
|---|
| 398 | if (igeom.eq.1) then ! geometry at t^{n-1}
|
|---|
| 399 |
|
|---|
| 400 | call makeq
|
|---|
| 401 | call lagscal
|
|---|
| 402 |
|
|---|
| 403 | else ! geometry at t^n
|
|---|
| 404 |
|
|---|
| 405 | if (nio.eq.0) write (6,*) istep,ifield,' explicit step'
|
|---|
| 406 |
|
|---|
| 407 |
|
|---|
| 408 | C New geometry
|
|---|
| 409 |
|
|---|
| 410 | isd = 1
|
|---|
| 411 | if (ifaxis.and.ifmhd) isd = 2 !This is a problem if T is to be T!
|
|---|
| 412 |
|
|---|
| 413 | intype = 0
|
|---|
| 414 | if (iftran) intype = -1
|
|---|
| 415 | call sethlm (h1,h2,intype)
|
|---|
| 416 |
|
|---|
| 417 | call bcneusc (ta,-1) ! Modify diagonal for Robin condition
|
|---|
| 418 | call add2 (h2,ta ,n)
|
|---|
| 419 | call col2 (h2,BM1,n)
|
|---|
| 420 |
|
|---|
| 421 | call bcneusc (tb,1) ! Modify rhs for flux bc
|
|---|
| 422 | call add2 (bq(1,1,1,1,ifield-1),tb,n)
|
|---|
| 423 |
|
|---|
| 424 | call dssum (bq(1,1,1,1,ifield-1),lx1,ly1,lz1)
|
|---|
| 425 | call dssum (h2,lx1,ly1,lz1)
|
|---|
| 426 |
|
|---|
| 427 | call invcol3 (t(1,1,1,1,ifield-1),bq(1,1,1,1,ifield-1),h2,n)
|
|---|
| 428 |
|
|---|
| 429 | call bcdirsc (t(1,1,1,1,ifield-1)) ! --> no mask needed
|
|---|
| 430 |
|
|---|
| 431 | endif ! geometry at t^n
|
|---|
| 432 |
|
|---|
| 433 | return
|
|---|
| 434 | end
|
|---|
| 435 | c-----------------------------------------------------------------------
|
|---|
| 436 | subroutine diffab ! explicit treatment of diffusion operator
|
|---|
| 437 | c
|
|---|
| 438 | c Eulerian scheme, add diffusion term to forcing function
|
|---|
| 439 | c at current time step.
|
|---|
| 440 | c
|
|---|
| 441 |
|
|---|
| 442 | include 'SIZE'
|
|---|
| 443 | include 'SOLN'
|
|---|
| 444 | include 'MASS'
|
|---|
| 445 | include 'TSTEP'
|
|---|
| 446 | include 'INPUT'
|
|---|
| 447 |
|
|---|
| 448 | common /scruz/ ta(lx1,ly1,lz1,lelt)
|
|---|
| 449 | $ ,h2(lx1,ly1,lz1,lelt)
|
|---|
| 450 |
|
|---|
| 451 | nel = nelfld(ifield)
|
|---|
| 452 | n = lx1*ly1*lz1*nel
|
|---|
| 453 |
|
|---|
| 454 | intype = 0
|
|---|
| 455 | if (iftran) intype = -1
|
|---|
| 456 |
|
|---|
| 457 | isd = 1
|
|---|
| 458 | if (ifaxis.and.ifmhd) isd = 2 !This is a problem if T is to be T!
|
|---|
| 459 |
|
|---|
| 460 | imesh = 1
|
|---|
| 461 | c if (iftmsh(ifield)) imesh=2
|
|---|
| 462 |
|
|---|
| 463 | call rzero (h2,n)
|
|---|
| 464 | call axhelm (ta,t(1,1,1,1,ifield-1),vdiff(1,1,1,1,ifield)
|
|---|
| 465 | $ ,h2,imesh,isd)
|
|---|
| 466 | call sub2 (bq(1,1,1,1,ifield-1),ta,n)
|
|---|
| 467 |
|
|---|
| 468 | return
|
|---|
| 469 | end
|
|---|
| 470 | c-----------------------------------------------------------------------
|
|---|
| 471 | subroutine set_eta_alpha2
|
|---|
| 472 |
|
|---|
| 473 | c Set up required dg terms, e.g., alpha, eta, etc.
|
|---|
| 474 | c Face weight: .5 interior, 1. boundary
|
|---|
| 475 |
|
|---|
| 476 | include 'SIZE'
|
|---|
| 477 | include 'TOTAL'
|
|---|
| 478 | common /mysedg/ eta
|
|---|
| 479 |
|
|---|
| 480 |
|
|---|
| 481 | integer e,f,pf
|
|---|
| 482 |
|
|---|
| 483 | nface = 2*ldim
|
|---|
| 484 | call dsset(lx1,ly1,lz1)
|
|---|
| 485 |
|
|---|
| 486 | eta = 5 ! Semi-optimized value, single domain
|
|---|
| 487 |
|
|---|
| 488 | do e=1,nelfld(ifield)
|
|---|
| 489 | do f=1,nface
|
|---|
| 490 |
|
|---|
| 491 | pf = eface1(f)
|
|---|
| 492 | js1 = skpdat(1,pf)
|
|---|
| 493 | jf1 = skpdat(2,pf)
|
|---|
| 494 | jskip1 = skpdat(3,pf)
|
|---|
| 495 | js2 = skpdat(4,pf)
|
|---|
| 496 | jf2 = skpdat(5,pf)
|
|---|
| 497 | jskip2 = skpdat(6,pf)
|
|---|
| 498 |
|
|---|
| 499 | i = 0
|
|---|
| 500 | do j2=js2,jf2,jskip2
|
|---|
| 501 | do j1=js1,jf1,jskip1
|
|---|
| 502 | i = i+1
|
|---|
| 503 | a = area(i,1,f,e) ! Check Fydkowski notes
|
|---|
| 504 | a = a*a*fw(f,e) ! For ds_avg used below, plus quad weight
|
|---|
| 505 | etalph(i,f,e) = eta*(a/bm1(j1,j2,1,e))
|
|---|
| 506 | c write(6,*) i,j1,j2,e,f,a,etalph(i,f,e)
|
|---|
| 507 | enddo
|
|---|
| 508 | enddo
|
|---|
| 509 | enddo
|
|---|
| 510 | enddo
|
|---|
| 511 |
|
|---|
| 512 | call fgslib_gs_op (dg_hndlx,etalph,1,1,0) ! 1 ==> +
|
|---|
| 513 |
|
|---|
| 514 | return
|
|---|
| 515 | end
|
|---|
| 516 | c-----------------------------------------------------------------------
|
|---|
| 517 | subroutine fwght(msk1,mult)
|
|---|
| 518 | include 'SIZE'
|
|---|
| 519 | include 'TOTAL'
|
|---|
| 520 | parameter (lx=lx1*ly1*lz1)
|
|---|
| 521 | real msk1(lx1,ly1,lz1),mult(lx1,ly1,lz1,1)
|
|---|
| 522 | integer e,f
|
|---|
| 523 |
|
|---|
| 524 | parameter(lf=lx1*lz1*2*ldim*lelt)
|
|---|
| 525 | common /scrdg/uf(lx1*lz1,2*ldim,lelt)
|
|---|
| 526 |
|
|---|
| 527 |
|
|---|
| 528 |
|
|---|
| 529 | do i=1,lf
|
|---|
| 530 | uf(i,1,1)=1.
|
|---|
| 531 | enddo
|
|---|
| 532 | call fgslib_gs_op (dg_hndlx,uf,1,1,0) ! 1 ==> +
|
|---|
| 533 |
|
|---|
| 534 | nface = 2*ldim
|
|---|
| 535 | do e=1,nelt
|
|---|
| 536 | do f=1,nface
|
|---|
| 537 | fw(f,e) = 1 ! Boundary
|
|---|
| 538 | if (uf(1,f,e).gt.1.1) fw(f,e)=0.5 ! Interior
|
|---|
| 539 | enddo
|
|---|
| 540 | enddo
|
|---|
| 541 |
|
|---|
| 542 | return
|
|---|
| 543 | end
|
|---|
| 544 | c-----------------------------------------------------------------------
|
|---|
| 545 | subroutine dg_setup
|
|---|
| 546 | include 'SIZE'
|
|---|
| 547 | include 'TOTAL'
|
|---|
| 548 |
|
|---|
| 549 | common /ivrtx/ vertex ((2**ldim)*lelt)
|
|---|
| 550 | common /ctmp1/ qs(lx1*ly1*lz1*lelt)
|
|---|
| 551 | integer vertex
|
|---|
| 552 |
|
|---|
| 553 | logical ifany
|
|---|
| 554 | save ifany
|
|---|
| 555 | data ifany /.false./
|
|---|
| 556 |
|
|---|
| 557 | do ifield=1,ldimt1
|
|---|
| 558 | if (ifdgfld(ifield)) ifany=.true.
|
|---|
| 559 | enddo
|
|---|
| 560 |
|
|---|
| 561 | if (ifany) then
|
|---|
| 562 |
|
|---|
| 563 | call setup_dg_gs(dg_hndlx,lx1,ly1,lz1,nelt,nelgt,vertex)
|
|---|
| 564 | call dg_set_fc_ptr
|
|---|
| 565 |
|
|---|
| 566 | param(59)=1
|
|---|
| 567 | call geom_reset(1)
|
|---|
| 568 | call set_unr
|
|---|
| 569 | endif
|
|---|
| 570 |
|
|---|
| 571 |
|
|---|
| 572 | n = lx1*ly1*lz1*nelt
|
|---|
| 573 | call invers2(binvdg,bm1,n)
|
|---|
| 574 |
|
|---|
| 575 | call set_dg_wgts
|
|---|
| 576 |
|
|---|
| 577 | return
|
|---|
| 578 | end
|
|---|
| 579 | c-----------------------------------------------------------------------
|
|---|
| 580 | subroutine dg_setup2(mask)
|
|---|
| 581 | include 'SIZE'
|
|---|
| 582 | include 'TOTAL'
|
|---|
| 583 |
|
|---|
| 584 | real mask(1)
|
|---|
| 585 | common /ctmp0/ qs(lx1*ly1*lz1*lelt)
|
|---|
| 586 |
|
|---|
| 587 | integer ifld_last
|
|---|
| 588 | save ifld_last
|
|---|
| 589 | data ifld_last /0/
|
|---|
| 590 |
|
|---|
| 591 | if (ifield.eq.ifld_last) return
|
|---|
| 592 | ifld_last = ifield
|
|---|
| 593 |
|
|---|
| 594 | call fwght (qs,mask) ! qs is work array
|
|---|
| 595 | call set_eta_alpha2 ! Set eta/h
|
|---|
| 596 |
|
|---|
| 597 | return
|
|---|
| 598 | end
|
|---|
| 599 | c-----------------------------------------------------------------------
|
|---|
| 600 | subroutine cdscal_dg (igeom)
|
|---|
| 601 | C
|
|---|
| 602 | C Solve the convection-diffusion equation for passive scalar IPSCAL
|
|---|
| 603 | C
|
|---|
| 604 | include 'SIZE'
|
|---|
| 605 | include 'TOTAL'
|
|---|
| 606 | common /cprint/ ifprint
|
|---|
| 607 | logical ifprint,ifconv
|
|---|
| 608 |
|
|---|
| 609 | parameter (lt=lx1*ly1*lz1*lelt)
|
|---|
| 610 | common /scrns/ ta(lt),tb(lt)
|
|---|
| 611 | common /scrvh/ h1(lt),h2(lt)
|
|---|
| 612 |
|
|---|
| 613 | include 'ORTHOT' ! This must be fixed
|
|---|
| 614 |
|
|---|
| 615 | call dg_setup2(tmask(1,1,1,1,ifield-1))
|
|---|
| 616 |
|
|---|
| 617 | nel = nelfld(ifield)
|
|---|
| 618 | n = lx1*ly1*lz1*nel
|
|---|
| 619 |
|
|---|
| 620 | if (igeom.eq.1) then ! old geometry at t^{n-1}
|
|---|
| 621 |
|
|---|
| 622 | call makeq
|
|---|
| 623 | call lagscal
|
|---|
| 624 |
|
|---|
| 625 | else ! new geometry at t^n
|
|---|
| 626 |
|
|---|
| 627 | if (ifprint.and.nio.eq.0)
|
|---|
| 628 | $ write (6,*) ' Temperature/Passive scalar solution',ifield
|
|---|
| 629 |
|
|---|
| 630 | if1=ifield-1
|
|---|
| 631 | write(name4t,1) if1-1
|
|---|
| 632 | 1 format('PS',i2)
|
|---|
| 633 | if (ifield.eq.2) write(name4t,'(A4)') 'TEMP'
|
|---|
| 634 |
|
|---|
| 635 | isd = 1
|
|---|
| 636 | if (ifaxis.and.ifaziv.and.ifield.eq.2) isd = 2
|
|---|
| 637 | c if (ifaxis.and.ifmhd) isd = 2 !This is a problem if T is to be T!
|
|---|
| 638 |
|
|---|
| 639 | intype = 0
|
|---|
| 640 | if (iftran) intype = -1
|
|---|
| 641 | call sethlm (h1,h2,intype)
|
|---|
| 642 |
|
|---|
| 643 | c call bcneusc (ta,-1) !! Not Yet supported for DG
|
|---|
| 644 | c call add2 (h2,ta,n) !! Not Yet supported for DG
|
|---|
| 645 |
|
|---|
| 646 | call rzero (tb,n)
|
|---|
| 647 | call bcdirsc ( t (1,1,1,1,ifield-1))
|
|---|
| 648 | call conv_bdry_dg_weak (tb,t (1,1,1,1,ifield-1))
|
|---|
| 649 | call hxdg_surfa (tb,t (1,1,1,1,ifield-1),h1,h2)
|
|---|
| 650 | call add2 (tb,bq(1,1,1,1,ifield-1),n)
|
|---|
| 651 |
|
|---|
| 652 | call hmholtz_dg(name4t,t(1,1,1,1,ifield-1),tb,h1,h2
|
|---|
| 653 | $ ,tmask(1,1,1,1,ifield-1)
|
|---|
| 654 | $ ,tolht(ifield),nmxt(ifield-1))
|
|---|
| 655 |
|
|---|
| 656 | endif ! End of IGEOM branch.
|
|---|
| 657 |
|
|---|
| 658 | return
|
|---|
| 659 | end
|
|---|
| 660 | c-----------------------------------------------------------------------
|
|---|
