subroutine help

[Akshay Patil]

Dear Nek-Nerds,

I was wondering if someone could help clarify what exactly the variables in the given subroutines mean and how I can use them?

planar_average_z(ua,u,w1,w2)

Here ua=planar average that the function returns

u = velocity component/parameter to be averaged (lets say I set it to vx)

however, I dont know what w1 and w2 mean? 

Additionally, when I need to define ua should be it something like ua(lz1*lelt) and/or do I need to define what ua is beforehand?

planar_average_s(uavg_pl,u,w1,w2)

Is this any different from the above subroutine?

avg3(avg,f,g, alpha,beta,n,name,ifverbose)

avg = the functions returns the result

f = the parameter that needs averaging?

g = can I leave this blank or do I need to specify this?

what is alpha, beta, n?

Another issue I had was when I use the code below and try to print it, it returns the results the number of mpi ranks active. How can I avoid doing this? (Please forgive my novice Fortran skills)

subroutine compute_ubar

      include 'SIZE'

      include 'TOTAL'

      include 'NEKUSE'

c     Initialize the parameters

      real ubar, vbar, wbar, U_mean

c     Total number of grid points

      n = nx1*ny1*nz1*nelv

c     Computing the mean velocity in the domain

      ubar = glsc2(vx,bm1,n)/volvm1

      vbar = glsc2(vy,bm1,n)/volvm1

      wbar = glsc2(vz,bm1,n)/volvm1

      U_mean = (ubar**2+vbar**2+wbar**2)**0.5

c     Print it to the screen

c If I comment out the if statement it prints out the results n times if I am using n mpi ranks.

c      if(nio.eq.0) then

       write (*,*) U_mean

c      endif

      return

      end

--

Regards,

Akshay Patil

[Tanmoy Chatterjee]

Hi Akshay,

Please see below.

For compute_ubar(), use if(nid .eq. 0) write(6,*) U_mean, 'mean'

avg3(avg,f,g, alpha,beta,n,name,ifverbose) 

You have to specify g, you cannot keep it blank. If you want to find the mean/average velocity you can use avg1() routines. Look for navier5.f 

source code if you want to create a cheat sheet for avg*() routines. 

planar_average_z(ua,u,w1,w2) and planar_average_s(ua,u,w1,w2) are farely similar. And if you want to call these planar average routines

you have to allocate memory to ua. w1's and w2's are work arrays. Its mode is input output. You can just set them as blank arrays. But you have to define them apriori.

w1 and w2 are same array size as vx,vy, etc.

In general, if you are also interested in looking at how these routines work, have a look at navier5.f in the source code and search for these subroutines. 

Best Regards,

Tanmoy

--
You received this message because you are subscribed to the Google Groups "Nek5000" group.
To unsubscribe from this group and stop receiving emails from it, send an email to nek5000+u...@googlegroups.com.
To view this discussion on the web visit https://groups.google.com/d/msgid/nek5000/CACqPjWG-3AfEDFSt_jpTHeqeGhc1XcGTPaQXzPGZ_O4ts1UqrA%40mail.gmail.com.

[Akshay Patil]

Hi Tanmoy,

Just to confirm, when using this if(nid .eq. 0) write(6,*) U_mean, 'mean', it does not just print the average value of U_mean on processor id 0 right? I was using this and got a bit confused with this.

With regards to the second part. I still dont seem to get this to work. I am trying to get planar average (plane x-y and z is vertical in my case) of the velocity component vx.

When i write this subroutine as follows and call it it throws me Program received signal SIGFPE: Floating-point exception - erroneous arithmetic operation. error.

      subroutine myZAverage

      include 'SIZE'

      include 'TOTAL'

      real ua(nz1),w1(nz1),w2(nz1)

      call planar_average_z(ua,vx,w1,w2)

      return

      end

When i use this below routine it says nelz cannot show up here error. So I am a little confused as to what the hell is happening (hopefully I am doing something stupid).

      subroutine myZAverage

      include 'SIZE'

      include 'TOTAL'

      real ua(nz1,nelz),w1(nz1,nelz),w2(nz1,nelz)

      call planar_average_z(ua,vx,w1,w2)

      return

      end

--

Regards,

Akshay Patil

[Akshay Patil]

Hi Guys,

