KTH framework for Nek5000 toolboxes; testing version  0.0.1
navier6.f
Go to the documentation of this file.
1 c===============================================================================
2 c pff@cfm.brown.edu 3/19/96
3 c
4 c
5 c This is a suite of routines for solving overlapping subdomain
6 c problems with finite elements on distributed memory machines.
7 c
8 c The overall hierarchy of data structures is as follows:
9 c
10 c System - index set denoted by _glob
11 c
12 c Processor - index set denoted by _pglob
13 c
14 c .Domain - index set denoted by _dloc (1,2,3,...,n_dloc)
15 c
16 c .Sp.Elem - index set denoted by _sloc (1,2,3,...,n_sloc)
17 c
18 c
19 c A critical component of the parallel DD solver is the gather-scatter
20 c operation. As there may be more than one domain on a processor, and
21 c communication must be minimized, it is critical that communication be
22 c processor oriented, not domain oriented. Hence domains will access data
23 c via the dloc_to_pglob interface, and the pglob indexed variables will
24 c be accessed via a gather-scatter which interacts via the pglob_glob
25 c interface. Noticed that, in a uni-processor application, the pglob_glob
26 c map will be simply the identity.
27 c
28 c===============================================================================
29 c
30  subroutine set_overlap
31 c
32 c Set up arrays for overlapping Schwartz algorithm *for pressure solver*
33 c
34  include 'SIZE'
35  include 'DOMAIN'
36  include 'ESOLV'
37  include 'INPUT'
38  include 'TSTEP'
39 c
40  real*8 dnekclock,t0
41 c
42  parameter( n_tri = 7*ltotd )
43  common /scrns/ tri(n_tri)
44  integer tri,elem
45 c
46  common /screv/ x(2*ltotd)
47  common /scrvh/ y(2*ltotd)
48  common /scrch/ z(2*ltotd)
49 c
50  common /ctmp0/ nv_to_t(2*ltotd)
51 c
52  parameter(lia = ltotd - 2 - 2*lelt)
53  common /scrcg/ ntri(lelt+1),nmask(lelt+1)
54  $ , ia(lia)
55 c
56  common /scruz/ color(4*ltotd)
57  common /scrmg/ ddmask(4*ltotd)
58  common /ctmp1/ mask(4*ltotd)
59 
60  parameter(lxx=lx1*lx1, levb=lelv+lbelv)
61  common /fastd/ df(lx1*ly1*lz1,levb)
62  $ , sr(lxx*2,levb),ss(lxx*2,levb),st(lxx*2,levb)
63 
64 
65  integer e
66 
67  if (lx1.eq.2) param(43)=1.
68  if (lx1.eq.2.and.nid.eq.0) write(6,*) 'No mgrid for lx1=2!'
69 
70  if (ifaxis) ifmgrid = .false.
71  if (param(43).ne.0) ifmgrid = .false.
72 
73  npass = 1
74  if (ifmhd) npass = 2
75  do ipass=1,npass
76  ifield = 1
77 
78  if (ifsplit.and.ifmgrid) then
79 
80  if (ipass.gt.1) ifield = ifldmhd
81 
82  call swap_lengths
83  call gen_fast_spacing(x,y,z)
84 
85  call hsmg_setup
86  call h1mg_setup
87 
88  elseif (.not.ifsplit) then ! Pn-Pn-2
89 
90  if (ipass.gt.1) ifield = ifldmhd
91 
92  if (param(44).eq.1) then ! Set up local overlapping solves
93  call set_fem_data_l2(nel_proc,ndom,n_o,x,y,z,tri)
94  else
95  call swap_lengths
96  endif
97 
98  e = 1
99  if (ifield.gt.1) e = nelv+1
100 
101  call gen_fast_spacing(x,y,z)
102  call gen_fast(df(1,e),sr(1,e),ss(1,e),st(1,e),x,y,z)
103 
104  call init_weight_op
105  if (param(43).eq.0) call hsmg_setup
106  endif
107 
108  call set_up_h1_crs
109 
110  enddo
111 
112  return
113  end
114 c-----------------------------------------------------------------------
115  subroutine overflow_ck(n_req,n_avail,signal)
116 c
117 c Check for buffer overflow
118 c
119  character*11 signal
120 c
121  nid = mynode()
122 c
123  if (n_req.gt.n_avail) then
124  write(6,9) nid,n_req,n_avail,nid,signal
125  9 format(i7,' ERROR: requested array space (',i12
126  $ ,') exceeds allocated amount (',i12,').'
127  $ ,/,i12,' ABORTING.',3x,a11
128  $ ,/,i12,' ABORTING.',3x,'from overflow_ck call.')
129  call exitt
130  endif
131  return
132  end
133 c-----------------------------------------------------------------------
134  subroutine iunswap(b,ind,n,temp)
135  integer b(1),ind(1),temp(1)
136 c
137 c sort associated elements by putting item(jj)
138 c into item(i), where jj=ind(i).
139 c
140  do 20 i=1,n
141  jj=ind(i)
142  temp(jj)=b(i)
143  20 continue
144  do 30 i=1,n
145  30 b(i)=temp(i)
146  return
147  end
148 c-----------------------------------------------------------------------
149  subroutine set_fem_data_l2(nep,nd,no,x,y,z,p)
150  include 'SIZE'
151  include 'DOMAIN'
152  include 'NONCON'
153  include 'TOTAL'
154 c
155  real x(lx1,ly1,lz1,1),y(lx1,ly1,lz1,1),z(lx1,ly1,lz1,1)
156  real p(lx1,ly1,lz1,1)
157  integer ce,cf
158 c
159  common /ctmp0/ w1(lx1,ly1),w2(lx1,ly1)
160 c
161 c First, copy local geometry to temporary, expanded, arrays
162 c
163  nxy1 = lx1*ly1
164  nxy2 = lx2*ly2
165  nxyz1 = lx1*ly1*lz1
166  nxyz2 = lx2*ly2*lz2
167  ntot1 = nelv*nxyz1
168  ntot2 = nelv*nxyz2
169 c
170  call rzero(x,ntot1)
171  call rzero(y,ntot1)
172  call rzero(z,ntot1)
173 c
174 c Take care of surface geometry
175 c
176  call map_face12(x,xm1,w1,w2)
177  call map_face12(y,ym1,w1,w2)
178  call map_face12(z,zm1,w1,w2)
179 c
180 c Take care of interior geometry
181 c
182  iz1 = 0
183  if (if3d) iz1=1
184  do ie=1,nelv
185  do iz=1,lz2
186  do iy=1,ly2
187  do ix=1,lx2
188  x(ix+1,iy+1,iz+iz1,ie) = xm2(ix,iy,iz,ie)
189  y(ix+1,iy+1,iz+iz1,ie) = ym2(ix,iy,iz,ie)
190  z(ix+1,iy+1,iz+iz1,ie) = zm2(ix,iy,iz,ie)
191  enddo
192  enddo
193  enddo
194  enddo
195 c
196 c
197 c Compute absolute value of distance between pressure and vel. planes
198 c
199  call dface_add1sa(x)
200  call dface_add1sa(y)
201  call dface_add1sa(z)
202 c
203 c Sum this distance
204 c
205  ifielt = ifield
206  ifield = 1
207 c call dssum(x,lx1,ly1,lz1)
208 c call dssum(y,lx1,ly1,lz1)
209 c call dssum(z,lx1,ly1,lz1)
210 c
211 c Scale children values by 0.5. This is a hack, based on assumption
212 c that number of children sharing a parent edge is 2
213 c
214  scale = 0.5
215  if (if3d) scale = 0.25
216  do im = 1 , mort_m
217  cf = noncon_f(im)
218  ce = noncon_e(im)
219  call faces(x,scale,ce,cf,lx1,ly1,lz1)
220  call faces(y,scale,ce,cf,lx1,ly1,lz1)
221  if (if3d) call faces(z,scale,ce,cf,lx1,ly1,lz1)
222  enddo
223  call fgslib_gs_op_many(gsh_fld(ifield), x,y,z,x,x,x,ldim, 1,1,0)
224 c
225  ifield = ifielt
226 c
227 c Dirichlet (pmask=0) faces all set...
228 c
229  nel_proc = nelv
230  ndom = nelv
231 c ndom = 1
232 c
233  n_o = 1
234 c
235 c These are the return values, since the calling routine doesn't
236 c know DOMAIN
237 c
238  nep = nel_proc
239  nd = ndom
240  no = n_o
241  ifield = ifielt
242  return
243  end
244 c-----------------------------------------------------------------------
245  subroutine map_face12(x2,x1,w1,w2)
246  include 'SIZE'
247  include 'INPUT'
248  include 'IXYZ'
249 c
250 c Interpolate iface of x1 (GLL pts) onto x2 (GL pts).
251 c Work arrays should be of size lx1*lx1 each.
252 c
253  real x2(lx1*ly1*lz1,1)
254  real x1(lx1*ly1*lz1,1)
255  real w1(1),w2(1)
256 c
257 c
258  do ie=1,nelv
259  if (if3d) then
260  call map_one_face12(x2(1,ie),x1(1,ie),1,ixm12,ixtm12,w1,w2)
261  call map_one_face12(x2(1,ie),x1(1,ie),2,ixm12,ixtm12,w1,w2)
262  call map_one_face12(x2(1,ie),x1(1,ie),3,ixm12,ixtm12,w1,w2)
263  call map_one_face12(x2(1,ie),x1(1,ie),4,ixm12,ixtm12,w1,w2)
264  call map_one_face12(x2(1,ie),x1(1,ie),5,ixm12,ixtm12,w1,w2)
265  call map_one_face12(x2(1,ie),x1(1,ie),6,ixm12,ixtm12,w1,w2)
266  else
267  call map_one_face12(x2(1,ie),x1(1,ie),1,ixm12,ixtm12,w1,w2)
268  call map_one_face12(x2(1,ie),x1(1,ie),2,ixm12,ixtm12,w1,w2)
269  call map_one_face12(x2(1,ie),x1(1,ie),3,ixm12,ixtm12,w1,w2)
270  call map_one_face12(x2(1,ie),x1(1,ie),4,ixm12,ixtm12,w1,w2)
271  endif
272  enddo
273 c
274  return
275  end
276 c-----------------------------------------------------------------------
277  subroutine map_one_face12(x2,x1,iface,i12,i12t,w1,w2)
278  include 'SIZE'
279  include 'INPUT'
280 c
281 c Interpolate iface of x1 (GLL pts) onto interior of face of x2 (GL pts).
282 c Work arrays should be of size lx1*lx1 each.
283 c
284  real x2(lx1,ly1,lz1)
285  real x1(lx1,ly1,lz1)
286  real w1(1),w2(1)
287  real i12(1),i12t(1)
288 c
289 c
290 c Extract surface data from x1
291 c
292  CALL facind (kx1,kx2,ky1,ky2,kz1,kz2,lx1,ly1,lz1,iface)
293  i=0
294  do iz=kz1,kz2
295  do iy=ky1,ky2
296  do ix=kx1,kx2
297  i = i+1
298  w1(i) = x1(ix,iy,iz)
299  enddo
300  enddo
301  enddo
302 c
303 c Interpolate from mesh 1 to 2
304 c
305  if (if3d) then
306  call mxm(i12 ,lx2,w1,lx1,w2,lx1)
307  call mxm(w2,lx2,i12t,lx1,w1,lx2)
308  else
309  call mxm(i12 ,lx2,w1,lx1,w2, 1)
310  call copy(w1,w2,lx2)
311  endif
312 c
313 c Write to interior points on face of x2
314 c
315  kx1=min(kx1+1,lx1,kx2)
316  kx2=max(kx2-1, 1,kx1)
317  ky1=min(ky1+1,ly1,ky2)
318  ky2=max(ky2-1, 1,ky1)
319  kz1=min(kz1+1,lz1,kz2)
320  kz2=max(kz2-1, 1,kz1)
321 c
322  i=0
323  do iz=kz1,kz2
324  do iy=ky1,ky2
325  do ix=kx1,kx2
326  i = i+1
327  x2(ix,iy,iz) = w1(i)
328  enddo
329  enddo
330  enddo
331 c
332  return
333  end
334 c-----------------------------------------------------------------------
335  subroutine dface_add1sa(x)
336 c Compute |face-interior|
337 c
338  include 'SIZE'
339  include 'INPUT'
340  real x(lx1,ly1,lz1,1)
341 c
342  do ie=1,nelv
343 c
344  if (if3d) then
345 c
346  do iz=2,lz1-1
347  do ix=2,lx1-1
348  x(ix,1 ,iz,ie)=abs(x(ix,1 ,iz,ie) - x(ix,2 ,iz,ie))
349  x(ix,ly1,iz,ie)=abs(x(ix,ly1,iz,ie) - x(ix,ly1-1,iz,ie))
350  enddo
351  enddo
352 c
353  do iz=2,lz1-1
354  do iy=2,ly1-1
355  x(1 ,iy,iz,ie)=abs(x(1 ,iy,iz,ie) - x(2 ,iy,iz,ie))
356  x(lx1,iy,iz,ie)=abs(x(lx1,iy,iz,ie) - x(lx1-1,iy,iz,ie))
357  enddo
358  enddo
359 c
360  do iy=2,ly1-1
361  do ix=2,lx1-1
362  x(ix,iy,1 ,ie)=abs(x(ix,iy,1 ,ie) - x(ix,iy,2 ,ie))
363  x(ix,iy,lz1,ie)=abs(x(ix,iy,lz1,ie) - x(ix,iy,lz1-1,ie))
364  enddo
365  enddo
366 c
367  else
368 c
369 c 2D
370 c
371  do ix=2,lx1-1
372  x(ix,1 ,1,ie)=abs(x(ix,1 ,1,ie) - x(ix,2 ,1,ie))
373  x(ix,ly1,1,ie)=abs(x(ix,ly1,1,ie) - x(ix,ly1-1,1,ie))
374  enddo
375  do iy=2,ly1-1
376  x(1 ,iy,1,ie)=abs(x(1 ,iy,1,ie) - x(2 ,iy,1,ie))
377  x(lx1,iy,1,ie)=abs(x(lx1,iy,1,ie) - x(lx1-1,iy,1,ie))
378  enddo
379 c
380  endif
381  enddo
382  return
383  end
384 c-----------------------------------------------------------------------
385  subroutine faces(a,s,ie,iface,nx,ny,nz)
386 C
387 C Scale face(IFACE,IE) of array A by s.
388 C IFACE is the input in the pre-processor ordering scheme.
389 C
390  include 'SIZE'
391  dimension a(nx,ny,nz,lelt)
392  CALL facind (kx1,kx2,ky1,ky2,kz1,kz2,nx,ny,nz,iface)
393  DO 100 iz=kz1,kz2
394  DO 100 iy=ky1,ky2
395  DO 100 ix=kx1,kx2
396  a(ix,iy,iz,ie)=s*a(ix,iy,iz,ie)
397  100 CONTINUE
398  RETURN
399  END
400 c-----------------------------------------------------------------------
exitt
void exitt()
Definition: comm_mpi.f:604
mynode
integer function mynode()
Definition: comm_mpi.f:382
dnekclock
real *8 function dnekclock()
Definition: comm_mpi.f:393
facind
subroutine facind(kx1, kx2, ky1, ky2, kz1, kz2, nx, ny, nz, iface)
Definition: connect1.f:1028
scale
subroutine scale(xyzl, nl)
Definition: connect2.f:602
swap_lengths
subroutine swap_lengths
Definition: fast3d.f:1541
gen_fast_spacing
subroutine gen_fast_spacing(x, y, z)
Definition: fast3d.f:143
gen_fast
subroutine gen_fast(df, sr, ss, st, x, y, z)
Definition: fast3d.f:3
init_weight_op
subroutine init_weight_op
Definition: fasts.f:311
hsmg_setup
subroutine hsmg_setup()
Definition: hsmg.f:23
h1mg_setup
subroutine h1mg_setup()
Definition: hsmg.f:2235
copy
subroutine copy(a, b, n)
Definition: math.f:260
rzero
subroutine rzero(a, n)
Definition: math.f:208
mxm
subroutine mxm(a, n1, b, n2, c, n3)
Definition: mxm_wrapper.f:2
iunswap
subroutine iunswap(b, ind, n, temp)
Definition: navier6.f:135
map_face12
subroutine map_face12(x2, x1, w1, w2)
Definition: navier6.f:246
set_overlap
subroutine set_overlap
Definition: navier6.f:31
set_fem_data_l2
subroutine set_fem_data_l2(nep, nd, no, x, y, z, p)
Definition: navier6.f:150
overflow_ck
subroutine overflow_ck(n_req, n_avail, signal)
Definition: navier6.f:116
dface_add1sa
subroutine dface_add1sa(x)
Definition: navier6.f:336
faces
subroutine faces(a, s, ie, iface, nx, ny, nz)
Definition: navier6.f:386
map_one_face12
subroutine map_one_face12(x2, x1, iface, i12, i12t, w1, w2)
Definition: navier6.f:278
set_up_h1_crs
subroutine set_up_h1_crs
Definition: navier8.f:84