Does solenoidal field orthogonal to dilatational field?

FMDenaro · April 11, 2019, 04:46

Quote:

Originally Posted by TurbJet

I don't understand. If u and kx is real, wouldn't this result in imaginary du/dx?

If you start from the Fourier series expressing u(x) as sum of U(k)*exp(i*k*x) you see the definition by using a single component

du/dx = U(k) d/dx exp(i*k*x) = i*k*U(k)*exp(i*k*x)= i*k*u(x)

you can think about the derivative of the spectral interpolant, see https://people.maths.ox.ac.uk/trefethen/7all.pdf

A clear explanation of the two approaches is provided in section 2.1.3 of the textbook https://www.springer.com/us/book/9783540307259

see also here, note 4 raised the point about the complex derivative

https://math.mit.edu/~stevenj/fft-deriv.pdf

However, you are working on the divergence-free function so I would use the equation like kx*U(kx)+ky*V(ky)+kz*W(kz)=0

TurbJet · April 11, 2019, 17:44

Quote:

Originally Posted by FMDenaro

If you start from the Fourier series expressing u(x) as sum of U(k)*exp(i*k*x) you see the definition by using a single component

du/dx = U(k) d/dx exp(i*k*x) = i*k*U(k)*exp(i*k*x)= i*k*u(x)

you can think about the derivative of the spectral interpolant, see https://people.maths.ox.ac.uk/trefethen/7all.pdf

A clear explanation of the two approaches is provided in section 2.1.3 of the textbook https://www.springer.com/us/book/9783540307259

see also here, note 4 raised the point about the complex derivative

https://math.mit.edu/~stevenj/fft-deriv.pdf

However, you are working on the divergence-free function so I would use the equation like kx*U(kx)+ky*V(ky)+kz*W(kz)=0

Um, this is new; different from what I learnt. Anyway, I'll take a look at it. Thanks.

But come back to my previous question: after I decompose the velocity field in spectral space and inverse FFT the velocity fields back to physical space, the imaginary parts are not close to zero. I tested my HHD code on Taylor-Green vortex, it works perfectly well; so I don't think it's the issue with my HHD code. Now I am thinking could it due to aliasing? Otherwise, I don't see other possibilities.

FMDenaro · April 11, 2019, 18:04

when you consider the derivative, you get [i*k*U(k)]*exp(i*k*x), the term in the square brackets is the Fourier transform of du/dx, have a look at Eq.(2.6) in the Trefethen notes.

Eifoehn4 · April 12, 2019, 03:21

Quote:

Originally Posted by FMDenaro

If you start from the Fourier series expressing u(x) as sum of U(k)*exp(i*k*x) you see the definition by using a single component

du/dx = U(k) d/dx exp(i*k*x) = i*k*U(k)*exp(i*k*x)= i*k*u(x)

you can think about the derivative of the spectral interpolant, see https://people.maths.ox.ac.uk/trefethen/7all.pdf

A clear explanation of the two approaches is provided in section 2.1.3 of the textbook https://www.springer.com/us/book/9783540307259

see also here, note 4 raised the point about the complex derivative

https://math.mit.edu/~stevenj/fft-deriv.pdf

However, you are working on the divergence-free function so I would use the equation like kx*U(kx)+ky*V(ky)+kz*W(kz)=0

I think you are talking about the same issue. Consider

$U(x)= \sum_{i=0}^N \hat{U}_i\phi_{i}(x)=\sum_{i=0}^N \hat{U}_i e^{ikx}$

being trigonometric expansion with $\phi_{i}$ being your modal basis functions and $\hat{U}_i$ are your modal coefficients (spectral).

You can also define the same polynomial in Lagrange form with

$U(x)= \sum_{i=0}^N \tilde{U}_i l_{i}(x)$

being the Lagrange expansion with $l_{i}$ being your nodal basis functions and $\tilde{U}_i$ are your nodal coefficients (physical)

Here

is the same polynomial in different representation.

You call $\mathcal{F}$ the discrete Fourier-transform and $\mathcal{F}^{-1}$ its inverse.

Formally you can write both

$\frac{dU(x)}{dx} = \mathcal{F}^{-1}(ik_x\mathcal{F}(U(x)))$

or

$\frac{dU(x)}{dx} = ik_x U(x)$

because a polynomial is a polynomial is a polynomial, or to put it another way, it is only a case of representation.

Perhaps it is better to write

$\hat{\boldsymbol{U}}= \mathcal{F}(\tilde{\boldsymbol{U}})$

or vice versa

$\tilde{\boldsymbol{U}} = \mathcal{F}^{-1}(\hat{\boldsymbol{U}})$

with