I am kind of lost with planform averaging routines. I updated the Nek5000 version from 17 to 19 and navier5.f no longer has the planar_average_z() subroutine. This is kind of a frustrating since I am trying to get the planform averaging for my channel case and I keep running into silly roadblock which seem highly cyptic. Either I am too stupid to understand the code layout or there is something nasty happening. Could someone tell me how I can do the following?

My geometric layout Channel Box: 0.7H x 0.3H x 0.1H, where H = 0.1m with 80 gridpoints in each direction, i.e. 80x80x80 gridpoints (for the test case)

I want to average (say vx) in x and y direction i.e. streamwise and spanwise direction in my case. This should give me about 80 values of vx averaged over x and y.

Thanks in advance! Apologise my failure to comprehend the simplest of the Nek5000 code bits.

--

Regards,

Akshay Patil

[Akshay Patil]

Just to clarify I was trying the planar_average() subroutine in Nek5000-v19. This subroutine I think assumes that the mesh is extruded in the idir.  Here is my code.

idir = 3 ! z

call gtpp_gs_setup(igs_z,nelx*nely,1,nelz,idir)

call planar_avg(uavg_z,vx,igs_z)

And here is the error it throws me

gtpp_gs_setup: invalid gtp mesh dimensions!

--

Regards,

Akshay Patil

[Fischer, Paul]

Does nelx*nely*nelz = nelgv ?

Paul

---

From: nek...@googlegroups.com <nek...@googlegroups.com> on behalf of Akshay Patil <patil.ak...@gmail.com>
Sent: Wednesday, June 24, 2020 6:34 PM
To: Tanmoy Chatterjee <t.chat...@asu.edu>
Cc: Nek5000 <nek...@googlegroups.com>
Subject: Re: [nek5000] subroutine help

 

To view this discussion on the web visit https://groups.google.com/d/msgid/nek5000/CACqPjWFG_ZwjFohjRN3Hz3VwsyHjALM4kde7uZuLR%2BJu0pNbhg%40mail.gmail.com.

[Akshay Patil]

Hi Paul,

Yes I have 10x10x10 elements in x,y,z directions. I just did a quick check while runtime and nelgv=1000. 

--

Regards,

Akshay Patil

[Fischer, Paul]

Hi Akshay,

cd

cd Nek5000/core

grep "invalid gtp" *

points to dssum.f

Looking at the code there, I see:

      subroutine gtpp_gs_setup(hndl,nelgx,nelgy,nelgz,idir)

      include 'SIZE'

      include 'TOTAL'

      integer hndl,nelgx,nelgy,nelgz,idir

      common /nekmpi/ mid,mp,nekcomm,nekgroup,nekreal

      common /c_is1/ glo_num(lx1,ly1,lz1,lelv)

      integer e,ex,ey,ez,eg

      integer*8 glo_num,ex_g

      nelgxyz = nelgx*nelgy*nelgz

      if (nelgxyz .ne. nelgv)

     $ call exitti('gtpp_gs_setup: invalid gtp mesh dimensions!$',

     $             nelgxyz)

hth,

Paul

---

From: Akshay Patil <patil.ak...@gmail.com>
Sent: Wednesday, June 24, 2020 7:23 PM
To: Fischer, Paul <fisc...@illinois.edu>
Cc: Tanmoy Chatterjee <t.chat...@asu.edu>; Nek5000 <nek...@googlegroups.com>

[Akshay Patil]

I have no idea why the error shows up.  When I tried to print it during the simulation run it indeed gives nelgxyz .ne. nelgv. I used the code below to confirm

      if(nid .eq. 0) then

            print *,  'nelgv:', nelgv, 'nelgxyz:', nelgxyz

      endif

The output is nelgv: 1000 nelgxyz: 21966

I have -10 -10 -10 in my channel.box file. 

--

Regards,

Akshay Patil

[Fischer, Paul]

Did you explicitly set

nelx=10

nely=10

nelz=10

in your .usr file ?  (The user file knows nothing about the box file...)

Best,

Paul

---

From: Akshay Patil <patil.ak...@gmail.com>
Sent: Wednesday, June 24, 2020 7:41 PM

[Akshay Patil]

Alright that was the issue, I had to prescribe what nelx, nely,nelz were. Thanks for that. So now when I try the following code without defining the size of uavg_z it gives me a segmentation fault (My understanding thats an invalid memory/pointer issue). Here is the code I used. 

      subroutine userchk()

      include 'SIZE'

      include 'TOTAL'

      include 'ZPER'

