KTH framework for Nek5000 toolboxes; testing version  0.0.1
fasts.f
Go to the documentation of this file.
1 c-----------------------------------------------------------------------
2  subroutine local_solves_fdm(u,v)
3 c
4 c Given an input vector v, this returns the additive Schwarz solution
5 c
6 c -1
7 c u = M v
8 c
9 c T -1
10 c where M = sum ( R_i A_i R_i )
11 c i
12 c
13 c The R_i's are simply index_set restriction operators.
14 c
15 c The local solves are performed with the fast diagonalization method.
16 c
17  include 'SIZE'
18  include 'INPUT'
19  include 'DOMAIN'
20  include 'PARALLEL'
21 c
22  include 'TSTEP'
23  include 'CTIMER'
24 c
25  real u(lx2,ly2,lz2,lelv),v(lx2,ly2,lz2,lelv)
26  common /scrpre/ v1(lx1,ly1,lz1,lelv)
27  $ ,w1(lx1,ly1,lz1),w2(lx1,ly1,lz1)
28  common /scrover/ ar(lelv)
29 
30  parameter(lxx=lx1*lx1, levb=lelv+lbelv)
31  common /fastd/ df(lx1*ly1*lz1,levb)
32  $ , sr(lxx*2,levb),ss(lxx*2,levb),st(lxx*2,levb)
33  integer e,eb,eoff
34 
35 c
36  if (icalld.eq.0) tsolv=0.0
37  icalld=icalld+1
38  nsolv=icalld
39 c
40  ntot1 = lx1*ly1*lz1*nelv
41  ntot2 = lx2*ly2*lz2*nelv
42 c
43 c Fill interiors
44  iz1 = 0
45  if (if3d) iz1=1
46  call rzero(v1,ntot1)
47  do e=1,nelv
48  do iz=1,lz2
49  do iy=1,ly2
50  do ix=1,lx2
51  v1(ix+1,iy+1,iz+iz1,e) = v(ix,iy,iz,e)
52  enddo
53  enddo
54  enddo
55  enddo
56  call dface_ext (v1)
57  call dssum (v1,lx1,ly1,lz1)
58  call dface_add1si (v1,-1.)
59 c
60 c Now solve each subdomain problem:
61 c
62  etime1=dnekclock()
63 
64  eoff = 0
65  if (ifield.gt.1) eoff = nelv
66 
67  do e = 1,nelv
68  eb = e + eoff
69  call fastdm1(v1(1,1,1,e),df(1,eb)
70  $ ,sr(1,eb),ss(1,eb),st(1,eb),w1,w2)
71  enddo
72  tsolv=tsolv+dnekclock()-etime1
73 c
74 c Exchange/add elemental solutions
75 c
76  call s_face_to_int (v1,-1.)
77  call dssum (v1,lx1,ly1,lz1)
78  call s_face_to_int (v1, 1.)
79  if(param(42).eq.0) call do_weight_op(v1)
80 c
81 c Map back to pressure grid (extract interior values)
82 c
83  do e=1,nelv
84  do iz=1,lz2
85  do iy=1,ly2
86  do ix=1,lx2
87  u(ix,iy,iz,e) = v1(ix+1,iy+1,iz+iz1,e)
88  enddo
89  enddo
90  enddo
91  enddo
92 c
93  return
94  end
95 c-----------------------------------------------------------------------
96  subroutine fastdm1(r,df,sr,ss,st,w1,w2)
97 c
98 c Fast diagonalization solver for FEM on mesh 1
99 c
100  include 'SIZE'
101  parameter(lxx=lx1*lx1,lxyz=lx1*ly1*lz1)
102 
103  real r(1),df(1),sr(lxx,2),ss(lxx,2),st(lxx,2),w1(1),w2(1)
104 c
105 c
106 c T
107 c S r
108  call tensr3 (w1,lx1,r ,lx1,sr(1,2),ss(1,1),st(1,1),w2)
109 c
110 c
111 c -1 T
112 c D S r
113 c
114  call col2 (w1,df,lxyz)
115 c
116 c
117 c -1 T
118 c S D S r
119 c
120  call tensr3 (r ,lx1,w1,lx1,sr(1,1),ss(1,2),st(1,2),w2)
121 c
122  return
123  end
124 c-----------------------------------------------------------------------
125  subroutine tensr3(v,nv,u,nu,A,Bt,Ct,w)
126 C
127 C - Tensor product application of v = (C x B x A) u
128 C NOTE -- the transpose of B & C must be input, rather than B & C.
129 C
130 C - scratch arrays: w(nu*nu*nv)
131 C
132 C
133  include 'SIZE'
134  include 'INPUT'
135  real v(nv,nv,nv),u(nu,nu,nu)
136  real A(1),Bt(1),Ct(1)
137  real w(1)
138 
139  if (nu.gt.nv) then
140  write(6,*) nid,nu,nv,' ERROR in tensr3. Contact P.Fischer.'
141  write(6,*) nid,nu,nv,' Memory problem.'
142  call exitt
143  endif
144 
145  if (if3d) then
146  nuv = nu*nv
147  nvv = nv*nv
148  call mxm(a,nv,u,nu,v,nu*nu)
149  k=1
150  l=1
151  do iz=1,nu
152  call mxm(v(k,1,1),nv,bt,nu,w(l),nv)
153  k=k+nuv
154  l=l+nvv
155  enddo
156  call mxm(w,nvv,ct,nu,v,nv)
157  else
158  call mxm(a,nv,u,nu,w,nu)
159  call mxm(w,nv,bt,nu,v,nv)
160  endif
161  return
162  end
163 c-----------------------------------------------------------------------
164  subroutine s_face_to_int(x,c)
165 c
166 c Scale face and add to interior of element
167 c
168  include 'SIZE'
169  include 'INPUT'
170  real x(lx1,ly1,lz1,1)
171 c
172  do ie=1,nelv
173 c
174  if (if3d) then
175 c
176  do iz=2,lz1-1
177  do ix=2,lx1-1
178  x(ix,2 ,iz,ie) = c*x(ix,1 ,iz,ie) + x(ix,2 ,iz,ie)
179  x(ix,ly1-1,iz,ie) = c*x(ix,ly1,iz,ie) + x(ix,ly1-1,iz,ie)
180  enddo
181  enddo
182 c
183  do iz=2,lz1-1
184  do iy=2,ly1-1
185  x(2 ,iy,iz,ie) = c*x(1 ,iy,iz,ie) + x(2 ,iy,iz,ie)
186  x(lx1-1,iy,iz,ie) = c*x(lx1,iy,iz,ie) + x(lx1-1,iy,iz,ie)
187  enddo
188  enddo
189 c
190  do iy=2,ly1-1
191  do ix=2,lx1-1
192  x(ix,iy,2 ,ie) = c*x(ix,iy,1 ,ie) + x(ix,iy,2 ,ie)
193  x(ix,iy,lz1-1,ie) = c*x(ix,iy,lz1,ie) + x(ix,iy,lz1-1,ie)
194  enddo
195  enddo
196 c
197  else
198 c 2D
199  do ix=2,lx1-1
200  x(ix,2 ,1,ie) = c*x(ix,1 ,1,ie) + x(ix,2 ,1,ie)
201  x(ix,ly1-1,1,ie) = c*x(ix,ly1,1,ie) + x(ix,ly1-1,1,ie)
202  enddo
203  do iy=2,ly1-1
204  x(2 ,iy,1,ie) = c*x(1 ,iy,1,ie) + x(2 ,iy,1,ie)
205  x(lx1-1,iy,1,ie) = c*x(lx1,iy,1,ie) + x(lx1-1,iy,1,ie)
206  enddo
207  endif
208  enddo
209  return
210  end
211 c-----------------------------------------------------------------------
212  subroutine dface_ext(x)
213 c Extend interior to face of element
214 c
215  include 'SIZE'
216  include 'INPUT'
217  real x(lx1,ly1,lz1,1)
218 c
219  do ie=1,nelv
220 c
221  if (if3d) then
222 c
223  do iz=2,lz1-1
224  do ix=2,lx1-1
225  x(ix,1 ,iz,ie) = x(ix,2 ,iz,ie)
226  x(ix,ly1,iz,ie) = x(ix,ly1-1,iz,ie)
227  enddo
228  enddo
229 c
230  do iz=2,lz1-1
231  do iy=2,ly1-1
232  x(1 ,iy,iz,ie) = x(2 ,iy,iz,ie)
233  x(lx1,iy,iz,ie) = x(lx1-1,iy,iz,ie)
234  enddo
235  enddo
236 c
237  do iy=2,ly1-1
238  do ix=2,lx1-1
239  x(ix,iy,1 ,ie) = x(ix,iy,2 ,ie)
240  x(ix,iy,lz1,ie) = x(ix,iy,lz1-1,ie)
241  enddo
242  enddo
243 c
244  else
245 c
246  do ix=2,lx1-1
247  x(ix,1 ,1,ie) = x(ix,2 ,1,ie)
248  x(ix,ly1,1,ie) = x(ix,ly1-1,1,ie)
249  enddo
250  do iy=2,ly1-1
251  x(1 ,iy,1,ie) = x(2 ,iy,1,ie)
252  x(lx1,iy,1,ie) = x(lx1-1,iy,1,ie)
253  enddo
254 c
255  endif
256  enddo
257  return
258  end
259 c-----------------------------------------------------------------------
260  subroutine dface_add1si(x,c)
261 c Scale interior and add to face of element
262 c
263  include 'SIZE'
264  include 'INPUT'
265  real x(lx1,ly1,lz1,1)
266 c
267  do ie=1,nelv
268 c
269  if (if3d) then
270 c
271  do iz=2,lz1-1
272  do ix=2,lx1-1
273  x(ix,1 ,iz,ie) = x(ix,1 ,iz,ie) + c*x(ix,2 ,iz,ie)
274  x(ix,ly1,iz,ie) = x(ix,ly1,iz,ie) + c*x(ix,ly1-1,iz,ie)
275  enddo
276  enddo
277 c
278  do iz=2,lz1-1
279  do iy=2,ly1-1
280  x(1 ,iy,iz,ie) = x(1 ,iy,iz,ie) + c*x(2 ,iy,iz,ie)
281  x(lx1,iy,iz,ie) = x(lx1,iy,iz,ie) + c*x(lx1-1,iy,iz,ie)
282  enddo
283  enddo
284 c
285  do iy=2,ly1-1
286  do ix=2,lx1-1
287  x(ix,iy,1 ,ie) = x(ix,iy,1 ,ie) + c*x(ix,iy,2 ,ie)
288  x(ix,iy,lz1,ie) = x(ix,iy,lz1,ie) + c*x(ix,iy,lz1-1,ie)
289  enddo
290  enddo
291 c
292  else
293 c
294 c 2D
295 c
296  do ix=2,lx1-1
297  x(ix,1 ,1,ie) = x(ix,1 ,1,ie) + c*x(ix,2 ,1,ie)
298  x(ix,ly1,1,ie) = x(ix,ly1,1,ie) + c*x(ix,ly1-1,1,ie)
299  enddo
300  do iy=2,ly1-1
301  x(1 ,iy,1,ie) = x(1 ,iy,1,ie) + c*x(2 ,iy,1,ie)
302  x(lx1,iy,1,ie) = x(lx1,iy,1,ie) + c*x(lx1-1,iy,1,ie)
303  enddo
304 c
305  endif
306  enddo
307  return
308  end
309 c-----------------------------------------------------------------------
310  subroutine init_weight_op
311 
312  include 'SIZE'
313  include 'INPUT'
314  include 'TSTEP'
315  parameter(levb=lelv+lbelv)
316  common /swaplengths/ l(lx1,ly1,lz1,lelv)
317  common /weightop/ w(lx2,lz2,2,3,levb)
318  real l
319  real w
320  integer e,e0,eb
321 
322  e0 = 0
323  if (ifield.gt.1) e0 = nelv
324 
325  n=lx2+1
326  if (if3d) then
327  do e=1,nelv
328  call rzero(l(1,1,1,e),lx1*ly1*lz1)
329  do k=2,n
330  do j=2,n
331  l(2,j,k,e)=1
332  l(n,j,k,e)=1
333  enddo
334  enddo
335  do k=2,n
336  do i=2,n
337  l(i,2,k,e)=1
338  l(i,n,k,e)=1
339  enddo
340  enddo
341  do j=2,n
342  do i=2,n
343  l(i,j,2,e)=1
344  l(i,j,n,e)=1
345  enddo
346  enddo
347  enddo
348  else
349  do e=1,nelv
350  call rzero(l(1,1,1,e),lx1*ly1*lz1)
351  do j=2,n
352  l(2,j,1,e)=1
353  l(n,j,1,e)=1
354  enddo
355  do i=2,n
356  l(i,2,1,e)=1
357  l(i,n,1,e)=1
358  enddo
359  enddo
360  endif
361 
362  call dface_ext(l)
363  call dssum(l,lx1,ly1,lz1)
364  call dface_add1si(l,-1.)
365  call s_face_to_int(l,1.)
366 c l now holds the count matrix C on the outer pressure nodes
367  if (if3d) then
368  do e=1,nelv
369  eb = e0+e
370  do k=1,lz2
371  do j=1,ly2
372  w(j,k,1,1,eb)=1.0/l(2,j+1,k+1,e)
373  w(j,k,2,1,eb)=1.0/l(n,j+1,k+1,e)
374  enddo
375  enddo
376  do k=1,lz2
377  do i=1,lx2
378  w(i,k,1,2,eb)=1.0/l(i+1,2,k+1,e)
379  w(i,k,2,2,eb)=1.0/l(i+1,n,k+1,e)
380  enddo
381  enddo
382  do j=1,ly2
383  do i=1,lx2
384  w(i,j,1,3,eb)=1.0/l(i+1,j+1,2,e)
385  w(i,j,2,3,eb)=1.0/l(i+1,j+1,n,e)
386  enddo
387  enddo
388  enddo
389  else
390  do e=1,nelv
391  eb = e0+e
392  do j=1,ly2
393  w(j,1,1,1,eb)=1.0/l(2,j+1,1,e)
394  w(j,1,2,1,eb)=1.0/l(n,j+1,1,e)
395  enddo
396  do i=1,lx2
397  w(i,1,1,2,eb)=1.0/l(i+1,2,1,e)
398  w(i,1,2,2,eb)=1.0/l(i+1,n,1,e)
399  enddo
400  enddo
401  endif
402  end
403 c-----------------------------------------------------------------------
404  subroutine do_weight_op(x)
405  include 'SIZE'
406  include 'INPUT'
407  include 'TSTEP'
408  parameter(levb=lelv+lbelv)
409  common /weightop/ w(lx2,lz2,2,3,levb)
410  real w
411 
412  real x(0:lx1-1,0:ly1-1,0:lz1-1,1)
413  integer e,e0,eb
414 
415  e0 = 0
416  if (ifield.gt.1) e0 = nelv
417 
418  if (if3d) then
419  do e=1,nelv
420  eb = e0 + e
421  do k=1,lz2
422  do j=1,ly2
423  x( 1,j,k,e)=w(j,k,1,1,eb)*x( 1,j,k,e)
424  x(lx2,j,k,e)=w(j,k,2,1,eb)*x(lx2,j,k,e)
425  enddo
426  enddo
427  do k=1,lz2
428  do i=2,lx2-1
429  x(i, 1,k,e)=w(i,k,1,2,eb)*x(i, 1,k,e)
430  x(i,ly2,k,e)=w(i,k,2,2,eb)*x(i,ly2,k,e)
431  enddo
432  enddo
433  do j=2,ly2-1
434  do i=2,lx2-1
435  x(i,j, 1,e)=w(i,j,1,3,eb)*x(i,j, 1,e)
436  x(i,j,lz2,e)=w(i,j,2,3,eb)*x(i,j,lz2,e)
437  enddo
438  enddo
439  enddo
440  else
441  do e=1,nelv
442  eb = e0 + e
443  do j=1,ly2
444  x( 1,j,0,e)=w(j,1,1,1,eb)*x( 1,j,0,e)
445  x(lx2,j,0,e)=w(j,1,2,1,eb)*x(lx2,j,0,e)
446  enddo
447  do i=2,lx2-1
448  x(i, 1,0,e)=w(i,1,1,2,eb)*x(i, 1,0,e)
449  x(i,ly2,0,e)=w(i,1,2,2,eb)*x(i,ly2,0,e)
450  enddo
451  enddo
452  endif
453  end
exitt
void exitt()
Definition: comm_mpi.f:604
dnekclock
real *8 function dnekclock()
Definition: comm_mpi.f:393
dssum
subroutine dssum(u, nx, ny, nz)
Definition: dssum.f:34
init_weight_op
subroutine init_weight_op
Definition: fasts.f:311
tensr3
subroutine tensr3(v, nv, u, nu, A, Bt, Ct, w)
Definition: fasts.f:126
s_face_to_int
subroutine s_face_to_int(x, c)
Definition: fasts.f:165
dface_ext
subroutine dface_ext(x)
Definition: fasts.f:213
dface_add1si
subroutine dface_add1si(x, c)
Definition: fasts.f:261
do_weight_op
subroutine do_weight_op(x)
Definition: fasts.f:405
local_solves_fdm
subroutine local_solves_fdm(u, v)
Definition: fasts.f:3
fastdm1
subroutine fastdm1(r, df, sr, ss, st, w1, w2)
Definition: fasts.f:97
col2
subroutine col2(a, b, n)
Definition: math.f:564
rzero
subroutine rzero(a, n)
Definition: math.f:208
mxm
subroutine mxm(a, n1, b, n2, c, n3)
Definition: mxm_wrapper.f:2