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

I have added an interpolation tool to the repo, int_tp, in
nek5_svn/trunk/tools/int_tp

There is also a README in that directory with instructions on the details
of using the tool.

(Note - int_tp should be compiled as any standard nek tool in the
trunk/tools/ directory with "maketools int_tp" or "maketools all")

Let me know if you encounter complications.

Thanks,
Katie

On Fri, May 4, 2012 at 7:51 AM, <nek5000-users at lists.mcs.anl.gov> wrote:

> Le 04/05/2012 13:55, nek5000-users at lists.mcs.anl.**gov<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<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<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 <Nek5000-users at lists.mcs.anl.gov>
>>>> https://lists.mcs.anl.gov/**mailman/listinfo/nek5000-users<https://lists.mcs.anl.gov/mailman/listinfo/nek5000-users>
>>>>
>>> ______________________________**_________________
>> Nek5000-users mailing list
>> Nek5000-users at lists.mcs.anl.**gov <Nek5000-users at lists.mcs.anl.gov>
>> https://lists.mcs.anl.gov/**mailman/listinfo/nek5000-users<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 <Nek5000-users at lists.mcs.anl.gov>
> https://lists.mcs.anl.gov/**mailman/listinfo/nek5000-users<https://lists.mcs.anl.gov/mailman/listinfo/nek5000-users>
>
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.mcs.anl.gov/pipermail/nek5000-users/attachments/20120507/cf8deff1/attachment.html>