c     Define the number of elements since userchk() wont know whats in channel.box

      nelx=10

      nely=10

      nelz=10

c     Define the direction over which averaging to be done

      idir = 3 ! z

      call gtpp_gs_setup(igs_z,nelx*nely,1,nelz,idir)

      call planar_avg(uavg_z,vx,igs_z)

      return

      end

--

Regards,

Akshay Patil

[Akshay Patil]

Oh I think I have not defined/prescribed the size of uavg_z which is why it does not know where to store it. I have it working now. If you dont mind can I ask another question with regards to the same averaging utility. 

What is the best way to compute velocity gradients and TKE budget terms? (I looked into what Philipp shared yesterday and could not really fathom the KTH utility). Additionally, I think personally its a good exercise to write the subroutines myself that way I can learn a litle bit of the Nek5000 framework better. 

Thanks for the help! Highly appreciate it!

--

Regards,

Akshay Patil

[Tanmoy Chatterjee]

You calculated uavg using glsc2(), which is a global summation using gop or gather operation, so all the processors have the global scalar value. the if(nid .eq. 0) condition ensures that you are just printing for rank = 0.

Regarding the planar average routines if you are using v19 instead of v17 (quite a bit of code change), then, yes please follow Dr. Fischer's advice since the planar avg routines are made more generic to incorporate extruded meshes rather than just box meshes. In general, if you are using planar avg routines its imperative that you explicitly specify nelx, nely, nelz.

Best Regards,

Tanmoy

To view this discussion on the web visit https://groups.google.com/d/msgid/nek5000/CACqPjWFF%3D0DgCeb6PEsJBFn8tGf-obeCX5p2Yd71CUp7wCXX3Q%40mail.gmail.com.

[Akshay Patil]

Hi Tanmoy,

I see, I did not expect the changes from v19 and v17. I realised only today that my local machine had v19 and my server had v17 which causes all the mess. I apologise for my rather frustrated tone. I had no idea of that this was the case. It was only at the end of the day that I realised this was the case. Again sorry about the fussy issues, thanks for understanding!

I think I have the hang of it now. One last thing, could you comment on computing TKE budgets (again planform averages such that I get a vertical profile of the TKE budgets). Thanks!

--

Regards,

Akshay Patil

[Tanmoy Chatterjee]

Hi Akshay, 

You can definitely write your own budget equations.

For gradients try using gradm1() which will give you the strong gradients. This should be fine if you are using DNS (highly resolved simulations)

call gradm1(vx1,vxx,vxy,vxz); it will give you gradients in all the three directions. TKE budget will invoke a lot of gradients you would want to design your code such that one array can be used for many gradients (by overwriting them carefully) and storing for the average. 

Additionally. you will need avg1(), avg2(), avg3(), and you can write one routine avg4() in userchk for average of triple products. You can look at the workflow at navier5.f

Best Regards,

Tanmoy

[Akshay Patil]

I will keep that in mind! Thanks a lot again. I hope its fine if I chime on this thread if I face issues? 

--

Regards,

Akshay Patil

[Tanmoy Chatterjee]

sure.

Best Regards,

Tanmoy

To view this discussion on the web visit https://groups.google.com/d/msgid/nek5000/CACqPjWEEdzWzYAmOO0pjS3KyOnh7MvHUHEzj2YSFNZNX%2Bwwrag%40mail.gmail.com.

[Akshay Patil]

Hi Tanmoy,

I have a followup question based on what you suggested yesterday.  I do not fully understand how the data structure looks like.

If we were to say average over the y direction is this how you do it?

     real uavg_y(lx1*ly1*lz1*lelt)

     idir = 2 ! y

     call gtpp_gs_setup(igs_y,nelx*nelz,nely,1,idir)

     call planar_avg(uavg_y,vx,igs_y)

This should give uavg_y which is size lx1*lz1*lelt (because we averaged over the y direction correct?). However, thats not the case since we need to initialise the uavg_y of size lx1*ly1*lz1*lelt or as you mentioned size(vx) = size(uavg_y). So I am a bit confused as to how this average routine works.

A followup question to the above one is that:

Now if I average uavg_y over the x direction it should give me a result (say) uavg_yx which is of size lz1*lelt. For that based on the discussion yesterday this is something I would use?

     real uavg_y(lx1*ly1*lz1*lelt), uavg_yx(lx1*ly1*lz1*lelt)

