Version 23 (modified by 10 years ago) ( diff ) | ,
---|
Spectra Objects
Spectra Objects can be created within the ProblemModuleInit routine in problem.f90. To create spectras you first need to add two USE
statements to your problem.f90
USE Spectras USE Fields
Then in ProblemModuleInit declare a variable pointer of type SpectraDef
TYPE(SpectraDef), POINTER :: Spectra
Then create the Spectra and set the various parameters as in the following example
CALL CreateSpectra(Spectra) ALLOCATE(Spectra%Fields(1)) Spectra%Fields(:)%id=(/Mass_Field/) Spectra%type=SCALAR_SPECT
Then at each process event (currently each frame) a curve file will be generated in the out directory ie. out/Spectra00021.curve
that will contain all of the spectra for that frame. See CurveFiles for more information on plotting curve files in visit.
The above example creates a spectrum of density (or rather density2) which will be 1 entry in the curve file.
You can also create spectra of vector fields like velocity2
CALL CreateSpectra(Spectra) ALLOCATE(Spectra%Fields(3)) Spectra%Fields(:)%id=(/vx_Field, vy_Field, vz_Field/) Spectra%type=VECTOR_SPECT
This will generate 6 curves: One for each component of the vector field, (vx2, vy2, vz2) as well as the total (v2) and the helmholtz decomposed fields (v_sol2 and v_div2). Note the total and the decomposed fields will be labeled vx_total, vx_sol, & vx_div
though they have nothing to do with the 'x' direction. Each curve represents the spectral energy density. So in the case of using the Mass_Field
which corresponds to density, we would have
And in the case of the velocity field
- For more information on the Field sub-object's properties see ProcessingFields
- If you are making several spectras, you can reuse the Spectra Pointer (with or without Nullifying it) by calling
CreateSpectra(Spectra)
for each new spectra. - The spectra are calculated using the PFFT module
There are some additional properties of the Spectra object that can be modified:
- Spectra%level — You can adjust the maximum level of data used to generate the spectra. Usually this is done when there are not enough computational resources to store a fixed grid at the maximum level… By default Spectra%level will use the maximum highest level available though you can manually adjust this
Spectra%level=MaxLevel
- Spectra%dk — By default the fourier transforms are binned with a radial bin size equal to the largest minimum wavenumber per dimension… If the domain is 16x64, then the minimum wavenumbers in x and y are 1/16 and 1/64 respectively. This means that the bin size will be 1/16 - or 4 times the minimum wave number. You can adjust this by setting the bin size dk in terms of the smallest resolvable wavenumber. So if dk = 1, the bin size would be 1/64, but there will be aliasing effects every 4 points which correspond to multiples of 1/16.
- Spectra%mB — By default the spectra is taken of the whole computational domain, however you can instead take a spectra of a subsection of the domain by specifying the bounds in the index space of the spectra's level. It is also a good idea to have the domain be multiples of 2 and to coincide with coarse cell boundaries. So typically it is good to first determine the size of the region in coarse cells
N
that you want to take the spectra of, and then assuming that the Spectra%level = MaxLevel, calculate the bounds asSpectra%mB(:,1)=levels(MaxLevel)%mX/2 - (N/2)*2**MaxLevel+1 Spectra%mB(:,2)=levels(MaxLevel)%mX/2 + (N/2)*2**MaxLevel
- Spectra%method — By default the interpolation method for producing the fixed grid data is constant interp. If the coarser grids have strong unresolved gradients, this can produce large signals in the resulting spectra at wavenumbers correposponding to grid sizes. To help mitigate this, consider setting the prolongation method to be PARABOLIC_INTERP or SPECTRAL_PROLONGATION
Spectra%method = SPECTRAL_PROLONGATION
- Spectra%WindowFunction — Finally if the data is not periodic, you can specify a window function to damp strong gradients that will appear at the periodic boundaries. Currently the only option available is a cosine window
Spectra%window=COSINE_WINDOW
Here is a full list of the various Spectra parameters with the default values in brackets:
TYPE SpectraDef TYPE(FieldDef), ALLOCATABLE :: Fields(:) REAL(KIND=qPREC), DIMENSION(:,:), ALLOCATABLE :: Data INTEGER :: Level=MAXIMUMSPECTRALEVEL INTEGER, DIMENSION(3,2) :: mB INTEGER :: WindowFunction = NO_WINDOW REAL(KIND=qPREC) :: dk = DEFAULT INTEGER :: Type = SCALAR_SPECT INTEGER :: method = CONSTANT_INTERP TYPE(SpectraDef), POINTER :: next END TYPE
What does the Spectra calculate?
First we must understand what the FFT calculates
What does the FFT calculate
The FFT performs a Discrete Fourier Transform (DFT) which is defined as
where
is any integer.We can see how this corresponds to a real fourier transform by making the substitutions
, , and and take the limit as
So that the DFT corresponds to a fourier transform where the physical wavelength is normalized by
and the magnitude is normalized byNow the FFT returns the array of transforms with
However,
is periodic in since
We can visually see this by plotting the discrete sampling of the real part of the continuous functions for
and for a grid whereSo we can interpret
as representing a signal with wavenumber or , or . Since we don't expect there to be signals in our data at frequencies above the Nyquist frequency - it is logical to assume that corresponds to a signal with .Mathematically this means that if
then or equivalentlyThis is why we call the resulting FFT array to be in wrap around order. That is, the array returned by the FFT corresponds to normalized wavenumbers of
Or if we index the FFT array starting at 1, and then plot the index of each component on the corresponding grid in k space we get
MultiDimensional FFTs
A multidimensional DFT can be thought of as a sequence of 1D transforms in the way that a multidimensional Fourier transform is a sequence of 1D integrals. That is
And again, the grid will be wrapped around in both x and y. So a 10x10 transform will have components that correspond to the following points in k space.
Now since the normalization of the wavenumber goes like
, if the region is not a square/cube, the spacing of the points in k space will not be uniform.For example, if
thenSo if our domain is 10x5, and we use
to normalize our wave vectors, then we would have the following components.And then each component of the transform is squared and added to the appropriate radial bin.
Attachments (14)
- SampleSpectra.jpeg (62.9 KB ) - added by 13 years ago.
-
Screen Shot 2015-01-06 at 3.31.14 PM.png
(191.0 KB
) - added by 10 years ago.
analytic
-
Screen Shot 2015-01-06 at 3.38.19 PM.png
(74.0 KB
) - added by 10 years ago.
discrete1
-
Screen Shot 2015-01-06 at 3.08.22 PM.png
(122.5 KB
) - added by 10 years ago.
discrete2
- Screen Shot 2015-01-06 at 4.28.36 PM.png (78.7 KB ) - added by 10 years ago.
- Screen Shot 2015-01-13 at 10.49.34 AM.png (53.1 KB ) - added by 10 years ago.
- Screen Shot 2015-01-13 at 10.57.02 AM.png (52.4 KB ) - added by 10 years ago.
- Screen Shot 2015-01-13 at 12.05.47 PM.png (15.8 KB ) - added by 10 years ago.
- Screen Shot 2015-01-13 at 2.44.52 PM.png (93.3 KB ) - added by 10 years ago.
- Screen Shot 2015-01-13 at 2.53.55 PM.png (63.4 KB ) - added by 10 years ago.
- Screen Shot 2015-01-13 at 2.58.13 PM.png (107.9 KB ) - added by 10 years ago.
- Screen Shot 2015-01-13 at 2.57.25 PM.png (75.6 KB ) - added by 10 years ago.
- Screen Shot 2015-01-13 at 3.07.07 PM.png (565.2 KB ) - added by 10 years ago.
- Screen Shot 2015-01-13 at 3.15.15 PM.png (485.0 KB ) - added by 10 years ago.
Download all attachments as: .zip