KTH framework for Nek5000 toolboxes; testing version  0.0.1
dssum.f
Go to the documentation of this file.
1  subroutine setupds(gs_handle,nx,ny,nz,nel,melg,vertex,glo_num)
2  include 'SIZE'
3  include 'INPUT'
4  include 'PARALLEL'
5  include 'NONCON'
6 
7  integer gs_handle
8  integer vertex(1)
9  integer*8 glo_num(1),ngv
10 
11  common /nekmpi/ mid,mp,nekcomm,nekgroup,nekreal
12 
13  t0 = dnekclock()
14 
15 c Global-to-local mapping for gs
16  call set_vert(glo_num,ngv,nx,nel,vertex,.false.)
17 
18 c Initialize gather-scatter code
19  ntot = nx*ny*nz*nel
20  call fgslib_gs_setup(gs_handle,glo_num,ntot,nekcomm,mp)
21 
22 c call gs_chkr(glo_num)
23 
24  t1 = dnekclock() - t0
25  if (nio.eq.0) then
26  write(6,1) t1,gs_handle,nx,ngv,melg
27  1 format(' setupds time',1pe11.4,' seconds ',2i3,2i12)
28  endif
29 c
30  return
31  end
32 c-----------------------------------------------------------------------
33  subroutine dssum(u,nx,ny,nz)
34  include 'SIZE'
35  include 'CTIMER'
36  include 'INPUT'
37  include 'NONCON'
38  include 'PARALLEL'
39  include 'TSTEP'
40  real u(1)
41 
42  parameter(lface=lx1*ly1)
43  common /nonctmp/ uin(lface,2*ldim),uout(lface)
44 
45  ifldt = ifield
46 c if (ifldt.eq.0) ifldt = 1
47  if (ifldt.eq.ifldmhd) ifldt = 1
48 c write(6,*) ifldt,ifield,gsh_fld(ifldt),imesh,' ifldt'
49 
50  if (ifsync) call nekgsync()
51 
52 #ifdef TIMER
53  if (icalld.eq.0) then
54  tdsmx=0.
55  tdsmn=0.
56  endif
57  icalld=icalld+1
58  etime1=dnekclock()
59 #endif
60 
61 c
62 c T ~ ~T T
63 c Implement QQ := J Q Q J
64 c
65 c
66 c T
67 c This is the J part, translating child data
68 c
69 c call apply_Jt(u,nx,ny,nz,nel)
70 c
71 c
72 c
73 c ~ ~T
74 c This is the Q Q part
75 c
76  if (gsh_fld(ifldt).ge.0) then
77  if (nio.eq.0.and.loglevel.gt.5)
78  $ write(6,*) 'dssum', ifldt
79  call fgslib_gs_op(gsh_fld(ifldt),u,1,1,0) ! 1 ==> +
80  endif
81 c
82 c
83 c
84 c This is the J part, interpolating parent solution onto child
85 c
86 c call apply_J(u,nx,ny,nz,nel)
87 c
88 c
89 #ifdef TIMER
90  timee=(dnekclock()-etime1)
91  tdsum=tdsum+timee
92  ndsum=icalld
93  tdsmx=max(timee,tdsmx)
94  tdsmn=min(timee,tdsmn)
95 #endif
96 c
97  return
98  end
99 c-----------------------------------------------------------------------
100  subroutine dsop(u,op,nx,ny,nz)
101  include 'SIZE'
102  include 'PARALLEL'
103  include 'INPUT'
104  include 'TSTEP'
105  include 'CTIMER'
106 
107  real u(1)
108  character*3 op
109  character*10 s1,s2
110 c
111 c o gs recognized operations:
112 c
113 c o "+" ==> addition.
114 c o "*" ==> multiplication.
115 c o "M" ==> maximum.
116 c o "m" ==> minimum.
117 c o "A" ==> (fabs(x)>fabs(y)) ? (x) : (y), ident=0.0.
118 c o "a" ==> (fabs(x)<fabs(y)) ? (x) : (y), ident=MAX_DBL
119 c o "e" ==> ((x)==0.0) ? (y) : (x), ident=0.0.
120 c
121 c o note: a binary function pointer flavor exists.
122 c
123 c
124 c o gs level:
125 c
126 c o level=0 ==> pure tree
127 c o level>=num_nodes-1 ==> pure pairwise
128 c o level = 1,...num_nodes-2 ==> mix tree/pairwise.
129 c
130 c
131  ifldt = ifield
132 c if (ifldt.eq.0) ifldt = 1
133  if (ifldt.eq.ifldmhd) ifldt = 1
134 
135 c if (nio.eq.0)
136 c $ write(6,*) istep,' dsop: ',op,ifield,ifldt,gsh_fld(ifldt)
137 
138  if(ifsync) call nekgsync()
139 
140  if (op.eq.'+ ') call fgslib_gs_op(gsh_fld(ifldt),u,1,1,0)
141  if (op.eq.'sum') call fgslib_gs_op(gsh_fld(ifldt),u,1,1,0)
142  if (op.eq.'SUM') call fgslib_gs_op(gsh_fld(ifldt),u,1,1,0)
143 
144  if (op.eq.'* ') call fgslib_gs_op(gsh_fld(ifldt),u,1,2,0)
145  if (op.eq.'mul') call fgslib_gs_op(gsh_fld(ifldt),u,1,2,0)
146  if (op.eq.'MUL') call fgslib_gs_op(gsh_fld(ifldt),u,1,2,0)
147 
148  if (op.eq.'m ') call fgslib_gs_op(gsh_fld(ifldt),u,1,3,0)
149  if (op.eq.'min') call fgslib_gs_op(gsh_fld(ifldt),u,1,3,0)
150  if (op.eq.'mna') call fgslib_gs_op(gsh_fld(ifldt),u,1,3,0)
151  if (op.eq.'MIN') call fgslib_gs_op(gsh_fld(ifldt),u,1,3,0)
152  if (op.eq.'MNA') call fgslib_gs_op(gsh_fld(ifldt),u,1,3,0)
153 
154  if (op.eq.'M ') call fgslib_gs_op(gsh_fld(ifldt),u,1,4,0)
155  if (op.eq.'max') call fgslib_gs_op(gsh_fld(ifldt),u,1,4,0)
156  if (op.eq.'mxa') call fgslib_gs_op(gsh_fld(ifldt),u,1,4,0)
157  if (op.eq.'MAX') call fgslib_gs_op(gsh_fld(ifldt),u,1,4,0)
158  if (op.eq.'MXA') call fgslib_gs_op(gsh_fld(ifldt),u,1,4,0)
159 c
160  return
161  end
162 c-----------------------------------------------------------------------
163  subroutine vec_dssum(u,v,w,nx,ny,nz)
164 c
165 c Direct stiffness summation of the face data, for field U.
166 c
167 c Boundary condition data corresponds to component IFIELD of
168 c the CBC array.
169 c
170  include 'SIZE'
171  include 'TOPOL'
172  include 'INPUT'
173  include 'PARALLEL'
174  include 'TSTEP'
175  include 'CTIMER'
176 
177  REAL U(1),V(1),W(1)
178 
179  if(ifsync) call nekgsync()
180 
181 #ifdef TIMER
182  if (icalld.eq.0) tvdss=0.0d0
183  if (icalld.eq.0) tgsum=0.0d0
184  icalld=icalld+1
185  nvdss=icalld
186  etime1=dnekclock()
187 #endif
188 
189 c
190 c============================================================================
191 c execution phase
192 c============================================================================
193 c
194  ifldt = ifield
195 c if (ifldt.eq.0) ifldt = 1
196  if (ifldt.eq.ifldmhd) ifldt = 1
197 
198  call fgslib_gs_op_many(gsh_fld(ifldt),u,v,w,u,u,u,ldim,1,1,0)
199 
200 #ifdef TIMER
201  timee=(dnekclock()-etime1)
202  tvdss=tvdss+timee
203  tdsmx=max(timee,tdsmx)
204  tdsmn=min(timee,tdsmn)
205 #endif
206 
207  return
208  end
209 
210 c-----------------------------------------------------------------------
211  subroutine vec_dsop(u,v,w,nx,ny,nz,op)
212 c
213 c Direct stiffness summation of the face data, for field U.
214 c
215 c Boundary condition data corresponds to component IFIELD of
216 c the CBC array.
217 c
218  include 'SIZE'
219  include 'TOPOL'
220  include 'INPUT'
221  include 'PARALLEL'
222  include 'TSTEP'
223  include 'CTIMER'
224 c
225  real u(1),v(1),w(1)
226  character*3 op
227 
228 c============================================================================
229 c execution phase
230 c============================================================================
231 
232  ifldt = ifield
233 c if (ifldt.eq.0) ifldt = 1
234  if (ifldt.eq.ifldmhd) ifldt = 1
235 
236 c write(6,*) 'opdsop: ',op,ifldt,ifield
237  if(ifsync) call nekgsync()
238 
239  if (op.eq.'+ ' .or. op.eq.'sum' .or. op.eq.'SUM')
240  $ call fgslib_gs_op_many(gsh_fld(ifldt),u,v,w,u,u,u,ldim,1,1,0)
241 
242 
243  if (op.eq.'* ' .or. op.eq.'mul' .or. op.eq.'MUL')
244  $ call fgslib_gs_op_many(gsh_fld(ifldt),u,v,w,u,u,u,ldim,1,2,0)
245 
246 
247  if (op.eq.'m ' .or. op.eq.'min' .or. op.eq.'mna'
248  $ .or. op.eq.'MIN' .or. op.eq.'MNA')
249  $ call fgslib_gs_op_many(gsh_fld(ifldt),u,v,w,u,u,u,ldim,1,3,0)
250 
251 
252  if (op.eq.'M ' .or. op.eq.'max' .or. op.eq.'mxa'
253  $ .or. op.eq.'MAX' .or. op.eq.'MXA')
254  $ call fgslib_gs_op_many(gsh_fld(ifldt),u,v,w,u,u,u,ldim,1,4,0)
255 
256 
257  return
258  end
259 c-----------------------------------------------------------------------
260  subroutine nvec_dssum(u,stride,n,gs_handle)
261 
262 c Direct stiffness summation of the array u for n fields
263 c
264  include 'SIZE'
265  include 'CTIMER'
266 
267  real u(1)
268  integer n,stride,gs_handle
269 
270  if(ifsync) call nekgsync()
271 
272 #ifdef TIMER
273  icalld=icalld+1
274  nvdss=icalld
275  etime1=dnekclock()
276 #endif
277  call fgslib_gs_op_fields(gs_handle,u,stride,n,1,1,0)
278 
279 #ifdef TIMER
280  timee=(dnekclock()-etime1)
281  tvdss=tvdss+timee
282  tdsmx=max(timee,tdsmx)
283  tdsmn=min(timee,tdsmn)
284 #endif
285 
286  return
287  end
288 
289 c----------------------------------------------------------------------
290  subroutine matvec3(uout,Jmat,uin,iftrsp,n1,n2)
291 c
292  include 'SIZE'
293 c
294  real Jmat (n1,n1,2)
295  real uin (1)
296  real uout (1)
297  logical iftrsp
298 c
299  common /matvtmp/ utmp(lx1,ly1)
300 c
301  if (ldim.eq.2) then
302  call mxm (jmat(1,1,1),n1,uin,n1,uout,n2)
303  else
304  if (iftrsp) then
305  call transpose(uout,n2,uin,n1)
306  else
307  call copy (uout,uin,n1*n2)
308  endif
309  call mxm (jmat(1,1,1),n1,uout,n1,utmp,n2)
310  call mxm (utmp,n2,jmat(1,1,2),n1,uout,n1)
311  endif
312 c
313  return
314  end
315 c-----------------------------------------------------------------------
316  subroutine matvec3t(uout,Jmat,uin,iftrsp,n1,n2)
317 c
318  include 'SIZE'
319 c
320  real Jmat (n1,n1,2)
321  real uin (n1,n2)
322  real uout (n1,n2)
323  logical iftrsp
324 c
325  common /matvtmp/ utmp(lx1*ly1)
326 c
327  call transpose(utmp,n2,uin,n1)
328  call mxm (jmat(1,1,2),n1,utmp,n1,uout,n2)
329  call mxm (uout,n2,jmat(1,1,1),n1,utmp,n1)
330  if (iftrsp) then
331  call copy (uout,utmp,n1*n2)
332  else
333  call transpose(uout,n2,utmp,n1)
334  endif
335 c
336  return
337  end
338 c-----------------------------------------------------------------------
339  subroutine matvect (out,d,vec,n1,n2)
340  dimension d(n1,n2),out(1),vec(1)
341 c
342 c handle non-square matrix in mat-vec mult -- TRANSPOSE
343 c N1 is still the number of rows
344 c N2 is still the number of cols
345 c
346 c
347  call mxm(vec,1,d,n1,out,n2)
348 c
349  return
350  end
351 c-----------------------------------------------------------------------
352 c subroutine opq_in_place(a,b,c)
353 c include 'SIZE'
354 c real a(1),b(1),c(1)
355 c
356 c call q_in_place(a)
357 c call q_in_place(b)
358 c if (ldim .eq.3) call q_in_place(c)
359 c
360 c return
361 c end
362 c-----------------------------------------------------------------------
363  subroutine vectof_add(b,a,ie,iface,nx,ny,nz)
364 C
365 C Copy vector A to the face (IFACE) of B
366 C IFACE is the input in the pre-processor ordering scheme.
367 C
368  dimension a(nx,ny)
369  dimension b(nx,ny,nz,1)
370  CALL facind (kx1,kx2,ky1,ky2,kz1,kz2,nx,ny,nz,iface)
371  k = 0
372  DO 100 iz=kz1,kz2
373  DO 100 iy=ky1,ky2
374  DO 100 ix=kx1,kx2
375  k = k + 1
376  b(ix,iy,iz,ie) = b(ix,iy,iz,ie) + a(k,1)
377  100 CONTINUE
378  return
379  END
380 c-----------------------------------------------------------------------
381  subroutine zero_f(b,ie,iface,nx,ny,nz)
382 C
383 C ZERO the face (IFACE) of B
384 C IFACE is the input in the pre-processor ordering scheme.
385 C
386  dimension b(nx,ny,nz,1)
387  CALL facind (kx1,kx2,ky1,ky2,kz1,kz2,nx,ny,nz,iface)
388 c
389  DO 100 iz=kz1,kz2
390  DO 100 iy=ky1,ky2
391  DO 100 ix=kx1,kx2
392  b(ix,iy,iz,ie) = 0.
393  100 CONTINUE
394  return
395  END
396 c-----------------------------------------------------------------------
397  subroutine ftovec_0(a,b,ie,iface,nx,ny,nz)
398 C
399 C Copy the face (IFACE) of B to vector A.
400 C IFACE is the input in the pre-processor ordering scheme.
401 C
402  dimension a(nx,ny)
403  dimension b(nx,ny,nz,1)
404  CALL facind (kx1,kx2,ky1,ky2,kz1,kz2,nx,ny,nz,iface)
405  k = 0
406  DO 100 iz=kz1,kz2
407  DO 100 iy=ky1,ky2
408  DO 100 ix=kx1,kx2
409  k = k + 1
410  a(k,1)=b(ix,iy,iz,ie)
411  b(ix,iy,iz,ie)=0.0
412  100 CONTINUE
413  return
414  END
415 c-----------------------------------------------------------------------
416  subroutine ftovec(a,b,ie,iface,nx,ny,nz)
417 C
418 C Copy the face (IFACE) of B to vector A.
419 C IFACE is the input in the pre-processor ordering scheme.
420 C
421  real A(NX,NY)
422  real B(NX,NY,NZ,1)
423  CALL facind (kx1,kx2,ky1,ky2,kz1,kz2,nx,ny,nz,iface)
424  k = 0
425  DO 100 iz=kz1,kz2
426  DO 100 iy=ky1,ky2
427  DO 100 ix=kx1,kx2
428  k = k + 1
429  a(k,1)=b(ix,iy,iz,ie)
430  100 CONTINUE
431  return
432  END
433 c-----------------------------------------------------------------------
434  subroutine vectof(b,a,ie,iface,nx,ny,nz)
435 C
436 C Copy vector A to the face (IFACE) of B
437 C IFACE is the input in the pre-processor ordering scheme.
438 C
439  real A(NX,NY)
440  real B(NX,NY,NZ,1)
441  CALL facind (kx1,kx2,ky1,ky2,kz1,kz2,nx,ny,nz,iface)
442  k = 0
443  DO 100 iz=kz1,kz2
444  DO 100 iy=ky1,ky2
445  DO 100 ix=kx1,kx2
446  k = k + 1
447  b(ix,iy,iz,ie) = a(k,1)
448  100 CONTINUE
449  return
450  END
451 c-----------------------------------------------------------------------
452  subroutine ftoveci(a,b,ie,iface,nx,ny,nz)
453 C
454 C Copy the face (IFACE) of B to vector A.
455 C IFACE is the input in the pre-processor ordering scheme.
456 C
457  integer A(NX,NY)
458  integer B(NX,NY,NZ,1)
459  CALL facind (kx1,kx2,ky1,ky2,kz1,kz2,nx,ny,nz,iface)
460  k = 0
461  DO 100 iz=kz1,kz2
462  DO 100 iy=ky1,ky2
463  DO 100 ix=kx1,kx2
464  k = k + 1
465  a(k,1)=b(ix,iy,iz,ie)
466  100 CONTINUE
467  return
468  END
469 c-----------------------------------------------------------------------
470  subroutine vectofi(b,a,ie,iface,nx,ny,nz)
471 C
472 C Copy vector A to the face (IFACE) of B
473 C IFACE is the input in the pre-processor ordering scheme.
474 C
475  integer A(NX,NY)
476  integer B(NX,NY,NZ,1)
477  CALL facind (kx1,kx2,ky1,ky2,kz1,kz2,nx,ny,nz,iface)
478  k = 0
479  DO 100 iz=kz1,kz2
480  DO 100 iy=ky1,ky2
481  DO 100 ix=kx1,kx2
482  k = k + 1
483  b(ix,iy,iz,ie) = a(k,1)
484  100 CONTINUE
485  return
486  END
487 c-----------------------------------------------------------------------
488  subroutine apply_jt(u,nx,ny,nz,nel)
489  include 'SIZE'
490  include 'CTIMER'
491  include 'INPUT'
492  include 'NONCON'
493  include 'PARALLEL'
494  include 'TSTEP'
495  real u(1)
496 c
497  parameter(lface=lx1*ly1)
498  common /nonctmp/ uin(lface,2*ldim),uout(lface)
499 c
500 c
501 c T
502 c This is the J part, translating child data
503 c
504  do ie = 1 , nel
505 c Note, we zero out u() on this face after extracting, for
506 c consistency reasons discovered during Jerry's thesis.
507 c Thus, "ftovec_0" rather than ftovec(). (iface -- Ed notation)
508  do iface = 1 , 2*ldim
509  im = mortar(iface,ie)
510  if (im.ne.0) then
511  call ftovec_0(uin(1,iface),u,ie,iface,nx,ny,nz)
512  endif
513  enddo
514  do iface=1,2*ldim
515  im = mortar(iface,ie)
516  if (im.ne.0) then
517  if (if3d) then
518  call matvec3t
519  $ (uout,jmat(1,1,1,im),uin(1,iface),ifjt(im),nx,nx)
520  else
521  call matvect (uout,jmat(1,1,1,im),uin(1,iface),nx,nx)
522  endif
523  call vectof_add(u,uout,ie,iface,nx,ny,nz)
524  endif
525  enddo
526  enddo
527 c
528  return
529  end
530 c-----------------------------------------------------------------------
531  subroutine apply_j(u,nx,ny,nz,nel)
532  include 'SIZE'
533  include 'CTIMER'
534  include 'INPUT'
535  include 'NONCON'
536  include 'PARALLEL'
537  include 'TSTEP'
538  real u(1)
539 c
540  parameter(lface=lx1*ly1)
541  common /nonctmp/ uin(lface,2*ldim),uout(lface)
542 c
543 c This is the J part, interpolating parent solution onto child
544 c
545 c
546  do ie = 1 , nel
547  do iface = 1 , 2*ldim
548  im = mortar(iface,ie)
549  if (im.ne.0) then
550  call ftovec(uin(1,iface),u,ie,iface,nx,ny,nz)
551  endif
552  enddo
553  do iface=1,2*ldim
554  im = mortar(iface,ie)
555  if (im.ne.0) then
556  call matvec3
557  $ (uout,jmat(1,1,1,im),uin(1,iface),ifjt(im),nx,nz)
558  call vectof (u,uout,ie,iface,nx,ny,nz)
559  endif
560  enddo
561  enddo
562 c
563  return
564  end
565 c-----------------------------------------------------------------------
566  subroutine h1_proj(u,nx,ny,nz)
567  include 'SIZE'
568  include 'CTIMER'
569  include 'INPUT'
570  include 'NONCON'
571  include 'PARALLEL'
572  include 'TSTEP'
573  real u(1)
574 c
575  parameter(lface=lx1*ly1)
576  common /nonctmp/ uin(lface,2*ldim),uout(lface)
577 
578  if(ifsync) call nekgsync()
579 
580 #ifdef TIMER
581  if (icalld.eq.0) then
582  tdsmx=0.
583  tdsmn=0.
584  endif
585  icalld=icalld+1
586  etime1=dnekclock()
587 #endif
588 c
589  ifldt = ifield
590 c if (ifldt.eq.0) ifldt = 1
591  nel = nelv
592  if (ifield.ge.2) nel=nelt
593  ntot = nx*ny*nz*nel
594 
595 
596 c
597 c
598 c ~ ~T
599 c Implement := J Q Q Mu
600 c
601 c
602 c T
603 c
604  call col2 (u,umult,ntot)
605 c
606 c ~ ~T
607 c This is the Q Q part
608 c
609  call fgslib_gs_op(gsh_fld(ifldt),u,1,1,0)
610 c
611 c
612 c This is the J part, interpolating parent solution onto child
613 c
614  call apply_j(u,nx,ny,nz,nel)
615 c
616 c
617 #ifdef TIMER
618  timee=(dnekclock()-etime1)
619  tdsum=tdsum+timee
620  ndsum=icalld
621  tdsmx=max(timee,tdsmx)
622  tdsmn=min(timee,tdsmn)
623 #endif
624 c
625  return
626  end
627 c-----------------------------------------------------------------------
628  subroutine dssum_msk(u,mask,nx,ny,nz)
629  include 'SIZE'
630  include 'CTIMER'
631  include 'INPUT'
632  include 'NONCON'
633  include 'PARALLEL'
634  include 'TSTEP'
635  real u(1),mask(1)
636 c
637  parameter(lface=lx1*ly1)
638  common /nonctmp/ uin(lface,2*ldim),uout(lface)
639 
640  if(ifsync) call nekgsync()
641 
642 #ifdef TIMER
643  if (icalld.eq.0) then
644  tdsmx=0.
645  tdsmn=0.
646  endif
647  icalld=icalld+1
648  etime1=dnekclock()
649 #endif
650 c
651  ifldt = ifield
652 c if (ifldt.eq.0) ifldt = 1
653  nel = nelv
654  if (ifield.ge.2) nel=nelt
655  ntot = nx*ny*nz*nel
656 
657 
658 c
659 c T ~ ~T T
660 c Implement Q M Q := J M Q Q J
661 c
662 c
663 c T
664 c This is the J part, translating child data
665 c
666  call apply_jt(u,nx,ny,nz,nel)
667 c
668 c
669 c
670 c ~ ~T
671 c This is the Q Q part
672 c
673  call fgslib_gs_op(gsh_fld(ifldt),u,1,1,0)
674  call col2 (u,mask,ntot)
675 c
676 c
677 c This is the J part, interpolating parent solution onto child
678 c
679  call apply_j(u,nx,ny,nz,nel)
680 c
681 c
682 #ifdef TIMER
683  timee=(dnekclock()-etime1)
684  tdsum=tdsum+timee
685  ndsum=icalld
686  tdsmx=max(timee,tdsmx)
687  tdsmn=min(timee,tdsmn)
688 #endif
689 c
690  return
691  end
692 c-----------------------------------------------------------------------
693  subroutine dssum_msk2(u,mask,binv,nx,ny,nz)
694  include 'SIZE'
695  include 'CTIMER'
696  include 'INPUT'
697  include 'NONCON'
698  include 'PARALLEL'
699  include 'TSTEP'
700  real u(1),mask(1),binv(1)
701 c
702  parameter(lface=lx1*ly1)
703  common /nonctmp/ uin(lface,2*ldim),uout(lface)
704 
705  if(ifsync) call nekgsync()
706 
707 #ifdef TIMER
708  if (icalld.eq.0) then
709  tdsmx=0.
710  tdsmn=0.
711  endif
712  icalld=icalld+1
713  etime1=dnekclock()
714 #endif
715 
716 c
717  ifldt = ifield
718 c if (ifldt.eq.0) ifldt = 1
719  nel = nelv
720  if (ifield.ge.2) nel=nelt
721  ntot = nx*ny*nz*nel
722 
723 
724 c
725 c
726 c T ~ ~T T
727 c Implement Q M Q := J M Q Q J
728 c
729 c
730 c T
731 c This is the J part, translating child data
732 c
733  call apply_jt(u,nx,ny,nz,nel)
734 c
735 c
736 c
737 c ~ ~T
738 c This is the Q Q part
739 c
740  call fgslib_gs_op(gsh_fld(ifldt),u,1,1,0)
741  call col3 (u,mask,binv,ntot)
742 c
743 c
744 c This is the J part, interpolating parent solution onto child
745 c
746  call apply_j(u,nx,ny,nz,nel)
747 c
748 c
749 #ifdef TIMER
750  timee=(dnekclock()-etime1)
751  tdsum=tdsum+timee
752  ndsum=icalld
753  tdsmx=max(timee,tdsmx)
754  tdsmn=min(timee,tdsmn)
755 #endif
756 c
757  return
758  end
759 c-----------------------------------------------------------------------
760  subroutine gtpp_gs_op(u,op,hndl)
761 c
762 c gather-scatter operation across global tensor product planes
763 c
764  include 'SIZE'
765  include 'TOTAL'
766 
767  real u(*)
768  character*3 op
769  integer hndl
770 
771  if (op.eq.'+ ' .or. op.eq.'sum' .or. op.eq.'SUM') then
772  call fgslib_gs_op(hndl,u,1,1,0)
773  elseif (op.eq.'* ' .or. op.eq.'mul' .or. op.eq.'MUL') then
774  call fgslib_gs_op(hndl,u,1,2,0)
775  elseif (op.eq.'m ' .or. op.eq.'min' .or. op.eq.'mna'
776  & .or. op.eq.'MIN' .or. op.eq.'MNA') then
777  call fgslib_gs_op(hndl,u,1,3,0)
778  elseif (op.eq.'M ' .or. op.eq.'max' .or. op.eq.'mxa'
779  & .or. op.eq.'MAX' .or. op.eq.'MXA') then
780  call fgslib_gs_op(hndl,u,1,4,0)
781  else
782  call exitti('gtpp_gs_op: invalid operation!$',1)
783  endif
784 
785  return
786  end
787 c-----------------------------------------------------------------------
788  subroutine gtpp_gs_setup(hndl,nelgx,nelgy,nelgz,idir)
789 
790  include 'SIZE'
791  include 'TOTAL'
792 
793  integer hndl,nelgx,nelgy,nelgz,idir
794 
795  common /nekmpi/ mid,mp,nekcomm,nekgroup,nekreal
796  common /c_is1/ glo_num(lx1,ly1,lz1,lelv)
797  integer e,ex,ey,ez,eg
798  integer*8 glo_num,ex_g
799 
800  nelgxyz = nelgx*nelgy*nelgz
801  if (nelgxyz .ne. nelgv)
802  $ call exitti('gtpp_gs_setup: invalid gtp mesh dimensions!$',
803  $ nelgxyz)
804 
805  nel = nelv
806  nelgxy = nelgx*nelgy
807  nelgyz = nelgy*nelgz
808  nelgzx = nelgz*nelgx
809 
810  if (idir.eq.1) then
811  ! x-direction
812  do e=1,nel
813  eg = lglel(e)
814  call get_exyz(ex,ey,ez,eg,nelgx,nelgyz,1)
815  ex_g = ey
816  do k=1,nz1 ! Enumerate points in the y-z plane
817  do j=1,ny1
818  do i=1,nx1
819  glo_num(i,j,k,e) = j+ny1*(k-1) + ny1*nz1*(ex_g-1)
820  enddo
821  enddo
822  enddo
823  enddo
824  elseif (idir.eq.2) then
825  ! y-direction
826  do e=1,nel
827  eg = lglel(e)
828  call get_exyz(ex,ey,ez,eg,nelgx,nelgy,nelgz)
829  ex_g = (ez-1)*nelgx+ex
830  do k=1,nz1 ! Enumerate points in the x-z plane
831  do j=1,ny1
832  do i=1,nx1
833  glo_num(i,j,k,e) = k+nz1*(i-1) + nx1*nz1*(ex_g-1)
834  enddo
835  enddo
836  enddo
837  enddo
838  elseif (idir.eq.3) then
839  ! z-direction
840  do e=1,nel
841  eg = lglel(e)
842  call get_exyz(ex,ey,ez,eg,nelgxy,1,1)
843  ex_g = ex
844  do k=1,nz1 ! Enumerate points in the x-y plane
845  do j=1,ny1
846  do i=1,nx1
847  glo_num(i,j,k,e) = i+nx1*(j-1) + nx1*ny1*(ex_g-1)
848  enddo
849  enddo
850  enddo
851  enddo
852  endif
853 
854  n = nel*nx1*ny1*nz1
855  call fgslib_gs_setup(hndl,glo_num,n,nekcomm,np)
856 
857  return
858  end
859 c-----------------------------------------------------------------------
860  subroutine gs_setup_ms(hndl,nel,nx,ny,nz)
861 
862  include 'SIZE'
863  include 'TOTAL'
864  include 'mpif.h'
865 
866  integer hndl
867  integer e,eg
868 
869  common /c_is1/ glo_num(lx1*ly1*lz1*lelt)
870  integer*8 glo_num
871 
872  do e=1,nel
873  eg = lglel(e)
874  do k=1,nz
875  do j=1,ny
876  do i=1,nx
877  ii = i + nx*(j-1) + nx*ny*(k-1) + nx*ny*nz*(e-1)
878  glo_num(ii) = i + nx*(j-1) + nx*ny*(k-1) +
879  $ nx*ny*nz*(eg-1)
880  enddo
881  enddo
882  enddo
883  enddo
884 
885  n = nel*nx*ny*nz
886  call fgslib_gs_setup(hndl,glo_num,n,iglobalcomm,np_global)
887 
888  return
889  end
890 c-----------------------------------------------------------------------
891  subroutine gs_op_ms(u,op,hndl)
892 c
893 c gather-scatter operation across sessions
894 c
895  include 'SIZE'
896  include 'PARALLEL'
897  include 'INPUT'
898  include 'TSTEP'
899  include 'CTIMER'
900 
901  real u(*)
902  character*3 op
903 
904  if(ifsync) call nekgsync()
905 
906  if (op.eq.'+ ' .or. op.eq.'sum' .or. op.eq.'SUM') then
907  call fgslib_gs_op(hndl,u,1,1,0)
908  else if (op.eq.'* ' .or. op.eq.'mul' .or. op.eq.'MUL') then
909  call fgslib_gs_op(hndl,u,1,2,0)
910  else if (op.eq.'m ' .or. op.eq.'min' .or. op.eq.'mna'
911  & .or. op.eq.'MIN' .or. op.eq.'MNA') then
912  call fgslib_gs_op(hndl,u,1,3,0)
913  else if (op.eq.'M ' .or. op.eq.'max' .or. op.eq.'mxa'
914  & .or. op.eq.'MAX' .or. op.eq.'MXA') then
915  call fgslib_gs_op(hndl,u,1,4,0)
916  else
917  call exitti('gs_op_ms: invalid operation!$',1)
918  endif
919 
920  return
921  end
nekgsync
subroutine nekgsync()
Definition: comm_mpi.f:502
exitti
subroutine exitti(stringi, idata)
Definition: comm_mpi.f:535
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
vec_dssum
subroutine vec_dssum(u, v, w, nx, ny, nz)
Definition: dssum.f:164
ftovec
subroutine ftovec(a, b, ie, iface, nx, ny, nz)
Definition: dssum.f:417
gs_op_ms
subroutine gs_op_ms(u, op, hndl)
Definition: dssum.f:892
dssum_msk2
subroutine dssum_msk2(u, mask, binv, nx, ny, nz)
Definition: dssum.f:694
apply_j
subroutine apply_j(u, nx, ny, nz, nel)
Definition: dssum.f:532
ftoveci
subroutine ftoveci(a, b, ie, iface, nx, ny, nz)
Definition: dssum.f:453
matvec3t
subroutine matvec3t(uout, Jmat, uin, iftrsp, n1, n2)
Definition: dssum.f:317
nvec_dssum
subroutine nvec_dssum(u, stride, n, gs_handle)
Definition: dssum.f:261
vectof_add
subroutine vectof_add(b, a, ie, iface, nx, ny, nz)
Definition: dssum.f:364
matvec3
subroutine matvec3(uout, Jmat, uin, iftrsp, n1, n2)
Definition: dssum.f:291
gtpp_gs_op
subroutine gtpp_gs_op(u, op, hndl)
Definition: dssum.f:761
zero_f
subroutine zero_f(b, ie, iface, nx, ny, nz)
Definition: dssum.f:382
ftovec_0
subroutine ftovec_0(a, b, ie, iface, nx, ny, nz)
Definition: dssum.f:398
gtpp_gs_setup
subroutine gtpp_gs_setup(hndl, nelgx, nelgy, nelgz, idir)
Definition: dssum.f:789
apply_jt
subroutine apply_jt(u, nx, ny, nz, nel)
Definition: dssum.f:489
h1_proj
subroutine h1_proj(u, nx, ny, nz)
Definition: dssum.f:567
dssum_msk
subroutine dssum_msk(u, mask, nx, ny, nz)
Definition: dssum.f:629
gs_setup_ms
subroutine gs_setup_ms(hndl, nel, nx, ny, nz)
Definition: dssum.f:861
vectof
subroutine vectof(b, a, ie, iface, nx, ny, nz)
Definition: dssum.f:435
dsop
subroutine dsop(u, op, nx, ny, nz)
Definition: dssum.f:101
setupds
subroutine setupds(gs_handle, nx, ny, nz, nel, melg, vertex, glo_num)
Definition: dssum.f:2
vec_dsop
subroutine vec_dsop(u, v, w, nx, ny, nz, op)
Definition: dssum.f:212
dssum
subroutine dssum(u, nx, ny, nz)
Definition: dssum.f:34
matvect
subroutine matvect(out, d, vec, n1, n2)
Definition: dssum.f:340
vectofi
subroutine vectofi(b, a, ie, iface, nx, ny, nz)
Definition: dssum.f:471
col3
subroutine col3(a, b, c, n)
Definition: math.f:598
transpose
subroutine transpose(a, lda, b, ldb)
Definition: math.f:1594
col2
subroutine col2(a, b, n)
Definition: math.f:564
copy
subroutine copy(a, b, n)
Definition: math.f:260
mxm
subroutine mxm(a, n1, b, n2, c, n3)
Definition: mxm_wrapper.f:2
get_exyz
subroutine get_exyz(ex, ey, ez, eg, nelx, nely, nelz)
Definition: navier5.f:1252
set_vert
subroutine set_vert(glo_num, ngv, nx, nel, vertex, ifcenter)
Definition: navier8.f:4