c First average in y and then use that result to average in x

c This way the only values that we get are a function of z (vertical coordinate in my case)

     idir = 2 ! y

     call gtpp_gs_setup(igs_y,nelx*nelz,nely,1,idir)

     call planar_avg(uavg_y,vx,igs_y)

     idir = 1 ! x

     call gtpp_gs_setup(igs_x,nelx,nely,nelz,idir)

     call planar_avg(uavg_yx,uavg_y,igs_x)

--

Regards,

Akshay Patil

[Akshay Patil]

Any thoughts on this? Thanks!

--

Regards,

Akshay Patil

[Akshay Patil]

I adapted the turbChannel 1D profiling to my case and only changed only the y direction averaging (because my coordinate system is z vertical). However, when I use the code below, I get an error! I have highlighted the stuff I changed  so its easy to isolate the issue.

c     Define the call for averaging

      integer icalld

      save    icalld

      data    icalld /0/

c     Define the averaging parameters

      real atime,timel

      save atime,timel

c     Define the data dumping parameters

      integer ntdump

      save    ntdump

c     Define the interpolation stuff

      real    rwk(INTP_NMAX,ldim+1) ! r, s, t, dist2

      integer iwk(INTP_NMAX,3)      ! code, proc, el 

      save    rwk, iwk

      integer nint, intp_h

      save    nint, intp_h

      logical iffpts

      save iffpts

c     Define the working array based on the interpolation point number

      real xint(INTP_NMAX),yint(INTP_NMAX),zint(INTP_NMAX)

      save xint, yint, zint

      save igs_x, igs_y

c     Define the parameters and the array sizes

      parameter(nstat=9)

      real ravg(lx1*ly1*lz1*lelt,nstat)

      real stat(lx1*ly1*lz1*lelt,nstat)

      real stat_z(INTP_NMAX*nstat)

      save ravg, stat, stat_z

c     Define the working array and logical parameter to save

      logical ifverbose

      common /gaaa/    wo1(lx1,ly1,lz1,lelv)

c     Define a fixed value of the shear velocity

      real u_tau

      u_tau = 5.3e-3

c     Define the number of points

      n     = nx1*ny1*nz1*nelv

c     Define the number of elements in each direction

      nelx  = NUM_X

      nely  = NUM_Y

      nelz  = NUM_Z

c     Computing the mean velocities

      ubar = glsc2(vx,bm1,n)/volvm1

      e2   = glsc3(vy,bm1,vy,n)+glsc3(vz,bm1,vz,n)

      e2   = e2/volvm1

c     Print the values at each time step

      if(nid.eq.0) write(6,*) 'Time:',time,'Ux:',ubar,'Uy+Uz:',e2

      

      if (time.lt.tSTATSTART) return

c

c     What follows computes some statistics ...

c

      if(icalld.eq.0) then

        if(nid.eq.0) write(6,*) 'Start collecting statistics ...', u_tau

        nxm = 1 ! mesh is linear

        call interp_setup(intp_h,0.0,nxm,nelt)

        nint = 0

        if (nid.eq.0) then

          nint = INTP_NMAX

          call cfill(xint,XCINT,size(xint))

          call cfill(yint,YCINT,size(yint))

          do i = 1,INTP_NMAX 

             zi = (i-1.)/(INTP_NMAX-1)

             zi = 0.1*zi  !Rescaling the z values

             zint(i) = tanh(BETAM*(2*zi-1))/tanh(BETAM)

          enddo

        endif

        iffpts = .true. ! dummy call to find points

        call interp_nfld(stat_z,ravg,1,xint,yint,zint,nint,

     $                        iwk,rwk,INTP_NMAX,iffpts,intp_h)

        iffpts = .false.

        call gtpp_gs_setup(igs_x,nelx     ,nely,nelz,1) ! x-avx

        call gtpp_gs_setup(igs_y,nelx*nelz,nely,1,2) ! y-avg

        call rzero(ravg,size(ravg))

        atime     = 0

        timel     = time

        ntdump    = int(time/tSTATFREQ)

        icalld = 1

      endif

c     Set the time averaging parameters

      dtime = time - timel

      atime = atime + dtime

