KTH framework for Nek5000 toolboxes; testing version  0.0.1
convect.f
Go to the documentation of this file.
1 c-----------------------------------------------------------------------
2  subroutine char_conv(p0,u,ulag,bm,bmlag,msk,c,cs,gsl)
3 c
4 c
5 c Convect over last NBD steps using characteristics scheme
6 c
7 c NOTE: Here, we assume that ulag is stored by time-slice first,
8 c then by field number (this is opposite to prior Nek5000)
9 c
10 c
11  include 'SIZE'
12  include 'TOTAL'
13  real p0(1),u(1),ulag(1),bm(1),bmlag(1),msk(1),c(1),cs(0:1)
14  integer gsl
15 
16  common /scrns/ ct(lxd*lyd*lzd*lelv*ldim)
17 
18  common /scrvh/ bmsk(lx1*ly1*lz1*lelv)
19  $ , bdwt(lx1*ly1*lz1*lelv)
20  $ , bmst(lx1*ly1*lz1*lelv)
21  $ , u1(lx1*ly1*lz1*lelv)
22 
23  common /scrmg/ r1(lx1*ly1*lz1*lelv)
24  $ , r2(lx1*ly1*lz1*lelv)
25  $ , r3(lx1*ly1*lz1*lelv)
26  $ , r4(lx1*ly1*lz1*lelv)
27 
28  nelc = nelv ! number of elements in convecting field
29  if (ifield.eq.ifldmhd) nelc = nelfld(ifield)
30 
31  nc = cs(0) ! number of stored convecting fields
32 
33  ln = lx1*ly1*lz1*lelt
34  n = lx1*ly1*lz1*nelfld(ifield)
35  m = lxd*lyd*lzd*nelc*ldim
36 
37 c if(nid.eq.0) write(*,*) 'going into char_conv1 '
38  call char_conv1 (p0,u,bmnv,n,ulag,ln,gsl,c,m,cs(1),nc,ct
39  $ ,u1,r1,r2,r3,r4,bmsk,bdivw,bdwt,bmass,bmst,bm,bmlag)
40 
41  return
42  end
43 c-----------------------------------------------------------------------
44  subroutine char_conv1 (p0,u,bmnv,n,ulag,ln,gsl,c,m,cs,nc,ct
45  $ ,u1,r1,r2,r3,r4,bmsk,bdivw,bdwt,bmass,bmst,bm,bmlag)
46 
47  include 'SIZE'
48  include 'INPUT'
49  include 'TSTEP'
50 
51  real p0(n),u(n),bmnv(n,1),ulag(ln,1),c(m,0:nc),cs(0:nc),bdivw(n,1)
52  $ ,bm(n), bmlag(ln,1)
53 
54  real ct(m),u1(n),r1(n),r2(n),r3(n),r4(n),bmsk(n),bdwt(n) ! work arrays
55 
56  integer gsl
57 
58 
59 ! Convect over last NBD steps using characteristics scheme
60 
61 ! n-q n-1
62 ! Given u(t , X ), q = 1,2,...,nbd, compute phi ( := p0 )
63 
64 ! n-1 nbd ~n-q
65 ! phi := sum u
66 ! q=1
67 
68 ! ~n-q ~n-q
69 ! each u satisfies u := v such that
70 
71 
72 ! dv
73 ! -- + C.grad v = 0 t \in [t^n-q,t^n], v(t^n-q,X) = u(t^n-q,X)
74 ! dt
75 
76 ! n = lx1*ly1*lz1*nelv
77 ! m = lxd*lyd*lzd*nelv
78 
79  tau = time-vlsum(dtlag,nbd) ! initialize time for u^n-k
80  call int_vel (ct ,tau,c ,m,nc,cs,nid) ! ct(t) = sum w_k c(.,k)
81  call int_vel (bmsk,tau,bmnv ,n,nc,cs,nid) ! B^-1(t^n-1)
82  call int_vel (bmst,tau,bmass,n,nc,cs,nid) ! B(t^n-1)
83  call int_vel (bdwt,tau,bdivw,n,nc,cs,nid) ! BdivW(t^n-1)
84 
85  call rzero(p0,n)
86 
87  do ilag = nbd,1,-1
88 
89  um = 0
90  if (ilag.eq.1) then
91  do i=1,n
92  p0(i) = p0(i)+bd(ilag+1)*u(i)*bm(i)
93  um=max(um,u(i))
94  enddo
95  else
96  if(ifmvbd) then
97  do i=1,n
98  p0(i) = p0(i)+bd(ilag+1)*ulag(i,ilag-1)*bmlag(i,ilag-1)
99  um=max(um,ulag(i,ilag-1))
100  enddo
101  else
102  do i=1,n
103  p0(i) = p0(i)+bd(ilag+1)*ulag(i,ilag-1)*bm(i)
104  um=max(um,ulag(i,ilag-1))
105  enddo
106  endif
107  endif
108 
109  dtau = dtlag(ilag)/ntaubd
110  do itau = 1,ntaubd ! ntaubd=number of RK4 substeps (typ. 1 or 2)
111 
112  tau1 = tau + dtau
113 
114  c1 = 1.
115  c2 = -dtau/2.
116  c3 = -dtau
117  th = tau+dtau/2.
118 
119  call invcol3 (u1,p0,bmst,n)
120  call conv_rhs(r1,u1,ct,bmsk,bmst,bdwt,gsl) ! STAGE 1
121  call col2 (r1,bmst,n) ! ! r1 = B(n-1)* r1
122 
123  call add3s12 (u1,p0,r1,c1,c2,n)
124  call int_vel (bmst,th,bmass,n,nc,cs,nid) ! B(n-1/2)
125  call invcol2 (u1,bmst,n) ! u2=B(n-1/2)
126 
127  call int_vel (ct ,th,c ,m,nc,cs,nid) ! STAGE 2
128  call int_vel (bmsk,th,bmnv ,n,nc,cs,nid) ! B^-1(n-1/2)
129  call int_vel (bdwt,th,bdivw,n,nc,cs,nid) ! BdivW(n-1/2)
130  call conv_rhs(r2,u1,ct,bmsk,bmst,bdwt,gsl)
131  call col2 (r2,bmst,n) ! du = B * du
132 
133  call add3s12 (u1,p0,r2,c1,c2,n) ! STAGE 3
134  call invcol2 (u1,bmst,n)
135  call conv_rhs(r3,u1,ct,bmsk,bmst,bdwt,gsl) ! B(n-1/2) (still)
136  call col2 (r3,bmst,n) ! du = B * du
137 
138  call add3s12 (u1,p0,r3,c1,c3,n)
139  call int_vel (bmst,tau1,bmass,n,nc,cs,nid) ! B^-1(n)
140  call invcol2 (u1,bmst,n) ! u2=B(n-1/2)
141 
142  call int_vel (ct ,tau1,c ,m,nc,cs,nid) ! STAGE 4
143  call int_vel (bmsk,tau1,bmnv ,n,nc,cs,nid) ! B^-1(n)
144  call int_vel (bdwt,tau1,bdivw,n,nc,cs,nid) ! BdivW(n)
145  call conv_rhs(r4,u1,ct,bmsk,bmst,bdwt,gsl)
146  call col2 (r4,bmst,n) ! du = B * du
147 
148  c1 = -dtau/6.
149  c2 = -dtau/3.
150  do i=1,n
151  p0(i) = p0(i)+c1*(r1(i)+r4(i))+c2*(r2(i)+r3(i))
152  enddo
153  tau = tau1
154  enddo
155  enddo
156 
157  return
158  end
159 c-----------------------------------------------------------------------
160  subroutine int_vel(c_t,t0,c,n,nc,ct,nid)
161 
162 c Interpolate convecting velocity field c(t_1,...,t_nconv) to
163 c time t0 and return result in c_t.
164 
165 c Ouput: c_t = sum wt_k * ct_i(k)
166 
167 c Here, t0 is the time of interest
168 
169  real c_t(n),c(n,0:nc),ct(0:nc)
170 
171  parameter(lwtmax=10)
172  real wt(0:lwtmax)
173 
174  if (nc.gt.lwtmax) then
175  write(6,*) nid,'ERROR int_vel: lwtmax too small',lwtmax,nc
176  call exitt
177  endif
178 
179  no = nc-1
180  call fd_weights_full(t0,ct(0),no,0,wt) ! interpolation weights
181 
182  call rzero(c_t,n)
183  do j=1,n
184  do i=0,no
185  c_t(j) = c_t(j) + wt(i)*c(j,i)
186  enddo
187  enddo
188 
189  return
190  end
191 c-----------------------------------------------------------------------
192  subroutine conv_rhs (du,u,c,bmsk,bmst,bdwt,gsl)
193 c
194  include 'SIZE'
195  include 'TOTAL'
196 c
197 c apply convecting field c(1,ldim) to scalar field u(1)
198 c
199  real du(1),u(1),c(1),bmsk(1),bdwt(1)
200  integer gsl
201 c
202  logical ifconv
203 c
204 c ifconv = .false.
205  ifconv = .true.
206 c
207  n = lx1*ly1*lz1*nelv
208 
209  if (ifdgfld(ifield)) then
210 
211  if (param(99).eq.1) call conv_rhs_dg (du,u,c)
212  if (param(99).eq.0) call conv_rhs_dg_aliased (du,u,c)
213 
214  elseif (ifconv) then
215 
216  if (ifcons) then
217  if (if3d ) call convop_cons_3d (du,u,c,lx1,lxd,nelv)
218  if (.not.if3d) call convop_cons_2d (du,u,c,lx1,lxd,nelv)
219  else
220  if (if3d ) call convop_fst_3d (du,u,c,lx1,lxd,nelv)
221  if (.not.if3d) call convop_fst_2d (du,u,c,lx1,lxd,nelv)
222  endif
223 
224  call subcol3(du,bdwt,u,n)
225  call fgslib_gs_op (gsl,du,1,1,0) ! +
226 
227  call col2 (du,bmsk,n) ! du = Binv * msk * du
228 
229  else
230  call rzero (du,n)
231  endif
232 
233  return
234  end
235 c-----------------------------------------------------------------------
236  subroutine convop_fst_3d(du,u,c,mx,md,nel)
237 c
238  include 'SIZE'
239 c
240 c apply convecting field c to scalar field u
241 c
242  real du(mx*mx*mx,nel)
243  real u(mx*mx*mx,nel)
244  real c(md*md*md,nel,3)
245  parameter(ldd=lxd*lyd*lzd)
246  common /ctmp1/ ur(ldd),us(ldd),ut(ldd),ud(ldd)
247 c
248  logical if3d,ifd
249  integer e
250 c
251  if3d = .true.
252  ifd = .false.
253  if (md.ne.mx) ifd=.true.
254  ifd=.true.
255 c
256  nrstd = md**3
257  call lim_chk(nrstd,ldd,'urus ','ldd ','convop_fst')
258 c
259  do e=1,nel
260  call grad_rstd(ur,us,ut,u(1,e),mx,md,if3d,ud)
261  if (ifd) then ! dealiased
262  do i=1,nrstd
263 c C has the mass matrix factored in per (4.8.5), p. 227, DFM.
264  ud(i) = c(i,e,1)*ur(i)+c(i,e,2)*us(i)+c(i,e,3)*ut(i)
265  enddo
266  call intp_rstd(du(1,e),ud,mx,md,if3d,1) ! 1 --> transpose
267  else
268  do i=1,nrstd
269 c C has the mass matrix factored in per (4.8.5), p. 227, DFM.
270  du(i,e) = c(i,e,1)*ur(i)+c(i,e,2)*us(i)+c(i,e,3)*ut(i)
271  enddo
272  endif
273  enddo
274 c
275  return
276  end
277 c-----------------------------------------------------------------------
278  subroutine convop_fst_2d(du,u,c,mx,md,nel)
279 c
280  include 'SIZE'
281 c
282 c apply convecting field c to scalar field u
283 c
284  real du(mx*mx,nel)
285  real u(mx*mx,nel)
286  real c(md*md,nel,2)
287  parameter(ldd=lxd*lyd*lzd)
288  common /ctmp1/ ur(ldd),us(ldd),ut(ldd),ud(ldd)
289 c
290  logical if3d,ifd
291  integer e
292 c
293  if3d = .false.
294  ifd = .false.
295  if (md.ne.mx) ifd=.true.
296  ifd=.true.
297 c
298  nrstd = md**2
299  call lim_chk(nrstd,ldd,'urus ','ldd ','convop_fst')
300 c
301  do e=1,nel
302  call grad_rstd(ur,us,ut,u(1,e),mx,md,if3d,ud)
303  if (ifd) then ! dealiased
304  do i=1,nrstd
305 c C has the mass matrix factored in per (4.8.5), p. 227, DFM.
306  ud(i) = c(i,e,1)*ur(i)+c(i,e,2)*us(i)
307  enddo
308  call intp_rstd(du(1,e),ud,mx,md,if3d,1) ! 1 --> transpose
309  else
310  do i=1,nrstd
311 c C has the mass matrix factored in per (4.8.5), p. 227, DFM.
312  du(i,e) = c(i,e,1)*ur(i)+c(i,e,2)*us(i)
313  enddo
314  endif
315  enddo
316 c
317  return
318  end
319 c-----------------------------------------------------------------------
320  subroutine grad_rstd(ur,us,ut,u,mx,md,if3d,ju) ! GLL->GL grad
321 
322  include 'SIZE'
323  include 'DXYZ'
324 
325  real ur(1),us(1),ut(1),u(1),ju(1)
326  logical if3d
327 
328  parameter(ldg=lxd**3,lwkd=4*lxd*lxd)
329  common /dgrad/ d(ldg),dt(ldg),dg(ldg),dgt(ldg),jgl(ldg),jgt(ldg)
330  $ , wkd(lwkd)
331  real jgl,jgt
332 
333  call intp_rstd(ju,u,mx,md,if3d,0) ! 0 = forward
334 
335  m0 = md-1
336  call get_dgl_ptr (ip,md,md)
337  if (if3d) then
338  call local_grad3(ur,us,ut,ju,m0,1,dg(ip),dgt(ip))
339  else
340  call local_grad2(ur,us ,ju,m0,1,dg(ip),dgt(ip))
341  endif
342 
343  return
344  end
345 c-----------------------------------------------------------------------
346  subroutine intp_rstd(ju,u,mx,md,if3d,idir) ! GLL->GL interpolation
347 
348 c GLL interpolation from mx to md.
349 
350 c If idir ^= 0, then apply transpose operator (md to mx)
351 
352  include 'SIZE'
353 
354  real ju(1),u(1)
355  logical if3d
356 
357  parameter(ldg=lxd**3,lwkd=4*lxd*lxd)
358  common /dgrad/ d(ldg),dt(ldg),dg(ldg),dgt(ldg),jgl(ldg),jgt(ldg)
359  $ , wkd(lwkd)
360  real jgl,jgt
361 
362  parameter(ld=2*lxd)
363  common /ctmp0/ w(ld**ldim,2)
364 
365  call lim_chk(md,ld,'md ','ld ','grad_rstd ')
366  call lim_chk(mx,ld,'mx ','ld ','grad_rstd ')
367 
368  ldw = 2*(ld**ldim)
369 
370  call get_int_ptr (i,mx,md)
371 c
372  if (idir.eq.0) then
373  call specmpn(ju,md,u,mx,jgl(i),jgt(i),if3d,w,ldw)
374  else
375  call specmpn(ju,mx,u,md,jgt(i),jgl(i),if3d,w,ldw)
376  endif
377 c
378  return
379  end
380 c-----------------------------------------------------------------------
381  subroutine gen_int(jgl,jgt,mp,np,w)
382 c
383 c Generate interpolation from np GLL points to mp GL points
384 c
385 c jgl = interpolation matrix, mapping from velocity nodes to pressure
386 c jgt = transpose of interpolation matrix
387 c w = work array of size (np+mp)
388 c
389 c np = number of points on GLL grid
390 c mp = number of points on GL grid
391 c
392 c
393  real jgl(mp,np),jgt(np*mp),w(1)
394 c
395  iz = 1
396  id = iz + np
397 c
398  call zwgll (w(iz),jgt,np)
399  call zwgl (w(id),jgt,mp)
400 c
401  n = np-1
402  do i=1,mp
403  call fd_weights_full(w(id+i-1),w(iz),n,0,jgt)
404  do j=1,np
405  jgl(i,j) = jgt(j) ! Interpolation matrix
406  enddo
407  enddo
408 c
409  call transpose(jgt,np,jgl,mp)
410 c
411  return
412  end
413 c-----------------------------------------------------------------------
414  subroutine gen_dgl(dgl,dgt,mp,np,w)
415 c
416 c Generate derivative from np GL points onto mp GL points
417 c
418 c dgl = interpolation matrix, mapping from velocity nodes to pressure
419 c dgt = transpose of interpolation matrix
420 c w = work array of size (3*np+mp)
421 c
422 c np = number of points on GLL grid
423 c mp = number of points on GL grid
424 c
425 c
426 c
427  real dgl(mp,np),dgt(np*mp),w(1)
428 c
429 c
430  iz = 1
431  id = iz + np
432 c
433  call zwgl (w(iz),dgt,np) ! GL points
434  call zwgl (w(id),dgt,mp) ! GL points
435 c
436  ndgt = 2*np
437  ldgt = mp*np
438  call lim_chk(ndgt,ldgt,'ldgt ','dgt ','gen_dgl ')
439 c
440  n = np-1
441  do i=1,mp
442  call fd_weights_full(w(id+i-1),w(iz),n,1,dgt) ! 1=1st deriv.
443  do j=1,np
444  dgl(i,j) = dgt(np+j) ! Derivative matrix
445  enddo
446  enddo
447 c
448  call transpose(dgt,np,dgl,mp)
449 c
450  return
451  end
452 c-----------------------------------------------------------------------
453  subroutine lim_chk(n,m,avar5,lvar5,sub_name10)
454  include 'SIZE' ! need nid
455  character*5 avar5,lvar5
456  character*10 sub_name10
457 
458  if (n.gt.m) then
459  write(6,1) nid,n,m,avar5,lvar5,sub_name10
460  1 format(i8,' ERROR: :',2i12,2(1x,a5),1x,a10)
461  call exitti('lim_chk problem. $',n)
462  endif
463 
464  return
465  end
466 c-----------------------------------------------------------------------
467  subroutine get_int_ptr (ip,mx,md) ! GLL-->GL pointer
468 
469 c Get pointer to jgl() for interpolation pair (mx,md)
470 
471  include 'SIZE'
472 
473  parameter(ldg=lxd**3,lwkd=4*lxd*lxd)
474  common /dgrad/ d(ldg),dt(ldg),dg(ldg),dgt(ldg),jgl(ldg),jgt(ldg)
475  $ , wkd(lwkd)
476  real jgl,jgt
477 c
478  parameter(ld=2*lxd)
479  common /igrad/ pd(0:ld*ld)
480  $ , pdg(0:ld*ld)
481  $ , pjgl(0:ld*ld)
482  integer pd , pdg , pjgl
483 c
484  ij = md + ld*(mx-1)
485  ip = pjgl(ij)
486 c
487  if (ip.eq.0) then
488 c
489  nstore = pjgl(0)
490  pjgl(ij) = nstore+1
491  nstore = nstore + md*mx
492  pjgl(0) = nstore
493  ip = pjgl(ij)
494 c
495  nwrkd = mx + md
496  call lim_chk(nstore,ldg ,'jgl ','ldg ','get_int_pt')
497  call lim_chk(nwrkd ,lwkd,'wkd ','lwkd ','get_int_pt')
498 c
499  call gen_int(jgl(ip),jgt(ip),md,mx,wkd)
500  endif
501 c
502  return
503  end
504 c-----------------------------------------------------------------------
505  subroutine get_dgl_ptr (ip,mx,md)
506 c
507 c Get pointer to GL-GL interpolation dgl() for pair (mx,md)
508 c
509  include 'SIZE'
510 c
511  parameter(ldg=lxd**3,lwkd=4*lxd*lxd)
512  common /dgrad/ d(ldg),dt(ldg),dg(ldg),dgt(ldg),jgl(ldg),jgt(ldg)
513  $ , wkd(lwkd)
514  real jgl,jgt
515 c
516  parameter(ld=2*lxd)
517  common /jgrad/ pd(0:ld*ld)
518  $ , pdg(0:ld*ld)
519  $ , pjgl(0:ld*ld)
520  integer pd , pdg , pjgl
521 c
522  ij = md + ld*(mx-1)
523  ip = pdg(ij)
524 
525  if (ip.eq.0) then
526 
527  nstore = pdg(0)
528  pdg(ij) = nstore+1
529  nstore = nstore + md*mx
530  pdg(0) = nstore
531  ip = pdg(ij)
532 c
533  nwrkd = mx + md
534  call lim_chk(nstore,ldg ,'dg ','ldg ','get_dgl_pt')
535  call lim_chk(nwrkd ,lwkd,'wkd ','lwkd ','get_dgl_pt')
536 c
537  call gen_dgl(dg(ip),dgt(ip),md,mx,wkd)
538  endif
539 c
540  return
541  end
542 c-----------------------------------------------------------------------
543  subroutine set_conv_char(ct,c,ux,uy,uz,nelc,tau,ifnew)
544  include 'SIZE'
545  include 'TSTEP'
546 
547  real ct(0:1) ! time stamps for saved field (0=#flds)
548  real c(1) ! saved vel. fields, dealiased etc.
549  real ux(1),uy(1),uz(1) ! input vel. field
550  integer nelc ! number of elements in conv. field
551  logical ifnew ! =true if shifting stack of fields
552 
553  numr = lxd*lyd*lzd*lelv*ldim*(lorder+1)
554  denr = lxd*lyd*lzd*nelv*ldim
555  nconv_max = numr/denr
556  if (nconv_max.lt.nbdinp+1)
557  $ call exitti(
558  $ 'ABORT: not enough memory for characteristics scheme!$',
559  $ nconv_max)
560 
561  nc = ct(0)
562 
563  m = lxd*lyd*lzd*nelc*ldim
564 
565  call set_ct_cvx
566  $ (ct,c,m,ux,uy,uz,tau,nc,nconv_max,nelc,ifnew)
567 
568  nc = min(nc,nbdinp)
569  ct(0) = nc ! store current count
570 
571  return
572  end
573 c-----------------------------------------------------------------------
574  subroutine set_ct_cvx(ct,c,m,u,v,w,tau,nc,mc,nelc,ifnew)
575  include 'SIZE'
576  include 'INPUT' ! ifcons
577 
578  real ct(0:1),c(m,1)
579  real u(1),v(1),w(1)
580  logical ifnew
581 
582  if (ifnew) then
583 
584 c Shift existing convecting fields
585 c Note: "1" entry is most recent
586 
587  nc = nc+1
588  nc = min(nc,mc)
589  ct(0) = nc
590 
591  do i=nc,2,-1
592  call copy(c(1,i),c(1,i-1),m)
593  ct(i) = ct(i-1)
594  enddo
595  endif
596 
597 c Save time and map the current velocity to rst coordinates.
598 
599  ix = 1
600  iy = ix + lxd*lyd*lzd*nelc
601  iz = iy + lxd*lyd*lzd*nelc
602 
603  if (ifcons) then
604  call set_convect_cons(c(ix,1),c(iy,1),c(iz,1),u,v,w)
605  else
606  call set_convect_new (c(ix,1),c(iy,1),c(iz,1),u,v,w)
607  endif
608 
609  ct(1) = tau
610 
611  return
612  end
613 c-----------------------------------------------------------------------
614  subroutine grad_rst(ur,us,ut,u,md,if3d) ! Gauss-->Gauss grad
615 
616  include 'SIZE'
617  include 'DXYZ'
618 
619  real ur(1),us(1),ut(1),u(1)
620  logical if3d
621 
622  parameter(ldg=lxd**3,lwkd=4*lxd*lxd)
623  common /dgrad/ d(ldg),dt(ldg),dg(ldg),dgt(ldg),jgl(ldg),jgt(ldg)
624  $ , wkd(lwkd)
625  real jgl,jgt
626 
627  m0 = md-1
628  call get_dgl_ptr (ip,md,md)
629  if (if3d) then
630  call local_grad3(ur,us,ut,u,m0,1,dg(ip),dgt(ip))
631  else
632  call local_grad2(ur,us ,u,m0,1,dg(ip),dgt(ip))
633  endif
634 
635  return
636  end
637 c-----------------------------------------------------------------------
638  subroutine convect_new(bdu,u,ifuf,cx,cy,cz,ifcf)
639 
640 C Compute dealiased form: J^T Bf *JC .grad Ju w/ correct Jacobians
641 C
642  include 'SIZE'
643  include 'TOTAL'
644 
645  real bdu(1),u(1),cx(1),cy(1),cz(1)
646  logical ifuf,ifcf ! u and/or c already on fine mesh?
647 
648  parameter(lxy=lx1*ly1*lz1,ltd=lxd*lyd*lzd)
649  common /scrcv/ fx(ltd),fy(ltd),fz(ltd)
650  $ , ur(ltd),us(ltd),ut(ltd)
651  $ , tr(ltd,3),uf(ltd)
652 
653  integer e
654 
655  call set_dealias_rx
656 
657  nxyz1 = lx1*ly1*lz1
658  nxyzd = lxd*lyd*lzd
659 
660  nxyzu = nxyz1
661  if (ifuf) nxyzu = nxyzd
662 
663  nxyzc = nxyz1
664  if (ifcf) nxyzc = nxyzd
665 
666  iu = 1 ! pointer to scalar field u
667  ic = 1 ! pointer to vector field C
668  ib = 1 ! pointer to scalar field Bdu
669 
670 
671  do e=1,nelv
672 
673  if (ifcf) then
674 
675  call copy(tr(1,1),cx(ic),nxyzd) ! already in rst form
676  call copy(tr(1,2),cy(ic),nxyzd)
677  if (if3d) call copy(tr(1,3),cz(ic),nxyzd)
678 
679  else ! map coarse velocity to fine mesh (C-->F)
680 
681  call intp_rstd(fx,cx(ic),lx1,lxd,if3d,0) ! 0 --> forward
682  call intp_rstd(fy,cy(ic),lx1,lxd,if3d,0) ! 0 --> forward
683  if (if3d) call intp_rstd(fz,cz(ic),lx1,lxd,if3d,0) ! 0 --> forward
684 
685  if (if3d) then ! Convert convector F to r-s-t coordinates
686 
687  do i=1,nxyzd
688  tr(i,1)=rx(i,1,e)*fx(i)+rx(i,2,e)*fy(i)+rx(i,3,e)*fz(i)
689  tr(i,2)=rx(i,4,e)*fx(i)+rx(i,5,e)*fy(i)+rx(i,6,e)*fz(i)
690  tr(i,3)=rx(i,7,e)*fx(i)+rx(i,8,e)*fy(i)+rx(i,9,e)*fz(i)
691  enddo
692 
693  else
694 
695  do i=1,nxyzd
696  tr(i,1)=rx(i,1,e)*fx(i)+rx(i,2,e)*fy(i)
697  tr(i,2)=rx(i,3,e)*fx(i)+rx(i,4,e)*fy(i)
698  enddo
699 
700  endif
701 
702  endif
703 
704  if (ifuf) then
705  call grad_rst(ur,us,ut,u(iu),lxd,if3d)
706  else
707  call intp_rstd(uf,u(iu),lx1,lxd,if3d,0) ! 0 --> forward
708  call grad_rst(ur,us,ut,uf,lxd,if3d)
709  endif
710 
711  if (if3d) then
712  do i=1,nxyzd ! mass matrix included, per DFM (4.8.5)
713  uf(i) = tr(i,1)*ur(i)+tr(i,2)*us(i)+tr(i,3)*ut(i)
714  enddo
715  else
716  do i=1,nxyzd ! mass matrix included, per DFM (4.8.5)
717  uf(i) = tr(i,1)*ur(i)+tr(i,2)*us(i)
718  enddo
719  endif
720  call intp_rstd(bdu(ib),uf,lx1,lxd,if3d,1) ! Project back to coarse
721 
722  ic = ic + nxyzc
723  iu = iu + nxyzu
724  ib = ib + nxyz1
725 
726  enddo
727 
728  return
729  end
730 c-----------------------------------------------------------------------
731  subroutine convect_cons(bdu,u,ifuf,cx,cy,cz,ifcf)
732 
733 c Compute dealiased form: J^T Bf *div. JC Ju w/ correct Jacobians
734 
735 c conservative form
736 
737 
738  include 'SIZE'
739  include 'TOTAL'
740 
741  real bdu(1),u(1),cx(1),cy(1),cz(1)
742 
743  logical ifuf,ifcf ! u and/or c already on fine mesh?
744 
745  parameter(lxy=lx1*ly1*lz1,ltd=lxd*lyd*lzd)
746  common /scrcv/ uf(ltd),cf(ltd),cu(ltd)
747  $ , cr(ltd),cs(ltd),ct(ltd)
748 
749 
750  integer e
751 
752  call set_dealias_rx
753 
754  nxyz1 = lx1*ly1*lz1
755  nxyzd = lxd*lyd*lzd
756 
757  nxyzu = nxyz1
758  if (ifuf) nxyzu = nxyzd
759 
760  nxyzc = nxyz1
761  if (ifcf) nxyzc = nxyzd
762 
763  iu = 1 ! pointer to scalar field u
764  ic = 1 ! pointer to vector field C
765  ib = 1 ! pointer to scalar field Bdu
766 
767  do e=1,nelv
768 
769  call intp_rstd(uf,u(iu),lx1,lxd,if3d,0) ! 0 --> forward
770 
771  call rzero(cu,nxyzd)
772  do i=1,ldim
773 
774  if (ifcf) then ! C is already on fine mesh
775 
776  call exitt ! exit for now
777 
778  else ! map coarse velocity to fine mesh (C-->F)
779 
780  if (i.eq.1) call intp_rstd(cf,cx(ic),lx1,lxd,if3d,0) ! 0 --> forward
781  if (i.eq.2) call intp_rstd(cf,cy(ic),lx1,lxd,if3d,0) ! 0 --> forward
782  if (i.eq.3) call intp_rstd(cf,cz(ic),lx1,lxd,if3d,0) ! 0 --> forward
783 
784  call col2(cf,uf,nxyzd) ! collocate C and u on fine mesh
785 
786  call grad_rst(cr,cs,ct,cf,lxd,if3d) ! d/dr (C_i*u)
787 
788  if (if3d) then
789 
790  do j=1,nxyzd
791  cu(j)=cu(j)
792  $ +cr(j)*rx(j,i,e)+cs(j)*rx(j,i+3,e)+ct(j)*rx(j,i+6,e)
793  enddo
794 
795  else ! 2D
796 
797  do j=1,nxyzd
798  cu(j)=cu(j)
799  $ +cr(j)*rx(j,i,e)+cs(j)*rx(j,i+2,e)
800  enddo
801 
802  endif
803  endif
804  enddo
805 
806  call intp_rstd(bdu(ib),cu,lx1,lxd,if3d,1) ! Project back to coarse
807 
808  ic = ic + nxyzc
809  iu = iu + nxyzu
810  ib = ib + nxyz1
811 
812  enddo
813 
814  return
815  end
816 c-----------------------------------------------------------------------
817  subroutine set_convect_cons(cx,cy,cz,ux,uy,uz)
818 
819 c Put vx,vy,vz on fine mesh (for conservation form)
820 
821 
822  include 'SIZE'
823  include 'TOTAL'
824 
825  parameter(lxy=lx1*ly1*lz1,ltd=lxd*lyd*lzd)
826 
827  real cx(ltd,1),cy(ltd,1),cz(ltd,1)
828  real ux(lxy,1),uy(lxy,1),uz(lxy,1)
829 
830  integer e
831 
832  call set_dealias_rx
833 
834  do e=1,nelv ! Map coarse velocity to fine mesh (C-->F)
835 
836  call intp_rstd(cx(1,e),ux(1,e),lx1,lxd,if3d,0) ! 0 --> forward
837  call intp_rstd(cy(1,e),uy(1,e),lx1,lxd,if3d,0) ! 0 --> forward
838  if (if3d) call intp_rstd(cz(1,e),uz(1,e),lx1,lxd,if3d,0) ! 0 --> forward
839 
840  enddo
841 
842  return
843  end
844 c-----------------------------------------------------------------------
845  subroutine set_convect_new(cr,cs,ct,ux,uy,uz)
846 C
847 C Put vxd,vyd,vzd into rst form on fine mesh
848 C
849 C For rst form, see eq. (4.8.5) in Deville, Fischer, Mund (2002).
850 C
851  include 'SIZE'
852  include 'TOTAL'
853 
854  parameter(lxy=lx1*ly1*lz1,ltd=lxd*lyd*lzd)
855 
856  real cr(ltd,1),cs(ltd,1),ct(ltd,1)
857  real ux(lxy,1),uy(lxy,1),uz(lxy,1)
858 
859  common /scrcv/ fx(ltd),fy(ltd),fz(ltd)
860  $ , ur(ltd),us(ltd),ut(ltd)
861  $ , tr(ltd,3),uf(ltd)
862 
863  integer e
864 
865  call set_dealias_rx
866 
867  nxyz1 = lx1*ly1*lz1
868  nxyzd = lxd*lyd*lzd
869 
870  ic = 1 ! pointer to vector field C
871 
872  do e=1,nelv
873 
874 c Map coarse velocity to fine mesh (C-->F)
875 
876  call intp_rstd(fx,ux(1,e),lx1,lxd,if3d,0) ! 0 --> forward
877  call intp_rstd(fy,uy(1,e),lx1,lxd,if3d,0) ! 0 --> forward
878  if (if3d) call intp_rstd(fz,uz(1,e),lx1,lxd,if3d,0) ! 0 --> forward
879 
880 c Convert convector F to r-s-t coordinates
881 
882  if (if3d) then
883 
884  do i=1,nxyzd
885  cr(i,e)=rx(i,1,e)*fx(i)+rx(i,2,e)*fy(i)+rx(i,3,e)*fz(i)
886  cs(i,e)=rx(i,4,e)*fx(i)+rx(i,5,e)*fy(i)+rx(i,6,e)*fz(i)
887  ct(i,e)=rx(i,7,e)*fx(i)+rx(i,8,e)*fy(i)+rx(i,9,e)*fz(i)
888  enddo
889 
890  else
891 
892  do i=1,nxyzd
893  cr(i,e)=rx(i,1,e)*fx(i)+rx(i,2,e)*fy(i)
894  cs(i,e)=rx(i,3,e)*fx(i)+rx(i,4,e)*fy(i)
895  enddo
896 
897  endif
898  enddo
899 
900  return
901  end
902 c-----------------------------------------------------------------------
903  subroutine advchar
904 c
905 c Compute convective contribution using
906 c operator-integrator-factor method (characteristics).
907 c
908  include 'SIZE'
909  include 'MASS'
910  include 'INPUT'
911  include 'SOLN'
912  include 'TSTEP'
913  include 'PARALLEL'
914 
915  include 'CTIMER'
916 
917  common /cchar/ ct_vx(0:lorder) ! time for each slice in c_vx()
918 
919  common /scruz/ phx(lx1*ly1*lz1*lelt)
920  $ , phy(lx1*ly1*lz1*lelt)
921  $ , phz(lx1*ly1*lz1*lelt)
922  $ , hmsk(lx1*ly1*lz1*lelt)
923 
924  if (icalld.eq.0) tadvc=0.0
925  icalld=icalld+1
926  nadvc=icalld
927  etime1=dnekclock()
928 
929  dti = 1./dt
930  n = lx1*ly1*lz1*nelv
931 
932  call char_conv(phx,vx,vxlag,bm1,bm1lag,hmsk,c_vx,ct_vx,gsh_fld(1))
933  call char_conv(phy,vy,vylag,bm1,bm1lag,hmsk,c_vx,ct_vx,gsh_fld(1))
934  if (if3d) call char_conv
935  $ (phz,vz,vzlag,bm1,bm1lag,hmsk,c_vx,ct_vx,gsh_fld(1))
936 
937  call cfill(hmsk,dti,n)
938  if(.not. iflomach) call col2(hmsk,vtrans,n)
939 
940  if (if3d) then
941 
942  do i=1,n
943  h2i = hmsk(i)
944  bfx(i,1,1,1) = bfx(i,1,1,1)+phx(i)*h2i
945  bfy(i,1,1,1) = bfy(i,1,1,1)+phy(i)*h2i
946  bfz(i,1,1,1) = bfz(i,1,1,1)+phz(i)*h2i
947  enddo
948 
949  else
950 
951  do i=1,n
952  h2i = hmsk(i)
953  bfx(i,1,1,1) = bfx(i,1,1,1)+phx(i)*h2i
954  bfy(i,1,1,1) = bfy(i,1,1,1)+phy(i)*h2i
955  enddo
956 
957  endif
958 
959  tadvc=tadvc+(dnekclock()-etime1)
960 
961  return
962  end
963 c-----------------------------------------------------------------------
964  subroutine convch
965 
966 c Compute convective contribution using
967 c operator-integrator-factor method (characteristics).
968 
969  include 'SIZE'
970  include 'MASS'
971  include 'INPUT'
972  include 'SOLN'
973  include 'TSTEP'
974  include 'PARALLEL'
975  include 'CTIMER'
976 
977  common /cchar/ ct_vx(0:lorder) ! time for each slice in c_vx()
978 
979  common /scruz/ phi(lx1*ly1*lz1*lelt)
980  $ , hmsk(lx1*ly1*lz1*lelt)
981 
982  if (icalld.eq.0) tadvc=0.0
983  icalld=icalld+1
984  nadvc=icalld
985  etime1=dnekclock()
986 
987  n = lx1*ly1*lz1*nelv
988  dti = 1./dt
989 
990  if(nid.eq.0 .and. loglevel.gt.2) write(6,*) 'convch', ifield
991  call char_conv(phi,t(1,1,1,1,ifield-1),tlag(1,1,1,1,1,ifield-1)
992  $ ,bm1,bm1lag,hmsk,c_vx,ct_vx,gsh_fld(1))
993 
994  do i=1,n
995  bq(i,1,1,1,ifield-1) = bq(i,1,1,1,ifield-1)
996  $ + phi(i)*vtrans(i,1,1,1,ifield)*dti
997  enddo
998 
999  tadvc=tadvc+(dnekclock()-etime1)
1000 
1001  return
1002  end
1003 c-----------------------------------------------------------------------
1004  subroutine convop_cons_3d(du,u,c,mx,md,nel) ! Conservation form
1005 
1006 c Apply convecting field c to scalar field u, conservation form d/dxj cj phi
1007 
1008 c Assumes that current convecting field is on dealias mesh, in c()
1009 
1010  include 'SIZE'
1011  include 'DEALIAS'
1012  include 'GEOM'
1013 
1014  real du(mx*mx*mx,nel)
1015  real u(mx*mx*mx,nel)
1016  real c(md*md*md,nel,3)
1017  parameter(ldd=lxd*lyd*lzd)
1018  common /ctmp1/ ur(ldd),us(ldd),ut(ldd),ju(ldd),ud(ldd),tu(ldd)
1019  real ju
1020 
1021  logical if3d,ifd
1022  integer e
1023 
1024  if3d = .true.
1025  ifd = .false.
1026  if (md.ne.mx) ifd=.true.
1027  ifd=.true.
1028 
1029  nrstd = md**3
1030  call lim_chk(nrstd,ldd,'urus ','ldd ','convp_cons')
1031 
1032  do e=1,nel
1033 
1034  call intp_rstd (ju,u(1,e),mx,md,if3d,0) ! 0 = forward; on Gauss points!
1035  call rzero (ud,nrstd)
1036 
1037  do j=1,ldim
1038  do i=1,nrstd
1039  tu(i)=c(i,e,j)*ju(i) ! C_j*T
1040  enddo
1041  call grad_rst(ur,us,ut,tu,md,if3d) ! Already on fine (Gauss) mesh
1042 
1043  j0 = j+0
1044  j3 = j+3
1045  j6 = j+6
1046  do i=1,nrstd ! rx has mass matrix and Jacobian on fine mesh
1047  ud(i)=ud(i)
1048  $ +rx(i,j0,e)*ur(i)+rx(i,j3,e)*us(i)+rx(i,j6,e)*ut(i)
1049  enddo
1050  enddo
1051 
1052  call intp_rstd(du(1,e),ud,mx,md,if3d,1) ! 1 --> transpose
1053 
1054  enddo
1055 
1056  return
1057  end
1058 c-----------------------------------------------------------------------
1059  subroutine convop_cons_2d(du,u,c,mx,md,nel) ! Conservation form
1060 
1061 c Apply convecting field c to scalar field u, conservation form d/dxj cj phi
1062 
1063 c Assumes that current convecting field is on dealias mesh, in c()
1064 
1065  include 'SIZE'
1066  include 'GEOM'
1067  include 'TSTEP'
1068 
1069 
1070  real du(mx*mx,nel)
1071  real u(mx*mx,nel)
1072  real c(md*md,nel,2)
1073  parameter(ldd=lxd*lyd*lzd)
1074  common /ctmp1/ ur(ldd),us(ldd),ut(ldd),ju(ldd),ud(ldd),tu(ldd)
1075  real ju
1076 
1077  logical if3d,ifd
1078  integer e
1079 
1080  if3d = .false.
1081  ifd = .false.
1082  if (md.ne.mx) ifd=.true.
1083 
1084  nrstd = md**2
1085  call lim_chk(nrstd,ldd,'urus ','ldd ','convp_cons')
1086 
1087  if (nio.eq.0.and.istep.lt.3) write(6,*) 'convp_cons',istep
1088 
1089  do e=1,nel
1090 
1091  call intp_rstd (ju,u(1,e),mx,md,if3d,0) ! 0 = forward; on Gauss points!
1092  call rzero (ud,nrstd)
1093 
1094 c call outmat(c(1,e,1),md,md,'fine u',e)
1095 c call outmat(c(1,e,2),md,md,'fine v',e)
1096 c call outmat(ju ,md,md,'fine T',e)
1097 
1098  do j=1,ldim
1099  do i=1,nrstd
1100  tu(i)=c(i,e,j)*ju(i) ! C_j*T
1101  enddo
1102  call grad_rst(ur,us,ut,tu,md,if3d) ! Already on fine (Gauss) mesh
1103 
1104  j0 = j+0
1105  j2 = j+2
1106  do i=1,nrstd ! rx has mass matrix and Jacobian on fine mesh
1107  ud(i)=ud(i)+rx(i,j0,e)*ur(i)+rx(i,j2,e)*us(i)
1108  enddo
1109  enddo
1110 
1111  call intp_rstd(du(1,e),ud,mx,md,if3d,1) ! 1 --> transpose
1112 
1113  enddo
1114 
1115 c call exitti('convop_cons_2d$',istep)
1116 
1117  return
1118  end
1119 c-----------------------------------------------------------------------
1120  subroutine iface_vert_int8(fa,va,jz0,jz1,nel)
1121  include 'SIZE'
1122  integer*8 fa(lx1*lz1,2*ldim,nel),va(0:lx1+1,0:ly1+1,jz0:jz1,nel)
1123  integer e,f
1124 
1125  n = lx1*lz1*2*ldim*nel
1126  call i8zero(fa,n)
1127 
1128  mx1 = lx1+2
1129  my1 = ly1+2
1130  mz1 = lz1+2
1131  if (ldim.eq.2) mz1=1
1132 
1133  nface = 2*ldim
1134  do e=1,nel
1135  do f=1,nface
1136  call facind (kx1,kx2,ky1,ky2,kz1,kz2,lx1,ly1,lz1,f)
1137 
1138  if (f.eq.1) then ! EB notation
1139  ky1=ky1-1
1140  ky2=ky1
1141  elseif (f.eq.2) then
1142  kx1=kx1+1
1143  kx2=kx1
1144  elseif (f.eq.3) then
1145  ky1=ky1+1
1146  ky2=ky1
1147  elseif (f.eq.4) then
1148  kx1=kx1-1
1149  kx2=kx1
1150  elseif (f.eq.5) then
1151  kz1=kz1-1
1152  kz2=kz1
1153  elseif (f.eq.6) then
1154  kz1=kz1+1
1155  kz2=kz1
1156  endif
1157 
1158  i = 0
1159  do iz=kz1,kz2
1160  do iy=ky1,ky2
1161  do ix=kx1,kx2
1162  i = i+1
1163  fa(i,f,e)=va(ix,iy,iz,e)
1164 c write(6,*) 'fa:',fa(i,f,e),i,f,e
1165  enddo
1166  enddo
1167  enddo
1168  enddo
1169  enddo
1170 
1171  return
1172  end
1173 c-----------------------------------------------------------------------
1174  subroutine set_binv(bmnv,hmsk,n) ! Store binvm1*(hyperbolic mask)
1175 
1176  include 'SIZE'
1177  include 'PARALLEL'
1178  include 'MASS'
1179 
1180  real bmnv(n,lorder),hmsk(n)
1181 
1182  do i=lorder,2,-1
1183  call copy(bmnv(1,i),bmnv(1,i-1),n)
1184  enddo
1185 
1186  call copy (bmnv,bm1,n) ! Fill bmnv(1,1)
1187 
1188  call fgslib_gs_op(gsh_fld(1),bmnv,1,1,0) ! 1 ==> +; gsh_fld(1) is velocity
1189 
1190  do i=1,n
1191  bmnv(i,1)=hmsk(i)/bmnv(i,1)
1192  enddo
1193 
1194  return
1195  end
1196 c-----------------------------------------------------------------------
1197  subroutine set_bdivw(bdivw,hmsk,n) ! Store binvm1*(hyperbolic mask)
1198 
1199  include 'SIZE'
1200  include 'PARALLEL'
1201  include 'MASS'
1202  include 'INPUT'
1203  include 'MVGEOM'
1204  common /scruz/ cx(lx1*ly1*lz1*lelt)
1205  $ , cy(lx1*ly1*lz1*lelt)
1206  $ , cz(lx1*ly1*lz1*lelt)
1207 
1208  real bdivw(n,lorder),hmsk(n)
1209 
1210  do i=lorder,2,-1
1211  call copy(bdivw(1,i),bdivw(1,i-1),n)
1212  enddo
1213 
1214  call gradm1 (bdivw,cy ,cz , wx )
1215  call gradm1 (cx ,cy ,cz , wy )
1216  call add2 (bdivw,cy , n )
1217  if (if3d) then
1218  call gradm1 (cx ,cy ,cz , wz )
1219  call add2 (bdivw,cz , n )
1220  endif
1221  call col2(bdivw,bm1,n)
1222 
1223  return
1224  end
1225 c-----------------------------------------------------------------------
1226  subroutine set_bmass(bmass,hmsk,n) ! Store bmass*(hyperbolic mask)
1227 
1228  include 'SIZE'
1229  include 'PARALLEL'
1230  include 'MASS'
1231 
1232  real bmass(n,lorder),hmsk(n)
1233 
1234  do i=lorder,2,-1
1235  call copy(bmass(1,i),bmass(1,i-1),n)
1236  enddo
1237 
1238  call copy (bmass,bm1,n) ! Fill bmass(1,1)
1239 
1240  return
1241  end
1242 c-----------------------------------------------------------------------
1243  subroutine setup_dg_gs(dgh,nx,ny,nz,nel,melg,vertex)
1244 
1245 c Global-to-local mapping for gs
1246 
1247  include 'SIZE'
1248  include 'TOTAL'
1249 
1250  integer dgh,vertex(1)
1251 
1252  parameter(lf=lx1*lz1*2*ldim*lelt)
1253  common /c_is1/ glo_num_face(lf)
1254  $ , glo_num_vol((lx1+2)*(ly1+2)*(lz1+2)*lelt)
1255  integer*8 glo_num_face,glo_num_vol,ngv
1256 
1257  common /nekmpi/ mid,mp,nekcomm,nekgroup,nekreal
1258 
1259  mx = nx+2
1260  call set_vert(glo_num_vol,ngv,mx,nel,vertex,.false.)
1261 
1262  mz0 = 1
1263  mz1 = 1
1264  if (if3d) mz0 = 0
1265  if (if3d) mz1 = lz1+1
1266  call iface_vert_int8 (glo_num_face,glo_num_vol,mz0,mz1,nelt)
1267 
1268  nf = lx1*lz1*2*ldim*nelt !total number of points on faces
1269  call fgslib_gs_setup(dgh,glo_num_face,nf,nekcomm,np)
1270 
1271  return
1272  end
1273 c-----------------------------------------------------------------------
1274  subroutine dg_set_fc_ptr
1275 
1276 c Set up pointer to restrict u to faces ! NOTE: compact
1277 
1278  include 'SIZE'
1279  include 'TOTAL'
1280 
1281  integer e,f,ef
1282 
1283  call dsset(lx1,ly1,lz1) ! set skpdat
1284 
1285  nxyz = lx1*ly1*lz1
1286  nxz = lx1*lz1
1287  nface = 2*ldim
1288  nxzf = lx1*lz1*nface ! red'd mod to area, unx, etc.
1289 
1290  k = 0
1291 
1292  do e=1,nelv
1293  do ef=1,nface ! EB notation
1294 
1295  f = eface1(ef)
1296  js1 = skpdat(1,f)
1297  jf1 = skpdat(2,f)
1298  jskip1 = skpdat(3,f)
1299  js2 = skpdat(4,f)
1300  jf2 = skpdat(5,f)
1301  jskip2 = skpdat(6,f)
1302 
1303  i = 0
1304  do j2=js2,jf2,jskip2
1305  do j1=js1,jf1,jskip1
1306 
1307  i = i+1
1308  k = i+nxz*(ef-1)+nxzf*(e-1) ! face numbering
1309  dg_face(k) = j1+lx1*(j2-1)+nxyz*(e-1) ! global numbering
1310 
1311  enddo
1312  enddo
1313 
1314  enddo
1315  enddo
1316  ndg_facex = nxzf*nelv
1317 
1318  return
1319  end
1320 c-----------------------------------------------------------------------
1321  subroutine full2face(faceary, vol_ary)
1322 
1323  include 'SIZE'
1324  include 'TOTAL'
1325 
1326  real faceary(lx1*lz1,2*ldim,lelt)
1327  real vol_ary(lx1,ly1,lz1,lelt)
1328  integer i,j
1329 
1330  do j=1,ndg_facex
1331  i=dg_face(j)
1332  faceary(j,1,1) = vol_ary(i,1,1,1)
1333  enddo
1334 
1335 
1336  return
1337  end
1338 c-----------------------------------------------------------------------
1339  subroutine face2full(vol_ary, faceary)
1340 
1341  include 'SIZE'
1342  include 'TOTAL'
1343 
1344  real faceary(lx1*lz1,2*ldim,lelt)
1345  real vol_ary(lx1,ly1,lz1,lelt)
1346  integer i,j
1347 
1348  n=lx1*ly1*lz1*nelfld(ifield)
1349  call rzero(vol_ary,n)
1350 
1351  do j=1,ndg_facex
1352  i=dg_face(j)
1353  vol_ary(i,1,1,1) = vol_ary(i,1,1,1)+faceary(j,1,1)
1354  enddo
1355 
1356  return
1357  end
1358 c-----------------------------------------------------------------------
1359  subroutine add_face2full(vol_ary, faceary)
1360 
1361  include 'SIZE'
1362  include 'TOTAL'
1363 
1364  real faceary(lx1*lz1,2*ldim,lelt)
1365  real vol_ary(lx1,ly1,lz1,lelt)
1366  integer i,j
1367 
1368  do j=1,ndg_facex
1369  i=dg_face(j)
1370  vol_ary(i,1,1,1) = vol_ary(i,1,1,1)+faceary(j,1,1)
1371  enddo
1372 
1373  return
1374  end
1375 c-----------------------------------------------------------------------
1376  subroutine conv_rhs_dg_aliased (du,u,c)
1377 c
1378  include 'SIZE'
1379  include 'TOTAL'
1380 
1381 c Apply convecting field c(1,ldim) to scalar field u(1).
1382 
1383  real du(1),u(1),c(1)
1384 
1385  parameter(lf=lx1*lz1*2*ldim*lelt)
1386  common /scrdg/ uf(lf),uxf(lf),uyf(lf),uzf(lf),upwind_wgt(lf)
1387 
1388  integer e,f
1389 
1390  n = lx1*ly1*lz1*nelv
1391  nf = lx1*lz1*2*ldim*nelt
1392 
1393 
1394  if (ifcons) then
1395  if (if3d ) call convop_cons_3d (du,u,c,lx1,lxd,nelv)
1396  if (.not.if3d) call convop_cons_2d (du,u,c,lx1,lxd,nelv)
1397  else
1398  if (if3d ) call convop_fst_3d (du,u,c,lx1,lxd,nelv)
1399  if (.not.if3d) call convop_fst_2d (du,u,c,lx1,lxd,nelv)
1400  endif
1401 
1402  call full2face(uf ,u )
1403  call full2face(uxf,vx)
1404  call full2face(uyf,vy)
1405  call full2face(uzf,vz)
1406  if (.not.if3d) call rzero(uzf,nf)
1407 
1408  beta_u = 0.00 ! 1=full upwind; 0=central flux
1409  beta_u = 0.25 ! 1=full upwind; 0=central flux
1410  beta_u = 1.00 ! 1=full upwind; 0=central flux
1411  beta_u = param(98)
1412  if (istep.le.5.and.nio.eq.0) write(6,*) beta_u,' dg upwind'
1413 
1414  nface = 2*ldim
1415  nxz = lx1*lz1
1416  k = 0
1417  do e=1,nelt
1418  do f=1,nface
1419  do i=1,nxz
1420 
1421  k=k+1
1422 
1423  beta = ( unx(i,1,f,e)*uxf(k)
1424  $ + uny(i,1,f,e)*uyf(k)
1425  $ + unz(i,1,f,e)*uzf(k))
1426 
1427  uf(k) = -beta*area(i,1,f,e)*uf(k)
1428 
1429  upwind_wgt(k) = 1.0
1430  if (beta.gt.0) upwind_wgt(k) = 0.0
1431  upwind_wgt(k) = 0.5*(1-beta_u) + upwind_wgt(k)*beta_u
1432  if (beta.eq.0) upwind_wgt(k) = 0.5
1433 
1434  enddo
1435  enddo
1436  enddo
1437 
1438  call fgslib_gs_op(dg_hndlx,uf,1,1,0) ! 1 ==> +
1439  call col2 (uf,upwind_wgt,nf) ! Inefficient, but ok for now.
1440  ! Should be combined with
1441  call add_face2full( du , uf) ! <--- this stmt.
1442 
1443  call fbinvert ( du ) ! Right now works only for undeformed
1444 
1445  return
1446  end
1447 c-----------------------------------------------------------------------
1448  subroutine map_faced(ju,u,mx,md,fdim,idir) ! GLL->GL interpolation
1449 
1450 c GLL interpolation from mx to md for a face array of size (nx,nz)
1451 
1452 c If idir ^= 0, then apply transpose operator (md to mx)
1453 
1454  include 'SIZE'
1455 
1456  real ju(1),u(1)
1457  integer fdim
1458 
1459  parameter(ldg=lxd**3,lwkd=4*lxd*lxd)
1460  common /dgrad/ d(ldg),dt(ldg),dg(ldg),dgt(ldg),jgl(ldg),jgt(ldg)
1461  $ , wkd(lwkd)
1462  real jgl,jgt
1463 
1464  parameter(ld=2*lxd)
1465  common /ctmp0/ w(ld**ldim,2)
1466 
1467  call lim_chk(md,ld,'md ','ld ','map_faced ')
1468  call lim_chk(mx,ld,'mx ','ld ','map_faced ')
1469 
1470 c call copy(ju,u,mx)
1471 c return
1472 
1473  call get_int_ptr (i,mx,md)
1474 
1475  if (idir.eq.0) then
1476  if (fdim.eq.2) then
1477  call mxm(jgl(i),md,u,mx,wkd,mx)
1478  call mxm(wkd,md,jgt(i),mx,ju,md)
1479  else
1480  call mxm(jgl(i),md,u,mx,ju,1)
1481  endif
1482  else
1483  if (fdim.eq.2) then
1484  call mxm(jgt(i),mx,u,md,wkd,md)
1485  call mxm(wkd,mx,jgl(i),md,ju,mx)
1486  else
1487  call mxm(jgt(i),mx,u,md,ju,1)
1488  endif
1489  endif
1490 
1491  return
1492  end
1493 c-----------------------------------------------------------------------
1494  subroutine fbinvert(rhs) ! Still in development. 10/10/15, pff.
1495 
1496  include 'SIZE'
1497  include 'TOTAL'
1498 
1499  real rhs(lx1,ly1,lz1,lelt)
1500 
1501  common /cfbinv/ qn(lx1),alpha_n,beta_n
1502  $ ,s1(ly1,lz1),bnv(lx1)
1503  $ ,tmp(lx1*ly1*lz1*lelt)
1504  integer icalld
1505  save icalld
1506  data icalld /0/
1507 
1508  integer e
1509 
1510  n = lx1*ly1*lz1*nelfld(ifield)
1511 
1512  do i=1,n ! FOR NOW, USE DIAGONAL MASS
1513  rhs(i,1,1,1)=rhs(i,1,1,1)/bm1(i,1,1,1) ! MATRIX. pff, 10/10/15
1514  enddo
1515 
1516  return
1517  end
1518 c-----------------------------------------------------------------------
1519  subroutine grad_rstd_ta(du,ur,us,ut,md,if3d) ! GL->GL gradt
1520 
1521  include 'SIZE'
1522  include 'DXYZ'
1523 
1524  real ur(1),us(1),ut(1),u(1)
1525 
1526  logical if3d
1527 
1528  parameter(ldg=lxd**3,lwkd=4*lxd*lxd)
1529  common /dgrad/ d(ldg),dt(ldg),dg(ldg),dgt(ldg),jgl(ldg),jgt(ldg)
1530  $ , wkd(lwkd)
1531  real jgl,jgt
1532 
1533  call get_dgl_ptr (ip,md,md)
1534  call gradrta (du,ur,us,ut,dgt(ip),dg(ip),dg(ip),md,md,md,if3d)
1535 
1536  return
1537  end
1538 c-----------------------------------------------------------------------
1539  subroutine set_dg_wgts
1540 
1541  include 'SIZE'
1542  include 'TOTAL'
1543 
1544  integer icalld
1545  save icalld
1546  data icalld /0/
1547 
1548  common /finewts/ zptf(lxd),wgtf(lxd),wghtf(lxd*lzd),wghtc(lx1*lz1)
1549 
1550  if (icalld.eq.0) then ! Set fine-scale surface weights
1551 
1552  icalld = 1
1553 
1554  call zwgl(zptf,wgtf,lxd)
1555  if (if3d) then
1556  k=0
1557  do j=1,ly1
1558  do i=1,lx1
1559  k=k+1
1560  wghtc(k)=wxm1(i)*wzm1(j)
1561  enddo
1562  enddo
1563  k=0
1564  do j=1,lyd
1565  do i=1,lxd
1566  k=k+1
1567  wghtf(k)=wgtf(i)*wgtf(j)
1568  enddo
1569  enddo
1570  else
1571  call copy(wghtc,wxm1,lx1)
1572  call copy(wghtf,wgtf,lxd)
1573  endif
1574  endif
1575 
1576  return
1577  end
1578 c-----------------------------------------------------------------------
1579  subroutine conv_rhs_dg (du,u,c)
1580 c
1581  include 'SIZE'
1582  include 'TOTAL'
1583 
1584 c Apply convecting field c(1,ldim) to scalar field u(1).
1585 
1586  parameter(ldd=lxd*lyd*lzd)
1587  real du(1),u(1),c(ldd*lelv,3)
1588 
1589  parameter(lf=lx1*lz1*2*ldim*lelt)
1590  common /scrdg/ uf(lf),uxf(lf),uyf(lf),uzf(lf),upwind_wgt(lf)
1591  $ , beta_c(lx1*lz1),jaco_c(lx1*lz1)
1592  $ , beta_f(lxd*lzd),jaco_f(lxd*lzd)
1593  $ , ufine(lxd*lzd)
1594  real jaco_c,jaco_f
1595  common /finewts/ zptf(lxd),wgtf(lxd),wghtf(lxd*lzd),wghtc(lx1*lz1)
1596 
1597  integer e,f,fdim
1598 
1599 
1600  n = lx1*ly1*lz1*lelv
1601  nf = lx1*lz1*2*ldim*lelt
1602 
1603  call conv_rhs_dg_weak (du,u,c(1,1),c(1,2),c(1,3))
1604 
1605  call full2face(uf ,u )
1606  call full2face(uxf,vx)
1607  call full2face(uyf,vy)
1608  call full2face(uzf,vz)
1609  if (.not.if3d) call rzero(uzf,nf)
1610 
1611  beta_u = 0.00 ! 1=full upwind; 0=central flux
1612  beta_u = 0.25 ! 1=full upwind; 0=central flux
1613  beta_u = 1.00 ! 1=full upwind; 0=central flux
1614  beta_u = param(98)
1615  if (istep.le.5.and.nio.eq.0) write(6,*) beta_u,' dg upwind'
1616 
1617  nface = 2*ldim
1618  nxz = lx1*lz1
1619  nxzd = lxd*lzd
1620  k = 0
1621  do e=1,nelt ! This formula for upwind weights appears to
1622  do f=1,nface ! assume that U is continuous, or at least of
1623 
1624  kface = k+1
1625  do i=1,nxz ! the same sign on either side of the interface.
1626 
1627  k=k+1
1628 
1629  beta = ( unx(i,1,f,e)*uxf(k)
1630  $ + uny(i,1,f,e)*uyf(k)
1631  $ + unz(i,1,f,e)*uzf(k))
1632 
1633  upwind_wgt(k) = 1.0
1634  if (beta.gt.0) upwind_wgt(k) = 0.0
1635  upwind_wgt(k) = 0.5*(1-beta_u) + upwind_wgt(k)*beta_u
1636 
1637  beta_c(i) = -beta
1638  jaco_c(i) = area(i,1,f,e)/wghtc(i)
1639 
1640  enddo
1641 
1642  fdim = ldim-1 ! Dimension of face
1643  call map_faced(beta_f,beta_c ,lx1,lxd,fdim,0) ! Dealiased quadrature,
1644  call map_faced(jaco_f,jaco_c ,lx1,lxd,fdim,0) ! 0 --> coarse to fine,
1645  call map_faced(ufine,uf(kface),lx1,lxd,fdim,0) ! ufine = J uf
1646 
1647  do i=1,nxzd
1648  ufine(i)=wghtf(i)*jaco_f(i)*beta_f(i)*ufine(i)
1649  enddo
1650  call map_faced(uf(kface),ufine,lx1,lxd,fdim,1) ! 1 --> uf = J^T ufine
1651 
1652  enddo
1653  enddo
1654 
1655  call fgslib_gs_op(dg_hndlx,uf,1,1,0) ! 1 ==> +
1656 
1657  call col2 (uf,upwind_wgt,nf) ! Inefficient, but ok for now.
1658  ! Should be combined with
1659  call add_face2full( du , uf) ! <--- this stmt.
1660 
1661  call fbinvert ( du ) ! Right now works only for undeformed
1662 
1663  return
1664  end
1665 c-----------------------------------------------------------------------
1666  subroutine convect_dg(du,u,ifuf,cr,cs,ct,ifcf)
1667  include 'SIZE'
1668  include 'TOTAL'
1669  real du(1),u(1),cr(1),cs(1),ct(1)
1670  logical ifuf,ifcf
1671 
1672  call conv_rhs_dg_weak (du,u,cr,cs,ct)
1673 
1674  return
1675  end
1676 c-----------------------------------------------------------------------
1677  subroutine convop_weak(du,u,cr,cs,ct,mx,md,nel) ! Weak Conservation form
1678 
1679 c Apply convecting field c to scalar field u, conservation form d/dxj cj phi
1680 
1681 c Assumes that current convecting field is on dealias mesh, in c()
1682 
1683  include 'SIZE'
1684  include 'TOTAL'
1685 
1686  parameter(lxx=lx1*ly1*lz1,ldd=lxd*lyd*lzd)
1687  real du(lxx,nel)
1688  real u(lxx,nel)
1689  real cr(ldd,nel),cs(ldd,nel),ct(ldd,nel)
1690  common /ctmp1/ ur(ldd),us(ldd),ut(ldd),ju(ldd),ud(ldd),tu(ldd)
1691  real ju
1692 
1693  integer e
1694 
1695  nxyz = lx1*ly1*lz1
1696  nrstd = md**ldim
1697 
1698  call lim_chk(nrstd,ldd,'urus5','ldd ','convp_cons')
1699 
1700  do e=1,nel
1701 
1702  call intp_rstd (ju,u(1,e),mx,md,if3d,0) ! 0 = forward; on Gauss points!
1703  if (ldim.eq.3) then
1704  do i=1,ldd
1705  ur(i)=ju(i)*cr(i,e) ! Already in r-s-t coordinates
1706  us(i)=ju(i)*cs(i,e)
1707  ut(i)=ju(i)*ct(i,e)
1708  ud(i)=0
1709  enddo
1710  else
1711  do i=1,ldd
1712  ur(i)=ju(i)*cr(i,e) ! Already in r-s-t coordinates
1713  us(i)=ju(i)*cs(i,e)
1714  ud(i)=0
1715  enddo
1716  endif
1717  call grad_rstd_ta (ud,ur,us,ut,lxd,if3d) ! GL->GL gradt
1718  call intp_rstd (du(1,e),ud,mx,md,if3d,1) ! 1 = backward; on Gauss points!
1719  do i=1,nxyz
1720  du(i,e) = -du(i,e)*binvdg(i,e)
1721  enddo
1722 
1723  enddo
1724 
1725  return
1726  end
1727 c-----------------------------------------------------------------------
1728  subroutine conv_bdry_dg_weak (du,u) ! THIS SHOULD HAVE: ,cr,cs,ct)
1729 c
1730 c Implement Cu = Div (cu) in weak form using DG
1731 c
1732  include 'SIZE'
1733  include 'TOTAL'
1734 
1735 c Apply convecting field c(1,ldim) to scalar field u(1).
1736 
1737  real du(1),u(1)
1738 
1739  parameter(lf=lx1*lz1*2*ldim*lelt)
1740  common /scrdg/uf(lf),uxf(lf),uyf(lf),uzf(lf),upwind_wgt(lf),us(lf)
1741  $ ,beta_c(lx1*lz1),jaco_c(lx1*lz1)
1742  $ ,beta_f(lxd*lzd),jaco_f(lxd*lzd)
1743  $ ,ufine(lxd*lzd)
1744  real jaco_c,jaco_f
1745 
1746  common /finewts/ zptf(lxd),wgtf(lxd),wghtf(lxd*lzd),wghtc(lx1*lz1)
1747 
1748  integer e,f,fdim
1749 
1750  n = lx1*ly1*lz1*nelv
1751  nf = lx1*lz1*2*ldim*nelt
1752 
1753  call full2face(uf ,u )
1754  call full2face(uxf,vx)
1755  call full2face(uyf,vy)
1756  call full2face(uzf,vz)
1757  if (.not.if3d) call rzero(uzf,nf)
1758 
1759  beta_u = 1.00 ! 1=full upwind; 0=central flux
1760 
1761  nface = 2*ldim
1762  nxz = lx1*lz1
1763  nxzd = lxd*lzd
1764  k = 0
1765  do e=1,nelt ! This formula for upwind weights appears to
1766  do f=1,nface ! assume that U is continuous, or at least of
1767 
1768  kface = k+1
1769  if (fw(f,e).gt.0.6) then
1770 
1771  do i=1,nxz ! the same sign on either side of the interface.
1772 
1773  k=k+1
1774 
1775  beta = ( unx(i,1,f,e)*uxf(k)
1776  $ + uny(i,1,f,e)*uyf(k)
1777  $ + unz(i,1,f,e)*uzf(k))
1778 
1779  upwind_wgt(k) = 0.0
1780  if (beta.lt.0) upwind_wgt(k) = 1.0
1781 
1782  beta_c(i) = beta
1783  jaco_c(i) = area(i,1,f,e)/wghtc(i)
1784 
1785  enddo
1786 
1787  fdim = ldim-1 ! Dimension of face
1788  call map_faced(beta_f,beta_c ,lx1,lxd,fdim,0) ! Dealiased quadrature,
1789  call map_faced(jaco_f,jaco_c ,lx1,lxd,fdim,0) ! 0 --> coarse to fine,
1790  call map_faced(ufine,uf(kface),lx1,lxd,fdim,0) ! ufine = J uf
1791 
1792  do i=1,nxzd
1793  ufine(i)=wghtf(i)*jaco_f(i)*beta_f(i)*ufine(i)
1794  enddo
1795  call map_faced(uf(kface),ufine,lx1,lxd,fdim,1) ! 1 --> uf = J^T ufine
1796 
1797  else
1798  do i=1,nxz
1799  k=k+1
1800  upwind_wgt(k) = 0.0
1801  enddo
1802  endif
1803 
1804  enddo
1805  enddo
1806 
1807  do j=1,ndg_facex
1808  i=dg_face(j)
1809  du(i) = du(i) - ( upwind_wgt(j)*uf(j) )
1810  enddo
1811 
1812  return
1813  end
1814 c-----------------------------------------------------------------------
1815  subroutine conv_rhs_dg_weak (du,u,cr,cs,ct)
1816 c
1817 c Implement Cu = Div (cu) in weak form using DG
1818 c
1819  include 'SIZE'
1820  include 'TOTAL'
1821 
1822 c Apply convecting field c(1,ldim) to scalar field u(1).
1823 
1824  real du(1),u(1),cr(1),cs(1),ct(1)
1825 
1826  parameter(lf=lx1*lz1*2*ldim*lelt)
1827  common /scrdg/uf(lf),uxf(lf),uyf(lf),uzf(lf),upwind_wgt(lf),us(lf)
1828  $ ,beta_c(lx1*lz1),jaco_c(lx1*lz1)
1829  $ ,beta_f(lxd*lzd),jaco_f(lxd*lzd)
1830  $ ,ufine(lxd*lzd)
1831  real jaco_c,jaco_f
1832 
1833  common /finewts/ zptf(lxd),wgtf(lxd),wghtf(lxd*lzd),wghtc(lx1*lz1)
1834 
1835  integer e,f,fdim
1836 
1837  n = lx1*ly1*lz1*lelv
1838  nf = lx1*lz1*2*ldim*lelt
1839 
1840 
1841  call convop_weak (du,u,cr,cs,ct,lx1,lxd,nelv) ! Volumetric term
1842 
1843  call full2face(uf ,u )
1844  call full2face(uxf,vx)
1845  call full2face(uyf,vy)
1846  call full2face(uzf,vz)
1847  if (.not.if3d) call rzero(uzf,nf)
1848 
1849  beta_u = 0.25 ! 1=full upwind; 0=central flux
1850  beta_u = 0.00 ! 1=full upwind; 0=central flux
1851  beta_u = 1.00 ! 1=full upwind; 0=central flux
1852  if (istep.le.5.and.nio.eq.0) write(6,*) beta_u,' dg upwind'
1853 
1854  nface = 2*ldim
1855  nxz = lx1*lz1
1856  nxzd = lxd*lzd
1857  k = 0
1858  do e=1,nelt ! This formula for upwind weights appears to
1859  do f=1,nface ! assume that U is continuous, or at least of
1860 
1861  kface = k+1
1862  do i=1,nxz ! the same sign on either side of the interface.
1863  k=k+1
1864  beta = ( unx(i,1,f,e)*uxf(k)
1865  $ + uny(i,1,f,e)*uyf(k)
1866  $ + unz(i,1,f,e)*uzf(k) )
1867 
1868  upwind_wgt(k) = 0.0
1869  if (beta.gt.0) upwind_wgt(k) = 1.0
1870  upwind_wgt(k) = 0.5*(1-beta_u) + beta_u*(1-upwind_wgt(k))
1871 
1872  if (fw(f,e).gt.0.6 .and. beta.lt.0) upwind_wgt(k)=1.
1873  if (fw(f,e).gt.0.6 .and. beta.gt.0) upwind_wgt(k)=0.
1874 
1875  beta_c(i) = beta
1876  jaco_c(i) = area(i,1,f,e)/wghtc(i)
1877  enddo
1878 
1879  fdim = ldim-1 ! Dimension of face
1880  call map_faced(beta_f,beta_c ,lx1,lxd,fdim,0) ! Dealiased quadrature,
1881  call map_faced(jaco_f,jaco_c ,lx1,lxd,fdim,0) ! 0 --> coarse to fine,
1882  call map_faced(ufine,uf(kface),lx1,lxd,fdim,0) ! ufine = J uf
1883 
1884  do i=1,nxzd
1885  ufine(i)=wghtf(i)*jaco_f(i)*beta_f(i)*ufine(i)
1886  enddo
1887  call map_faced(uf(kface),ufine,lx1,lxd,fdim,1) ! 1 --> uf = J^T ufine
1888  call copy (us(kface),uf(kface),lx1*lz1) ! Save uf for later recombination
1889 
1890  enddo
1891  enddo
1892 
1893  call fgslib_gs_op(dg_hndlx,uf,1,1,0) ! 1 ==> + : uf <-- uf^- + uf^+
1894 
1895  do j=1,ndg_facex
1896  i=dg_face(j)
1897  du(i) = du(i) + ( us(j)-upwind_wgt(j)*uf(j) )*binvdg(i,1)
1898  enddo
1899 
1900  return
1901  end
1902 c-----------------------------------------------------------------------
exitt
void exitt()
Definition: comm_mpi.f:604
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
dsset
subroutine dsset(nx, ny, nz)
Definition: connect1.f:553
set_dealias_rx
subroutine set_dealias_rx
Definition: convect2.f:121
conv_rhs
subroutine conv_rhs(du, u, c, bmsk, bmst, bdwt, gsl)
Definition: convect.f:193
iface_vert_int8
subroutine iface_vert_int8(fa, va, jz0, jz1, nel)
Definition: convect.f:1121
set_dg_wgts
subroutine set_dg_wgts
Definition: convect.f:1540
convop_cons_3d
subroutine convop_cons_3d(du, u, c, mx, md, nel)
Definition: convect.f:1005
convch
subroutine convch
Definition: convect.f:965
set_conv_char
subroutine set_conv_char(ct, c, ux, uy, uz, nelc, tau, ifnew)
Definition: convect.f:544
char_conv1
subroutine char_conv1(p0, u, bmnv, n, ulag, ln, gsl, c, m, cs, nc, ct, u1, r1, r2, r3, r4, bmsk, bdivw, bdwt, bmass, bmst, bm, bmlag)
Definition: convect.f:46
set_binv
subroutine set_binv(bmnv, hmsk, n)
Definition: convect.f:1175
convop_fst_3d
subroutine convop_fst_3d(du, u, c, mx, md, nel)
Definition: convect.f:237
dg_set_fc_ptr
subroutine dg_set_fc_ptr
Definition: convect.f:1275
set_bdivw
subroutine set_bdivw(bdivw, hmsk, n)
Definition: convect.f:1198
conv_rhs_dg
subroutine conv_rhs_dg(du, u, c)
Definition: convect.f:1580
get_dgl_ptr
subroutine get_dgl_ptr(ip, mx, md)
Definition: convect.f:506
gen_int
subroutine gen_int(jgl, jgt, mp, np, w)
Definition: convect.f:382
int_vel
subroutine int_vel(c_t, t0, c, n, nc, ct, nid)
Definition: convect.f:161
set_convect_new
subroutine set_convect_new(cr, cs, ct, ux, uy, uz)
Definition: convect.f:846
conv_rhs_dg_weak
subroutine conv_rhs_dg_weak(du, u, cr, cs, ct)
Definition: convect.f:1816
full2face
subroutine full2face(faceary, vol_ary)
Definition: convect.f:1322
face2full
subroutine face2full(vol_ary, faceary)
Definition: convect.f:1340
get_int_ptr
subroutine get_int_ptr(ip, mx, md)
Definition: convect.f:468
intp_rstd
subroutine intp_rstd(ju, u, mx, md, if3d, idir)
Definition: convect.f:347
convop_fst_2d
subroutine convop_fst_2d(du, u, c, mx, md, nel)
Definition: convect.f:279
gen_dgl
subroutine gen_dgl(dgl, dgt, mp, np, w)
Definition: convect.f:415
grad_rstd
subroutine grad_rstd(ur, us, ut, u, mx, md, if3d, ju)
Definition: convect.f:321
convop_cons_2d
subroutine convop_cons_2d(du, u, c, mx, md, nel)
Definition: convect.f:1060
add_face2full
subroutine add_face2full(vol_ary, faceary)
Definition: convect.f:1360
fbinvert
subroutine fbinvert(rhs)
Definition: convect.f:1495
set_convect_cons
subroutine set_convect_cons(cx, cy, cz, ux, uy, uz)
Definition: convect.f:818
convect_dg
subroutine convect_dg(du, u, ifuf, cr, cs, ct, ifcf)
Definition: convect.f:1667
set_ct_cvx
subroutine set_ct_cvx(ct, c, m, u, v, w, tau, nc, mc, nelc, ifnew)
Definition: convect.f:575
convect_cons
subroutine convect_cons(bdu, u, ifuf, cx, cy, cz, ifcf)
Definition: convect.f:732
conv_rhs_dg_aliased
subroutine conv_rhs_dg_aliased(du, u, c)
Definition: convect.f:1377
set_bmass
subroutine set_bmass(bmass, hmsk, n)
Definition: convect.f:1227
setup_dg_gs
subroutine setup_dg_gs(dgh, nx, ny, nz, nel, melg, vertex)
Definition: convect.f:1244
map_faced
subroutine map_faced(ju, u, mx, md, fdim, idir)
Definition: convect.f:1449
char_conv
subroutine char_conv(p0, u, ulag, bm, bmlag, msk, c, cs, gsl)
Definition: convect.f:3
grad_rstd_ta
subroutine grad_rstd_ta(du, ur, us, ut, md, if3d)
Definition: convect.f:1520
convop_weak
subroutine convop_weak(du, u, cr, cs, ct, mx, md, nel)
Definition: convect.f:1678
lim_chk
subroutine lim_chk(n, m, avar5, lvar5, sub_name10)
Definition: convect.f:454
convect_new
subroutine convect_new(bdu, u, ifuf, cx, cy, cz, ifcf)
Definition: convect.f:639
grad_rst
subroutine grad_rst(ur, us, ut, u, md, if3d)
Definition: convect.f:615
advchar
subroutine advchar
Definition: convect.f:904
conv_bdry_dg_weak
subroutine conv_bdry_dg_weak(du, u)
Definition: convect.f:1729
fd_weights_full
subroutine fd_weights_full(xx, x, n, m, c)
Definition: fast3d.f:1293
gradrta
subroutine gradrta(u, ur, us, ut, Drt, Ds, Dt, nr, ns, nt, if3d)
Definition: hmholtz.f:1743
invcol2
subroutine invcol2(a, b, n)
Definition: math.f:73
transpose
subroutine transpose(a, lda, b, ldb)
Definition: math.f:1594
add3s12
subroutine add3s12(x, y, z, c1, c2, n)
Definition: math.f:1534
i8zero
subroutine i8zero(a, n)
Definition: math.f:939
add2
subroutine add2(a, b, n)
Definition: math.f:622
col2
subroutine col2(a, b, n)
Definition: math.f:564
copy
subroutine copy(a, b, n)
Definition: math.f:260
vlsum
real function vlsum(vec, n)
Definition: math.f:417
subcol3
subroutine subcol3(a, b, c, n)
Definition: math.f:186
cfill
subroutine cfill(a, b, n)
Definition: math.f:244
invcol3
subroutine invcol3(a, b, c, n)
Definition: math.f:98
rzero
subroutine rzero(a, n)
Definition: math.f:208
mxm
subroutine mxm(a, n1, b, n2, c, n3)
Definition: mxm_wrapper.f:2
local_grad2
subroutine local_grad2(ur, us, u, N, e, D, Dt)
Definition: navier5.f:448
local_grad3
subroutine local_grad3(ur, us, ut, u, N, e, D, Dt)
Definition: navier5.f:429
gradm1
subroutine gradm1(ux, uy, uz, u)
Definition: navier5.f:463
set_vert
subroutine set_vert(glo_num, ngv, nx, nel, vertex, ifcenter)
Definition: navier8.f:4
specmpn
subroutine specmpn(b, nb, a, na, ba, ab, if3d, w, ldw)
Definition: navier8.f:941
zwgl
subroutine zwgl(Z, W, NP)
Definition: speclib.f:91
zwgll
subroutine zwgll(Z, W, NP)
Definition: speclib.f:108