$\hat{\boldsymbol{U}}=[\hat{U}_0,\hat{U}_1,...,\hat{U}_N]^T$

being the the vector of modal coefficients and

$\tilde{\boldsymbol{U}}=[\tilde{U}_0,\tilde{U}_1,...,\tilde{U}_N]^T$

being the the vector of nodal coefficients

The discrete Fourier-transform only changes the coefficients.

sbaffini · April 12, 2019, 10:17

This is how I would get the spectral derivative dy/dx of a periodic function y(x) in MATLAB for a given set of coordinates x with constant spacing:

%MATLAB SECTION
nx=length(x);
Lx=x(end)-x(1);

y=y(1:end-1); nx=nx-1; %Remove periodicity
Nc=(nx+mod(nx,2))/2-1;
kx1=2*pi*(0:Nc)/Lx;
kx2=2*pi*(-Nc:-1)/Lx;

if mod(nx,2)==1

kx=[kx1 kx2];

else

kx=[kx1 0 kx2];

end

y_hat = fft(y);
dydx_hat = 1i*kx.*y_hat;
dydx = real(ifft(dydx_hat));
dydx = [dydx dydx(1)]; %Add back periodicity
%END OF MATLAB SECTION

where the remove and add-back periodicity are NOT needed if the end points are NOT actually the same. For example, use it for finite difference computations but not finite volume ones.

TurbJet · April 12, 2019, 19:52

Quote:

Originally Posted by Eifoehn4

I think you are talking about the same issue. Consider

$U(x)= \sum_{i=0}^N \hat{U}_i\phi_{i}(x)=\sum_{i=0}^N \hat{U}_i e^{ikx}$

being trigonometric expansion with $\phi_{i}$ being your modal basis functions and $\hat{U}_i$ are your modal coefficients (spectral).

You can also define the same polynomial in Lagrange form with

$U(x)= \sum_{i=0}^N \tilde{U}_i l_{i}(x)$

being the Lagrange expansion with $l_{i}$ being your nodal basis functions and $\tilde{U}_i$ are your nodal coefficients (physical)

Here

is the same polynomial in different representation.

You call $\mathcal{F}$ the discrete Fourier-transform and $\mathcal{F}^{-1}$ its inverse.

Formally you can write both

$\frac{dU(x)}{dx} = \mathcal{F}^{-1}(ik_x\mathcal{F}(U(x)))$

or

$\frac{dU(x)}{dx} = ik_x U(x)$

because a polynomial is a polynomial is a polynomial, or to put it another way, it is only a case of representation.

Perhaps it is better to write

$\hat{\boldsymbol{U}}= \mathcal{F}(\tilde{\boldsymbol{U}})$

or vice versa

$\tilde{\boldsymbol{U}} = \mathcal{F}^{-1}(\hat{\boldsymbol{U}})$

with

$\hat{\boldsymbol{U}}=[\hat{U}_0,\hat{U}_1,...,\hat{U}_N]^T$

being the the vector of modal coefficients and

$\tilde{\boldsymbol{U}}=[\tilde{U}_0,\tilde{U}_1,...,\tilde{U}_N]^T$

being the the vector of nodal coefficients

The discrete Fourier-transform only changes the coefficients.

Yeah, this is what I learned using Fourier transform to compute derivative.

TurbJet · April 12, 2019, 19:54

Quote:

Originally Posted by sbaffini

This is how I would get the spectral derivative dy/dx of a periodic function y(x) in MATLAB for a given set of coordinates x with constant spacing:

%MATLAB SECTION
nx=length(x);
Lx=x(end)-x(1);

y=y(1:end-1); nx=nx-1; %Remove periodicity
Nc=(nx+mod(nx,2))/2-1;
kx1=2*pi*(0:Nc)/Lx;
kx2=2*pi*(-Nc:-1)/Lx;

if mod(nx,2)==1

kx=[kx1 kx2];

else

kx=[kx1 0 kx2];

end

y_hat = fft(y);
dydx_hat = 1i*kx.*y_hat;
dydx = real(ifft(dydx_hat));
dydx = [dydx dydx(1)]; %Add back periodicity
%END OF MATLAB SECTION

where the remove and add-back periodicity are NOT needed if the end points are NOT actually the same. For example, use it for finite difference computations but not finite volume ones.

Yes, this is pretty much how I did the derivative.

By my question was: the solution I got back from iFFT is complex, and the imaginary part is significant. I believe it's due to the aliasing errors. So I am wondering is there anyway I can separate these aliasing errors, analyze them quantitatively?

Eifoehn4 · April 13, 2019, 04:14

Quote:

Originally Posted by TurbJet

Hi,

