KTH framework for Nek5000 toolboxes; testing version  0.0.1
hsmg.f
Go to the documentation of this file.
1 c-----------------------------------------------------------------------
2 c
3 c Some relevant parameters
4 c
5 c param(41):
6 c 0 - use additive SEMG
7 c 1 - use hybrid SEMG (not yet working... but coming soon!)
8 c
9 c param(42): navier0.f, fasts.f
10 c 0 - use GMRES for iterative solver, use non-symmetric weighting
11 c 1 - use PCG for iterative solver, do not use weighting
12 c
13 c param(43): uzawa_gmres.f, navier6.f
14 c 0 - use additive multilevel scheme (requires param(42).eq.0)
15 c 1 - use original 2 level scheme
16 c
17 c param(44): fast3d.f, navier6.f
18 c 0 - base top-level additive Schwarz on restrictions of E
19 c 1 - base top-level additive Schwarz on restrictions of A
20 c
21 c----------------------------------------------------------------------
22  subroutine hsmg_setup()
23  include 'SIZE'
24  include 'INPUT'
25  include 'PARALLEL'
26  include 'HSMG'
27  include 'SEMHAT'
28  include 'TSTEP'
29 
30  integer nf,nc,nr
31  integer nx,ny,nz
32 
33  mg_fld = 1
34  if (ifield.gt.1) mg_fld = 2
35  if (ifield.eq.1) call hsmg_index_0 ! initialize index sets
36 
37  call hsmg_setup_mg_nx ! set nx values for each level of multigrid
38  call hsmg_setup_semhat ! set spectral element hat matrices
39  call hsmg_setup_intp
40  call hsmg_setup_dssum ! set direct stiffness summation handles
41  call hsmg_setup_wtmask ! set restriction weight matrices and bc masks
42  call hsmg_setup_fdm ! set up fast diagonalization method
43  call hsmg_setup_schwarz_wt(.false.)
44  call hsmg_setup_solve ! set up the solver
45 c call hsmg_setup_dbg
46 
47  return
48  end
49 c----------------------------------------------------------------------
50  subroutine hsmg_setup_semhat
51  include 'SIZE'
52  include 'INPUT'
53  include 'HSMG'
54  include 'SEMHAT'
55  integer n,l
56 c generate the SEM hat matrices for each level
57 c top level
58  n = mg_nx(mg_lmax)
59  call semhat(ah,bh,ch,dh,zh,dph,jph,bgl,zglhat,dgl,jgl,n,wh)
60  call copy(mg_zh(1,mg_lmax),zglhat,n-1) !top level based on gl points
61  mg_nh(mg_lmax)=n-1
62  mg_nhz(mg_lmax)=n-1
63  if(.not.if3d) mg_nhz(mg_lmax)=1
64 c lower levels
65  do l=1,mg_lmax-1
66  n = mg_nx(l)
67  if(n.gt.1) then
68  call semhat(ah,bh,ch,dh,zh,dph,jph,bgl,zglhat,dgl,jgl,n,wh)
69  call copy(mg_ah(1,l),ah,(n+1)*(n+1))
70  call copy(mg_bh(1,l),bh,n+1)
71  call copy(mg_dh(1,l),dh,(n+1)*(n+1))
72  call transpose(mg_dht(1,l),n+1,dh,n+1)
73  call copy(mg_zh(1,l),zh,n+1)
74  else
75  mg_zh(1,l) = -1.
76  mg_zh(2,l) = 1.
77  endif
78  mg_nh(l)=n+1
79  mg_nhz(l)=mg_nz(l)+1
80  enddo
81  end
82 c----------------------------------------------------------------------
83  subroutine hsmg_setup_intp
84  include 'SIZE'
85  include 'HSMG'
86  include 'SEMHAT'
87  integer l,nf,nc
88 
89  do l=1,mg_lmax-1
90 
91  nf=mg_nh(l+1)
92  nc=mg_nh(l)
93 
94 ! Standard multigrid coarse-to-fine interpolation
95  call hsmg_setup_intpm(
96  $ mg_jh(1,l),mg_zh(1,l+1),mg_zh(1,l),nf,nc)
97  call transpose(mg_jht(1,l),nc,mg_jh(1,l),nf)
98 
99 ! Fine-to-coarse interpolation for variable-coefficient operators
100  call hsmg_setup_intpm(
101  $ mg_jhfc(1,l),mg_zh(1,l),mg_zh(1,l+1),nc,nf)
102  call transpose(mg_jhfct(1,l),nf,mg_jhfc(1,l),nc)
103 c call outmat(mg_jhfc(1,l),nc,nf,'MG_JHFC',l)
104 
105  enddo
106  end
107 c----------------------------------------------------------------------
108  subroutine hsmg_setup_intpm(jh,zf,zc,nf,nc)
109  integer nf,nc
110  real jh(nf,nc),zf(1),zc(1)
111  include 'SIZE'
112  real w(2*lx1+2)
113  do i=1,nf
114  call fd_weights_full(zf(i),zc,nc-1,1,w)
115  do j=1,nc
116  jh(i,j)=w(j)
117  enddo
118  enddo
119  return
120  end
121 c----------------------------------------------------------------------
122  subroutine hsmg_setup_dssum
123  include 'SIZE'
124  include 'INPUT'
125  include 'PARALLEL'
126  include 'HSMG'
127  parameter(lxyz=(lx1+2)*(ly1+2)*(lz1+2))
128  common /c_is1/ glo_num(lxyz*lelv)
129  common /ivrtx/ vertex((2**ldim)*lelt)
130 
131  integer*8 glo_num
132  integer vertex
133  integer nx,ny,nz
134  integer l
135 
136 c set up direct stiffness summation for each level
137  ncrnr = 2**ldim
138  call get_vert
139 
140 c++ write(6,*) mg_fld,' mgfld in hsmg_setup_dssum'
141 
142  do l=1,mg_lmax-1
143  nx=mg_nh(l)
144  ny=mg_nh(l)
145  nz=mg_nhz(l)
146  call setupds(mg_gsh_handle(l,mg_fld),nx,ny,nz
147  $ ,nelv,nelgv,vertex,glo_num)
148  nx=nx+2
149  ny=ny+2
150  nz=nz+2
151  if(.not.if3d) nz=1
152  call setupds(mg_gsh_schwarz_handle(l,mg_fld),nx,ny,nz
153  $ ,nelv,nelgv,vertex,glo_num)
154  enddo
155  end
156 c----------------------------------------------------------------------
157  subroutine h1mg_setup_wtmask
158  include 'SIZE'
159  include 'HSMG'
160  integer i,l
161  i = mg_mask_index(mg_lmax,mg_fld-1)
162  do l=1,mg_lmax
163  mg_rstr_wt_index(l,mg_fld)=i
164  mg_mask_index(l,mg_fld)=i
165  i=i+mg_nh(l)*mg_nhz(l)*2*ldim*nelv
166  if(i .gt. lmgs*lmg_rwt*2*ldim*lelv) then
167  itmp = i/(2*ldim*lelv)
168  write(6,*) 'parameter lmg_rwt too small',i,itmp,lmg_rwt
169  call exitt
170  endif
171  call hsmg_setup_rstr_wt(
172  $ mg_rstr_wt(mg_rstr_wt_index(l,mg_fld))
173  $ ,mg_nh(l),mg_nh(l),mg_nhz(l),l,mg_work)
174  call hsmg_setup_mask(
175  $ mg_mask(mg_mask_index(l,mg_fld))
176  $ ,mg_nh(l),mg_nh(l),mg_nhz(l),l,mg_work)
177  enddo
178  mg_mask_index(l,mg_fld)=i
179  end
180 c----------------------------------------------------------------------
181  subroutine hsmg_setup_wtmask
182  include 'SIZE'
183  include 'HSMG'
184  integer i,l
185  i = mg_mask_index(mg_lmax,mg_fld-1)
186  do l=1,mg_lmax-1
187  mg_rstr_wt_index(l,mg_fld)=i
188  mg_mask_index(l,mg_fld)=i
189  i=i+mg_nh(l)*mg_nhz(l)*2*ldim*nelv
190  if(i .gt. lmgs*lmg_rwt*2*ldim*lelv) then
191  itmp = i/(2*ldim*lelv)
192  write(6,*) 'parameter lmg_rwt too small',i,itmp,lmg_rwt
193  call exitt
194  endif
195  call hsmg_setup_rstr_wt(
196  $ mg_rstr_wt(mg_rstr_wt_index(l,mg_fld))
197  $ ,mg_nh(l),mg_nh(l),mg_nhz(l),l,mg_work)
198  call hsmg_setup_mask(
199  $ mg_mask(mg_mask_index(l,mg_fld))
200  $ ,mg_nh(l),mg_nh(l),mg_nhz(l),l,mg_work)
201  enddo
202  mg_mask_index(l,mg_fld)=i
203  end
204 c----------------------------------------------------------------------
205  subroutine hsmg_intp(uf,uc,l) ! l is coarse level
206  real uf(1),uc(1)
207  integer l
208  include 'SIZE'
209  include 'HSMG'
210  call hsmg_tnsr(uf,mg_nh(l+1),uc,mg_nh(l),mg_jh(1,l),mg_jht(1,l))
211  return
212  end
213 c----------------------------------------------------------------------
214  subroutine hsmg_rstr(uc,uf,l) ! l is coarse level
215  real uf(1),uc(1)
216  integer l
217  include 'SIZE'
218  include 'HSMG'
219  if(l.ne.mg_lmax-1)
220  $ call hsmg_do_wt(uf,mg_rstr_wt(mg_rstr_wt_index(l+1,mg_fld))
221  $ ,mg_nh(l+1),mg_nh(l+1),mg_nhz(l+1))
222  call hsmg_tnsr(uc,mg_nh(l),uf,mg_nh(l+1),mg_jht(1,l),mg_jh(1,l))
223  call hsmg_dssum(uc,l)
224  return
225  end
226 c----------------------------------------------------------------------
227  subroutine hsmg_rstr_no_dssum(uc,uf,l) ! l is coarse level
228  real uf(1),uc(1)
229  integer l
230  include 'SIZE'
231  include 'HSMG'
232  if(l.ne.mg_lmax-1)
233  $ call hsmg_do_wt(uf,mg_rstr_wt(mg_rstr_wt_index(l+1,mg_fld))
234  $ ,mg_nh(l+1),mg_nh(l+1),mg_nhz(l+1))
235  call hsmg_tnsr(uc,mg_nh(l),uf,mg_nh(l+1),mg_jht(1,l),mg_jh(1,l))
236  return
237  end
238 c----------------------------------------------------------------------
239 c computes
240 c v = [A (x) A] u or
241 c v = [A (x) A (x) A] u
242  subroutine hsmg_tnsr(v,nv,u,nu,A,At)
243  integer nv,nu
244  real v(1),u(1),A(1),At(1)
245  include 'SIZE'
246  include 'INPUT'
247  if (.not. if3d) then
248  call hsmg_tnsr2d(v,nv,u,nu,a,at)
249  else
250  call hsmg_tnsr3d(v,nv,u,nu,a,at,at)
251  endif
252  return
253  end
254 c----------------------------------------------------------------------
255 c computes
256 c T
257 c v = A u B
258  subroutine hsmg_tnsr2d(v,nv,u,nu,A,Bt)
259  integer nv,nu
260  real v(nv*nv,nelv),u(nu*nu,nelv),A(1),Bt(1)
261  include 'SIZE'
262  common /hsmgw/ work((lx1+2)*(lx1+2))
263  integer ie
264  do ie=1,nelv
265  call mxm(a,nv,u(1,ie),nu,work,nu)
266  call mxm(work,nv,bt,nu,v(1,ie),nv)
267  enddo
268  return
269  end
270 c----------------------------------------------------------------------
271 c computes
272 c
273 c v = [C (x) B (x) A] u
274  subroutine hsmg_tnsr3d(v,nv,u,nu,A,Bt,Ct)
275  integer nv,nu
276  real v(nv*nv*nv,nelv),u(nu*nu*nu,nelv),A(1),Bt(1),Ct(1)
277  include 'SIZE'
278  parameter(lwk=(lx1+2)*(ly1+2)*(lz1+2))
279  common /hsmgw/ work(0:lwk-1),work2(0:lwk-1)
280  integer ie, i
281  do ie=1,nelv
282  call mxm(a,nv,u(1,ie),nu,work,nu*nu)
283  do i=0,nu-1
284  call mxm(work(nv*nu*i),nv,bt,nu,work2(nv*nv*i),nv)
285  enddo
286  call mxm(work2,nv*nv,ct,nu,v(1,ie),nv)
287  enddo
288  return
289  end
290 c----------------------------------------------------------------------
291 c computes
292 c T
293 c v = A u B
294  subroutine hsmg_tnsr2d_el(v,nv,u,nu,A,Bt)
295  integer nv,nu
296  real v(nv*nv),u(nu*nu),A(1),Bt(1)
297  include 'SIZE'
298  common /hsmgw/ work((lx1+2)*(lx1+2))
299 c
300  call mxm(a,nv,u,nu,work,nu)
301  call mxm(work,nv,bt,nu,v,nv)
302 c
303  return
304  end
305 c----------------------------------------------------------------------
306 c computes
307 c
308 c v = [C (x) B (x) A] u
309  subroutine hsmg_tnsr3d_el(v,nv,u,nu,A,Bt,Ct)
310  integer nv,nu
311  real v(nv*nv*nv),u(nu*nu*nu),A(1),Bt(1),Ct(1)
312  include 'SIZE'
313  parameter(lwk=(lx1+2)*(ly1+2)*(lz1+2))
314  common /hsmgw/ work(0:lwk-1),work2(0:lwk-1)
315  integer i
316 c
317  call mxm(a,nv,u,nu,work,nu*nu)
318  do i=0,nu-1
319  call mxm(work(nv*nu*i),nv,bt,nu,work2(nv*nv*i),nv)
320  enddo
321  call mxm(work2,nv*nv,ct,nu,v,nv)
322 c
323  return
324  end
325 c----------------------------------------------------------------------
326  subroutine hsmg_dssum(u,l)
327  include 'SIZE'
328  include 'HSMG'
329  include 'CTIMER'
330  real u(1)
331 
332  if (ifsync) call nekgsync()
333  etime1=dnekclock()
334 
335  call fgslib_gs_op(mg_gsh_handle(l,mg_fld),u,1,1,0)
336  tdadd =tdadd + dnekclock()-etime1
337 
338 
339  return
340  end
341 c----------------------------------------------------------------------
342  subroutine hsmg_dsprod(u,l)
343  include 'SIZE'
344  include 'HSMG'
345  include 'CTIMER'
346  real u(1)
347 
348  if (ifsync) call nekgsync()
349 
350  call fgslib_gs_op(mg_gsh_handle(l,mg_fld),u,1,2,0)
351  return
352  end
353 c----------------------------------------------------------------------
354  subroutine hsmg_schwarz_dssum(u,l)
355  include 'SIZE'
356  include 'HSMG'
357  include 'CTIMER'
358 
359  if (ifsync) call nekgsync()
360  etime1=dnekclock()
361 
362  call fgslib_gs_op(mg_gsh_schwarz_handle(l,mg_fld),u,1,1,0)
363  tdadd =tdadd + dnekclock()-etime1
364 
365  return
366  end
367 c----------------------------------------------------------------------
368  subroutine hsmg_extrude(arr1,l1,f1,arr2,l2,f2,nx,ny,nz)
369  include 'SIZE'
370  include 'INPUT'
371  integer l1,l2,nx,ny,nz
372  real arr1(nx,ny,nz,nelv),arr2(nx,ny,nz,nelv)
373  real f1,f2
374 
375  integer i,j,k,ie,i0,i1
376  i0=2
377  i1=nx-1
378 
379  if(.not.if3d) then
380  do ie=1,nelv
381  do j=i0,i1
382  arr1(l1+1 ,j,1,ie) = f1*arr1(l1+1 ,j,1,ie)
383  $ +f2*arr2(l2+1 ,j,1,ie)
384  arr1(nx-l1,j,1,ie) = f1*arr1(nx-l1,j,1,ie)
385  $ +f2*arr2(nx-l2,j,1,ie)
386  enddo
387  do i=i0,i1
388  arr1(i,l1+1 ,1,ie) = f1*arr1(i,l1+1 ,1,ie)
389  $ +f2*arr2(i,l2+1 ,1,ie)
390  arr1(i,ny-l1,1,ie) = f1*arr1(i,ny-l1,1,ie)
391  $ +f2*arr2(i,nx-l2,1,ie)
392  enddo
393  enddo
394  else
395  do ie=1,nelv
396  do k=i0,i1
397  do j=i0,i1
398  arr1(l1+1 ,j,k,ie) = f1*arr1(l1+1 ,j,k,ie)
399  $ +f2*arr2(l2+1 ,j,k,ie)
400  arr1(nx-l1,j,k,ie) = f1*arr1(nx-l1,j,k,ie)
401  $ +f2*arr2(nx-l2,j,k,ie)
402  enddo
403  enddo
404  do k=i0,i1
405  do i=i0,i1
406  arr1(i,l1+1 ,k,ie) = f1*arr1(i,l1+1 ,k,ie)
407  $ +f2*arr2(i,l2+1 ,k,ie)
408  arr1(i,nx-l1,k,ie) = f1*arr1(i,nx-l1,k,ie)
409  $ +f2*arr2(i,nx-l2,k,ie)
410  enddo
411  enddo
412  do j=i0,i1
413  do i=i0,i1
414  arr1(i,j,l1+1 ,ie) = f1*arr1(i,j,l1+1 ,ie)
415  $ +f2*arr2(i,j,l2+1 ,ie)
416  arr1(i,j,nx-l1,ie) = f1*arr1(i,j,nx-l1,ie)
417  $ +f2*arr2(i,j,nx-l2,ie)
418  enddo
419  enddo
420  enddo
421  endif
422  return
423  end
424 c----------------------------------------------------------------------
425  subroutine h1mg_schwarz(e,r,sigma,l)
426  include 'SIZE'
427  include 'HSMG'
428 
429  real e(1),r(1)
430 
431  n = mg_h1_n(l,mg_fld)
432 
433  call h1mg_schwarz_part1 (e,r,l)
434  call hsmg_schwarz_wt (e,l) ! e := W e
435  call cmult (e,sigma,n) ! l l
436 
437  return
438  end
439 c----------------------------------------------------------------------
440  subroutine h1mg_schwarz_part1 (e,r,l)
441  include 'SIZE'
442  include 'INPUT' ! if3d
443  include 'TSTEP' ! ifield
444  include 'HSMG'
445 
446  real e(1),r(1)
447 
448  integer enx,eny,enz,pm
449 
450  zero = 0
451  one = 1
452  onem = -1
453 
454  n = mg_h1_n(l,mg_fld)
455  pm = p_mg_msk(l,mg_fld)
456 
457  call h1mg_mask (r,mg_imask(pm),nelfld(ifield)) ! Zero Dirichlet nodes
458 
459  if (if3d) then ! extended array
460  call hsmg_schwarz_toext3d(mg_work,r,mg_nh(l))
461  else
462  call hsmg_schwarz_toext2d(mg_work,r,mg_nh(l))
463  endif
464 
465  enx=mg_nh(l)+2
466  eny=mg_nh(l)+2
467  enz=mg_nh(l)+2
468  if(.not.if3d) enz=1
469  i = enx*eny*enz*nelv+1
470 
471 c exchange interior nodes
472  call hsmg_extrude(mg_work,0,zero,mg_work,2,one,enx,eny,enz)
473  call hsmg_schwarz_dssum(mg_work,l)
474  call hsmg_extrude(mg_work,0,one ,mg_work,2,onem,enx,eny,enz)
475 
476  call hsmg_fdm(mg_work(i),mg_work,l) ! Do the local solves
477 
478 c Sum overlap region (border excluded)
479  call hsmg_extrude(mg_work,0,zero,mg_work(i),0,one ,enx,eny,enz)
480  call hsmg_schwarz_dssum(mg_work(i),l)
481  call hsmg_extrude(mg_work(i),0,one ,mg_work,0,onem,enx,eny,enz)
482  call hsmg_extrude(mg_work(i),2,one,mg_work(i),0,one,enx,eny,enz)
483 
484  if(.not.if3d) then ! Go back to regular size array
485  call hsmg_schwarz_toreg2d(e,mg_work(i),mg_nh(l))
486  else
487  call hsmg_schwarz_toreg3d(e,mg_work(i),mg_nh(l))
488  endif
489 
490  call hsmg_dssum(e,l) ! sum border nodes
491  call h1mg_mask (e,mg_imask(pm),nelfld(ifield)) ! apply mask
492 
493  return
494  end
495 c----------------------------------------------------------------------
496  subroutine hsmg_schwarz(e,r,l)
497  include 'SIZE'
498  include 'INPUT'
499  include 'HSMG'
500  real e(1),r(1)
501  integer l
502  integer enx,eny,enz
503  integer i
504 
505  real zero,one,onem
506  zero = 0
507  one = 1
508  onem = -1
509 
510 c apply mask (zeros Dirichlet nodes)
511  !!!!! uncommenting
512  call hsmg_do_wt(r,mg_mask(mg_mask_index(l,mg_fld)),
513  $ mg_nh(l),mg_nh(l),mg_nhz(l))
514 
515 c go to extended size array (room for overlap)
516  if (if3d) then
517  call hsmg_schwarz_toext3d(mg_work,r,mg_nh(l))
518  else
519  call hsmg_schwarz_toext2d(mg_work,r,mg_nh(l))
520  endif
521 
522  enx=mg_nh(l)+2
523  eny=mg_nh(l)+2
524  enz=mg_nh(l)+2
525  if(.not.if3d) enz=1
526  i = enx*eny*enz*nelv+1
527 
528 c exchange interior nodes
529  call hsmg_extrude(mg_work,0,zero,mg_work,2,one,enx,eny,enz)
530  call hsmg_schwarz_dssum(mg_work,l)
531  call hsmg_extrude(mg_work,0,one ,mg_work,2,onem,enx,eny,enz)
532 
533 c do the local solves
534  call hsmg_fdm(mg_work(i),mg_work,l)
535 c sum overlap region (border excluded)
536  call hsmg_extrude(mg_work,0,zero,mg_work(i),0,one ,enx,eny,enz)
537  call hsmg_schwarz_dssum(mg_work(i),l)
538  call hsmg_extrude(mg_work(i),0,one ,mg_work,0,onem,enx,eny,enz)
539  call hsmg_extrude(mg_work(i),2,one,mg_work(i),0,one,enx,eny,enz)
540 c go back to regular size array
541  if(.not.if3d) then
542  call hsmg_schwarz_toreg2d(e,mg_work(i),mg_nh(l))
543  else
544  call hsmg_schwarz_toreg3d(e,mg_work(i),mg_nh(l))
545  endif
546 c sum border nodes
547  call hsmg_dssum(e,l)
548 c apply mask (zeros Dirichlet nodes)
549  !!!!!! changing r to e
550  call hsmg_do_wt(e,mg_mask(mg_mask_index(l,mg_fld)),
551  $ mg_nh(l),mg_nh(l),mg_nhz(l))
552  return
553  end
554 c----------------------------------------------------------------------
555  subroutine hsmg_schwarz_toext2d(a,b,n)
556  include 'SIZE'
557  integer n
558  real a(0:n+1,0:n+1,nelv),b(n,n,nelv)
559 
560  integer i,j,ie
561 c call rzero(a,(n+2)*(n+2)*nelv)
562  do ie=1,nelv
563  do i=0,n+1
564  a(i,0,ie)=0.
565  enddo
566  do j=1,n
567  a(0 ,j,ie)=0.
568  do i=1,n
569  a(i,j,ie)=b(i,j,ie)
570  enddo
571  a(n+1,j,ie)=0.
572  enddo
573  do i=0,n+1
574  a(i,n+1,ie)=0.
575  enddo
576  enddo
577  return
578  end
579 c----------------------------------------------------------------------
580  subroutine hsmg_schwarz_toext3d(a,b,n)
581  include 'SIZE'
582  integer n
583  real a(0:n+1,0:n+1,0:n+1,nelv),b(n,n,n,nelv)
584 
585  integer i,j,k,ie
586  call rzero(a,(n+2)*(n+2)*(n+2)*nelv)
587  do ie=1,nelv
588  do k=1,n
589  do j=1,n
590  do i=1,n
591  a(i,j,k,ie)=b(i,j,k,ie)
592  enddo
593  enddo
594  enddo
595  enddo
596  return
597  end
598 c----------------------------------------------------------------------
599  subroutine hsmg_schwarz_toreg2d(b,a,n)
600  include 'SIZE'
601  integer n
602  real a(0:n+1,0:n+1,nelv),b(n,n,nelv)
603 
604  integer i,j,ie
605  do ie=1,nelv
606  do j=1,n
607  do i=1,n
608  b(i,j,ie)=a(i,j,ie)
609  enddo
610  enddo
611  enddo
612  return
613  end
614 c----------------------------------------------------------------------
615  subroutine hsmg_schwarz_toreg3d(b,a,n)
616  include 'SIZE'
617  integer n
618  real a(0:n+1,0:n+1,0:n+1,nelv),b(n,n,n,nelv)
619 
620  integer i,j,k,ie
621  do ie=1,nelv
622  do k=1,n
623  do j=1,n
624  do i=1,n
625  b(i,j,k,ie)=a(i,j,k,ie)
626  enddo
627  enddo
628  enddo
629  enddo
630  return
631  end
632 c----------------------------------------------------------------------
633  subroutine h1mg_setup_fdm()
634  include 'SIZE'
635  include 'INPUT'
636  include 'HSMG'
637 
638  integer l,i,j,nl
639  i = mg_fast_s_index(mg_lmax,mg_fld-1)
640  j = mg_fast_d_index(mg_lmax,mg_fld-1)
641  do l=2,mg_lmax
642  mg_fast_s_index(l,mg_fld)=i
643  nl = mg_nh(l)+2
644  i=i+nl*nl*2*ldim*nelv
645  if(i .gt. lmg_fasts*2*ldim*lelv) then
646  itmp = i/(2*ldim*lelv)
647  write(6,*) 'lmg_fasts too small',i,itmp,lmg_fasts,l
648  call exitt
649  endif
650  mg_fast_d_index(l,mg_fld)=j
651  j=j+(nl**ldim)*nelv
652  if(j .gt. lmg_fastd*lelv) then
653  itmp = i/(2*ldim*lelv)
654  write(6,*) 'lmg_fastd too small',i,itmp,lmg_fastd,l
655  call exitt
656  endif
657  call hsmg_setup_fast(
658  $ mg_fast_s(mg_fast_s_index(l,mg_fld))
659  $ ,mg_fast_d(mg_fast_d_index(l,mg_fld))
660  $ ,mg_nh(l)+2,mg_ah(1,l),mg_bh(1,l),mg_nx(l))
661  enddo
662  mg_fast_s_index(l,mg_fld)=i
663  mg_fast_d_index(l,mg_fld)=j
664  return
665  end
666 c----------------------------------------------------------------------
667  subroutine hsmg_setup_fdm()
668  include 'SIZE'
669  include 'INPUT'
670  include 'HSMG'
671 
672  integer l,i,j,nl
673  i = mg_fast_s_index(mg_lmax,mg_fld-1)
674  j = mg_fast_d_index(mg_lmax,mg_fld-1)
675  do l=2,mg_lmax-1
676  mg_fast_s_index(l,mg_fld)=i
677  nl = mg_nh(l)+2
678  i=i+nl*nl*2*ldim*nelv
679  if(i .gt. lmg_fasts*2*ldim*lelv) then
680  itmp = i/(2*ldim*lelv)
681  write(6,*) 'lmg_fasts too small',i,itmp,lmg_fasts,l
682  call exitt
683  endif
684  mg_fast_d_index(l,mg_fld)=j
685  j=j+(nl**ldim)*nelv
686  if(j .gt. lmg_fastd*lelv) then
687  itmp = i/(2*ldim*lelv)
688  write(6,*) 'lmg_fastd too small',i,itmp,lmg_fastd,l
689  call exitt
690  endif
691  call hsmg_setup_fast(
692  $ mg_fast_s(mg_fast_s_index(l,mg_fld))
693  $ ,mg_fast_d(mg_fast_d_index(l,mg_fld))
694  $ ,mg_nh(l)+2,mg_ah(1,l),mg_bh(1,l),mg_nx(l))
695  enddo
696  mg_fast_s_index(l,mg_fld)=i
697  mg_fast_d_index(l,mg_fld)=j
698  return
699  end
700 c----------------------------------------------------------------------
701  subroutine hsmg_setup_fast(s,d,nl,ah,bh,n)
702  include 'SIZE'
703  include 'INPUT'
704  include 'HSMG'
705  real s(nl*nl,2,ldim,nelv)
706  real d(nl**ldim,nelv)
707  real ah(1),bh(1)
708  common /ctmpf/ lr(2*lx1+4),ls(2*lx1+4),lt(2*lx1+4)
709  $ , llr(lelt),lls(lelt),llt(lelt)
710  $ , lmr(lelt),lms(lelt),lmt(lelt)
711  $ , lrr(lelt),lrs(lelt),lrt(lelt)
712  real lr ,ls ,lt
713  real llr,lls,llt
714  real lmr,lms,lmt
715  real lrr,lrs,lrt
716 
717  integer i,j,k
718  integer ie,il,nr,ns,nt
719  integer lbr,rbr,lbs,rbs,lbt,rbt,two
720  real eps,diag
721 
722  two = 2
723  ierr = 0
724  do ie=1,nelv
725  call get_fast_bc(lbr,rbr,lbs,rbs,lbt,rbt,ie,two,ierr)
726  nr=nl
727  ns=nl
728  nt=nl
729  call hsmg_setup_fast1d(s(1,1,1,ie),lr,nr,lbr,rbr
730  $ ,llr(ie),lmr(ie),lrr(ie),ah,bh,n,ie)
731  call hsmg_setup_fast1d(s(1,1,2,ie),ls,ns,lbs,rbs
732  $ ,lls(ie),lms(ie),lrs(ie),ah,bh,n,ie)
733  if(if3d) call hsmg_setup_fast1d(s(1,1,3,ie),lt,nt,lbt,rbt
734  $ ,llt(ie),lmt(ie),lrt(ie),ah,bh,n,ie)
735  il=1
736  if(.not.if3d) then
737  eps = 1.e-5*(vlmax(lr(2),nr-2) + vlmax(ls(2),ns-2))
738  do j=1,ns
739  do i=1,nr
740  diag = lr(i)+ls(j)
741  if (diag.gt.eps) then
742  d(il,ie) = 1.0/diag
743  else
744 c write(6,2) ie,'Reset Eig in hsmg setup fast:',i,j,l
745 c $ ,eps,diag,lr(i),ls(j)
746  2 format(i6,1x,a21,3i5,1p4e12.4)
747  d(il,ie) = 0.0
748  endif
749  il=il+1
750  enddo
751  enddo
752  else
753  eps = 1.e-5 * (vlmax(lr(2),nr-2)
754  $ + vlmax(ls(2),ns-2) + vlmax(lt(2),nt-2))
755  do k=1,nt
756  do j=1,ns
757  do i=1,nr
758  diag = lr(i)+ls(j)+lt(k)
759  if (diag.gt.eps) then
760  d(il,ie) = 1.0/diag
761  else
762 c write(6,3) ie,'Reset Eig in hsmg setup fast:',i,j,k,l
763 c $ ,eps,diag,lr(i),ls(j),lt(k)
764  3 format(i6,1x,a21,4i5,1p5e12.4)
765  d(il,ie) = 0.0
766  endif
767  il=il+1
768  enddo
769  enddo
770  enddo
771  endif
772  enddo
773 
774  ierrmx = iglmax(ierr,1)
775  if (ierrmx.gt.0) then
776  if (ierr.gt.0) write(6,*) nid,ierr,' BC FAIL'
777  call exitti('A INVALID BC FOUND in genfast$',ierrmx)
778  endif
779 
780  return
781  end
782 c----------------------------------------------------------------------
783  subroutine hsmg_setup_fast1d(s,lam,nl,lbc,rbc,ll,lm,lr,ah,bh,n,ie)
784  integer nl,lbc,rbc,n
785  real s(nl,nl,2),lam(nl),ll,lm,lr
786  real ah(0:n,0:n),bh(0:n)
787 
788  include 'SIZE'
789  parameter(lxm=lx1+2)
790  common /ctmp0/ b(2*lxm*lxm),w(2*lxm*lxm)
791 
792  call hsmg_setup_fast1d_a(s,lbc,rbc,ll,lm,lr,ah,n)
793  call hsmg_setup_fast1d_b(b,lbc,rbc,ll,lm,lr,bh,n)
794 
795 c if (nid.eq.0) write(6,*) 'THIS is generalev call',nl,lbc
796  call generalev(s,b,lam,nl,w)
797  if(lbc.gt.0) call row_zero(s,nl,nl,1)
798  if(lbc.eq.1) call row_zero(s,nl,nl,2)
799  if(rbc.gt.0) call row_zero(s,nl,nl,nl)
800  if(rbc.eq.1) call row_zero(s,nl,nl,nl-1)
801 
802  call transpose(s(1,1,2),nl,s,nl)
803  return
804  end
805 c----------------------------------------------------------------------
806  subroutine hsmg_setup_fast1d_a(a,lbc,rbc,ll,lm,lr,ah,n)
807  integer lbc,rbc,n
808  real a(0:n+2,0:n+2),ll,lm,lr
809  real ah(0:n,0:n)
810 
811  real fac
812  integer i,j,i0,i1
813  i0=0
814  if(lbc.eq.1) i0=1
815  i1=n
816  if(rbc.eq.1) i1=n-1
817 
818  call rzero(a,(n+3)*(n+3))
819  fac = 2.0/lm
820  a(1,1)=1.0
821  a(n+1,n+1)=1.0
822  do j=i0,i1
823  do i=i0,i1
824  a(i+1,j+1)=fac*ah(i,j)
825  enddo
826  enddo
827  if(lbc.eq.0) then
828  fac = 2.0/ll
829  a(0,0)=fac*ah(n-1,n-1)
830  a(1,0)=fac*ah(n ,n-1)
831  a(0,1)=fac*ah(n-1,n )
832  a(1,1)=a(1,1)+fac*ah(n ,n )
833  else
834  a(0,0)=1.0
835  endif
836  if(rbc.eq.0) then
837  fac = 2.0/lr
838  a(n+1,n+1)=a(n+1,n+1)+fac*ah(0,0)
839  a(n+2,n+1)=fac*ah(1,0)
840  a(n+1,n+2)=fac*ah(0,1)
841  a(n+2,n+2)=fac*ah(1,1)
842  else
843  a(n+2,n+2)=1.0
844  endif
845  return
846  end
847 c----------------------------------------------------------------------
848  subroutine hsmg_setup_fast1d_b(b,lbc,rbc,ll,lm,lr,bh,n)
849  integer lbc,rbc,n
850  real b(0:n+2,0:n+2),ll,lm,lr
851  real bh(0:n)
852 
853  real fac
854  integer i,j,i0,i1
855  i0=0
856  if(lbc.eq.1) i0=1
857  i1=n
858  if(rbc.eq.1) i1=n-1
859 
860  call rzero(b,(n+3)*(n+3))
861  fac = 0.5*lm
862  b(1,1)=1.0
863  b(n+1,n+1)=1.0
864  do i=i0,i1
865  b(i+1,i+1)=fac*bh(i)
866  enddo
867  if(lbc.eq.0) then
868  fac = 0.5*ll
869  b(0,0)=fac*bh(n-1)
870  b(1,1)=b(1,1)+fac*bh(n )
871  else
872  b(0,0)=1.0
873  endif
874  if(rbc.eq.0) then
875  fac = 0.5*lr
876  b(n+1,n+1)=b(n+1,n+1)+fac*bh(0)
877  b(n+2,n+2)=fac*bh(1)
878  else
879  b(n+2,n+2)=1.0
880  endif
881  return
882  end
883 c----------------------------------------------------------------------
884 c clobbers r
885  subroutine hsmg_fdm(e,r,l)
886  include 'SIZE'
887  include 'INPUT'
888  include 'HSMG'
889  call hsmg_do_fast(e,r,
890  $ mg_fast_s(mg_fast_s_index(l,mg_fld)),
891  $ mg_fast_d(mg_fast_d_index(l,mg_fld)),
892  $ mg_nh(l)+2)
893  return
894  end
895 c----------------------------------------------------------------------
896 c clobbers r
897  subroutine hsmg_do_fast(e,r,s,d,nl)
898  include 'SIZE'
899  include 'INPUT'
900  real e(nl**ldim,nelv)
901  real r(nl**ldim,nelv)
902  real s(nl*nl,2,ldim,nelv)
903  real d(nl**ldim,nelv)
904 
905  integer ie,nn,i
906  nn=nl**ldim
907  if(.not.if3d) then
908  do ie=1,nelv
909  call hsmg_tnsr2d_el(e(1,ie),nl,r(1,ie),nl
910  $ ,s(1,2,1,ie),s(1,1,2,ie))
911  do i=1,nn
912  r(i,ie)=d(i,ie)*e(i,ie)
913  enddo
914  call hsmg_tnsr2d_el(e(1,ie),nl,r(1,ie),nl
915  $ ,s(1,1,1,ie),s(1,2,2,ie))
916  enddo
917  else
918  do ie=1,nelv
919  call hsmg_tnsr3d_el(e(1,ie),nl,r(1,ie),nl
920  $ ,s(1,2,1,ie),s(1,1,2,ie),s(1,1,3,ie))
921  do i=1,nn
922  r(i,ie)=d(i,ie)*e(i,ie)
923  enddo
924  call hsmg_tnsr3d_el(e(1,ie),nl,r(1,ie),nl
925  $ ,s(1,1,1,ie),s(1,2,2,ie),s(1,2,3,ie))
926  enddo
927  endif
928  return
929  end
930 c----------------------------------------------------------------------
931 c u = wt .* u
932  subroutine hsmg_do_wt(u,wt,nx,ny,nz)
933  include 'SIZE'
934  include 'INPUT'
935  integer nx,ny,nz
936  real u(nx,ny,nz,nelv)
937  real wt(nx,nz,2,ldim,nelv)
938 
939  integer e
940 
941 c if (nx.eq.2) then
942 c do e=1,nelv
943 c call outmat(wt(1,1,1,1,e),nx,nz,'wt 1-1',e)
944 c call outmat(wt(1,1,2,1,e),nx,nz,'wt 2-1',e)
945 c call outmat(wt(1,1,1,2,e),nx,nz,'wt 1-2',e)
946 c call outmat(wt(1,1,2,2,e),nx,nz,'wt 2-2',e)
947 c enddo
948 c call exitti('hsmg_do_wt quit$',nelv)
949 c endif
950 
951  if (.not. if3d) then
952  do ie=1,nelv
953  do j=1,ny
954  u( 1,j,1,ie)=u( 1,j,1,ie)*wt(j,1,1,1,ie)
955  u(nx,j,1,ie)=u(nx,j,1,ie)*wt(j,1,2,1,ie)
956  enddo
957  do i=2,nx-1
958  u(i, 1,1,ie)=u(i, 1,1,ie)*wt(i,1,1,2,ie)
959  u(i,ny,1,ie)=u(i,ny,1,ie)*wt(i,1,2,2,ie)
960  enddo
961  enddo
962  else
963  do ie=1,nelv
964  do k=1,nz
965  do j=1,ny
966  u( 1,j,k,ie)=u( 1,j,k,ie)*wt(j,k,1,1,ie)
967  u(nx,j,k,ie)=u(nx,j,k,ie)*wt(j,k,2,1,ie)
968  enddo
969  enddo
970  do k=1,nz
971  do i=2,nx-1
972  u(i, 1,k,ie)=u(i, 1,k,ie)*wt(i,k,1,2,ie)
973  u(i,ny,k,ie)=u(i,ny,k,ie)*wt(i,k,2,2,ie)
974  enddo
975  enddo
976  do j=2,ny-1
977  do i=2,nx-1
978  u(i,j, 1,ie)=u(i,j, 1,ie)*wt(i,j,1,3,ie)
979  u(i,j,nz,ie)=u(i,j,nz,ie)*wt(i,j,2,3,ie)
980  enddo
981  enddo
982  enddo
983  endif
984  return
985  end
986 c----------------------------------------------------------------------
987  subroutine hsmg_setup_rstr_wt(wt,nx,ny,nz,l,w)
988  include 'SIZE'
989  include 'INPUT'
990  integer nx,ny,nz,l
991  real w(nx,ny,nz,nelv)
992  real wt(nx,nz,2,ldim,nelv)
993 
994  integer ie
995  !init border nodes to 1
996  call rzero(w,nx*ny*nz*nelv)
997 c print *, 'Setup rstr wt: ',nx,ny,nz,nelv
998  if (.not.if3d) then
999  do ie=1,nelv
1000  do i=1,nx
1001  w(i,1,1,ie)=1.0
1002  w(i,ny,1,ie)=1.0
1003  enddo
1004  do j=1,ny
1005  w(1,j,1,ie)=1.0
1006  w(nx,j,1,ie)=1.0
1007  enddo
1008  enddo
1009  else
1010  do ie=1,nelv
1011  do j=1,ny
1012  do i=1,nx
1013  w(i,j,1,ie)=1.0
1014  w(i,j,nz,ie)=1.0
1015  enddo
1016  enddo
1017  do k=1,nz
1018  do i=1,nx
1019  w(i,1,k,ie)=1.0
1020  w(i,ny,k,ie)=1.0
1021  enddo
1022  enddo
1023  do k=1,nz
1024  do j=1,ny
1025  w(1,j,k,ie)=1.0
1026  w(nx,j,k,ie)=1.0
1027  enddo
1028  enddo
1029  enddo
1030  endif
1031  call hsmg_dssum(w,l)
1032  !invert the count w to get the weight wt
1033  if (.not. if3d) then
1034  do ie=1,nelv
1035  do j=1,ny
1036  wt(j,1,1,1,ie)=1.0/w(1,j,1,ie)
1037  wt(j,1,2,1,ie)=1.0/w(nx,j,1,ie)
1038  enddo
1039  do i=1,nx
1040  wt(i,1,1,2,ie)=1.0/w(i,1,1,ie)
1041  wt(i,1,2,2,ie)=1.0/w(i,ny,1,ie)
1042  enddo
1043  enddo
1044  else
1045  do ie=1,nelv
1046  do k=1,nz
1047  do j=1,ny
1048  wt(j,k,1,1,ie)=1.0/w(1,j,k,ie)
1049  wt(j,k,2,1,ie)=1.0/w(nx,j,k,ie)
1050  enddo
1051  enddo
1052  do k=1,nz
1053  do i=1,nx
1054  wt(i,k,1,2,ie)=1.0/w(i,1,k,ie)
1055  wt(i,k,2,2,ie)=1.0/w(i,ny,k,ie)
1056  enddo
1057  enddo
1058  do j=1,ny
1059  do i=1,nx
1060  wt(i,j,1,3,ie)=1.0/w(i,j,1,ie)
1061  wt(i,j,2,3,ie)=1.0/w(i,j,nz,ie)
1062  enddo
1063  enddo
1064  enddo
1065  endif
1066  return
1067  end
1068 c----------------------------------------------------------------------
1069  subroutine hsmg_setup_mask(wt,nx,ny,nz,l,w)
1070  include 'SIZE'
1071  include 'INPUT'
1072  integer nx,ny,nz,l
1073  real w(nx,ny,nz,nelv)
1074  real wt(nx,nz,2,ldim,nelv)
1075 
1076  integer ie
1077  integer lbr,rbr,lbs,rbs,lbt,rbt,two
1078 c init everything to 1
1079 
1080  n = nx*ny*nz*nelv
1081  call rone(w,n)
1082 
1083 c set dirichlet nodes to zero
1084  ierr = 0
1085  two = 2
1086  do ie=1,nelv
1087  call get_fast_bc(lbr,rbr,lbs,rbs,lbt,rbt,ie,two,ierr)
1088  if (ierr.ne.0) then
1089  ierr = -1
1090  call get_fast_bc(lbr,rbr,lbs,rbs,lbt,rbt,ie,two,ierr)
1091  endif
1092 
1093  if(lbr.eq.1) then
1094  do k=1,nz
1095  do j=1,ny
1096  w(1,j,k,ie)=0.0
1097  enddo
1098  enddo
1099  endif
1100  if(rbr.eq.1) then
1101  do k=1,nz
1102  do j=1,ny
1103  w(nx,j,k,ie)=0.0
1104  enddo
1105  enddo
1106  endif
1107  if(lbs.eq.1) then
1108  do k=1,nz
1109  do i=1,nx
1110  w(i,1,k,ie)=0.0
1111  enddo
1112  enddo
1113  endif
1114  if(rbs.eq.1) then
1115  do k=1,nz
1116  do i=1,nx
1117  w(i,ny,k,ie)=0.0
1118  enddo
1119  enddo
1120  endif
1121  if(if3d) then
1122  if(lbt.eq.1) then
1123  do j=1,ny
1124  do i=1,nx
1125  w(i,j,1,ie)=0.0
1126  enddo
1127  enddo
1128  endif
1129  if(rbt.eq.1) then
1130  do j=1,ny
1131  do i=1,nx
1132  w(i,j,nz,ie)=0.0
1133  enddo
1134  enddo
1135  endif
1136  endif
1137  enddo
1138 c do direct stiffness multiply
1139 
1140  call hsmg_dsprod(w,l)
1141 
1142 
1143 c store weight
1144  if (.not. if3d) then
1145  do ie=1,nelv
1146  do j=1,ny
1147  wt(j,1,1,1,ie)=w(1,j,1,ie)
1148  wt(j,1,2,1,ie)=w(nx,j,1,ie)
1149  enddo
1150  do i=1,nx
1151  wt(i,1,1,2,ie)=w(i,1,1,ie)
1152  wt(i,1,2,2,ie)=w(i,ny,1,ie)
1153  enddo
1154  enddo
1155  else
1156  do ie=1,nelv
1157  do k=1,nz
1158  do j=1,ny
1159  wt(j,k,1,1,ie)=w(1,j,k,ie)
1160  wt(j,k,2,1,ie)=w(nx,j,k,ie)
1161  enddo
1162  enddo
1163  do k=1,nz
1164  do i=1,nx
1165  wt(i,k,1,2,ie)=w(i,1,k,ie)
1166  wt(i,k,2,2,ie)=w(i,ny,k,ie)
1167  enddo
1168  enddo
1169  do k=1,nz
1170  do j=1,ny
1171  wt(j,k,1,3,ie)=w(i,j,1,ie)
1172  wt(j,k,2,3,ie)=w(i,j,nz,ie)
1173  enddo
1174  enddo
1175  enddo
1176  endif
1177 
1178  ierrmx = iglmax(ierr,1)
1179  if (ierrmx.gt.0) then
1180  if (ierr.gt.0) write(6,*) nid,ierr,' BC FAIL'
1181  call exitti('B INVALID BC FOUND in genfast$',ierrmx)
1182  endif
1183 
1184  return
1185  end
1186 c----------------------------------------------------------------------
1187  subroutine hsmg_setup_schwarz_wt(ifsqrt)
1188  logical ifsqrt
1189  include 'SIZE'
1190  include 'INPUT'
1191  include 'HSMG'
1192 
1193  integer l,i,nl,nlz
1194 
1195  i = mg_schwarz_wt_index(mg_lmax,mg_fld-1)
1196  do l=2,mg_lmax-1
1197  mg_schwarz_wt_index(l,mg_fld)=i
1198  nl = mg_nh(l)
1199  nlz = mg_nh(l)
1200  if(.not.if3d) nlz=1
1201  i=i+nl*nlz*4*ldim*nelv
1202  if(i .gt. lmg_swt*4*ldim*lelv) then
1203  itmp = i/(4*ldim*lelv)
1204  write(6,*) 'lmg_swt too small',i,itmp,lmg_swt,l
1205  call exitt
1206  endif
1207 
1208  call h1mg_setup_schwarz_wt_1(
1209  $ mg_schwarz_wt(mg_schwarz_wt_index(l,mg_fld)),l,ifsqrt)
1210 
1211  enddo
1212  mg_schwarz_wt_index(l,mg_fld)=i
1213 
1214  return
1215  end
1216 c----------------------------------------------------------------------
1217  subroutine h1mg_setup_schwarz_wt(ifsqrt)
1218  logical ifsqrt
1219  include 'SIZE'
1220  include 'INPUT'
1221  include 'HSMG'
1222 
1223  integer l,i,nl,nlz
1224 
1225  i = mg_schwarz_wt_index(mg_lmax,mg_fld-1)
1226  do l=2,mg_lmax
1227 
1228  mg_schwarz_wt_index(l,mg_fld)=i
1229  nl = mg_nh(l)
1230  nlz = mg_nhz(l)
1231  i = i+nl*nlz*4*ldim*nelv
1232 
1233  if (i .gt. lmg_swt*4*ldim*lelv) then
1234  itmp = i/(4*ldim*lelv)
1235  write(6,*) 'lmg_swt too small',i,itmp,lmg_swt,l
1236  call exitt
1237  endif
1238 
1239  call h1mg_setup_schwarz_wt_1(
1240  $ mg_schwarz_wt(mg_schwarz_wt_index(l,mg_fld)),l,ifsqrt)
1241 
1242  enddo
1243 
1244  mg_schwarz_wt_index(l,mg_fld)=i
1245 
1246  return
1247  end
1248 c----------------------------------------------------------------------
1249  subroutine hsmg_schwarz_wt(e,l)
1250  include 'SIZE'
1251  include 'INPUT'
1252  include 'HSMG'
1253 
1254  if(.not.if3d) call hsmg_schwarz_wt2d(
1255  $ e,mg_schwarz_wt(mg_schwarz_wt_index(l,mg_fld)),mg_nh(l))
1256  if(if3d) call hsmg_schwarz_wt3d(
1257  $ e,mg_schwarz_wt(mg_schwarz_wt_index(l,mg_fld)),mg_nh(l))
1258  return
1259  end
1260 c----------------------------------------------------------------------
1261  subroutine hsmg_schwarz_wt2d(e,wt,n)
1262  include 'SIZE'
1263  integer n
1264  real e(n,n,nelv)
1265  real wt(n,4,2,nelv)
1266 
1267  integer ie,i,j
1268  do ie=1,nelv
1269  do j=1,n
1270  e(1 ,j,ie)=e(1 ,j,ie)*wt(j,1,1,ie)
1271  e(2 ,j,ie)=e(2 ,j,ie)*wt(j,2,1,ie)
1272  e(n-1,j,ie)=e(n-1,j,ie)*wt(j,3,1,ie)
1273  e(n ,j,ie)=e(n ,j,ie)*wt(j,4,1,ie)
1274  enddo
1275  do i=3,n-2
1276  e(i,1 ,ie)=e(i,1 ,ie)*wt(i,1,2,ie)
1277  e(i,2 ,ie)=e(i,2 ,ie)*wt(i,2,2,ie)
1278  e(i,n-1,ie)=e(i,n-1,ie)*wt(i,3,2,ie)
1279  e(i,n ,ie)=e(i,n ,ie)*wt(i,4,2,ie)
1280  enddo
1281  enddo
1282  return
1283  end
1284 c----------------------------------------------------------------------
1285  subroutine hsmg_schwarz_wt3d(e,wt,n)
1286  include 'SIZE'
1287  integer n
1288  real e(n,n,n,nelv)
1289  real wt(n,n,4,3,nelv)
1290 
1291  integer ie,i,j,k
1292  do ie=1,nelv
1293  do k=1,n
1294  do j=1,n
1295  e(1 ,j,k,ie)=e(1 ,j,k,ie)*wt(j,k,1,1,ie)
1296  e(2 ,j,k,ie)=e(2 ,j,k,ie)*wt(j,k,2,1,ie)
1297  e(n-1,j,k,ie)=e(n-1,j,k,ie)*wt(j,k,3,1,ie)
1298  e(n ,j,k,ie)=e(n ,j,k,ie)*wt(j,k,4,1,ie)
1299  enddo
1300  enddo
1301  do k=1,n
1302  do i=3,n-2
1303  e(i,1 ,k,ie)=e(i,1 ,k,ie)*wt(i,k,1,2,ie)
1304  e(i,2 ,k,ie)=e(i,2 ,k,ie)*wt(i,k,2,2,ie)
1305  e(i,n-1,k,ie)=e(i,n-1,k,ie)*wt(i,k,3,2,ie)
1306  e(i,n ,k,ie)=e(i,n ,k,ie)*wt(i,k,4,2,ie)
1307  enddo
1308  enddo
1309  do j=3,n-2
1310  do i=3,n-2
1311  e(i,j,1 ,ie)=e(i,j,1 ,ie)*wt(i,j,1,3,ie)
1312  e(i,j,2 ,ie)=e(i,j,2 ,ie)*wt(i,j,2,3,ie)
1313  e(i,j,n-1,ie)=e(i,j,n-1,ie)*wt(i,j,3,3,ie)
1314  e(i,j,n ,ie)=e(i,j,n ,ie)*wt(i,j,4,3,ie)
1315  enddo
1316  enddo
1317  enddo
1318  return
1319  end
1320 c----------------------------------------------------------------------
1321  subroutine hsmg_coarse_solve(e,r)
1322  include 'SIZE'
1323  include 'DOMAIN'
1324  include 'ESOLV'
1325  include 'GEOM'
1326  include 'SOLN'
1327  include 'PARALLEL'
1328  include 'HSMG'
1329  include 'CTIMER'
1330  include 'INPUT'
1331  include 'TSTEP'
1332  real e(1),r(1)
1333 c
1334  integer n_crs_tot
1335  save n_crs_tot
1336  data n_crs_tot /0/
1337 c
1338  if (icalld.eq.0) then ! timer info
1339  ncrsl=0
1340  tcrsl=0.0
1341  endif
1342  icalld = 1
1343 
1344  if (ifsync) call nekgsync()
1345 
1346  ncrsl = ncrsl + 1
1347  etime1=dnekclock()
1348 
1349  call fgslib_crs_solve(xxth(ifield),e,r)
1350 
1351  tcrsl=tcrsl+dnekclock()-etime1
1352 
1353  return
1354  end
1355 c----------------------------------------------------------------------
1356  subroutine hsmg_setup_solve
1357  include 'SIZE'
1358  include 'HSMG'
1359 
1360  integer l,i,nl,nlz
1361  i = mg_solve_index(mg_lmax+1,mg_fld-1)
1362  do l=1,mg_lmax
1363  mg_solve_index(l,mg_fld)=i
1364  i=i+mg_nh(l)*mg_nh(l)*mg_nhz(l)*nelv
1365  if(i .gt. lmg_solve*lelv) then
1366  itmp = i/lelv
1367  write(6,*) 'lmg_solve too small',i,itmp,lmg_solve,l
1368  call exitt
1369  endif
1370  enddo
1371  mg_solve_index(l,mg_fld)=i
1372 
1373  return
1374  end
1375 c----------------------------------------------------------------------
1376  subroutine hsmg_solve(e,r)
1377  real e(1),r(1)
1378  include 'SIZE'
1379  include 'HSMG'
1380  include 'GEOM'
1381  include 'INPUT'
1382  include 'MASS'
1383  include 'SOLN'
1384  include 'TSTEP'
1385  include 'CTIMER'
1386  include 'PARALLEL'
1387 
1388  common /quick/ ecrs(2) ! quick work array
1389  $ , ecrs2(2) ! quick work array
1390 c common /quick/ ecrs (lx2*ly2*lz2*lelv) ! quick work array
1391 c $ , ecrs2 (lx2*ly2*lz2*lelv) ! quick work array
1392 
1393  common /scrhi/ h2inv(lx1,ly1,lz1,lelv)
1394  common /scrvh/ h1(lx1,ly1,lz1,lelv),
1395  $ h2(lx1,ly1,lz1,lelv)
1396 
1397  integer ilstep,iter
1398  save ilstep,iter
1399  data ilstep,iter /0,0/
1400 
1401  real rhoavg,copt(2),copw(2)
1402  save rhoavg,copt1,copt2
1403  data rhoavg,copt1,copt2 /3*1./ ! Default copt = 1 for additive
1404 
1405  integer l,nt
1406  integer*8 ntotg,nxyz2
1407 
1408  logical if_hybrid
1409 
1410  mg_fld = 1
1411  if (ifield.gt.1) mg_fld = 2
1412 
1413  if (istep.ne.ilstep) then
1414  ilstep = istep
1415  ntot1 = lx1*ly1*lz1*nelv
1416  rhoavg = glsc2(vtrans,bm1,ntot1)/volvm1
1417  endif
1418 
1419  n = lx2*ly2*lz2*nelv
1420 c call copy(e,r,n)
1421 c return
1422 
1423  if (icalld.eq.0) then
1424 
1425  tddsl=0.0
1426  nddsl=0
1427 
1428  icalld = 1
1429  taaaa = 0
1430  tbbbb = 0
1431  tcccc = 0
1432  tdddd = 0
1433  teeee = 0
1434  endif
1435 
1436  nddsl = nddsl + 1
1437  etime1 = dnekclock()
1438 
1439 
1440 c n = lx2*ly2*lz2*nelv
1441 c rmax = glmax(r,n)
1442 c if (nid.eq.0) write(6,*) istep,n,rmax,' rmax1'
1443 
1444  iter = iter + 1
1445 
1446  l = mg_lmax
1447  nt = mg_nh(l)*mg_nh(l)*mg_nhz(l)*nelv
1448  ! e := W M r
1449  ! Schwarz
1450  time_0 = dnekclock()
1451  call local_solves_fdm(e,r)
1452 
1453  time_1 = dnekclock()
1454 
1455 c if (param(41).eq.1) if_hybrid = .true.
1456  if_hybrid = .false.
1457 
1458  if (if_hybrid) then
1459  ! w := E e
1460  rbd1dt = rhoavg*bd(1)/dt ! Assumes constant density!!!
1461  call cdabdtp(mg_work2,e,h1,h2,h2inv,1)
1462  call cmult (mg_work2,rbd1dt,nt)
1463  time_2 = dnekclock()
1464  if (istep.eq.1) then
1465  copt(1) = vlsc2(r ,mg_work2,nt)
1466  copt(2) = vlsc2(mg_work2,mg_work2,nt)
1467  call gop(copt,copw,'+ ', 2)
1468  copt(1) = copt(1)/copt(2)
1469  avg2 = 1./iter
1470  avg1 = 1.-avg2
1471  copt1 = avg1*copt1 + avg2*copt(1)
1472  if(nio.eq.0)write(6,1)istep,iter,rbd1dt,copt(1),copt1,'cpt1'
1473  1 format(2i6,1p3e14.5,2x,a4)
1474  endif
1475  ! w := r - w
1476  do i = 1,nt
1477  mg_work2(i) = r(i) - copt1*mg_work2(i)
1478  e(i) = copt1*e(i)
1479  ecrs2(i) = mg_work2(i)
1480  enddo
1481 
1482  else ! Additive
1483  ! w := r - w
1484  do i = 1,nt
1485  mg_work2(i) = r(i)
1486  enddo
1487  time_2 = dnekclock()
1488  endif
1489 
1490  do l = mg_lmax-1,2,-1
1491 
1492 c rmax = glmax(mg_work2,nt)
1493 c if (nid.eq.0) write(6,*) l,nt,rmax,' rmax2'
1494 
1495  nt = mg_nh(l)*mg_nh(l)*mg_nhz(l)*nelv
1496  ! T
1497  ! r := J w
1498  ! l
1499  call hsmg_rstr(mg_solve_r(mg_solve_index(l,mg_fld)),mg_work2,l)
1500 
1501  ! w := r
1502  ! l
1503  call copy(mg_work2,mg_solve_r(mg_solve_index(l,mg_fld)),nt)
1504  ! e := M w
1505  ! l Schwarz
1506  call hsmg_schwarz(
1507  $ mg_solve_e(mg_solve_index(l,mg_fld)),mg_work2,l)
1508 
1509  ! e := W e
1510  ! l l
1511  call hsmg_schwarz_wt(mg_solve_e(mg_solve_index(l,mg_fld)),l)
1512 
1513 c call exitti('quit in mg$',l)
1514 
1515  ! w := r - w
1516  ! l
1517  do i = 0,nt-1
1518  mg_work2(i+1) = mg_solve_r(mg_solve_index(l,mg_fld)+i)
1519  $ !-alpha*mg_work2(i+1)
1520  enddo
1521  enddo
1522 
1523  call hsmg_rstr_no_dssum(
1524  $ mg_solve_r(mg_solve_index(1,mg_fld)),mg_work2,1)
1525 
1526  nzw = ldim-1
1527 
1528  call hsmg_do_wt(mg_solve_r(mg_solve_index(1,mg_fld)),
1529  $ mg_mask(mg_mask_index(1,mg_fld)),2,2,nzw)
1530 
1531  ! -1
1532  ! e := A r
1533  ! 1 1
1534  call hsmg_coarse_solve(mg_solve_e(mg_solve_index(1,mg_fld)),
1535  $ mg_solve_r(mg_solve_index(1,mg_fld)))
1536 
1537  call hsmg_do_wt(mg_solve_e(mg_solve_index(1,mg_fld)),
1538  $ mg_mask(mg_mask_index(1,mg_fld)),2,2,nzw)
1539  time_3 = dnekclock()
1540  do l = 2,mg_lmax-1
1541  nt = mg_nh(l)*mg_nh(l)*mg_nhz(l)*nelv
1542  ! w := J e
1543  ! l-1
1544  call hsmg_intp
1545  $ (mg_work2,mg_solve_e(mg_solve_index(l-1,mg_fld)),l-1)
1546 
1547  ! e := e + w
1548  ! l l
1549  do i = 0,nt-1
1550  mg_solve_e(mg_solve_index(l,mg_fld)+i) =
1551  $ + mg_solve_e(mg_solve_index(l,mg_fld)+i) + mg_work2(i+1)
1552  enddo
1553  enddo
1554  l = mg_lmax
1555  nt = mg_nh(l)*mg_nh(l)*mg_nhz(l)*nelv
1556  ! w := J e
1557  ! m-1
1558 
1559  call hsmg_intp(mg_work2,
1560  $ mg_solve_e(mg_solve_index(l-1,mg_fld)),l-1)
1561 
1562  if (if_hybrid.and.istep.eq.1) then
1563  ! ecrs := E e_c
1564  call cdabdtp(ecrs,mg_work2,h1,h2,h2inv,1)
1565  call cmult (ecrs,rbd1dt,nt)
1566  copt(1) = vlsc2(ecrs2,ecrs,nt)
1567  copt(2) = vlsc2(ecrs ,ecrs,nt)
1568  call gop(copt,copw,'+ ', 2)
1569  copt(1) = copt(1)/copt(2)
1570  avg2 = 1./iter
1571  avg1 = 1.-avg2
1572  copt2 = avg1*copt2 + avg2*copt(1)
1573  if(nio.eq.0)write(6,1)istep,iter,rbd1dt,copt(1),copt2,'cpt2'
1574  endif
1575  ! e := e + w
1576 
1577  do i = 1,nt
1578  e(i) = e(i) + copt2*mg_work2(i)
1579  enddo
1580  time_4 = dnekclock()
1581 c print *, 'Did an MG iteration'
1582 c
1583  taaaa = taaaa + (time_1 - time_0)
1584  tbbbb = tbbbb + (time_2 - time_1)
1585  tcccc = tcccc + (time_3 - time_2)
1586  tdddd = tdddd + (time_4 - time_3)
1587  teeee = teeee + (time_4 - time_0)
1588 c
1589 c A typical time breakdown:
1590 c
1591 c 1.3540E+01 5.4390E+01 1.1440E+01 1.2199E+00 8.0590E+01 HSMG time
1592 c
1593 c ==> 54/80 = 67 % of preconditioner time is in residual evaluation!
1594 c
1595  call ortho (e)
1596 
1597  tddsl = tddsl + ( dnekclock()-etime1 )
1598 
1599 
1600 
1601  return
1602  end
1603 c----------------------------------------------------------------------
1604  subroutine hsmg_setup_mg_nx()
1605  include 'SIZE'
1606  include 'INPUT'
1607  include 'PARALLEL'
1608  include 'HSMG'
1609  include 'SEMHAT'
1610  real w(lx1+2)
1611  integer nf,nc,nr
1612  integer nx,ny,nz
1613 
1614  integer mgn2(10)
1615  save mgn2
1616  data mgn2 / 1, 2, 2, 2, 2, 3, 3, 5, 5, 5/
1617 c data mgn2 / 1, 2, 3, 4, 5, 6, 7, 8, 9, 0
1618 
1619 c if (param(82).eq.0) param(82)=2 ! nek default
1620 c if (np.eq.1) param(82)=2 ! single proc. too slow
1621  p82 = 2 ! potentially variable nxc
1622  mg_lmax = 3
1623 c mg_lmax = 4
1624  if (lx1.eq.4) mg_lmax = 2
1625 c if (param(79).ne.0) mg_lmax = param(79)
1626  mglx1 = p82-1 !1
1627  mg_nx(1) = mglx1
1628  mg_ny(1) = mglx1
1629  mg_nz(1) = mglx1
1630  if (.not.if3d) mg_nz(1) = 0
1631 
1632  mglx2 = 2*(lx2/4) + 1
1633  if (lx1.eq.5) mglx2 = 3
1634 c if (lx1.eq.6) mglx2 = 3
1635  if (lx1.le.10) mglx2 = mgn2(lx1)
1636  if (lx1.eq.8) mglx2 = 4
1637  if (lx1.eq.8) mglx2 = 3
1638 
1639 c mglx2 = min(3,mglx2)
1640 
1641 
1642  mg_nx(2) = mglx2
1643  mg_ny(2) = mglx2
1644  mg_nz(2) = mglx2
1645  if (.not.if3d) mg_nz(2) = 0
1646 
1647  mg_nx(3) = mglx2+1
1648  mg_ny(3) = mglx2+1
1649  mg_nz(3) = mglx2+1
1650  if (.not.if3d) mg_nz(3) = 0
1651 
1652  mg_nx(mg_lmax) = lx1-1
1653  mg_ny(mg_lmax) = ly1-1
1654  mg_nz(mg_lmax) = lz1-1
1655 
1656  if (nio.eq.0) write(*,*) 'mg_nx:',(mg_nx(i),i=1,mg_lmax)
1657  if (nio.eq.0) write(*,*) 'mg_ny:',(mg_ny(i),i=1,mg_lmax)
1658  if (nio.eq.0) write(*,*) 'mg_nz:',(mg_nz(i),i=1,mg_lmax)
1659 
1660  return
1661  end
1662 c----------------------------------------------------------------------
1663  subroutine hsmg_index_0 ! initialize index sets
1664  include 'SIZE'
1665  include 'HSMG'
1666 
1667  n = lmgn*(lmgs+1)
1668 
1669  call izero( mg_rstr_wt_index , n )
1670  call izero( mg_mask_index , n )
1671  call izero( mg_solve_index , n )
1672  call izero( mg_fast_s_index , n )
1673  call izero( mg_fast_d_index , n )
1674  call izero( mg_schwarz_wt_index , n )
1675 
1676  return
1677  end
1678 c----------------------------------------------------------------------
1679  subroutine outfldn (x,n,txt10,ichk) ! writes into unit=40+ifiled
1680  include 'SIZE'
1681  include 'TSTEP'
1682  real x(n,n,1,lelt)
1683  character*10 txt10
1684 c
1685  integer idum,e
1686  save idum
1687  data idum /3/
1688  if (idum.lt.0) return
1689  m = 40 + ifield ! unit #
1690 c
1691 C
1692  mtot = n*n*nelv
1693  if (n.gt.7.or.nelv.gt.16) return
1694  xmin = glmin(x,mtot)
1695  xmax = glmax(x,mtot)
1696 c
1697  rnel = nelv
1698  snel = sqrt(rnel)+.1
1699  ne = snel
1700  ne1 = nelv-ne+1
1701  do ie=ne1,1,-ne
1702  l=ie-1
1703  do k=1,1
1704  if (ie.eq.ne1) write(m,116) txt10,k,ie,xmin,xmax,ichk,time
1705  write(m,117)
1706  do j=n,1,-1
1707  if (n.eq.2) write(m,102) ((x(i,j,k,e+l),i=1,n),e=1,ne)
1708  if (n.eq.3) write(m,103) ((x(i,j,k,e+l),i=1,n),e=1,ne)
1709  if (n.eq.4) write(m,104) ((x(i,j,k,e+l),i=1,n),e=1,ne)
1710  if (n.eq.5) write(m,105) ((x(i,j,k,e+l),i=1,n),e=1,ne)
1711  if (n.eq.6) write(m,106) ((x(i,j,k,e+l),i=1,n),e=1,ne)
1712  if (n.eq.7) write(m,107) ((x(i,j,k,e+l),i=1,n),e=1,ne)
1713  if (n.eq.8) write(m,108) ((x(i,j,k,e+l),i=1,n),e=1,ne)
1714  enddo
1715  enddo
1716  enddo
1717 
1718 C
1719  102 FORMAT(4(2f9.5,2x))
1720  103 FORMAT(4(3f9.5,2x))
1721  104 FORMAT(4(4f7.3,2x))
1722  105 FORMAT(5f9.5,10x,5f9.5)
1723  106 FORMAT(6f9.5,5x,6f9.5)
1724  107 FORMAT(7f8.4,5x,7f8.4)
1725  108 FORMAT(8f8.4,4x,8f8.4)
1726 c
1727  116 FORMAT( /,5x,' ^ ',/,
1728  $ 5x,' Y | ',/,
1729  $ 5x,' | ',a10,/,
1730  $ 5x,' +----> ','Plane = ',i2,'/',i2,2x,2e12.4,/,
1731  $ 5x,' X ','Step =',i9,f15.5)
1732  117 FORMAT(' ')
1733 c
1734 c if (ichk.eq.1.and.idum.gt.0) call checkit(idum)
1735  return
1736  end
1737 c-----------------------------------------------------------------------
1738  subroutine outfldn0 (x,n,txt10,ichk)
1739  include 'SIZE'
1740  include 'TSTEP'
1741  real x(n,n,1,lelt)
1742  character*10 txt10
1743 c
1744  integer idum,e
1745  save idum
1746  data idum /3/
1747  if (idum.lt.0) return
1748 c
1749 C
1750  mtot = n*n*nelv
1751  if (n.gt.7.or.nelv.gt.16) return
1752  xmin = glmin(x,mtot)
1753  xmax = glmax(x,mtot)
1754 c
1755  rnel = nelv
1756  snel = sqrt(rnel)+.1
1757  ne = snel
1758  ne1 = nelv-ne+1
1759  do ie=ne1,1,-ne
1760  l=ie-1
1761  do k=1,1
1762  if (ie.eq.ne1) write(6,116) txt10,k,ie,xmin,xmax,ichk,time
1763  write(6,117)
1764  do j=n,1,-1
1765  if (n.eq.2) write(6,102) ((x(i,j,k,e+l),i=1,n),e=1,ne)
1766  if (n.eq.3) write(6,103) ((x(i,j,k,e+l),i=1,n),e=1,ne)
1767  if (n.eq.4) write(6,104) ((x(i,j,k,e+l),i=1,n),e=1,ne)
1768  if (n.eq.5) write(6,105) ((x(i,j,k,e+l),i=1,n),e=1,ne)
1769  if (n.eq.6) write(6,106) ((x(i,j,k,e+l),i=1,n),e=1,ne)
1770  if (n.eq.7) write(6,107) ((x(i,j,k,e+l),i=1,n),e=1,ne)
1771  if (n.eq.8) write(6,108) ((x(i,j,k,e+l),i=1,n),e=1,ne)
1772  enddo
1773  enddo
1774  enddo
1775 
1776 C
1777  102 FORMAT(4(2f9.5,2x))
1778  103 FORMAT(4(3f9.5,2x))
1779  104 FORMAT(4(4f7.3,2x))
1780  105 FORMAT(5f9.5,10x,5f9.5)
1781  106 FORMAT(6f9.5,5x,6f9.5)
1782  107 FORMAT(7f8.4,5x,7f8.4)
1783  108 FORMAT(8f8.4,4x,8f8.4)
1784 c
1785  116 FORMAT( /,5x,' ^ ',/,
1786  $ 5x,' Y | ',/,
1787  $ 5x,' | ',a10,/,
1788  $ 5x,' +----> ','Plane = ',i2,'/',i2,2x,2e12.4,/,
1789  $ 5x,' X ','Step =',i9,f15.5)
1790  117 FORMAT(' ')
1791 c
1792 c if (ichk.eq.1.and.idum.gt.0) call checkit(idum)
1793  return
1794  end
1795 c-----------------------------------------------------------------------
1796  subroutine outflda (x,n,txt10,ichk) ! writes into unit=p130+ifiled
1797  include 'SIZE' ! or into std. output for p130<9
1798  include 'TSTEP' ! truncated below eps=p131
1799  include 'INPUT' ! param(130)
1800  real x(1)
1801  character*10 txt10 ! note: n is not used
1802 c parameter (eps=1.e-7)
1803 C
1804  p130 = param(130)
1805  eps = param(131)
1806  if (p130.le.0) return
1807  m = 6
1808  if (p130.gt.9) m = p130 + ifield
1809 
1810  ntot = lx1*ly1*lz1*nelfld(ifield)
1811 
1812  xmin = glmin(x,ntot)
1813  xmax = glmax(x,ntot)
1814  xavg = glsum(x,ntot)/ntot
1815 
1816  if (abs(xavg).lt.eps) xavg = 0. ! truncation
1817 
1818  if (nid.eq.0) write(m,10) txt10,ichk,ntot,xavg,xmin,xmax
1819 
1820  10 format(3x,a10,2i8,' pts, avg,min,max = ',1p3g14.6)
1821 c
1822  return
1823  end
1824 c-----------------------------------------------------------------------
1825  subroutine outfldan(x,n,txt10,ichk) ! writes x(1:n) into unit=p130+ifiled
1826  include 'SIZE' ! or into std. output for 0<p130<9
1827  include 'TSTEP' ! truncated below eps=p131
1828  include 'INPUT'
1829  real x(1)
1830  character*10 txt10
1831 c parameter (eps=1.e-7)
1832 C
1833  p130 = param(130)
1834  eps = param(131)
1835  if (p130.le.0) return
1836  m = 6
1837  if (p130.gt.9) m = p130 + ifield
1838 
1839  ntot = n
1840 
1841  xmin = glmin(x,ntot)
1842  xmax = glmax(x,ntot)
1843  xavg = glsum(x,ntot)/ntot
1844 
1845  if (abs(xavg).lt.eps) xavg = 0. ! truncation
1846 
1847  if (nid.eq.0) write(m,10) txt10,ichk,ntot,xavg,xmin,xmax
1848 
1849  10 format(3x,a10,2i8,' pts, avg,min,max = ',1p3g11.3)
1850 c 10 format(3X,a10,2i8,' pts, avg,min,max = ',1p3g14.6)
1851 c
1852  return
1853  end
1854 c-----------------------------------------------------------------------
1855  subroutine h1mg_solve(z,rhs,if_hybrid) ! Solve preconditioner: Mz=rhs
1856  real z(1),rhs(1)
1857 
1858 c Assumes that preprocessing has been completed via h1mg_setup()
1859 
1860 
1861  include 'SIZE'
1862  include 'HSMG' ! Same array space as HSMG
1863  include 'GEOM'
1864  include 'INPUT'
1865  include 'MASS'
1866  include 'SOLN'
1867  include 'TSTEP'
1868  include 'CTIMER'
1869  include 'PARALLEL'
1870 
1871  common /scrhi/ h2inv(lx1,ly1,lz1,lelv)
1872  common /scrvh/ h1(lx1,ly1,lz1,lelv),
1873  $ h2(lx1,ly1,lz1,lelv)
1874  parameter(lt=lx1*ly1*lz1*lelt)
1875  common /scrmg/ e(2*lt),w(lt),r(lt)
1876  integer p_msk,p_b
1877  logical if_hybrid
1878 
1879 c if_hybrid = .true. ! Control this from gmres, according
1880 c if_hybrid = .false. ! to convergence efficiency
1881 
1882  nel = nelfld(ifield)
1883 
1884  op = 1. ! Coefficients for h1mg_ax
1885  om = -1.
1886  sigma = 1.
1887  if (if_hybrid) sigma = 2./3.
1888 
1889  l = mg_h1_lmax
1890  n = mg_h1_n(l,mg_fld)
1891  is = 1 ! solve index
1892 
1893  call h1mg_schwarz(z,rhs,sigma,l) ! z := sigma W M rhs
1894  ! Schwarz
1895  call copy(r,rhs,n) ! r := rhs
1896  if (if_hybrid) call h1mg_axm(r,z,op,om,l,w) ! r := rhs - A z
1897  ! l
1898 
1899  do l = mg_h1_lmax-1,2,-1 ! DOWNWARD Leg of V-cycle
1900  is = is + n
1901  n = mg_h1_n(l,mg_fld)
1902  ! T
1903  call h1mg_rstr(r,l,.true.) ! r := J r
1904  ! l l+1
1905 ! OVERLAPPING Schwarz exchange and solve:
1906  call h1mg_schwarz(e(is),r,sigma,l) ! e := sigma W M r
1907  ! l Schwarz l
1908 
1909  if(if_hybrid)call h1mg_axm(r,e(is),op,om,l,w)! r := r - A e
1910  ! l l
1911  enddo
1912  is = is+n
1913  ! T
1914  call h1mg_rstr(r,1,.false.) ! r := J r
1915  ! l l+1
1916  p_msk = p_mg_msk(l,mg_fld)
1917  call h1mg_mask(r,mg_imask(p_msk),nel) ! -1
1918  call hsmg_coarse_solve ( e(is) , r ) ! e := A r
1919  call h1mg_mask(e(is),mg_imask(p_msk),nel) ! 1 1 1
1920 
1921 c nx = mg_nh(1)
1922 c call outnxfld (e(is),nx,nelv,'ecrsb4',is)
1923 c call h1mg_mask(e(is),mg_imask(p_msk),nel) ! 1 1 1
1924 c call outnxfld (e(is),nx,nelv,'ecrsaf',is)
1925 c call exitt
1926 
1927  do l = 2,mg_h1_lmax-1 ! UNWIND. No smoothing.
1928  im = is
1929  is = is - n
1930  n = mg_h1_n(l,mg_fld)
1931  call hsmg_intp (w,e(im),l-1) ! w := J e
1932  i1=is-1 ! l-1
1933  do i=1,n
1934  e(i1+i) = e(i1+i) + w(i) ! e := e + w
1935  enddo ! l l
1936  enddo
1937 
1938  l = mg_h1_lmax
1939  n = mg_h1_n(l,mg_fld)
1940  im = is ! solve index
1941  call hsmg_intp(w,e(im),l-1) ! w := J e
1942  do i = 1,n ! l-1
1943  z(i) = z(i) + w(i) ! z := z + w
1944  enddo
1945 
1946  call dsavg(z) ! Emergency hack --- to ensure continuous z!
1947 
1948  return
1949  end
1950 c-----------------------------------------------------------------------
1951  subroutine h1mg_axm(w,p,aw,ap,l,wk)
1952 c
1953 c w := aw*w + ap*H*p, level l, with mask and dssum
1954 c
1955 c Hu := div. h1 grad u + h2 u
1956 c
1957 c ~= h1 A u + h2 B u
1958 c
1959 c Here, we assume that pointers into g() and h1() and h2() have
1960 c been established
1961 c
1962  include 'SIZE'
1963  include 'HSMG'
1964  include 'TSTEP' ! nelfld
1965 
1966  real w(1),p(1),wk(1)
1967 
1968  integer p_h1,p_h2,p_g,p_b,p_msk
1969  logical ifh2
1970 
1971  p_h1 = p_mg_h1(l,mg_fld)
1972  p_h2 = p_mg_h2(l,mg_fld)
1973  p_g = p_mg_g(l,mg_fld)
1974  p_b = p_mg_b(l,mg_fld)
1975  p_msk = p_mg_msk(l,mg_fld)
1976 
1977  if (p_h1 .eq.0) call mg_set_h1 (p_h1 ,l)
1978  if (p_h2 .eq.0) call mg_set_h2 (p_h2 ,l)
1979  if (p_g .eq.0) call mg_set_gb (p_g,p_b,l)
1980  if (p_msk.eq.0) call mg_set_msk (p_msk,l)
1981 
1982  ifh2 = .true.
1983  if (p_h2.eq.0) ifh2 = .false. ! Need a much better mech.
1984  ifh2 = .false.
1985 
1986  nx = mg_nh(l)
1987  ny = mg_nh(l)
1988  nz = mg_nhz(l)
1989  ng = 3*ldim-3
1990 
1991 
1992  call h1mg_axml (wk,p
1993  $ ,mg_h1(p_h1),mg_h2(p_h2),nx,ny,nz,nelfld(ifield)
1994  $ ,mg_g(p_g) , ng ,mg_b(p_b), mg_imask(p_msk),ifh2)
1995 
1996 
1997  call hsmg_dssum (wk,l)
1998 
1999  n = nx*ny*nz*nelfld(ifield)
2000  call add2sxy (w,aw,wk,ap,n)
2001 
2002  return
2003  end
2004 c-----------------------------------------------------------------------
2005  subroutine h1mg_axml
2006  $ (w,p,h1,h2,nx,ny,nz,nel,g,ng,b,mask,ifh2)
2007 c
2008 c w := aw*w + ap*H*p, level l, with mask and dssum
2009 c
2010 c Hu := div. h1 grad u + h2 u
2011 c
2012 c ~= h1 A u + h2 B u
2013 c
2014 
2015  include 'SIZE'
2016  include 'TOTAL'
2017  include 'HSMG'
2018 
2019  real w (nx*ny*nz,nel), p (nx*ny*nz,nel)
2020  $ , h1(nx*ny*nz,nel), h2(nx*ny*nz,nel)
2021  $ , b (nx*ny*nz,nel), g (ng*nx*ny*nz,nel)
2022  integer mask(1)
2023 
2024  logical ifh2
2025 
2026  parameter(lxyz=lx1*ly1*lz1)
2027  common /ctmp0/ ur(lxyz),us(lxyz),ut(lxyz)
2028 
2029  integer e
2030 
2031  do e=1,nel
2032 
2033  call axe(w(1,e),p(1,e),h1(1,e),h2(1,e),g(1,e),ng,b(1,e)
2034  $ ,nx,ny,nz,ur,us,ut,ifh2,ifrzer(e),e)
2035 
2036  im = mask(e)
2037  call mg_mask_e(w,mask(im)) ! Zero out Dirichlet conditions
2038 
2039  enddo
2040 
2041  return
2042  end
2043 c-----------------------------------------------------------------------
2044  subroutine h1mg_mask(w,mask,nel)
2045  include 'SIZE'
2046 
2047  real w (1)
2048  integer mask(1) ! Pointer to Dirichlet BCs
2049  integer e
2050 
2051  do e=1,nel
2052  im = mask(e)
2053  call mg_mask_e(w,mask(im)) ! Zero out Dirichlet conditions
2054  enddo
2055 
2056  return
2057  end
2058 c----------------------------------------------------------------------
2059  subroutine mg_mask_e(w,mask) ! Zero out Dirichlet conditions
2060  include 'SIZE'
2061  real w(1)
2062  integer mask(0:1)
2063 
2064  n=mask(0)
2065  do i=1,n
2066 c write(6,*) i,mask(i),n,' MG_MASK'
2067  w(mask(i)) = 0.
2068  enddo
2069 
2070  return
2071  end
2072 c-----------------------------------------------------------------------
2073  subroutine axe
2074  $ (w,p,h1,h2,g,ng,b,nx,ny,nz,ur,us,ut,ifh2,ifrz,e)
2075 
2076  include 'SIZE'
2077  include 'INPUT' ! if3d
2078  logical ifh2,ifrz
2079 
2080  real w (nx*ny*nz), p (nx*ny*nz)
2081  $ , h1(nx*ny*nz), h2(nx*ny*nz)
2082  $ , b (nx*ny*nz), g (ng,nx*ny*nz)
2083  $ , ur(nx*ny*nz), us(nx*ny*nz), ut(nx*ny*nz)
2084  integer e
2085 
2086  nxyz = nx*ny*nz
2087 
2088  call gradl_rst(ur,us,ut,p,nx,if3d)
2089 c if (e.eq.1) then
2090 c call outmat(p ,nx,ny,'ur A p',e)
2091 c call outmat(ur,nx,ny,'ur A r',e)
2092 c call outmat(us,nx,ny,'ur A s',e)
2093 c endif
2094 
2095  if (if3d) then
2096  do i=1,nxyz
2097  wr = g(1,i)*ur(i) + g(4,i)*us(i) + g(5,i)*ut(i)
2098  ws = g(4,i)*ur(i) + g(2,i)*us(i) + g(6,i)*ut(i)
2099  wt = g(5,i)*ur(i) + g(6,i)*us(i) + g(3,i)*ut(i)
2100  ur(i) = wr*h1(i)
2101  us(i) = ws*h1(i)
2102  ut(i) = wt*h1(i)
2103  enddo
2104  elseif (ifaxis) then
2105  call exitti('Blame Paul for no gradl_rst support yet$',nx)
2106  else
2107  do i=1,nxyz
2108  wr = g(1,i)*ur(i) + g(3,i)*us(i)
2109  ws = g(3,i)*ur(i) + g(2,i)*us(i)
2110 c write(6,3) i,ur(i),wr,g(1,i)/b(i),b(i)
2111 c 3 format(i5,1p4e12.4,' wrws')
2112  ur(i) = wr*h1(i)
2113  us(i) = ws*h1(i)
2114  enddo
2115  endif
2116 
2117  call gradl_rst_t(w,ur,us,ut,nx,if3d)
2118 
2119 c if (e.eq.1) then
2120 c call outmat(w ,nx,ny,'ur B w',e)
2121 c call outmat(ur,nx,ny,'ur B r',e)
2122 c call outmat(us,nx,ny,'ur B s',e)
2123 c endif
2124 
2125  if (ifh2) then
2126  do i=1,nxyz
2127  w(i)=w(i)+h2(i)*b(i)*p(i)
2128  enddo
2129  endif
2130 
2131  return
2132  end
2133 c----------------------------------------------------------------------
2134  subroutine hsmg_tnsr1(v,nv,nu,A,At)
2135 c
2136 c v = [A (x) A] u or
2137 c v = [A (x) A (x) A] u
2138 c
2139  integer nv,nu
2140  real v(1),A(1),At(1)
2141  include 'SIZE'
2142  include 'INPUT'
2143  if (.not. if3d) then
2144  call hsmg_tnsr1_2d(v,nv,nu,a,at)
2145  else
2146  call hsmg_tnsr1_3d(v,nv,nu,a,at,at)
2147  endif
2148  return
2149  end
2150 c-------------------------------------------------------T--------------
2151  subroutine hsmg_tnsr1_2d(v,nv,nu,A,Bt) ! u = A u B
2152  integer nv,nu
2153  real v(1),A(1),Bt(1)
2154  include 'SIZE'
2155  common /hsmgw/ work(lx1*lx1)
2156  integer e
2157 
2158  nv2 = nv*nv
2159  nu2 = nu*nu
2160 
2161  if (nv.le.nu) then
2162  iv=1
2163  iu=1
2164  do e=1,nelv
2165  call mxm(a,nv,v(iu),nu,work,nu)
2166  call mxm(work,nv,bt,nu,v(iv),nv)
2167  iv = iv + nv2
2168  iu = iu + nu2
2169  enddo
2170  else
2171  do e=nelv,1,-1
2172  iu=1+nu2*(e-1)
2173  iv=1+nv2*(e-1)
2174  call mxm(a,nv,v(iu),nu,work,nu)
2175  call mxm(work,nv,bt,nu,v(iv),nv)
2176  enddo
2177  endif
2178 
2179  return
2180  end
2181 c----------------------------------------------------------------------
2182  subroutine hsmg_tnsr1_3d(v,nv,nu,A,Bt,Ct) ! v = [C (x) B (x) A] u
2183  integer nv,nu
2184  real v(1),A(1),Bt(1),Ct(1)
2185  include 'SIZE'
2186  parameter(lwk=(lx1+2)*(ly1+2)*(lz1+2))
2187  common /hsmgw/ work(0:lwk-1),work2(0:lwk-1)
2188  integer e,e0,ee,es
2189 
2190  e0=1
2191  es=1
2192  ee=nelv
2193 
2194  if (nv.gt.nu) then
2195  e0=nelv
2196  es=-1
2197  ee=1
2198  endif
2199 
2200  nu3 = nu**3
2201  nv3 = nv**3
2202 
2203  do e=e0,ee,es
2204  iu = 1 + (e-1)*nu3
2205  iv = 1 + (e-1)*nv3
2206  call mxm(a,nv,v(iu),nu,work,nu*nu)
2207  do i=0,nu-1
2208  call mxm(work(nv*nu*i),nv,bt,nu,work2(nv*nv*i),nv)
2209  enddo
2210  call mxm(work2,nv*nv,ct,nu,v(iv),nv)
2211  enddo
2212 
2213  return
2214  end
2215 c------------------------------------------ T -----------------------
2216  subroutine h1mg_rstr(r,l,ifdssum) ! r =J r, l is coarse level
2217  include 'SIZE'
2218  include 'HSMG'
2219  logical ifdssum
2220 
2221  real r(1)
2222  integer l
2223 
2224  call hsmg_do_wt(r,mg_rstr_wt(mg_rstr_wt_index(l+1,mg_fld))
2225  $ ,mg_nh(l+1),mg_nh(l+1),mg_nhz(l+1))
2226 
2227  call hsmg_tnsr1(r,mg_nh(l),mg_nh(l+1),mg_jht(1,l),mg_jh(1,l))
2228 
2229  if (ifdssum) call hsmg_dssum(r,l)
2230 
2231  return
2232  end
2233 c----------------------------------------------------------------------
2234  subroutine h1mg_setup()
2235  include 'SIZE'
2236  include 'TOTAL'
2237  include 'HSMG'
2238 
2239  common /scrhi/ h2inv(lx1,ly1,lz1,lelt)
2240  common /scrvh/ h1(lx1,ly1,lz1,lelt),
2241  $ h2(lx1,ly1,lz1,lelt)
2242 
2243  integer p_h1,p_h2,p_g,p_b,p_msk
2244 
2245 
2246  param(59) = 1
2247  call geom_reset(1) ! Recompute g1m1 etc. with deformed only
2248 
2249  n = lx1*ly1*lz1*nelt
2250  call rone (h1 ,n)
2251  call rzero(h2 ,n)
2252  call rzero(h2inv,n)
2253 
2254  call h1mg_setup_mg_nx
2255  call h1mg_setup_semhat ! SEM hat matrices for each level
2256  call hsmg_setup_intp ! Interpolation operators
2257  call h1mg_setup_dssum ! set direct stiffness summation handles
2258  call h1mg_setup_wtmask ! set restriction weight matrices and bc masks
2259  call h1mg_setup_fdm ! set up fast diagonalization method
2260  call h1mg_setup_schwarz_wt(.false.)
2261  call hsmg_setup_solve ! set up the solver
2262 
2263  l=mg_h1_lmax
2264  call mg_set_h1 (p_h1 ,l)
2265  call mg_set_h2 (p_h2 ,l)
2266  call mg_set_gb (p_g,p_b,l)
2267  call mg_set_msk (p_msk,l)
2268 
2269  return
2270  end
2271 c-----------------------------------------------------------------------
2272  subroutine h1mg_setup_mg_nx()
2273  include 'SIZE'
2274  include 'INPUT'
2275  include 'PARALLEL'
2276  include 'HSMG'
2277  include 'SEMHAT'
2278  include 'TSTEP' ! nelfld
2279  real w(lx1+2)
2280  integer nf,nc,nr
2281  integer nx,ny,nz
2282 
2283  integer mgn2(10)
2284  save mgn2
2285  data mgn2 / 1, 2, 2, 2, 2, 3, 3, 5, 5, 5/
2286 c data mgn2 / 1, 2, 3, 4, 5, 6, 7, 8, 9, 0
2287 
2288 c if (param(82).eq.0) param(82)=2 ! nek default
2289 c if (np.eq.1) param(82)=2 ! single proc. too slow
2290  p82 = 2 ! potentially variable nxc
2291  mg_h1_lmax = 3
2292 c mg_h1_lmax = 4
2293  if (lx1.eq.4) mg_h1_lmax = 2
2294 c if (param(79).ne.0) mg_h1_lmax = param(79)
2295  mg_lmax = mg_h1_lmax
2296  mglx1 = p82-1 !1
2297  mg_nx(1) = mglx1
2298  mg_ny(1) = mglx1
2299  mg_nz(1) = mglx1
2300  if (.not.if3d) mg_nz(1) = 0
2301 
2302  mglx2 = 2*(lx2/4) + 1
2303  if (lx1.eq.5) mglx2 = 3
2304 c if (lx1.eq.6) mglx2 = 3
2305  if (lx1.le.10) mglx2 = mgn2(lx1)
2306  if (lx1.eq.8) mglx2 = 4
2307  if (lx1.eq.8) mglx2 = 3
2308 
2309  mglx2 = min(3,mglx2) ! This choice seems best (9/24/12)
2310 
2311  mg_nx(2) = mglx2
2312  mg_ny(2) = mglx2
2313  mg_nz(2) = mglx2
2314  if (.not.if3d) mg_nz(2) = 0
2315 
2316  mg_nx(3) = mglx2+1
2317  mg_ny(3) = mglx2+1
2318  mg_nz(3) = mglx2+1
2319  if (.not.if3d) mg_nz(3) = 0
2320 
2321  mg_nx(mg_h1_lmax) = lx1-1
2322  mg_ny(mg_h1_lmax) = ly1-1
2323  mg_nz(mg_h1_lmax) = lz1-1
2324 
2325  if (nio.eq.0) write(*,*) 'h1_mg_nx:',(mg_nx(i),i=1,mg_h1_lmax)
2326  if (nio.eq.0) write(*,*) 'h1_mg_ny:',(mg_ny(i),i=1,mg_h1_lmax)
2327  if (nio.eq.0) write(*,*) 'h1_mg_nz:',(mg_nz(i),i=1,mg_h1_lmax)
2328 
2329  do ifld=1,ldimt1
2330  do l=1,mg_lmax
2331  mg_h1_n(l,ifld)=(mg_nx(l)+1)
2332  $ *(mg_ny(l)+1)
2333  $ *(mg_nz(l)+1)*nelfld(ifld)
2334  enddo
2335  enddo
2336 
2337  return
2338  end
2339 c----------------------------------------------------------------------
2340  subroutine h1mg_setup_semhat ! SEM hat matrices for each level
2341  include 'SIZE'
2342  include 'INPUT'
2343  include 'HSMG'
2344  include 'SEMHAT'
2345 
2346  do l=1,mg_h1_lmax
2347  n = mg_nx(l) ! polynomial order
2348  call semhat(ah,bh,ch,dh,zh,dph,jph,bgl,zglhat,dgl,jgl,n,wh)
2349  call copy(mg_ah(1,l),ah,(n+1)*(n+1))
2350  call copy(mg_bh(1,l),bh,n+1)
2351  call copy(mg_dh(1,l),dh,(n+1)*(n+1))
2352  call transpose(mg_dht(1,l),n+1,dh,n+1)
2353  call copy(mg_zh(1,l),zh,n+1)
2354 
2355  mg_nh(l)=n+1
2356  mg_nhz(l)=mg_nz(l)+1
2357 
2358  enddo
2359  end
2360 c----------------------------------------------------------------------
2361  subroutine h1mg_setup_dssum
2362  include 'SIZE'
2363  include 'INPUT'
2364  include 'PARALLEL'
2365  include 'HSMG'
2366  parameter(lxyz=(lx1+2)*(ly1+2)*(lz1+2))
2367  common /c_is1/ glo_num(lxyz*lelt)
2368  common /ivrtx/ vertex((2**ldim)*lelt)
2369 
2370  integer*8 glo_num
2371  integer vertex
2372  integer nx,ny,nz
2373  integer l
2374 
2375  ncrnr = 2**ldim
2376  call get_vert
2377 
2378 
2379  do l=1,mg_lmax ! set up direct stiffness summation for each level
2380  nx=mg_nh(l)
2381  ny=mg_nh(l)
2382  nz=mg_nhz(l)
2383  call setupds(mg_gsh_handle(l,mg_fld),nx,ny,nz
2384  $ ,nelv,nelgv,vertex,glo_num)
2385  nx=nx+2
2386  ny=ny+2
2387  nz=nz+2
2388  if(.not.if3d) nz=1
2389  call setupds(mg_gsh_schwarz_handle(l,mg_fld),nx,ny,nz
2390  $ ,nelv,nelgv,vertex,glo_num)
2391  enddo
2392 
2393  return
2394  end
2395 c----------------------------------------------------------------------
2396  subroutine mg_set_msk(p_msk ,l0)
2397  include 'SIZE'
2398  include 'HSMG'
2399  include 'TSTEP'
2400  integer p_msk
2401 
2402  l = mg_h1_lmax
2403  p_mg_msk(l,mg_fld) = 0
2404  n = mg_h1_n(l,mg_fld)
2405 
2406 
2407  do l=mg_h1_lmax,1,-1
2408  nx = mg_nh(l)
2409  ny = mg_nh(l)
2410  nz = mg_nhz(l)
2411 
2412  p_msk = p_mg_msk(l,mg_fld)
2413 
2414  call h1mg_setup_mask
2415  $ (mg_imask(p_msk),nm,nx,ny,nz,nelfld(ifield),l,mg_work)
2416 
2417  if (l.gt.1) p_mg_msk(l-1,mg_fld)=p_mg_msk(l,mg_fld)+nm
2418 
2419  enddo
2420 
2421  p_msk = p_mg_msk(l0,mg_fld)
2422 
2423  return
2424  end
2425 c----------------------------------------------------------------------
2426  subroutine h1mg_setup_mask(mask,nm,nx,ny,nz,nel,l,w)
2427  include 'SIZE'
2428  include 'INPUT' ! if3d
2429 
2430  integer mask(1) ! Pointer to Dirichlet BCs
2431  integer nx,ny,nz,l
2432  real w(nx,ny,nz,nel)
2433 
2434  integer e,count,ptr
2435  integer lbr,rbr,lbs,rbs,lbt,rbt,two
2436 
2437  zero = 0
2438  nxyz = nx*ny*nz
2439  n = nx*ny*nz*nel
2440 
2441  call rone(w,n) ! Init everything to 1
2442 
2443  ierrmx = 0 ! BC verification
2444  two = 2
2445  do e=1,nel ! Set dirichlet nodes to zero
2446 
2447  call get_fast_bc(lbr,rbr,lbs,rbs,lbt,rbt,e,two,ierr)
2448 c write(6,6) e,lbr,rbr,lbs,rbs,ierr,nx
2449 c 6 format(i5,2x,4i3,2x,i2,3x,i5,' lbr,rbr,lbs')
2450 
2451  if (lbr.eq.1) call facev(w,e,4,zero,nx,ny,nz)
2452  if (rbr.eq.1) call facev(w,e,2,zero,nx,ny,nz)
2453  if (lbs.eq.1) call facev(w,e,1,zero,nx,ny,nz)
2454  if (rbs.eq.1) call facev(w,e,3,zero,nx,ny,nz)
2455  if (if3d) then
2456  if (lbt.eq.1) call facev(w,e,5,zero,nx,ny,nz)
2457  if (rbt.eq.1) call facev(w,e,6,zero,nx,ny,nz)
2458  endif
2459  ierrmx = max(ierrmx,ierr)
2460  enddo
2461 
2462  call hsmg_dsprod(w,l) ! direct stiffness multiply
2463 
2464 c
2465 c Prototypical mask layout, nel=5:
2466 c
2467 c e=1 ... 10
2468 c 1 2 3 4 5 ... 10 | 11 12 13 14 | 15 | 16 |
2469 c 11 15 16 ... | 3 p1 p2 p3 | 0 | 0 | ...
2470 c ^
2471 c |
2472 c |_count for e=1
2473 c
2474 
2475  nm = 1 ! store mask
2476  do e=1,nel
2477 
2478  mask(e) = nel+nm
2479  count = 0 ! # Dirchlet points on element e
2480  ptr = mask(e)
2481 
2482  do i=1,nxyz
2483  if (w(i,1,1,e).eq.0) then
2484  nm = nm +1
2485  count = count+1
2486  ptr = ptr +1
2487  mask(ptr) = i + nxyz*(e-1) ! where I mask on element e
2488  endif
2489  enddo
2490 
2491 
2492  ptr = mask(e)
2493  mask(ptr) = count
2494 
2495  nm = nm+1 ! bump pointer to hold next count
2496 
2497  enddo
2498 
2499  nm = nel + nm-1 ! Return total number of mask pointers/counters
2500 
2501  ierrmx = iglmax(ierrmx,1)
2502  if (ierrmx.gt.0) then
2503  if (ierr.gt.0) write(6,*) nid,ierr,' BC FAIL h1'
2504  call exitti('D INVALID BC FOUND in h1mg_setup_mask$',ierrmx)
2505  endif
2506 
2507  return
2508  end
2509 c----------------------------------------------------------------------
2510  subroutine mg_set_h1 (p_h1 ,l0)
2511  include 'SIZE'
2512  include 'HSMG'
2513  integer pf,pc
2514 
2515 c As a first pass, rely on the cheesy common-block interface to get h1
2516 
2517  common /scrvh/ h1(lx1,ly1,lz1,lelv)
2518  $ , h2(lx1,ly1,lz1,lelv)
2519  $ , h2inv(lx1,ly1,lz1,lelv)
2520 
2521  integer p_h1
2522 
2523  l = mg_h1_lmax
2524  p_mg_h1(l,mg_fld) = 0
2525  n = mg_h1_n(l,mg_fld)
2526 
2527  call copy (mg_h1,h1,n) ! Fine grid is just original h1
2528 
2529  nx = mg_nh(l)
2530  ny = mg_nh(l)
2531  nz = mg_nhz(l)
2532 
2533  do l=mg_h1_lmax-1,1,-1
2534 
2535  p_mg_h1(l,mg_fld) = p_mg_h1(l+1,mg_fld) + n
2536  n = mg_h1_n(l ,mg_fld)
2537 
2538  pf = p_mg_h1(l+1,mg_fld)
2539  pc = p_mg_h1(l ,mg_fld)
2540 
2541  call hsmg_intp_fc (mg_h1(pc),mg_h1(pf),l)
2542 
2543  enddo
2544 
2545  p_h1 = p_mg_h1(l0,mg_fld)
2546 
2547  return
2548  end
2549 c-----------------------------------------------------------------------
2550  subroutine mg_set_h2 (p_h2 ,l0)
2551  include 'SIZE'
2552  include 'HSMG'
2553 
2554 c As a first pass, rely on the cheesy common-block interface to get h2
2555 
2556  common /scrvh/ h1(lx1,ly1,lz1,lelv)
2557  $ , h2(lx1,ly1,lz1,lelv)
2558  $ , h2inv(lx1,ly1,lz1,lelv)
2559 
2560  integer p_h2,pf,pc
2561 
2562  l = mg_h1_lmax
2563  p_mg_h2(l,mg_fld) = 0
2564  n = mg_h1_n(l,mg_fld)
2565 
2566  call copy (mg_h2,h2,n) ! Fine grid is just original h2
2567 
2568  nx = mg_nh(l)
2569  ny = mg_nh(l)
2570  nz = mg_nhz(l)
2571 
2572  do l=mg_h1_lmax-1,1,-1
2573 
2574  p_mg_h2(l,mg_fld) = p_mg_h2(l+1,mg_fld) + n
2575  n = mg_h1_n(l ,mg_fld)
2576 
2577  pf = p_mg_h2(l+1,mg_fld)
2578  pc = p_mg_h2(l ,mg_fld)
2579 
2580  call hsmg_intp_fc (mg_h2(pc),mg_h2(pf),l)
2581 
2582  enddo
2583 
2584  p_h2 = p_mg_h2(l0,mg_fld)
2585 
2586  return
2587  end
2588 c-----------------------------------------------------------------------
2589  subroutine hsmg_intp_fc(uc,uf,l) ! l is coarse level
2590 
2591  include 'SIZE'
2592  include 'HSMG'
2593 
2594  real uc(1),uf(1)
2595 
2596 
2597  nc = mg_nh(l)
2598  nf = mg_nh(l+1)
2599  call hsmg_tnsr(uc,nc,uf,nf,mg_jhfc(1,l),mg_jhfct(1,l))
2600 
2601  return
2602  end
2603 c-----------------------------------------------------------------------
2604  subroutine mg_intp_fc_e(uc,uf,nxc,nyc,nzc,nxf,nyf,nzf,e,l,w)
2605  include 'SIZE'
2606  include 'INPUT' ! if3d
2607  include 'HSMG'
2608 
2609  real uf(nxf,nyf,nzf),uc(nxc,nyc,nzc),w(1)
2610  integer e
2611 
2612  if (if3d) then
2613 
2614  n1=nxf*nyf
2615  n2=nzf
2616  n3=nzc
2617  call mxm(uf,n1,mg_jhfct(1,l),n2,w,n3)
2618 
2619  lf=1 ! pointers into work array w()
2620  lc=1 + n1*n3
2621  lc0=lc
2622 
2623  n1=nxf
2624  n2=nyf
2625  n3=nyc
2626 
2627  do k=1,nzc
2628  call mxm(w(lf),n1,mg_jhfct(1,l),n2,w(lc),n3)
2629  lf = lf + n1*n2
2630  lc = lc + n1*n3
2631  enddo
2632 
2633  lf=lc0 ! Rewind fine pointer to start of coarse data
2634  n1=nxc
2635  n2=nxf
2636  n3=nyc*nzc
2637  call mxm(mg_jhfc(1,l),n1,w(lf),n2,uc,n3)
2638 
2639  else ! 2D
2640 
2641  n1=nxf
2642  n2=nyf
2643  n3=nyc
2644  call mxm(uf,n1,mg_jhfct(1,l),n2,w,n3)
2645 
2646  n1=nxc
2647  n2=nxf
2648  n3=nyc
2649  call mxm(mg_jhfc(1,l),n1,w,n2,uc,n3)
2650 
2651  endif
2652 
2653  return
2654  end
2655 c-----------------------------------------------------------------------
2656  subroutine mg_intp_gfc_e(gc,gf,ng,nxc,nyc,nzc,nxf,nyf,nzf,e,l,w)
2657  include 'SIZE'
2658  include 'INPUT' ! if3d
2659  include 'HSMG'
2660 
2661  real gf(ng,nxf,nyf,nzf),gc(ng,nxc,nyc,nzc),w(1)
2662  integer e
2663 
2664 
2665  if (if3d) then
2666 
2667  n1=ng*nxf*nyf
2668  n2=nzf
2669  n3=nzc
2670  call mxm(gf,n1,mg_jhfct(1,l),n2,w,n3)
2671 
2672  lf=1 ! pointers into work array w()
2673  lc=1 + n1*n3
2674  lc0=lc
2675 
2676  n1=ng*nxf
2677  n2=nyf
2678  n3=nyc
2679 
2680  do k=1,nzc
2681  call mxm(w(lf),n1,mg_jhfct(1,l),n2,w(lc),n3)
2682  lf = lf + n1*n2
2683  lc = lc + n1*n3
2684  enddo
2685 
2686  lf=lc0 ! Rewind fine pointer to start of coarse data
2687  n1=ng
2688  n2=nxf
2689  n3=nxc
2690 
2691  do k=1,nyc*nzc
2692  call mxm(w(lf),n1,mg_jhfct(1,l),n2,gc(1,1,k,1),n3)
2693  lf = lf + n1*n2
2694  enddo
2695 
2696  else ! 2D
2697 
2698  n1=ng*nxf
2699  n2=nyf
2700  n3=nyc
2701  call mxm(gf,n1,mg_jhfct(1,l),n2,w,n3)
2702 
2703  lf=1 ! pointers into work array w()
2704 
2705  n1=ng
2706  n2=nxf
2707  n3=nxc
2708 
2709  do k=1,nyc
2710  call mxm(w(lf),n1,mg_jhfct(1,l),n2,gc(1,1,k,1),n3)
2711  lf = lf + n1*n2
2712  enddo
2713 
2714  endif
2715 
2716  return
2717  end
2718 c-----------------------------------------------------------------------
2719  subroutine mg_scale_mass (b,g,wt,ng,nx,ny,nz,wk,ifinv)
2720  include 'SIZE'
2721  include 'INPUT' ! if3d
2722  include 'HSMG'
2723 
2724  real b(1),g(ng,1),wt(1),wk(1)
2725  logical ifinv
2726 
2727  common /ctmp0/ wi(2*lx1+4)
2728 
2729  n = nx*ny*nz
2730 
2731  if (nx.le.2*lx1) then
2732  if (ifinv) then
2733  call invers2(wi,wt,nx)
2734  else
2735  call copy(wi,wt,nx)
2736  endif
2737  else
2738  call exitti('mg_scale_mass: wi too small$',nx)
2739  endif
2740 
2741  if (if3d) then
2742  l=0
2743  do k=1,nz
2744  do j=1,ny
2745  wjk=wi(j)*wi(k)
2746  do i=1,nx
2747  l=l+1
2748  wk(l) = wjk*wi(i)
2749  enddo
2750  enddo
2751  enddo
2752 
2753  do k=1,n
2754  b(k) = wk(k)*b(k)
2755  g(1,k) = wk(k)*g(1,k)
2756  g(2,k) = wk(k)*g(2,k)
2757  g(3,k) = wk(k)*g(3,k)
2758  g(4,k) = wk(k)*g(4,k)
2759  g(5,k) = wk(k)*g(5,k)
2760  g(6,k) = wk(k)*g(6,k)
2761  enddo
2762 
2763  else ! 2D
2764  l=0
2765  do j=1,ny
2766  do i=1,nx
2767  l=l+1
2768  wk(l) = wi(i)*wi(j)
2769  enddo
2770  enddo
2771 
2772  do k=1,n
2773  b(k) = wk(k)*b(k)
2774  g(1,k) = wk(k)*g(1,k)
2775  g(2,k) = wk(k)*g(2,k)
2776  g(3,k) = wk(k)*g(3,k)
2777  enddo
2778 
2779  endif
2780 
2781  return
2782  end
2783 c-----------------------------------------------------------------------
2784  subroutine mg_set_gb (p_g,p_b,l0)
2785  include 'SIZE'
2786  include 'HSMG'
2787  include 'MASS' ! bm1
2788  include 'TSTEP' ! nelfld
2789 
2790  integer p_g,p_b,e
2791  common /ctmp1/ w(lx1*ly1*lz1*lelt*2)
2792 
2793  l = mg_h1_lmax
2794  p_mg_b(l,mg_fld) = 0
2795  p_mg_g(l,mg_fld) = 0
2796  n = mg_h1_n(l,mg_fld)
2797 
2798 
2799  ng = 3*(ldim-1) ! 3 or 6 elements to symm dxd tensor
2800 
2801  do l=mg_h1_lmax-1,1,-1
2802 
2803  p_mg_b(l,mg_fld) = p_mg_b(l+1,mg_fld) + n
2804  p_mg_g(l,mg_fld) = p_mg_g(l+1,mg_fld) + n*ng
2805  n = mg_h1_n(l ,mg_fld)
2806 
2807  enddo
2808 
2809  do e=1,nelfld(ifield)
2810  do l=mg_h1_lmax,1,-1
2811 
2812  nx = mg_nh(l)
2813  ny = mg_nh(l)
2814  nz = mg_nhz(l)
2815  nxyz = nx*ny*nz
2816 
2817  p_g = p_mg_g(l,mg_fld) + ng*nx*ny*nz*(e-1)
2818  p_b = p_mg_b(l,mg_fld) + nx*ny*nz*(e-1)
2819 
2820  if (l.eq.mg_h1_lmax) then
2821  call gxfer_e (mg_g(p_g) ,ng,e ) ! Fine grid=original G
2822  call copy (mg_b(p_b) ,bm1(1,1,1,e),nxyz) ! Fine grid=original B
2823  call mg_scale_mass ! Divide out Wghts
2824  $ (mg_b(p_b),mg_g(p_g),mg_bh(1,l),ng,nx,ny,nz,w,.true.)
2825  else
2826 
2827 c Generate G and B by interpolating their continous counterparts onto
2828 c the coarse grid and collocating with coarse-grid quadrature weights
2829 
2830  call mg_intp_gfc_e
2831  $ (mg_g(p_g),mg_g(l_g),ng,nx,ny,nz,nxl,nyl,nzl,e,l,w)
2832 
2833  call mg_intp_fc_e
2834  $ (mg_b(p_b),mg_b(l_b) ,nx,ny,nz,nxl,nyl,nzl,e,l,w)
2835 
2836  call mg_scale_mass ! Reinstate weights
2837  $ (mg_b(l_b),mg_g(l_g),mg_bh(1,l+1),ng,nxl,nyl,nzl,w,.false.)
2838 
2839  endif
2840 
2841  l_b = p_b
2842  l_g = p_g
2843 
2844  nxl = nx
2845  nyl = ny
2846  nzl = nz
2847 
2848  enddo
2849 
2850  call mg_scale_mass ! Reinstate weights
2851  $ (mg_b(l_b),mg_g(l_g),mg_bh(1,1),ng,nxl,nyl,nzl,w,.false.)
2852 
2853 
2854  enddo
2855 
2856  p_b = p_mg_b(l0,mg_fld)
2857  p_g = p_mg_g(l0,mg_fld)
2858 
2859  return
2860  end
2861 c-----------------------------------------------------------------------
2862  subroutine gxfer_e (g,ng,e)
2863  include 'SIZE'
2864  include 'TOTAL'
2865 
2866  real g(ng,1)
2867  integer e
2868 
2869  nxyz = lx1*ly1*lz1
2870 
2871 c ifdfrm(e) = .true. ! TOO LATE
2872 
2873  if (if3d) then
2874  do i=1,nxyz
2875  g(1,i) = g1m1(i,1,1,e)
2876  g(2,i) = g2m1(i,1,1,e)
2877  g(3,i) = g3m1(i,1,1,e)
2878  g(4,i) = g4m1(i,1,1,e)
2879  g(5,i) = g5m1(i,1,1,e)
2880  g(6,i) = g6m1(i,1,1,e)
2881  enddo
2882  else
2883  do i=1,nxyz
2884  g(1,i) = g1m1(i,1,1,e)
2885  g(2,i) = g2m1(i,1,1,e)
2886  g(3,i) = g4m1(i,1,1,e)
2887  enddo
2888  endif
2889 
2890  return
2891  end
2892 c-----------------------------------------------------------------------
2893  subroutine chkr(name3,ii)
2894  include 'SIZE'
2895  include 'TOTAL'
2896  include 'HSMG'
2897  character*3 name3
2898 
2899  write(6,*) mg_h1_lmax,ii,' ',name3,' CHKR'
2900 
2901  return
2902  end
2903 c-----------------------------------------------------------------------
2904  subroutine outgmat(a,ng,nx,ny,name6,k,e)
2905 
2906  integer e
2907  real a(ng,nx,ny)
2908  common /ctmp0/ w(100000)
2909  character*6 name6
2910 
2911 c do i=1,ng
2912  do i=1,1
2913  sum = 0.
2914  do ii=1,nx*ny
2915  w(ii)=a(i,ii,1)
2916  sum = sum + a(i,ii,1)
2917  enddo
2918 
2919  write(6,1) name6,i,k,e,nx,ny,ng,sum
2920  1 format(a6,6i5,f12.5,' outgmat')
2921 
2922  call outmatz(w,nx,ny,name6,i)
2923  enddo
2924 
2925  return
2926  end
2927 c-----------------------------------------------------------------------
2928  subroutine outmatz(a,m,n,name6,ie)
2929  real a(m,n)
2930  character*6 name6
2931 
2932  sum=0.
2933  sua=0.
2934  do i=1,m*n
2935  sum=sum+ a(i,1)
2936  sua=sua+abs(a(i,1))
2937  enddo
2938  sum=sum/(m*n)
2939  sua=sua/(m*n)
2940 
2941  write(6,*)
2942  write(6,1) ie,name6,m,n,sum,sua
2943  1 format(i8,' matrix: ',a6,2i5,1p2e12.4)
2944 
2945  n12 = min(m,12)
2946  do j=m,1,-1
2947  write(6,6) ie,name6,(a(i,j),i=1,n12)
2948  enddo
2949  6 format(i3,1x,a6,12f9.5)
2950 c write(6,*)
2951  return
2952  end
2953 c-----------------------------------------------------------------------
2954  subroutine h1mg_setup_schwarz_wt_2(wt,ie,n,work,ifsqrt)
2955  include 'SIZE'
2956  real wt(1),work(1)
2957  logical ifsqrt
2958 
2959  if (ldim.eq.2) call h1mg_setup_schwarz_wt2d_2(wt,ie,n,work,ifsqrt)
2960  if (ldim.eq.3) call h1mg_setup_schwarz_wt3d_2(wt,ie,n,work,ifsqrt)
2961 
2962  return
2963  end
2964 c----------------------------------------------------------------------
2965  subroutine h1mg_setup_schwarz_wt2d_2(wt,ie,n,work,ifsqrt)
2966  include 'SIZE'
2967  logical ifsqrt
2968  integer n
2969  real wt(n,4,2,nelv)
2970  real work(n,n)
2971 
2972  integer ie,i,j
2973  do j=1,n
2974  wt(j,1,1,ie)=1.0/work(1,j)
2975  wt(j,2,1,ie)=1.0/work(2,j)
2976  wt(j,3,1,ie)=1.0/work(n-1,j)
2977  wt(j,4,1,ie)=1.0/work(n,j)
2978  enddo
2979  do i=1,n
2980  wt(i,1,2,ie)=1.0/work(i,1)
2981  wt(i,2,2,ie)=1.0/work(i,2)
2982  wt(i,3,2,ie)=1.0/work(i,n-1)
2983  wt(i,4,2,ie)=1.0/work(i,n)
2984  enddo
2985  if(ifsqrt) then
2986  do ii=1,2
2987  do j=1,4
2988  do i=1,n
2989  wt(i,j,ii,ie)=sqrt(wt(i,j,ii,ie))
2990  enddo
2991  enddo
2992  enddo
2993  endif
2994 
2995  return
2996  end
2997 c----------------------------------------------------------------------
2998  subroutine h1mg_setup_schwarz_wt3d_2(wt,ie,n,work,ifsqrt)
2999  include 'SIZE'
3000  logical ifsqrt
3001  integer n
3002  real wt(n,n,4,3,nelv)
3003  real work(n,n,n)
3004 
3005  integer ie,i,j,k
3006  integer lbr,rbr,lbs,rbs,lbt,rbt
3007 
3008  ierr = 0
3009  do k=1,n
3010  do j=1,n
3011  wt(j,k,1,1,ie)=1.0/work(1,j,k)
3012  wt(j,k,2,1,ie)=1.0/work(2,j,k)
3013  wt(j,k,3,1,ie)=1.0/work(n-1,j,k)
3014  wt(j,k,4,1,ie)=1.0/work(n,j,k)
3015  enddo
3016  enddo
3017  do k=1,n
3018  do i=1,n
3019  wt(i,k,1,2,ie)=1.0/work(i,1,k)
3020  wt(i,k,2,2,ie)=1.0/work(i,2,k)
3021  wt(i,k,3,2,ie)=1.0/work(i,n-1,k)
3022  wt(i,k,4,2,ie)=1.0/work(i,n,k)
3023  enddo
3024  enddo
3025  do j=1,n
3026  do i=1,n
3027  wt(i,j,1,3,ie)=1.0/work(i,j,1)
3028  wt(i,j,2,3,ie)=1.0/work(i,j,2)
3029  wt(i,j,3,3,ie)=1.0/work(i,j,n-1)
3030  wt(i,j,4,3,ie)=1.0/work(i,j,n)
3031  enddo
3032  enddo
3033  if(ifsqrt) then
3034  do ii=1,3
3035  do k=1,4
3036  do j=1,4
3037  do i=1,n
3038  wt(i,j,k,ii,ie)=sqrt(wt(i,j,k,ii,ie))
3039  enddo
3040  enddo
3041  enddo
3042  enddo
3043  endif
3044 
3045  return
3046  end
3047 c----------------------------------------------------------------------
3048  subroutine h1mg_setup_schwarz_wt_1(wt,l,ifsqrt)
3049  include 'SIZE'
3050  include 'INPUT' ! if3d
3051  include 'TSTEP' ! ifield
3052  include 'HSMG'
3053 
3054  real wt(1),work(1)
3055  logical ifsqrt
3056 
3057  integer enx,eny,enz,pm
3058 
3059  zero = 0
3060  one = 1
3061  onem = -1
3062 
3063  n = mg_h1_n(l,mg_fld)
3064  pm = p_mg_msk(l,mg_fld)
3065 
3066  enx=mg_nh(l)+2
3067  eny=mg_nh(l)+2
3068  enz=mg_nh(l)+2
3069  if(.not.if3d) enz=1
3070  ns = enx*eny*enz*nelfld(ifield)
3071  i = ns+1
3072 
3073  call rone(mg_work(i),ns)
3074 
3075 c Sum overlap region (border excluded)
3076  call hsmg_extrude(mg_work,0,zero,mg_work(i),0,one ,enx,eny,enz)
3077  call hsmg_schwarz_dssum(mg_work(i),l)
3078  call hsmg_extrude(mg_work(i),0,one ,mg_work,0,onem,enx,eny,enz)
3079  call hsmg_extrude(mg_work(i),2,one,mg_work(i),0,one,enx,eny,enz)
3080 
3081  if(.not.if3d) then ! Go back to regular size array
3082  call hsmg_schwarz_toreg2d(mg_work,mg_work(i),mg_nh(l))
3083  else
3084  call hsmg_schwarz_toreg3d(mg_work,mg_work(i),mg_nh(l))
3085  endif
3086 
3087  call hsmg_dssum(mg_work,l) ! sum border nodes
3088 
3089 
3090  nx = mg_nh(l)
3091  ny = mg_nh(l)
3092  nz = mg_nh(l)
3093  if (.not.if3d) nz=1
3094  nxyz = nx*ny*nz
3095  k = 1
3096  do ie=1,nelfld(ifield)
3097 c call outmat(mg_work(k),nx,ny,'NEW WT',ie)
3098  call h1mg_setup_schwarz_wt_2(wt,ie,nx,mg_work(k),ifsqrt)
3099  k = k+nxyz
3100  enddo
3101 c stop
3102 
3103  return
3104  end
3105 c----------------------------------------------------------------------
exitt
void exitt()
Definition: comm_mpi.f:604
gop
subroutine gop(x, w, op, n)
Definition: comm_mpi.f:201
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
facev
subroutine facev(a, ie, iface, val, nx, ny, nz)
Definition: connect1.f:1077
setupds
subroutine setupds(gs_handle, nx, ny, nz, nel, melg, vertex, glo_num)
Definition: dssum.f:2
row_zero
subroutine row_zero(a, m, n, e)
Definition: fast3d.f:1617
semhat
subroutine semhat(a, b, c, d, z, dgll, jgll, bgl, zgl, dgl, jgl, n, w)
Definition: fast3d.f:1216
get_fast_bc
subroutine get_fast_bc(lbr, rbr, lbs, rbs, lbt, rbt, e, bsym, ierr)
Definition: fast3d.f:803
fd_weights_full
subroutine fd_weights_full(xx, x, n, m, c)
Definition: fast3d.f:1293
local_solves_fdm
subroutine local_solves_fdm(u, v)
Definition: fasts.f:3
generalev
subroutine generalev(a, b, lam, n, w)
Definition: hmholtz.f:1369
outfldan
subroutine outfldan(x, n, txt10, ichk)
Definition: hsmg.f:1826
h1mg_setup_schwarz_wt3d_2
subroutine h1mg_setup_schwarz_wt3d_2(wt, ie, n, work, ifsqrt)
Definition: hsmg.f:2999
hsmg_tnsr2d
subroutine hsmg_tnsr2d(v, nv, u, nu, A, Bt)
Definition: hsmg.f:259
hsmg_rstr_no_dssum
subroutine hsmg_rstr_no_dssum(uc, uf, l)
Definition: hsmg.f:228
h1mg_setup_fdm
subroutine h1mg_setup_fdm()
Definition: hsmg.f:634
mg_set_gb
subroutine mg_set_gb(p_g, p_b, l0)
Definition: hsmg.f:2785
h1mg_setup_wtmask
subroutine h1mg_setup_wtmask
Definition: hsmg.f:158
hsmg_extrude
subroutine hsmg_extrude(arr1, l1, f1, arr2, l2, f2, nx, ny, nz)
Definition: hsmg.f:369
hsmg_rstr
subroutine hsmg_rstr(uc, uf, l)
Definition: hsmg.f:215
hsmg_schwarz_toreg3d
subroutine hsmg_schwarz_toreg3d(b, a, n)
Definition: hsmg.f:616
mg_set_msk
subroutine mg_set_msk(p_msk, l0)
Definition: hsmg.f:2397
outmatz
subroutine outmatz(a, m, n, name6, ie)
Definition: hsmg.f:2929
hsmg_tnsr
subroutine hsmg_tnsr(v, nv, u, nu, A, At)
Definition: hsmg.f:243
h1mg_setup_mg_nx
subroutine h1mg_setup_mg_nx()
Definition: hsmg.f:2273
hsmg_tnsr3d_el
subroutine hsmg_tnsr3d_el(v, nv, u, nu, A, Bt, Ct)
Definition: hsmg.f:310
hsmg_intp
subroutine hsmg_intp(uf, uc, l)
Definition: hsmg.f:206
h1mg_setup_schwarz_wt
subroutine h1mg_setup_schwarz_wt(ifsqrt)
Definition: hsmg.f:1218
mg_intp_gfc_e
subroutine mg_intp_gfc_e(gc, gf, ng, nxc, nyc, nzc, nxf, nyf, nzf, e, l, w)
Definition: hsmg.f:2657
h1mg_setup_schwarz_wt_1
subroutine h1mg_setup_schwarz_wt_1(wt, l, ifsqrt)
Definition: hsmg.f:3049
hsmg_tnsr3d
subroutine hsmg_tnsr3d(v, nv, u, nu, A, Bt, Ct)
Definition: hsmg.f:275
h1mg_schwarz_part1
subroutine h1mg_schwarz_part1(e, r, l)
Definition: hsmg.f:441
hsmg_schwarz_wt3d
subroutine hsmg_schwarz_wt3d(e, wt, n)
Definition: hsmg.f:1286
hsmg_setup_rstr_wt
subroutine hsmg_setup_rstr_wt(wt, nx, ny, nz, l, w)
Definition: hsmg.f:988
gxfer_e
subroutine gxfer_e(g, ng, e)
Definition: hsmg.f:2863
mg_mask_e
subroutine mg_mask_e(w, mask)
Definition: hsmg.f:2060
mg_scale_mass
subroutine mg_scale_mass(b, g, wt, ng, nx, ny, nz, wk, ifinv)
Definition: hsmg.f:2720
hsmg_tnsr1_2d
subroutine hsmg_tnsr1_2d(v, nv, nu, A, Bt)
Definition: hsmg.f:2152
hsmg_fdm
subroutine hsmg_fdm(e, r, l)
Definition: hsmg.f:886
hsmg_dssum
subroutine hsmg_dssum(u, l)
Definition: hsmg.f:327
hsmg_schwarz_toreg2d
subroutine hsmg_schwarz_toreg2d(b, a, n)
Definition: hsmg.f:600
hsmg_setup_fast1d
subroutine hsmg_setup_fast1d(s, lam, nl, lbc, rbc, ll, lm, lr, ah, bh, n, ie)
Definition: hsmg.f:784
hsmg_setup
subroutine hsmg_setup()
Definition: hsmg.f:23
hsmg_tnsr1_3d
subroutine hsmg_tnsr1_3d(v, nv, nu, A, Bt, Ct)
Definition: hsmg.f:2183
hsmg_schwarz_dssum
subroutine hsmg_schwarz_dssum(u, l)
Definition: hsmg.f:355
h1mg_rstr
subroutine h1mg_rstr(r, l, ifdssum)
Definition: hsmg.f:2217
outgmat
subroutine outgmat(a, ng, nx, ny, name6, k, e)
Definition: hsmg.f:2905
hsmg_schwarz
subroutine hsmg_schwarz(e, r, l)
Definition: hsmg.f:497
h1mg_axm
subroutine h1mg_axm(w, p, aw, ap, l, wk)
Definition: hsmg.f:1952
h1mg_setup
subroutine h1mg_setup()
Definition: hsmg.f:2235
mg_set_h1
subroutine mg_set_h1(p_h1, l0)
Definition: hsmg.f:2511
hsmg_setup_fast1d_a
subroutine hsmg_setup_fast1d_a(a, lbc, rbc, ll, lm, lr, ah, n)
Definition: hsmg.f:807
chkr
subroutine chkr(name3, ii)
Definition: hsmg.f:2894
mg_intp_fc_e
subroutine mg_intp_fc_e(uc, uf, nxc, nyc, nzc, nxf, nyf, nzf, e, l, w)
Definition: hsmg.f:2605
outfldn0
subroutine outfldn0(x, n, txt10, ichk)
Definition: hsmg.f:1739
hsmg_setup_solve
subroutine hsmg_setup_solve
Definition: hsmg.f:1357
outfldn
subroutine outfldn(x, n, txt10, ichk)
Definition: hsmg.f:1680
h1mg_solve
subroutine h1mg_solve(z, rhs, if_hybrid)
Definition: hsmg.f:1856
hsmg_setup_mask
subroutine hsmg_setup_mask(wt, nx, ny, nz, l, w)
Definition: hsmg.f:1070
h1mg_setup_mask
subroutine h1mg_setup_mask(mask, nm, nx, ny, nz, nel, l, w)
Definition: hsmg.f:2427
axe
subroutine axe(w, p, h1, h2, g, ng, b, nx, ny, nz, ur, us, ut, ifh2, ifrz, e)
Definition: hsmg.f:2075
hsmg_intp_fc
subroutine hsmg_intp_fc(uc, uf, l)
Definition: hsmg.f:2590
h1mg_mask
subroutine h1mg_mask(w, mask, nel)
Definition: hsmg.f:2045
hsmg_setup_fdm
subroutine hsmg_setup_fdm()
Definition: hsmg.f:668
hsmg_schwarz_wt
subroutine hsmg_schwarz_wt(e, l)
Definition: hsmg.f:1250
hsmg_tnsr2d_el
subroutine hsmg_tnsr2d_el(v, nv, u, nu, A, Bt)
Definition: hsmg.f:295
outflda
subroutine outflda(x, n, txt10, ichk)
Definition: hsmg.f:1797
hsmg_do_wt
subroutine hsmg_do_wt(u, wt, nx, ny, nz)
Definition: hsmg.f:933
mg_set_h2
subroutine mg_set_h2(p_h2, l0)
Definition: hsmg.f:2551
hsmg_setup_semhat
subroutine hsmg_setup_semhat
Definition: hsmg.f:51
h1mg_setup_schwarz_wt2d_2
subroutine h1mg_setup_schwarz_wt2d_2(wt, ie, n, work, ifsqrt)
Definition: hsmg.f:2966
h1mg_axml
subroutine h1mg_axml(w, p, h1, h2, nx, ny, nz, nel, g, ng, b, mask, ifh2)
Definition: hsmg.f:2007
hsmg_setup_intp
subroutine hsmg_setup_intp
Definition: hsmg.f:84
h1mg_schwarz
subroutine h1mg_schwarz(e, r, sigma, l)
Definition: hsmg.f:426
hsmg_solve
subroutine hsmg_solve(e, r)
Definition: hsmg.f:1377
hsmg_setup_wtmask
subroutine hsmg_setup_wtmask
Definition: hsmg.f:182
hsmg_coarse_solve
subroutine hsmg_coarse_solve(e, r)
Definition: hsmg.f:1322
hsmg_setup_mg_nx
subroutine hsmg_setup_mg_nx()
Definition: hsmg.f:1605
hsmg_tnsr1
subroutine hsmg_tnsr1(v, nv, nu, A, At)
Definition: hsmg.f:2135
hsmg_schwarz_wt2d
subroutine hsmg_schwarz_wt2d(e, wt, n)
Definition: hsmg.f:1262
hsmg_index_0
subroutine hsmg_index_0
Definition: hsmg.f:1664
hsmg_setup_intpm
subroutine hsmg_setup_intpm(jh, zf, zc, nf, nc)
Definition: hsmg.f:109
hsmg_setup_schwarz_wt
subroutine hsmg_setup_schwarz_wt(ifsqrt)
Definition: hsmg.f:1188
hsmg_do_fast
subroutine hsmg_do_fast(e, r, s, d, nl)
Definition: hsmg.f:898
hsmg_setup_fast1d_b
subroutine hsmg_setup_fast1d_b(b, lbc, rbc, ll, lm, lr, bh, n)
Definition: hsmg.f:849
hsmg_setup_fast
subroutine hsmg_setup_fast(s, d, nl, ah, bh, n)
Definition: hsmg.f:702
hsmg_schwarz_toext2d
subroutine hsmg_schwarz_toext2d(a, b, n)
Definition: hsmg.f:556
h1mg_setup_semhat
subroutine h1mg_setup_semhat
Definition: hsmg.f:2341
hsmg_setup_dssum
subroutine hsmg_setup_dssum
Definition: hsmg.f:123
h1mg_setup_schwarz_wt_2
subroutine h1mg_setup_schwarz_wt_2(wt, ie, n, work, ifsqrt)
Definition: hsmg.f:2955
hsmg_dsprod
subroutine hsmg_dsprod(u, l)
Definition: hsmg.f:343
h1mg_setup_dssum
subroutine h1mg_setup_dssum
Definition: hsmg.f:2362
hsmg_schwarz_toext3d
subroutine hsmg_schwarz_toext3d(a, b, n)
Definition: hsmg.f:581
dsavg
subroutine dsavg(u)
Definition: ic.f:1878
geom_reset
subroutine geom_reset(icall)
Definition: ic.f:1802
get_vert
subroutine get_vert
Definition: map2.f:126
invers2
subroutine invers2(a, b, n)
Definition: math.f:51
glmin
function glmin(a, n)
Definition: math.f:973
transpose
subroutine transpose(a, lda, b, ldb)
Definition: math.f:1594
glsc2
function glsc2(x, y, n)
Definition: math.f:794
vlmax
real function vlmax(vec, n)
Definition: math.f:396
copy
subroutine copy(a, b, n)
Definition: math.f:260
add2sxy
subroutine add2sxy(x, a, y, b, n)
Definition: math.f:1574
iglmax
function iglmax(a, n)
Definition: math.f:913
izero
subroutine izero(a, n)
Definition: math.f:215
vlsc2
real function vlsc2(x, y, n)
Definition: math.f:739
glsum
function glsum(x, n)
Definition: math.f:861
cmult
subroutine cmult(a, const, n)
Definition: math.f:315
rone
subroutine rone(a, n)
Definition: math.f:230
rzero
subroutine rzero(a, n)
Definition: math.f:208
glmax
function glmax(a, n)
Definition: math.f:960
mxm
subroutine mxm(a, n1, b, n2, c, n3)
Definition: mxm_wrapper.f:2
gradl_rst
subroutine gradl_rst(ur, us, ut, u, md, if3d)
Definition: navier1.f:4643
ortho
subroutine ortho(respr)
Definition: navier1.f:224
gradl_rst_t
subroutine gradl_rst_t(u, ur, us, ut, md, if3d)
Definition: navier1.f:4614
cdabdtp
subroutine cdabdtp(ap, wp, h1, h2, h2inv, intype)
Definition: navier1.f:259
avg1
subroutine avg1(avg, f, alpha, beta, n, name, ifverbose)
Definition: navier5.f:1019
avg2
subroutine avg2(avg, f, alpha, beta, n, name, ifverbose)
Definition: navier5.f:1042