KTH framework for Nek5000 toolboxes; testing version  0.0.1
pstat2D.f
Go to the documentation of this file.
1 
6 !=======================================================================
10  subroutine pstat2d_register()
11  implicit none
12 
13  include 'SIZE'
14  include 'INPUT'
15  include 'FRAMELP'
16  include 'PSTAT2D'
17 
18  ! local variables
19  integer lpmid, il
20  real ltim
21  character*2 str
22 
23  ! functions
24  real dnekclock
25 !-----------------------------------------------------------------------
26  ! timing
27  ltim = dnekclock()
28 
29  ! check if the current module was already registered
30  call mntr_mod_is_name_reg(lpmid,pstat_name)
31  if (lpmid.gt.0) then
32  call mntr_warn(lpmid,
33  $ 'module ['//trim(pstat_name)//'] already registered')
34  return
35  endif
36 
37  ! find parent module
38  call mntr_mod_is_name_reg(lpmid,'FRAME')
39  if (lpmid.le.0) then
40  lpmid = 1
41  call mntr_abort(lpmid,
42  $ 'parent module ['//'FRAME'//'] not registered')
43  endif
44 
45  ! register module
46  call mntr_mod_reg(pstat_id,lpmid,pstat_name,
47  $ 'Post processing for statistics')
48 
49  ! register timers
50  call mntr_tmr_is_name_reg(lpmid,'FRM_TOT')
51  ! total time
52  call mntr_tmr_reg(pstat_tmr_tot_id,lpmid,pstat_id,
53  $ 'PSTAT_TOT','Pstat total time',.false.)
54  lpmid = pstat_tmr_tot_id
55  ! initialisation
56  call mntr_tmr_reg(pstat_tmr_ini_id,lpmid,pstat_id,
57  $ 'PSTAT_INI','Pstat initialisation time',.true.)
58  ! averaging
59  call mntr_tmr_reg(pstat_tmr_avg_id,lpmid,pstat_id,
60  $ 'PSTAT_AVG','Pstat averaging time',.true.)
61  ! derivative calculation
62  call mntr_tmr_reg(pstat_tmr_der_id,lpmid,pstat_id,
63  $ 'PSTAT_DER','Pstat derivative calculation time',.true.)
64  ! interpolation
65  call mntr_tmr_reg(pstat_tmr_int_id,lpmid,pstat_id,
66  $ 'PSTAT_INT','Pstat interpolation time',.true.)
67 
68  ! register and set active section
69  call rprm_sec_reg(pstat_sec_id,pstat_id,'_'//adjustl(pstat_name),
70  $ 'Runtime paramere section for pstat module')
71  call rprm_sec_set_act(.true.,pstat_sec_id)
72 
73  ! register parameters
74  call rprm_rp_reg(pstat_amr_irnr_id,pstat_sec_id,'AMR_NREF',
75  $ 'Nr. of initial refinemnt (AMR only)',rpar_int,1,0.0,.false.,' ')
76  call rprm_rp_reg(pstat_ffile_id,pstat_sec_id,'STS_FFILE',
77  $ 'First stat file number',rpar_int,1,0.0,.false.,' ')
78  call rprm_rp_reg(pstat_nfile_id,pstat_sec_id,'STS_NFILE',
79  $ 'Last stat file number',rpar_int,1,0.0,.false.,' ')
80  call rprm_rp_reg(pstat_stime_id,pstat_sec_id,'STS_STIME',
81  $ 'Statistics starting time',rpar_real,1,0.0,.false.,' ')
82  call rprm_rp_reg(pstat_nstep_id,pstat_sec_id,'STS_NSTEP',
83  $ 'Number of steps between averaging (in sts file)',
84  $ rpar_int,10,0.0,.false.,' ')
85 
86  ! set initialisation flag
87  pstat_ifinit=.false.
88 
89  ! timing
90  ltim = dnekclock() - ltim
91  call mntr_tmr_add(pstat_tmr_tot_id,1,ltim)
92 
93  return
94  end subroutine
95 !=======================================================================
99  subroutine pstat2d_init()
100  implicit none
101 
102  include 'SIZE'
103  include 'FRAMELP'
104  include 'PSTAT2D'
105 
106  ! local variables
107  integer itmp, il
108  real rtmp, ltim
109  logical ltmp
110  character*20 ctmp
111 
112  ! functions
113  real dnekclock
114 !-----------------------------------------------------------------------
115 ! check if the module was already initialised
116  if (pstat_ifinit) then
117  call mntr_warn(pstat_id,
118  $ 'module ['//trim(pstat_name)//'] already initiaised.')
119  return
120  endif
121 
122  ! timing
123  ltim = dnekclock()
124 
125  ! get runtime parameters
126  call rprm_rp_get(itmp,rtmp,ltmp,ctmp,pstat_amr_irnr_id,rpar_int)
127  pstat_amr_irnr = abs(itmp)
128  call rprm_rp_get(itmp,rtmp,ltmp,ctmp,pstat_ffile_id,rpar_int)
129  pstat_ffile = abs(itmp)
130  call rprm_rp_get(itmp,rtmp,ltmp,ctmp,pstat_nfile_id,rpar_int)
131  pstat_nfile = abs(itmp)
132  call rprm_rp_get(itmp,rtmp,ltmp,ctmp,pstat_stime_id,rpar_real)
133  pstat_stime = rtmp
134  call rprm_rp_get(itmp,rtmp,ltmp,ctmp,pstat_nstep_id,rpar_int)
135  pstat_nstep = abs(itmp)
136 
137  ! set field swapping array
138  do il = 1,26
139  pstat_swfield(il) = il
140  enddo
141  do il = 27,32
142  pstat_swfield(il) = il+1
143  enddo
144  pstat_swfield(33) = 27
145  pstat_swfield(34) = 38
146  do il = 35,38
147  pstat_swfield(il) = il-1
148  enddo
149  do il = 39,pstat_svar
150  pstat_swfield(il) = il
151  enddo
152 #ifdef AMR
153  ! in case we are dealing with amr mesh first generate the mesh
154  ! read element centre data
155  call mntr_log(pstat_id,lp_inf,'Updating 2D mesh structure')
156  call pstat2d_mfi_crd2d(pstat_ffile)
157  ! perform element ordering to get 2D mesh corresponing to 3D projected one
158  ! in case of AMR this performs refinemnt of 2D mesh as well
159  call pstat2d_mesh_amr()
160 #endif
161  ! everything is initialised
162  pstat_ifinit=.true.
163 
164  ! timing
165  ltim = dnekclock() - ltim
166  call mntr_tmr_add(pstat_tmr_ini_id,1,ltim)
167 
168  return
169  end subroutine
170 !=======================================================================
174  logical function pstat2d_is_initialised()
175  implicit none
176 
177  include 'SIZE'
178  include 'PSTAT2D'
179 !-----------------------------------------------------------------------
180  pstat2d_is_initialised = pstat_ifinit
181 
182  return
183  end function
184 !=======================================================================
187  subroutine pstat2d_main
188  implicit none
189 
190  include 'SIZE'
191  include 'FRAMELP'
192  include 'PSTAT2D'
193 
194  ! local variables
195  integer il
196 
197  ! simple timing
198  real ltim
199 
200  ! functions
201  real dnekclock
202 !-----------------------------------------------------------------------
203  ! read and average fields
204  ltim = dnekclock()
205  call mntr_log(pstat_id,lp_inf,'Field averaging')
206  call pstat2d_sts_avg
207  ltim = dnekclock() - ltim
208  call mntr_tmr_add(pstat_tmr_avg_id,1,ltim)
209 
210  ! calculate derivatives
211  ltim = dnekclock()
212  call mntr_log(pstat_id,lp_inf,'Derivative calculation')
213  call pstat2d_deriv
214  ltim = dnekclock() - ltim
215  call mntr_tmr_add(pstat_tmr_der_id,1,ltim)
216 
217  ! interpolate into the set of points
218  ltim = dnekclock()
219  call mntr_log(pstat_id,lp_inf,'Point interpolation')
220  call pstat2d_interp
221  ltim = dnekclock() - ltim
222  call mntr_tmr_add(pstat_tmr_int_id,1,ltim)
223 
224  return
225  end subroutine
226 !=======================================================================
229  subroutine pstat2d_mesh_amr
230  implicit none
231 #ifdef AMR
232  include 'SIZE'
233  include 'PARALLEL'
234  include 'FRAMELP'
235  include 'PSTAT2D'
236 
237  ! local variables
238  integer ierr
239  integer il, inf_cnt
240 
241  ! functions
242  integer iglsum
243 !-----------------------------------------------------------------------
244  ! reset amr level falg
245  do il = 1, lelt
246  pstat_refl(il) = 0
247  enddo
248 
249  ! initial refinement
250  do il=1,pstat_amr_irnr
251  call amr_refinement()
252  enddo
253 
254  ! infinit loop
255  inf_cnt = 0
256  do
257  ! get element mapping
258  call pstat2d_mesh_map(.true.)
259 
260  ! perform refinement
261  call amr_refinement()
262 
263  inf_cnt = inf_cnt + 1
264  if (inf_cnt.gt.99) call mntr_abort(pstat_id,
265  $ 'Infinit loop too long; possible problem with mark check.')
266 
267  pstat_elmod = iglsum(pstat_elmod,1)
268  if (nelgv.eq.pstat_nelg.and.pstat_elmod.eq.0) then
269  call mntr_logi(pstat_id,lp_prd,
270  $ 'Finished refinement; cycles number :', inf_cnt)
271  exit
272  endif
273 
274  enddo
275 #endif
276  return
277  end subroutine
278 !=======================================================================
282  subroutine pstat2d_mesh_map(ifrref)
283  implicit none
284 
285  include 'SIZE'
286  include 'GEOM'
287  include 'PARALLEL'
288  include 'INPUT'
289  include 'FRAMELP'
290  include 'PSTAT2D'
291 
292  ! global data structures
293  integer mid,mp,nekcomm,nekgroup,nekreal
294  common /nekmpi/ mid,mp,nekcomm,nekgroup,nekreal
295 
296  ! argument list
297  logical ifrref
298 
299  ! local variables
300  integer ierr
301  integer il, jl, kl, inf_cnt
302  integer ifpts ! findpts flag
303  real tol ! interpolation tolerance
304  integer nt
305  integer npt_max ! max communication set
306  integer nxf, nyf, nzf ! fine mesh for bb-test
307  real bb_t ! relative size to expand bounding boxes by
308  real toldist
309  parameter(toldist = 5e-6)
310 
311  integer rcode(lelt), proc(lelt),elid(lelt)
312  real dist(lelt), rst(ldim*lelt)
313 
314  integer nfail, npass
315 
316  ! for sorting
317  integer idim, isdim
318  parameter(idim=4, isdim = 4*lelt)
319  integer isort(idim,isdim), ind(isdim), iwork(idim)
320  integer ipass, iseg, nseg
321  integer ninseg(lelt) ! elements in segment
322  logical ifseg(lelt) ! segment borders
323 
324  integer ibuf(2)
325 
326  ! functions
327  integer iglsum
328 
329 !#define DEBUG
330 #ifdef DEBUG
331  character*3 str1, str2
332  integer iunit
333  ! call number
334  integer icalld
335  save icalld
336  data icalld /0/
337 #endif
338 !-----------------------------------------------------------------------
339  ! find the mapping of the 3D projected elements on the 2D mesh and refine
340  tol = 5e-13
341  nt = lx1*ly1*lz1*lelt
342  npt_max = 256
343  nxf = 2*lx1
344  nyf = 2*ly1
345  nzf = 2*lz1
346  bb_t = 0.01
347 
348  ! start interpolation tool on given mesh
349  call fgslib_findpts_setup(ifpts,nekcomm,mp,ldim,xm1,ym1,zm1,
350  & lx1,ly1,lz1,nelt,nxf,nyf,nzf,bb_t,nt,nt,npt_max,tol)
351  ! identify elements
352  call fgslib_findpts(ifpts,rcode,1,proc,1,elid,1,rst,ldim,
353  & dist,1,pstat_cnt(1,1),ldim,pstat_cnt(2,1),ldim,
354  & pstat_cnt(ldim,1),ldim,pstat_nel)
355  ! close interpolation tool
356  call fgslib_findpts_free(ifpts)
357 
358  ! find problems with interpolation
359  nfail = 0
360  do il = 1,pstat_nel
361  ! check return code
362  if (rcode(il).eq.1) then
363  if (sqrt(dist(il)).gt.toldist) nfail = nfail + 1
364  elseif(rcode(il).eq.2) then
365  nfail = nfail + 1
366  endif
367  enddo
368  nfail = iglsum(nfail,1)
369 
370 #ifdef DEBUG
371  ! for testing
372  ! to output refinement
373  icalld = icalld+1
374  call io_file_freeid(iunit, ierr)
375  write(str1,'(i3.3)') nid
376  write(str2,'(i3.3)') icalld
377  open(unit=iunit,file='CRDfpts.txt'//str1//'i'//str2)
378 
379  write(iunit,*) pstat_nelg, pstat_nel, nfail
380  do il=1,pstat_nel
381  write(iunit,*) il, pstat_gnel(il), proc(il), elid(il)+1,
382  & rcode(il), dist(il), (rst(jl+(il-1)*ldim),jl=1,ldim)
383  enddo
384 
385  close(iunit)
386 #endif
387 
388  if (nfail.gt.0) call mntr_abort(pstat_id,
389  $ 'Elements not identified in pstat_mesh_map')
390 
391  ! sort data
392  ! pack it first
393  do il=1,pstat_nel
394  isort(1,il) = proc(il)
395  isort(2,il) = elid(il) + 1 ! go from c to fortran count
396  isort(3,il) = pstat_gnel(il)
397 #ifdef AMR
398  isort(4,il) = pstat_lev(il) ! max amr level value
399 #else
400  isort(4,il) = 0
401 #endif
402  enddo
403 
404  !tupple sort
405  ! mark no section boundaries
406  do il=1,pstat_nel
407  ifseg(il) = .false.
408  enddo
409  ! perform local sorting to identify unique set sorting by directions
410  ! first run => whole set is one segment
411  nseg = 1
412  ifseg(1) = .true.
413  ninseg(1) = pstat_nel
414  ! Multiple passes eliminates false positives
415  do ipass=1,2
416  do jl=1,2 ! Sort within each segment (proc, element nr)
417  il=1
418  do iseg=1,nseg
419  call ituple_sort(isort(1,il),idim,ninseg(iseg),jl,1,
420  $ ind,iwork) ! key = jl
421  il = il + ninseg(iseg)
422  enddo
423 
424  do il=2,pstat_nel
425  ! find segments borders
426  if (isort(jl,il).ne.isort(jl,il-1)) ifseg(il)=.true.
427  enddo
428 
429  ! Count up number of different segments
430  nseg = 0
431  do il=1,pstat_nel
432  if (ifseg(il)) then
433  nseg = nseg+1
434  ninseg(nseg) = 1
435  else
436  ninseg(nseg) = ninseg(nseg) + 1
437  endif
438  enddo
439  enddo ! jl=1,2
440  enddo ! ipass=1,2
441  ! sorting end
442 
443  ! remove multiplicities
444  jl = 1
445  do iseg=1,nseg
446  do kl=1,idim
447  isort(kl,iseg) = isort(kl,jl)
448  enddo
449 #ifdef AMR
450  do il=jl+1,jl + ninseg(iseg)
451  if(isort(4,iseg).lt.isort(4,il)) then
452  isort(3,iseg) = isort(3,il)
453  isort(4,iseg) = isort(4,il)
454  endif
455  enddo
456 #endif
457  jl = jl + ninseg(iseg)
458  enddo
459 
460 #ifdef DEBUG
461  ! for testing
462  ! to output refinement
463  call io_file_freeid(iunit, ierr)
464  write(str1,'(i3.3)') nid
465  write(str2,'(i3.3)') icalld
466  open(unit=iunit,file='CRDsort.txt'//str1//'i'//str2)
467 
468  write(iunit,*) nseg
469 
470  do il=1,nseg
471  write(iunit,*) il, ninseg(il), (isort(jl,il),jl=1,idim)
472  enddo
473 
474  close(iunit)
475 #endif
476 
477  ! transfer and sort data
478  call fgslib_crystal_ituple_transfer(cr_h,isort,idim,nseg,isdim,1)
479  il = 2
480  call fgslib_crystal_ituple_sort(cr_h,isort,idim,nseg,il,1)
481 
482  ! remove duplicates
483  iseg = 1
484  do il = 2, nseg
485  if (isort(2,iseg).eq.isort(2,il)) then
486 #ifdef AMR
487  if (isort(4,iseg).lt.isort(4,il)) then
488  isort(3,iseg) = isort(3,il)
489  isort(4,iseg) = isort(4,il)
490  endif
491 #endif
492  else
493  iseg = iseg + 1
494  if (iseg.ne.il) then
495  do kl=1,idim
496  isort(kl,iseg) = isort(kl,il)
497  enddo
498  endif
499  endif
500  enddo
501  nseg = iseg
502 
503 #ifdef DEBUG
504  ! for testing
505  ! to output refinement
506  call io_file_freeid(iunit, ierr)
507  write(str1,'(i3.3)') nid
508  write(str2,'(i3.3)') icalld
509  open(unit=iunit,file='CRDtrans.txt'//str1//'i'//str2)
510 
511  write(iunit,*) nseg, nelv
512 
513  do il=1,nseg
514  write(iunit,*) il, (isort(jl,il),jl=1,idim)
515  enddo
516 
517  close(iunit)
518 #endif
519 
520  ! check consistency
521  if (nseg.ne.nelt) then
522  ierr = 1
523  else
524  ierr = 0
525  endif
526  call mntr_check_abort(pstat_id,ierr,
527  $ 'Inconsistent segment number in pstat_mesh_map')
528 
529 #ifdef AMR
530  if (ifrref) then
531  ! transfer data for refinement mark
532  do il = 1, nseg
533  pstat_refl(isort(2,il)) = isort(4,il)
534  enddo
535 
536  ! for refinement return here
537  return
538  endif
539 #endif
540 
541  ! set global element number for transfer
542  do il = 1, nelt
543  pstat_gnel(isort(2,il)) = isort(3,il)
544  enddo
545 
546  ! notify element owners
547  do il = 1, nelt
548  isort(1,il) = gllnid(pstat_gnel(il))
549  isort(2,il) = gllel(pstat_gnel(il))
550  isort(3,il) = il
551  isort(4,il) = lglel(il)
552  enddo
553 
554  ! transfer and sort data
555  iseg = nelt
556  call fgslib_crystal_ituple_transfer(cr_h,isort,idim,iseg,isdim,1)
557  ! sanity check
558  ierr = 0
559  if (iseg.ne.nelt) ierr = 1 !!! test for good elements?
560  call mntr_check_abort(pstat_id,ierr,
561  $ 'Inconsistent element number in pstat_mesh_map')
562  il = 2
563  call fgslib_crystal_ituple_sort(cr_h,isort,idim,iseg,il,1)
564 
565 
566 #ifdef DEBUG
567  ! sanity check
568  ierr = 0
569  do il=1,iseg
570  if (isort(3,il).ne.gllel(isort(4,il))) ierr = 1
571  enddo
572  call mntr_check_abort(pstat_id,ierr,
573  $ 'Inconsistent element transfer number in pstat_mesh_map')
574  ierr = 0
575  do il=1,iseg
576  if (isort(1,il).ne.gllnid(isort(4,il))) ierr = 1
577  enddo
578  call mntr_check_abort(pstat_id,ierr,
579  $ 'Inconsistent process transfer number in pstat_mesh_map')
580 #endif
581 
582  ! save transfer data
583  do il=1,nelt
584  pstat_gnel(il) = isort(4,il)
585  enddo
586 
587 #ifdef DEBUG
588  ! for testing
589  ! to output refinement
590  call io_file_freeid(iunit, ierr)
591  write(str1,'(i3.3)') nid
592  write(str2,'(i3.3)') icalld
593  open(unit=iunit,file='CRDmap.txt'//str1//'i'//str2)
594 
595  write(iunit,*) nelgt,nelt
596  do il=1,nelt
597  write(iunit,*) il, pstat_gnel(il)
598  enddo
599 
600  close(iunit)
601 #endif
602 
603 #undef DEBUG
604 
605  return
606  end subroutine
607 !=======================================================================
610  subroutine pstat2d_transfer()
611  implicit none
612 
613  include 'SIZE'
614  include 'GEOM'
615  include 'SOLN'
616  include 'PARALLEL'
617  include 'FRAMELP'
618  include 'PSTAT2D'
619 
620  ! local variables
621  integer il, jl
622  integer lnelt, itmp
623  integer isw
624  parameter(isw=2)
625 
626  ! transfer arrays
627  integer vi(isw,LELT)
628  integer*8 vl(1) ! required by crystal rauter
629 
630  integer key(1) ! required by crystal rauter; for sorting
631 !-----------------------------------------------------------------------
632  ! element size
633  itmp = lx1*ly1*lz1
634 
635  ! coordinates
636  ! number of local elements
637  lnelt = nelt
638 
639  ! pack transfer data
640  do il=1,lnelt
641  vi(1,il) = gllel(pstat_gnel(il))
642  vi(2,il) = gllnid(pstat_gnel(il))
643  enddo
644 
645  ! transfer arrays
646  call fgslib_crystal_tuple_transfer
647  $ (cr_h,lnelt,lelt,vi,isw,vl,0,xm1,itmp,2)
648 
649  ! test local element number
650  if (lnelt.ne.nelt) then
651  call mntr_abort('Error: pstat_transfer; lnelt /= nelt')
652  endif
653 
654  ! sort elements acording to their global number
655  key(1) = 1
656  call fgslib_crystal_tuple_sort
657  $ (cr_h,lnelt,vi,isw,vl,0,xm1,itmp,key,1)
658 
659  ! number of local elements
660  lnelt = nelt
661 
662  ! pack transfer data
663  do il=1,lnelt
664  vi(1,il) = gllel(pstat_gnel(il))
665  vi(2,il) = gllnid(pstat_gnel(il))
666  enddo
667 
668  ! transfer arrays
669  call fgslib_crystal_tuple_transfer
670  $ (cr_h,lnelt,lelt,vi,isw,vl,0,ym1,itmp,2)
671 
672  ! test local element number
673  if (lnelt.ne.nelt) then
674  call mntr_abort('Error: pstat_transfer; lnelt /= nelt')
675  endif
676 
677  ! sort elements acording to their global number
678  key(1) = 1
679  call fgslib_crystal_tuple_sort
680  $ (cr_h,lnelt,vi,isw,vl,0,ym1,itmp,key,1)
681 
682  ! for each variable
683  do il=1, pstat_svar
684  ! number of local elements
685  lnelt = nelt
686 
687  ! pack transfer data
688  do jl=1,lnelt
689  vi(1,jl) = gllel(pstat_gnel(jl))
690  vi(2,jl) = gllnid(pstat_gnel(jl))
691  enddo
692 
693  ! transfer array
694  call fgslib_crystal_tuple_transfer
695  $ (cr_h,lnelt,lelt,vi,isw,vl,0,t(1,1,1,1,il+1),itmp,2)
696 
697  ! test local element number
698  if (lnelt.ne.nelt) then
699  call mntr_abort('Error: pstat_transfer; lnelt /= nelt')
700  endif
701 
702  ! sort elements acording to their global number
703  key(1) = 1
704  call fgslib_crystal_tuple_sort
705  $ (cr_h,lnelt,vi,isw,vl,0,t(1,1,1,1,il+1),itmp,key,1)
706  enddo
707 
708  return
709  end subroutine
710 !=======================================================================
713  subroutine pstat2d_sts_avg
714  implicit none
715 
716  include 'SIZE'
717  include 'INPUT'
718  include 'RESTART'
719  include 'TSTEP'
720  include 'SOLN'
721  include 'FRAMELP'
722  include 'PSTAT2D'
723 
724  ! local variables
725  integer il, jl ! loop index
726  integer ierr ! error mark
727  integer nvec ! single field length
728  character*3 prefix ! file prefix
729  character*132 fname ! file name
730  character*132 bname ! base name
731  character*6 str ! file number
732  integer nps1,nps0,npsr ! number of passive scalars in the file
733  integer istepr ! number of time step in restart files
734  real ltime, dtime ! simulation time and time update
735  real rtmp
736 
737 !-----------------------------------------------------------------------
738  ! no regular mesh; important for file name generation
739  ifreguo = .false.
740 
741  call io_init
742 
743  ierr=0
744  ! open files on i/o nodes
745  prefix='sts'
746  ! get base name (SESSION)
747  bname = trim(adjustl(session))
748 
749  ! mark variables to be read
750  ifgetx=.true.
751  ifgetz=.false.
752  ifgetu=.false.
753  ifgetw=.false.
754  ifgetp=.false.
755  ifgett=.false.
756  do il=1,pstat_svar
757  ifgtps(il)=.true.
758  enddo
759 
760  ! initial time and step count
761  ltime = pstat_stime
762  istepr = 0
763 
764  ! initilise vectors
765  call rzero(pstat_ruavg,lx1**ldim*lelt*pstat_svar)
766  nvec = lx1*ly1*nelt
767 
768  ! loop over stat files
769  do il = pstat_ffile, pstat_nfile
770 
771  call mntr_logi(pstat_id,lp_inf,'Opening sts file nr=',il)
772 
773  ! read cenre data to generate proper sorting
774  call pstat2d_mfi_crd2d(il)
775  ! get element mapping
776  call pstat2d_mesh_map(.false.)
777 
778  call io_mfo_fname(fname,bname,prefix,ierr)
779  write(str,'(i5.5)') il
780  fname = trim(fname)//trim(str)
781 
782  fid0 = 0
783  call addfid(fname,fid0)
784 
785  ! add directory name
786  fname = 'DATA/'//trim(fname)
787 
788  !call load_fld(fname)
789  call mfi(fname,il)
790 
791  ! extract number of passive scalars in the file and check consistency
792  do jl=1,10
793  if (rdcode1(jl).eq.'S') then
794  read(rdcode1(jl+1),'(i1)') nps1
795  read(rdcode1(jl+2),'(i1)') nps0
796  npsr = 10*nps1+nps0
797  endif
798  enddo
799  if (npsr.ne.pstat_svar) call mntr_abort(pstat_id,
800  $ 'Inconsistent number of fielsd in sts file.')
801 
802  ! calculate interval and update time
803  dtime = timer - ltime
804  ltime = timer
805 
806  ! sum number of time steps
807  istepr = istepr + istpr
808 
809  ! reshuffle lelements
810  call pstat2d_transfer()
811 
812  ! accumulate fileds
813  do jl = 1,pstat_svar
814 ! call add2s2(pstat_ruavg(1,1,pstat_swfield(jl)),
815 ! $ t(1,1,1,1,jl+1),dtime,nvec)
816  call add2s2(pstat_ruavg(1,1,jl),t(1,1,1,1,jl+1),dtime,nvec)
817  enddo
818 
819  enddo
820 
821  ! save data for file header
822  pstat_etime = ltime
823  pstat_istepr = istepr
824 
825  ! divide by time span
826  if (ltime.ne.pstat_stime) then
827  rtmp = 1.0/(ltime-pstat_stime)
828  else
829  rtmp = 1.0
830  endif
831 
832  do il = 1,pstat_svar
833  call cmult(pstat_ruavg(1,1,il),rtmp,nvec)
834  enddo
835 
836  ! save all averaged fields
837  ifvo = .false.
838  ifpo = .false.
839  ifto = .false.
840  do il=1,pstat_svar
841  ifpsco(il)=.true.
842  enddo
843  ! put variables back to temperature
844  do il=1,pstat_svar
845  call copy(t(1,1,1,1,il+1),pstat_ruavg(1,1,il),nvec)
846  enddo
847  call outpost2(vx,vy,vz,pr,t,pstat_svar+1,'st1')
848  ! sort variables according to pstat_swfield
849  do il = 1,pstat_svar
850  call copy(pstat_ruavg(1,1,pstat_swfield(il)),t(1,1,1,1,il+1),
851  $ nvec)
852  end do
853 
854  return
855  end subroutine
856 !=======================================================================
859  subroutine pstat2d_deriv
860  implicit none
861 
862  include 'SIZE'
863  include 'INPUT'
864  include 'SOLN'
865  include 'FRAMELP'
866  include 'PSTAT2D'
867 
868  ! local variables
869  integer il ! loop index
870  integer nvec ! single field length
871  real dudx(lx1*ly1*lz1,lelt,3) ! field derivatives
872 !-----------------------------------------------------------------------
873  ! update derivative arrays
874  call geom_reset(1)
875 
876  ! initilise vectors
877  call rzero(pstat_ruder,lx1**ldim*lelt*pstat_dvar)
878  nvec = lx1*ly1*nelt
879 
880  ! Tensor 1. Compute 2D derivative tensor with U
881  call gradm1(dudx(1,1,1),dudx(1,1,2),dudx(1,1,3),
882  $ pstat_ruavg(1,1,1))
883  call copy(pstat_ruder(1,1,1),dudx(1,1,1),nvec) ! dU/dx
884  call copy(pstat_ruder(1,1,2),dudx(1,1,2),nvec) ! dU/dy
885 
886  ! Tensor 1. Compute 2D derivative tensor with V
887  call gradm1(dudx(1,1,1),dudx(1,1,2),dudx(1,1,3),
888  $ pstat_ruavg(1,1,2))
889  call copy(pstat_ruder(1,1,3),dudx(1,1,1),nvec) ! dV/dx
890  call copy(pstat_ruder(1,1,4),dudx(1,1,2),nvec) ! dV/dy
891 
892  ! Tensor 2. Compute 2D derivative tensor with W
893  call gradm1(dudx(1,1,1),dudx(1,1,2),dudx(1,1,3),
894  $ pstat_ruavg(1,1,3))
895  call copy(pstat_ruder(1,1,5),dudx(1,1,1),nvec) ! dW/dx
896  call copy(pstat_ruder(1,1,6),dudx(1,1,2),nvec) ! dW/dy
897 
898  ! Tensor 2. Compute 2D derivative tensor with P
899  call gradm1(dudx(1,1,1),dudx(1,1,2),dudx(1,1,3),
900  $ pstat_ruavg(1,1,4))
901  call copy(pstat_ruder(1,1,7),dudx(1,1,1),nvec) ! dP/dx
902  call copy(pstat_ruder(1,1,8),dudx(1,1,2),nvec) ! dP/dy
903 
904  ! Tensor 3. Compute 2D derivative tensor with <uu>
905  call gradm1(dudx(1,1,1),dudx(1,1,2),dudx(1,1,3),
906  $ pstat_ruavg(1,1,5))
907  call copy(pstat_ruder(1,1,9),dudx(1,1,1),nvec) ! d<uu>/dx
908  call copy(pstat_ruder(1,1,10),dudx(1,1,2),nvec) ! d<uu>/dy
909 
910  ! Tensor 3. Compute 2D derivative tensor with <vv>
911  call gradm1(dudx(1,1,1),dudx(1,1,2),dudx(1,1,3),
912  $ pstat_ruavg(1,1,6))
913  call copy(pstat_ruder(1,1,11),dudx(1,1,1),nvec) ! d<vv>/dx
914  call copy(pstat_ruder(1,1,12),dudx(1,1,2),nvec) ! d<vv>/dy
915 
916  ! Tensor 4. Compute 2D derivative tensor with <ww>
917  call gradm1(dudx(1,1,1),dudx(1,1,2),dudx(1,1,3),
918  $ pstat_ruavg(1,1,7))
919  call copy(pstat_ruder(1,1,13),dudx(1,1,1),nvec) ! d<ww>/dx
920  call copy(pstat_ruder(1,1,14),dudx(1,1,2),nvec) ! d<ww>/dy
921 
922  ! Tensor 4. Compute 2D derivative tensor with <pp>
923  call gradm1(dudx(1,1,1),dudx(1,1,2),dudx(1,1,3),
924  $ pstat_ruavg(1,1,8))
925  call copy(pstat_ruder(1,1,15),dudx(1,1,1),nvec) ! d<pp>/dx
926  call copy(pstat_ruder(1,1,16),dudx(1,1,2),nvec) ! d<pp>/dy
927 
928  ! Tensor 5. Compute 2D derivative tensor with <uv>
929  call gradm1(dudx(1,1,1),dudx(1,1,2),dudx(1,1,3),
930  $ pstat_ruavg(1,1,9))
931  call copy(pstat_ruder(1,1,17),dudx(1,1,1),nvec) ! d<uv>/dx
932  call copy(pstat_ruder(1,1,18),dudx(1,1,2),nvec) ! d<uv>/dy
933 
934  ! Tensor 5. Compute 2D derivative tensor with <vw>
935  call gradm1(dudx(1,1,1),dudx(1,1,2),dudx(1,1,3),
936  $ pstat_ruavg(1,1,10))
937  call copy(pstat_ruder(1,1,19),dudx(1,1,1),nvec) ! d<vw>/dx
938  call copy(pstat_ruder(1,1,20),dudx(1,1,2),nvec) ! d<vw>/dy
939 
940  ! Tensor 6. Compute 2D derivative tensor with <uw>
941  call gradm1(dudx(1,1,1),dudx(1,1,2),dudx(1,1,3),
942  $ pstat_ruavg(1,1,11))
943  call copy(pstat_ruder(1,1,21),dudx(1,1,1),nvec) ! d<uw>/dx
944  call copy(pstat_ruder(1,1,22),dudx(1,1,2),nvec) ! d<uw>/dy
945 
946  ! Tensor 6. Compute 2D derivative tensor with <uuu>
947  call gradm1(dudx(1,1,1),dudx(1,1,2),dudx(1,1,3),
948  $ pstat_ruavg(1,1,24))
949  call copy(pstat_ruder(1,1,23),dudx(1,1,1),nvec) ! d<uuu>/dx
950  call copy(pstat_ruder(1,1,24),dudx(1,1,2),nvec) ! d<uuu>/dy
951 
952  ! Tensor 7. Compute 2D derivative tensor with <vvv>
953  call gradm1(dudx(1,1,1),dudx(1,1,2),dudx(1,1,3),
954  $ pstat_ruavg(1,1,25))
955  call copy(pstat_ruder(1,1,25),dudx(1,1,1),nvec) ! d<vvv>/dx
956  call copy(pstat_ruder(1,1,26),dudx(1,1,2),nvec) ! d<vvv>/dy
957 
958  ! Tensor 7. Compute 2D derivative tensor with <www>
959  call gradm1(dudx(1,1,1),dudx(1,1,2),dudx(1,1,3),
960  $ pstat_ruavg(1,1,26))
961  call copy(pstat_ruder(1,1,27),dudx(1,1,1),nvec) ! d<www>/dx
962  call copy(pstat_ruder(1,1,28),dudx(1,1,2),nvec) ! d<www>/dy
963 
964  ! Tensor 8. Compute 2D derivative tensor with <ppp>
965  call gradm1(dudx(1,1,1),dudx(1,1,2),dudx(1,1,3),
966  $ pstat_ruavg(1,1,27))
967  call copy(pstat_ruder(1,1,29),dudx(1,1,1),nvec) ! d<ppp>/dx
968  call copy(pstat_ruder(1,1,30),dudx(1,1,2),nvec) ! d<ppp>/dy
969 
970  ! Tensor 8. Compute 2D derivative tensor with <uuv>
971  call gradm1(dudx(1,1,1),dudx(1,1,2),dudx(1,1,3),
972  $ pstat_ruavg(1,1,28))
973  call copy(pstat_ruder(1,1,31),dudx(1,1,1),nvec) ! d<uuv>/dx
974  call copy(pstat_ruder(1,1,32),dudx(1,1,2),nvec) ! d<uuv>/dy
975 
976  ! Tensor 9. Compute 2D derivative tensor with <uuw>
977  call gradm1(dudx(1,1,1),dudx(1,1,2),dudx(1,1,3),
978  $ pstat_ruavg(1,1,29))
979  call copy(pstat_ruder(1,1,33),dudx(1,1,1),nvec) ! d<uuw>/dx
980  call copy(pstat_ruder(1,1,34),dudx(1,1,2),nvec) ! d<uuw>/dy
981 
982  ! Tensor 9. Compute 2D derivative tensor with <vvu>
983  call gradm1(dudx(1,1,1),dudx(1,1,2),dudx(1,1,3),
984  $ pstat_ruavg(1,1,30))
985  call copy(pstat_ruder(1,1,35),dudx(1,1,1),nvec) ! d<vvu>/dx
986  call copy(pstat_ruder(1,1,36),dudx(1,1,2),nvec) ! d<vvu>/dy
987 
988  ! Tensor 10. Compute 2D derivative tensor with <vvw>
989  call gradm1(dudx(1,1,1),dudx(1,1,2),dudx(1,1,3),
990  $ pstat_ruavg(1,1,31))
991  call copy(pstat_ruder(1,1,37),dudx(1,1,1),nvec) ! d<vvw>/dx
992  call copy(pstat_ruder(1,1,38),dudx(1,1,2),nvec) ! d<vvw>/dy
993 
994  ! Tensor 10. Compute 2D derivative tensor with <wwu>
995  call gradm1(dudx(1,1,1),dudx(1,1,2),dudx(1,1,3),
996  $ pstat_ruavg(1,1,32))
997  call copy(pstat_ruder(1,1,39),dudx(1,1,1),nvec) ! d<wwu>/dx
998  call copy(pstat_ruder(1,1,40),dudx(1,1,2),nvec) ! d<wwu>/dy
999 
1000  ! Tensor 11. Compute 2D derivative tensor with <wwv>
1001  call gradm1(dudx(1,1,1),dudx(1,1,2),dudx(1,1,3),
1002  $ pstat_ruavg(1,1,33))
1003  call copy(pstat_ruder(1,1,41),dudx(1,1,1),nvec) ! d<wwu>/dx
1004  call copy(pstat_ruder(1,1,42),dudx(1,1,2),nvec) ! d<wwu>/dy
1005 
1006  ! Tensor 11. Compute 2D derivative tensor with <uvw>
1007  call gradm1(dudx(1,1,1),dudx(1,1,2),dudx(1,1,3),
1008  $ pstat_ruavg(1,1,34))
1009  call copy(pstat_ruder(1,1,43),dudx(1,1,1),nvec) ! d<uvw>/dx
1010  call copy(pstat_ruder(1,1,44),dudx(1,1,2),nvec) ! d<uvw>/dy
1011 
1012  ! Tensor 12. Compute 2D derivative tensor with dU/dx
1013  call gradm1(dudx(1,1,1),dudx(1,1,2),dudx(1,1,3),
1014  $ pstat_ruder(1,1,1))
1015  call copy(pstat_ruder(1,1,45),dudx(1,1,1),nvec) ! d2U/dx2
1016 
1017  ! Tensor 12. Compute 2D derivative tensor with dU/dy
1018  call gradm1(dudx(1,1,1),dudx(1,1,2),dudx(1,1,3),
1019  $ pstat_ruder(1,1,2))
1020  call copy(pstat_ruder(1,1,46),dudx(1,1,2),nvec) ! d2U/dy2
1021 
1022  ! Tensor 13. Compute 2D derivative tensor with dV/dx
1023  call gradm1(dudx(1,1,1),dudx(1,1,2),dudx(1,1,3),
1024  $ pstat_ruder(1,1,3))
1025  call copy(pstat_ruder(1,1,47),dudx(1,1,1),nvec) ! d2V/dx2
1026 
1027  ! Tensor 13. Compute 2D derivative tensor with dV/dy
1028  call gradm1(dudx(1,1,1),dudx(1,1,2),dudx(1,1,3),
1029  $ pstat_ruder(1,1,4))
1030  call copy(pstat_ruder(1,1,48),dudx(1,1,2),nvec) ! d2V/dy2
1031 
1032  ! Tensor 14. Compute 2D derivative tensor with dW/dx
1033  call gradm1(dudx(1,1,1),dudx(1,1,2),dudx(1,1,3),
1034  $ pstat_ruder(1,1,5))
1035  call copy(pstat_ruder(1,1,49),dudx(1,1,1),nvec) ! d2W/dx2
1036 
1037  ! Tensor 14. Compute 2D derivative tensor with dW/dy
1038  call gradm1(dudx(1,1,1),dudx(1,1,2),dudx(1,1,3),
1039  $ pstat_ruder(1,1,6))
1040  call copy(pstat_ruder(1,1,50),dudx(1,1,2),nvec) ! d2W/dy2
1041 
1042  ! Tensor 15. Compute 2D derivative tensor with d<uu>/dx
1043  call gradm1(dudx(1,1,1),dudx(1,1,2),dudx(1,1,3),
1044  $ pstat_ruder(1,1,9))
1045  call copy(pstat_ruder(1,1,51),dudx(1,1,1),nvec) ! d2<uu>/dx2
1046 
1047  ! Tensor 15. Compute 2D derivative tensor with d<uu>/dy
1048  call gradm1(dudx(1,1,1),dudx(1,1,2),dudx(1,1,3),
1049  $ pstat_ruder(1,1,10))
1050  call copy(pstat_ruder(1,1,52),dudx(1,1,2),nvec) ! d2<uu>/dy2
1051 
1052  ! Tensor 16. Compute 2D derivative tensor with d<vv>/dx
1053  call gradm1(dudx(1,1,1),dudx(1,1,2),dudx(1,1,3),
1054  $ pstat_ruder(1,1,11))
1055  call copy(pstat_ruder(1,1,53),dudx(1,1,1),nvec) ! d2<vv>/dx2
1056 
1057  ! Tensor 16. Compute 2D derivative tensor with d<vv>/dy
1058  call gradm1(dudx(1,1,1),dudx(1,1,2),dudx(1,1,3),
1059  $ pstat_ruder(1,1,12))
1060  call copy(pstat_ruder(1,1,54),dudx(1,1,2),nvec) ! d2<vv>/dy2
1061 
1062  ! Tensor 17. Compute 2D derivative tensor with d<ww>/dx
1063  call gradm1(dudx(1,1,1),dudx(1,1,2),dudx(1,1,3),
1064  $ pstat_ruder(1,1,13))
1065  call copy(pstat_ruder(1,1,55),dudx(1,1,1),nvec) ! d2<ww>/dx2
1066 
1067  ! Tensor 17. Compute 2D derivative tensor with d<ww>/dy
1068  call gradm1(dudx(1,1,1),dudx(1,1,2),dudx(1,1,3),
1069  $ pstat_ruder(1,1,14))
1070  call copy(pstat_ruder(1,1,56),dudx(1,1,2),nvec) ! d2<ww>/dy2
1071 
1072  ! Tensor 18. Compute 2D derivative tensor with d<uv>/dx
1073  call gradm1(dudx(1,1,1),dudx(1,1,2),dudx(1,1,3),
1074  $ pstat_ruder(1,1,17))
1075  call copy(pstat_ruder(1,1,57),dudx(1,1,1),nvec) ! d2<uv>/dx2
1076 
1077  ! Tensor 18. Compute 2D derivative tensor with d<uv>/dy
1078  call gradm1(dudx(1,1,1),dudx(1,1,2),dudx(1,1,3),
1079  $ pstat_ruder(1,1,18))
1080  call copy(pstat_ruder(1,1,58),dudx(1,1,2),nvec) ! d2<uv>/dy2
1081 
1082  ! Tensor 19. Compute 2D derivative tensor with d<uw>/dx
1083  call gradm1(dudx(1,1,1),dudx(1,1,2),dudx(1,1,3),
1084  $ pstat_ruder(1,1,21))
1085  call copy(pstat_ruder(1,1,59),dudx(1,1,1),nvec) ! d2<uw>/dx2
1086 
1087  ! Tensor 19. Compute 2D derivative tensor with d<uw>/dy
1088  call gradm1(dudx(1,1,1),dudx(1,1,2),dudx(1,1,3),
1089  $ pstat_ruder(1,1,22))
1090  call copy(pstat_ruder(1,1,60),dudx(1,1,2),nvec) ! d2<uw>/dy2
1091 
1092  ! Tensor 20. Compute 2D derivative tensor with d<vw>/dx
1093  call gradm1(dudx(1,1,1),dudx(1,1,2),dudx(1,1,3),
1094  $ pstat_ruder(1,1,19))
1095  call copy(pstat_ruder(1,1,61),dudx(1,1,1),nvec) ! d2<vw>/dx2
1096 
1097  ! Tensor 20. Compute 2D derivative tensor with d<vw>/dy
1098  call gradm1(dudx(1,1,1),dudx(1,1,2),dudx(1,1,3),
1099  $ pstat_ruder(1,1,20))
1100  call copy(pstat_ruder(1,1,62),dudx(1,1,2),nvec) ! d2<vw>/dy2
1101 
1102  ! Tensor 21. Compute 2D derivative tensor with <pu>
1103  call gradm1(dudx(1,1,1),dudx(1,1,2),dudx(1,1,3),
1104  $ pstat_ruavg(1,1,12))
1105  call copy(pstat_ruder(1,1,63),dudx(1,1,1),nvec) ! d<pu>/dx
1106  call copy(pstat_ruder(1,1,64),dudx(1,1,2),nvec) ! d<pu>/dy
1107 
1108  ! Tensor 21. Compute 2D derivative tensor with <pv>
1109  call gradm1(dudx(1,1,1),dudx(1,1,2),dudx(1,1,3),
1110  $ pstat_ruavg(1,1,13))
1111  call copy(pstat_ruder(1,1,65),dudx(1,1,1),nvec) ! d<pv>/dx
1112  call copy(pstat_ruder(1,1,66),dudx(1,1,2),nvec) ! d<pv>/dy
1113 
1114  ! Tensor 22. Compute 2D derivative tensor with <pw>
1115  call gradm1(dudx(1,1,1),dudx(1,1,2),dudx(1,1,3),
1116  $ pstat_ruavg(1,1,14))
1117  call copy(pstat_ruder(1,1,67),dudx(1,1,1),nvec) ! d<pw>/dx
1118  call copy(pstat_ruder(1,1,68),dudx(1,1,2),nvec) ! d<pw>/dy
1119 
1120  ! save all derivatives
1121  ifvo = .false.
1122  ifpo = .false.
1123  ifto = .false.
1124  do il=1,pstat_dvar
1125  ifpsco(il)=.true.
1126  enddo
1127  ! put variables back to temperature
1128  do il=1,pstat_dvar
1129  call copy(t(1,1,1,1,il+1),pstat_ruder(1,1,il),nvec)
1130  enddo
1131  call outpost2(vx,vy,vz,pr,t,pstat_dvar+1,'st2')
1132 
1133  return
1134  end subroutine
1135 !=======================================================================
1138  subroutine pstat2d_interp
1139  implicit none
1140 
1141  include 'SIZE'
1142  include 'GEOM'
1143  include 'FRAMELP'
1144  include 'PSTAT2D'
1145 
1146  ! global data structures
1147  integer mid,mp,nekcomm,nekgroup,nekreal
1148  common /nekmpi/ mid,mp,nekcomm,nekgroup,nekreal
1149 
1150  ! local variables
1151  integer ifpts ! findpts flag
1152  real tol ! interpolation tolerance
1153  integer nt
1154  integer npt_max ! max communication set
1155  integer nxf, nyf, nzf ! fine mesh for bb-test
1156  real bb_t ! relative size to expand bounding boxes by
1157  real toldist
1158  parameter(toldist = 5e-6)
1159 
1160  integer rcode(lhis), proc(lhis),elid(lhis)
1161  real dist(lhis), rst(ldim*lhis)
1162 
1163  integer nfail, npass
1164 
1165  integer il ! loop index
1166  !integer nvec ! single field length
1167 
1168  ! functions
1169  integer iglsum
1170 
1171 !#define DEBUG
1172 #ifdef DEBUG
1173  character*3 str1, str2
1174  integer iunit, ierr, jl
1175  ! call number
1176  integer icalld
1177  save icalld
1178  data icalld /0/
1179 #endif
1180 !-----------------------------------------------------------------------
1181  ! read point position
1182  call pstat2d_mfi_interp
1183 
1184  ! initialise interpolation tool
1185  tol = 5e-13
1186  nt = lx1*ly1*lz1*lelt
1187  npt_max = 256
1188  nxf = 2*lx1
1189  nyf = 2*ly1
1190  nzf = 2*lz1
1191  bb_t = 0.01
1192 
1193  ! start interpolation tool on given mesh
1194  call fgslib_findpts_setup(ifpts,nekcomm,mp,ldim,xm1,ym1,zm1,
1195  & lx1,ly1,lz1,nelt,nxf,nyf,nzf,bb_t,nt,nt,npt_max,tol)
1196 
1197  ! identify points
1198  call fgslib_findpts(ifpts,rcode,1,proc,1,elid,1,rst,ldim,dist,1,
1199  & pstat_int_pts(1,1),ldim,pstat_int_pts(2,1),ldim,
1200  & pstat_int_pts(ldim,1),ldim,pstat_npt)
1201 
1202  ! find problems with interpolation
1203  nfail = 0
1204  do il = 1,pstat_npt
1205  ! check return code
1206  if (rcode(il).eq.1) then
1207  if (sqrt(dist(il)).gt.toldist) nfail = nfail + 1
1208  elseif(rcode(il).eq.2) then
1209  nfail = nfail + 1
1210  endif
1211  enddo
1212  nfail = iglsum(nfail,1)
1213 
1214 #ifdef DEBUG
1215  ! for testing
1216  ! to output refinement
1217  icalld = icalld+1
1218  call io_file_freeid(iunit, ierr)
1219  write(str1,'(i3.3)') nid
1220  write(str2,'(i3.3)') icalld
1221  open(unit=iunit,file='INTfpts.txt'//str1//'i'//str2)
1222 
1223  write(iunit,*) pstat_nptot, pstat_npt, nfail
1224  do il=1,pstat_npt
1225  write(iunit,*) il, proc(il), elid(il), rcode(il), dist(il),
1226  $ (rst(jl+(il-1)*ldim),jl=1,ldim)
1227  enddo
1228 
1229  close(iunit)
1230 #endif
1231 
1232  if (nfail.gt.0) call mntr_abort(pstat_id,
1233  $ 'pstat_interp: Points not mapped')
1234 
1235  ! Interpolate averaged fields
1236  do il=1,pstat_svar
1237  call fgslib_findpts_eval(ifpts,pstat_int_avg(1,il),1,
1238  & rcode,1,proc,1,elid,1,rst,ndim,pstat_npt,
1239  & pstat_ruavg(1,1,il))
1240  enddo
1241 
1242 #ifdef DEBUG
1243  ! for testing
1244  ! to output refinement
1245  call io_file_freeid(iunit, ierr)
1246  write(str1,'(i3.3)') nid
1247  write(str2,'(i3.3)') icalld
1248  open(unit=iunit,file='INTavg.txt'//str1//'i'//str2)
1249 
1250  write(iunit,*) pstat_nptot, pstat_npt
1251  do il=1,pstat_npt
1252  write(iunit,*) il, (pstat_int_avg(il,jl),jl=1,4)
1253  enddo
1254 
1255  close(iunit)
1256 #endif
1257 
1258  ! Interpolate fields derivatives
1259  do il=1,pstat_dvar
1260  call fgslib_findpts_eval(ifpts,pstat_int_der(1,il),1,
1261  & rcode,1,proc,1,elid,1,rst,ndim,pstat_npt,
1262  & pstat_ruder(1,1,il))
1263  enddo
1264 
1265 #ifdef DEBUG
1266  ! for testing
1267  ! to output refinement
1268  call io_file_freeid(iunit, ierr)
1269  write(str1,'(i3.3)') nid
1270  write(str2,'(i3.3)') icalld
1271  open(unit=iunit,file='INTder.txt'//str1//'i'//str2)
1272 
1273  write(iunit,*) pstat_nptot, pstat_npt
1274  do il=1,pstat_npt
1275  write(iunit,*) il, (pstat_int_der(il,jl),jl=1,4)
1276  enddo
1277 
1278  close(iunit)
1279 #endif
1280 
1281  ! finalise interpolation tool
1282  call fgslib_findpts_free(ifpts)
1283 
1284  ! write down interpolated values
1285  call pstat2d_mfo_interp
1286 
1287 #undef DEBUG
1288  return
1289  end subroutine
1290 !=======================================================================
gllel
integer function gllel(ieg)
Definition: dprocmap.f:183
gllnid
integer function gllnid(ieg)
Definition: dprocmap.f:161
io_file_freeid
subroutine io_file_freeid(iunit, ierr)
Get free file unit number and store max unit value.
Definition: io_tools.f:47
io_mfo_fname
subroutine io_mfo_fname(fname, bname, prefix, ierr)
Generate file name according to nek rulles without opening the file.
Definition: io_tools.f:109
mntr_logi
subroutine mntr_logi(mid, priority, logs, ivar)
Write log message adding single integer.
Definition: mntrlog.f:709
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
mntr_check_abort
subroutine mntr_check_abort(mid, ierr, logs)
Abort simulation.
Definition: mntrlog.f:856
pstat2d_init
subroutine pstat2d_init()
Initilise pstat module.
Definition: pstat2D.f:100
pstat2d_is_initialised
logical function pstat2d_is_initialised()
Check if module was initialised.
Definition: pstat2D.f:175
pstat2d_mfi_crd2d
subroutine pstat2d_mfi_crd2d(fidl)
Read nonconforming data from the file.
Definition: pstat2D_IO.f:11
pstat2d_interp
subroutine pstat2d_interp
Interpolate int the set of points.
Definition: pstat2D.f:1139
pstat2d_transfer
subroutine pstat2d_transfer()
Reshuffle elements between sts and current ordering.
Definition: pstat2D.f:611
pstat2d_deriv
subroutine pstat2d_deriv
Calculate derivatives.
Definition: pstat2D.f:860
pstat2d_mfo_interp
subroutine pstat2d_mfo_interp
Geather data and write it down.
Definition: pstat2D_IO.f:392
pstat2d_mesh_amr
subroutine pstat2d_mesh_amr
Manipulate mesh to get proper AMR refinement.
Definition: pstat2D.f:230
pstat2d_register
subroutine pstat2d_register()
Register post processing statistics module.
Definition: pstat2D.f:11
pstat2d_main
subroutine pstat2d_main
Main interface of pstat module.
Definition: pstat2D.f:188
pstat2d_sts_avg
subroutine pstat2d_sts_avg
Read in fields and average them.
Definition: pstat2D.f:714
pstat2d_mesh_map
subroutine pstat2d_mesh_map(ifrref)
Find mapping of sts mesh to the existing in the run.
Definition: pstat2D.f:283
pstat2d_mfi_interp
subroutine pstat2d_mfi_interp
Read interpolation points position and redistribute them.
Definition: pstat2D_IO.f:199
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
mfi
subroutine mfi(fname_in, ifile)
Definition: ic.f:2355
addfid
subroutine addfid(fname, fid)
Definition: ic.f:2512
geom_reset
subroutine geom_reset(icall)
Definition: ic.f:1802
add2s2
subroutine add2s2(a, b, c1, n)
Definition: math.f:690
copy
subroutine copy(a, b, n)
Definition: math.f:260
cmult
subroutine cmult(a, const, n)
Definition: math.f:315
rzero
subroutine rzero(a, n)
Definition: math.f:208
gradm1
subroutine gradm1(ux, uy, uz, u)
Definition: navier5.f:463
ituple_sort
subroutine ituple_sort(a, lda, n, key, nkey, ind, aa)
Definition: navier8.f:327
outpost2
subroutine outpost2(v1, v2, v3, vp, vt, nfldt, name3)
Definition: prepost.f:1394
io_init
subroutine io_init
Definition: prepost.f:1024