Currently I am reading a paper about compressible homogeneous isotropic turbulence. In this paper, the authors claim that

$\langle u^s_iu^d_i\rangle = 0$

where <> stands for volume averaging, and superscript s & d represent solenoidal/dilatational velocity fields, which can be obtained from Helmholtz decomposition.

I really don't see how this identity is valid. Solenoidal and dilatational components are orthogonal in spectral space, that's for sure. But I don't think it is equivalent to be orthogonal in physical space. Or I am misunderstanding this identity; it's equal to zero because some other reason?

Can anyone give me some hints?

Appreciate it.

Could you share the Paper you are talking about? On the other hand, we probably talk past each other.

sbaffini · April 13, 2019, 07:42

Quote:

Originally Posted by TurbJet

Yes, this is pretty much how I did the derivative.

By my question was: the solution I got back from iFFT is complex, and the imaginary part is significant. I believe it's due to the aliasing errors. So I am wondering is there anyway I can separate these aliasing errors, analyze them quantitatively?

Well, first of all, I would check that your grid can actually support the functions you are differentiating. You need at least 4 points per wavelength.

In general, I would use a test field to see that everything works as expected, before moving to the actual data you have no control over.

FMDenaro · April 13, 2019, 09:39

Quote:

Originally Posted by TurbJet

Yes, this is pretty much how I did the derivative.

By my question was: the solution I got back from iFFT is complex, and the imaginary part is significant. I believe it's due to the aliasing errors. So I am wondering is there anyway I can separate these aliasing errors, analyze them quantitatively?

Have you checked your FFT and iFFT processes first on a known derivative? What happens if you consider a function like sin(2*pi*k/L) for several values of k?

FMDenaro · April 13, 2019, 09:41

Have also a look to the HHD in spectral space here
http://bohr.physics.berkeley.edu/cla...amclassemf.pdf

TurbJet · April 13, 2019, 19:19

Quote:

Originally Posted by FMDenaro

Have you checked your FFT and iFFT processes first on a known derivative? What happens if you consider a function like sin(2*pi*k/L) for several values of k?

I am pretty sure my HHD function as well as my derivative code work perfectly well. I've tested them on sin/cos function as you mentioned, I've also tested them on Taylor-Green vortex field. The results are all perfectly well.

And I checked the note you attached. It is the same as what I did.

So, I believe the complex issue stems from the aliasing error.

TurbJet · April 13, 2019, 19:21

Quote:

Originally Posted by sbaffini

Well, first of all, I would check that your grid can actually support the functions you are differentiating. You need at least 4 points per wavelength.

In general, I would use a test field to see that everything works as expected, before moving to the actual data you have no control over.

I am pretty sure my HHD and derivative codes work perfectly fine. I've tested them on various analytical functions, such as sin/cos, Taylor-Green vortex, etc.. So the only problem that I can think of is the aliasing errors.

TurbJet · April 13, 2019, 19:23

Quote:

Originally Posted by Eifoehn4

Could you share the Paper you are talking about? On the other hand, we probably talk past each other.

The paper is
Peterson, Livescu, Forcing for statistically stationary compressible isotropic turbulence. PoF (2010).

And I made a mistake; it actually should be

$\langle w_i^sw_i^d\rangle = 0$

where $w_i^\alpha = \sqrt{\rho}u_i^\alpha \;\;,\; \alpha = s,d$

And yeah, the topic kinda drifted away....

sbaffini · April 13, 2019, 19:33

Quote:

Originally Posted by TurbJet

I pretty sure my HHD and derivative codes work perfectly fine. I've tested them on varies analytically functions, such as sin/cos, Taylor-Green vortex, etc.. So the only problem that I can think of is the aliasing errors.

As final check you could purpotedly feed your routine with an analytical function that you know would be aliased.

However, Taylor solution, if I recall well, is soleinodal. You should try something that is dilatational as well, if you haven't already.

TurbJet · April 13, 2019, 19:34

Quote:

Originally Posted by sbaffini

As final check you could purpotedly feed your routine with an analytical function that you know would be aliased.

However, Taylor solution, if I recall well, is soleinodal. You should try something that is dilatational as well, if you haven't already.

I also constructed a purely dilatational field, and it worked well.

As for the aliased field, do you happen to know any?

sbaffini · April 13, 2019, 19:42

In 2d the taylor solution can be parameterized with respect to the wavelength. You can check my phd thesis here on my blog for it (but it is easy to derive as well). Just use a wavelength higher than what your resolution can sustain

FMDenaro · April 14, 2019, 04:53

Quote:

Originally Posted by TurbJet

