KTH framework for Nek5000 toolboxes; testing version  0.0.1
hpf.f
Go to the documentation of this file.
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
ident
subroutine ident(a, n)
Definition: calcz.f:147
build_hpf_fld
subroutine build_hpf_fld(v, u, f, nx, nz)
Definition: hpf.f:179
hpf_trns_fcn
subroutine hpf_trns_fcn(diag, kut)
Definition: hpf.f:247
build_hpf_mat
subroutine build_hpf_mat(op_mat, f_filter, ifboyd)
Definition: hpf.f:130
make_hpf
subroutine make_hpf
Definition: hpf.f:2
spec_coeff_init
subroutine spec_coeff_init(ref_xmap, ifboyd)
Definition: hpf.f:284
transpose
subroutine transpose(a, lda, b, ldb)
Definition: math.f:1594
addcol3
subroutine addcol3(a, b, c, n)
Definition: math.f:654
copy
subroutine copy(a, b, n)
Definition: math.f:260
sub3
subroutine sub3(a, b, c, n)
Definition: math.f:175
cmult
subroutine cmult(a, const, n)
Definition: math.f:315
mxm
subroutine mxm(a, n1, b, n2, c, n3)
Definition: mxm_wrapper.f:2
opadd2col
subroutine opadd2col(a1, a2, a3, b1, b2, b3, c)
Definition: navier1.f:2751
gaujordf
subroutine gaujordf(a, m, n, indr, indc, ipiv, ierr, rmult)
Definition: navier5.f:668
legendre_poly
subroutine legendre_poly(L, x, N)
Definition: navier5.f:785