| 1 | real function avm_vdiff(ix,iy,iz,e,c1,ncut)
|
|---|
| 2 | c
|
|---|
| 3 | c c1 and ncut a user tuneable control parameters.
|
|---|
| 4 | c Set c1 = 1.0 and reduce/increase as much possible/required,
|
|---|
| 5 | c depending on your application.
|
|---|
| 6 | c Typically ncut 1 or 2 works well.
|
|---|
| 7 | c Note, avoid using lx1 < 6!
|
|---|
| 8 | c
|
|---|
| 9 | include 'SIZE'
|
|---|
| 10 | include 'TOTAL'
|
|---|
| 11 |
|
|---|
| 12 | integer ix, iy, iz, e
|
|---|
| 13 | real c1, c2
|
|---|
| 14 | integer ncut
|
|---|
| 15 |
|
|---|
| 16 | logical ifcont
|
|---|
| 17 | parameter (ifcont=.false.)
|
|---|
| 18 |
|
|---|
| 19 | real visc(lx1,ly1,lz1,lelt)
|
|---|
| 20 | save visc
|
|---|
| 21 |
|
|---|
| 22 | common /SCRMG/ r (lx1*ly1*lz1,lelt),
|
|---|
| 23 | $ tx(lx1*ly1*lz1,lelt),
|
|---|
| 24 | $ ty(lx1*ly1*lz1,lelt),
|
|---|
| 25 | $ tz(lx1*ly1*lz1,lelt)
|
|---|
| 26 |
|
|---|
| 27 | parameter (lm=40)
|
|---|
| 28 | parameter (lm2=lm*lm)
|
|---|
| 29 | real hpf_filter(lm2)
|
|---|
| 30 | real hpf_op(lx1*lx1,ldimt1)
|
|---|
| 31 | save hpf_op
|
|---|
| 32 |
|
|---|
| 33 | integer ibuild(ldimt1)
|
|---|
| 34 | save ibuild
|
|---|
| 35 |
|
|---|
| 36 | save icalld
|
|---|
| 37 | data icalld / 0 /
|
|---|
| 38 |
|
|---|
| 39 | real h0,h0max
|
|---|
| 40 | real viscc(8,lelt)
|
|---|
| 41 |
|
|---|
| 42 | if (ix*iy*iz*e .ne. 1) then ! use cache
|
|---|
| 43 | avm_vdiff = visc(ix,iy,iz,e)
|
|---|
| 44 | return
|
|---|
| 45 | endif
|
|---|
| 46 |
|
|---|
| 47 | if (icalld.eq.0) then
|
|---|
| 48 | iffilter(ifield) = .false.
|
|---|
| 49 | do i = 1,ldimt1
|
|---|
| 50 | ibuild(i) = 0
|
|---|
| 51 | enddo
|
|---|
| 52 | icalld = 1
|
|---|
| 53 | endif
|
|---|
| 54 |
|
|---|
| 55 | nxyz = lx1*ly1*lz1
|
|---|
| 56 | n = nxyz*nelv
|
|---|
| 57 | c2 = 0.5
|
|---|
| 58 |
|
|---|
| 59 | ! compute residual
|
|---|
| 60 | if (ifield.eq.1) then
|
|---|
| 61 | if (nid.eq.0) write(6,*) 'avm not supported for ifield=1 !'
|
|---|
| 62 | call exitt
|
|---|
| 63 | else
|
|---|
| 64 | if (ibuild(ifield).eq.0) then
|
|---|
| 65 | call hpf_trns_fcn(hpf_filter,ncut)
|
|---|
| 66 | call build_hpf_mat(hpf_op(1,ifield),hpf_filter,.false.)
|
|---|
| 67 | ibuild(ifield) = ibuild(ifield) + 1
|
|---|
| 68 | endif
|
|---|
| 69 | call build_hpf_fld(r,t(1,1,1,1,ifield-1),hpf_op(1,ifield),
|
|---|
| 70 | $ lx1,lz1)
|
|---|
| 71 |
|
|---|
| 72 | psave = param(99)
|
|---|
| 73 | param(99) = 0
|
|---|
| 74 | call copy(tx,r,n)
|
|---|
| 75 | call convop(r,tx)
|
|---|
| 76 | param(99) = psave
|
|---|
| 77 |
|
|---|
| 78 | ! normalize
|
|---|
| 79 | uavg = gl2norm(t(1,1,1,1,ifield-1),n)
|
|---|
| 80 | call cfill(tx,uavg,n)
|
|---|
| 81 | call sub2(tx,t(1,1,1,1,ifield-1),n)
|
|---|
| 82 | uinf = 1.
|
|---|
| 83 | if(uavg.gt.0) uinf = glamax(tx, n)
|
|---|
| 84 | endif
|
|---|
| 85 |
|
|---|
| 86 | ! evaluate arificial viscosity
|
|---|
| 87 | uinf = 1./uinf
|
|---|
| 88 | do ie = 1,nelv
|
|---|
| 89 | h0max = vlmax(deltaf(1,1,1,ie),nxyz)
|
|---|
| 90 | do i = 1,nxyz
|
|---|
| 91 | h0 = deltaf(i,1,1,ie)
|
|---|
| 92 | vmax = sqrt(vx(i,1,1,ie)**2 + vy(i,1,1,ie)**2
|
|---|
| 93 | $ + vz(i,1,1,ie)**2)
|
|---|
| 94 | vismax = c2 * h0max * vmax
|
|---|
| 95 | visc(i,1,1,ie) = min(vismax, c1*h0**2 * abs(r(i,ie))*uinf)
|
|---|
| 96 | enddo
|
|---|
| 97 | enddo
|
|---|
| 98 |
|
|---|
| 99 | ! make it piecewise constant
|
|---|
| 100 | do ie = 1,nelv
|
|---|
| 101 | vmax = vlmax(visc(1,1,1,ie),nxyz)
|
|---|
| 102 | call cfill(visc(1,1,1,ie),vmax,nxyz)
|
|---|
| 103 | enddo
|
|---|
| 104 |
|
|---|
| 105 | ! make it P1 continuous
|
|---|
| 106 | if (ifcont) then
|
|---|
| 107 | call dsop (visc,'max',lx1,ly1,lz1)
|
|---|
| 108 | do ie = 1,nelv
|
|---|
| 109 | viscc(1,ie) = visc(1 ,1 ,1 ,ie)
|
|---|
| 110 | viscc(2,ie) = visc(lx1,1 ,1 ,ie)
|
|---|
| 111 | viscc(3,ie) = visc(1 ,ly1,1 ,ie)
|
|---|
| 112 | viscc(4,ie) = visc(lx1,ly1,1 ,ie)
|
|---|
| 113 |
|
|---|
| 114 | viscc(5,ie) = visc(1 ,1 ,lz1,ie)
|
|---|
| 115 | viscc(6,ie) = visc(lx1,1 ,lz1,ie)
|
|---|
| 116 | viscc(7,ie) = visc(1 ,ly1,lz1,ie)
|
|---|
| 117 | viscc(8,ie) = visc(lx1,ly1,lz1,ie)
|
|---|
| 118 | enddo
|
|---|
| 119 | call map_c_to_f_h1_bilin(visc,viscc)
|
|---|
| 120 | endif
|
|---|
| 121 |
|
|---|
| 122 | if (mod(istep,10).eq.0) then
|
|---|
| 123 | vismx = glamax(visc,n)
|
|---|
| 124 | vismn = glamin(visc,n)
|
|---|
| 125 | visav = glsc2(visc,bm1,n)/volvm1
|
|---|
| 126 | if (nio.eq.0) write(6,10) time,vismx,vismn,visav,ifield
|
|---|
| 127 | 10 format(1p4e12.4,' AVM',i6)
|
|---|
| 128 | endif
|
|---|
| 129 |
|
|---|
| 130 | avm_vdiff = visc(ix,iy,iz,e)
|
|---|
| 131 |
|
|---|
| 132 | return
|
|---|
| 133 | end
|
|---|
| 134 | c---------------------------------------------------------------
|
|---|
| 135 | real function deltaf(ix,iy,iz,iel)
|
|---|
| 136 | c
|
|---|
| 137 | c compute characteristic sgs length scale
|
|---|
| 138 | c defined as min of GLL nodes within an elements but many
|
|---|
| 139 | c others possible.
|
|---|
| 140 | c
|
|---|
| 141 | include 'SIZE'
|
|---|
| 142 | include 'TOTAL'
|
|---|
| 143 |
|
|---|
| 144 | real dx(lx1,ly1,lz1,lelt)
|
|---|
| 145 | save dx
|
|---|
| 146 |
|
|---|
| 147 | data icalld /0/
|
|---|
| 148 | save icalld
|
|---|
| 149 |
|
|---|
| 150 | nxyz = nx1*ny1*nz1
|
|---|
| 151 | n = nxyz*nelv
|
|---|
| 152 |
|
|---|
| 153 | if (icalld.eq.0 .or. ifmvbd) then
|
|---|
| 154 | dinv = 1./ldim
|
|---|
| 155 | do ie = 1,nelv
|
|---|
| 156 | c volavg = 0
|
|---|
| 157 | c do i = 1,nxyz
|
|---|
| 158 | c volavg = volavg + bm1(i,1,1,ie)
|
|---|
| 159 | c enddo
|
|---|
| 160 | c dd = (volavg**dinv)/(lx1-1)
|
|---|
| 161 |
|
|---|
| 162 | dd = dxmin_e(ie)
|
|---|
| 163 |
|
|---|
| 164 | call cfill(dx(1,1,1,ie),dd,nxyz)
|
|---|
| 165 | enddo
|
|---|
| 166 | icalld = 1
|
|---|
| 167 | endif
|
|---|
| 168 |
|
|---|
| 169 | deltaf = dx(ix,iy,iz,iel)
|
|---|
| 170 |
|
|---|
| 171 | end
|
|---|
