[Nek5000-users] Transfer data from an nelx*nely*nelz array of Nth-order to an n*n*n array for FFT analysis

[nek5000-users at lists.mcs.anl.gov]

Hi All,

In regards to the discussion on this issue. I have been lately  
interpolating  from the SEM nek5000 grid onto an equidistant grid and  
performing FFT to obtain two-point correlations and energy spectrum  
for pipe flow simulations.

During the simulation, I outpost fields from nek5000 containing the  
velocity and pressure. After the simulation ends, I would use a  
simplified version of nek5000 where I loop on these files and  
interpolate in parallel. Upon reading each file, I interpolate to a  
cylindrical mesh that is equidistant in the axial and azimuthal  
direction, and non-uniform in the radial direction.

Since my grid is huge, I do the interpolation in tubes. I loop in the
radial direction, and for each radial position, I generate the
equidistant (z,theta)-grid then interpolate the velocities and
pressure etc.. onto that grid using subroutines used in hpts. Each
data in a tube is sent to a separate processors where I perform FFT
analysis. The FFT analysis is performed using FFT subroutines included
in a .f file which I added to nek5000.

This procedure continues until I finish looping on all my files.

Best
George


Quoting nek5000-users at lists.mcs.anl.gov:

> Le 04/05/2012 13:55, nek5000-users at lists.mcs.anl.gov a écrit :
>> Hi Neks,
>> another point I'm (even more) interested in is: how do you do your  
>> FFT analysis?
>> I guess in a first step the grid is mapped onto a regular grid
>> (ifreguo = .true.) and in a second step the FFT is computed at that
>> grid.
>> Then the question arises: What is a good choice for "nrg" (with nrg as
>> the dimension of regular grid (nrg**ndim)...)? nrg=lx1 or lx1-1?
>> Regards
>> Jan
>>
>> 2012/4/13 Jan Frielinghausen<jan.frielinghausen at googlemail.com>:
>>> Hi Neks, Hi Can,
>>> I too would be very interested in getting more information about  
>>> these routines.
>>> regards
>>> Jan
>>>
>>> Am 5. April 2012 18:10 schrieb<nek5000-users at lists.mcs.anl.gov>:
>>>> Hi neks !
>>>>
>>>> For a given velocity field
>>>> (vx(1:lx1,1:ly1,1:lz1,1:lelv),vy(1:lx1,1:ly1,1:lz1,1:lelv),vz(1:lx1,1:ly1,1:lz1,1:lelv))
>>>> i would like to calculate the energy spectrum, this implies FFT
>>>> transformation on regular arrays
>>>> vx(1:n,1:n,1:n),vy(1:n,1:n,1:n),vz(1:n,1:n,1:n) (need
>>>> deltax=deltay=deltaz=cst).
>>>> On the emailing-list it has been said :
>>>> "We do have off-line serial routines
>>>> that can transfer data from an nelx x nely x nelz array of Nth-order
>>>> elements to an nxnxn array - this allows for analysis by FFTs etc"
>>>>
>>>>
>>>> Can you please tell me more about this ? What are these routines ??
>>>> B.regards
>>>> Can
>>>>
>>>> _______________________________________________
>>>> Nek5000-users mailing list
>>>> Nek5000-users at lists.mcs.anl.gov
>>>> https://lists.mcs.anl.gov/mailman/listinfo/nek5000-users
>> _______________________________________________
>> Nek5000-users mailing list
>> Nek5000-users at lists.mcs.anl.gov
>> https://lists.mcs.anl.gov/mailman/listinfo/nek5000-users
>>
> The solution (i think its more an hack than a real solution) i  
> found, was to create an hpts.in file with the coordinates of every  
> points of the regular mesh
> i wrote the fallowing routine for this :
> (My domain is a cube of length 2pi*2pi*2pi, you can surely adapt the  
> fallowing routine for other simple geometries)
>
> subroutine generehptsin(ndummy,ldummy,uprn)
>     implicit none
> c  ndummy = number of grid points per dimension
> c  ldummy = length in one dimension
> c  uprn = dummy variable for output
> c  you need nx=ny=nz and deltax=deltay=deltaz
> c  For fft analysis you have to ensure that ndummy**3 is equal to  
> 2^n (ex : 128 points = 2^7)
>     integer ndummy,i,j,k,uprn
>     real*8    ldummy,deltax,
> &    x(0:ndummy-1),y(0:ndummy-1),z(0:ndummy-1)
>     open(uprn,file='hpts.in',status='unknown')
> c    compute deltax
>     deltax=0
>     deltax=ldummy/(ndummy-1)
> c    fil vectors : x,y,z
>     x=0
>     y=0
>     z=0
>     do i=0, ndummy-1
>     x(i)=i*deltax
>     y(i)=i*deltax
>     z(i)=i*deltax
>     enddo
>     write(uprn,*)ndummy*ndummy*ndummy    ! total number of grid  
> points in the first line
>     do k=0, ndummy-1
>     do j=0, ndummy-1
>     do i=0, ndummy-1
>     write(uprn,'(1p20E15.7)')x(i),y(j),z(k)
>     enddo
>     enddo
>     enddo
>     close(uprn)
>     return
>     end
>
> Then modify the hpts() subroutine slightly (just to be able to  
> access to the "fieldout" matrix  when you call hpts() in userchk())
> Eg :  modify "subroutine hpts()" to "subroutine hpts(fieldout)"  
> (it's the first line of hpts routine in postpro.f)
> You will also probably want to comment the fallowing lines to avoid  
> hpts to write in the hpts.out.
>
>       ! write interpolation results to hpts.out
> c      if(nid.eq.0) then
> c        do ip = 1,npoints
> c           write(50,'(1p20E15.7)') time,
> c &      (fieldout(i,ip), i=1,nflds)
> c        enddo
> c      endif
> c      if(nid.eq.0) write(6,*) 'done :: dump history points'
>
>
> Then, to express the velocity field in a v(i,j,k) layout, i wrote  
> the fallowing routine :
>
>     subroutine reformatfield(nfldm,fieldout,
> &    ndummy,vxrdummy,vyrdummy,vzrdummy)
>
> c  You have to call this in userchk routine without changing the  
> first two input arguments
>     implicit none
> c  fieldout is now the output matrix of the hpts routine, we need this matrix
> c  nfldm is the number of field you resolve
> c  ndummy is the number of grid points per dimension
> c  vxrdummx is the x component of the velocity field in a regular grid
> c  vxrdummy is the y component of the velocity field in a regular grid
> c  vxrdummz is the z component of the velocity field in a regular grid
>     integer ndummy,i,j,k,z,nfldm
>     real     field1(ndummy*ndummy*ndummy),
> &        field2(ndummy*ndummy*ndummy),
> &        field3(ndummy*ndummy*ndummy),
> &        fieldout(nfldm,ndummy*ndummy*ndummy),
> &        vxrdummy(1:ndummy,1:ndummy,1:ndummy),
> &        vyrdummy(1:ndummy,1:ndummy,1:ndummy),
> &        vzrdummy(1:ndummy,1:ndummy,1:ndummy)
>
>
>
>     do i=1,ndummy*ndummy*ndummy
>     field1(i)=fieldout(1,i)    ! vx
>     field2(i)=fieldout(2,i)    ! vy
>     field3(i)=fieldout(3,i)    ! vz
>     enddo
>
>     z=0
>     do k=1,ndummy
>         do j=1,ndummy
>             do i=1,ndummy
>             vxrdummy(i,j,k)=field1(i+z)
>             vyrdummy(i,j,k)=field2(i+z)
>             vzrdummy(i,j,k)=field3(i+z)
>             enddo
>         z=z+ndummy
>         enddo
>     enddo
>     return
>     end
>
> This routine will output a velocity field expressed in a regular  
> grid with v=v(i,j,k) layout. You can, from there, use any FFT  
> routine for fft analysis.
> I used the "four1" routine that you can find on the "numerical recipies" tome
> Sorry for my poor english and  hope that will help
> For your question i think nrg=nelx*lx1 (1D) because, the number of  
> grid points on a 1D gll mesh is nelx * lx1 with lx1=N+1
> B.regards
>
>
>
> _______________________________________________
> Nek5000-users mailing list
> Nek5000-users at lists.mcs.anl.gov
> https://lists.mcs.anl.gov/mailman/listinfo/nek5000-users
>



----------------------------------------------------------------
This message was sent using IMP, the Internet Messaging Program.