[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]

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