| 1 | subroutine MAKE_HPF
|
|---|
| 2 |
|
|---|
| 3 | include 'SIZE'
|
|---|
| 4 | include 'SOLN'
|
|---|
| 5 | include 'INPUT' ! param(110),(111)
|
|---|
| 6 | include 'TSTEP' ! ifield
|
|---|
| 7 | include 'MASS' ! BM1
|
|---|
| 8 |
|
|---|
| 9 | integer nxyz
|
|---|
| 10 | parameter (nxyz=lx1*ly1*lz1)
|
|---|
| 11 | integer n
|
|---|
| 12 |
|
|---|
| 13 | integer lm,lm2
|
|---|
| 14 | parameter (lm=40)
|
|---|
| 15 | parameter (lm2=lm*lm)
|
|---|
| 16 | real hpf_filter(lm2)
|
|---|
| 17 |
|
|---|
| 18 | integer hpf_kut
|
|---|
| 19 | real hpf_chi
|
|---|
| 20 | logical hpf_ifboyd
|
|---|
| 21 |
|
|---|
| 22 | integer nel
|
|---|
| 23 |
|
|---|
| 24 | integer icalld
|
|---|
| 25 | save icalld
|
|---|
| 26 | data icalld /0/
|
|---|
| 27 |
|
|---|
| 28 | real TA1,TA2,TA3
|
|---|
| 29 | COMMON /SCRUZ/ TA1 (LX1,LY1,LZ1,LELV)
|
|---|
| 30 | $ , TA2 (LX1,LY1,LZ1,LELV)
|
|---|
| 31 | $ , TA3 (LX1,LY1,LZ1,LELV)
|
|---|
| 32 |
|
|---|
| 33 | real hpf_op(lx1,lx1)
|
|---|
| 34 | save hpf_op
|
|---|
| 35 |
|
|---|
| 36 | c----------------------------------------
|
|---|
| 37 | if(.not. iffilter(ifield)) return
|
|---|
| 38 |
|
|---|
| 39 | hpf_kut = int(param(101))+1
|
|---|
| 40 | hpf_chi = -1.0*abs(param(103))
|
|---|
| 41 | c Boyd transform to preserve element boundary values is
|
|---|
| 42 | c linearly unstable when used as forcing.
|
|---|
| 43 | hpf_ifboyd = .false.
|
|---|
| 44 |
|
|---|
| 45 | nel = nelv
|
|---|
| 46 | n = nxyz*nel
|
|---|
| 47 |
|
|---|
| 48 | if (hpf_chi.eq.0) return
|
|---|
| 49 | if(nid.eq.0 .and. loglevel.gt.2) write(6,*) 'apply hpf ',
|
|---|
| 50 | $ ifield, hpf_kut, hpf_chi
|
|---|
| 51 |
|
|---|
| 52 | if (icalld.eq.0) then
|
|---|
| 53 | c Create the filter transfer function
|
|---|
| 54 | call hpf_trns_fcn(hpf_filter,hpf_kut)
|
|---|
| 55 |
|
|---|
| 56 | c Build the matrix to apply the filter function
|
|---|
| 57 | c to an input field
|
|---|
| 58 | call build_hpf_mat(hpf_op,hpf_filter,hpf_ifboyd)
|
|---|
| 59 |
|
|---|
| 60 | c Only initialize once
|
|---|
| 61 | icalld=icalld+1
|
|---|
| 62 | endif
|
|---|
| 63 |
|
|---|
| 64 | if (ifield.eq.1) then
|
|---|
| 65 | c Apply the filter
|
|---|
| 66 | c to velocity fields
|
|---|
| 67 | if (jp.eq.0) then
|
|---|
| 68 | call build_hpf_fld(ta1,vx,hpf_op,lx1,lz1)
|
|---|
| 69 | call build_hpf_fld(ta2,vy,hpf_op,lx1,lz1)
|
|---|
| 70 | if (if3d) call build_hpf_fld(ta3,vz,hpf_op,lx1,lz1)
|
|---|
| 71 |
|
|---|
| 72 | c Multiply by filter weight (chi)
|
|---|
| 73 | call cmult(ta1,hpf_chi,n)
|
|---|
| 74 | call cmult(ta2,hpf_chi,n)
|
|---|
| 75 | if (if3d) call cmult(ta3,hpf_chi,n)
|
|---|
| 76 |
|
|---|
| 77 | c Multiply by Mass matrix
|
|---|
| 78 | c and add to forcing term
|
|---|
| 79 | call opadd2col (bfx,bfy,bfz,ta1,ta2,ta3,bm1)
|
|---|
| 80 | else
|
|---|
| 81 | c Apply filter on velocity perturbation fields
|
|---|
| 82 | call build_hpf_fld(ta1,vxp(1,jp),hpf_op,lx1,lz1)
|
|---|
| 83 | call build_hpf_fld(ta2,vyp(1,jp),hpf_op,lx1,lz1)
|
|---|
| 84 | if (if3d) call build_hpf_fld(ta3,vzp(1,jp),hpf_op,lx1,lz1)
|
|---|
| 85 |
|
|---|
| 86 | c Multiply by filter weight (chi)
|
|---|
| 87 | call cmult(ta1,hpf_chi,n)
|
|---|
| 88 | call cmult(ta2,hpf_chi,n)
|
|---|
| 89 | if (if3d) call cmult(ta3,hpf_chi,n)
|
|---|
| 90 |
|
|---|
| 91 | c Multiply by Mass matrix
|
|---|
| 92 | c and add to forcing term
|
|---|
| 93 | call opadd2col (bfxp(1,jp),bfyp(1,jp),bfzp(1,jp),
|
|---|
| 94 | $ ta1,ta2,ta3,bm1)
|
|---|
| 95 | endif
|
|---|
| 96 |
|
|---|
| 97 | else
|
|---|
| 98 |
|
|---|
| 99 | c Apply filter to temp/passive scalar fields
|
|---|
| 100 | if (jp.eq.0) then
|
|---|
| 101 | call build_hpf_fld(ta1,t(1,1,1,1,ifield-1),
|
|---|
| 102 | $ hpf_op,lx1,lz1)
|
|---|
| 103 |
|
|---|
| 104 | c Multiply by filter weight (chi)
|
|---|
| 105 | call cmult(ta1,hpf_chi,n)
|
|---|
| 106 |
|
|---|
| 107 | c Multiply by Mass matrix
|
|---|
| 108 | c and add to source term
|
|---|
| 109 | call addcol3(bq(1,1,1,1,ifield-1),ta1,bm1,n)
|
|---|
| 110 | else
|
|---|
| 111 | c Apply filter on scalar perturbation field
|
|---|
| 112 | call build_hpf_fld(ta1,tp(1,ifield-1,jp),hpf_op,lx1,lz1)
|
|---|
| 113 |
|
|---|
| 114 | c Multiply by filter weight (chi)
|
|---|
| 115 | call cmult(ta1,hpf_chi,n)
|
|---|
| 116 |
|
|---|
| 117 | c Multiply by Mass matrix
|
|---|
| 118 | c and add to source term
|
|---|
| 119 | call addcol3(bqp(1,ifield-1,jp),ta1,bm1,n)
|
|---|
| 120 | endif
|
|---|
| 121 |
|
|---|
| 122 | endif
|
|---|
| 123 |
|
|---|
| 124 | return
|
|---|
| 125 | end
|
|---|
| 126 |
|
|---|
| 127 | c----------------------------------------------------------------------
|
|---|
| 128 |
|
|---|
| 129 | subroutine build_hpf_mat(op_mat,f_filter,ifboyd)
|
|---|
| 130 |
|
|---|
| 131 | c Builds the operator for high pass filtering
|
|---|
| 132 | c Transformation matrix from nodal to modal space.
|
|---|
| 133 | c Applies f_filter to the the legendre coefficients
|
|---|
| 134 | c Transforms back to nodal space
|
|---|
| 135 | c Operation: V * f_filter * V^(-1)
|
|---|
| 136 | c Where V is the transformation matrix from modal to nodal space
|
|---|
| 137 |
|
|---|
| 138 | c implicit none
|
|---|
| 139 |
|
|---|
| 140 | include 'SIZE'
|
|---|
| 141 |
|
|---|
| 142 | logical IFBOYD
|
|---|
| 143 | integer n
|
|---|
| 144 | parameter (n=lx1*lx1)
|
|---|
| 145 | integer lm, lm2
|
|---|
| 146 | parameter (lm=40)
|
|---|
| 147 | parameter (lm2=lm*lm)
|
|---|
| 148 |
|
|---|
| 149 | real f_filter(lm2)
|
|---|
| 150 | real op_mat(lx1,lx1)
|
|---|
| 151 |
|
|---|
| 152 | real ref_xmap(lm2)
|
|---|
| 153 | real wk_xmap(lm2)
|
|---|
| 154 |
|
|---|
| 155 | real wk1(lm2),wk2(lm2)
|
|---|
| 156 | real indr(lm),ipiv(lm),indc(lm)
|
|---|
| 157 |
|
|---|
| 158 | real rmult(lm)
|
|---|
| 159 | integer ierr
|
|---|
| 160 |
|
|---|
| 161 | integer i,j
|
|---|
| 162 |
|
|---|
| 163 | call spec_coeff_init(ref_xmap,ifboyd)
|
|---|
| 164 |
|
|---|
| 165 | call copy(wk_xmap,ref_xmap,lm2)
|
|---|
| 166 | call copy(wk1,wk_XMAP,lm2)
|
|---|
| 167 |
|
|---|
| 168 | call gaujordf (wk1,lx1,lx1,indr,indc,ipiv,ierr,rmult) ! xmap inverse
|
|---|
| 169 |
|
|---|
| 170 | call mxm (f_filter,lx1,wk1,lx1,wk2,lx1) ! -1
|
|---|
| 171 | call mxm (wk_xmap,lx1,wk2,lx1,op_mat,lx1) ! V D V
|
|---|
| 172 |
|
|---|
| 173 | return
|
|---|
| 174 | end
|
|---|
| 175 |
|
|---|
| 176 | c----------------------------------------------------------------------
|
|---|
| 177 |
|
|---|
| 178 | subroutine build_hpf_fld(v,u,f,nx,nz)
|
|---|
| 179 |
|
|---|
| 180 | c Appies the operator f to field u
|
|---|
| 181 | c using tensor operations
|
|---|
| 182 | c v = f*u
|
|---|
| 183 |
|
|---|
| 184 | c implicit none
|
|---|
| 185 |
|
|---|
| 186 | include 'SIZE'
|
|---|
| 187 | include 'TSTEP' ! ifield
|
|---|
| 188 |
|
|---|
| 189 | integer nxyz
|
|---|
| 190 | parameter (nxyz=lx1*ly1*lz1)
|
|---|
| 191 |
|
|---|
| 192 | real w1(nxyz),w2(nxyz) ! work arrays
|
|---|
| 193 | real v(nxyz,lelv) ! output fld
|
|---|
| 194 | real u(nxyz,lelv) ! input fld
|
|---|
| 195 | c
|
|---|
| 196 | integer nx,nz
|
|---|
| 197 |
|
|---|
| 198 | real f(nx,nx),ft(nx,nx) ! operator f and its transpose
|
|---|
| 199 | c
|
|---|
| 200 | integer e,i,j,k
|
|---|
| 201 | integer nel
|
|---|
| 202 |
|
|---|
| 203 | nel = nelv
|
|---|
| 204 |
|
|---|
| 205 | call copy(v,u,nxyz*nel)
|
|---|
| 206 | c
|
|---|
| 207 | call transpose(ft,nx,f,nx)
|
|---|
| 208 | c
|
|---|
| 209 | if (ldim.eq.3) then
|
|---|
| 210 | do e=1,nel
|
|---|
| 211 | c Filter
|
|---|
| 212 | call copy(w2,v(1,e),nxyz)
|
|---|
| 213 | call mxm(f,nx,w2,nx,w1,nx*nx)
|
|---|
| 214 | i=1
|
|---|
| 215 | j=1
|
|---|
| 216 | do k=1,nx
|
|---|
| 217 | call mxm(w1(i),nx,ft,nx,w2(j),nx)
|
|---|
| 218 | i = i+nx*nx
|
|---|
| 219 | j = j+nx*nx
|
|---|
| 220 | enddo
|
|---|
| 221 | call mxm (w2,nx*nx,ft,nx,w1,nx)
|
|---|
| 222 |
|
|---|
| 223 | call sub3(w2,u(1,e),w1,nxyz)
|
|---|
| 224 | call copy(v(1,e),w2,nxyz)
|
|---|
| 225 | c call copy(v(1,e),w1,nxyz)
|
|---|
| 226 |
|
|---|
| 227 | enddo
|
|---|
| 228 | else
|
|---|
| 229 | do e=1,nel
|
|---|
| 230 | c Filter
|
|---|
| 231 | call copy(w1,v(1,e),nxyz)
|
|---|
| 232 | call mxm(f ,nx,w1,nx,w2,nx)
|
|---|
| 233 | call mxm(w2,nx,ft,nx,w1,nx)
|
|---|
| 234 |
|
|---|
| 235 | call sub3(w2,u(1,e),w1,nxyz)
|
|---|
| 236 | call copy(v(1,e),w2,nxyz)
|
|---|
| 237 | c call copy(v(1,e),w1,nxyz)
|
|---|
| 238 | enddo
|
|---|
| 239 | endif
|
|---|
| 240 | c
|
|---|
| 241 | return
|
|---|
| 242 | end
|
|---|
| 243 |
|
|---|
| 244 | c----------------------------------------------------------------------
|
|---|
| 245 |
|
|---|
| 246 | subroutine hpf_trns_fcn(diag,kut)
|
|---|
| 247 |
|
|---|
| 248 | c implicit none
|
|---|
| 249 |
|
|---|
| 250 | include 'SIZE'
|
|---|
| 251 | include 'PARALLEL'
|
|---|
| 252 |
|
|---|
| 253 | real diag(lx1*lx1)
|
|---|
| 254 | integer nx,k0,kut,kk,k
|
|---|
| 255 |
|
|---|
| 256 | real amp
|
|---|
| 257 |
|
|---|
| 258 | c Set up transfer function
|
|---|
| 259 | c
|
|---|
| 260 | nx = lx1
|
|---|
| 261 | call ident (diag,nx)
|
|---|
| 262 | c
|
|---|
| 263 | kut=kut ! kut=additional modes
|
|---|
| 264 | k0 = nx-kut
|
|---|
| 265 | do k=k0+1,nx
|
|---|
| 266 | kk = k+nx*(k-1)
|
|---|
| 267 | amp = ((k-k0)*(k-k0)+0.)/(kut*kut+0.) ! Normalized amplitude. quadratic growth
|
|---|
| 268 | diag(kk) = 1.-amp
|
|---|
| 269 | enddo
|
|---|
| 270 |
|
|---|
| 271 | c Output normalized transfer function
|
|---|
| 272 | k0 = lx1+1
|
|---|
| 273 | if (nio.eq.0) then
|
|---|
| 274 | write(6,6) 'HPF :',((1.-diag(k)), k=1,lx1*lx1,k0)
|
|---|
| 275 | 6 format(a8,16f9.6,6(/,8x,16f9.6))
|
|---|
| 276 | endif
|
|---|
| 277 |
|
|---|
| 278 | return
|
|---|
| 279 | end
|
|---|
| 280 |
|
|---|
| 281 | c----------------------------------------------------------------------
|
|---|
| 282 |
|
|---|
| 283 | subroutine spec_coeff_init(ref_xmap,ifboyd)
|
|---|
| 284 | c Initialise spectral coefficients
|
|---|
| 285 | c For legendre transform
|
|---|
| 286 |
|
|---|
| 287 | c implicit none
|
|---|
| 288 |
|
|---|
| 289 | include 'SIZE'
|
|---|
| 290 | include 'WZ'
|
|---|
| 291 |
|
|---|
| 292 | integer lm, lm2
|
|---|
| 293 | parameter (lm=40)
|
|---|
| 294 | parameter (lm2=lm*lm)
|
|---|
| 295 |
|
|---|
| 296 | c local variables
|
|---|
| 297 | integer i, j, k, n, nx, kj
|
|---|
| 298 | c Legendre polynomial
|
|---|
| 299 | real plegx(lm)
|
|---|
| 300 | real z
|
|---|
| 301 | real ref_xmap(lm2)
|
|---|
| 302 | real pht(lm2)
|
|---|
| 303 |
|
|---|
| 304 | c Change of basis
|
|---|
| 305 | logical IFBOYD
|
|---|
| 306 | c----------------------------------------
|
|---|
| 307 |
|
|---|
| 308 | nx = LX1
|
|---|
| 309 | kj = 0
|
|---|
| 310 | n = nx-1
|
|---|
| 311 | do j=1,nx
|
|---|
| 312 | z = ZGM1(j,1)
|
|---|
| 313 | call legendre_poly(plegx,z,n)
|
|---|
| 314 | kj = kj+1
|
|---|
| 315 | pht(kj) = plegx(1)
|
|---|
| 316 | kj = kj+1
|
|---|
| 317 | pht(kj) = plegx(2)
|
|---|
| 318 |
|
|---|
| 319 | if (IFBOYD) then ! change basis to preserve element boundary values
|
|---|
| 320 | do k=3,nx
|
|---|
| 321 | kj = kj+1
|
|---|
| 322 | pht(kj) = plegx(k)-plegx(k-2)
|
|---|
| 323 | enddo
|
|---|
| 324 | else ! legendre basis
|
|---|
| 325 | do k=3,nx
|
|---|
| 326 | kj = kj+1
|
|---|
| 327 | pht(kj) = plegx(k)
|
|---|
| 328 | enddo
|
|---|
| 329 | endif
|
|---|
| 330 | enddo
|
|---|
| 331 |
|
|---|
| 332 | call transpose (ref_xmap,nx,pht,nx)
|
|---|
| 333 |
|
|---|
| 334 | return
|
|---|
| 335 | end
|
|---|