I am pretty sure my HHD function as well as my derivative code work perfectly well. I've tested them on sin/cos function as you mentioned, I've also tested them on Taylor-Green vortex field. The results are all perfectly well.

And I checked the note you attached. It is the same as what I did.

So, I believe the complex issue stems from the aliasing error.

To assess if the problem is in the aliasing you should replicate a controlled test. On a grid of size h and Nyquist frequency Kc=pi/h, sample the function sin (k*x) for k<Kc and k>Kc and check the resulting FFT/iFFT for the function and derivative

TurbJet · April 14, 2019, 22:40

Quote:

Originally Posted by FMDenaro

To assess if the problem is in the aliasing you should replicate a controlled test. On a grid of size h and Nyquist frequency Kc=pi/h, sample the function sin (k*x) for k<Kc and k>Kc and check the resulting FFT/iFFT for the function and derivative

I've tried this. This is what I did. For example, I have the following MATLAB code

Code:

clearvars;  clc;    close all

N = 64;
L = 2*pi;
dx = L / N;

x = 0 : dx : L;

w = 33;
y = sin(w * x);


% FFT/iFFT
yk = fft(y(1:end-1));
yp = ifft(yk);


% derivative
k = [0:N/2-1, -N/2:-1];
dydxhat = 1i * k .* yk;

dydx = ifft(dydxhat);

The wavenumber of y is larger than the Nyquist. But the resulting y is practically zero. I can't even build the original function.

FMDenaro · April 15, 2019, 04:34

"But the resulting y is practically zero"

what do you mean? I run your code at w=1,16,33 and checked the y and yp vectors and they are not zero

April 12, 2019, 10:17		#25
sbaffini Senior Member Paolo Lampitella Join Date: Mar 2009 Location: Italy Posts: 2,195 Blog Entries: 29 Rep Power: 39	This is how I would get the spectral derivative dy/dx of a periodic function y(x) in MATLAB for a given set of coordinates x with constant spacing: %MATLAB SECTION nx=length(x); Lx=x(end)-x(1); y=y(1:end-1); nx=nx-1; %Remove periodicity Nc=(nx+mod(nx,2))/2-1; kx1=2pi(0:Nc)/Lx; kx2=2pi(-Nc:-1)/Lx; if mod(nx,2)==1 kx=[kx1 kx2]; else kx=[kx1 0 kx2]; end y_hat = fft(y); dydx_hat = 1ikx.y_hat; dydx = real(ifft(dydx_hat)); dydx = [dydx dydx(1)]; %Add back periodicity %END OF MATLAB SECTION where the remove and add-back periodicity are NOT needed if the end points are NOT actually the same. For example, use it for finite difference computations but not finite volume ones.

Similar Threads
Thread	Thread Starter	Forum	Replies	Last Post
simpleDyMFoam on v1806	gamemakerh	OpenFOAM Programming & Development	0	November 8, 2018 09:15
potential flows, helmholtz decomposition and other stuffs	pigna	Main CFD Forum	1	October 26, 2017 09:34
Moving mesh	Niklas Wikstrom (Wikstrom)	OpenFOAM Running, Solving & CFD	122	June 15, 2014 07:20
Zero size field	taranov	OpenFOAM Bugs	2	April 20, 2010 05:51
Problem with rhoSimpleFoam	matteo_gautero	OpenFOAM Running, Solving & CFD	0	February 28, 2008 07:51

April 11, 2019, 18:04		#23
FMDenaro Senior Member Filippo Maria Denaro Join Date: Jul 2010 Posts: 6,896 Rep Power: 73	when you consider the derivative, you get [ikU(k)]exp(ik*x), the term in the square brackets is the Fourier transform of du/dx, have a look at Eq.(2.6) in the Trefethen notes.

April 13, 2019, 09:41		#31
FMDenaro Senior Member Filippo Maria Denaro Join Date: Jul 2010 Posts: 6,896 Rep Power: 73	Have also a look to the HHD in spectral space here http://bohr.physics.berkeley.edu/cla...amclassemf.pdf

April 13, 2019, 19:42		#37
sbaffini Senior Member Paolo Lampitella Join Date: Mar 2009 Location: Italy Posts: 2,195 Blog Entries: 29 Rep Power: 39	In 2d the taylor solution can be parameterized with respect to the wavelength. You can check my phd thesis here on my blog for it (but it is easy to derive as well). Just use a wavelength higher than what your resolution can sustain

April 15, 2019, 04:34		#40
FMDenaro Senior Member Filippo Maria Denaro Join Date: Jul 2010 Posts: 6,896 Rep Power: 73	"But the resulting y is practically zero" what do you mean? I run your code at w=1,16,33 and checked the y and yp vectors and they are not zero