KTH framework for Nek5000 toolboxes; testing version  0.0.1
stat.f
Go to the documentation of this file.
1 
8 !=======================================================================
12  subroutine stat_register()
13  implicit none
14 
15  include 'SIZE'
16  include 'INPUT'
17  include 'FRAMELP'
18  include 'STATD'
19 
20  ! local variables
21  integer lpmid, il
22  real ltim
23  character*2 str
24 
25  ! functions
26  real dnekclock
27 !-----------------------------------------------------------------------
28  ! timing
29  ltim = dnekclock()
30 
31  ! check if the current module was already registered
32  call mntr_mod_is_name_reg(lpmid,stat_name)
33  if (lpmid.gt.0) then
34  call mntr_warn(lpmid,
35  $ 'module ['//trim(stat_name)//'] already registered')
36  return
37  endif
38 
39  ! find parent module
40  call mntr_mod_is_name_reg(lpmid,'FRAME')
41  if (lpmid.le.0) then
42  lpmid = 1
43  call mntr_abort(lpmid,
44  $ 'parent module ['//'FRAME'//'] not registered')
45  endif
46 
47  ! register module
48  call mntr_mod_reg(stat_id,lpmid,stat_name,
49  $ '2D and 3D statistics')
50 
51  ! check if 2D mapping module is registered
52  if (stat_rdim.eq.1) then
53  call mntr_mod_is_name_reg(lpmid,'MAP2D')
54  ! if not, register module
55  if (lpmid.le.0) then
56  call map2d_register()
57  endif
58  endif
59 
60  ! register timers
61  call mntr_tmr_is_name_reg(lpmid,'FRM_TOT')
62  ! total time
63  call mntr_tmr_reg(stat_tmr_tot_id,lpmid,stat_id,
64  $ 'STAT_TOT','Statistics total time',.false.)
65  lpmid = stat_tmr_tot_id
66  ! initialisation
67  call mntr_tmr_reg(stat_tmr_ini_id,lpmid,stat_id,
68  $ 'STAT_INI','Statistics initialisation time',.true.)
69  ! averagign
70  call mntr_tmr_reg(stat_tmr_avg_id,lpmid,stat_id,
71  $ 'STAT_AVG','Statistics averaging time',.true.)
72 
73  if (stat_rdim.eq.1) then
74  ! communication
75  call mntr_tmr_reg(stat_tmr_cmm_id,lpmid,stat_id,
76  $ 'STAT_CMM','Statistics communication time',.true.)
77  endif
78  ! IO
79  call mntr_tmr_reg(stat_tmr_io_id,lpmid,stat_id,
80  $ 'STAT_IO','Statistics IO time',.true.)
81 
82  ! register and set active section
83  call rprm_sec_reg(stat_sec_id,stat_id,'_'//adjustl(stat_name),
84  $ 'Runtime paramere section for statistics module')
85  call rprm_sec_set_act(.true.,stat_sec_id)
86 
87  ! register parameters
88  call rprm_rp_reg(stat_avstep_id,stat_sec_id,'AVSTEP',
89  $ 'Frequency of averaging',rpar_int,10,0.0,.false.,' ')
90 
91  call rprm_rp_reg(stat_skstep_id,stat_sec_id,'SKSTEP',
92  $ 'Skipped initial steps',rpar_int,10,0.0,.false.,' ')
93 
94  call rprm_rp_reg(stat_iostep_id,stat_sec_id,'IOSTEP',
95  $ 'Frequency of filed saving',rpar_int,100,0.0,.false.,' ')
96 
97  ! set initialisation flag
98  stat_ifinit=.false.
99 
100  ! timing
101  ltim = dnekclock() - ltim
102  call mntr_tmr_add(stat_tmr_tot_id,1,ltim)
103 
104  return
105  end subroutine
106 !=======================================================================
110  subroutine stat_init()
111  implicit none
112 
113  include 'SIZE'
114  include 'FRAMELP'
115  include 'TSTEP'
116  include 'MAP2D'
117  include 'STATD'
118 
119  ! local variables
120  integer itmp, il
121  real rtmp, ltim
122  logical ltmp
123  character*20 ctmp
124 
125  ! functions
126  real dnekclock
127 !-----------------------------------------------------------------------
128  call mntr_log(stat_id,lp_inf,'Initialisation started')
129 
130  ! check if the module was already initialised
131  if (stat_ifinit) then
132  call mntr_warn(stat_id,
133  $ 'module ['//trim(stat_name)//'] already initiaised.')
134  return
135  endif
136 
137  ! check if map2d module is initialised
138  if (stat_rdim.eq.1) then
139  if (.not.map2d_ifinit) call map2d_init()
140  endif
141  call mntr_log(stat_id,lp_inf,'Map 2D finalised')
142 
143  ! timing
144  ltim = dnekclock()
145 
146  ! get runtime parameters
147  call rprm_rp_get(itmp,rtmp,ltmp,ctmp,stat_avstep_id,rpar_int)
148  stat_avstep = itmp
149  call rprm_rp_get(itmp,rtmp,ltmp,ctmp,stat_skstep_id,rpar_int)
150  stat_skstep = itmp
151  call rprm_rp_get(itmp,rtmp,ltmp,ctmp,stat_iostep_id,rpar_int)
152  stat_iostep = itmp
153 
154  ! initialise time averaging variables for a given statistics file
155  stat_atime = 0.0
156  stat_tstart = time
157 
158  if (stat_rdim.eq.1) then
159  ! finish 3D => 2D mapping operations
160  ! get local integration coefficients
161  call mntr_log(stat_id,lp_vrb,'Getting loacal int. coeff.')
162  call stat_init_int1d
163  endif
164  call mntr_log(stat_id,lp_inf,'1D local finalised')
165 
166  ! reset statistics variables
167  itmp = lx1**(ldim-stat_rdim)*lelt*stat_lvar
168  call rzero(stat_ruavg,itmp)
169 
170  ! everything is initialised
171  stat_ifinit=.true.
172 
173  call mntr_log(stat_id,lp_inf,'Initialisation finalised')
174 
175  ! timing
176  ltim = dnekclock() - ltim
177  call mntr_tmr_add(stat_tmr_ini_id,1,ltim)
178 
179  return
180  end subroutine
181 !=======================================================================
185  subroutine stat_end()
186  implicit none
187 
188  include 'SIZE'
189  include 'FRAMELP'
190  include 'STATD'
191 
192  ! local variables
193  real ltim
194 
195  ! functions
196  real dnekclock
197 !-----------------------------------------------------------------------
198  ! make sure all data in the buffer is saved
199  if (stat_atime.gt.0.0) then
200  if (stat_rdim.eq.1) then
201  ltim = dnekclock()
202  call mntr_log(stat_id,lp_inf,'Global sum')
203  call stat_gs_sum
204  ltim = dnekclock() - ltim
205  call mntr_tmr_add(stat_tmr_cmm_id,1,ltim)
206  endif
207 
208  ltim = dnekclock()
209  call mntr_log(stat_id,lp_inf,'Writing stat file')
210  call stat_mfo
211  ltim = dnekclock() - ltim
212  call mntr_tmr_add(stat_tmr_io_id,1,ltim)
213  endif
214 
215  return
216  end subroutine
217 !=======================================================================
221  logical function stat_is_initialised()
222  implicit none
223 
224  include 'SIZE'
225  include 'STATD'
226 !-----------------------------------------------------------------------
227  stat_is_initialised = stat_ifinit
228 
229  return
230  end function
231 !=======================================================================
236  subroutine stat_avg
237  implicit none
238 
239  include 'SIZE'
240  include 'FRAMELP'
241  include 'TSTEP'
242  include 'STATD'
243 
244  ! local variables
245  integer itmp
246 
247  ! simple timing
248  real ltim
249 
250  ! functions
251  real dnekclock
252 !-----------------------------------------------------------------------
253  ! skip initial steps
254  if (istep.gt.stat_skstep) then
255  itmp = istep - stat_skstep
256 
257  ! average
258  if (mod(itmp,stat_avstep).eq.0) then
259  ! simple timing
260  ltim = dnekclock()
261  call mntr_log(stat_id,lp_inf,'Average compute')
262  call stat_compute
263  ! timing
264  ltim = dnekclock() - ltim
265  call mntr_tmr_add(stat_tmr_avg_id,1,ltim)
266  endif
267 
268  ! save statistics file and restart statistics variables
269  if (mod(itmp,stat_iostep).eq.0) then
270  if (stat_rdim.eq.1) then
271 
272  ltim = dnekclock()
273  call mntr_log(stat_id,lp_inf,'Global sum')
274  call stat_gs_sum
275  ltim = dnekclock() - ltim
276  call mntr_tmr_add(stat_tmr_cmm_id,1,ltim)
277  endif
278 
279  ltim = dnekclock()
280  call mntr_log(stat_id,lp_inf,'Writing stat file')
281  call stat_mfo
282  ltim = dnekclock() - ltim
283  call mntr_tmr_add(stat_tmr_io_id,1,ltim)
284 
285 
286  ! clean up array
287  itmp = lx1**(ldim-stat_rdim)*lelt*stat_lvar
288  call rzero(stat_ruavg,itmp)
289 
290  ! reset averaging parameters
291  stat_atime = 0.0
292  stat_tstart = time
293  endif
294  else
295  ! set averaging parameters
296  stat_atime = 0.0
297  stat_tstart = time
298  endif
299 
300  return
301  end subroutine
302 !======================================================================
309  subroutine stat_init_int1d()
310  implicit none
311 
312  include 'SIZE'
313  include 'FRAMELP'
314  include 'WZ' ! W?M1
315  include 'GEOM' ! ?M1
316  include 'INPUT' ! IFAXIS, IF3D
317  include 'DXYZ' ! D?M1, D?TM1
318  include 'MAP2D'
319  include 'STATD'
320 
321  ! global variables
322  integer mid,mp,nekcomm,nekgroup,nekreal
323  common /nekmpi/ mid,mp,nekcomm,nekgroup,nekreal
324  ! scratch space
325  real lxyzd(lx1,ly1,lz1,lelt,3)
326  common /scrsf/ lxyzd ! coordinate derivatives
327 
328  ! local variables
329  integer il, jl, kl, el ! loop index
330  integer el2 ! index of 2D element
331  real lwm1(lx1) ! wieghts for 1D integration
332 
333  ! global communication
334  integer gs_handle ! gather-scatter handle
335  integer*8 unodes(lx1*lz1*lelt) ! unique local nodes
336 !-----------------------------------------------------------------------
337  if(ifaxis) call mntr_abort(stat_id,
338  $ 'stat_init_int1D; IFAXIS not supported')
339 
340  ! copy wieghts depending on the uniform direction
341  if (map2d_idir.eq.1) then
342  call copy(lwm1,wxm1,nx1)
343  stat_nm1 = nx1
344  stat_nm2 = ny1
345  stat_nm3 = nz1
346  ! get coordinates derivatives d[XYZ]/dr
347  il = ny1*nz1
348  do el = 1, nelt
349  if(map2d_lmap(el).ne.-1) then
350  call mxm(dxm1,nx1,xm1(1,1,1,el),nx1,lxyzd(1,1,1,el,1),il)
351  call mxm(dxm1,nx1,ym1(1,1,1,el),nx1,lxyzd(1,1,1,el,2),il)
352  if (if3d) call mxm(dxm1,nx1,zm1(1,1,1,el),nx1,
353  $ lxyzd(1,1,1,el,3),il)
354  endif
355  enddo
356  elseif (map2d_idir.eq.2) then
357  call copy(lwm1,wym1,ny1)
358  stat_nm1 = ny1
359  stat_nm2 = nx1
360  stat_nm3 = nz1
361  ! get coordinates derivatives d[XYZ]/ds
362  do el = 1, nelt
363  if(map2d_lmap(el).ne.-1) then
364  do il=1, nz1
365  call mxm(xm1(1,1,il,el),nx1,dytm1,ny1,
366  $ lxyzd(1,1,il,el,1),ny1)
367  call mxm(ym1(1,1,il,el),nx1,dytm1,ny1,
368  $ lxyzd(1,1,il,el,2),ny1)
369  if (if3d) call mxm(zm1(1,1,il,el),nx1,dytm1,ny1,
370  $ lxyzd(1,1,il,el,3),ny1)
371  enddo
372  endif
373  enddo
374  else
375  if (if3d) then
376  call copy(lwm1,wzm1,nz1)
377  stat_nm1 = nz1
378  stat_nm2 = nx1
379  stat_nm3 = ny1
380  ! get coordinates derivatives d[XYZ]/dt
381  il = nx1*ny1
382  do el = 1, nelt
383  if(map2d_lmap(el).ne.-1) then
384  call mxm(xm1(1,1,1,el),il,dztm1,nz1,
385  $ lxyzd(1,1,1,el,1),nz1)
386  call mxm(ym1(1,1,1,el),il,dztm1,nz1,
387  $ lxyzd(1,1,1,el,2),nz1)
388  call mxm(zm1(1,1,1,el),il,dztm1,nz1,
389  $ lxyzd(1,1,1,el,3),nz1)
390  endif
391  enddo
392  else
393  call mntr_abort(stat_id,'2D run cannot be z averaged.')
394  endif
395  endif
396 
397  ! for now I assume lx1=stat_nm1=stat_nm2=stat_nm3
398  ! check if that is true
399  if (if3d) then
400  if(lx1.ne.stat_nm1.or.lx1.ne.stat_nm2.or.lx1.ne.stat_nm3) then
401  call mntr_abort(stat_id,
402  $ 'stat_init_int1D; unequal array sizes')
403  endif
404  else
405  if(lx1.ne.stat_nm1.or.lx1.ne.stat_nm2.or.1.ne.stat_nm3) then
406  call mntr_abort(stat_id,
407  $ 'stat_init_int1D; unequal array sizes')
408  endif
409  endif
410 
411  ! get 1D mass matrix ordering directions in such a way that
412  ! the uniform direction corresponds to the the first index
413  il = stat_nm1*stat_nm2*stat_nm3
414  ! get arc length
415  do el = 1, nelt
416  if(map2d_lmap(el).ne.-1) then
417  call vsq(lxyzd(1,1,1,el,1),il)
418  call vsq(lxyzd(1,1,1,el,2),il)
419  if (if3d) call vsq(lxyzd(1,1,1,el,3),il)
420 
421  call add2(lxyzd(1,1,1,el,1),lxyzd(1,1,1,el,2),il)
422  if (if3d) call add2(lxyzd(1,1,1,el,1),lxyzd(1,1,1,el,3),il)
423 
424  call vsqrt(lxyzd(1,1,1,el,1),il)
425  endif
426  enddo
427 
428  il=il*nelt
429  call rzero(stat_bm1d,il)
430 
431  ! reshuffle array
432  call stat_reshufflev(stat_bm1d,lxyzd,nelt)
433 
434  ! multiply by wieghts
435  do el=1, nelt
436  if(map2d_lmap(el).ne.-1) then
437  do kl=1, stat_nm3
438  do jl=1, stat_nm2
439  do il=1, stat_nm1
440  stat_bm1d(il,jl,kl,el) =
441  $ lwm1(il)*stat_bm1d(il,jl,kl,el)
442  enddo
443  enddo
444  enddo
445  endif
446  enddo
447 
448  ! get total line length
449  ! sum contributions from different 3D elements to get
450  ! local arc length
451 
452  il = stat_nm2*stat_nm3*nelt
453  call rzero(stat_abm1d,il)
454 
455  do el = 1, nelv
456  el2 = map2d_lmap(el)
457  if(el2.gt.0) then
458  do kl=1, stat_nm3
459  do jl=1, stat_nm2
460  do il=1, stat_nm1
461  stat_abm1d(jl,kl,el2) = stat_abm1d(jl,kl,el2) +
462  $ stat_bm1d(il,jl,kl,el)
463  enddo
464  enddo
465  enddo
466  endif
467  enddo
468 
469  ! Global communication to sum local contributions for arc lenght
470  ! set up communicator
471  el = stat_nm2*stat_nm3
472  do il = 1,map2d_lnum
473  kl = map2d_gmap(il) - 1
474  do jl=1,el
475  unodes(el*(il-1) + jl) = int(el,8)*int(kl,8) + jl
476  enddo
477  enddo
478  kl = el*map2d_lnum
479  call fgslib_gs_setup(gs_handle,unodes,kl,nekcomm,mp)
480 
481  call fgslib_gs_op(gs_handle,stat_abm1d,1,1,0)
482 
483  ! destroy communicator
484  call fgslib_gs_free (gs_handle)
485 
486  return
487  end subroutine
488 !=======================================================================
496  subroutine stat_reshufflev(rvar, var, nl)
497  implicit none
498 
499  include 'SIZE'
500  include 'MAP2D'
501  include 'STATD'
502 
503  ! argument list
504  real rvar(lx1,ly1,lz1,lelt) !reshuffled array
505  real var(lx1,ly1,lz1,lelt) ! input array
506  integer nl ! element number to reshuffle
507 
508  ! local variables
509  integer il, jl, kl, el ! loop index
510 !-----------------------------------------------------------------------
511  ! if no space averaging copy
512  if (stat_rdim.eq.0) then
513  el = lx1*ly1*lz1*lelt
514  call copy(rvar,var,el)
515  else
516  ! if space averagign swap data
517  if (map2d_idir.eq.1) then
518  do el=1, nl
519  if(map2d_lmap(el).ne.-1) then
520  do kl=1, stat_nm3
521  do jl=1, stat_nm2
522  do il=1, stat_nm1
523  rvar(il,jl,kl,el) = var(il,jl,kl,el)
524  enddo
525  enddo
526  enddo
527  endif
528  enddo
529  elseif (map2d_idir.eq.2) then
530  do el=1, nl
531  if(map2d_lmap(el).ne.-1) then
532  do kl=1, stat_nm3
533  do jl=1, stat_nm2
534  do il=1, stat_nm1
535  rvar(il,jl,kl,el) = var(jl,il,kl,el)
536  enddo
537  enddo
538  enddo
539  endif
540  enddo
541  else
542  do el=1, nl
543  if(map2d_lmap(el).ne.-1) then
544  do kl=1, stat_nm3
545  do jl=1, stat_nm2
546  do il=1, stat_nm1
547  rvar(il,jl,kl,el) = var(jl,kl,il,el)
548  enddo
549  enddo
550  enddo
551  endif
552  enddo
553  endif
554  endif
555 
556  return
557  end subroutine
558 !=======================================================================
565  subroutine stat_compute_1dav1(lvar,npos,alpha,beta)
566  implicit none
567 
568  include 'SIZE'
569  include 'FRAMELP'
570  include 'STATD'
571 
572  ! argument list
573  real lvar(lx1,ly1,lz1,lelt)
574  integer npos
575  real alpha, beta
576 
577  ! global variables
578  real rtmp(lx1,lz1,lelt) ! dummy array
579  common /ctmp0/ rtmp
580 
581  ! local variables
582  integer il, jl, kl, el ! loop index
583 !-----------------------------------------------------------------------
584  ! consistency check
585  if(npos.gt.stat_lvar)
586  $ call mntr_abort(stat_id,'inconsistent npos ')
587 
588  if (stat_rdim.eq.1) then
589  ! zero work array
590  el = lx1*lz1*lelt
591  call rzero(rtmp,el)
592 
593  ! perform 1D integral
594  do el = 1, nelv
595  do kl=1, stat_nm3
596  do jl=1, stat_nm2
597  do il=1, stat_nm1
598  rtmp(jl,kl,el) = rtmp(jl,kl,el) +
599  $ stat_bm1d(il,jl,kl,el)*lvar(il,jl,kl,el)
600  enddo
601  enddo
602  enddo
603  enddo
604 
605  ! time average
606  el = stat_nm2*stat_nm3*nelv
607  call add2sxy(stat_ruavg(1,1,npos),alpha,rtmp,beta,el)
608  else
609  el = lx1**(ldim)*lelt
610  call add2sxy(stat_ruavg(1,1,npos),alpha,lvar,beta,el)
611  endif
612 
613  return
614  end subroutine
615 !=======================================================================
622  subroutine stat_compute_1dav2(lvar1,lvar2,npos,alpha,beta)
623  implicit none
624 
625  include 'SIZE'
626  include 'FRAMELP'
627  include 'STATD'
628 
629  ! argument list
630  real lvar1(lx1,ly1,lz1,lelt), lvar2(lx1,ly1,lz1,lelt)
631  integer npos
632  real alpha, beta
633 
634  ! global variables
635  real rtmp(lx1,lz1,lelt) ! dummy array
636  real rtmp2(lx1,ly1,lz1,lelt) ! dummy array
637  common /ctmp0/ rtmp, rtmp2
638 
639  ! local variables
640  integer il, jl, kl, el ! loop index
641 !-----------------------------------------------------------------------
642  ! consistency check
643  if(npos.gt.stat_lvar)
644  $ call mntr_abort(stat_id,'inconsistent npos ')
645 
646  if (stat_rdim.eq.1) then
647  ! zero work array
648  el = lx1*lz1*lelt
649  call rzero(rtmp,el)
650 
651  ! perform 1D integral
652  do el = 1, nelv
653  do kl=1, stat_nm3
654  do jl=1, stat_nm2
655  do il=1, stat_nm1
656  rtmp(jl,kl,el) = rtmp(jl,kl,el) +
657  $ stat_bm1d(il,jl,kl,el)*
658  $ lvar1(il,jl,kl,el)*lvar2(il,jl,kl,el)
659  enddo
660  enddo
661  enddo
662  enddo
663 
664  ! time average
665  el = stat_nm2*stat_nm3*nelv
666  call add2sxy(stat_ruavg(1,1,npos),alpha,rtmp,beta,el)
667  else
668  el = lx1**(ldim)*lelt
669  call col3(rtmp2,lvar1,lvar2,el)
670  call add2sxy(stat_ruavg(1,1,npos),alpha,rtmp2,beta,el)
671  endif
672 
673  return
674  end subroutine
675 !=======================================================================
678  subroutine stat_gs_sum
679  implicit none
680 
681  include 'SIZE'
682  include 'FRAMELP'
683  include 'PARALLEL'
684  include 'INPUT' ! if3d
685  include 'MAP2D'
686  include 'STATD'
687 
688  ! global variables
689  integer mid,mp,nekcomm,nekgroup,nekreal
690  common /nekmpi/ mid,mp,nekcomm,nekgroup,nekreal
691 
692  ! local variables
693  integer gs_handle ! gather-scatter handle
694  integer*8 unodes(lx1*lz1*lelt) ! unique local nodes
695  integer il, jl, el ! loop index
696  integer el2 ! index of 2D element
697  integer itmp1, itmp2
698  real rtmp_ruavg(lx1*lz1,lelt) ! tmp array for local 2D element aggragation
699 !-----------------------------------------------------------------------
700  ! if no space averaging return
701  if (stat_rdim.eq.0) return
702 
703  ! stamp logs
704  call mntr_log(stat_id,lp_vrb,'Global statistics summation.')
705 
706  ! perform local summation of 2D contributions; not required for 3D version
707  do il = 1, stat_lvar
708  el = lx1*lz1*lelt
709  call rzero(rtmp_ruavg,el)
710  do el=1,nelv
711  el2 = map2d_lmap(el)
712  if(el2.gt.0) then
713  do jl = 1,stat_nm2*stat_nm3
714  rtmp_ruavg(jl,el2) = rtmp_ruavg(jl,el2) +
715  $ stat_ruavg(jl,el,il)
716  enddo
717  endif
718  enddo
719  ! copy data back
720  el = stat_nm2*stat_nm3*map2d_lnum
721  call copy(stat_ruavg(1,1,il),rtmp_ruavg,el)
722  enddo
723 
724  ! set up communicator
725  itmp1 = stat_nm2*stat_nm3
726  do il = 1, map2d_lnum
727  itmp2 = map2d_gmap(il) - 1
728  do jl=1,itmp1
729  unodes(itmp1*(il-1) + jl) = int(itmp1,8)*int(itmp2,8) + jl
730  enddo
731  enddo
732  itmp2 = itmp1*map2d_lnum
733  call fgslib_gs_setup(gs_handle,unodes,itmp2,nekcomm,mp)
734 
735  ! sum variables
736  do il=1,stat_lvar
737  call fgslib_gs_op(gs_handle,stat_ruavg(1,1,il),1,1,0)
738  enddo
739 
740  ! destroy communicator
741  call fgslib_gs_free (gs_handle)
742 
743  ! divide data by arc length
744  if (if3d) then
745  itmp1=stat_nm2*stat_nm3
746  do il=1, map2d_lnum
747  if (map2d_own(il).eq.nid) then
748  do jl=1, stat_lvar
749  call invcol2(stat_ruavg(1,il,jl),
750  $ stat_abm1d(1,1,il),itmp1)
751  enddo
752  endif
753  enddo
754  endif
755 
756  return
757  end subroutine
758 !=======================================================================
762  subroutine stat_compute()
763  implicit none
764 
765  include 'SIZE'
766  include 'FRAMELP'
767  include 'SOLN'
768  include 'TSTEP'
769  include 'STATD' ! Variables from the statistics
770  include 'INPUT' ! if3d
771 
772  ! global variables
773  ! work arrays
774  real slvel(LX1,LY1,LZ1,LELT,3), slp(LX1,LY1,LZ1,LELT)
775  common /scrmg/ slvel, slp
776  real tmpvel(LX1,LY1,LZ1,LELT,3), tmppr(LX1,LY1,LZ1,LELT)
777  common /scruz/ tmpvel, tmppr
778  real dudx(LX1,LY1,LZ1,LELT,3) ! du/dx, du/dy and du/dz
779  real dvdx(LX1,LY1,LZ1,LELT,3) ! dv/dx, dv/dy and dv/dz
780  real dwdx(LX1,LY1,LZ1,LELT,3) ! dw/dx, dw/dy and dw/dz
781  common /scrns/ dudx, dvdx
782  common /scrsf/ dwdx
783 
784  ! local variables
785  integer npos ! position in STAT_RUAVG
786  real alpha, beta,dtime ! time averaging parameters
787  integer lnvar ! count number of variables
788  integer i ! loop index
789  integer itmp ! dummy variable
790  real rtmp ! dummy variable
791  integer ntot
792 
793 !-----------------------------------------------------------------------
794  ! stamp logs
795  call mntr_log(stat_id,lp_vrb,'Average fields.')
796 
797  ! Calculate time span of current statistical sample
798  dtime=time-stat_atime-stat_tstart
799 
800  ! Update total time over which the current stat file is averaged
801  stat_atime=time-stat_tstart
802 
803  ! Time average is compuated as:
804  ! Accumulated=alpha*Accumulated+beta*New
805  ! Calculate alpha and beta
806  beta=dtime/stat_atime
807  alpha=1.0-beta
808 
809  ! Map pressure to velocity mesh
810  call mappr(tmppr,pr,tmpvel(1,1,1,1,2),tmpvel(1,1,1,1,3))
811 
812  ! Compute derivative tensor and normalise pressure
813  call user_stat_trnsv(tmpvel,dudx,dvdx,dwdx,slvel,tmppr)
814 
815  ! reset varaible counter
816  lnvar = 0
817 
818  ! reshuffle arrays
819  ! velocity
820  call stat_reshufflev(slvel(1,1,1,1,1),tmpvel(1,1,1,1,1),nelv)
821  call stat_reshufflev(slvel(1,1,1,1,2),tmpvel(1,1,1,1,2),nelv)
822  if (if3d) call stat_reshufflev(slvel(1,1,1,1,3),
823  $ tmpvel(1,1,1,1,3),nelv)
824 
825  ! pressure
826  call stat_reshufflev(slp,tmppr,nelv)
827 
828  ! reshuffle velocity derivatives
829  ! VX
830  call stat_reshufflev(tmpvel(1,1,1,1,1),dudx(1,1,1,1,1),nelv)
831  call stat_reshufflev(tmpvel(1,1,1,1,2),dudx(1,1,1,1,2),nelv)
832  if (if3d) call stat_reshufflev(tmpvel(1,1,1,1,3),
833  $ dudx(1,1,1,1,3),nelv)
834  ! copy
835  itmp = lx1*ly1*lz1*lelt*ldim
836  call copy(dudx,tmpvel,itmp)
837 
838  ! VY
839  call stat_reshufflev(tmpvel(1,1,1,1,1),dvdx(1,1,1,1,1),nelv)
840  call stat_reshufflev(tmpvel(1,1,1,1,2),dvdx(1,1,1,1,2),nelv)
841  if (if3d) call stat_reshufflev(tmpvel(1,1,1,1,3),
842  $ dvdx(1,1,1,1,3),nelv)
843  ! copy
844  itmp = lx1*ly1*lz1*lelt*ldim
845  call copy(dvdx,tmpvel,itmp)
846 
847  ! VZ
848  if (if3d) then
849  call stat_reshufflev(tmpvel(1,1,1,1,1),dwdx(1,1,1,1,1),nelv)
850  call stat_reshufflev(tmpvel(1,1,1,1,2),dwdx(1,1,1,1,2),nelv)
851  call stat_reshufflev(tmpvel(1,1,1,1,3),dwdx(1,1,1,1,3),nelv)
852  ! copy
853  itmp = lx1*ly1*lz1*lelt*ldim
854  call copy(dwdx,tmpvel,itmp)
855  endif
856 
857 !=======================================================================
858  ! Computation of statistics
859 
860  ! <u>t
861  lnvar = lnvar + 1
862  npos = lnvar
863  call stat_compute_1dav1(slvel(1,1,1,1,1),npos,alpha,beta)
864 
865  ! <v>t
866  lnvar = lnvar + 1
867  npos = lnvar
868  call stat_compute_1dav1(slvel(1,1,1,1,2),npos,alpha,beta)
869 
870  ! <w>t
871  lnvar = lnvar + 1
872  npos = lnvar
873  call stat_compute_1dav1(slvel(1,1,1,1,3),npos,alpha,beta)
874 
875  ! <p>t
876  lnvar = lnvar + 1
877  npos = lnvar
878  call stat_compute_1dav1(slp(1,1,1,1),npos,alpha,beta)
879 
880 !-----------------------------------------------------------------------
881  ! <uu>t
882  lnvar = lnvar + 1
883  npos = lnvar
884  call stat_compute_1dav2(slvel(1,1,1,1,1),slvel(1,1,1,1,1),
885  $ npos,alpha,beta)
886 
887  ! <vv>t
888  lnvar = lnvar + 1
889  npos = lnvar
890  call stat_compute_1dav2(slvel(1,1,1,1,2),slvel(1,1,1,1,2),
891  $ npos,alpha,beta)
892 
893  ! <ww>t
894  lnvar = lnvar + 1
895  npos = lnvar
896  call stat_compute_1dav2(slvel(1,1,1,1,3),slvel(1,1,1,1,3),
897  $ npos,alpha,beta)
898 
899  ! <pp>t
900  lnvar = lnvar + 1
901  npos = lnvar
902  call stat_compute_1dav2(slp(1,1,1,1),slp(1,1,1,1),
903  $ npos,alpha,beta)
904 
905 !-----------------------------------------------------------------------
906  ! <uv>t
907  lnvar = lnvar + 1
908  npos = lnvar
909  call stat_compute_1dav2(slvel(1,1,1,1,1),slvel(1,1,1,1,2),
910  $ npos,alpha,beta)
911 
912  ! <vw>t
913  lnvar = lnvar + 1
914  npos = lnvar
915  call stat_compute_1dav2(slvel(1,1,1,1,2),slvel(1,1,1,1,3),
916  $ npos,alpha,beta)
917 
918  ! <uw>t
919  lnvar = lnvar + 1
920  npos = lnvar
921  call stat_compute_1dav2(slvel(1,1,1,1,1),slvel(1,1,1,1,3),
922  $ npos,alpha,beta)
923 
924 !-----------------------------------------------------------------------
925  ! <pu>t
926  lnvar = lnvar + 1
927  npos = lnvar
928  call stat_compute_1dav2(slp(1,1,1,1),slvel(1,1,1,1,1),
929  $ npos,alpha,beta)
930 
931  ! <pv>t
932  lnvar = lnvar + 1
933  npos = lnvar
934  call stat_compute_1dav2(slp(1,1,1,1),slvel(1,1,1,1,2),
935  $ npos,alpha,beta)
936 
937  ! <pw>t
938  lnvar = lnvar + 1
939  npos = lnvar
940  call stat_compute_1dav2(slp(1,1,1,1),slvel(1,1,1,1,3),
941  $ npos,alpha,beta)
942 
943 !-----------------------------------------------------------------------
944  ! <pdudx>t
945  lnvar = lnvar + 1
946  npos = lnvar
947  call stat_compute_1dav2(slp(1,1,1,1),dudx(1,1,1,1,1),
948  $ npos,alpha,beta)
949 
950  ! <pdudy>t
951  lnvar = lnvar + 1
952  npos = lnvar
953  call stat_compute_1dav2(slp(1,1,1,1),dudx(1,1,1,1,2),
954  $ npos,alpha,beta)
955 
956  ! <pdudz>t
957  lnvar = lnvar + 1
958  npos = lnvar
959  call stat_compute_1dav2(slp(1,1,1,1),dudx(1,1,1,1,3),
960  $ npos,alpha,beta)
961 
962 !-----------------------------------------------------------------------
963  ! <pdvdx>t
964  lnvar = lnvar + 1
965  npos = lnvar
966  call stat_compute_1dav2(slp(1,1,1,1),dvdx(1,1,1,1,1),
967  $ npos,alpha,beta)
968 
969  ! <pdvdy>t
970  lnvar = lnvar + 1
971  npos = lnvar
972  call stat_compute_1dav2(slp(1,1,1,1),dvdx(1,1,1,1,2),
973  $ npos,alpha,beta)
974 
975  ! <pdvdz>t
976  lnvar = lnvar + 1
977  npos = lnvar
978  call stat_compute_1dav2(slp(1,1,1,1),dvdx(1,1,1,1,3),
979  $ npos,alpha,beta)
980 
981 !-----------------------------------------------------------------------
982  ! <pdwdx>t
983  lnvar = lnvar + 1
984  npos = lnvar
985  call stat_compute_1dav2(slp(1,1,1,1),dwdx(1,1,1,1,1),
986  $ npos,alpha,beta)
987 
988  ! <pdwdy>t
989  lnvar = lnvar + 1
990  npos = lnvar
991  call stat_compute_1dav2(slp(1,1,1,1),dwdx(1,1,1,1,2),
992  $ npos,alpha,beta)
993 
994  ! <pdwdz>t
995  lnvar = lnvar + 1
996  npos = lnvar
997  call stat_compute_1dav2(slp(1,1,1,1),dwdx(1,1,1,1,3),
998  $ npos,alpha,beta)
999 
1000 !-----------------------------------------------------------------------
1001  ! UU, VV, WW
1002  itmp = lx1*ly1*lz1*lelt*ldim
1003  call col3(tmpvel(1,1,1,1,1),slvel(1,1,1,1,1),slvel(1,1,1,1,1),
1004  $ itmp)
1005 
1006  ! <uuu>t
1007  lnvar = lnvar + 1
1008  npos = lnvar
1009  call stat_compute_1dav2(slvel(1,1,1,1,1),tmpvel(1,1,1,1,1),
1010  $ npos,alpha,beta)
1011 
1012  ! <vvv>t
1013  lnvar = lnvar + 1
1014  npos = lnvar
1015  call stat_compute_1dav2(slvel(1,1,1,1,2),tmpvel(1,1,1,1,2),
1016  $ npos,alpha,beta)
1017 
1018  ! <www>t
1019  lnvar = lnvar + 1
1020  npos = lnvar
1021  call stat_compute_1dav2(slvel(1,1,1,1,3),tmpvel(1,1,1,1,3),
1022  $ npos,alpha,beta)
1023 
1024 !-----------------------------------------------------------------------
1025  ! <uuv>t
1026  lnvar = lnvar + 1
1027  npos = lnvar
1028  call stat_compute_1dav2(tmpvel(1,1,1,1,1),slvel(1,1,1,1,2),
1029  $ npos,alpha,beta)
1030 
1031  ! <uuw>t
1032  lnvar = lnvar + 1
1033  npos = lnvar
1034  call stat_compute_1dav2(tmpvel(1,1,1,1,1),slvel(1,1,1,1,3),
1035  $ npos,alpha,beta)
1036 
1037  ! <vvu>t
1038  lnvar = lnvar + 1
1039  npos = lnvar
1040  call stat_compute_1dav2(tmpvel(1,1,1,1,2),slvel(1,1,1,1,1),
1041  $ npos,alpha,beta)
1042 
1043  ! <vvw>t
1044  lnvar = lnvar + 1
1045  npos = lnvar
1046  call stat_compute_1dav2(tmpvel(1,1,1,1,2),slvel(1,1,1,1,3),
1047  $ npos,alpha,beta)
1048 
1049 !-----------------------------------------------------------------------
1050  ! <wwu>t
1051  lnvar = lnvar + 1
1052  npos = lnvar
1053  call stat_compute_1dav2(tmpvel(1,1,1,1,3),slvel(1,1,1,1,1),
1054  $ npos,alpha,beta)
1055 
1056  ! <wwv>t
1057  lnvar = lnvar + 1
1058  npos = lnvar
1059  call stat_compute_1dav2(tmpvel(1,1,1,1,3),slvel(1,1,1,1,2),
1060  $ npos,alpha,beta)
1061 
1062 !-----------------------------------------------------------------------
1063  ! <ppp>t
1064  lnvar = lnvar + 1
1065  npos = lnvar
1066  itmp = lx1*ly1*lz1*lelt
1067  call col3(tmppr(1,1,1,1),slp(1,1,1,1),slp(1,1,1,1),
1068  $ itmp)
1069  call stat_compute_1dav2(tmppr(1,1,1,1),slp(1,1,1,1),
1070  $ npos,alpha,beta)
1071 
1072 !-----------------------------------------------------------------------
1073  ! <pppp>t
1074  lnvar = lnvar + 1
1075  npos = lnvar
1076  call stat_compute_1dav2(tmppr(1,1,1,1),tmppr(1,1,1,1),
1077  $ npos,alpha,beta)
1078 
1079 !-----------------------------------------------------------------------
1080  ! <uvw>t
1081  lnvar = lnvar + 1
1082  npos = lnvar
1083  ! copy uv to tmppr (do not need pp anymore)
1084  itmp = lx1*ly1*lz1*lelt
1085  call col3(tmppr(1,1,1,1),slvel(1,1,1,1,1),slvel(1,1,1,1,2),
1086  $ itmp)
1087  call stat_compute_1dav2(tmppr(1,1,1,1),slvel(1,1,1,1,3),
1088  $ npos,alpha,beta)
1089 
1090 !-----------------------------------------------------------------------
1091  ! <uuuu>t
1092  lnvar = lnvar + 1
1093  npos = lnvar
1094  call stat_compute_1dav2(tmpvel(1,1,1,1,1),tmpvel(1,1,1,1,1),
1095  $ npos,alpha,beta)
1096 
1097  ! <vvvv>t
1098  lnvar = lnvar + 1
1099  npos = lnvar
1100  call stat_compute_1dav2(tmpvel(1,1,1,1,2),tmpvel(1,1,1,1,2),
1101  $ npos,alpha,beta)
1102 
1103  ! <wwww>t
1104  lnvar = lnvar + 1
1105  npos = lnvar
1106  call stat_compute_1dav2(tmpvel(1,1,1,1,3),tmpvel(1,1,1,1,3),
1107  $ npos,alpha,beta)
1108 
1109 !-----------------------------------------------------------------------
1110  ! <e11>t : (du/dx)^2 + (du/dy)^2 + (du/dz)^2
1111  itmp = lx1*ly1*lz1*lelt*ldim
1112  call col3(tmpvel(1,1,1,1,1),dudx(1,1,1,1,1),dudx(1,1,1,1,1),
1113  $ itmp)
1114  itmp = lx1*ly1*lz1*lelt
1115  call add2(tmpvel(1,1,1,1,1),tmpvel(1,1,1,1,2),itmp)
1116  call add2(tmpvel(1,1,1,1,1),tmpvel(1,1,1,1,3),itmp)
1117  lnvar = lnvar + 1
1118  npos = lnvar
1119  call stat_compute_1dav1(tmpvel(1,1,1,1,1),npos,alpha,beta)
1120 
1121  ! <e22>t: (dv/dx)^2 + (dv/dy)^2 + (dv/dz)^2
1122  itmp = lx1*ly1*lz1*lelt*ldim
1123  call col3(tmpvel(1,1,1,1,1),dvdx(1,1,1,1,1),dvdx(1,1,1,1,1),
1124  $ itmp)
1125  itmp = lx1*ly1*lz1*lelt
1126  call add2(tmpvel(1,1,1,1,1),tmpvel(1,1,1,1,2),itmp)
1127  call add2(tmpvel(1,1,1,1,1),tmpvel(1,1,1,1,3),itmp)
1128  lnvar = lnvar + 1
1129  npos = lnvar
1130  call stat_compute_1dav1(tmpvel(1,1,1,1,1),npos,alpha,beta)
1131 
1132  ! <e33>t: (dw/dx)^2 + (dw/dy)^2 + (dw/dz)^2
1133  itmp = lx1*ly1*lz1*lelt*ldim
1134  call col3(tmpvel(1,1,1,1,1),dwdx(1,1,1,1,1),dwdx(1,1,1,1,1),
1135  $ itmp)
1136  itmp = lx1*ly1*lz1*lelt
1137  call add2(tmpvel(1,1,1,1,1),tmpvel(1,1,1,1,2),itmp)
1138  call add2(tmpvel(1,1,1,1,1),tmpvel(1,1,1,1,3),itmp)
1139  lnvar = lnvar + 1
1140  npos = lnvar
1141  call stat_compute_1dav1(tmpvel(1,1,1,1,1),npos,alpha,beta)
1142 
1143 !-----------------------------------------------------------------------
1144  ! <e12>t: (du/dx)*(dv/dx) + (du/dy)*(dv/dy) + (du/dz)*(dv/dz)
1145  itmp = lx1*ly1*lz1*lelt*ldim
1146  call col3(tmpvel(1,1,1,1,1),dudx(1,1,1,1,1),dvdx(1,1,1,1,1),
1147  $ itmp)
1148  itmp = lx1*ly1*lz1*lelt
1149  call add2(tmpvel(1,1,1,1,1),tmpvel(1,1,1,1,2),itmp)
1150  call add2(tmpvel(1,1,1,1,1),tmpvel(1,1,1,1,3),itmp)
1151  lnvar = lnvar + 1
1152  npos = lnvar
1153  call stat_compute_1dav1(tmpvel(1,1,1,1,1),npos,alpha,beta)
1154 
1155  ! <e13>t: (du/dx)*(dw/dx) + (du/dy)*(dw/dy) + (du/dz)*(dw/dz)
1156  itmp = lx1*ly1*lz1*lelt*ldim
1157  call col3(tmpvel(1,1,1,1,1),dudx(1,1,1,1,1),dwdx(1,1,1,1,1),
1158  $ itmp)
1159  itmp = lx1*ly1*lz1*lelt
1160  call add2(tmpvel(1,1,1,1,1),tmpvel(1,1,1,1,2),itmp)
1161  call add2(tmpvel(1,1,1,1,1),tmpvel(1,1,1,1,3),itmp)
1162  lnvar = lnvar + 1
1163  npos = lnvar
1164  call stat_compute_1dav1(tmpvel(1,1,1,1,1),npos,alpha,beta)
1165 
1166  ! <e23>t: (dv/dx)*(dw/dx) + (dv/dy)*(dw/dy) + (dv/dz)*(dw/dz)
1167  itmp = lx1*ly1*lz1*lelt*ldim
1168  call col3(tmpvel(1,1,1,1,1),dvdx(1,1,1,1,1),dwdx(1,1,1,1,1),
1169  $ itmp)
1170  itmp = lx1*ly1*lz1*lelt
1171  call add2(tmpvel(1,1,1,1,1),tmpvel(1,1,1,1,2),itmp)
1172  call add2(tmpvel(1,1,1,1,1),tmpvel(1,1,1,1,3),itmp)
1173  lnvar = lnvar + 1
1174  npos = lnvar
1175  call stat_compute_1dav1(tmpvel(1,1,1,1,1),npos,alpha,beta)
1176 
1177 !=======================================================================
1178  !End of local compute
1179 
1180  ! save number of variables
1181  stat_nvar = lnvar
1182 
1183  return
1184  end subroutine
1185 !=======================================================================
map2d_register
subroutine map2d_register()
Register 2D mapping routines.
Definition: map2D.f:13
map2d_init
subroutine map2d_init()
Initilise map2d module.
Definition: map2D.f:70
mntr_tmr_is_name_reg
subroutine mntr_tmr_is_name_reg(mid, mname)
Check if timer name is registered and return its id.
Definition: mntrtmr.f:146
mntr_warn
subroutine mntr_warn(mid, logs)
Write warning message.
Definition: mntrlog.f:803
mntr_tmr_add
subroutine mntr_tmr_add(mid, icount, time)
Check if timer id is registered. This operation is performed locally.
Definition: mntrtmr.f:237
mntr_mod_is_name_reg
subroutine mntr_mod_is_name_reg(mid, mname)
Check if module name is registered and return its id.
Definition: mntrlog.f:459
mntr_abort
subroutine mntr_abort(mid, logs)
Abort simulation.
Definition: mntrlog.f:837
mntr_log
subroutine mntr_log(mid, priority, logs)
Write log message.
Definition: mntrlog.f:600
mntr_mod_reg
subroutine mntr_mod_reg(mid, pmid, mname, mdscr)
Register new module.
Definition: mntrlog.f:346
mntr_tmr_reg
subroutine mntr_tmr_reg(mid, pmid, modid, mname, mdscr, ifsum)
Register new timer.
Definition: mntrtmr.f:16
rprm_rp_get
subroutine rprm_rp_get(ipval, rpval, lpval, cpval, rpid, ptype)
Get runtime parameter form active section. This operation is performed locally.
Definition: rprm.f:883
rprm_rp_reg
subroutine rprm_rp_reg(rpid, mid, pname, pdscr, ptype, ipval, rpval, lpval, cpval)
Register new runtime parameter.
Definition: rprm.f:484
rprm_sec_set_act
subroutine rprm_sec_set_act(ifact, rpid)
Set section's activation flag. Master value is broadcasted.
Definition: rprm.f:422
rprm_sec_reg
subroutine rprm_sec_reg(rpid, mid, pname, pdscr)
Register new parameter section.
Definition: rprm.f:165
stat_avg
subroutine stat_avg
Main interface of statistics module.
Definition: stat.f:237
stat_gs_sum
subroutine stat_gs_sum
Global statistics summation.
Definition: stat.f:679
stat_init
subroutine stat_init()
Initilise statistics module.
Definition: stat.f:111
stat_register
subroutine stat_register()
Register 2D and 3D statistics module.
Definition: stat.f:13
stat_is_initialised
logical function stat_is_initialised()
Check if module was initialised.
Definition: stat.f:222
stat_compute_1dav2
subroutine stat_compute_1dav2(lvar1, lvar2, npos, alpha, beta)
Perform local 1D integration on multiplication of 2 variables.
Definition: stat.f:623
stat_mfo
subroutine stat_mfo()
Main interface for saving statistics.
Definition: stat_IO.f:13
stat_compute
subroutine stat_compute()
Compute statistics.
Definition: stat.f:763
stat_init_int1d
subroutine stat_init_int1d()
Get local integration coefficients.
Definition: stat.f:310
stat_reshufflev
subroutine stat_reshufflev(rvar, var, nl)
Array reshuffle.
Definition: stat.f:497
stat_end
subroutine stat_end()
Finalise statistics module.
Definition: stat.f:186
stat_compute_1dav1
subroutine stat_compute_1dav1(lvar, npos, alpha, beta)
Perform local 1D integration on 1 variable.
Definition: stat.f:566
col3
subroutine col3(a, b, c, n)
Definition: math.f:598
invcol2
subroutine invcol2(a, b, n)
Definition: math.f:73
add2
subroutine add2(a, b, n)
Definition: math.f:622
copy
subroutine copy(a, b, n)
Definition: math.f:260
add2sxy
subroutine add2sxy(x, a, y, b, n)
Definition: math.f:1574
vsqrt
subroutine vsqrt(A, N)
Definition: math.f:41
vsq
subroutine vsq(A, N)
Definition: math.f:31
rzero
subroutine rzero(a, n)
Definition: math.f:208
mxm
subroutine mxm(a, n1, b, n2, c, n3)
Definition: mxm_wrapper.f:2
mappr
subroutine mappr(pm1, pm2, pa, pb)
Definition: navier5.f:237