| 1 | c-----------------------------------------------------------------------
|
|---|
| 2 | c
|
|---|
| 3 | c Stability limits:
|
|---|
| 4 | c
|
|---|
| 5 | c AB3: .7236 w/safety (1.2): .603
|
|---|
| 6 | c
|
|---|
| 7 | c RK3: 1.73 (sqrt 3) w/safety (1.2): 1.44
|
|---|
| 8 | c
|
|---|
| 9 | c RK4: 2.828 w/safety (1.2): 2.36
|
|---|
| 10 | c
|
|---|
| 11 | c SEM Safety factor: 1.52 for N=3
|
|---|
| 12 | c < 1.20 for N=16
|
|---|
| 13 | c ~ 1.16 for N=256
|
|---|
| 14 | c
|
|---|
| 15 | c-----------------------------------------------------------------------
|
|---|
| 16 | subroutine setup_convect(igeom)
|
|---|
| 17 | include 'SIZE'
|
|---|
| 18 | include 'TOTAL'
|
|---|
| 19 | logical ifnew
|
|---|
| 20 |
|
|---|
| 21 | common /cchar/ ct_vx(0:lorder+1) ! time for each slice in c_vx()
|
|---|
| 22 |
|
|---|
| 23 | common /scruz/ cx (lx1*ly1*lz1*lelt)
|
|---|
| 24 | $ , cy (lx1*ly1*lz1*lelt)
|
|---|
| 25 | $ , cz (lx1*ly1*lz1*lelt)
|
|---|
| 26 | $ , hmsk(lx1*ly1*lz1*lelt)
|
|---|
| 27 |
|
|---|
| 28 |
|
|---|
| 29 | if (igeom.eq.1) return
|
|---|
| 30 | if (param(99).lt.0) return ! no dealiasing
|
|---|
| 31 |
|
|---|
| 32 | if (ifchar) then
|
|---|
| 33 |
|
|---|
| 34 | nelc = nelv
|
|---|
| 35 | if (ifmhd) nelc = max(nelv,nelfld(ifldmhd))
|
|---|
| 36 | if (ifmhd) call exitti('no characteristics for mhd yet$',istep)
|
|---|
| 37 |
|
|---|
| 38 | ifnew = .true.
|
|---|
| 39 | if (igeom.gt.2) ifnew = .false.
|
|---|
| 40 |
|
|---|
| 41 | if (ifgeom) then ! Moving mesh
|
|---|
| 42 | call opsub3(cx,cy,cz,vx,vy,vz,wx,wy,wz)
|
|---|
| 43 | call set_conv_char(ct_vx,c_vx,cx,cy,cz,nelc,time,ifnew)
|
|---|
| 44 | c call set_char_mask(hmsk,cx,cy,cz) ! mask for hyperbolic system
|
|---|
| 45 | else
|
|---|
| 46 | call set_conv_char(ct_vx,c_vx,vx,vy,vz,nelc,time,ifnew)
|
|---|
| 47 | c call set_char_mask(hmsk,vx,vy,vz) ! mask for hyperbolic system
|
|---|
| 48 | endif
|
|---|
| 49 |
|
|---|
| 50 | n=lx1*ly1*lz1*nelv
|
|---|
| 51 | call rone(hmsk,n) ! TURN OFF MASK
|
|---|
| 52 | call set_binv (bmnv, hmsk,n) ! Store binvm1*(hyperbolic mask)
|
|---|
| 53 | if(ifmvbd) then
|
|---|
| 54 | call set_bdivw(bdivw,hmsk,n) ! Store Bdivw *(hyperbolic mask)
|
|---|
| 55 | else
|
|---|
| 56 | call rzero(bdivw,n*lorder)
|
|---|
| 57 | endif
|
|---|
| 58 | call set_bmass(bmass,hmsk,n) ! Store binvm1*(hyperbolic mask)
|
|---|
| 59 |
|
|---|
| 60 | else
|
|---|
| 61 |
|
|---|
| 62 | if (.not.ifpert) then
|
|---|
| 63 | if (ifcons) then
|
|---|
| 64 | call set_convect_cons (vxd,vyd,vzd,vx,vy,vz)
|
|---|
| 65 | else
|
|---|
| 66 | call set_convect_new (vxd,vyd,vzd,vx,vy,vz)
|
|---|
| 67 | endif
|
|---|
| 68 | endif
|
|---|
| 69 |
|
|---|
| 70 | endif
|
|---|
| 71 |
|
|---|
| 72 | return
|
|---|
| 73 | end
|
|---|
| 74 | c-----------------------------------------------------------------------
|
|---|
| 75 | subroutine set_char_mask(mask,u,v,w) ! mask for hyperbolic system
|
|---|
| 76 |
|
|---|
| 77 | include 'SIZE'
|
|---|
| 78 | include 'INPUT'
|
|---|
| 79 | include 'GEOM'
|
|---|
| 80 | include 'SOLN'
|
|---|
| 81 | include 'TSTEP'
|
|---|
| 82 | integer msk(0:1)
|
|---|
| 83 | character cb*3
|
|---|
| 84 | parameter (lxyz1=lx1*ly1*lz1)
|
|---|
| 85 | common /ctmp1/ work(lxyz1,lelt)
|
|---|
| 86 |
|
|---|
| 87 | real mask(lxyz1,1),u(lxyz1,1),v(lxyz1,1),w(lxyz1,1)
|
|---|
| 88 |
|
|---|
| 89 | integer e,f
|
|---|
| 90 |
|
|---|
| 91 | nfaces= 2*ldim
|
|---|
| 92 | ntot1 = lx1*ly1*lz1*nelv
|
|---|
| 93 | call rzero (work,ntot1)
|
|---|
| 94 | call rone (mask,NTOT1)
|
|---|
| 95 |
|
|---|
| 96 | ifldv = 1
|
|---|
| 97 | do 100 e=1,nelv
|
|---|
| 98 | do 100 f=1,nfaces
|
|---|
| 99 | cb=cbc(f,e,ifldv)
|
|---|
| 100 | if (cb(1:1).eq.'v' .or. cb(1:1).eq.'V' .or.
|
|---|
| 101 | $ cb.eq.'mv ' .or. cb.eq.'MV ') then
|
|---|
| 102 |
|
|---|
| 103 | call faccl3 (work(1,e),u(1,e),unx(1,1,f,e),f)
|
|---|
| 104 | call faddcl3(work(1,e),v(1,e),uny(1,1,f,e),f)
|
|---|
| 105 | if (if3d)
|
|---|
| 106 | $ call faddcl3(work(1,e),w(1,e),unz(1,1,f,e),f)
|
|---|
| 107 |
|
|---|
| 108 | call fcaver (vaver,work,e,f)
|
|---|
| 109 |
|
|---|
| 110 | if (vaver.lt.0) call facev (mask,e,f,0.0,lx1,ly1,lz1)
|
|---|
| 111 | endif
|
|---|
| 112 | if (cb(1:2).eq.'ws' .or. cb(1:2).eq.'WS')
|
|---|
| 113 | $ call facev (mask,e,f,0.0,lx1,ly1,lz1)
|
|---|
| 114 | 100 continue
|
|---|
| 115 | call dsop(mask,'MUL',lx1,ly1,lz1)
|
|---|
| 116 |
|
|---|
| 117 | return
|
|---|
| 118 | end
|
|---|
| 119 | c-----------------------------------------------------------------------
|
|---|
| 120 | subroutine set_dealias_rx
|
|---|
| 121 | C
|
|---|
| 122 | C Eulerian scheme, add convection term to forcing function
|
|---|
| 123 | C at current time step.
|
|---|
| 124 | C
|
|---|
| 125 | include 'SIZE'
|
|---|
| 126 | include 'INPUT'
|
|---|
| 127 | include 'GEOM'
|
|---|
| 128 | include 'TSTEP' ! for istep
|
|---|
| 129 |
|
|---|
| 130 | common /dealias1/ zd(lxd),wd(lxd)
|
|---|
| 131 | integer e
|
|---|
| 132 |
|
|---|
| 133 | integer ilstep
|
|---|
| 134 | save ilstep
|
|---|
| 135 | data ilstep /-1/
|
|---|
| 136 |
|
|---|
| 137 | if (.not.ifgeom.and.ilstep.gt.1) return ! already computed
|
|---|
| 138 | if (ifgeom.and.ilstep.eq.istep) return ! already computed
|
|---|
| 139 | ilstep = istep
|
|---|
| 140 |
|
|---|
| 141 | nxyz1 = lx1*ly1*lz1
|
|---|
| 142 | nxyzd = lxd*lyd*lzd
|
|---|
| 143 |
|
|---|
| 144 | call zwgl (zd,wd,lxd) ! zwgl -- NOT zwgll!
|
|---|
| 145 |
|
|---|
| 146 | if (if3d) then
|
|---|
| 147 | c
|
|---|
| 148 | do e=1,nelv
|
|---|
| 149 |
|
|---|
| 150 | c Interpolate z+ and z- into fine mesh, translate to r-s-t coords
|
|---|
| 151 |
|
|---|
| 152 | call intp_rstd(rx(1,1,e),rxm1(1,1,1,e),lx1,lxd,if3d,0) ! 0 --> fwd
|
|---|
| 153 | call intp_rstd(rx(1,2,e),rym1(1,1,1,e),lx1,lxd,if3d,0) ! 0 --> fwd
|
|---|
| 154 | call intp_rstd(rx(1,3,e),rzm1(1,1,1,e),lx1,lxd,if3d,0) ! 0 --> fwd
|
|---|
| 155 | call intp_rstd(rx(1,4,e),sxm1(1,1,1,e),lx1,lxd,if3d,0) ! 0 --> fwd
|
|---|
| 156 | call intp_rstd(rx(1,5,e),sym1(1,1,1,e),lx1,lxd,if3d,0) ! 0 --> fwd
|
|---|
| 157 | call intp_rstd(rx(1,6,e),szm1(1,1,1,e),lx1,lxd,if3d,0) ! 0 --> fwd
|
|---|
| 158 | call intp_rstd(rx(1,7,e),txm1(1,1,1,e),lx1,lxd,if3d,0) ! 0 --> fwd
|
|---|
| 159 | call intp_rstd(rx(1,8,e),tym1(1,1,1,e),lx1,lxd,if3d,0) ! 0 --> fwd
|
|---|
| 160 | call intp_rstd(rx(1,9,e),tzm1(1,1,1,e),lx1,lxd,if3d,0) ! 0 --> fwd
|
|---|
| 161 |
|
|---|
| 162 | l = 0
|
|---|
| 163 | do k=1,lzd
|
|---|
| 164 | do j=1,lyd
|
|---|
| 165 | do i=1,lxd
|
|---|
| 166 | l = l+1
|
|---|
| 167 | w = wd(i)*wd(j)*wd(k)
|
|---|
| 168 | do ii=1,9
|
|---|
| 169 | rx(l,ii,e) = w*rx(l,ii,e)
|
|---|
| 170 | enddo
|
|---|
| 171 | enddo
|
|---|
| 172 | enddo
|
|---|
| 173 | enddo
|
|---|
| 174 | enddo
|
|---|
| 175 |
|
|---|
| 176 | else ! 2D
|
|---|
| 177 | c
|
|---|
| 178 | do e=1,nelv
|
|---|
| 179 |
|
|---|
| 180 | c Interpolate z+ and z- into fine mesh, translate to r-s-t coords
|
|---|
| 181 |
|
|---|
| 182 | call intp_rstd(rx(1,1,e),rxm1(1,1,1,e),lx1,lxd,if3d,0) ! 0 --> fwd
|
|---|
| 183 | call intp_rstd(rx(1,2,e),rym1(1,1,1,e),lx1,lxd,if3d,0) ! 0 --> fwd
|
|---|
| 184 | call intp_rstd(rx(1,3,e),sxm1(1,1,1,e),lx1,lxd,if3d,0) ! 0 --> fwd
|
|---|
| 185 | call intp_rstd(rx(1,4,e),sym1(1,1,1,e),lx1,lxd,if3d,0) ! 0 --> fwd
|
|---|
| 186 |
|
|---|
| 187 | l = 0
|
|---|
| 188 | do j=1,lyd
|
|---|
| 189 | do i=1,lxd
|
|---|
| 190 | l = l+1
|
|---|
| 191 | w = wd(i)*wd(j)
|
|---|
| 192 | do ii=1,4
|
|---|
| 193 | rx(l,ii,e) = w*rx(l,ii,e)
|
|---|
| 194 | enddo
|
|---|
| 195 | enddo
|
|---|
| 196 | enddo
|
|---|
| 197 | enddo
|
|---|
| 198 |
|
|---|
| 199 | endif
|
|---|
| 200 |
|
|---|
| 201 | return
|
|---|
| 202 | end
|
|---|
| 203 | c-----------------------------------------------------------------------
|
|---|
| 204 |
|
|---|