c     Averaging over time

      if (atime.ne.0. .and. dtime.ne.0.) then

        beta      = dtime / atime

        alpha     = 1. - beta

        ifverbose = .false.

        call avg1(ravg(1,1),vx   ,alpha,beta,n,'uavg',ifverbose)

        call avg2(ravg(1,2),vx   ,alpha,beta,n,'urms',ifverbose)

        call avg2(ravg(1,3),vy   ,alpha,beta,n,'vrms',ifverbose)

        call avg2(ravg(1,4),vz   ,alpha,beta,n,'wrms',ifverbose)

        call avg3(ravg(1,5),vx,vy,alpha,beta,n,'uvmm',ifverbose)

        call avg1(ravg(1,6),t    ,alpha,beta,n,'tavg',ifverbose)

        call avg2(ravg(1,7),t    ,alpha,beta,n,'trms',ifverbose)

        call avg3(ravg(1,8),vx,t ,alpha,beta,n,'utmm',ifverbose)

        call avg3(ravg(1,9),vy,t ,alpha,beta,n,'vtmm',ifverbose)

      endif

c     Reset old time as new time

      timel = time

      ! write statistics to file

      if(istep.gt.0 .and. time.gt.(ntdump+1)*tSTATFREQ) then

         ! averaging over statistical homogeneous directions (x-y)

         do i = 1,nstat

            call planar_avg(wo1      ,ravg(1,i),igs_x)

            call planar_avg(stat(1,i),wo1      ,igs_y)

         enddo

         ! extract data along wall normal direction (1D profile)

         call interp_nfld(stat_z,stat,nstat,xint,yint,zint,nint,

     $                    iwk,rwk,INTP_NMAX,iffpts,intp_h)

         ntdump = ntdump + 1

         if (nid.ne.0) goto 998 

         write(6,*) 'Dumping statistics ...'

 

         open(unit=56,file='vel_fluc_prof.dat')

         write(56,'(A,1pe14.7)') '#time = ', time

         write(56,'(A)') 

     $    '#  z    uu    vv    ww    uv'

         open(unit=57,file='mean_prof.dat')

         write(57,'(A,1pe14.7)') '#time = ', time

         write(57,'(A)') 

     $    '#  z    Umean'

         do i = 1,nint

            zz = 1+zint(i)

            write(56,3) 

     &           zz,

     &           (stat_z(1*nint+i)-(stat_z(0*nint+i))**2),

     &           stat_z(2*nint+i),

     &           stat_z(3*nint+i),

     &           stat_z(4*nint+i)

            write(57,3) 

     &           zz,

     &           stat_z(0*nint+i)

  3         format(1p15e17.9)

         enddo

         close(56)

         close(57)

 998  endif

When I use this, I get the error as follows:

call userchk

 Time:   0.0000000000000000      Ux:   5.7054317786125811E-002 Uy+Uz:   3.3048747827478064E-003

 Start collecting statistics ...   5.3000000000000000E-003

 call intp_setup tol=   4.9999999999999999E-013

 Nmsh for findpts:           1

 WARNING: point not within mesh xy[z]: !  4.2439916-314  4.2439916-314 -1.0000000E+00

 WARNING: point not within mesh xy[z]: !  4.2439916-314  4.2439916-314 -9.9894341E-01

 WARNING: point not within mesh xy[z]: !  4.2439916-314  4.2439916-314 -9.9788331E-01

 WARNING: point not within mesh xy[z]: !  4.2439916-314  4.2439916-314 -9.9681969E-01

 WARNING: point not within mesh xy[z]: !  4.2439916-314  4.2439916-314 -9.9575256E-01

   total number of points =          100

   failed =          100

--

Regards,

Akshay Patil

[Tanmoy Chatterjee]

seems there is a problem with your interp_nfld(). Are your points defined within the domain?

Best Regards,

Tanmoy

To view this discussion on the web visit https://groups.google.com/d/msgid/nek5000/CACqPjWFtj3u9v%3Dj4VwgQh00OnFzw-pPkJY5MDDZnpXFioz-mXQ%40mail.gmail.com.

[Akshay Patil]

Hi Guys,

For a given channel flow with 10x10x10 elements in x, y, and z-direction the planform averaging can be done as follows right?

Given that u(i,j,k,e) is the data structure with FORTRAN column format. With N=8 and Nelx=Nely=Nelz=10, we can say that the data structure is laid out as

z_i = u(80*80,1,1,processor)

