| 1 | c-----------------------------------------------------------------------
|
|---|
| 2 | subroutine char_conv(p0,u,ulag,bm,bmlag,msk,c,cs,gsl)
|
|---|
| 3 | c
|
|---|
| 4 | c
|
|---|
| 5 | c Convect over last NBD steps using characteristics scheme
|
|---|
| 6 | c
|
|---|
| 7 | c NOTE: Here, we assume that ulag is stored by time-slice first,
|
|---|
| 8 | c then by field number (this is opposite to prior Nek5000)
|
|---|
| 9 | c
|
|---|
| 10 | c
|
|---|
| 11 | include 'SIZE'
|
|---|
| 12 | include 'TOTAL'
|
|---|
| 13 | real p0(1),u(1),ulag(1),bm(1),bmlag(1),msk(1),c(1),cs(0:1)
|
|---|
| 14 | integer gsl
|
|---|
| 15 |
|
|---|
| 16 | common /scrns/ ct (lxd*lyd*lzd*lelv*ldim)
|
|---|
| 17 |
|
|---|
| 18 | common /scrvh/ bmsk(lx1*ly1*lz1*lelv)
|
|---|
| 19 | $ , bdwt(lx1*ly1*lz1*lelv)
|
|---|
| 20 | $ , bmst(lx1*ly1*lz1*lelv)
|
|---|
| 21 | $ , u1 (lx1*ly1*lz1*lelv)
|
|---|
| 22 |
|
|---|
| 23 | common /scrmg/ r1 (lx1*ly1*lz1*lelv)
|
|---|
| 24 | $ , r2 (lx1*ly1*lz1*lelv)
|
|---|
| 25 | $ , r3 (lx1*ly1*lz1*lelv)
|
|---|
| 26 | $ , r4 (lx1*ly1*lz1*lelv)
|
|---|
| 27 |
|
|---|
| 28 | nelc = nelv ! number of elements in convecting field
|
|---|
| 29 | if (ifield.eq.ifldmhd) nelc = nelfld(ifield)
|
|---|
| 30 |
|
|---|
| 31 | nc = cs(0) ! number of stored convecting fields
|
|---|
| 32 |
|
|---|
| 33 | ln = lx1*ly1*lz1*lelt
|
|---|
| 34 | n = lx1*ly1*lz1*nelfld(ifield)
|
|---|
| 35 | m = lxd*lyd*lzd*nelc*ldim
|
|---|
| 36 |
|
|---|
| 37 | c if(nid.eq.0) write(*,*) 'going into char_conv1 '
|
|---|
| 38 | call char_conv1 (p0,u,bmnv,n,ulag,ln,gsl,c,m,cs(1),nc,ct
|
|---|
| 39 | $ ,u1,r1,r2,r3,r4,bmsk,bdivw,bdwt,bmass,bmst,bm,bmlag)
|
|---|
| 40 |
|
|---|
| 41 | return
|
|---|
| 42 | end
|
|---|
| 43 | c-----------------------------------------------------------------------
|
|---|
| 44 | subroutine char_conv1 (p0,u,bmnv,n,ulag,ln,gsl,c,m,cs,nc,ct
|
|---|
| 45 | $ ,u1,r1,r2,r3,r4,bmsk,bdivw,bdwt,bmass,bmst,bm,bmlag)
|
|---|
| 46 |
|
|---|
| 47 | include 'SIZE'
|
|---|
| 48 | include 'INPUT'
|
|---|
| 49 | include 'TSTEP'
|
|---|
| 50 |
|
|---|
| 51 | real p0(n),u(n),bmnv(n,1),ulag(ln,1),c(m,0:nc),cs(0:nc),bdivw(n,1)
|
|---|
| 52 | $ ,bm(n), bmlag(ln,1)
|
|---|
| 53 |
|
|---|
| 54 | real ct(m),u1(n),r1(n),r2(n),r3(n),r4(n),bmsk(n),bdwt(n) ! work arrays
|
|---|
| 55 |
|
|---|
| 56 | integer gsl
|
|---|
| 57 |
|
|---|
| 58 |
|
|---|
| 59 | ! Convect over last NBD steps using characteristics scheme
|
|---|
| 60 |
|
|---|
| 61 | ! n-q n-1
|
|---|
| 62 | ! Given u(t , X ), q = 1,2,...,nbd, compute phi ( := p0 )
|
|---|
| 63 |
|
|---|
| 64 | ! n-1 nbd ~n-q
|
|---|
| 65 | ! phi := sum u
|
|---|
| 66 | ! q=1
|
|---|
| 67 |
|
|---|
| 68 | ! ~n-q ~n-q
|
|---|
| 69 | ! each u satisfies u := v such that
|
|---|
| 70 |
|
|---|
| 71 |
|
|---|
| 72 | ! dv
|
|---|
| 73 | ! -- + C.grad v = 0 t \in [t^n-q,t^n], v(t^n-q,X) = u(t^n-q,X)
|
|---|
| 74 | ! dt
|
|---|
| 75 |
|
|---|
| 76 | ! n = lx1*ly1*lz1*nelv
|
|---|
| 77 | ! m = lxd*lyd*lzd*nelv
|
|---|
| 78 |
|
|---|
| 79 | tau = time-vlsum(dtlag,nbd) ! initialize time for u^n-k
|
|---|
| 80 | call int_vel (ct ,tau,c ,m,nc,cs,nid) ! ct(t) = sum w_k c(.,k)
|
|---|
| 81 | call int_vel (bmsk,tau,bmnv ,n,nc,cs,nid) ! B^-1(t^n-1)
|
|---|
| 82 | call int_vel (bmst,tau,bmass,n,nc,cs,nid) ! B(t^n-1)
|
|---|
| 83 | call int_vel (bdwt,tau,bdivw,n,nc,cs,nid) ! BdivW(t^n-1)
|
|---|
| 84 |
|
|---|
| 85 | call rzero(p0,n)
|
|---|
| 86 |
|
|---|
| 87 | do ilag = nbd,1,-1
|
|---|
| 88 |
|
|---|
| 89 | um = 0
|
|---|
| 90 | if (ilag.eq.1) then
|
|---|
| 91 | do i=1,n
|
|---|
| 92 | p0(i) = p0(i)+bd(ilag+1)*u(i)*bm(i)
|
|---|
| 93 | um=max(um,u(i))
|
|---|
| 94 | enddo
|
|---|
| 95 | else
|
|---|
| 96 | if(ifmvbd) then
|
|---|
| 97 | do i=1,n
|
|---|
| 98 | p0(i) = p0(i)+bd(ilag+1)*ulag(i,ilag-1)*bmlag(i,ilag-1)
|
|---|
| 99 | um=max(um,ulag(i,ilag-1))
|
|---|
| 100 | enddo
|
|---|
| 101 | else
|
|---|
| 102 | do i=1,n
|
|---|
| 103 | p0(i) = p0(i)+bd(ilag+1)*ulag(i,ilag-1)*bm(i)
|
|---|
| 104 | um=max(um,ulag(i,ilag-1))
|
|---|
| 105 | enddo
|
|---|
| 106 | endif
|
|---|
| 107 | endif
|
|---|
| 108 |
|
|---|
| 109 | dtau = dtlag(ilag)/ntaubd
|
|---|
| 110 | do itau = 1,ntaubd ! ntaubd=number of RK4 substeps (typ. 1 or 2)
|
|---|
| 111 |
|
|---|
| 112 | tau1 = tau + dtau
|
|---|
| 113 |
|
|---|
| 114 | c1 = 1.
|
|---|
| 115 | c2 = -dtau/2.
|
|---|
| 116 | c3 = -dtau
|
|---|
| 117 | th = tau+dtau/2.
|
|---|
| 118 |
|
|---|
| 119 | call invcol3 (u1,p0,bmst,n)
|
|---|
| 120 | call conv_rhs(r1,u1,ct,bmsk,bmst,bdwt,gsl) ! STAGE 1
|
|---|
| 121 | call col2 (r1,bmst,n) ! ! r1 = B(n-1)* r1
|
|---|
| 122 |
|
|---|
| 123 | call add3s12 (u1,p0,r1,c1,c2,n)
|
|---|
| 124 | call int_vel (bmst,th,bmass,n,nc,cs,nid) ! B(n-1/2)
|
|---|
| 125 | call invcol2 (u1,bmst,n) ! u2=B(n-1/2)
|
|---|
| 126 |
|
|---|
| 127 | call int_vel (ct ,th,c ,m,nc,cs,nid) ! STAGE 2
|
|---|
| 128 | call int_vel (bmsk,th,bmnv ,n,nc,cs,nid) ! B^-1(n-1/2)
|
|---|
| 129 | call int_vel (bdwt,th,bdivw,n,nc,cs,nid) ! BdivW(n-1/2)
|
|---|
| 130 | call conv_rhs(r2,u1,ct,bmsk,bmst,bdwt,gsl)
|
|---|
| 131 | call col2 (r2,bmst,n) ! du = B * du
|
|---|
| 132 |
|
|---|
| 133 | call add3s12 (u1,p0,r2,c1,c2,n) ! STAGE 3
|
|---|
| 134 | call invcol2 (u1,bmst,n)
|
|---|
| 135 | call conv_rhs(r3,u1,ct,bmsk,bmst,bdwt,gsl) ! B(n-1/2) (still)
|
|---|
| 136 | call col2 (r3,bmst,n) ! du = B * du
|
|---|
| 137 |
|
|---|
| 138 | call add3s12 (u1,p0,r3,c1,c3,n)
|
|---|
| 139 | call int_vel (bmst,tau1,bmass,n,nc,cs,nid) ! B^-1(n)
|
|---|
| 140 | call invcol2 (u1,bmst,n) ! u2=B(n-1/2)
|
|---|
| 141 |
|
|---|
| 142 | call int_vel (ct ,tau1,c ,m,nc,cs,nid) ! STAGE 4
|
|---|
| 143 | call int_vel (bmsk,tau1,bmnv ,n,nc,cs,nid) ! B^-1(n)
|
|---|
| 144 | call int_vel (bdwt,tau1,bdivw,n,nc,cs,nid) ! BdivW(n)
|
|---|
| 145 | call conv_rhs(r4,u1,ct,bmsk,bmst,bdwt,gsl)
|
|---|
| 146 | call col2 (r4,bmst,n) ! du = B * du
|
|---|
| 147 |
|
|---|
| 148 | c1 = -dtau/6.
|
|---|
| 149 | c2 = -dtau/3.
|
|---|
| 150 | do i=1,n
|
|---|
| 151 | p0(i) = p0(i)+c1*(r1(i)+r4(i))+c2*(r2(i)+r3(i))
|
|---|
| 152 | enddo
|
|---|
| 153 | tau = tau1
|
|---|
| 154 | enddo
|
|---|
| 155 | enddo
|
|---|
| 156 |
|
|---|
| 157 | return
|
|---|
| 158 | end
|
|---|
| 159 | c-----------------------------------------------------------------------
|
|---|
| 160 | subroutine int_vel(c_t,t0,c,n,nc,ct,nid)
|
|---|
| 161 |
|
|---|
| 162 | c Interpolate convecting velocity field c(t_1,...,t_nconv) to
|
|---|
| 163 | c time t0 and return result in c_t.
|
|---|
| 164 |
|
|---|
| 165 | c Ouput: c_t = sum wt_k * ct_i(k)
|
|---|
| 166 |
|
|---|
| 167 | c Here, t0 is the time of interest
|
|---|
| 168 |
|
|---|
| 169 | real c_t(n),c(n,0:nc),ct(0:nc)
|
|---|
| 170 |
|
|---|
| 171 | parameter (lwtmax=10)
|
|---|
| 172 | real wt(0:lwtmax)
|
|---|
| 173 |
|
|---|
| 174 | if (nc.gt.lwtmax) then
|
|---|
| 175 | write(6,*) nid,'ERROR int_vel: lwtmax too small',lwtmax,nc
|
|---|
| 176 | call exitt
|
|---|
| 177 | endif
|
|---|
| 178 |
|
|---|
| 179 | no = nc-1
|
|---|
| 180 | call fd_weights_full(t0,ct(0),no,0,wt) ! interpolation weights
|
|---|
| 181 |
|
|---|
| 182 | call rzero(c_t,n)
|
|---|
| 183 | do j=1,n
|
|---|
| 184 | do i=0,no
|
|---|
| 185 | c_t(j) = c_t(j) + wt(i)*c(j,i)
|
|---|
| 186 | enddo
|
|---|
| 187 | enddo
|
|---|
| 188 |
|
|---|
| 189 | return
|
|---|
| 190 | end
|
|---|
| 191 | c-----------------------------------------------------------------------
|
|---|
| 192 | subroutine conv_rhs (du,u,c,bmsk,bmst,bdwt,gsl)
|
|---|
| 193 | c
|
|---|
| 194 | include 'SIZE'
|
|---|
| 195 | include 'TOTAL'
|
|---|
| 196 | c
|
|---|
| 197 | c apply convecting field c(1,ldim) to scalar field u(1)
|
|---|
| 198 | c
|
|---|
| 199 | real du(1),u(1),c(1),bmsk(1),bdwt(1)
|
|---|
| 200 | integer gsl
|
|---|
| 201 | c
|
|---|
| 202 | logical ifconv
|
|---|
| 203 | c
|
|---|
| 204 | c ifconv = .false.
|
|---|
| 205 | ifconv = .true.
|
|---|
| 206 | c
|
|---|
| 207 | n = lx1*ly1*lz1*nelv
|
|---|
| 208 |
|
|---|
| 209 | if (ifdgfld(ifield)) then
|
|---|
| 210 |
|
|---|
| 211 | if (param(99).eq.1) call conv_rhs_dg (du,u,c)
|
|---|
| 212 | if (param(99).eq.0) call conv_rhs_dg_aliased (du,u,c)
|
|---|
| 213 |
|
|---|
| 214 | elseif (ifconv) then
|
|---|
| 215 |
|
|---|
| 216 | if (ifcons) then
|
|---|
| 217 | if (if3d ) call convop_cons_3d (du,u,c,lx1,lxd,nelv)
|
|---|
| 218 | if (.not.if3d) call convop_cons_2d (du,u,c,lx1,lxd,nelv)
|
|---|
| 219 | else
|
|---|
| 220 | if (if3d ) call convop_fst_3d (du,u,c,lx1,lxd,nelv)
|
|---|
| 221 | if (.not.if3d) call convop_fst_2d (du,u,c,lx1,lxd,nelv)
|
|---|
| 222 | endif
|
|---|
| 223 |
|
|---|
| 224 | call subcol3(du,bdwt,u,n)
|
|---|
| 225 | call fgslib_gs_op (gsl,du,1,1,0) ! +
|
|---|
| 226 |
|
|---|
| 227 | call col2 (du,bmsk,n) ! du = Binv * msk * du
|
|---|
| 228 |
|
|---|
| 229 | else
|
|---|
| 230 | call rzero (du,n)
|
|---|
| 231 | endif
|
|---|
| 232 |
|
|---|
| 233 | return
|
|---|
| 234 | end
|
|---|
| 235 | c-----------------------------------------------------------------------
|
|---|
| 236 | subroutine convop_fst_3d(du,u,c,mx,md,nel)
|
|---|
| 237 | c
|
|---|
| 238 | include 'SIZE'
|
|---|
| 239 | c
|
|---|
| 240 | c apply convecting field c to scalar field u
|
|---|
| 241 | c
|
|---|
| 242 | real du(mx*mx*mx,nel)
|
|---|
| 243 | real u(mx*mx*mx,nel)
|
|---|
| 244 | real c(md*md*md,nel,3)
|
|---|
| 245 | parameter (ldd=lxd*lyd*lzd)
|
|---|
| 246 | common /ctmp1/ ur(ldd),us(ldd),ut(ldd),ud(ldd)
|
|---|
| 247 | c
|
|---|
| 248 | logical if3d,ifd
|
|---|
| 249 | integer e
|
|---|
| 250 | c
|
|---|
| 251 | if3d = .true.
|
|---|
| 252 | ifd = .false.
|
|---|
| 253 | if (md.ne.mx) ifd=.true.
|
|---|
| 254 | ifd=.true.
|
|---|
| 255 | c
|
|---|
| 256 | nrstd = md**3
|
|---|
| 257 | call lim_chk(nrstd,ldd,'urus ','ldd ','convop_fst')
|
|---|
| 258 | c
|
|---|
| 259 | do e=1,nel
|
|---|
| 260 | call grad_rstd(ur,us,ut,u(1,e),mx,md,if3d,ud)
|
|---|
| 261 | if (ifd) then ! dealiased
|
|---|
| 262 | do i=1,nrstd
|
|---|
| 263 | c C has the mass matrix factored in per (4.8.5), p. 227, DFM.
|
|---|
| 264 | ud(i) = c(i,e,1)*ur(i)+c(i,e,2)*us(i)+c(i,e,3)*ut(i)
|
|---|
| 265 | enddo
|
|---|
| 266 | call intp_rstd(du(1,e),ud,mx,md,if3d,1) ! 1 --> transpose
|
|---|
| 267 | else
|
|---|
| 268 | do i=1,nrstd
|
|---|
| 269 | c C has the mass matrix factored in per (4.8.5), p. 227, DFM.
|
|---|
| 270 | du(i,e) = c(i,e,1)*ur(i)+c(i,e,2)*us(i)+c(i,e,3)*ut(i)
|
|---|
| 271 | enddo
|
|---|
| 272 | endif
|
|---|
| 273 | enddo
|
|---|
| 274 | c
|
|---|
| 275 | return
|
|---|
| 276 | end
|
|---|
| 277 | c-----------------------------------------------------------------------
|
|---|
| 278 | subroutine convop_fst_2d(du,u,c,mx,md,nel)
|
|---|
| 279 | c
|
|---|
| 280 | include 'SIZE'
|
|---|
| 281 | c
|
|---|
| 282 | c apply convecting field c to scalar field u
|
|---|
| 283 | c
|
|---|
| 284 | real du(mx*mx,nel)
|
|---|
| 285 | real u(mx*mx,nel)
|
|---|
| 286 | real c(md*md,nel,2)
|
|---|
| 287 | parameter (ldd=lxd*lyd*lzd)
|
|---|
| 288 | common /ctmp1/ ur(ldd),us(ldd),ut(ldd),ud(ldd)
|
|---|
| 289 | c
|
|---|
| 290 | logical if3d,ifd
|
|---|
| 291 | integer e
|
|---|
| 292 | c
|
|---|
| 293 | if3d = .false.
|
|---|
| 294 | ifd = .false.
|
|---|
| 295 | if (md.ne.mx) ifd=.true.
|
|---|
| 296 | ifd=.true.
|
|---|
| 297 | c
|
|---|
| 298 | nrstd = md**2
|
|---|
| 299 | call lim_chk(nrstd,ldd,'urus ','ldd ','convop_fst')
|
|---|
| 300 | c
|
|---|
| 301 | do e=1,nel
|
|---|
| 302 | call grad_rstd(ur,us,ut,u(1,e),mx,md,if3d,ud)
|
|---|
| 303 | if (ifd) then ! dealiased
|
|---|
| 304 | do i=1,nrstd
|
|---|
| 305 | c C has the mass matrix factored in per (4.8.5), p. 227, DFM.
|
|---|
| 306 | ud(i) = c(i,e,1)*ur(i)+c(i,e,2)*us(i)
|
|---|
| 307 | enddo
|
|---|
| 308 | call intp_rstd(du(1,e),ud,mx,md,if3d,1) ! 1 --> transpose
|
|---|
| 309 | else
|
|---|
| 310 | do i=1,nrstd
|
|---|
| 311 | c C has the mass matrix factored in per (4.8.5), p. 227, DFM.
|
|---|
| 312 | du(i,e) = c(i,e,1)*ur(i)+c(i,e,2)*us(i)
|
|---|
| 313 | enddo
|
|---|
| 314 | endif
|
|---|
| 315 | enddo
|
|---|
| 316 | c
|
|---|
| 317 | return
|
|---|
| 318 | end
|
|---|
| 319 | c-----------------------------------------------------------------------
|
|---|
| 320 | subroutine grad_rstd(ur,us,ut,u,mx,md,if3d,ju) ! GLL->GL grad
|
|---|
| 321 |
|
|---|
| 322 | include 'SIZE'
|
|---|
| 323 | include 'DXYZ'
|
|---|
| 324 |
|
|---|
| 325 | real ur(1),us(1),ut(1),u(1),ju(1)
|
|---|
| 326 | logical if3d
|
|---|
| 327 |
|
|---|
| 328 | parameter (ldg=lxd**3,lwkd=4*lxd*lxd)
|
|---|
| 329 | common /dgrad/ d(ldg),dt(ldg),dg(ldg),dgt(ldg),jgl(ldg),jgt(ldg)
|
|---|
| 330 | $ , wkd(lwkd)
|
|---|
| 331 | real jgl,jgt
|
|---|
| 332 |
|
|---|
| 333 | call intp_rstd(ju,u,mx,md,if3d,0) ! 0 = forward
|
|---|
| 334 |
|
|---|
| 335 | m0 = md-1
|
|---|
| 336 | call get_dgl_ptr (ip,md,md)
|
|---|
| 337 | if (if3d) then
|
|---|
| 338 | call local_grad3(ur,us,ut,ju,m0,1,dg(ip),dgt(ip))
|
|---|
| 339 | else
|
|---|
| 340 | call local_grad2(ur,us ,ju,m0,1,dg(ip),dgt(ip))
|
|---|
| 341 | endif
|
|---|
| 342 |
|
|---|
| 343 | return
|
|---|
| 344 | end
|
|---|
| 345 | c-----------------------------------------------------------------------
|
|---|
| 346 | subroutine intp_rstd(ju,u,mx,md,if3d,idir) ! GLL->GL interpolation
|
|---|
| 347 |
|
|---|
| 348 | c GLL interpolation from mx to md.
|
|---|
| 349 |
|
|---|
| 350 | c If idir ^= 0, then apply transpose operator (md to mx)
|
|---|
| 351 |
|
|---|
| 352 | include 'SIZE'
|
|---|
| 353 |
|
|---|
| 354 | real ju(1),u(1)
|
|---|
| 355 | logical if3d
|
|---|
| 356 |
|
|---|
| 357 | parameter (ldg=lxd**3,lwkd=4*lxd*lxd)
|
|---|
| 358 | common /dgrad/ d(ldg),dt(ldg),dg(ldg),dgt(ldg),jgl(ldg),jgt(ldg)
|
|---|
| 359 | $ , wkd(lwkd)
|
|---|
| 360 | real jgl,jgt
|
|---|
| 361 |
|
|---|
| 362 | parameter (ld=2*lxd)
|
|---|
| 363 | common /ctmp0/ w(ld**ldim,2)
|
|---|
| 364 |
|
|---|
| 365 | call lim_chk(md,ld,'md ','ld ','grad_rstd ')
|
|---|
| 366 | call lim_chk(mx,ld,'mx ','ld ','grad_rstd ')
|
|---|
| 367 |
|
|---|
| 368 | ldw = 2*(ld**ldim)
|
|---|
| 369 |
|
|---|
| 370 | call get_int_ptr (i,mx,md)
|
|---|
| 371 | c
|
|---|
| 372 | if (idir.eq.0) then
|
|---|
| 373 | call specmpn(ju,md,u,mx,jgl(i),jgt(i),if3d,w,ldw)
|
|---|
| 374 | else
|
|---|
| 375 | call specmpn(ju,mx,u,md,jgt(i),jgl(i),if3d,w,ldw)
|
|---|
| 376 | endif
|
|---|
| 377 | c
|
|---|
| 378 | return
|
|---|
| 379 | end
|
|---|
| 380 | c-----------------------------------------------------------------------
|
|---|
| 381 | subroutine gen_int(jgl,jgt,mp,np,w)
|
|---|
| 382 | c
|
|---|
| 383 | c Generate interpolation from np GLL points to mp GL points
|
|---|
| 384 | c
|
|---|
| 385 | c jgl = interpolation matrix, mapping from velocity nodes to pressure
|
|---|
| 386 | c jgt = transpose of interpolation matrix
|
|---|
| 387 | c w = work array of size (np+mp)
|
|---|
| 388 | c
|
|---|
| 389 | c np = number of points on GLL grid
|
|---|
| 390 | c mp = number of points on GL grid
|
|---|
| 391 | c
|
|---|
| 392 | c
|
|---|
| 393 | real jgl(mp,np),jgt(np*mp),w(1)
|
|---|
| 394 | c
|
|---|
| 395 | iz = 1
|
|---|
| 396 | id = iz + np
|
|---|
| 397 | c
|
|---|
| 398 | call zwgll (w(iz),jgt,np)
|
|---|
| 399 | call zwgl (w(id),jgt,mp)
|
|---|
| 400 | c
|
|---|
| 401 | n = np-1
|
|---|
| 402 | do i=1,mp
|
|---|
| 403 | call fd_weights_full(w(id+i-1),w(iz),n,0,jgt)
|
|---|
| 404 | do j=1,np
|
|---|
| 405 | jgl(i,j) = jgt(j) ! Interpolation matrix
|
|---|
| 406 | enddo
|
|---|
| 407 | enddo
|
|---|
| 408 | c
|
|---|
| 409 | call transpose(jgt,np,jgl,mp)
|
|---|
| 410 | c
|
|---|
| 411 | return
|
|---|
| 412 | end
|
|---|
| 413 | c-----------------------------------------------------------------------
|
|---|
| 414 | subroutine gen_dgl(dgl,dgt,mp,np,w)
|
|---|
| 415 | c
|
|---|
| 416 | c Generate derivative from np GL points onto mp GL points
|
|---|
| 417 | c
|
|---|
| 418 | c dgl = interpolation matrix, mapping from velocity nodes to pressure
|
|---|
| 419 | c dgt = transpose of interpolation matrix
|
|---|
| 420 | c w = work array of size (3*np+mp)
|
|---|
| 421 | c
|
|---|
| 422 | c np = number of points on GLL grid
|
|---|
| 423 | c mp = number of points on GL grid
|
|---|
| 424 | c
|
|---|
| 425 | c
|
|---|
| 426 | c
|
|---|
| 427 | real dgl(mp,np),dgt(np*mp),w(1)
|
|---|
| 428 | c
|
|---|
| 429 | c
|
|---|
| 430 | iz = 1
|
|---|
| 431 | id = iz + np
|
|---|
| 432 | c
|
|---|
| 433 | call zwgl (w(iz),dgt,np) ! GL points
|
|---|
| 434 | call zwgl (w(id),dgt,mp) ! GL points
|
|---|
| 435 | c
|
|---|
| 436 | ndgt = 2*np
|
|---|
| 437 | ldgt = mp*np
|
|---|
| 438 | call lim_chk(ndgt,ldgt,'ldgt ','dgt ','gen_dgl ')
|
|---|
| 439 | c
|
|---|
| 440 | n = np-1
|
|---|
| 441 | do i=1,mp
|
|---|
| 442 | call fd_weights_full(w(id+i-1),w(iz),n,1,dgt) ! 1=1st deriv.
|
|---|
| 443 | do j=1,np
|
|---|
| 444 | dgl(i,j) = dgt(np+j) ! Derivative matrix
|
|---|
| 445 | enddo
|
|---|
| 446 | enddo
|
|---|
| 447 | c
|
|---|
| 448 | call transpose(dgt,np,dgl,mp)
|
|---|
| 449 | c
|
|---|
| 450 | return
|
|---|
| 451 | end
|
|---|
| 452 | c-----------------------------------------------------------------------
|
|---|
| 453 | subroutine lim_chk(n,m,avar5,lvar5,sub_name10)
|
|---|
| 454 | include 'SIZE' ! need nid
|
|---|
| 455 | character*5 avar5,lvar5
|
|---|
| 456 | character*10 sub_name10
|
|---|
| 457 |
|
|---|
| 458 | if (n.gt.m) then
|
|---|
| 459 | write(6,1) nid,n,m,avar5,lvar5,sub_name10
|
|---|
| 460 | 1 format(i8,' ERROR: :',2i12,2(1x,a5),1x,a10)
|
|---|
| 461 | call exitti('lim_chk problem. $',n)
|
|---|
| 462 | endif
|
|---|
| 463 |
|
|---|
| 464 | return
|
|---|
| 465 | end
|
|---|
| 466 | c-----------------------------------------------------------------------
|
|---|
| 467 | subroutine get_int_ptr (ip,mx,md) ! GLL-->GL pointer
|
|---|
| 468 |
|
|---|
| 469 | c Get pointer to jgl() for interpolation pair (mx,md)
|
|---|
| 470 |
|
|---|
| 471 | include 'SIZE'
|
|---|
| 472 |
|
|---|
| 473 | parameter (ldg=lxd**3,lwkd=4*lxd*lxd)
|
|---|
| 474 | common /dgrad/ d(ldg),dt(ldg),dg(ldg),dgt(ldg),jgl(ldg),jgt(ldg)
|
|---|
| 475 | $ , wkd(lwkd)
|
|---|
| 476 | real jgl,jgt
|
|---|
| 477 | c
|
|---|
| 478 | parameter (ld=2*lxd)
|
|---|
| 479 | common /igrad/ pd (0:ld*ld)
|
|---|
| 480 | $ , pdg (0:ld*ld)
|
|---|
| 481 | $ , pjgl (0:ld*ld)
|
|---|
| 482 | integer pd , pdg , pjgl
|
|---|
| 483 | c
|
|---|
| 484 | ij = md + ld*(mx-1)
|
|---|
| 485 | ip = pjgl(ij)
|
|---|
| 486 | c
|
|---|
| 487 | if (ip.eq.0) then
|
|---|
| 488 | c
|
|---|
| 489 | nstore = pjgl(0)
|
|---|
| 490 | pjgl(ij) = nstore+1
|
|---|
| 491 | nstore = nstore + md*mx
|
|---|
| 492 | pjgl(0) = nstore
|
|---|
| 493 | ip = pjgl(ij)
|
|---|
| 494 | c
|
|---|
| 495 | nwrkd = mx + md
|
|---|
| 496 | call lim_chk(nstore,ldg ,'jgl ','ldg ','get_int_pt')
|
|---|
| 497 | call lim_chk(nwrkd ,lwkd,'wkd ','lwkd ','get_int_pt')
|
|---|
| 498 | c
|
|---|
| 499 | call gen_int(jgl(ip),jgt(ip),md,mx,wkd)
|
|---|
| 500 | endif
|
|---|
| 501 | c
|
|---|
| 502 | return
|
|---|
| 503 | end
|
|---|
| 504 | c-----------------------------------------------------------------------
|
|---|
| 505 | subroutine get_dgl_ptr (ip,mx,md)
|
|---|
| 506 | c
|
|---|
| 507 | c Get pointer to GL-GL interpolation dgl() for pair (mx,md)
|
|---|
| 508 | c
|
|---|
| 509 | include 'SIZE'
|
|---|
| 510 | c
|
|---|
| 511 | parameter (ldg=lxd**3,lwkd=4*lxd*lxd)
|
|---|
| 512 | common /dgrad/ d(ldg),dt(ldg),dg(ldg),dgt(ldg),jgl(ldg),jgt(ldg)
|
|---|
| 513 | $ , wkd(lwkd)
|
|---|
| 514 | real jgl,jgt
|
|---|
| 515 | c
|
|---|
| 516 | parameter (ld=2*lxd)
|
|---|
| 517 | common /jgrad/ pd (0:ld*ld)
|
|---|
| 518 | $ , pdg (0:ld*ld)
|
|---|
| 519 | $ , pjgl (0:ld*ld)
|
|---|
| 520 | integer pd , pdg , pjgl
|
|---|
| 521 | c
|
|---|
| 522 | ij = md + ld*(mx-1)
|
|---|
| 523 | ip = pdg (ij)
|
|---|
| 524 |
|
|---|
| 525 | if (ip.eq.0) then
|
|---|
| 526 |
|
|---|
| 527 | nstore = pdg (0)
|
|---|
| 528 | pdg (ij) = nstore+1
|
|---|
| 529 | nstore = nstore + md*mx
|
|---|
| 530 | pdg (0) = nstore
|
|---|
| 531 | ip = pdg (ij)
|
|---|
| 532 | c
|
|---|
| 533 | nwrkd = mx + md
|
|---|
| 534 | call lim_chk(nstore,ldg ,'dg ','ldg ','get_dgl_pt')
|
|---|
| 535 | call lim_chk(nwrkd ,lwkd,'wkd ','lwkd ','get_dgl_pt')
|
|---|
| 536 | c
|
|---|
| 537 | call gen_dgl(dg (ip),dgt(ip),md,mx,wkd)
|
|---|
| 538 | endif
|
|---|
| 539 | c
|
|---|
| 540 | return
|
|---|
| 541 | end
|
|---|
| 542 | c-----------------------------------------------------------------------
|
|---|
| 543 | subroutine set_conv_char(ct,c,ux,uy,uz,nelc,tau,ifnew)
|
|---|
| 544 | include 'SIZE'
|
|---|
| 545 | include 'TSTEP'
|
|---|
| 546 |
|
|---|
| 547 | real ct(0:1) ! time stamps for saved field (0=#flds)
|
|---|
| 548 | real c(1) ! saved vel. fields, dealiased etc.
|
|---|
| 549 | real ux(1),uy(1),uz(1) ! input vel. field
|
|---|
| 550 | integer nelc ! number of elements in conv. field
|
|---|
| 551 | logical ifnew ! =true if shifting stack of fields
|
|---|
| 552 |
|
|---|
| 553 | numr = lxd*lyd*lzd*lelv*ldim*(lorder+1)
|
|---|
| 554 | denr = lxd*lyd*lzd*nelv*ldim
|
|---|
| 555 | nconv_max = numr/denr
|
|---|
| 556 | if (nconv_max.lt.nbdinp+1)
|
|---|
| 557 | $ call exitti(
|
|---|
| 558 | $ 'ABORT: not enough memory for characteristics scheme!$',
|
|---|
| 559 | $ nconv_max)
|
|---|
| 560 |
|
|---|
| 561 | nc = ct(0)
|
|---|
| 562 |
|
|---|
| 563 | m = lxd*lyd*lzd*nelc*ldim
|
|---|
| 564 |
|
|---|
| 565 | call set_ct_cvx
|
|---|
| 566 | $ (ct,c,m,ux,uy,uz,tau,nc,nconv_max,nelc,ifnew)
|
|---|
| 567 |
|
|---|
| 568 | nc = min (nc,nbdinp)
|
|---|
| 569 | ct(0) = nc ! store current count
|
|---|
| 570 |
|
|---|
| 571 | return
|
|---|
| 572 | end
|
|---|
| 573 | c-----------------------------------------------------------------------
|
|---|
| 574 | subroutine set_ct_cvx(ct,c,m,u,v,w,tau,nc,mc,nelc,ifnew)
|
|---|
| 575 | include 'SIZE'
|
|---|
| 576 | include 'INPUT' ! ifcons
|
|---|
| 577 |
|
|---|
| 578 | real ct(0:1),c(m,1)
|
|---|
| 579 | real u(1),v(1),w(1)
|
|---|
| 580 | logical ifnew
|
|---|
| 581 |
|
|---|
| 582 | if (ifnew) then
|
|---|
| 583 |
|
|---|
| 584 | c Shift existing convecting fields
|
|---|
| 585 | c Note: "1" entry is most recent
|
|---|
| 586 |
|
|---|
| 587 | nc = nc+1
|
|---|
| 588 | nc = min(nc,mc)
|
|---|
| 589 | ct(0) = nc
|
|---|
| 590 |
|
|---|
| 591 | do i=nc,2,-1
|
|---|
| 592 | call copy(c(1,i),c(1,i-1),m)
|
|---|
| 593 | ct(i) = ct(i-1)
|
|---|
| 594 | enddo
|
|---|
| 595 | endif
|
|---|
| 596 |
|
|---|
| 597 | c Save time and map the current velocity to rst coordinates.
|
|---|
| 598 |
|
|---|
| 599 | ix = 1
|
|---|
| 600 | iy = ix + lxd*lyd*lzd*nelc
|
|---|
| 601 | iz = iy + lxd*lyd*lzd*nelc
|
|---|
| 602 |
|
|---|
| 603 | if (ifcons) then
|
|---|
| 604 | call set_convect_cons(c(ix,1),c(iy,1),c(iz,1),u,v,w)
|
|---|
| 605 | else
|
|---|
| 606 | call set_convect_new (c(ix,1),c(iy,1),c(iz,1),u,v,w)
|
|---|
| 607 | endif
|
|---|
| 608 |
|
|---|
| 609 | ct(1) = tau
|
|---|
| 610 |
|
|---|
| 611 | return
|
|---|
| 612 | end
|
|---|
| 613 | c-----------------------------------------------------------------------
|
|---|
| 614 | subroutine grad_rst(ur,us,ut,u,md,if3d) ! Gauss-->Gauss grad
|
|---|
| 615 |
|
|---|
| 616 | include 'SIZE'
|
|---|
| 617 | include 'DXYZ'
|
|---|
| 618 |
|
|---|
| 619 | real ur(1),us(1),ut(1),u(1)
|
|---|
| 620 | logical if3d
|
|---|
| 621 |
|
|---|
| 622 | parameter (ldg=lxd**3,lwkd=4*lxd*lxd)
|
|---|
| 623 | common /dgrad/ d(ldg),dt(ldg),dg(ldg),dgt(ldg),jgl(ldg),jgt(ldg)
|
|---|
| 624 | $ , wkd(lwkd)
|
|---|
| 625 | real jgl,jgt
|
|---|
| 626 |
|
|---|
| 627 | m0 = md-1
|
|---|
| 628 | call get_dgl_ptr (ip,md,md)
|
|---|
| 629 | if (if3d) then
|
|---|
| 630 | call local_grad3(ur,us,ut,u,m0,1,dg(ip),dgt(ip))
|
|---|
| 631 | else
|
|---|
| 632 | call local_grad2(ur,us ,u,m0,1,dg(ip),dgt(ip))
|
|---|
| 633 | endif
|
|---|
| 634 |
|
|---|
| 635 | return
|
|---|
| 636 | end
|
|---|
| 637 | c-----------------------------------------------------------------------
|
|---|
| 638 | subroutine convect_new(bdu,u,ifuf,cx,cy,cz,ifcf)
|
|---|
| 639 |
|
|---|
| 640 | C Compute dealiased form: J^T Bf *JC .grad Ju w/ correct Jacobians
|
|---|
| 641 | C
|
|---|
| 642 | include 'SIZE'
|
|---|
| 643 | include 'TOTAL'
|
|---|
| 644 |
|
|---|
| 645 | real bdu(1),u(1),cx(1),cy(1),cz(1)
|
|---|
| 646 | logical ifuf,ifcf ! u and/or c already on fine mesh?
|
|---|
| 647 |
|
|---|
| 648 | parameter (lxy=lx1*ly1*lz1,ltd=lxd*lyd*lzd)
|
|---|
| 649 | common /scrcv/ fx(ltd),fy(ltd),fz(ltd)
|
|---|
| 650 | $ , ur(ltd),us(ltd),ut(ltd)
|
|---|
| 651 | $ , tr(ltd,3),uf(ltd)
|
|---|
| 652 |
|
|---|
| 653 | integer e
|
|---|
| 654 |
|
|---|
| 655 | call set_dealias_rx
|
|---|
| 656 |
|
|---|
| 657 | nxyz1 = lx1*ly1*lz1
|
|---|
| 658 | nxyzd = lxd*lyd*lzd
|
|---|
| 659 |
|
|---|
| 660 | nxyzu = nxyz1
|
|---|
| 661 | if (ifuf) nxyzu = nxyzd
|
|---|
| 662 |
|
|---|
| 663 | nxyzc = nxyz1
|
|---|
| 664 | if (ifcf) nxyzc = nxyzd
|
|---|
| 665 |
|
|---|
| 666 | iu = 1 ! pointer to scalar field u
|
|---|
| 667 | ic = 1 ! pointer to vector field C
|
|---|
| 668 | ib = 1 ! pointer to scalar field Bdu
|
|---|
| 669 |
|
|---|
| 670 |
|
|---|
| 671 | do e=1,nelv
|
|---|
| 672 |
|
|---|
| 673 | if (ifcf) then
|
|---|
| 674 |
|
|---|
| 675 | call copy(tr(1,1),cx(ic),nxyzd) ! already in rst form
|
|---|
| 676 | call copy(tr(1,2),cy(ic),nxyzd)
|
|---|
| 677 | if (if3d) call copy(tr(1,3),cz(ic),nxyzd)
|
|---|
| 678 |
|
|---|
| 679 | else ! map coarse velocity to fine mesh (C-->F)
|
|---|
| 680 |
|
|---|
| 681 | call intp_rstd(fx,cx(ic),lx1,lxd,if3d,0) ! 0 --> forward
|
|---|
| 682 | call intp_rstd(fy,cy(ic),lx1,lxd,if3d,0) ! 0 --> forward
|
|---|
| 683 | if (if3d) call intp_rstd(fz,cz(ic),lx1,lxd,if3d,0) ! 0 --> forward
|
|---|
| 684 |
|
|---|
| 685 | if (if3d) then ! Convert convector F to r-s-t coordinates
|
|---|
| 686 |
|
|---|
| 687 | do i=1,nxyzd
|
|---|
| 688 | tr(i,1)=rx(i,1,e)*fx(i)+rx(i,2,e)*fy(i)+rx(i,3,e)*fz(i)
|
|---|
| 689 | tr(i,2)=rx(i,4,e)*fx(i)+rx(i,5,e)*fy(i)+rx(i,6,e)*fz(i)
|
|---|
| 690 | tr(i,3)=rx(i,7,e)*fx(i)+rx(i,8,e)*fy(i)+rx(i,9,e)*fz(i)
|
|---|
| 691 | enddo
|
|---|
| 692 |
|
|---|
| 693 | else
|
|---|
| 694 |
|
|---|
| 695 | do i=1,nxyzd
|
|---|
| 696 | tr(i,1)=rx(i,1,e)*fx(i)+rx(i,2,e)*fy(i)
|
|---|
| 697 | tr(i,2)=rx(i,3,e)*fx(i)+rx(i,4,e)*fy(i)
|
|---|
| 698 | enddo
|
|---|
| 699 |
|
|---|
| 700 | endif
|
|---|
| 701 |
|
|---|
| 702 | endif
|
|---|
| 703 |
|
|---|
| 704 | if (ifuf) then
|
|---|
| 705 | call grad_rst(ur,us,ut,u(iu),lxd,if3d)
|
|---|
| 706 | else
|
|---|
| 707 | call intp_rstd(uf,u(iu),lx1,lxd,if3d,0) ! 0 --> forward
|
|---|
| 708 | call grad_rst(ur,us,ut,uf,lxd,if3d)
|
|---|
| 709 | endif
|
|---|
| 710 |
|
|---|
| 711 | if (if3d) then
|
|---|
| 712 | do i=1,nxyzd ! mass matrix included, per DFM (4.8.5)
|
|---|
| 713 | uf(i) = tr(i,1)*ur(i)+tr(i,2)*us(i)+tr(i,3)*ut(i)
|
|---|
| 714 | enddo
|
|---|
| 715 | else
|
|---|
| 716 | do i=1,nxyzd ! mass matrix included, per DFM (4.8.5)
|
|---|
| 717 | uf(i) = tr(i,1)*ur(i)+tr(i,2)*us(i)
|
|---|
| 718 | enddo
|
|---|
| 719 | endif
|
|---|
| 720 | call intp_rstd(bdu(ib),uf,lx1,lxd,if3d,1) ! Project back to coarse
|
|---|
| 721 |
|
|---|
| 722 | ic = ic + nxyzc
|
|---|
| 723 | iu = iu + nxyzu
|
|---|
| 724 | ib = ib + nxyz1
|
|---|
| 725 |
|
|---|
| 726 | enddo
|
|---|
| 727 |
|
|---|
| 728 | return
|
|---|
| 729 | end
|
|---|
| 730 | c-----------------------------------------------------------------------
|
|---|
| 731 | subroutine convect_cons(bdu,u,ifuf,cx,cy,cz,ifcf)
|
|---|
| 732 |
|
|---|
| 733 | c Compute dealiased form: J^T Bf *div. JC Ju w/ correct Jacobians
|
|---|
| 734 |
|
|---|
| 735 | c conservative form
|
|---|
| 736 |
|
|---|
| 737 |
|
|---|
| 738 | include 'SIZE'
|
|---|
| 739 | include 'TOTAL'
|
|---|
| 740 |
|
|---|
| 741 | real bdu(1),u(1),cx(1),cy(1),cz(1)
|
|---|
| 742 |
|
|---|
| 743 | logical ifuf,ifcf ! u and/or c already on fine mesh?
|
|---|
| 744 |
|
|---|
| 745 | parameter (lxy=lx1*ly1*lz1,ltd=lxd*lyd*lzd)
|
|---|
| 746 | common /scrcv/ uf(ltd),cf(ltd),cu(ltd)
|
|---|
| 747 | $ , cr(ltd),cs(ltd),ct(ltd)
|
|---|
| 748 |
|
|---|
| 749 |
|
|---|
| 750 | integer e
|
|---|
| 751 |
|
|---|
| 752 | call set_dealias_rx
|
|---|
| 753 |
|
|---|
| 754 | nxyz1 = lx1*ly1*lz1
|
|---|
| 755 | nxyzd = lxd*lyd*lzd
|
|---|
| 756 |
|
|---|
| 757 | nxyzu = nxyz1
|
|---|
| 758 | if (ifuf) nxyzu = nxyzd
|
|---|
| 759 |
|
|---|
| 760 | nxyzc = nxyz1
|
|---|
| 761 | if (ifcf) nxyzc = nxyzd
|
|---|
| 762 |
|
|---|
| 763 | iu = 1 ! pointer to scalar field u
|
|---|
| 764 | ic = 1 ! pointer to vector field C
|
|---|
| 765 | ib = 1 ! pointer to scalar field Bdu
|
|---|
| 766 |
|
|---|
| 767 | do e=1,nelv
|
|---|
| 768 |
|
|---|
| 769 | call intp_rstd(uf,u(iu),lx1,lxd,if3d,0) ! 0 --> forward
|
|---|
| 770 |
|
|---|
| 771 | call rzero(cu,nxyzd)
|
|---|
| 772 | do i=1,ldim
|
|---|
| 773 |
|
|---|
| 774 | if (ifcf) then ! C is already on fine mesh
|
|---|
| 775 |
|
|---|
| 776 | call exitt ! exit for now
|
|---|
| 777 |
|
|---|
| 778 | else ! map coarse velocity to fine mesh (C-->F)
|
|---|
| 779 |
|
|---|
| 780 | if (i.eq.1) call intp_rstd(cf,cx(ic),lx1,lxd,if3d,0) ! 0 --> forward
|
|---|
| 781 | if (i.eq.2) call intp_rstd(cf,cy(ic),lx1,lxd,if3d,0) ! 0 --> forward
|
|---|
| 782 | if (i.eq.3) call intp_rstd(cf,cz(ic),lx1,lxd,if3d,0) ! 0 --> forward
|
|---|
| 783 |
|
|---|
| 784 | call col2(cf,uf,nxyzd) ! collocate C and u on fine mesh
|
|---|
| 785 |
|
|---|
| 786 | call grad_rst(cr,cs,ct,cf,lxd,if3d) ! d/dr (C_i*u)
|
|---|
| 787 |
|
|---|
| 788 | if (if3d) then
|
|---|
| 789 |
|
|---|
| 790 | do j=1,nxyzd
|
|---|
| 791 | cu(j)=cu(j)
|
|---|
| 792 | $ +cr(j)*rx(j,i,e)+cs(j)*rx(j,i+3,e)+ct(j)*rx(j,i+6,e)
|
|---|
| 793 | enddo
|
|---|
| 794 |
|
|---|
| 795 | else ! 2D
|
|---|
| 796 |
|
|---|
| 797 | do j=1,nxyzd
|
|---|
| 798 | cu(j)=cu(j)
|
|---|
| 799 | $ +cr(j)*rx(j,i,e)+cs(j)*rx(j,i+2,e)
|
|---|
| 800 | enddo
|
|---|
| 801 |
|
|---|
| 802 | endif
|
|---|
| 803 | endif
|
|---|
| 804 | enddo
|
|---|
| 805 |
|
|---|
| 806 | call intp_rstd(bdu(ib),cu,lx1,lxd,if3d,1) ! Project back to coarse
|
|---|
| 807 |
|
|---|
| 808 | ic = ic + nxyzc
|
|---|
| 809 | iu = iu + nxyzu
|
|---|
| 810 | ib = ib + nxyz1
|
|---|
| 811 |
|
|---|
| 812 | enddo
|
|---|
| 813 |
|
|---|
| 814 | return
|
|---|
| 815 | end
|
|---|
| 816 | c-----------------------------------------------------------------------
|
|---|
| 817 | subroutine set_convect_cons(cx,cy,cz,ux,uy,uz)
|
|---|
| 818 |
|
|---|
| 819 | c Put vx,vy,vz on fine mesh (for conservation form)
|
|---|
| 820 |
|
|---|
| 821 |
|
|---|
| 822 | include 'SIZE'
|
|---|
| 823 | include 'TOTAL'
|
|---|
| 824 |
|
|---|
| 825 | parameter (lxy=lx1*ly1*lz1,ltd=lxd*lyd*lzd)
|
|---|
| 826 |
|
|---|
| 827 | real cx(ltd,1),cy(ltd,1),cz(ltd,1)
|
|---|
| 828 | real ux(lxy,1),uy(lxy,1),uz(lxy,1)
|
|---|
| 829 |
|
|---|
| 830 | integer e
|
|---|
| 831 |
|
|---|
| 832 | call set_dealias_rx
|
|---|
| 833 |
|
|---|
| 834 | do e=1,nelv ! Map coarse velocity to fine mesh (C-->F)
|
|---|
| 835 |
|
|---|
| 836 | call intp_rstd(cx(1,e),ux(1,e),lx1,lxd,if3d,0) ! 0 --> forward
|
|---|
| 837 | call intp_rstd(cy(1,e),uy(1,e),lx1,lxd,if3d,0) ! 0 --> forward
|
|---|
| 838 | if (if3d) call intp_rstd(cz(1,e),uz(1,e),lx1,lxd,if3d,0) ! 0 --> forward
|
|---|
| 839 |
|
|---|
| 840 | enddo
|
|---|
| 841 |
|
|---|
| 842 | return
|
|---|
| 843 | end
|
|---|
| 844 | c-----------------------------------------------------------------------
|
|---|
| 845 | subroutine set_convect_new(cr,cs,ct,ux,uy,uz)
|
|---|
| 846 | C
|
|---|
| 847 | C Put vxd,vyd,vzd into rst form on fine mesh
|
|---|
| 848 | C
|
|---|
| 849 | C For rst form, see eq. (4.8.5) in Deville, Fischer, Mund (2002).
|
|---|
| 850 | C
|
|---|
| 851 | include 'SIZE'
|
|---|
| 852 | include 'TOTAL'
|
|---|
| 853 |
|
|---|
| 854 | parameter (lxy=lx1*ly1*lz1,ltd=lxd*lyd*lzd)
|
|---|
| 855 |
|
|---|
| 856 | real cr(ltd,1),cs(ltd,1),ct(ltd,1)
|
|---|
| 857 | real ux(lxy,1),uy(lxy,1),uz(lxy,1)
|
|---|
| 858 |
|
|---|
| 859 | common /scrcv/ fx(ltd),fy(ltd),fz(ltd)
|
|---|
| 860 | $ , ur(ltd),us(ltd),ut(ltd)
|
|---|
| 861 | $ , tr(ltd,3),uf(ltd)
|
|---|
| 862 |
|
|---|
| 863 | integer e
|
|---|
| 864 |
|
|---|
| 865 | call set_dealias_rx
|
|---|
| 866 |
|
|---|
| 867 | nxyz1 = lx1*ly1*lz1
|
|---|
| 868 | nxyzd = lxd*lyd*lzd
|
|---|
| 869 |
|
|---|
| 870 | ic = 1 ! pointer to vector field C
|
|---|
| 871 |
|
|---|
| 872 | do e=1,nelv
|
|---|
| 873 |
|
|---|
| 874 | c Map coarse velocity to fine mesh (C-->F)
|
|---|
| 875 |
|
|---|
| 876 | call intp_rstd(fx,ux(1,e),lx1,lxd,if3d,0) ! 0 --> forward
|
|---|
| 877 | call intp_rstd(fy,uy(1,e),lx1,lxd,if3d,0) ! 0 --> forward
|
|---|
| 878 | if (if3d) call intp_rstd(fz,uz(1,e),lx1,lxd,if3d,0) ! 0 --> forward
|
|---|
| 879 |
|
|---|
| 880 | c Convert convector F to r-s-t coordinates
|
|---|
| 881 |
|
|---|
| 882 | if (if3d) then
|
|---|
| 883 |
|
|---|
| 884 | do i=1,nxyzd
|
|---|
| 885 | cr(i,e)=rx(i,1,e)*fx(i)+rx(i,2,e)*fy(i)+rx(i,3,e)*fz(i)
|
|---|
| 886 | cs(i,e)=rx(i,4,e)*fx(i)+rx(i,5,e)*fy(i)+rx(i,6,e)*fz(i)
|
|---|
| 887 | ct(i,e)=rx(i,7,e)*fx(i)+rx(i,8,e)*fy(i)+rx(i,9,e)*fz(i)
|
|---|
| 888 | enddo
|
|---|
| 889 |
|
|---|
| 890 | else
|
|---|
| 891 |
|
|---|
| 892 | do i=1,nxyzd
|
|---|
| 893 | cr(i,e)=rx(i,1,e)*fx(i)+rx(i,2,e)*fy(i)
|
|---|
| 894 | cs(i,e)=rx(i,3,e)*fx(i)+rx(i,4,e)*fy(i)
|
|---|
| 895 | enddo
|
|---|
| 896 |
|
|---|
| 897 | endif
|
|---|
| 898 | enddo
|
|---|
| 899 |
|
|---|
| 900 | return
|
|---|
| 901 | end
|
|---|
| 902 | c-----------------------------------------------------------------------
|
|---|
| 903 | subroutine advchar
|
|---|
| 904 | c
|
|---|
| 905 | c Compute convective contribution using
|
|---|
| 906 | c operator-integrator-factor method (characteristics).
|
|---|
| 907 | c
|
|---|
| 908 | include 'SIZE'
|
|---|
| 909 | include 'MASS'
|
|---|
| 910 | include 'INPUT'
|
|---|
| 911 | include 'SOLN'
|
|---|
| 912 | include 'TSTEP'
|
|---|
| 913 | include 'PARALLEL'
|
|---|
| 914 |
|
|---|
| 915 | include 'CTIMER'
|
|---|
| 916 |
|
|---|
| 917 | common /cchar/ ct_vx(0:lorder) ! time for each slice in c_vx()
|
|---|
| 918 |
|
|---|
| 919 | common /scruz/ phx (lx1*ly1*lz1*lelt)
|
|---|
| 920 | $ , phy (lx1*ly1*lz1*lelt)
|
|---|
| 921 | $ , phz (lx1*ly1*lz1*lelt)
|
|---|
| 922 | $ , hmsk (lx1*ly1*lz1*lelt)
|
|---|
| 923 |
|
|---|
| 924 | if (icalld.eq.0) tadvc=0.0
|
|---|
| 925 | icalld=icalld+1
|
|---|
| 926 | nadvc=icalld
|
|---|
| 927 | etime1=dnekclock()
|
|---|
| 928 |
|
|---|
| 929 | dti = 1./dt
|
|---|
| 930 | n = lx1*ly1*lz1*nelv
|
|---|
| 931 |
|
|---|
| 932 | call char_conv(phx,vx,vxlag,bm1,bm1lag,hmsk,c_vx,ct_vx,gsh_fld(1))
|
|---|
| 933 | call char_conv(phy,vy,vylag,bm1,bm1lag,hmsk,c_vx,ct_vx,gsh_fld(1))
|
|---|
| 934 | if (if3d) call char_conv
|
|---|
| 935 | $ (phz,vz,vzlag,bm1,bm1lag,hmsk,c_vx,ct_vx,gsh_fld(1))
|
|---|
| 936 |
|
|---|
| 937 | call cfill(hmsk,dti,n)
|
|---|
| 938 | if(.not. iflomach) call col2(hmsk,vtrans,n)
|
|---|
| 939 |
|
|---|
| 940 | if (if3d) then
|
|---|
| 941 |
|
|---|
| 942 | do i=1,n
|
|---|
| 943 | h2i = hmsk(i)
|
|---|
| 944 | bfx(i,1,1,1) = bfx(i,1,1,1)+phx(i)*h2i
|
|---|
| 945 | bfy(i,1,1,1) = bfy(i,1,1,1)+phy(i)*h2i
|
|---|
| 946 | bfz(i,1,1,1) = bfz(i,1,1,1)+phz(i)*h2i
|
|---|
| 947 | enddo
|
|---|
| 948 |
|
|---|
| 949 | else
|
|---|
| 950 |
|
|---|
| 951 | do i=1,n
|
|---|
| 952 | h2i = hmsk(i)
|
|---|
| 953 | bfx(i,1,1,1) = bfx(i,1,1,1)+phx(i)*h2i
|
|---|
| 954 | bfy(i,1,1,1) = bfy(i,1,1,1)+phy(i)*h2i
|
|---|
| 955 | enddo
|
|---|
| 956 |
|
|---|
| 957 | endif
|
|---|
| 958 |
|
|---|
| 959 | tadvc=tadvc+(dnekclock()-etime1)
|
|---|
| 960 |
|
|---|
| 961 | return
|
|---|
| 962 | end
|
|---|
| 963 | c-----------------------------------------------------------------------
|
|---|
| 964 | subroutine convch
|
|---|
| 965 |
|
|---|
| 966 | c Compute convective contribution using
|
|---|
| 967 | c operator-integrator-factor method (characteristics).
|
|---|
| 968 |
|
|---|
| 969 | include 'SIZE'
|
|---|
| 970 | include 'MASS'
|
|---|
| 971 | include 'INPUT'
|
|---|
| 972 | include 'SOLN'
|
|---|
| 973 | include 'TSTEP'
|
|---|
| 974 | include 'PARALLEL'
|
|---|
| 975 | include 'CTIMER'
|
|---|
| 976 |
|
|---|
| 977 | common /cchar/ ct_vx(0:lorder) ! time for each slice in c_vx()
|
|---|
| 978 |
|
|---|
| 979 | common /scruz/ phi (lx1*ly1*lz1*lelt)
|
|---|
| 980 | $ , hmsk (lx1*ly1*lz1*lelt)
|
|---|
| 981 |
|
|---|
| 982 | if (icalld.eq.0) tadvc=0.0
|
|---|
| 983 | icalld=icalld+1
|
|---|
| 984 | nadvc=icalld
|
|---|
| 985 | etime1=dnekclock()
|
|---|
| 986 |
|
|---|
| 987 | n = lx1*ly1*lz1*nelv
|
|---|
| 988 | dti = 1./dt
|
|---|
| 989 |
|
|---|
| 990 | if(nid.eq.0 .and. loglevel.gt.2) write(6,*) 'convch', ifield
|
|---|
| 991 | call char_conv(phi,t(1,1,1,1,ifield-1),tlag(1,1,1,1,1,ifield-1)
|
|---|
| 992 | $ ,bm1,bm1lag,hmsk,c_vx,ct_vx,gsh_fld(1))
|
|---|
| 993 |
|
|---|
| 994 | do i=1,n
|
|---|
| 995 | bq(i,1,1,1,ifield-1) = bq(i,1,1,1,ifield-1)
|
|---|
| 996 | $ + phi(i)*vtrans(i,1,1,1,ifield)*dti
|
|---|
| 997 | enddo
|
|---|
| 998 |
|
|---|
| 999 | tadvc=tadvc+(dnekclock()-etime1)
|
|---|
| 1000 |
|
|---|
| 1001 | return
|
|---|
| 1002 | end
|
|---|
| 1003 | c-----------------------------------------------------------------------
|
|---|
| 1004 | subroutine convop_cons_3d(du,u,c,mx,md,nel) ! Conservation form
|
|---|
| 1005 |
|
|---|
| 1006 | c Apply convecting field c to scalar field u, conservation form d/dxj cj phi
|
|---|
| 1007 |
|
|---|
| 1008 | c Assumes that current convecting field is on dealias mesh, in c()
|
|---|
| 1009 |
|
|---|
| 1010 | include 'SIZE'
|
|---|
| 1011 | include 'DEALIAS'
|
|---|
| 1012 | include 'GEOM'
|
|---|
| 1013 |
|
|---|
| 1014 | real du(mx*mx*mx,nel)
|
|---|
| 1015 | real u(mx*mx*mx,nel)
|
|---|
| 1016 | real c(md*md*md,nel,3)
|
|---|
| 1017 | parameter (ldd=lxd*lyd*lzd)
|
|---|
| 1018 | common /ctmp1/ ur(ldd),us(ldd),ut(ldd),ju(ldd),ud(ldd),tu(ldd)
|
|---|
| 1019 | real ju
|
|---|
| 1020 |
|
|---|
| 1021 | logical if3d,ifd
|
|---|
| 1022 | integer e
|
|---|
| 1023 |
|
|---|
| 1024 | if3d = .true.
|
|---|
| 1025 | ifd = .false.
|
|---|
| 1026 | if (md.ne.mx) ifd=.true.
|
|---|
| 1027 | ifd=.true.
|
|---|
| 1028 |
|
|---|
| 1029 | nrstd = md**3
|
|---|
| 1030 | call lim_chk(nrstd,ldd,'urus ','ldd ','convp_cons')
|
|---|
| 1031 |
|
|---|
| 1032 | do e=1,nel
|
|---|
| 1033 |
|
|---|
| 1034 | call intp_rstd (ju,u(1,e),mx,md,if3d,0) ! 0 = forward; on Gauss points!
|
|---|
| 1035 | call rzero (ud,nrstd)
|
|---|
| 1036 |
|
|---|
| 1037 | do j=1,ldim
|
|---|
| 1038 | do i=1,nrstd
|
|---|
| 1039 | tu(i)=c(i,e,j)*ju(i) ! C_j*T
|
|---|
| 1040 | enddo
|
|---|
| 1041 | call grad_rst(ur,us,ut,tu,md,if3d) ! Already on fine (Gauss) mesh
|
|---|
| 1042 |
|
|---|
| 1043 | j0 = j+0
|
|---|
| 1044 | j3 = j+3
|
|---|
| 1045 | j6 = j+6
|
|---|
| 1046 | do i=1,nrstd ! rx has mass matrix and Jacobian on fine mesh
|
|---|
| 1047 | ud(i)=ud(i)
|
|---|
| 1048 | $ +rx(i,j0,e)*ur(i)+rx(i,j3,e)*us(i)+rx(i,j6,e)*ut(i)
|
|---|
| 1049 | enddo
|
|---|
| 1050 | enddo
|
|---|
| 1051 |
|
|---|
| 1052 | call intp_rstd(du(1,e),ud,mx,md,if3d,1) ! 1 --> transpose
|
|---|
| 1053 |
|
|---|
| 1054 | enddo
|
|---|
| 1055 |
|
|---|
| 1056 | return
|
|---|
| 1057 | end
|
|---|
| 1058 | c-----------------------------------------------------------------------
|
|---|
| 1059 | subroutine convop_cons_2d(du,u,c,mx,md,nel) ! Conservation form
|
|---|
| 1060 |
|
|---|
| 1061 | c Apply convecting field c to scalar field u, conservation form d/dxj cj phi
|
|---|
| 1062 |
|
|---|
| 1063 | c Assumes that current convecting field is on dealias mesh, in c()
|
|---|
| 1064 |
|
|---|
| 1065 | include 'SIZE'
|
|---|
| 1066 | include 'GEOM'
|
|---|
| 1067 | include 'TSTEP'
|
|---|
| 1068 |
|
|---|
| 1069 |
|
|---|
| 1070 | real du(mx*mx,nel)
|
|---|
| 1071 | real u(mx*mx,nel)
|
|---|
| 1072 | real c(md*md,nel,2)
|
|---|
| 1073 | parameter (ldd=lxd*lyd*lzd)
|
|---|
| 1074 | common /ctmp1/ ur(ldd),us(ldd),ut(ldd),ju(ldd),ud(ldd),tu(ldd)
|
|---|
| 1075 | real ju
|
|---|
| 1076 |
|
|---|
| 1077 | logical if3d,ifd
|
|---|
| 1078 | integer e
|
|---|
| 1079 |
|
|---|
| 1080 | if3d = .false.
|
|---|
| 1081 | ifd = .false.
|
|---|
| 1082 | if (md.ne.mx) ifd=.true.
|
|---|
| 1083 |
|
|---|
| 1084 | nrstd = md**2
|
|---|
| 1085 | call lim_chk(nrstd,ldd,'urus ','ldd ','convp_cons')
|
|---|
| 1086 |
|
|---|
| 1087 | if (nio.eq.0.and.istep.lt.3) write(6,*) 'convp_cons',istep
|
|---|
| 1088 |
|
|---|
| 1089 | do e=1,nel
|
|---|
| 1090 |
|
|---|
| 1091 | call intp_rstd (ju,u(1,e),mx,md,if3d,0) ! 0 = forward; on Gauss points!
|
|---|
| 1092 | call rzero (ud,nrstd)
|
|---|
| 1093 |
|
|---|
| 1094 | c call outmat(c(1,e,1),md,md,'fine u',e)
|
|---|
| 1095 | c call outmat(c(1,e,2),md,md,'fine v',e)
|
|---|
| 1096 | c call outmat(ju ,md,md,'fine T',e)
|
|---|
| 1097 |
|
|---|
| 1098 | do j=1,ldim
|
|---|
| 1099 | do i=1,nrstd
|
|---|
| 1100 | tu(i)=c(i,e,j)*ju(i) ! C_j*T
|
|---|
| 1101 | enddo
|
|---|
| 1102 | call grad_rst(ur,us,ut,tu,md,if3d) ! Already on fine (Gauss) mesh
|
|---|
| 1103 |
|
|---|
| 1104 | j0 = j+0
|
|---|
| 1105 | j2 = j+2
|
|---|
| 1106 | do i=1,nrstd ! rx has mass matrix and Jacobian on fine mesh
|
|---|
| 1107 | ud(i)=ud(i)+rx(i,j0,e)*ur(i)+rx(i,j2,e)*us(i)
|
|---|
| 1108 | enddo
|
|---|
| 1109 | enddo
|
|---|
| 1110 |
|
|---|
| 1111 | call intp_rstd(du(1,e),ud,mx,md,if3d,1) ! 1 --> transpose
|
|---|
| 1112 |
|
|---|
| 1113 | enddo
|
|---|
| 1114 |
|
|---|
| 1115 | c call exitti('convop_cons_2d$',istep)
|
|---|
| 1116 |
|
|---|
| 1117 | return
|
|---|
| 1118 | end
|
|---|
| 1119 | c-----------------------------------------------------------------------
|
|---|
| 1120 | subroutine iface_vert_int8(fa,va,jz0,jz1,nel)
|
|---|
| 1121 | include 'SIZE'
|
|---|
| 1122 | integer*8 fa(lx1*lz1,2*ldim,nel),va(0:lx1+1,0:ly1+1,jz0:jz1,nel)
|
|---|
| 1123 | integer e,f
|
|---|
| 1124 |
|
|---|
| 1125 | n = lx1*lz1*2*ldim*nel
|
|---|
| 1126 | call i8zero(fa,n)
|
|---|
| 1127 |
|
|---|
| 1128 | mx1 = lx1+2
|
|---|
| 1129 | my1 = ly1+2
|
|---|
| 1130 | mz1 = lz1+2
|
|---|
| 1131 | if (ldim.eq.2) mz1=1
|
|---|
| 1132 |
|
|---|
| 1133 | nface = 2*ldim
|
|---|
| 1134 | do e=1,nel
|
|---|
| 1135 | do f=1,nface
|
|---|
| 1136 | call facind (kx1,kx2,ky1,ky2,kz1,kz2,lx1,ly1,lz1,f)
|
|---|
| 1137 |
|
|---|
| 1138 | if (f.eq.1) then ! EB notation
|
|---|
| 1139 | ky1=ky1-1
|
|---|
| 1140 | ky2=ky1
|
|---|
| 1141 | elseif (f.eq.2) then
|
|---|
| 1142 | kx1=kx1+1
|
|---|
| 1143 | kx2=kx1
|
|---|
| 1144 | elseif (f.eq.3) then
|
|---|
| 1145 | ky1=ky1+1
|
|---|
| 1146 | ky2=ky1
|
|---|
| 1147 | elseif (f.eq.4) then
|
|---|
| 1148 | kx1=kx1-1
|
|---|
| 1149 | kx2=kx1
|
|---|
| 1150 | elseif (f.eq.5) then
|
|---|
| 1151 | kz1=kz1-1
|
|---|
| 1152 | kz2=kz1
|
|---|
| 1153 | elseif (f.eq.6) then
|
|---|
| 1154 | kz1=kz1+1
|
|---|
| 1155 | kz2=kz1
|
|---|
| 1156 | endif
|
|---|
| 1157 |
|
|---|
| 1158 | i = 0
|
|---|
| 1159 | do iz=kz1,kz2
|
|---|
| 1160 | do iy=ky1,ky2
|
|---|
| 1161 | do ix=kx1,kx2
|
|---|
| 1162 | i = i+1
|
|---|
| 1163 | fa(i,f,e)=va(ix,iy,iz,e)
|
|---|
| 1164 | c write(6,*) 'fa:',fa(i,f,e),i,f,e
|
|---|
| 1165 | enddo
|
|---|
| 1166 | enddo
|
|---|
| 1167 | enddo
|
|---|
| 1168 | enddo
|
|---|
| 1169 | enddo
|
|---|
| 1170 |
|
|---|
| 1171 | return
|
|---|
| 1172 | end
|
|---|
| 1173 | c-----------------------------------------------------------------------
|
|---|
| 1174 | subroutine set_binv(bmnv,hmsk,n) ! Store binvm1*(hyperbolic mask)
|
|---|
| 1175 |
|
|---|
| 1176 | include 'SIZE'
|
|---|
| 1177 | include 'PARALLEL'
|
|---|
| 1178 | include 'MASS'
|
|---|
| 1179 |
|
|---|
| 1180 | real bmnv(n,lorder),hmsk(n)
|
|---|
| 1181 |
|
|---|
| 1182 | do i=lorder,2,-1
|
|---|
| 1183 | call copy(bmnv(1,i),bmnv(1,i-1),n)
|
|---|
| 1184 | enddo
|
|---|
| 1185 |
|
|---|
| 1186 | call copy (bmnv,bm1,n) ! Fill bmnv(1,1)
|
|---|
| 1187 |
|
|---|
| 1188 | call fgslib_gs_op(gsh_fld(1),bmnv,1,1,0) ! 1 ==> +; gsh_fld(1) is velocity
|
|---|
| 1189 |
|
|---|
| 1190 | do i=1,n
|
|---|
| 1191 | bmnv(i,1)=hmsk(i)/bmnv(i,1)
|
|---|
| 1192 | enddo
|
|---|
| 1193 |
|
|---|
| 1194 | return
|
|---|
| 1195 | end
|
|---|
| 1196 | c-----------------------------------------------------------------------
|
|---|
| 1197 | subroutine set_bdivw(bdivw,hmsk,n) ! Store binvm1*(hyperbolic mask)
|
|---|
| 1198 |
|
|---|
| 1199 | include 'SIZE'
|
|---|
| 1200 | include 'PARALLEL'
|
|---|
| 1201 | include 'MASS'
|
|---|
| 1202 | include 'INPUT'
|
|---|
| 1203 | include 'MVGEOM'
|
|---|
| 1204 | common /scruz/ cx (lx1*ly1*lz1*lelt)
|
|---|
| 1205 | $ , cy (lx1*ly1*lz1*lelt)
|
|---|
| 1206 | $ , cz (lx1*ly1*lz1*lelt)
|
|---|
| 1207 |
|
|---|
| 1208 | real bdivw(n,lorder),hmsk(n)
|
|---|
| 1209 |
|
|---|
| 1210 | do i=lorder,2,-1
|
|---|
| 1211 | call copy(bdivw(1,i),bdivw(1,i-1),n)
|
|---|
| 1212 | enddo
|
|---|
| 1213 |
|
|---|
| 1214 | call gradm1 (bdivw,cy ,cz , wx )
|
|---|
| 1215 | call gradm1 (cx ,cy ,cz , wy )
|
|---|
| 1216 | call add2 (bdivw,cy , n )
|
|---|
| 1217 | if (if3d) then
|
|---|
| 1218 | call gradm1 (cx ,cy ,cz , wz )
|
|---|
| 1219 | call add2 (bdivw,cz , n )
|
|---|
| 1220 | endif
|
|---|
| 1221 | call col2(bdivw,bm1,n)
|
|---|
| 1222 |
|
|---|
| 1223 | return
|
|---|
| 1224 | end
|
|---|
| 1225 | c-----------------------------------------------------------------------
|
|---|
| 1226 | subroutine set_bmass(bmass,hmsk,n) ! Store bmass*(hyperbolic mask)
|
|---|
| 1227 |
|
|---|
| 1228 | include 'SIZE'
|
|---|
| 1229 | include 'PARALLEL'
|
|---|
| 1230 | include 'MASS'
|
|---|
| 1231 |
|
|---|
| 1232 | real bmass(n,lorder),hmsk(n)
|
|---|
| 1233 |
|
|---|
| 1234 | do i=lorder,2,-1
|
|---|
| 1235 | call copy(bmass(1,i),bmass(1,i-1),n)
|
|---|
| 1236 | enddo
|
|---|
| 1237 |
|
|---|
| 1238 | call copy (bmass,bm1,n) ! Fill bmass(1,1)
|
|---|
| 1239 |
|
|---|
| 1240 | return
|
|---|
| 1241 | end
|
|---|
| 1242 | c-----------------------------------------------------------------------
|
|---|
| 1243 | subroutine setup_dg_gs(dgh,nx,ny,nz,nel,melg,vertex)
|
|---|
| 1244 |
|
|---|
| 1245 | c Global-to-local mapping for gs
|
|---|
| 1246 |
|
|---|
| 1247 | include 'SIZE'
|
|---|
| 1248 | include 'TOTAL'
|
|---|
| 1249 |
|
|---|
| 1250 | integer dgh,vertex(1)
|
|---|
| 1251 |
|
|---|
| 1252 | parameter(lf=lx1*lz1*2*ldim*lelt)
|
|---|
| 1253 | common /c_is1/ glo_num_face(lf)
|
|---|
| 1254 | $ , glo_num_vol((lx1+2)*(ly1+2)*(lz1+2)*lelt)
|
|---|
| 1255 | integer*8 glo_num_face,glo_num_vol,ngv
|
|---|
| 1256 |
|
|---|
| 1257 | common /nekmpi/ mid,mp,nekcomm,nekgroup,nekreal
|
|---|
| 1258 |
|
|---|
| 1259 | mx = nx+2
|
|---|
| 1260 | call set_vert(glo_num_vol,ngv,mx,nel,vertex,.false.)
|
|---|
| 1261 |
|
|---|
| 1262 | mz0 = 1
|
|---|
| 1263 | mz1 = 1
|
|---|
| 1264 | if (if3d) mz0 = 0
|
|---|
| 1265 | if (if3d) mz1 = lz1+1
|
|---|
| 1266 | call iface_vert_int8 (glo_num_face,glo_num_vol,mz0,mz1,nelt)
|
|---|
| 1267 |
|
|---|
| 1268 | nf = lx1*lz1*2*ldim*nelt !total number of points on faces
|
|---|
| 1269 | call fgslib_gs_setup(dgh,glo_num_face,nf,nekcomm,np)
|
|---|
| 1270 |
|
|---|
| 1271 | return
|
|---|
| 1272 | end
|
|---|
| 1273 | c-----------------------------------------------------------------------
|
|---|
| 1274 | subroutine dg_set_fc_ptr
|
|---|
| 1275 |
|
|---|
| 1276 | c Set up pointer to restrict u to faces ! NOTE: compact
|
|---|
| 1277 |
|
|---|
| 1278 | include 'SIZE'
|
|---|
| 1279 | include 'TOTAL'
|
|---|
| 1280 |
|
|---|
| 1281 | integer e,f,ef
|
|---|
| 1282 |
|
|---|
| 1283 | call dsset(lx1,ly1,lz1) ! set skpdat
|
|---|
| 1284 |
|
|---|
| 1285 | nxyz = lx1*ly1*lz1
|
|---|
| 1286 | nxz = lx1*lz1
|
|---|
| 1287 | nface = 2*ldim
|
|---|
| 1288 | nxzf = lx1*lz1*nface ! red'd mod to area, unx, etc.
|
|---|
| 1289 |
|
|---|
| 1290 | k = 0
|
|---|
| 1291 |
|
|---|
| 1292 | do e=1,nelv
|
|---|
| 1293 | do ef=1,nface ! EB notation
|
|---|
| 1294 |
|
|---|
| 1295 | f = eface1(ef)
|
|---|
| 1296 | js1 = skpdat(1,f)
|
|---|
| 1297 | jf1 = skpdat(2,f)
|
|---|
| 1298 | jskip1 = skpdat(3,f)
|
|---|
| 1299 | js2 = skpdat(4,f)
|
|---|
| 1300 | jf2 = skpdat(5,f)
|
|---|
| 1301 | jskip2 = skpdat(6,f)
|
|---|
| 1302 |
|
|---|
| 1303 | i = 0
|
|---|
| 1304 | do j2=js2,jf2,jskip2
|
|---|
| 1305 | do j1=js1,jf1,jskip1
|
|---|
| 1306 |
|
|---|
| 1307 | i = i+1
|
|---|
| 1308 | k = i+nxz*(ef-1)+nxzf*(e-1) ! face numbering
|
|---|
| 1309 | dg_face(k) = j1+lx1*(j2-1)+nxyz*(e-1) ! global numbering
|
|---|
| 1310 |
|
|---|
| 1311 | enddo
|
|---|
| 1312 | enddo
|
|---|
| 1313 |
|
|---|
| 1314 | enddo
|
|---|
| 1315 | enddo
|
|---|
| 1316 | ndg_facex = nxzf*nelv
|
|---|
| 1317 |
|
|---|
| 1318 | return
|
|---|
| 1319 | end
|
|---|
| 1320 | c-----------------------------------------------------------------------
|
|---|
| 1321 | subroutine full2face(faceary, vol_ary)
|
|---|
| 1322 |
|
|---|
| 1323 | include 'SIZE'
|
|---|
| 1324 | include 'TOTAL'
|
|---|
| 1325 |
|
|---|
| 1326 | real faceary(lx1*lz1,2*ldim,lelt)
|
|---|
| 1327 | real vol_ary(lx1,ly1,lz1,lelt)
|
|---|
| 1328 | integer i,j
|
|---|
| 1329 |
|
|---|
| 1330 | do j=1,ndg_facex
|
|---|
| 1331 | i=dg_face(j)
|
|---|
| 1332 | faceary(j,1,1) = vol_ary(i,1,1,1)
|
|---|
| 1333 | enddo
|
|---|
| 1334 |
|
|---|
| 1335 |
|
|---|
| 1336 | return
|
|---|
| 1337 | end
|
|---|
| 1338 | c-----------------------------------------------------------------------
|
|---|
| 1339 | subroutine face2full(vol_ary, faceary)
|
|---|
| 1340 |
|
|---|
| 1341 | include 'SIZE'
|
|---|
| 1342 | include 'TOTAL'
|
|---|
| 1343 |
|
|---|
| 1344 | real faceary(lx1*lz1,2*ldim,lelt)
|
|---|
| 1345 | real vol_ary(lx1,ly1,lz1,lelt)
|
|---|
| 1346 | integer i,j
|
|---|
| 1347 |
|
|---|
| 1348 | n=lx1*ly1*lz1*nelfld(ifield)
|
|---|
| 1349 | call rzero(vol_ary,n)
|
|---|
| 1350 |
|
|---|
| 1351 | do j=1,ndg_facex
|
|---|
| 1352 | i=dg_face(j)
|
|---|
| 1353 | vol_ary(i,1,1,1) = vol_ary(i,1,1,1)+faceary(j,1,1)
|
|---|
| 1354 | enddo
|
|---|
| 1355 |
|
|---|
| 1356 | return
|
|---|
| 1357 | end
|
|---|
| 1358 | c-----------------------------------------------------------------------
|
|---|
| 1359 | subroutine add_face2full(vol_ary, faceary)
|
|---|
| 1360 |
|
|---|
| 1361 | include 'SIZE'
|
|---|
| 1362 | include 'TOTAL'
|
|---|
| 1363 |
|
|---|
| 1364 | real faceary(lx1*lz1,2*ldim,lelt)
|
|---|
| 1365 | real vol_ary(lx1,ly1,lz1,lelt)
|
|---|
| 1366 | integer i,j
|
|---|
| 1367 |
|
|---|
| 1368 | do j=1,ndg_facex
|
|---|
| 1369 | i=dg_face(j)
|
|---|
| 1370 | vol_ary(i,1,1,1) = vol_ary(i,1,1,1)+faceary(j,1,1)
|
|---|
| 1371 | enddo
|
|---|
| 1372 |
|
|---|
| 1373 | return
|
|---|
| 1374 | end
|
|---|
| 1375 | c-----------------------------------------------------------------------
|
|---|
| 1376 | subroutine conv_rhs_dg_aliased (du,u,c)
|
|---|
| 1377 | c
|
|---|
| 1378 | include 'SIZE'
|
|---|
| 1379 | include 'TOTAL'
|
|---|
| 1380 |
|
|---|
| 1381 | c Apply convecting field c(1,ldim) to scalar field u(1).
|
|---|
| 1382 |
|
|---|
| 1383 | real du(1),u(1),c(1)
|
|---|
| 1384 |
|
|---|
| 1385 | parameter(lf=lx1*lz1*2*ldim*lelt)
|
|---|
| 1386 | common /scrdg/ uf(lf),uxf(lf),uyf(lf),uzf(lf),upwind_wgt(lf)
|
|---|
| 1387 |
|
|---|
| 1388 | integer e,f
|
|---|
| 1389 |
|
|---|
| 1390 | n = lx1*ly1*lz1*nelv
|
|---|
| 1391 | nf = lx1*lz1*2*ldim*nelt
|
|---|
| 1392 |
|
|---|
| 1393 |
|
|---|
| 1394 | if (ifcons) then
|
|---|
| 1395 | if (if3d ) call convop_cons_3d (du,u,c,lx1,lxd,nelv)
|
|---|
| 1396 | if (.not.if3d) call convop_cons_2d (du,u,c,lx1,lxd,nelv)
|
|---|
| 1397 | else
|
|---|
| 1398 | if (if3d ) call convop_fst_3d (du,u,c,lx1,lxd,nelv)
|
|---|
| 1399 | if (.not.if3d) call convop_fst_2d (du,u,c,lx1,lxd,nelv)
|
|---|
| 1400 | endif
|
|---|
| 1401 |
|
|---|
| 1402 | call full2face(uf ,u )
|
|---|
| 1403 | call full2face(uxf,vx)
|
|---|
| 1404 | call full2face(uyf,vy)
|
|---|
| 1405 | call full2face(uzf,vz)
|
|---|
| 1406 | if (.not.if3d) call rzero(uzf,nf)
|
|---|
| 1407 |
|
|---|
| 1408 | beta_u = 0.00 ! 1=full upwind; 0=central flux
|
|---|
| 1409 | beta_u = 0.25 ! 1=full upwind; 0=central flux
|
|---|
| 1410 | beta_u = 1.00 ! 1=full upwind; 0=central flux
|
|---|
| 1411 | beta_u = param(98)
|
|---|
| 1412 | if (istep.le.5.and.nio.eq.0) write(6,*) beta_u,' dg upwind'
|
|---|
| 1413 |
|
|---|
| 1414 | nface = 2*ldim
|
|---|
| 1415 | nxz = lx1*lz1
|
|---|
| 1416 | k = 0
|
|---|
| 1417 | do e=1,nelt
|
|---|
| 1418 | do f=1,nface
|
|---|
| 1419 | do i=1,nxz
|
|---|
| 1420 |
|
|---|
| 1421 | k=k+1
|
|---|
| 1422 |
|
|---|
| 1423 | beta = ( unx (i,1,f,e)*uxf(k)
|
|---|
| 1424 | $ + uny (i,1,f,e)*uyf(k)
|
|---|
| 1425 | $ + unz (i,1,f,e)*uzf(k))
|
|---|
| 1426 |
|
|---|
| 1427 | uf(k) = -beta*area(i,1,f,e)*uf(k)
|
|---|
| 1428 |
|
|---|
| 1429 | upwind_wgt(k) = 1.0
|
|---|
| 1430 | if (beta.gt.0) upwind_wgt(k) = 0.0
|
|---|
| 1431 | upwind_wgt(k) = 0.5*(1-beta_u) + upwind_wgt(k)*beta_u
|
|---|
| 1432 | if (beta.eq.0) upwind_wgt(k) = 0.5
|
|---|
| 1433 |
|
|---|
| 1434 | enddo
|
|---|
| 1435 | enddo
|
|---|
| 1436 | enddo
|
|---|
| 1437 |
|
|---|
| 1438 | call fgslib_gs_op(dg_hndlx,uf,1,1,0) ! 1 ==> +
|
|---|
| 1439 | call col2 (uf,upwind_wgt,nf) ! Inefficient, but ok for now.
|
|---|
| 1440 | ! Should be combined with
|
|---|
| 1441 | call add_face2full( du , uf) ! <--- this stmt.
|
|---|
| 1442 |
|
|---|
| 1443 | call fbinvert ( du ) ! Right now works only for undeformed
|
|---|
| 1444 |
|
|---|
| 1445 | return
|
|---|
| 1446 | end
|
|---|
| 1447 | c-----------------------------------------------------------------------
|
|---|
| 1448 | subroutine map_faced(ju,u,mx,md,fdim,idir) ! GLL->GL interpolation
|
|---|
| 1449 |
|
|---|
| 1450 | c GLL interpolation from mx to md for a face array of size (nx,nz)
|
|---|
| 1451 |
|
|---|
| 1452 | c If idir ^= 0, then apply transpose operator (md to mx)
|
|---|
| 1453 |
|
|---|
| 1454 | include 'SIZE'
|
|---|
| 1455 |
|
|---|
| 1456 | real ju(1),u(1)
|
|---|
| 1457 | integer fdim
|
|---|
| 1458 |
|
|---|
| 1459 | parameter (ldg=lxd**3,lwkd=4*lxd*lxd)
|
|---|
| 1460 | common /dgrad/ d(ldg),dt(ldg),dg(ldg),dgt(ldg),jgl(ldg),jgt(ldg)
|
|---|
| 1461 | $ , wkd(lwkd)
|
|---|
| 1462 | real jgl,jgt
|
|---|
| 1463 |
|
|---|
| 1464 | parameter (ld=2*lxd)
|
|---|
| 1465 | common /ctmp0/ w(ld**ldim,2)
|
|---|
| 1466 |
|
|---|
| 1467 | call lim_chk(md,ld,'md ','ld ','map_faced ')
|
|---|
| 1468 | call lim_chk(mx,ld,'mx ','ld ','map_faced ')
|
|---|
| 1469 |
|
|---|
| 1470 | c call copy(ju,u,mx)
|
|---|
| 1471 | c return
|
|---|
| 1472 |
|
|---|
| 1473 | call get_int_ptr (i,mx,md)
|
|---|
| 1474 |
|
|---|
| 1475 | if (idir.eq.0) then
|
|---|
| 1476 | if (fdim.eq.2) then
|
|---|
| 1477 | call mxm(jgl(i),md,u,mx,wkd,mx)
|
|---|
| 1478 | call mxm(wkd,md,jgt(i),mx,ju,md)
|
|---|
| 1479 | else
|
|---|
| 1480 | call mxm(jgl(i),md,u,mx,ju,1)
|
|---|
| 1481 | endif
|
|---|
| 1482 | else
|
|---|
| 1483 | if (fdim.eq.2) then
|
|---|
| 1484 | call mxm(jgt(i),mx,u,md,wkd,md)
|
|---|
| 1485 | call mxm(wkd,mx,jgl(i),md,ju,mx)
|
|---|
| 1486 | else
|
|---|
| 1487 | call mxm(jgt(i),mx,u,md,ju,1)
|
|---|
| 1488 | endif
|
|---|
| 1489 | endif
|
|---|
| 1490 |
|
|---|
| 1491 | return
|
|---|
| 1492 | end
|
|---|
| 1493 | c-----------------------------------------------------------------------
|
|---|
| 1494 | subroutine fbinvert(rhs) ! Still in development. 10/10/15, pff.
|
|---|
| 1495 |
|
|---|
| 1496 | include 'SIZE'
|
|---|
| 1497 | include 'TOTAL'
|
|---|
| 1498 |
|
|---|
| 1499 | real rhs(lx1,ly1,lz1,lelt)
|
|---|
| 1500 |
|
|---|
| 1501 | common /cfbinv/ qn(lx1),alpha_n,beta_n
|
|---|
| 1502 | $ ,s1(ly1,lz1),bnv(lx1)
|
|---|
| 1503 | $ ,tmp(lx1*ly1*lz1*lelt)
|
|---|
| 1504 | integer icalld
|
|---|
| 1505 | save icalld
|
|---|
| 1506 | data icalld /0/
|
|---|
| 1507 |
|
|---|
| 1508 | integer e
|
|---|
| 1509 |
|
|---|
| 1510 | n = lx1*ly1*lz1*nelfld(ifield)
|
|---|
| 1511 |
|
|---|
| 1512 | do i=1,n ! FOR NOW, USE DIAGONAL MASS
|
|---|
| 1513 | rhs(i,1,1,1)=rhs(i,1,1,1)/bm1(i,1,1,1) ! MATRIX. pff, 10/10/15
|
|---|
| 1514 | enddo
|
|---|
| 1515 |
|
|---|
| 1516 | return
|
|---|
| 1517 | end
|
|---|
| 1518 | c-----------------------------------------------------------------------
|
|---|
| 1519 | subroutine grad_rstd_ta(du,ur,us,ut,md,if3d) ! GL->GL gradt
|
|---|
| 1520 |
|
|---|
| 1521 | include 'SIZE'
|
|---|
| 1522 | include 'DXYZ'
|
|---|
| 1523 |
|
|---|
| 1524 | real ur(1),us(1),ut(1),u(1)
|
|---|
| 1525 |
|
|---|
| 1526 | logical if3d
|
|---|
| 1527 |
|
|---|
| 1528 | parameter (ldg=lxd**3,lwkd=4*lxd*lxd)
|
|---|
| 1529 | common /dgrad/ d(ldg),dt(ldg),dg(ldg),dgt(ldg),jgl(ldg),jgt(ldg)
|
|---|
| 1530 | $ , wkd(lwkd)
|
|---|
| 1531 | real jgl,jgt
|
|---|
| 1532 |
|
|---|
| 1533 | call get_dgl_ptr (ip,md,md)
|
|---|
| 1534 | call gradrta (du,ur,us,ut,dgt(ip),dg(ip),dg(ip),md,md,md,if3d)
|
|---|
| 1535 |
|
|---|
| 1536 | return
|
|---|
| 1537 | end
|
|---|
| 1538 | c-----------------------------------------------------------------------
|
|---|
| 1539 | subroutine set_dg_wgts
|
|---|
| 1540 |
|
|---|
| 1541 | include 'SIZE'
|
|---|
| 1542 | include 'TOTAL'
|
|---|
| 1543 |
|
|---|
| 1544 | integer icalld
|
|---|
| 1545 | save icalld
|
|---|
| 1546 | data icalld /0/
|
|---|
| 1547 |
|
|---|
| 1548 | common /finewts/ zptf(lxd),wgtf(lxd),wghtf(lxd*lzd),wghtc(lx1*lz1)
|
|---|
| 1549 |
|
|---|
| 1550 | if (icalld.eq.0) then ! Set fine-scale surface weights
|
|---|
| 1551 |
|
|---|
| 1552 | icalld = 1
|
|---|
| 1553 |
|
|---|
| 1554 | call zwgl(zptf,wgtf,lxd)
|
|---|
| 1555 | if (if3d) then
|
|---|
| 1556 | k=0
|
|---|
| 1557 | do j=1,ly1
|
|---|
| 1558 | do i=1,lx1
|
|---|
| 1559 | k=k+1
|
|---|
| 1560 | wghtc(k)=wxm1(i)*wzm1(j)
|
|---|
| 1561 | enddo
|
|---|
| 1562 | enddo
|
|---|
| 1563 | k=0
|
|---|
| 1564 | do j=1,lyd
|
|---|
| 1565 | do i=1,lxd
|
|---|
| 1566 | k=k+1
|
|---|
| 1567 | wghtf(k)=wgtf(i)*wgtf(j)
|
|---|
| 1568 | enddo
|
|---|
| 1569 | enddo
|
|---|
| 1570 | else
|
|---|
| 1571 | call copy(wghtc,wxm1,lx1)
|
|---|
| 1572 | call copy(wghtf,wgtf,lxd)
|
|---|
| 1573 | endif
|
|---|
| 1574 | endif
|
|---|
| 1575 |
|
|---|
| 1576 | return
|
|---|
| 1577 | end
|
|---|
| 1578 | c-----------------------------------------------------------------------
|
|---|
| 1579 | subroutine conv_rhs_dg (du,u,c)
|
|---|
| 1580 | c
|
|---|
| 1581 | include 'SIZE'
|
|---|
| 1582 | include 'TOTAL'
|
|---|
| 1583 |
|
|---|
| 1584 | c Apply convecting field c(1,ldim) to scalar field u(1).
|
|---|
| 1585 |
|
|---|
| 1586 | parameter(ldd=lxd*lyd*lzd)
|
|---|
| 1587 | real du(1),u(1),c(ldd*lelv,3)
|
|---|
| 1588 |
|
|---|
| 1589 | parameter(lf=lx1*lz1*2*ldim*lelt)
|
|---|
| 1590 | common /scrdg/ uf(lf),uxf(lf),uyf(lf),uzf(lf),upwind_wgt(lf)
|
|---|
| 1591 | $ , beta_c(lx1*lz1),jaco_c(lx1*lz1)
|
|---|
| 1592 | $ , beta_f(lxd*lzd),jaco_f(lxd*lzd)
|
|---|
| 1593 | $ , ufine (lxd*lzd)
|
|---|
| 1594 | real jaco_c,jaco_f
|
|---|
| 1595 | common /finewts/ zptf(lxd),wgtf(lxd),wghtf(lxd*lzd),wghtc(lx1*lz1)
|
|---|
| 1596 |
|
|---|
| 1597 | integer e,f,fdim
|
|---|
| 1598 |
|
|---|
| 1599 |
|
|---|
| 1600 | n = lx1*ly1*lz1*lelv
|
|---|
| 1601 | nf = lx1*lz1*2*ldim*lelt
|
|---|
| 1602 |
|
|---|
| 1603 | call conv_rhs_dg_weak (du,u,c(1,1),c(1,2),c(1,3))
|
|---|
| 1604 |
|
|---|
| 1605 | call full2face(uf ,u )
|
|---|
| 1606 | call full2face(uxf,vx)
|
|---|
| 1607 | call full2face(uyf,vy)
|
|---|
| 1608 | call full2face(uzf,vz)
|
|---|
| 1609 | if (.not.if3d) call rzero(uzf,nf)
|
|---|
| 1610 |
|
|---|
| 1611 | beta_u = 0.00 ! 1=full upwind; 0=central flux
|
|---|
| 1612 | beta_u = 0.25 ! 1=full upwind; 0=central flux
|
|---|
| 1613 | beta_u = 1.00 ! 1=full upwind; 0=central flux
|
|---|
| 1614 | beta_u = param(98)
|
|---|
| 1615 | if (istep.le.5.and.nio.eq.0) write(6,*) beta_u,' dg upwind'
|
|---|
| 1616 |
|
|---|
| 1617 | nface = 2*ldim
|
|---|
| 1618 | nxz = lx1*lz1
|
|---|
| 1619 | nxzd = lxd*lzd
|
|---|
| 1620 | k = 0
|
|---|
| 1621 | do e=1,nelt ! This formula for upwind weights appears to
|
|---|
| 1622 | do f=1,nface ! assume that U is continuous, or at least of
|
|---|
| 1623 |
|
|---|
| 1624 | kface = k+1
|
|---|
| 1625 | do i=1,nxz ! the same sign on either side of the interface.
|
|---|
| 1626 |
|
|---|
| 1627 | k=k+1
|
|---|
| 1628 |
|
|---|
| 1629 | beta = ( unx (i,1,f,e)*uxf(k)
|
|---|
| 1630 | $ + uny (i,1,f,e)*uyf(k)
|
|---|
| 1631 | $ + unz (i,1,f,e)*uzf(k))
|
|---|
| 1632 |
|
|---|
| 1633 | upwind_wgt(k) = 1.0
|
|---|
| 1634 | if (beta.gt.0) upwind_wgt(k) = 0.0
|
|---|
| 1635 | upwind_wgt(k) = 0.5*(1-beta_u) + upwind_wgt(k)*beta_u
|
|---|
| 1636 |
|
|---|
| 1637 | beta_c(i) = -beta
|
|---|
| 1638 | jaco_c(i) = area(i,1,f,e)/wghtc(i)
|
|---|
| 1639 |
|
|---|
| 1640 | enddo
|
|---|
| 1641 |
|
|---|
| 1642 | fdim = ldim-1 ! Dimension of face
|
|---|
| 1643 | call map_faced(beta_f,beta_c ,lx1,lxd,fdim,0) ! Dealiased quadrature,
|
|---|
| 1644 | call map_faced(jaco_f,jaco_c ,lx1,lxd,fdim,0) ! 0 --> coarse to fine,
|
|---|
| 1645 | call map_faced(ufine,uf(kface),lx1,lxd,fdim,0) ! ufine = J uf
|
|---|
| 1646 |
|
|---|
| 1647 | do i=1,nxzd
|
|---|
| 1648 | ufine(i)=wghtf(i)*jaco_f(i)*beta_f(i)*ufine(i)
|
|---|
| 1649 | enddo
|
|---|
| 1650 | call map_faced(uf(kface),ufine,lx1,lxd,fdim,1) ! 1 --> uf = J^T ufine
|
|---|
| 1651 |
|
|---|
| 1652 | enddo
|
|---|
| 1653 | enddo
|
|---|
| 1654 |
|
|---|
| 1655 | call fgslib_gs_op(dg_hndlx,uf,1,1,0) ! 1 ==> +
|
|---|
| 1656 |
|
|---|
| 1657 | call col2 (uf,upwind_wgt,nf) ! Inefficient, but ok for now.
|
|---|
| 1658 | ! Should be combined with
|
|---|
| 1659 | call add_face2full( du , uf) ! <--- this stmt.
|
|---|
| 1660 |
|
|---|
| 1661 | call fbinvert ( du ) ! Right now works only for undeformed
|
|---|
| 1662 |
|
|---|
| 1663 | return
|
|---|
| 1664 | end
|
|---|
| 1665 | c-----------------------------------------------------------------------
|
|---|
| 1666 | subroutine convect_dg(du,u,ifuf,cr,cs,ct,ifcf)
|
|---|
| 1667 | include 'SIZE'
|
|---|
| 1668 | include 'TOTAL'
|
|---|
| 1669 | real du(1),u(1),cr(1),cs(1),ct(1)
|
|---|
| 1670 | logical ifuf,ifcf
|
|---|
| 1671 |
|
|---|
| 1672 | call conv_rhs_dg_weak (du,u,cr,cs,ct)
|
|---|
| 1673 |
|
|---|
| 1674 | return
|
|---|
| 1675 | end
|
|---|
| 1676 | c-----------------------------------------------------------------------
|
|---|
| 1677 | subroutine convop_weak(du,u,cr,cs,ct,mx,md,nel) ! Weak Conservation form
|
|---|
| 1678 |
|
|---|
| 1679 | c Apply convecting field c to scalar field u, conservation form d/dxj cj phi
|
|---|
| 1680 |
|
|---|
| 1681 | c Assumes that current convecting field is on dealias mesh, in c()
|
|---|
| 1682 |
|
|---|
| 1683 | include 'SIZE'
|
|---|
| 1684 | include 'TOTAL'
|
|---|
| 1685 |
|
|---|
| 1686 | parameter (lxx=lx1*ly1*lz1,ldd=lxd*lyd*lzd)
|
|---|
| 1687 | real du(lxx,nel)
|
|---|
| 1688 | real u(lxx,nel)
|
|---|
| 1689 | real cr(ldd,nel),cs(ldd,nel),ct(ldd,nel)
|
|---|
| 1690 | common /ctmp1/ ur(ldd),us(ldd),ut(ldd),ju(ldd),ud(ldd),tu(ldd)
|
|---|
| 1691 | real ju
|
|---|
| 1692 |
|
|---|
| 1693 | integer e
|
|---|
| 1694 |
|
|---|
| 1695 | nxyz = lx1*ly1*lz1
|
|---|
| 1696 | nrstd = md**ldim
|
|---|
| 1697 |
|
|---|
| 1698 | call lim_chk(nrstd,ldd,'urus5','ldd ','convp_cons')
|
|---|
| 1699 |
|
|---|
| 1700 | do e=1,nel
|
|---|
| 1701 |
|
|---|
| 1702 | call intp_rstd (ju,u(1,e),mx,md,if3d,0) ! 0 = forward; on Gauss points!
|
|---|
| 1703 | if (ldim.eq.3) then
|
|---|
| 1704 | do i=1,ldd
|
|---|
| 1705 | ur(i)=ju(i)*cr(i,e) ! Already in r-s-t coordinates
|
|---|
| 1706 | us(i)=ju(i)*cs(i,e)
|
|---|
| 1707 | ut(i)=ju(i)*ct(i,e)
|
|---|
| 1708 | ud(i)=0
|
|---|
| 1709 | enddo
|
|---|
| 1710 | else
|
|---|
| 1711 | do i=1,ldd
|
|---|
| 1712 | ur(i)=ju(i)*cr(i,e) ! Already in r-s-t coordinates
|
|---|
| 1713 | us(i)=ju(i)*cs(i,e)
|
|---|
| 1714 | ud(i)=0
|
|---|
| 1715 | enddo
|
|---|
| 1716 | endif
|
|---|
| 1717 | call grad_rstd_ta (ud,ur,us,ut,lxd,if3d) ! GL->GL gradt
|
|---|
| 1718 | call intp_rstd (du(1,e),ud,mx,md,if3d,1) ! 1 = backward; on Gauss points!
|
|---|
| 1719 | do i=1,nxyz
|
|---|
| 1720 | du(i,e) = -du(i,e)*binvdg(i,e)
|
|---|
| 1721 | enddo
|
|---|
| 1722 |
|
|---|
| 1723 | enddo
|
|---|
| 1724 |
|
|---|
| 1725 | return
|
|---|
| 1726 | end
|
|---|
| 1727 | c-----------------------------------------------------------------------
|
|---|
| 1728 | subroutine conv_bdry_dg_weak (du,u) ! THIS SHOULD HAVE: ,cr,cs,ct)
|
|---|
| 1729 | c
|
|---|
| 1730 | c Implement Cu = Div (cu) in weak form using DG
|
|---|
| 1731 | c
|
|---|
| 1732 | include 'SIZE'
|
|---|
| 1733 | include 'TOTAL'
|
|---|
| 1734 |
|
|---|
| 1735 | c Apply convecting field c(1,ldim) to scalar field u(1).
|
|---|
| 1736 |
|
|---|
| 1737 | real du(1),u(1)
|
|---|
| 1738 |
|
|---|
| 1739 | parameter(lf=lx1*lz1*2*ldim*lelt)
|
|---|
| 1740 | common /scrdg/uf(lf),uxf(lf),uyf(lf),uzf(lf),upwind_wgt(lf),us(lf)
|
|---|
| 1741 | $ ,beta_c(lx1*lz1),jaco_c(lx1*lz1)
|
|---|
| 1742 | $ ,beta_f(lxd*lzd),jaco_f(lxd*lzd)
|
|---|
| 1743 | $ ,ufine (lxd*lzd)
|
|---|
| 1744 | real jaco_c,jaco_f
|
|---|
| 1745 |
|
|---|
| 1746 | common /finewts/ zptf(lxd),wgtf(lxd),wghtf(lxd*lzd),wghtc(lx1*lz1)
|
|---|
| 1747 |
|
|---|
| 1748 | integer e,f,fdim
|
|---|
| 1749 |
|
|---|
| 1750 | n = lx1*ly1*lz1*nelv
|
|---|
| 1751 | nf = lx1*lz1*2*ldim*nelt
|
|---|
| 1752 |
|
|---|
| 1753 | call full2face(uf ,u )
|
|---|
| 1754 | call full2face(uxf,vx)
|
|---|
| 1755 | call full2face(uyf,vy)
|
|---|
| 1756 | call full2face(uzf,vz)
|
|---|
| 1757 | if (.not.if3d) call rzero(uzf,nf)
|
|---|
| 1758 |
|
|---|
| 1759 | beta_u = 1.00 ! 1=full upwind; 0=central flux
|
|---|
| 1760 |
|
|---|
| 1761 | nface = 2*ldim
|
|---|
| 1762 | nxz = lx1*lz1
|
|---|
| 1763 | nxzd = lxd*lzd
|
|---|
| 1764 | k = 0
|
|---|
| 1765 | do e=1,nelt ! This formula for upwind weights appears to
|
|---|
| 1766 | do f=1,nface ! assume that U is continuous, or at least of
|
|---|
| 1767 |
|
|---|
| 1768 | kface = k+1
|
|---|
| 1769 | if (fw(f,e).gt.0.6) then
|
|---|
| 1770 |
|
|---|
| 1771 | do i=1,nxz ! the same sign on either side of the interface.
|
|---|
| 1772 |
|
|---|
| 1773 | k=k+1
|
|---|
| 1774 |
|
|---|
| 1775 | beta = ( unx (i,1,f,e)*uxf(k)
|
|---|
| 1776 | $ + uny (i,1,f,e)*uyf(k)
|
|---|
| 1777 | $ + unz (i,1,f,e)*uzf(k))
|
|---|
| 1778 |
|
|---|
| 1779 | upwind_wgt(k) = 0.0
|
|---|
| 1780 | if (beta.lt.0) upwind_wgt(k) = 1.0
|
|---|
| 1781 |
|
|---|
| 1782 | beta_c(i) = beta
|
|---|
| 1783 | jaco_c(i) = area(i,1,f,e)/wghtc(i)
|
|---|
| 1784 |
|
|---|
| 1785 | enddo
|
|---|
| 1786 |
|
|---|
| 1787 | fdim = ldim-1 ! Dimension of face
|
|---|
| 1788 | call map_faced(beta_f,beta_c ,lx1,lxd,fdim,0) ! Dealiased quadrature,
|
|---|
| 1789 | call map_faced(jaco_f,jaco_c ,lx1,lxd,fdim,0) ! 0 --> coarse to fine,
|
|---|
| 1790 | call map_faced(ufine,uf(kface),lx1,lxd,fdim,0) ! ufine = J uf
|
|---|
| 1791 |
|
|---|
| 1792 | do i=1,nxzd
|
|---|
| 1793 | ufine(i)=wghtf(i)*jaco_f(i)*beta_f(i)*ufine(i)
|
|---|
| 1794 | enddo
|
|---|
| 1795 | call map_faced(uf(kface),ufine,lx1,lxd,fdim,1) ! 1 --> uf = J^T ufine
|
|---|
| 1796 |
|
|---|
| 1797 | else
|
|---|
| 1798 | do i=1,nxz
|
|---|
| 1799 | k=k+1
|
|---|
| 1800 | upwind_wgt(k) = 0.0
|
|---|
| 1801 | enddo
|
|---|
| 1802 | endif
|
|---|
| 1803 |
|
|---|
| 1804 | enddo
|
|---|
| 1805 | enddo
|
|---|
| 1806 |
|
|---|
| 1807 | do j=1,ndg_facex
|
|---|
| 1808 | i=dg_face(j)
|
|---|
| 1809 | du(i) = du(i) - ( upwind_wgt(j)*uf(j) )
|
|---|
| 1810 | enddo
|
|---|
| 1811 |
|
|---|
| 1812 | return
|
|---|
| 1813 | end
|
|---|
| 1814 | c-----------------------------------------------------------------------
|
|---|
| 1815 | subroutine conv_rhs_dg_weak (du,u,cr,cs,ct)
|
|---|
| 1816 | c
|
|---|
| 1817 | c Implement Cu = Div (cu) in weak form using DG
|
|---|
| 1818 | c
|
|---|
| 1819 | include 'SIZE'
|
|---|
| 1820 | include 'TOTAL'
|
|---|
| 1821 |
|
|---|
| 1822 | c Apply convecting field c(1,ldim) to scalar field u(1).
|
|---|
| 1823 |
|
|---|
| 1824 | real du(1),u(1),cr(1),cs(1),ct(1)
|
|---|
| 1825 |
|
|---|
| 1826 | parameter(lf=lx1*lz1*2*ldim*lelt)
|
|---|
| 1827 | common /scrdg/uf(lf),uxf(lf),uyf(lf),uzf(lf),upwind_wgt(lf),us(lf)
|
|---|
| 1828 | $ ,beta_c(lx1*lz1),jaco_c(lx1*lz1)
|
|---|
| 1829 | $ ,beta_f(lxd*lzd),jaco_f(lxd*lzd)
|
|---|
| 1830 | $ ,ufine (lxd*lzd)
|
|---|
| 1831 | real jaco_c,jaco_f
|
|---|
| 1832 |
|
|---|
| 1833 | common /finewts/ zptf(lxd),wgtf(lxd),wghtf(lxd*lzd),wghtc(lx1*lz1)
|
|---|
| 1834 |
|
|---|
| 1835 | integer e,f,fdim
|
|---|
| 1836 |
|
|---|
| 1837 | n = lx1*ly1*lz1*lelv
|
|---|
| 1838 | nf = lx1*lz1*2*ldim*lelt
|
|---|
| 1839 |
|
|---|
| 1840 |
|
|---|
| 1841 | call convop_weak (du,u,cr,cs,ct,lx1,lxd,nelv) ! Volumetric term
|
|---|
| 1842 |
|
|---|
| 1843 | call full2face(uf ,u )
|
|---|
| 1844 | call full2face(uxf,vx)
|
|---|
| 1845 | call full2face(uyf,vy)
|
|---|
| 1846 | call full2face(uzf,vz)
|
|---|
| 1847 | if (.not.if3d) call rzero(uzf,nf)
|
|---|
| 1848 |
|
|---|
| 1849 | beta_u = 0.25 ! 1=full upwind; 0=central flux
|
|---|
| 1850 | beta_u = 0.00 ! 1=full upwind; 0=central flux
|
|---|
| 1851 | beta_u = 1.00 ! 1=full upwind; 0=central flux
|
|---|
| 1852 | if (istep.le.5.and.nio.eq.0) write(6,*) beta_u,' dg upwind'
|
|---|
| 1853 |
|
|---|
| 1854 | nface = 2*ldim
|
|---|
| 1855 | nxz = lx1*lz1
|
|---|
| 1856 | nxzd = lxd*lzd
|
|---|
| 1857 | k = 0
|
|---|
| 1858 | do e=1,nelt ! This formula for upwind weights appears to
|
|---|
| 1859 | do f=1,nface ! assume that U is continuous, or at least of
|
|---|
| 1860 |
|
|---|
| 1861 | kface = k+1
|
|---|
| 1862 | do i=1,nxz ! the same sign on either side of the interface.
|
|---|
| 1863 | k=k+1
|
|---|
| 1864 | beta = ( unx (i,1,f,e)*uxf(k)
|
|---|
| 1865 | $ + uny (i,1,f,e)*uyf(k)
|
|---|
| 1866 | $ + unz (i,1,f,e)*uzf(k) )
|
|---|
| 1867 |
|
|---|
| 1868 | upwind_wgt(k) = 0.0
|
|---|
| 1869 | if (beta.gt.0) upwind_wgt(k) = 1.0
|
|---|
| 1870 | upwind_wgt(k) = 0.5*(1-beta_u) + beta_u*(1-upwind_wgt(k))
|
|---|
| 1871 |
|
|---|
| 1872 | if (fw(f,e).gt.0.6 .and. beta.lt.0) upwind_wgt(k)=1.
|
|---|
| 1873 | if (fw(f,e).gt.0.6 .and. beta.gt.0) upwind_wgt(k)=0.
|
|---|
| 1874 |
|
|---|
| 1875 | beta_c(i) = beta
|
|---|
| 1876 | jaco_c(i) = area(i,1,f,e)/wghtc(i)
|
|---|
| 1877 | enddo
|
|---|
| 1878 |
|
|---|
| 1879 | fdim = ldim-1 ! Dimension of face
|
|---|
| 1880 | call map_faced(beta_f,beta_c ,lx1,lxd,fdim,0) ! Dealiased quadrature,
|
|---|
| 1881 | call map_faced(jaco_f,jaco_c ,lx1,lxd,fdim,0) ! 0 --> coarse to fine,
|
|---|
| 1882 | call map_faced(ufine,uf(kface),lx1,lxd,fdim,0) ! ufine = J uf
|
|---|
| 1883 |
|
|---|
| 1884 | do i=1,nxzd
|
|---|
| 1885 | ufine(i)=wghtf(i)*jaco_f(i)*beta_f(i)*ufine(i)
|
|---|
| 1886 | enddo
|
|---|
| 1887 | call map_faced(uf(kface),ufine,lx1,lxd,fdim,1) ! 1 --> uf = J^T ufine
|
|---|
| 1888 | call copy (us(kface),uf(kface),lx1*lz1) ! Save uf for later recombination
|
|---|
| 1889 |
|
|---|
| 1890 | enddo
|
|---|
| 1891 | enddo
|
|---|
| 1892 |
|
|---|
| 1893 | call fgslib_gs_op(dg_hndlx,uf,1,1,0) ! 1 ==> + : uf <-- uf^- + uf^+
|
|---|
| 1894 |
|
|---|
| 1895 | do j=1,ndg_facex
|
|---|
| 1896 | i=dg_face(j)
|
|---|
| 1897 | du(i) = du(i) + ( us(j)-upwind_wgt(j)*uf(j) )*binvdg(i,1)
|
|---|
| 1898 | enddo
|
|---|
| 1899 |
|
|---|
| 1900 | return
|
|---|
| 1901 | end
|
|---|
| 1902 | c-----------------------------------------------------------------------
|
|---|