So if I want to carryout planform averages I could simply average the values every (N*Nelx)*(N*Nely) data points and obtain the profile for u(z) i.e. obtain a 1D profile for the channel flow, right?

I wrote a code that looks like the one below. This is assuming the above case is true.

c	Define the uavg array which is of a predefined size (in this case N*Nez = 8*10 = 80)

	real uavg(80)
	real numProc

c	Set the number of processor #define NUM_PROC <userValue> defined earlier

	numProc = NUM_PROC

c	Define the total number of data points (on one processor/mpi rank?)

	n = nx1*ny1*nz1*nelv

c	Define the counters which will be used

	zCounter = 1
	counter = 1

c	Main loop which loops over the number of processors (mpi ranks)

	do proc=1,numProc

c		Reset the startt and endd variables to 1 and 80*80 so that it starts at the right values		

		startt = 1
		endd = 80*80

c		Loop over the points

		do i=1,n
			uavg(zCounter) = sum(vx(startt:endd,1,1,numProc))/endd
			if(counter == 80*80) then
c				print *, uavg(zCounter), zCounter	!NOTE
				counter = 1                             !Reset the counter to 1   
				zCounter = zCounter + 1                 !Increment the z position counter 
				startt = i                              !Reset the startt of the array sum   
				endd = startt + 80*80                          
			endif

			counter = counter + 1
		enddo
	enddo 

I have a two-part question about this:

1. This code seems to give an array that repeats after index 40 and I am not sure why that is the case.  What I was expecting it would do is continue sampling the velocity profile for z(40+i) until 40+i == 80. My preliminary guess is that this is happening because I use the sum() function and not the global sum. However, I am not entirely sure how to tackle this particular issue.

2. How can I collect all the values in uavg array and write it to a single file say every Nth time step? I tried doing

if(nid.eq.0) then

      open(unit=10,file='z.dat')

      write(10,*) uavg

endif

However, this seems to only write the first 40 points. My understanding is this is happening because it's only writing the points on the master node. Since I am using 2 mpi ranks to do the test case, that makes sense. However, how can I make it write the entire array and not just the first N*Nelz/mpiRank points?

--

Regards,

Akshay Patil

[Akshay Patil]

Hi guys,

I ended up recycling the turbChannel code which uses interpolation. Here is the code in case anyone needs is using the z-direction as the vertical coordinate instead of the conventional y-direction. Note: That this setup is in dimensional units with H=0.1 m hence a scaling factor was introduced to scale the geometry from 1 m to 0.1m. These definitions are included at the top of the *.usr file. All capital variables are in the #define format (not included below).

c     Define the call for averaging
      integer icalld
      save    icalld
      data    icalld /0/

c     Define the averaging parameters
      real atime,timel
      save atime,timel

c     Define the data dumping parameters
      integer ntdump
      save    ntdump

c     Define the interpolation stuff
      real    rwk(INTP_NMAX,ldim+1) ! r, s, t, dist2
      integer iwk(INTP_NMAX,3)      ! code, proc, el 
      save    rwk, iwk

      integer nint, intp_h
      save    nint, intp_h

      logical iffpts
      save iffpts

c     Define the working array based on the interpolation point number
      real xint(INTP_NMAX),yint(INTP_NMAX),zint(INTP_NMAX)
      save xint, yint, zint
      save igs_x, igs_y


c     Define the parameters and the array sizes
      parameter(nstat=9)
      real ravg(lx1*ly1*lz1*lelt,nstat)
      real stat(lx1*ly1*lz1*lelt,nstat)
      real stat_z(INTP_NMAX*nstat)
      save ravg, stat, stat_z

c     Define the working array and logical parameter to save
      logical ifverbose

      common /gaaa/    wo1(lx1,ly1,lz1,lelv)

c     Define the number of points
      n     = nx1*ny1*nz1*nelv

c     Define the number of elements in each direction
      nelx  = NUM_X
      nely  = NUM_Y
      nelz  = NUM_Z

c     Define the scaling factor
      scalingFactor = SCALEFAC
      
c     Computing the mean velocities

      ubar = glsc2(vx,bm1,n)/volvm1
      e2   = glsc3(vy,bm1,vy,n)+glsc3(vz,bm1,vz,n)
      e2   = e2/volvm1

c     Print the values at each time step
      if(nid.eq.0) write(6,*) 'Time:',time,'Ux:',ubar,'Uy+Uz:',e2
      
      if (time.lt.tSTATSTART) return

c
c     What follows computes some statistics ...
c

      if(icalld.eq.0) then
        if(nid.eq.0) write(6,*) 'Start collecting statistics ...'

        nxm = 1 ! mesh is linear
        call interp_setup(intp_h,0.0,nxm,nelt)
        nint = 0
        if (nid.eq.0) then
          nint = INTP_NMAX
          call cfill(xint,XCINT,size(xint))
          call cfill(yint,YCINT,size(yint))
          do i = 1,INTP_NMAX 
             zi = (i-1.)/(INTP_NMAX-1)
c             zi = 0.1*zi  !Rescaling the z values
c             zint(i) = tanh(BETAM*(2*zi-1))/tanh(BETAM)
             zint(i) = scalingFactor*zi
          enddo
        endif


        iffpts = .true. ! dummy call to find points
        call interp_nfld(stat_z,ravg,1,xint,yint,zint,nint,
     $                        iwk,rwk,INTP_NMAX,iffpts,intp_h)

        iffpts = .false.
        call gtpp_gs_setup(igs_x,nelx*nelz,nely,1,1) ! x-avx
        call gtpp_gs_setup(igs_y,nelx,nely,nelz,2) ! y-avg

        call rzero(ravg,size(ravg))
        atime     = 0
        timel     = time
        ntdump    = int(time/tSTATFREQ)

        icalld = 1
      endif

c     Set the time averaging parameters
      dtime = time - timel
      atime = atime + dtime

c     Averaging over time

      if (atime.ne.0. .and. dtime.ne.0.) then
        beta      = dtime / atime
        alpha     = 1. - beta

        ifverbose = .false.
        call avg1(ravg(1,1),vx   ,alpha,beta,n,'uavg',ifverbose)
        call avg2(ravg(1,2),vx   ,alpha,beta,n,'urms',ifverbose)
        call avg2(ravg(1,3),vy   ,alpha,beta,n,'vrms',ifverbose)
        call avg2(ravg(1,4),vz   ,alpha,beta,n,'wrms',ifverbose)
        call avg3(ravg(1,5),vx,vy,alpha,beta,n,'uvmm',ifverbose)

        call avg1(ravg(1,6),t    ,alpha,beta,n,'tavg',ifverbose)
        call avg2(ravg(1,7),t    ,alpha,beta,n,'trms',ifverbose)
        call avg3(ravg(1,8),vx,t ,alpha,beta,n,'utmm',ifverbose)
        call avg3(ravg(1,9),vy,t ,alpha,beta,n,'vtmm',ifverbose)
      endif

c     Reset old time as new time
      timel = time

      ! write statistics to file
      if(istep.gt.0 .and. time.gt.(ntdump+1)*tSTATFREQ) then
         ! averaging over statistical homogeneous directions (x-y)
         do i = 1,nstat
            call planar_avg(wo1      ,ravg(1,i),igs_x)
            call planar_avg(stat(1,i),wo1      ,igs_y)
         enddo

         ! extract data along wall normal direction (1D profile)
         call interp_nfld(stat_z,stat,nstat,xint,yint,zint,nint,
     $                    iwk,rwk,INTP_NMAX,iffpts,intp_h)

         ntdump = ntdump + 1

         if (nid.ne.0) goto 998 

         write(6,*) 'Dumping statistics ...'
 
         open(unit=56,file='vel_fluc_prof.dat')
         write(56,'(A,1pe14.7)') '#time = ', time
         write(56,'(A)') 
     $    '#  z    uu    vv    ww    uv'

         open(unit=57,file='mean_prof.dat')
         write(57,'(A,1pe14.7)') '#time = ', time
         write(57,'(A)') 
     $    '#  z    Umean'

         do i = 1,nint
c            zz = 1+zint(i)
             zz = zint(i)
            write(56,3) 
     &           zz,
     &           (stat_z(1*nint+i)-(stat_z(0*nint+i))**2),
     &           stat_z(2*nint+i),
     &           stat_z(3*nint+i),
     &           stat_z(4*nint+i)

            write(57,3) 
     &           zz,
     &           stat_z(0*nint+i)

  3         format(1p15e17.9)
         enddo
         close(56)
         close(57)

 998  endif

--

Regards,

Akshay Patil